Weak Field Collapse in AdS: Introducing a Charge Density

We study the effect of a non-vanishing chemical potential on the thermalization time of a strongly coupled large $N_c$ gauge theory in $(2+1)$-dimensions, using a specific bottom-up gravity model in asymptotically AdS space. We first construct a perturbative solution to the gravity-equations, which dynamically interpolates between two AdS black hole backgrounds with different temperatures and chemical potentials, in a perturbative expansion of a bulk neutral scalar field. In the dual field theory, this corresponds to a quench dynamics by a marginal operator, where the corresponding coupling serves as the small parameter in which the perturbation is carried out. The evolution of non-local observables, such as the entanglement entropy, suggests that thermalization time decreases with increasing chemical potential. We also comment on the validity of our perturbative analysis.


Introduction
The AdS/CFT correspondence describes an equivalence between a classical gravitational dynamics and a large N c gauge theory. This remarkable correspondence has proved to be very useful in addressing aspects of strongly coupled dynamics in various models, ranging from understanding aspects of strongly coupled Quantum chromodynamics (QCD) to condensed matter-inspired systems. See e.g. [1,2] for recent reviews on some of these attempts.
Most of such endeavours usually discuss equilibrium properties of a class of strongly coupled SU(N c ) gauge theories at large N c . However, since this is a correspondence between the path integrals of the corresponding Quantum Field Theory (QFT) and the classical Gravity description, it is natural to assume that it extends to time-dependent dynamical situations as well. In fact, dynamical processes in a prototypical field theory model are extremely interesting to explore and learn about, since we do not have a good understanding of the governing rules and laws for systems completely out-of-equilibrium, specially at strong coupling. In this regime, AdS/CFT correspondence is potentially a very useful tool.
Such time-dependent issues fall under two broad categories: the study of quantum quenches, where a system is prepared in an energy eigenstate of a given Hamiltonian. The Hamiltonian is then perturbed by a time-dependent external parameter. Recent developments in cold atoms experiments have initiated a very active research where this perturbation occurs abruptly, see e.g. [3] for a condensed-mater approach and [4,5,6] for a review of the holographic attempts. The other broad issue is to understand the physics of thermalization for strongly coupled system. See e.g. [7,8,9] for earlier works on this. More recently, there has been a renewed interest to understand the physics of thermalization at strong coupling to shed light on the physics of the quark-gluon plasma (QGP) at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC). Most of these works rely heavily on numerical efforts that study a black hole formation process in an asymptotically anti de-Sitter (AdS) space, see e.g. [10,11,12,13,14,15,16,17,18,19,20,21,22].
A generic thermalization process typically describes the dynamical evolution from a "low temperature" phase to a "high temperature" thermal state, where the evolution is highly non-trivial and, in the case of a black hole formation process, highly non-linear as well. Holographically, such a process in asymptotically AdS space can be set up by turning on a non-normalizable (or a normalizable) mode of a bulk scalar field; as a result a shell of the corresponding field collapses in AdS space. It was shown in [23], using a weak-field perturbation method, that if the boundary non-normalizable mode is chosen to be coordinate independent and only have support over a brief time interval, the collapse of the corresponding homogeneous wave will always lead to black hole formation. 1 On the other hand, an alternate approach is to phenomenologically model the black hole formation process with as much simplicity as possible, such that the corresponding geometry can be probed to learn further physics. The hope is to learn at least qualitatively useful lessons which are presumably not heavily dependent on the details of the model. In the present context this can be achieved by exploring the AdS-Vaidya background, which describes a smooth evolving geometry from an empty AdS to an AdS-Schwarzschild background. Gravitationally, this geometry describes the collapse of a null dust in an asymptotically AdS-space. Probing this geometry has already led to numerous interesting results, see e.g. [29,30,31,32,33,34,35], where the 1 If the non-normalizable mode is not turned on, there is a large class of homogeneous, spherically symmetric, initial conditions whose temporal evolution does not lead to black hole formation [25,26,27]. These bulk solutions correspond to states of the boundary conformal field theory (CFT) that fail to thermalize at late times. These solutions have recently been proposed as the analogue-dual of the "quantum revivals" observed in finite size isolated quantum systems and widely studied in the condensed matter literature [28].
behavior of various non-local observables in such a dynamical geometry has been explored.
In [23], an interesting bridge between these two approaches has been established. The authors studied a collapse process for a massless scalar field in a so called "weak field approximation" limit, where the amplitude of the perturbation was chosen small and a perturbative solution of the Einstein's equations was obtained. At the leading order, this solution takes the form of an AdS-Vaidya background which is characterized by one mass function that interpolates between an AdS vacuum to an AdS-Schwarzschild geometry. In the dual field theory side, this corresponds to a dynamical evolution from zero temperature ground state to a thermal state of a certain large N gauge theory (such as the N = 4 super Yang-Mills or the ABJM theory). Thus, at least at the leading order, the analysis of [23] validates the AdS-Vaidya-based phenomenology from a first principle gravity computation, see [24] for more details.
Motivated by this observation, we will explore the possibility of introducing a conserved charge in the boundary theory -the simplest case of an U(1)-charge in this article -in an analogue of the "weak field approximation" limit starting from an effective gravity action. Our motivation is to address how thermalization time is affected in the presence of global conserved charges in a strongly coupled system, with as much analytical control as possible. This may be relevant for understanding the effect of a non-vanishing chemical potential (or a finite-density) on the strong coupling dynamics of out-of-equilibrium QCD at RHIC or at LHC. 2 We will thus generalize the construction of [23] introducing a chemical potential, which -albeit in a suitable approximation -will provide an analytical control on the background. This also provides us with a model in which non-local observables can subsequently be studied to explore the behaviour of thermalization time, in the spirit of [33,34].
Our results may be of general interest, beyond the QGP physics. For example, a qualitative behaviour, if it is universal, may shed light in condensed matter systems which are typically accompanied by a non-vanishing chemical potential. Moreover, to the best of our knowledge, there is no known field theoretic result about how thermalization time scales in a strongly coupled finite density system. Thus it is useful to explore models where this possibility can be realized.
This article is divided in the following parts: We begin with a summary of our results and a brief discussion on those in the next section. In section 3 we present the effective gravity model and provide the details of the perturbative solution. We then discuss the initial and the final states in details in section 4, and subsequently discuss the behaviour of the thermalization time in section 5. Finally we conclude in section 6.

Summary and Discussion
Let us begin with a more specific description. We will discuss a thermal quench, and for simplicity we will restrict ourselves to a (2 + 1)-dimensional large N c gauge theory in the strong coupling regime. Now, consider the Hamiltonian (or the Lagrangian) of the system, denoted by H 0 (or L 0 ), which is perturbed by a time-dependent perturbation of the form where H 0 (or L 0 ) describes the Hamiltonian (or the Lagrangian) of the original quantum field theory, λ(t) is an external parameter and δH ∆ (or O ∆ ) corresponds to the deformation of the QFT by an operator of dimension ∆.
Here we will restrict ourselves completely on asymptotically locally AdS-spaces, which implies that there is an underlying conformal field theory (CFT) that governs the physics. As a first attempt, we restrict ourselves to the case of ∆ = d, which here becomes ∆ = 3, i.e. an exactly marginal deformation. In principle, it is possible to study the quench dynamics by a relevant operator as well. However, we know that relevant operators can trigger an RG-flow and there may be a new CFT -or perhaps a non-relativistic cousin of it -that emerges in the infrared. In such a situation, unless we consider a temperature scale much larger than the scale set by the relevant operator, the black hole formation process may be governed by this infrared geometry instead. 3 To avoid this possible subtlety for relevant operators, we will consider an exactly marginal operator, which does not require a hierarchy between the RG-scale and the temperature-scale. Thus the underlying CFT will remain the same and in the gravitational dual it will suffice for us to specify the asymptotically locally AdS condition with a given radius of curvature as the boundary condition. Now we need to specify the initial conditions. Typically this is specified at t → −∞ (in which limit λ → 0) as a particular energy eigenstate of the QFT Hamiltonian H 0 . As the new coupling is turned on, depending on whether the process is adiabatic or abrupt, the system is expected to evolve differently with the Hamiltonian H λ . In a quantum mechanical system, under an adiabatic process the system remains in an energy eigenstate whose energy evolves with time following the response of the time-evolution of λ(t). On the other hand, for an abrupt quench, the system evolves in a linear superposition of eigenstates of the new Hamiltonian H λ . Here we will focus on the fast quench only 4 , in which the initial state is macroscopically characterized by {E, O φ , µ, T } initial -with E, O φ , µ and T respectively representing the energy, VEV for the marginal operator, chemical potential and the temperature of the state in consideration -and the final state is macroscopically characterized by {E, O φ , µ, T } final .
To properly account for the scale of measuring dimensionful quantities, we define the following dimensionless parameters: which means that we measure all dimensionful quantities in units of temperature of the corresponding state. Evidently, there can be several hierarchy of scales depending on how κ E , κ O φ and κ µ are parametrically separated. Furthermore, there are a couple of dimensionful parameters associated with the quench process itself: the energy injected (denoted by ∆E), and the duration (denoted by δt) or the rapidity of the quench. Thus we can further define Now, our initial state is characterized by the following parametric regime The quench process, to be amenable to a perturbative analysis, is characterized by Finally, the final state is characterized by We remind the reader that the conditions in (2.4)-(2.6) are specific to our perturbative analysis. We further note that, the perturbative solution that we construct is not similar to the AdS-RN-Vaidya geometry described and analyzed in [34] and thus provides a more generic case-study. It is amusing to further note that the regime of parameters outlined in (2.6) physically implies that we are considering a small chemical potential limit, which is qualitatively similar to the QGP-phase at the LHC. The geometric data describing the corresponding evolution is given by a metric, a gauge field, and a scalar field: {G, A, φ} of the following form

7)
where U EF denotes the ingoing Eddington-Finkelstein (EF) patch, g(z, v) and f (z, v) are two functions that are determined by solving the equations, z is the radial coordinate (in which the boundary is located at z → 0) and v is the EF-ingoing null direction. We obtain a dynamical evolution that can be summarized as: Let us now comment on the evolution of the event-horizon and the apparent horizon. We begin with the notion of the apparent horizon. Following [34], let us define the tangent vectors to the ingoing and outgoing null geodesics The co-dimension two spacelike surface, which is orthogonal to the tangent vectors above, has the following volume element The expansion of this volume element along the ingoing and outgoing null directions are given by (2.14) Here L ± denotes the Lie derivatives along the null direction corresponding to µ ± . Now, we can define the invariant expansion by Θ = θ + θ − which is given by Here z aH denotes the location of the apparent horizon.
On the other hand, the event horizon is a null surface in the background (2.7): which gives Thus, solving (2.15) and (2.17) gives the time-evolution of the apparent and the event horizon, respectively. Evidently, at the initial and the final states they coincide: z H = z aH . During the evolution, since the collapse is sourced by a physically reasonable matter field, the apparent horizon lies behind the event horizon, i.e. z aH > z H in our choice of the radial coordinate. One way to summarize our perturbative solution is to state that the evolution of e.g. the event-horizon can be obtained in a series expansion as follows: where Υ 2n can be determined from a first order differential equation of the form where Ξ is a functional of 20) and, finally, ε ∼ κ 1/2 quench 1. In section 5, we will present a pictorial representation of how the apparent and the event-horizons evole.
It is clear from (2.18) that the asymptotic series captures the physics as long as ||ε 2n Υ (2n) || 1. It will be shown In section 3.3, that this imposes a constraint and our perturbative analysis is valid up to a time-scale t pert ∼ 1/ (∆E) 1/3 . For t t pert , we will need to solve the system of equations numerically, which we will not pursue here. Written in terms of the duration of the pulse, the perturbative treatment is valid up to a time-scale t pert which is given by Thus, by tuning ε 1, we parametrically separate t pert and δt at will, and in this sense our approach is equivalent to considering an "fast quench". In the strict fast quench limit, i.e. t pert /(δt) → ∞, the t > t pert dynamics is completely frozen and we are left with the perturbative solution for all times.
To measure thermalization time, we need to identify suitable observables, which primarily fall in two classes: local and non-local. Note that, unlike the Vaidya-construction, for gaugeinvariant local observables thermalization is not instantaneous in this case. 5 This is buried in the details of the solution in (2.7), since as time varies, the scalar field dynamically acquires a nonzero expectation value. However, such local operators thermalize over a time-scale of O(δt), and does not contain the information about long-range correlations. On the other hand, non-local operators provide a more global information and we will explore the evolution of entanglement entropy in this article. Now, let us comment on a naìve expectation about the scaling of the thermalization time. It is known that for integrable systems, which contains an infinite number of conserved charges, thermalization does not happen. This is intuitively clear, since it becomes unlikely to populate the entire phase space. Furthermore, for a single U(1)-charge, if we consider a bosonic system, increasing the chemical potential will enhance Bose-Einstein condensation and thus will inhibit thermalization. For fermionic degrees of freedom, a higher value of chemical potential is associated with a higher Fermi surface which will subsequently need to be populated to achieve a thermal state. Thus, introducing a conserved charge is likely to have a slowing down of the thermalization time.
However, contrary to the expectation outlined in the above paragraph, we gather evidence that thermalization speeds up for increasing chemical potential in the regime µ/T 1. Thus, modulo the caveats of an effective gravity description and the approximate measures of the thermalization time, we are lead to think that the strong coupling dynamics perhaps gives rise to qualitatively new physics. Furthermore, since the thermalization time is unlikely to vanish for arbitrarily high chemical potential, we expect it to either saturate or turn back. In both cases, the scaling of the thermalization time seems to group the dynamics in two qualitatively different regimes: µ/T 1 and µ/T 1, much like what was observed in [34]. Now we will turn to discussing the details of our model.

Einstein-Maxwell-Dilaton system
We begin with the following action which leads to the following equations of motion

2)
where Let us now specify our ansatz which consists of the metric field, the Maxwell field and the scalar field: {G, A, φ}. We will assume translational invariance and we choose the ingoing Eddington-Finkelstein patch to represent our ansatz data: We need two more sets of data: the boundary conditions and the initial conditions. The boundary conditions simply impose that the geometry is asymptotically locally AdS, which is represented by

13)
A v (z, v) = const + O(z) . (3.14) This choice fixes the gauge completely. 6 Let us specify the initial condition now. Our initial state in the dual field theory corresponds to a thermal state with a non-vanishing chemical potential. This is represented by Our goal now is to find a solution of the system of equations in (3.2)-(3.4) subject to the initial conditions in (3.11)-(3.13) and with the asymptotically locally AdS boundary condition in (3.15)- (3.18). We want to introduce dynamics in the system by exciting a time-dependent nonnormalizable mode for the scalar field near the boundary, which can be represented by where φ 1 (v) is now a function 7 of O(1) and the dimensionless parameter 1 will eventually serve as the expansion parameter. To connect with the discussion in section 2, we note that: = κ 1/2 quench . Note that, the energy-momentum tensor in (3.6), (3.7) evaluated on (3.19) is not of a null-dust-type, as considered in i.e. [34], and thus we will not encounter potential subtleties associated with violating null energy condition [41].
Before proceeding further, let us comment on the particular coordinate patch. To incorporate the dynamics, we have chosen the Eddington-Finkelstein coordinates, collectively denoted by U EF , which is well-defined everywhere. On the other hand, in order to read off the stress-energy tensor of the dual field theory, it is very useful to express all data in terms of the Fefferman-Graham patch, which we denote by U FG . We can define a map ϕ : U EF → U FG , such that where ≡ ∂ r and˙≡ ∂ t . Near the boundary we have Note that the symmetry and boundary behavior requirements allow the scalar φ to be an arbitrary function of time at the boundary.
After appropriate holographic renormalization, and using the GKPW recipe, the stress-energy tensor of the dual field theory is obtained to be [42,43,44] Evidently, once we obtain a solution in the U EF -patch, we can use the map ϕ : U EF → U FG by solving equations in (3.22)-(3.24) and finally using (3.25), (3.26) we can read off the field theory stress-tensor of the corresponding state. In practice though, for the initial and the final states, the boundary energy-momentum tensor can be obtained by analyzing the thermodynamics of the corresponding state.

Asymptotic, z → 0, expansion
Let us first investigate the near boundary behavior of the solution to (3.2)-(3.4) to ensure the asymptotically locally AdS criterion. As z → 0, we introduce the formal expansion

27)
A where φ b (v), as introduced in (3.19) before, denotes the source which we are turning on at the boundary and the set of functions {f n , g n , p n , a n } can be systematically determined from the equations of motion at each order in z-expansion.
For illustrative purpose, we provide below explicit formulae up to O(z 5 )-term The two time-dependent functions M (v) and L(v), which are not determined by the asymptotic expansion above, physically correspond to the mass of the black hole and the VEV of the marginal operator, respectively. Note that, the solutions in (3.28)-(3.31) represent an asymptotic solution of the full geometry. So far, we have imposed the asymptotically locally AdS condition in details. In order to carry out a perturbative treatment, we need to identify the "correct" initial state around which it is meaningful to carry out a perturbation order by order. We will answer this question a priori: it turns out that one choice for which such a perturbative treatment works is to start from an AdS-RN geometry with a "small" mass and a "small" charge. Thus, instead of arbitrary M and Q in (3.15), we rewrite them as where 1 1 is another small parameter. We can identify 1 = κ 3 µ to connect to the discussion in section 2. Clearly, and 1 are hitherto independent parameters.
Thus, the initial state can be written as

35)
A Here µ i can be fixed by demanding regularity of the gauge field at the event horizon.

Expansion in amplitude of
To solve the equations of motion, we now work with the following formal expansion where the data {f n , g n , Φ n , A vn } are to be systematically determined from the equations of motion.
In order to perturbatively solve the equations of motion and motivated by convenience, we will treat 1 as a parameter of O( ). Thus, Note that the above relation is akin to the condition O(κ µ ) = O(κ quench ). It is clear that in this regime we regroup the dynamical contributions at various unrelated orders involving and 1 .
Physically this corresponds to linearly combining processes which occur at e.g. O( 2 ) and O( 2 1 ), and hence there is only one expansion parameter at the end.
We will be able to check two independent limits of setting → 0 and 1 → 0 using (3.41). Solving the dilaton equation of motion requiresḣ(0) = 0. Similarly, from Einstein field equation we obtain h(0) = 1. To present the explicit solution up to, e.g. O( 4 ), we first define the following functions In order to properly account for the powers of , let us write the -dependence of C 4 (v), P (v) and a 4 (v) explicitly. From (3.43), (3.45) and (3.46) we have, where now neither c 4 , p nor c 4 , p, a 4 depend on . With these definitions the solution can now be written as: 8 (3.53) We can now relate the boundary quantities {M (v), L(v), c} appearing in (3.28), (3.30) to the amplitude expansion: Before going further, let us check a couple of trivial limits: (i) First, note that setting = 0 keeping 1 = 0, we kill off the entire dynamics and get back the initial AdS-RN state in (3.34)- (3.36). In other words, this limit is equivalent to taking v → −∞. (ii) Secondly, if we set 1 = 0 keeping = 0, our initial state reduces to empty AdS. It is straightforward to check that this case reduces to the one with a vanishing chemical potential [23]. However, the non-trivial dynamics remain. (iii) Third, if we set both 1 = 0 = , then we are left in the empty AdS-background with no dynamics. This is the most trivial limit of the solution described above. As alluded in (3.41), we will consider 1 = 0 and = 0, with the constraint that k = 1 / ∼ O(1).
At late times, i.e. v δt, the final state is given by Here the tilde denotes the function evaluated at any v δt. Since φ 1 is of compact support, this is equal to the value of the function at infinity, Note that to fourth order, the mass of the black hole has increased by an amount C 2 2 + c 4 4 as compared to the initial state, similarly to [23]. Exploring the solution at higher order we find that at sixth order g(z, v) has a dependence onḧ(0) that survives at late time and gives a subleading contribution to the mass of the final black hole. Note also that the gauge field at late times differs from the original gauge field only by a shift in the chemical potential. This is consistent with our expectations since the charge of the black hole remains the same. However, the position of the horizon changes and thus -as we will see in section 4 -the chemical potential also changes due to the boundary condition on the gauge field. We will analyze the resulting thermodynamics in Section 4.

Analytic structure and regime of validity
Let us briefly investigate the analytic structure and regime of validity of the perturbative solution for v > δt. We can systematically find the solution to arbitrary order in . Then, following a similar approach as in [23], we can inductively show that, among the functions that appear in (3.37)-(3.40), Thus the non-trivial information about the dynamics is contained in the set of functions: G ≡ {G n } = {g 2n , f 2n , φ 2n+1 , A vn } and they take the following general form In general G is a function of v and a functional of φ 1 (v) and its derivatives: Note that even powers of are absent in φ(z, v) and odd powers are absent in g(z, v), f (z, v) and The fact that the functionsG n are polynomials in v whose degree grows with n implies that the series expansion will break down at late times. In order to characterize the regime of validity lets focus on φ(z, v), which has the maximum degree in v: It can further be checked that for late times It would suffice for our purposes, if the perturbative solution is valid up to the event-horizon of the geometry. Recall that the event-horizon is given by, up to the leading order in , (3.66) 9 As stated previously, the equations of motion demand that h(0) =ḣ(0) = 0. The requirement that the odd ... = 0) comes from assuming a series of the form 3.62 for the dilaton and the gauge field. 10 The late time behavior of A v is different; all terms of order 4 or higher go to zero at late times. Thus, at v > δt we recover the original gauge field and A v does not enter in the analysis of the late time validity of the solution.
Also recall that k and m are O(1) numbers whileC 2 ∼ 1 δt 3 . Thus we get: z H ∼ δt 2/3 . Therefore, close to the horizon and for late times, the term with labels n, k in (3.64) will scale approximately as which to leading order in corresponds to 1/ (∆E) 1/3 (see eq. (4.33)).

Thermodynamics of the states
We will now briefly discuss the initial and the final states which are interpolated by the evolution described in (3.50)-(3.53). Both our initial and the final states are characterized by a nonvanishing temperature and a chemical potential, both of which have dimension length −1 . It is also accompanied by the VEV of the marginal operator, denoted by O φ ∼ length −3 , which is being quenched via the dynamics. Thus the state is specified by the following data: Furthermore, note that from the full solution in (3.53) it is clear that the normalizable mode of the gauge field does not have any dynamics associated and thus lim z→0 which implies that the charge density in the dual field theory is kept fixed. Thus, we are considering the dynamical evolution in a "canonical ensemble". We will now discuss the initial and the final states in more detail. The regularized on-shell Euclidean action corresponds to the Gibbs free energy, which characterizes the grand-canonical ensemble. If we denote the Gibbs free energy by W , then the Helmholtz free energy, which we denote by F , characterizes the canonical ensemble and is obtained by a Legendre transformation of the Gibbs potential where µ and Q are respectively the chemical potential and the charge density.

Initial state
Here we will reinstate 1 where they originally appear. According to (3.34)-(3.36), in the limit v → −∞, we get: With the scalar field turned off, the geometry reduces to the usual AdS-RN black hole. For the above solution, the corresponding mass M and charge Q are given by Alternatively, the mass can also be given as On the other hand, the one-form A v must be regular at the horizon such that ||A|| remains finite at the bifurcation point of the Kruskal-extended patch. This imposes a constraint, relating the chemical potential to the charge and the mass: The subscript i in all subsequent physical quantities will stand for the initial state. To calculate the Hawking temperature of the initial state T i , we first perform a Wick rotation obtained as usual by the replacement t → iτ . Since the Euclidean time direction shrinks to zero size at z = z H , we must require that τ be periodically identified with appropriate period β i , i.e. τ ∼ τ + β i . A simple calculation shows that 11 For concreteness, we will evaluate all physical quantities in two leading order terms in 1 . or equivalently, It is convenient to invert the relations µ i (m, q) and β i (m, q) given in (4.11)-(4.13) to obtain m(µ i , β i ) and q(µ i , β i ). However, we have to proceed with some care given that µ i ∼ O( 1 β i such that they both are of order O ( 0 1 ). Then, an expansion in 1 is reliable and the inversions can be found perturbatively. To our order of approximation we find, and To study the thermodynamics of these solutions, we first evaluate the Euclidean action I onshell which defines the grand canonical (Gibbs) potential W = I/β i (see for instance [61]). With our conventions, the full Euclidean action is given by analytically continuing (3.1). Moreover, when the space is asymptotically AdS the Gibbons-Hawking boundary term gives a vanishing contribution.
As usual, the action diverges upon integration, given that the volume of any asymptotically AdS geometry goes to infinity near the boundary. These divergences can be eliminated by subtracting the pure AdS contribution, obtaining This may be rewritten entirely in terms of β i and µ i as The grand canonical potential is given by , (4.20) Together, they indeed satisfy the first law of thermodynamics, dE i = T i dS i + µ i dQ i . Moreover, the free energy W = I/β i is always negative, indicating stability of the solutions. Notice that we could have alternatively worked in the canonical ensemble, which is characterized by the Helmholtz free energy F = E − ST . A brief computation shows that, at the same order of approximation, Let us now comment on the identification of the field theory quantity corresponding to the putative small parameter 1 . From (4.12), (4.11) and (4.20), it is clear that parametrically Thus, we are considering an initial state in which the chemical potential is small compared to the temperature. It also becomes clear that, in obtaining the perturbative solution outlined in . By analyzing the final state, we will now observe that corresponds to an otherwise independent parameter in the dual field theory.

Final state
Let us now consider the final state. For v → ∞ (in practice for v δt) the bulk solution in (3.57)-(3.60) can be written as:

26)
A v (z) = µ f + 2 k 2 qz , (4.27) where k = ( 1 / ) and all the tildes denote the corresponding functions evaluated at v → ∞, which are of order O ( 0 ). The mass and charge of the final background are now which gives The new horizon lies at For convenience, we have defined m 2 = k 2 m + C 2 . Regularity of the gauge field at the horizon now yields On the other hand, the inverse temperature is given by After subtracting the AdS contribution, the Euclidean on-shell action evaluates to: Finally, we can compute the state variables for the final state. At leading order we get: Again, these quantities satisfy the first law of thermodynamics, dE f = T f dS f + µ f dQ f . Note that, upon integration by parts, we get that C 2 = 1 2 ∞ −∞ φ 1 (x) 2 dx > 0 so we always have m 2 > m. Both the energy and entropy increase but the charge remains the same, as expected. Moreover, for this final state, the Helmholtz free energy is given by Before concluding this section, let us comment on the parameter . It is straightforward to check that Note that, this scaling is in keeping with the scaling behaviour obtained in [45,46,47], which follows from dimensional analysis in our case. Evidently, (4.39) provides us with a natural meaning for the expansion parameter in gravity purely in terms of the field theory data. What is more, we observe that the perturbative solution is consistent as long as we impose

Thermalization time
Given the background that we obtained in the previous section, we will now explore how the thermalization time behaves with relevant parameters for the system. We will measure thermalization time by measuring non-local observables, specially entanglement entropy and define our thermalization time according to the behavior of this non-local probe. Via the AdS/CFT correspondence, this exercise amounts to computing an extremal area surface in the bulk geometry. For this purpose, we will only require the solution for the metric (3.8), which is given explicitly in (3.57)-(3.60). The specific form of f (z, v) and g(z, v) will not be important for now. The metric (3.8) is asymptotically AdS. To make this manifest, let us rescale the metric functions such that where we have definedf (z, v) ≡ zf (z, v) andg(z, v) ≡ z 2 g(z, v). These new functions satisfy thatf → 1 andg → 1 as z → 0. Also, the boundary time coordinate is related to the Eddington-Finkelstein coordinates in (5.1) through 13 dv = dt − dz g(z, v) .
This relation will be important below, for the computation of the thermalization time.

Entanglement entropy
We will compute entanglement entropy in this background for a particular shape: namely a "rectangular strip" which can be parametrized by {x ∈ (− /2, /2)} ∪ {y ∈ (− ⊥ /2, ⊥ /2)}. In a quantum system, entanglement (or geometric) entropy of a region A with its complement B is defined using the reduced density matrix a la von Neuman: where ρ A is obtained by tracing over the degrees of freedom in B, ρ A = tr B ρ.
In AdS d+1 /CFT d , the covariant prescription to compute entanglement entropy for a given region A was proposed by Hubeny, Rangamani and Takayanagi in [48], generalizing the originally proposed Ryu-Takayanagi formula [49]: where G N is the bulk Newton's constant, and γ A is the (d − 1)-dimensional extremal-area surface such that ∂γ A = ∂A. 14 Before delving into computation, let us flesh out the approximation in which we will be working for the rest of this article. The perturbative solution that we have obtained in (3.50)-(3.53) can be numerically subtle to handle, specially since we have a small parameter . Instead, we make use of the small parameter in the following manner: The metric data can be summarized as follows: Subsequently the geodesics and the entanglement entropy can also be determined by a similar expansion A + . . . , (5.5) A + . . . . (5.6) 13 Here we neglect the redshift effect indicated in [35]. In the thin-shell limit for v 0 → 0, the redshift factor will reduce to one for the regime outside the shell. Note that the t crit introduced in (5.30) will be slightly longer when incorporating the redshift effect. 14 We note that at this point the evidence in favor of the Ryu-Takayanagi proposal is much stronger than the covariant HRT proposal. However, till this issue is settled, we are assuming that HRT proposal is the correct prescription in our case.
It can be checked that upon using the equations of motion at the first non-trivial order in , which in this case happens to be at O( 2 ), γ (0) A completely [50,51]. Thus, if we limit ourselves to results at the first non-trivial order in , the task of determining geodesics simplifies significantly. On the other hand, it is clear from the metric in e.g. (3.50), the effect of the charge enters at O( 4 ), at which order there is no such simplification. In this article, for simplicity, we will keep our discussions limited to O( 2 ). We therefore emphasize that, our numerical results should be taken as indications of a certain physics rather than an observation. On the other hand, it is possible to go beyond O( 2 ) systematically, which we leave for future work.
Also note that, as we are measuring the response of the system at the first non-trivial order in , so we will measure the external dial. In our case, the external dial is a combination of e.g. (µ/T ) ∼ O( 2/3 ). Thus, effectively we are using a different precision for measuring the response-observables as compared to the parameters of the system. In a very precise sense, the subsequent dynamics that we analyze is primarily determined by the temperature-scale of the system. Now let us discuss the details. We will now compute S A in the limit ⊥ → ∞, so that the construction becomes invariant under translations in y. We can parameterize the extremal surfaces with functions z(x), v(x) subject to the boundary conditions z(± /2) = z 0 and v(± /2) = t , where z 0 is the usual UV cut-off needed to regularize the on-shell acion. These boundary conditions impose that the boundary of γ A coincides with the boundary of A along the temporal evolution. The area functional is given by where ≡ d/dx. Since there is no explicit x-dependence, there is a corresponding conservation equation given byf In this expression, z * is defined through z(0) = z * . Now, the right hand side of (5.2) in empty AdS is identically zero. In our case it is of order O( 4 ), which we can safely neglect since we are working to order O( 2 ). Thus, we have Combining (5.9) and (5.10) we get In practice we solve the above system as follows. First, we rewrite (5.10) aṡ v + 1 g(z, v) = 0 , (5.12) where˙≡ d/dz and we solve for v(z) subject to the boundary condition: v(z 0 ) = t boundary . With this function at hand and given a value of z * , we can then obtain a solution for z(x) by direct integration of (5.11). However, note that this last step is not necessary if we only want to extract the values of and A for a given z * . More specifically, these values can be obtained from respectively. In particular, the last equation in (5.13) arises upon substituting (5.10)-(5.11) in (5.8) and then changing the integration variable dx → dz/z . On the other hand, the area in (5.13) is divergent as the UV cut-off z 0 → 0 and must be regularized. The divergence comes from the fact that the volume of any asymptotically AdS background is infinite and the spatial surface γ A reaches the boundary. In particular, near the boundary it is clear thatf (z → 0) → 1, z (z → 0) = z 2 * /z 2 and therefore Subtracting this divergence, we obtain the finite term of the area which is the quantity we are interested in:

Toy model and choice of parameters
Before proceeding to the numerical results, we must specify the system we want to study. First, notice that at O( 4 ) the computation of entanglement entropy does not rely on the specific form of the coupling h(φ). Then, for the purposes of this section it is enough to set h(0) = 1 anḋ h(0) = 0, as required by the perturbative solution. For the scalar profile, we choose for simplicity a Gaussian function of the form where λ and v 0 are numerical parameters that control the amplitude and the width. With this choice, the leading-order amplitude of the dilaton scales as φ ∼ λ. Note that the perturbative expansion requires each order in to be at most of order O(1) in . Here we have introduced the extra parameter λ for convenience, which could yield (C 2 m ) ∼ O(1) with a proper choice of its numerical value. We will come back to this point below. One advantage of the above profile is that it can be integrated analytically to obtain the explicit form of C 2 (v), C 4 (v) and P (v). A brief computation leads to: Here erf(x) denotes the error function. Furthermore, we only need derivatives of φ 1 (v) and P (v) to compute the analytic forms off (z, v) andg(z, v) according to (3.57)-(3.60). Next, we have to select values for the various parameters in order to be consistent with the perturbative expansion. First, we set m = 1 and v 0 = 0.01. The first choice fixes a reference scale for the initial energy whereas the second guarantees that we are dealing with the thin shell limit. Then, we must fix λ in order to have m 2 ∼ O(1) for the final state. Note thatC 2 evaluates toC Then, a good choice for λ is so thatC 2 = 1. We further set k = 1 so = 1 and m 2 = k 2 m +C 2 = 2. In Figure 1 we plot both the gaussian profile for the scalar field and the mass function m 2 (v) ≡ k 2 m +C 2 (v) for the parameters given above. The remaining functions C 4 (v) and P (v) are of order O (10 −7 ) and O (10 −4 ), respectively. Note that the mass function does not increase monotonically in time, in stark contrast with the usual behavior of collapsing Vaidya geometries. Nevertheless, our field content is physically sensible and all the energy conditions are satisfied. Following the arguments of [41], then, we expect reasonable results for the thermalization process in the boundary theory. Another reason to argue that this must be true is that, although the apparent horizon (2.15) also shows this non-monotonic behavior, the event horizon (2.17), on the other hand, always increases along the temporal evolution. 16 Figure 2 shows representative behaviors of these two quantities. At any rate, it is worth pointing out that such non-monotonic behavior disappears as we take the thin shell limit, which is the case we are focusing on. In this limit both m 2 (v) and z aH (v) take the form of a step function. The expansion parameter should be a small number. We observe that the first and second corrections to the functionsf (z, v) andg(z, v) over AdS are of order O( 2 ) and O( 4 ), respectively. A reasonable choice for is then = 0.1. Finally, the value of q can be tuned in order to vary the temperature and the chemical potential of the solutions. There are two requirements for this quantity: (i) we must impose that q ∼ O(1) to be consistent with the perturbative expansion and (ii) we cannot exceed the maximum value allowed in order to avoid a naked singularity at early times. Regarding this last point, recall that to our order of approximation, the initial state is characterized by Therefore, from (5.22) it follows that the condition (ii) for not having a naked singularity sets a 16 If we truncate the metric at order O( 4 ), we find that for v 0 ∼ 10 or larger the event horizon also presents signs of non-monotonic behavior. This issue is corrected as we consider higher order corrections in . above which T i < 0 and with no real roots of g(z H ) = 0. Fortunately, notice that (5.23) is also within the acceptable range for satisfying item (i). The final state, on the other hand, is characterized by where It is easy to check that for values of q in the range allowed by (5.23), the final states are also free of naked singularities. Of course, this is directly related to the fact that the mass of the black hole is increased while the charge is kept fixed. Now, we have both temperature and chemical potential in our system and we need to construct a dimensionless ratio which would be our tunable parameter. Notice that both temperature and chemical potential has inverse length dimensions. Then, we can consider [34] χ ≡ 1 4π to be the relevant parameter that we will vary. In practice, we can vary χ by tuning the parameter q ∈ [−8.04, 8.04]. In particular, it is found that for this range of values, χ(q) ∈ [−0.42, 0.42] and increases almost linearly with q. Thus, within the perturbative approximation and for the choices of the various parameters, we are restricted to small values of χ as compared to [34].

Regimes of thermalization
Now we will discuss the numerical results. In Figure 3 (a) we plot sample solutions for the embedding functions z(x) for a fixed as the boundary times t is varied. Some of them cross the shell, located at v = 0, and refract. This refraction is suppressed by a factor of 2 given that the energy-momentum of the shell is itself of order O( 2 ). Nevertheless, the aforementioned effect is noticeable to the naked eye for large enough distances, (4π∆T ) 1. We also show in (b) the behavior of entanglement entropy as a function of time, for a fixed distance . In order to compare the results for various values of χ we subtract the entropy of the initial state ∆S(t) = S(t) − S(−∞), and focus on the entanglement growth over time. 17 We note some general properties of the behavior of ∆S(t) as we change the value of χ. Qualitatively, our results agree with those of [52,53] obtained in the context of Vaidya geometries. At early times, i.e. the "pre-local-equilibration" regime, the evolution is almost quadratic in time and is weakly dependent on the size . This stage is followed by a linear growth phase at intermediate times and, finally, a saturation is reached at late times, t ≥ t sat , where the entropy abruptly flattens out.
Let us explain these regimes in more details. First, in Figure 4, we have shown the functional dependence of entanglement entropy growth in various regimes. Following [52] we can introduce a "local equilibrium scale", denoted by eq , which means that a local thermodynamic description applies at length scales∼ O( eq ) even though the global description is out-of-equilibirum. For early times, which can be represented by the regime t eq , the rate of entanglement entropy growth is expected to be proportional to the area of the entangling surface. Furthermore, we can assume that the growth is proportional to a characteristic energy scale of the system. Note that, at the initial stage, we have the energy of the initial state, denoted by E initial , and the energy which is being injected to induce the dynamics, denoted by ∆E. Thus the typical energy-scale should be identified as A priori, E initial and ∆E are independent. However, in our case, we already demanded E initial ∼ ∆E and thus, using dimensional analysis, we can arrive at where α 1,2 are two constants and A denotes the area of the entangling surface.
On the other hand, in the regime t eq , we have a notion of a thermal entropy density which is denoted by S thermal . If we further assume that the entanglement entropy is proportional to the area of the entangling surface and the local thermal entropy density, then by dimensional analysis we get where v E is the entanglement production rate which has been analyzed in [52,53]. Evidently, our numerical data agrees very well with the intuition outlined in (5.28) and (5.29). It is intriguing to further note that, although we are not starting from a vacuum state, the analysis of [52,53] continues to hold. 18 Finally, since we are considering the "rectangular" shaped entangling surface, we expect that the saturation will be accompanied by an abrupt jump in the corresponding extremal area geodesic. Now we will discuss the scaling of thermalization time. From the time-evolution of entanglement entropy, we can extract t sat for a given length . Alternatively, we can define another time-scale t crit as a measure of the thermalization time. Recall that the shell is densely peaked around v = 0 for v 0 1. Thus, we can define t crit to be the time at which the corresponding extremal surface grazes the shell at v = 0. By definition from (5.2), we get In practice, it is easier to extract t crit rather than t sat . 19 Furthermore, in the limit v 0 → 0 these two quantities are found to agree. For the case of finite thickness, they only differ by a factor of order O(v 0 ) so we expect similar results for the thermalization time as long as v 0 is not too large. We thus focus on the critical time t crit . Let us, however, clarify a further caveat regarding t crit . Precisely speaking, t crit measure the 19 Note, however, the approach to equilibrium is expected to be abrupt in the case of a rectangular entangling region [52]. Thus, strictly speaking, t crit may not be the correct measure of thermalization time. However, in our case, we have checked that t sat and t crit exhibit qualitatively similar behaviour. Since we are only concerned with qualitative features, we do not attempt to make this more precise here. time a null ray takes to reach the AdS-boundary starting from the bulk point z * . Thus, for a given boundary length , the corresponding extremal area geodesic will indeed graze the shell at t crit provided the shell propagates at the speed of light. However, as can be checked from (3.6) and (3.7), the matter field is non-null and thus propagates slower than the speed of light. Thus, in reality t crit will serve as a lower bound for the actual thermalization time.
The dependence of t crit with is shown in Figure 5 (a), which approximates a linear growth similar to the one observed in [30,31]. Generally, we can represent the dependence as where A(χ) represents the slope of the linear regime, i.e. the velocity at which thermalization propagates in the system. Numerically we find A(χ) = 0.83 for blue (5.32) = 0.79 for red , (5.33) which implies that thermalization is super-critical (faster than speed of light). At larger lengthscales, the deviation from linearity is characterized by the second term above, where α(χ) is an index which can, in principle, depend on χ.
On the other hand t crit monotonically decreases with χ for fixed (in the allowed range for χ), which means that thermalization is faster in the presence of a chemical potential. We have shown a representative behaviors in Figure 5 (b). Our findings agree qualitatively with the results of [34], in the regime of validity of our solutions.
Let us now comment on what may happen beyond the validity of our perturbative solutions. It is unlikely that the thermalization time keep decreasing with increasing chemical potential. Hence there are two possibilities: (i) thermalization time plateaus or (ii) thermalization time eventually turns around and starts increasing with increasing chemical potential. Either way, it implies a qualitatively different scaling behaviour of thermalization time in two regimes: when χ 1 and when χ 1, i.e. which has an obvious interpretation as a "classical" and a "quantum" regime respectively. The latter observation is similar to the one made in [34].

A remark on the scrambling time
Before closing this section, we wish to make a few comments. 20 Our initial state is thermal, and thus is represented in the gravity description by an AdS black hole. It is generally believed that black holes are endowed with the special property of "fast scrambling" [55,56]. In particular, this implies that the time scale associated to the process of thermalization grows logarithmically with the number of degrees of freedom of the system where β is independent of N c . In our case, we can count N c by evaluating the thermal entropy of the initial state S i Thermal , or of the final state S f Thermal , or of their difference ∆S Thermal . Now, in order to make contact with the scrambling property of the black hole, we have to focus on the regime ∆T 1 where entanglement entropy approaches the thermal entropy, S(∆T 1) ∼ S Thermal . In this same limit, we can identify t crit ∼ t * since t crit (∆T 1) serves as a measure of "global" equilibration. In Figure 6 we show the behaviour of log ∆S as a function of the ∆T t crit for a fixed value of chemical potential χ. To generate the plot, we vary the length to generate pair of points {∆S( ), ∆T t crit ( )} and then we join them. Larger values of correspond to larger t crit ( ) for a fixed ∆T (see Figure 5), so we are interested in the rightmost part of Figure 6. It is worth pointing out that, due to limitations of the numerical accuracy and the validity of our perturbative solution, we can only go as far as ∆T ∼ 1. This is because in order to increase the length of the boundary strip we have to reduce z H − z * accordingly, which needs to be fine-tuned to high precision in order to increase ∆T . However, it is intriguing that in the region around ∆T ∼ 1 the curve smoothly approaches a straight line, in concordance with (5.34). It would be remarkable if this statement holds true as ∆T is increased. If it does, this may lead to important clues on how the effective degrees of freedom interact towards the process of thermalization [57,58]. It will be interesting to investigate this issue further and make contact with other approaches within the framework of AdS/CFT (see, e.g. [59,60]).

Evolution of L(v): a local observable
In section 2, we remarked that we can observe a non-trivial evolution of local observables in our case. A natural local observable is the vacuum expectation value (VEV) of the marginal operator which we are quenching here. To extract the precise VEV, we need to carefully perform holographic renormalization of our perturbative solution in (3.50)-(3.53). However, here we merely want to illustrate this time-evolution, and the entire machinery of holographic renormalization is redundant for our purposes. Thus, in figure 7 we show the time-evolution of L(v), which is given in equation (3.55). 21 Note that, the relatively small magnitude of the vertical axis is a consequence of the choice of the parameters that we have previously made. It is clear that the evolution occurs at a time-scale of O(δt) = O(v 0 ). First, we can generalize our perturbative analysis when the quench is being carried out for a relevant operator. This will correspond to introducing a non-trivial scalar potential and it will be interesting to check to what extent our observations are universal. In turn, it is easier to embed such effective gravity descriptions within 11-dimensional supergravity. Subsequently it will also be interesting to understand the specific embedding of our effective gravity description, in which the scalar potential is vanishing, in string theory. Note that, our model is unlikely to be realized within ABJM, since the latter does not have any marginal scalar-deformation. 22 To test the robustness of our qualitative observation, it will be interesting to model a similar dynamical background in the presence of more than one global U(1)-charge, by e.g. coupling the STU-model with a neutral massless scalar. 23 The quench of this scalar field will correspond to the quench of a marginal or relevant operator in the dual field theory. It will be very interesting to check how thermalization behaves in this case.
We have considered only (3 + 1)-dimensional bulk theory which corresponds to a (2 + 1)dimensional boundary theory. It is well-known that the dynamics in asymptotically locally AdS d+1 differs qualitatively for even vs. odd d, see e.g. [23]. Thus generalizing our results for the asymptotically AdS 5 -background will be an interesting avenue to pursue.
Intuitively a generic charged-system is much richer than the neutral case and one can consider various possibilities. In this article, we considered the quench by a neutral scalar. One can also consider a charged-scalar field. In the latter, depending on the temperature of the system, the ground state of the system can be described by a Reissner-Nordstorm black hole or a hairy black hole, corresponding to a "normal" phase and a superconducting phase of the field theory. Clearly, this corresponds to a more complicated dynamical process and it will be worthwhile to check how much physics can be accessed via a perturbative approach analogous to what we have adopted here.
Generalizing on this theme, it will be interesting to catalogue the possibilities in which Einstein equations admit such perturbative analyses. A bottom-up approach can be followed with various matter content in an Einstein-gravity theory with a negative cosmological constant. We can further relax the asymptotically AdS-condition, by using an asymptotically Lifshitz or Hyperscaling-Violating geometries.
The importance of a complete numerical evolution can hardly be overestimated. It is very crucial that eventually we have access to the entire evolution process in the full parameter space. As outlined above, there are numerous physically inequivalent real-time phenomena that can be captured within Einstein gravity with Maxwell and (charged or uncharged) scalar fields in various dimensions, which needs extensive numerical explorations. This is a long-term goal which we leave for future. hospitality and support, where significant parts of this work were done.

A Series solution up to 8th order in amplitude
In section 3.2 we presented the perturbative solution up to order 4 . In this appendix we will present the explicit form of the solution at higher orders of , where specially the non-trivial contributions coming from the coupling of the neutral scalar with the gauge field become transparent.
Recall that in (3.62) we gave the formal structure of the functions entering the expansion (3.37) -(3.40), The solution up to order 4 involves the functions C 2 (v), P (v), C 4 (v) and a 4 (v) whose explicit form was given in (3.42) - (3.46). Since these functions also enter the higher order expansion we quote them here again for sake of completness.
. We now proceed to write the solution to order 5 , 6 , 7 and 8 .