Holographic Thermalization with Weyl Corrections

We consider holographic thermalization in the presence of a Weyl correction in five dimensional AdS space. We first obtain the Weyl corrected black brane solution perturbatively, up to first order in the coupling. The corresponding AdS-Vaidya like solution is then constructed. This is then used to numerically analyze the time dependence of the two point correlation functions and the expectation values of rectangular Wilson loops in the boundary field theory, and we discuss how the Weyl correction can modify the thermalization time scales in the dual field theory. In this context, the subtle interplay between the Weyl coupling constant and the chemical potential is studied in detail.


Introduction
The AdS/CFT correspondence [1], [2], [3] or the gauge/gravity duality is one of the most striking aspects of string theory and is being extensively used over the past few years to study strongly coupled condensed matter systems. The duality relates a classical theory of gravity in a (d + 1) dimensional anti-de Sitter (AdS) spacetime to a strongly coupled conformal field theory living on the boundary of the AdS space in d spacetime dimensions. The strongly coupled nature of the boundary theory does not allow the usage of standard perturbative techniques. However, using the gauge/gravity duality, the computation becomes simpler to handle, because the dual gravity theory is classical. This is the primary motivation to explore phenomena in strongly coupled quantum systems from a holographic point of view. It is probably fair to say that by now, a clear understanding of the near-equilibrium physics of strongly coupled quantum field theories arising from the dual gravity sector has emerged. For example, one can calculate the bulk correlation functions [4] and compute different observables holographically and understand the linear response of the system to perturbations from equilibrium [5], [6]. However, it is quite difficult to understand the physics of a strongly coupled system, even from the dual gravity sector, when the system is out of equilibrium, because now one can not apply linear response theory. It is in fact very interesting to analyze how such a system reaches thermal equilibrium, once it is out of equilibrium, and calculate the "thermalization time." In a class of examples, this issue has been resolved holographically by constructing a time-dependent gravity solution in AdS space which describes the formation of a black hole at late times.
The other fact that motivates the study of non-equilibrium dynamics from the gauge/gravity duality, is the experimental input from the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC). When two large energy heavy ions collide in RHIC, some of their kinetic energy is transformed into heat energy. Because of the large amount of heat produced, the quarks and the gluons form a plasma-like state, known as quark gluon plasma (QGP). The process of forming QGP is known as thermalization. After the QGP is formed, i.e., after the thermalization process is over, it reaches the thermal equilibrium where one can apply linear response theory to understand the physics of QGP. The experimental result shows that the QGP formed in RHIC behaves like an ideal fluid, with a very small shear viscosity to entropy density ratio (η/s) indicating the strong coupling nature of the QGP [7]. Because of the strong coupling constant it is appropriate to use the AdS/CFT correspondence to compute η/s, and check whether they match with the experimental results. A lot of work [5], [8], [9], [10], [11] has been done in this direction concluding that there exists a small lower bound of the ratio which may depend on the coupling constant of the theory. While it is easy to compute different observables after the QGP is formed, it is difficult to compute the thermalization time, since, as mentioned before, thermalization is a non-equilibrium process. Also the observed thermalization time in RHIC is much shorter than what is predicted by calculations via perturbative techniques [12]. This indicates that the thermalization process also takes place within a strong coupling regime of QCD. These experimental results provide further important motivation to analyze thermalization, using the gauge/gravity duality. We also point out that there are large number of articles [13]- [22] where the thermalization process has been described as the dual of a black hole formation by a gravitational collapse in the bulk.
In [23] and [24], the authors considered an interesting model for thermalization in the context of AdS/CFT, and examined the thermalization after a sudden injection of energy to the boundary field theory by three different nonlocal probes : two-point correlator, Wilson loop and entanglement entropy. All of these three probes are well-defined in the dual gravity theory, and they are described by different geometric quantities. For example, the two-point correlator on the boundary corresponds to a geodesic connecting the two points and extending into the bulk. Similarly, the expectation value of the Wilson-loop operator and the entanglement entropy correspond to a minimal area surface and minimal volume, respectively, extending into the bulk. The model they considered is known as the AdS-Vaidya metric which describes the collapse of a thin shell of matter from the boundary to the bulk. As the shell collapses, it divides the spacetime into two region: the outer region of the shell represents a black brane while the inner region corresponds to pure AdS spacetime. Hence, at the early time, the AdS-Vaidya metric corresponds to a pure AdS space (representing a vacuum state of the boundary QFT), while it represents a black brane metric (representing a thermal state of the boundary QFT) at late times after the shell collapses.
For all kinds of probes, it was shown in [23], [24] that the UV degrees of freedom thermalize first and the IR degrees of freedom thermalize later, i.e., the thermalization process is top-down. 1 While a standard perturbation technique in QCD predicts that the thermalization process should be bottom-up [26], these papers get an opposite behavior. Note that the bottom-up behavior of thermalization in heavy-ion-collisions in a perturbative QCD can be realized as follows : when the thermalization initiates, a large number of soft gluons are emitted because of the large collision energy. These soft gluons collide amongst themselves, and equilibrate quickly to form a thermal bath. Hence, it is the low-energy modes or, the IR modes which thermalize first. Then the thermal bath absorbs energy from the hard gluons and when the hard gluons lose all their energy, the whole system thermalizes. However, top-down holographic thermalization is sensible from the dual gravity perspective. Simply put, since the IR modes probe more deeply into the bulk, the corresponding thermalization time should also be larger. It was also found that the thermalization time scales with the length of the probe l as τ ∼ l/2.
The work of [27] and [28] extended this model to study the thermalization in the presence of a chemical potential. These authors modelled the dual gravity theory in such a way that at late time when the shell collapses, the corresponding AdS-Vaidya metric would represent a Reissner-Nordström AdS black brane. They found that for larger probes, as one increases the ratio of the chemical potential to the temperature, the thermalization time increases. Then the idea was generalized to investigate the non-trivial corrections in the thermalization time due to the consideration of the higher derivative gravity [29], [30], Born-Infeld electrodynamics [31], Lifshitz and hyperscaling violating geometries [32], [33]. In a previous article [34], the authors studied the time evolution of holographic entanglement entropy in thermal and electromagnetic quenches using similar kind of model. In [35], the author has studied the late-time behavior of different nonlocal observables in an expanding boost-invariant plasma and shown how the fluid parameters (e.g., shear viscosity) affect the relaxation of these observables. In a recent paper [36], the time evolution of holographic n-partite information has also been investigated.
In a phenomenological i.e., bottom-up approach, it is of substantial interest to understand the process of thermalization in strongly coupled field theories upon the inclusion of general four derivative terms in the dual gravity, apart from the leading Maxwell term. Such terms are known to give rise to interesting effects -for example they non-trivially affect the η/s ratio [37]. In [38], the authors introduced these class of terms in the effective action and showed that they can lead to a violation of causality which can be prevented by the possibility of pair production. In a top-down approach, such terms are expected to arise in a string theory as quantum corrections to the low energy effective action. These terms are expected to be suppressed in a perturbative sense, and on the CFT side should represent terms suppressed by inverse powers of the 'tHooft coupling. For example, in the context of five dimensional AdS theories, one can generically think of adding all possible four derivative interactions to a usual Einstein-Maxwell action. As explained in [37], and as we review in the next section, there are a large number of such terms, but the action can be considerably simplified by choosing particular linear combination of the coupling constants. Understanding thermalization in a field theory dual to five dimensional AdS with a generic four derivative action seems a daunting task. In particular, the presence of five different coefficients that appear in such a theory is likely to make a general scan of the parameter space tedious, since the physics there is likely to depend on (fine tuned) values of the coefficients. In this paper, we consider one particular simplified situation. We consider the two derivative Einstein-Maxwell action corrected by a Weyl coupling.
In the context of holographic thermalization, one of the main motivations for considering such a coupling is the fact that it introduces an extra control parameter in the theory which might non-trivially affect the thermalization process. As we have said, a theory with all possible higher derivative couplings might be complicated, and in this paper we will see that a Weyl corrected theory already indicates non-trivial effects, which points to features that might be valid for such a generic theory. Indeed, this type of correction has previously appeared in [39], [40] who argued that as far as charge transport is concerned, starting with a general four derivative term, it is enough to consider only a linear combination of those terms, which involves a coupling of the Maxwell field to the bulk Weyl tensor. They computed the correction in the conductivity and the diffusion constant due to the Weyl coupling constant γ and predicted a bound in γ from the physical consistency conditions. This kind of four derivative interaction terms were also encountered before in [41] and [42]. In fact, QED in a general curved background leads to the Weyl coupling term at 1-loop [43].
In this paper, we will consider a Weyl correction term along with the two derivative Maxwell term in a five dimensional bulk AdS. Our purpose here would be to treat this model phenomenologically and compute the effects on thermalization due to the Weyl correction. For this purpose, we construct an AdS-Vaidya metric which interpolates between a pure AdS space at early time and a black brane solution with Weyl corrections at late time. In contrast to the previous works, the AdS-Vaidya spacetime we consider here will not be modelled by a collapsing thin shell of pressureless null dust. Rather, it is more appropriate to be modelled by a collapsing thin shell of charged fluid with some non-zero pressure.
Here, we will compute the two-point correlation function and the expectation value of the Wilson loop operator on the boundary field theory by probing the bulk with the geodesic and minimal area surfaces respectively. We find that the thermalization is always top-down and for a fixed value of the characteristic probe length l, the thermalization time decreases as one increases the value of γ, i.e., the QGP with a higher value of the Weyl coupling would thermalize faster. Further, we elaborate upon several interesting properties of the theory with a Weyl correction, for example the appearance of a swallow-tail pattern in the thermalization curve that can be controlled by γ. The paper is organized as follows: In section 2 we construct the black brane solution in linear order in the Weyl coupling constant. In section 3 we start our study of the holographic thermalization and set up the corresponding AdS-Vaidya solution by modelling the dynamical gravity with a thin shell of charged matter. In section 4 we discuss in detail about the two non-local observables we would probe: the two point correlation function and the Wilson loop. In section 5 we give a detailed description of the numerical procedure and explain the effect of the Weyl coupling constant and the chemical potential on the thermalization. In section 6 we summarize our main results. Three appendices at the end of the paper provide material supplementary to the main text.

Black Brane Solution with Weyl Corrections
In this section, we first write down the model action and then construct the black brane solution solving the Einstein and Maxwell equations. We consider the following action where a five-dimensional gravity with a negative cosmological constant is coupled to a U(1) gauge field A by the following two and four-derivative interactions : where F = dA is the usual Faraday 2-form. The Maxwell term F µν F µν represents the familiar two derivative interaction whereas the coefficients c 1 , c 2 and c 3 represent the coupling constants for the four derivative interaction terms which couple two derivatives of the gauge field to the spacetime curvature. L symbolizes the AdS length which is related to the cosmological constant Λ by Λ = − 6 L 2 . We will work in a unit where 16πG 5 = 1, G 5 being the five dimensional gravitational constant. The factor L 2 in the four-derivative interaction terms is brought in to make the coefficients c 1 , c 2 and c 3 dimensionless.
The action above is phenomenological in nature and let us briefly elaborate on this. As pointed out in the introduction, the work of [37] started with the most general four derivative action of gravity with a single U(1) gauge field. It was shown in that paper that by choosing proper field redefinition, the action can in fact be written in terms of a fewer number of terms. The final action in that paper contained five four-derivative interaction terms proportional to R µνρλ R µνρλ , R µνρλ F µν F ρλ , (F 2 ) 2 , F 4 and ǫ µνρλσ A µ R νραβ R λσ αβ along with the standard two derivative terms for five dimensional gravity along with a Chern-Simons term, where, In [44], an alternative field redefinition was used, which in turn retained the R 2 and RF 2 in the action. One can, in principle, understand thermalization with all five generic terms turned on in the action, but this will be complicated, especially since we expect the physics to depend strongly on relative values of the coefficients. We will rather focus on one particular type of correction. Our approach here is to write an action involving the coupling of the curvature to the gauge field. As pointed out earlier, such terms are sufficient to study charge transport properties of the dual CFT. We treat this as a phenomenological model to study thermalization in the dual field theory.
We consider a linear combination [39], [40], [45] of the three four-derivative interaction terms of (1) to express it into a simple form : where the five dimensional Weyl tensor C µνρλ is given by Here, γ represents the effective coupling for these higher derivative interaction terms. We will refer γ as the 'Weyl coupling' throughout the text. We will make a couple of remarks on the Weyl coupling here. In a five dimensional bulk AdS spacetime, it was shown by [39] that γ is bounded, namely − L 2 16 ≤ γ ≤ L 2 24 . While the lower bound arises to avoid the the possibility of superluminal propagation in the CFT by metastable quasi-particles, the upper bound appears to avoid the creation of certain ghost-like modes near the horizon. Hence in the probe limit, both signs of γ are feasible. Unfortunately, in the present analysis which includes backreaction, we cannot establish such a rigorous bound, as we will be working in a perturbative approximation up to first order in γ (as we elaborate upon shortly). However, in the spirit of [39], we will consider both signs of γ in our numerical analysis.
Before deriving the equations of motion one should note that the term C µνρλ F µν F ρλ can be written in the following form [37], Using the above form and making use of the Palatini identities (see Appendix A) we can construct the Einstein equation, where T µν represents the energy-momentum tensor and has the following expression, On the other hand, it is straightforward to write down the Maxwell equation, Now, taking into account the backreaction of the U(1) gauge field on the spacetime, we wish to solve the above equations (5) and (7) and try to construct a planar black brane solution. We consider the following ansatz for the metric 2 and the gauge field, A = (φ(r), 0, 0, 0, 0) .
If we take into account the backreaction of the U(1) gauge field, it is difficult to obtain an exact analytical expression for the metric that is a solution to the Einstein and Maxwell equations. Hence, we will try to perturbatively solve those equations up to linear order in γ. As we have mentioned, in a string theory, four derivative interactions are expected to arise as quantum corrections to a two derivative action, and our assumption is thus reasonable. There is however a caveat here which we should elaborate upon. We found that obtaining a controlled perturbative expansion in γ in the context of holographic thermalization is difficult. In our numerical analysis, choosing appropriate small values of γ might therefore seem a little arbitrary. Currently, we do not have a clear answer to this question. However, we have checked that the essential qualitative features of our analysis remains unaltered if numerical values of γ are chosen to be smaller than those considered in this paper. We thus expect our results to capture the essential physics up to first order in γ, with the chosen numerical values of the Weyl coupling. This is a drawback of our analysis and we will comment more on this towards the end of the paper.
In order to obtain the black brane metric, we will follow the by now standard procedure in the literature. We consider the following forms for f (r), χ(r) and φ(r) where f 0 (r), χ 0 (r) and φ 0 (r) are the zeroth order solutions representing a Reissner-Nordström AdS black brane and have the following exprssions, Here, M and Q are integration constants related to the ADM mass and the charge of the black brane. q = ( * F ) xyη = 2 √ 3 Q L 3 represents the charge density and r h denotes the position of the event horizon of the black brane.
We deonte by F (r), χ 1 (r) and φ 1 (r), the O(γ) corrections which we get by solving (5) and (7) keeping the terms up to linear order in γ. These are given as where k 1 , k 2 , k 3 and k 4 are dimensionless integration constants. Our next task would be to determine these constants. Following [37], we evaluate them in a standard fashion by imposing several constraints on the above equations. First, we note the asymptotic behaviour of the black brane metric, where, (f e −2χ ) ∞ = lim r→∞ f (r)e −2χ(r) . This metric at r → ∞ represents the background metric where the dual boundary CFT lives. Now, to fix the speed of light in the CFT to be unity, we then demand that (f e −2χ ) ∞ = 1, which in turn gives k 2 = 0. Now our second requirement is that the charge density q remains unchanged. Note that, one can write the Maxwell equation (7) in the form ∇ µ X µλ = 0, where, X µλ is an antisymmetric tensor. Hence, its dual ( * X) xyη is a constant and it is convenient to choose this constant to be the fixed charge density q, i.e., ( * X) xyη = q. Since the quantity ( * X) xyη does not depend on r, we demand lim r→∞ ( * X) xyη = q.
On the other hand, we can always compute this quantity in the asymptotic limit as, Comparing (14) and (15), we obtain k 4 = 0.
The third constraint is that we want to fix the position of the event horizon at r = r h for simplicity. Hence, we need f 0 (r)F (r)| r=r h = 0 which sets, The final requirement is that we need A t has to vanish at the horizon, in order to have a well defined one-form for the gauge field A. This implies φ 1 (r h ) = 0, which in turn fixes the constant k 3 , Since we have determined all the integration constants we can write down the final expressions for F (r), χ 1 (r) and φ 1 (r) which we show here for completeness, The Hawking temperature of this black brane is given by According to the gauge/gravity duality, it represents the temperature of the boundary field theory. Note that, for γ = 0, it reduces to the Hawking tem- between the charge (Q) and the mass parameter (M), Demanding T = 0, one can calculate the extremal charge of the black brane from (18) as, The gauge/gravity duality suggests that the chemical potential µ of the boundary field theory should be identified with the asymptotic value of the time component of the bulk gauge field, φ(r). On the boundary field theory µ has the dimension of energy (i.e., Length −1 ). So, we redefine the gauge field asÃ t = A t /L with some relevant scale 'L' such that the chemical potential has the unit of energy, wherẽ L has the dimension of length. Hence, the chemical potential on the boundary field theory is given by Now, in the boundary field theory we define a uesful dimensionless quantity, We can explore the full range of µ T , i.e., from µ T = 0 to µ T = ∞. The case µ T = 0 corresponds to Q = 0 and the other case µ T = ∞ corresponds to Q = Q ext when the Hawking temperature T vanishes. Figure 1 shows the variation of µ T as a function of the black brane charge parameter Q for different values of the Weyl coupling constant γ. Here we have set the radius of the event horizon r h = 1, the AdS length L = 1 and the relevant scaleL = 1. The green curve corresponds to γ = 0 when the Weyl corrected black brane reduces to a RNAdS black brane. The red and blue curves correspond to γ = 0.02 and −0.02 respectively. For a small value of the charge parameter (say Q < 0.15) the Weyl coupling γ does not affect the µ T ratio. But for a sufficiently large value of the charge parameter the µ T ratio increases as γ increases from a negative to a positive value. Since we are considering only the O(γ) correction to the metric and the gauge field, we will restrict ourselves to small values of γ in all our numerical calculations.

Holographic Thermalization with Weyl Corrections
In this section we briefly recapitulate the thermalization of the strongly coupled quantum field theory dual to the gravitational model introduced in the previous section. In particular, we will discuss how the new Weyl coupling term would affect such a thermalization process. We begin with a zero temperature state of a quantum field theory in four spacetime dimensions. According to the gauge/gravity duality, it is dual to a five dimensional pure AdS space. Now, if we inject some energy to this zero-temperature state, it would evolve, and after a certain time reach the thermal equilibrium at some non-zero finite temperature. This phenomenon is known as the thermalization. After the state thermalizes at some non-zero temperature, using the gauge/gravity duality, it would be identified with a charged black brane in AdS space. The Hawking temperature of this black brane represents the temperature of the final state in the dual field theory. So to describe the thermalization process in the field theory, we need to create an AdS black brane from a pure AdS space in the dual gravity sector. This formation of black branes is well-described in literature [23], [27], [28] by modeling the spacetime with a AdS-Vaidya metric which desribes the collapse of a thin shell of charged matter from the boundary, into the bulk interior. The outer region of this collapsing shell represents a black brane metric while the inner one corresponds to a pure AdS spacetime. For our purpose, to study the effect of the Weyl coupling on thermalization, we need to construct a similar AdS-Vaidya type metric which at early time would correspond to a pure AdS space, and at late times would merge to the Weyl-corrected black brane metric after the shell collapses. Hence, we now focus on constructing the Weyl corrected AdS-Vaidya metric: First, we introduce a new radial coordinate z = L 2 r . Note that the boundary is at z = 0 and the horizon is at z = L 2 r h in terms of this new coordinate. Writing down the metric and the gauge field in the usual Eddington-Finkelstein coordinate yields the following form, Although we have a z-component of the gauge field in the Eddington-Finkelstein coordinate, we can set A z = 0 through a gauge transformation. So, in this new coordinate system the gauge field becomes A = A t dv = φ(z)dv. Now, instead of dealing with a constant mass (M) and charge(Q) parameter, if we assume them to depend on the advanced time coordinate v, i.e., M = M(v) and Q = Q(v), the functions f (z) and χ(z) in the metric would assume the form f (v, z) and χ(v, z) respectively. Also, the gauge field φ(z) should be now written as φ(v, z). After considering M and Q as functions of the advanced time coordinate, the metric would no longer satisfy the Einstein and Maxwell equation. We then need an external matter source to vary M(v) and Q(v), with the advanced time v. Considering this external matter source, the Einstein and Maxwell equation can be written as, where the external matter source must have the following expression for its nonzero components in order to obey the Einstein and Maxwell equation, and Here, the prime denotes a derivative with respect to v. Equation (28) implies that J λ (ext) J λ(ext) = 0. Hence, the matter current, which sources the gauge field, is null which is permissible for a physical matter current (for related discussions, see [46]).
The stress-energy tensor T (ext) µν can be diagonalized to determine the energy density and pressure. Solving the eigenvalue problem, T ν(ext) µ n ν = λ n µ we get the energy density, ρ = 72γ ,z) , 0, 0, 0) as the corresponding eigenvector, which is null. This implies that, the stress energy tensor has support along the lightlike direction. In fact, the stress energy tensor can be written using two linearly independent null vectors n µ Hence, the characteristic surfaces of the Weyl-corrected Einstein and Maxwell equations are lightcones, which supports our analysis considering lightlike collapse.
However, (27) reveals that the infalling matter we are dealing with does not represent a shell of pressure-less null dust as considered in the works of [23] - [31]. Rather, we now have a thin shell of charged null fluid having finite pressure. 3 Therefore, we will consider the collapse of this null shell of fluid into a black brane to study the thermalization. Here, we should mention that the authors of [48] considered the evolution of thin shells made of different kinds of degrees of freedom, to study the dynamics of thermalization. These degrees of freedom are governed by different equations of state. It was shown there that the shells move and collapse with different velocities depending on the equation of state.
In a similar manner one can in principle start with our Weyl-corrected Vaidya geometry, and taking a particular equation of state (e.g., p = aρ with p and ρ being the pressure and energy density, respectively, with in the shell and a is a constant), one can, compute the shell velocity and study the collapse using the Israel junction condition [49], [50]. However, in this paper, we are interested in the effect of the Weyl-coupling constant, rather than different equations of state, or the degrees of freedom that constitute the shell. Now, we must make sure that the energy-momentum tensor supporting the time dependent bulk solution satisfies appropriate energy conditions, namely the null energy condition (NEC): T (ext) µν n µ i n ν i ≥ 0 (i = 1, 2), where n µ 1 and n µ 2 are the null vectors that we have discussed above. 4 It can be checked that the NEC with n µ 1 is satisfied trivially. On the other hand, replacing n µ 2 into the NEC, we have the following constraint condition on the functions f (v, z), χ(v, z) and φ(v, z), where, T are defined in (27). Now, we would choose the mass and the charge function in such a way that it can describe the thermalization process holographically, i.e., at early time (v → −∞), M(v) and Q(v) would represent a pure AdS geometry while at late time (v → ∞) they would imply a black brane geometry. In particular, we would take two smooth functions which are often used in literature [23], [27], [28] for a numerical study of thermalization: where v 0 is the thickness of the null shell. Notice that, as v → −∞, M(v) = Q(v) = 0 and the spacetime represents a pure AdS geometry while as v → ∞, M(v) = M, Q(v) = Q and the spacetime reduces to a finite temperature black brane given by (24). An useful physical situation would be the collapse of a shell having zero thickness (v 0 → 0). Then one can assume that is a step function. v = 0 would represent the position of the shell : for v < 0 the spacetime would be a pure AdS space whereas for v > 0 it would represent a black brane geometry. Substituting these forms of M(v) and Q(v) in (29) one can check that the NEC is valid for any value of v in this case. But one should be very careful on the validation of the NEC while performing numerical computations with the functions given in (30) i.e., with a shell of finite thickness v 0 . We should choose the parameters Q, M and γ in such a way that the NEC is satisfied and in order to satisfy the NEC one can always tune the values of L and r h . However, since the most physical situation would be the case with zero shell thickness, following [27], we scan the full parameter range i.e., from

Non-local Observables
In this section, we will pick out a set of observables to probe the thermalization in the four dimensional boundary field theory. Using the holographic principle, we would relate the particular observables to extended geometric objects in the bulk, representing the gravity dual for those. We would use the gravity duals to compute the particular observables and understand the thermalization. However, we cannot choose any local observable (e.g., expectation value of the energymomentum tensor) to probe the thermalization. We have to compute some nonlocal observables (e.g., correlation function of the local operators) to understand the details of the thermalization process. In particular, we would like to calculate the two-point correlation function of the local gauge invariant operators at a fixed time and the expectation value of a rectangular Wilson loop on the boundary field theory. The gauge/gravity duality provides a nice way to evaluate these non-local observables. The two-point correlation function can be interpreted as a renormalized geodesic length connecting the two points on the boundary. On the other hand, the expectation value of the rectangular Wilson loop corresponds to a minimal area surface extended into the bulk.

Two-point Correlation Function
We want to compute the equal time correlation function of local gauge invariant scalar operators O(t, x) with conformal dimension ∆ and study its evolution with time. As is well known, using the AdS/CFT correspondence and employing the saddle-point approximation for large values of the conformal dimension (i.e., ∆ ≫ 1) one can evaluate this as, where L is the length of the geodesic connecting the two points (t, x) and (t, x ′ ) on the boundary of the AdS space. Hence, for a very large value of the conformal dimension, L is proportional to the logarithm of the two-point correlation function.
We consider two points P (t, − l 2 , y, η) and Q(t, l 2 , y, η) separated by a distance l on the boundary field theory and connect them to construct a spacelike geodesic which would extend into the bulk. From symmetry, it is clear that, the shape of the geodesic would depend only on the boundary spatial coordinate x, but not on y or η. So we would treat x as the geodesic parameter. Then the solution to the geodesic equation is given by two functions v(x) and z(x) with the following boundary conditions : where t is the time at the end point of the interval on the boundary and z 0 is an UV cutoff in the theory. Now, using (24), we can write down the length of the curve connecting the two points on the boundary as, Here, the prime denotes a derivative with respect to x. The two functions v(x) and z(x) would minimize the length L curve and yield the geodesic. Hence, we can think of the integrand of the above equation as the Lagrangian and L curve as the action in the sense of classical mechanics. Note that, the integrand is not an explicit function of x, so there will be a conserved quantity, where, in analogy to classical mechanics, H corresponds to the Hamiltonian. Also note that, x = 0 is the turning point of the geodesics in a sense that the geodesic is symmetric around x = 0. Now, we impose the initial conditions at x = 0 : Then (34) gets simplified, since By employing the Euler-Lagrange equation, we can write down the equations for v(x) and z(x) using (33) and (36), This is a pair of second order non-linear coupled differential equations and difficult to solve analytically. However, we can solve these numerically for different pairs of (v * , z * ), where we are given the initial conditions (35). Finally, we use (36) to write the on-shell geodesic length in a simple form, Note that, L is a function of the separation l between the two points and the time t at those points. The l dependence is clear from the integral, where as the t dependence appears because of the factor z * . For a particular initial value (v * , z * ), we get a particular time from the condition v(± l 2 ) = t. Now, if we vary z * we would get a different value of t.
Notice that, the geodesic length L(l, t) has a divergent piece due to the contribution of the AdS boundary (z = 0). So one needs an UV cut-off (denoted by z 0 before) to get rid of this divergent piece. We define the relevant finite part of the geodesic length by taking out the divergent part in pure AdS geometry as where, 'L ren ' is called the 'renormalized geodesic length' which would represent the renormalized two-point correlation function on the boundary field theory.

Wilson Loop
We will use the Wilson loop as our second tool to probe the thermalization. We will explicitly calculate the expectation value of the Wilson loop operator in the boundary field theory and study its time evolution. The Wilson loop is a non-local gauge invariant observable which is defined as, where, the notation P denotes that we have a path-ordered integral of the gauge field A over a closed loop C and N stands for the rank of the gauge group.
Although it is in general difficult to compute the expectation value of the Wilson loop operator in a quantum field theory, it turns out to be easy using the AdS/CFT duality. It is related to the partition function in a string theory as where, Σ denotes the world sheet extending into the bulk and having the closed loop C as its boundary, i.e., ∂Σ = C. S N G (Σ) represents the string action, known as the Nambu-Gotto action. If the boundary field theory is strongly coupled, we can make use of a saddle-point approximation and write the above partition function as where, α ′ denotes the inverse string tension and A(Σ 0 ) represents the area of the minimal surface world-sheet (Σ 0 ) having the same boundary C. Now we construct a spacelike rectangular Wilson loop C on the x − y plane of the AdS boundary having sides l along the x-axis and R along the y-axis and center at the origin. Further, if we assume translational invariance along the y-axis, the shape of the bulk surface would depend only on x. Then the solution to the minimal area surface can again be given by two functions v(x) and z(x) parameterized only by x. The same boundary conditions as in the earlier case will be applicable here also, v(±l/2) = t , z(±l/2) = z 0 .
Using (24) we can write down the area of the surface Σ as, Like in the previous case of the geodesic, the integrand of the above equation has no explicit dependence on x. Hence, a conserved quantity would be associated with it, namely, Since the turning point of the minimal area surface is again at x = 0 because of the symmetry of the problem, we can impose the initial conditions, Using these initial conditions, (45) simplifies to, We now minimize the area functional given by (44) and get the equations for v(x) and z(x). These are similar to (37) but we show them here for completeness, We will solve these equations numerically the same way we do for the the twopoint correlator. Now we use the conservation equation (47) to express the area of the minimal surface as, Note that the minimal area surface is also divergent because of the z = 0 contribution to the integral. So, we will work with a 'renormalized' notion of the minimal area surface which we get by subtracting the divergent part. We define the 'renormalized minimal area surface' as, where, z 0 is the UV cut-off of the theory we are already familiar with.

Numerical Results
In this section, we provide the numerical details for computing the thermalization time in our Weyl corrected gravity model, and explain the results we get by probing the two-point correlation function and the Wilson loop on the boundary. For this purpose, it is necessary to separately solve the set of coupled differential equations (37) and (48). From now on, we would set the radius of the event horizon r h = 1, the AdS radius L = 1, the UV cut-off z 0 = 0.01 and the shell thickness v 0 = 0.01 in all of our numerical calculations. For a particular value of the charge parameter Q, from (19), we have the mass parameter, M = 1 + Q 2 . Now, for this fixed value of M and Q, we first solve the pair of equations (37) subject to the initial conditions of (35). To extract the boundary time t, we fix the value of v * and tune the value of z * until we get z = z 0 = 0.01 at the end point of the geodesic. For example, setting Q = 0.5 and γ = 0.02, if we take the separation of the two points on the boundary to be l = 3 and fix the initial time v * = −1.42, we get z(± l 2 ) = 0.01 just when the tip of the geodesic reaches the position z * = 1.273155. Thus we get the geodesic profiles z(x) and v(x) for this particular value of (v * , z * ) and the corresponding boundary time would be computed as, t = v(± l 2 ) = 0.0117278. Fixing the boundary separation to be l = 3, we can now take different values of v * and repeat this procedure to visualize the geodesic profiles at different stages of time. Thus we obtain a number of geodesic profiles v(x) and z(x) at different boundary times corresponding to different initial times v * . We can calculate the renormalized geodesic length L ren (l, t) from (39) for all these geodesic profiles, and can study the behaviour of the renormalized geodesic length as a function of the boundary time. At this stage, it is convenient to introduce a dimensionless l-independent quantity Lren We use the same technique to solve the pair of equations (48) and view the profiles of v(x) and z(x) characterizing the minimal area surface for the rectangular Wilson loop. Like the previous case, we calculate the boundary time by using the same boundary conditions (43) and study the time evolution of the renormalized minimal area surface. As before, it is convenient here to talk about a dimensionless renormalized quantity which does not depend on the area of the rectangular boundary Wilson loop. Hence, we will plot ∆A = (A ren − A thermal )/Rl as a function of the boundary time t, where A thermal is the thermal value of the renormalized minimal area.

Two-point Correlation Function and the Renormalized Geodesic Length
Here we study the time evolution of the two-point correlation function by probing the geodesic connecting those two points into the bulk AdS space. Figure 2 shows how the Weyl coupling parameter γ affects the geodesic profiles at different stages of time and thus affects the thermalization. We have fixed the charge Q = 1 and the separation between the boundary points l = 3 and compared the geodesic profiles that appear at different times for different values of γ. The left column represents the time evolution of the geodesics for γ = −0.01, the middle one corresponds to γ = 0 and the right column represents the time evolution for γ = 0.02. Notice that as time elapses, the shell, shown by the dashed green line, approaches z = 1 (shown by the dashed red line) where the horizon of the black brane would be formed at late times. Since we have two different metrics on either side of the shell, whenever the geodesic crosses the shell, it gets refracted by the shell. In the outer region of the shell, the geodesic propagates through an AdS black brane geometry whereas in the inner region, it propagates through a pure AdS geometry. It is the refraction of the geodesic at the shell which makes the dual field theory stray from thermality. As time elapses, the shell comes closer to the position where the event horizon would be formed and the geodesic refraction at the shell becomes more prominent. The time when the geodesic no longer penetrates the shell defines the thermalization time since in this case the geodesic only propagates through the black brane geometry and hence on the dual field theory we have a thermal correlator. Note that the geodesic with γ = −0.01 always penetrates 'more' into the bulk than the other two geodesics. Now compare the three figures of the last row at a fixed boundary time t ≈ 1.35. With γ = 0.02 the geodesic no longer crosses the shell. This suggests that the boundary system has thermalized at this value of γ. But it still needs a small amount of time for γ = 0 and an even larger amount of time for γ = −0.01 for the geodesic not to penetrate the shell. Hence, fixing the charge parameter Q, as we increase the Weyl coupling from a negative value to a positive value, we see that the boundary field theory takes less amount of time to thermalize. This is indicative of how a Weyl correction term can affect the thermalization time in the strongly coupled boundary field theory. Now, we compute the dimensionless l-independent renormalized geodesic length ∆L and plot it as a function of the boundary time t and generate the thermalization curves. We have already discussed the numerical method to calculate ∆L(t) extensively and so here we directly show our results in figure 3. This figure shows how the renormalized geodesic length and hence the two-point correlation function evolves with the boundary time t with the Weyl coupling γ as a parameter.
In figures 3(a), 3(b), and 3(c) we have fixed Q = 0.5, 1 and √ 2 respectively, while the separation between the two boundary points is taken to be l = 3. These figures show that, starting with a negative value, ∆L increases with time and at a certain time it saturates ending up at zero. But, at the beginning of the thermalization process, there is a delay, which was also reported by [23] and [27]. They argued that the delay was because the boundary field theory experiences the sudden injection of energy only at a distance of the order of the thermal wavelength ∼ 1 T . The time, when all the curves reach their corresponding thermal value, i.e., ∆L(t) reaches zero, is referred to the thermalization time. Clearly, it sets a time scale for the Weyl corrected black brane to form. We have zoomed the regions near the thermalization time to analyze the precise effect of γ on the thermalization time.
These are shown in the insets of the corresponding figures. Notice that, when Q is small, γ has little effect on the thermalization time. But for a sufficiently large value of Q, as one tunes γ from a negative value to a positive value, the thermalization time decreases which is also expected from the time evolution of the geodesics as shown in figure 2. It is important to point that when Q is very large, e.g., Q = √ 2, a swallow-tail pattern appears at the end of the thermalization curve with γ = −0.01, whereas, with γ = 0.02, 0.01 and 0 we have no such behaviour in the thermalization curve. One can also check that as γ becomes more negative, the swallow-tail pattern becomes more prominent. This swallow-tail kind of behaviour has been reported in [23], [27], [30] and [34] in the context of holographic thermalization, thermal and electromagnetic quenches. In [23] it was discussed that the emergence of the swallow-tail pattern depends on the dimension of the system while [27] argued that it is rather a universal phenomenon which does not depend on the dimensionality of the AdS space. We should point out here that the swallow-tail behaviour appears for γ = 0 for higher values of the charge Q, for large values of the boundary separation l. This seems to imply that this behaviour may not be related to potential causality violating issues as one switches on the Weyl coupling.
In figure 3(d), we present a zoomed view of the swallow-tail pattern. The pattern emerges because of the presence of three different geodesic profiles at a certain time before the thermalization and these three geodesics simultaneously extremize the bulk action at that particular time. Figure 4(a) shows the time evolution of z * , which is zoomed and shown in figure 4(b). It clearly shows that between t ≈ 1.4315 and t = 1.44, at any particular time, z * has three different values corresponding to the three different geodesics at that time. Because of the multivaluedness of z * (t) one should be careful before applying the saddle-point approximation at the late time. Now consider the dashed magenta line in figure 4(b) connecting the three points at t = 1.4332. We have shown the three geodesics corresponding to these three points in figure 5, where the left figure corresponds to the top point on the magenta line in figure 4(b), the middle one corresponds to the middle point on the magenta line while the right figure corresponds to the bottom point on the magenta line. Note that in figure 5(a) the geodesic crosses the shell and propagates inside the horizon while the shell is just inside the horizon. In figure  5(b) the shell is outside the horizon, the geodesic penetrates the shell but does not cross the horizon. In figure 5(c), the shell lies well outside the horizon and in this case the geodesic does not cross the shell.  In general, the appearance of a swallow-tail is usually an indicator of different scales present in the problem. In this case, if we follow the main curve, ∆L has three branches before and beyond t = 1.4357 (the position of the kink in figure 3(d)). For t = 1.4332, as we have just discussed, two of the geodesics do not penetrate the horizon. This feature in fact continues beyond t = 1.4357, up to t ≈ 1.439, as is evident from figure 4(b). If we assume that the saddle point approximation remains valid in these regions, then this would seem to indicate two physical geodesic solutions corresponding to two different scales in the problem. From a field theory perspective however, the two point correlator should be single valued, and this would suggest the breakdown of the saddle point approximation at late times. This issue is not resolved, and requires further investigation.
For completeness, we also provide the thermalization curves with Q as a parameter at a fixed value of γ. These are shown in figure 6. From figure  6(a), one can see that with γ = −0.01, the thermalization time increases as one increases the charge Q. In other words, if we change the µ T ratio from zero to ∞, the thermalization time enhances. If we zoom the figure to view the behavior of  the curves near the thermalization time, we notice that for the extremal charge Q = √ 2, which corresponds to µ T = ∞, we get a swallow-tail kind of behaviour before the thermalization. For a larger value of γ the change in the thermalization time with Q is negligible. Even in the inset of figure 6(b) we see a negligible change of the thermalization time within a very short region for γ = 0.02. So, we would comment on it on after probing it by another non-local observable, namely, the expectation value of the Wilson loop operator. But it is important to mention that we do not get any swallow-tail appearance in the thermalization curve for γ = 0.02 even with the extremal charge Q = √ 2. At this point, we define a time scale for the thermalization, τ crit , following [24]. It is the critical time when the peak of the geodesic touches the middle of the shell at v = 0. This can be computed as, where, z 0 is the UV cut-off and z * is the value of the z coordinate at the peak of the geodesic. Note that, for a particular boundary separation l, z * would pick a particular value and substituting into the above formula we can determine τ crit for that particular l. In figure 7 the critical thermalization time τ crit is plotted as a function of the boundary separation length l, which reveals that the thermalization is always topdown. A similar result was reported in [23], [27], [28]. This is not surprising and is a natural outcome of the dual geometrical probes we are using. If the boundary separation is small, the geodesic cannot cross the shell and always propagate in the black brane geometry, yielding a thermal correlator on the boundary. But, if the boundary separation is large enough, the geodesic would penetrate the shell and enter into the pure AdS geometry, thus the correlator would not be thermal and it would take a sufficient amount of time to be thermal. Hence τ crit would be larger for a larger l. It is interesting to note that for small values of l, τ crit ∼ l 2 , but as l increases there is a deviation from the linearity. From figure 7(b) it is also clear that, for a sufficiently large value of the charge parameter, the deviation occurs in τ > l 2 for negative values of γ and in τ < l 2 for positive values of γ. The figure also shows that for very small values of l, τ crit has a negligible dependence on γ, whereas, for a larger value of l, τ crit decreases as one increases γ at a fixed value of Q. Now fixing the value of γ, if one increases the charge parameter from Q = 1 to Q = √ 2, τ crit enhances.

Wilson Loop and the Renormalized Minimal Area Surfaces
In this subsection, we consider the time evolution of the minimal area surfaces by probing the Wilson loop into the bulk AdS space. Solving the set of equations (48) with proper initial conditions as prescribed earlier, we generate a sequence of minimal area surface profiles at different times for different values of γ. This has been shown in figure 8 where we have choosen a rectangular Wilson loop of length l = 1.5 and width R = 2 on the boundary of the AdS space and we have fixed the charge Q = 1. The left column corresponds to the time evolution of the minimal area surfaces for γ = −0.01, the middle one represents γ = 0 and the right column corresponds to the time evolution with γ = 0.02. Note that as expected, the time evolution profiles are almost the same as the geodesic profiles in the previous subsection. As time elapses, the shell approaches towards the surface z = 1 where the horizon of the black brane would be formed at late time. Also note that there is a refraction of the minimal area surfaces on the shell when they cross the shell. As γ becomes more negative, the minimal area surfaces penetrate more into the bulk. Now consider the three figures of the last row at a fixed boundary time t ≈ 1.125. With γ = 0.02 the surface propagates only in the black brane and does not cross the shell suggesting that the boundary system has been thermalized. But with γ = 0 and with γ = −0.01 the minimal area surfaces are still crossing the shell and propagating into the pure AdS geometry and hence we do not yet have a thermal Wilson loop on the boundary. Hence, as in the previous subsection, we expect that as we increase the Weyl coupling from a negative value to a positive value, the thermalization in the boundary field theory would also be easier. Now we compute the dimensionless renormalized minimal area δA and generate the thermalization curves to see whether they are consistent with our expectation. This has been shown in figure 9 which describes the behaviour of the renormalized minimal area and hence the boundary Wilson loop with time for different values of the Weyl coupling constant γ at fixed Q. We see a delay at the beginning of the thermalization time as in the case with two-point correlators. Starting with a negative value δA increases with time and reaches zero at the thermalization time.
In each figure, at a fixed value of Q, the thermalization time decreases as we increase the Weyl coupling from a negative to a positive value which agrees with our result for the time evolution of the minimal area surfaces. The insets in each figure contain more details about the swallow-tail appearance. Notice that with Q = 1 we have the swallow-tail appearance only for negative values of γ. Remember that we did not get any swallow-tail emergence in the thermalization curve for the two-point correlators with Q = 1. Further, if we increase the charge to the extremal value Q = √ 2 even γ = 0 and γ = 0.01 shows the swallow-tail behaviour before the thermalization. Hence, we conclude that the swallow-tail behaviour is more prominent with negative values of γ and large values of the charge parameter Q. Also for a given value of γ and Q probing the Wilson loop gives a clearer picture of the swallow-tail emergence than probing the two-point correlator. Figure 10 shows the thermalization curves with Q as a parameter and keeping γ fixed. Notice that for γ = −0.01, as Q increases there is a delay in the thermalization time as in the case with the two-point correlator. We also get a swallow-tail pattern for Q = 1 and Q = √ 2. But as we fix γ to a positive value, γ = 0.02, we see again a negligible change in the thermalization time. But if we zoom the figure in the region just before the thermalization time, we see that Q = √ 2 takes the minimum time to thermalize. Hence with Q = √ 2 even if the initial state of the dual field theory is much away from thermal equilibrium than the states with Q = 1 and 0.5, but a shorter delay time at the beginning of the thermalization and a faster growth of δA than the other two make the thermalization time for Q = √ 2 lesser. Hence for high positive values of γ, it is the large values of the charge parameter, which make the thermalization faster, although the situation is different with negative values of γ. This behavior was also noticed in [28] where for small values of the boundary separation, the boundary field theory thermalizes faster with large values of the charge parameter, while for large separation, the boundary theory thermalizes later with large values of the charge parameter. These authors argued that with fixed value of the charge parameter there may exist two different regimes in the boundary field theory : for small l the boundary theory is in quantum regime and for large l it lies in the classical regime. We have checked that for large value of the boundary separation, it is the smaller values of the charge parameter which makes the thermalization faster. But it is unclear whether there exists a notion of classical and quantum regime, since the system is out of equilibrium. Notice that there is no appearance of the swallow-tail pattern with larger positive values of γ even with the extremal charge (we are assuming here that our perturbative analysis is valid for such values of γ). Now we plot the critical time τ crit as a function of the length of the rectangular Wilson loop l while we fix the Width R = 2. Figure 11 shows that for very small values of l, τ crit is almost independent of γ. The thermalization time no longer behaves as τ crit ∼ l 2 in this case. It is clear that as we increase γ, τ crit decreases and this is more prominent when Q is sufficiently large.

Discussions and Conclusions
In this section we summarize the main results of this paper. We have considered five dimensional AdS gravity coupled to a U(1) gauge field by a combination of two and four derivative interactions, where, the four derivative interaction couples two powers of the Maxwell field to the bulk Weyl tensor. We call the coefficient of this four derivative interaction as the Weyl coupling constant (γ) and have studied how this coupling along with the chemical potential affect thermalization in the dual gauge theory. This gives an extra control parameter that might be important to construct models of thermalization in realistic situations. For this purpose, first we constructed the black brane metric which solves the Einstein and Maxwell equations up to linear order in γ, since an exact analytical solution seems to be intractable. Then we constructed the Vaidya-like dynamical black brane metric which represents a pure AdS space at early times and a Weyl-corrected black brane at late times. We used two non-local observables on the boundary field theory to study the thermalization: equal time two-point correlation functions and the expectation value of the Wilson loop operator. Using the gauge/gravity duality, these two observables were identified to be the dual of geometric quantities in the bulk gravity: the geodesic length and the minimal area surface, respectively, extending into the bulk from the boundary. These were then used to compute several physical quantities associated with the thermalization of the strongly coupled boundary theory.
Broadly, we have presented all necessary physical details of thermalization in Weyl corrected five dimensional AdS gravity duals. Namely, we have analyzed the time evolution of geodesics and the time dependence of the geodesic length for two point correlators. Correspondingly, we have studied the time dependence of minimal area surfaces for rectangular Wilson loops. For both these cases, we analyzed the thermalization time scale. The subtle interplay between the Weyl constant and the chemical potential has been elaborated upon, and we have seen that the thermalization times for strongly coupled field theories may increase or decrease depending on the relative values of the two. An outcome of our analysis was the appearance of a swallow-tail behaviour in the thermalization curve, whose onset is affected by the Weyl coupling, and we have seen evidence that this might indicate distinct physical possibilities relating to different scales in the problem, assuming that the saddle point approximation in computing the thermalization time continues to be valid.
One can also use another non-local observable, the entanglement entropy, to probe the thermalization in the boundary field theory. Using the gauge/gravity duality, it would be dual to the minimal volume extending into the bulk. In [23], it was shown that it is the entanglement entropy which sets the relevant time scale in the problem, since it thermalizes the last. In the present case we found it difficult to generate the thermalization curves for holographic entanglement entropy, since the numerical computations there need a very high working precision. Let us elaborate on this in some details. In Appendix C, we have provided the setup for calculating the Holographic entanglement entropy in the framework of Weyl corrected gravity (up to first order in γ).
In principle, we can compute the holographic entanglement entropy numerically following exactly same procedure as we did for the two-point correlator and the Wilson loop. However, we find that this requires an enhancement of the working precision in MATHEMATICA by a large amount in order for the numerical values to be reliable. For example, as we have already mentioned in section 5, to produce the thermalization curves for the two-point correlator we have to first solve the pair of equations (37) subject to the initial conditions of (35) for a fixed value of γ, Q and M. To extract the boundary time t, we fix the value of v * and tune the value of z * until we get z = z 0 = 0.01 (as chosen in section 5) at the end point of the geodesic. Here, we adjust the value of z * in such a way that we get z = z 0 = 0.01 with l 2 ∈ (1.499999, 1.500001), i.e., we have a tolerance of 10 −6 determining the length of the boundary separation. With this tolerance, we produce sensible results for the time evolution of the two-point correlators.
However, while dealing with the entanglement entropy, we have checked that a tolerance of 10 −6 does not give trustable numerical results (which can be compared for example with the entanglement entropy after setting γ = 0, where a tolerance of 10 −6 produces existing results in the literature). In fact, due to the complexity of the equation, we have to increase the tolerance upto 10 −15 to rely on numerical values produced. Tuning the value of z * with this amount of tolerance is a daunting task. Although we have performed a limited analysis in this regard, and seen indication that the behaviour of the entanglement entropy produces qualitatively similar results as from those of the other probes, a complete analysis seems rather tedious. So, although the procedure to evaluate the entanglement entropy is similar to the computation of the other two non-local probes, we did not establish in details the numerical values of the entanglement entropy here.
As pointed out in section 2 (see discussion following (9)), our numerical analysis of section 5 has a possible caveat. Namely, we have chosen certain small values of the Weyl coupling γ in the absence of a controlled perturbative expansion in powers of the same. Indeed, while a second order γ-corrected metric can be obtained, analysis of thermalization becomes difficult due to the complexity of the equations. However, we have checked that in the present analysis, smaller numerical values of γ (than those chosen here) does not change the qualitative aspects of our analysis, indicating that our numerical results up to O(γ) are trustable for the values of the coupling used in this paper. However we admit that this is a drawback of our numerical analysis.
Arguably, our analysis has been limited to the fact that we have used a bottom-up phenomenological approach. Although understanding thermalization with the most general four derivative action in five dimensional AdS gravity might be difficult, it should be interesting to see if progress can be made on this.
Acknowledgements AD and SM would like to thank Sayantani Bhattacharya for helpful discussions. AD also wishes to thank Damian Galante and Arnab Kundu for helpful email correspondence. SM would like to thank Joydeep Chakrabortty and Indian Institute of Technology, Kanpur for hospitality and financial support. The work of SM is supported by project No. SPO/DST/PHY/20130115.

A Energy-Momentum Tensor
For completeness here we briefly discuss how to derive the bulk energy-momentum tensor making use of the Palatini identities. Using (4), we can write down the action (2) in the following form, Because of the non-zero Weyl coupling constant γ the energy-momentum tensor would get a finite contribution from the variation of the corresponding part in the action, The variation of which gives, Now to evaluate the variations of R µνρλ , R µν and R, we use the Paalatini identities, Also we use the fact that the variation of the Christoffel symbol is a tensor given by Then after few steps one can derive the following expressions, Now having all the expressions for the variations of R µνρλ , R µν and R, substituting them in (A.3) we can derive the expression for the energy-momentum tensor T µν given in (6) and hence write down the Einstein equation.

B Solution with a different metric ansatz
Following [37] one can choose a metric ansatz different from the ansatz chosen in section 2, i.e., one can choose, dr 2 + r 2 L 2 (dx 2 + dy 2 + dη 2 ) , A = (φ(r), 0, 0, 0, 0) . (B.1) Again we have to solve it up to linear order in γ since an exact analytical solution for the metric seems to be imposible. We consider the same form as [37] for f (r), g(r) and φ(r) where f 0 (r) and φ 0 (r) are the zeroth order solutions representing a Reissner Nordström black brane given by where q = ( * F ) xyη = 2 √ 3 Q L 3 represents the charge density and 'r h ' denotes the position of the event horizon.
Again we solve the equations (5) and (7) to linear order in γ as we did in section 2 and get the perturbations F (r), G(r) and ψ(r) representing the O(γ) corrections, where k 1 , k 2 , k 3 and k 4 are dimensionless integration constants. Imposing the same constraints on the above equations as in section 2, one can evaluate those integration constants and write down the final form of the metric perturbations, Now introducing z = L 2 r and writing down this metric and gauge field in the Eddington-Finkelstein coordinate dv = dt − dz f (z)g(z) , (B.6) we have, Again using a proper gauge, we set A z = 0 and the gauge field becomes A = A t dv = φ(z)dv. Now if we consider the mass and charge parameter to depend on the advanced time coordinate v, the function F , G and ψ would also explicitly depend on v. So we need to introduce an external matter source to satisfy the Einstein and Maxwell equation given by (26) where the external matter source would satisfy, Now constructing a null shell of charged fluid with this energy-momentum tensor and checking the null energy conditions with the mass and charge functions given in (30) we can write down a similar set of equations like (37) and (48) for the geodesic and the minimal area surfaces. We noticed that when the boundary time t is small, we can solve the the two coupled equations and calculate δL and δA. But for large value of t the differential equations exhibit some form of stiffness and we could not get a stable solution for z(x) and v(x). So we could not get a complete thermalization curve both for the two-point correlator and the Wilson loop. Hence we switched to a different metric ansatz as explained in section 2 and resolved this issue. With that metric ansatz the differential equations (37) and (48) did not exhibit any kind of stiffness problem. However, we have checked that, both the metric ansatzs give the same results (same up to five decimal places) for small value of t, as expected.

C Holographic Entanglement Entropy
We can use entanglement entropy as another tool for probing the thermalization. For the sake of completeness, we discuss the basic features of the same in this appendix. If our boundary system is divided into two subsystems A and its complement B, the entanglement entropy of the subsystem A is defined as, where ρ A is the reduced density matrix of A, obtained by considering the trace over the degrees of freedom of B, i.e., ρ A = T r B (ρ), where ρ is the density matrix of the full quantum system. It is known that direct calculation of entanglement entropy in a quantum field theory is difficult beyond 1 + 1 dimensions. However, it becomes tractable if one uses the Ryu-Takayanagi formula [53]. Using this formula, the holographic entanglement entropy of the subsystem A that lives on the boundary of the AdS space is Here, G N is the Newton's constant of the bulk theory and Γ A is a codimension-2 minimal-area hypersurface that extends into the AdS bulk, and it shares the same boundary ∂A as that of the subsystem A. However, as is known, this formula is only applicable to static backgrounds in the absence of higher derivative terms in the bulk action. In presence of such higher derivative corrections, the Ryu-Takayanagi conjecture no longer holds. In [54], a formula for holographic entanglement entropy for a general higher derivative gravity theory was derived. This was shown to consist of the Wald entropy as the leading term, with subleading corrections due to the extrinsic curvature. However, one can explicitly check that, the extra four derivative interaction term in our Lagrangian, C µνρλ F µν F ρλ , does not give rise to any subleading term in the covariant expression of the holographic entanglement entropy given in [54]. Hence, for our purposes, we can use the expression of holographic entanglement entropy as where,L is the Lagrangian obtained from (2) and ε µν is the binormal Killing vector, normalised as ε µν ε µν = −2. Thus we have Like the previous cases, we will again have a conserved quantity, given bỹ Since the turning point of the codimension-2 hypersurface is at x = 0 because of the symmetry of the problem, we can impose the initial conditions, v(0) = v * , z(0) = z * , v ′ (0) = 0, z ′ (0) = 0. (C.6) Using these initial conditions, (C.5) simplifies to We now minimize the functional given in (C.4) and get the equations for v(x) and z(x) : We have to solve these equations numerically the same way we did for the the two-point correlator and the Wilson loop. Using the conservation equation (C.7), we express the extremized area as, 6 1 − 12γ z 6 L 10 q(v) 2 1 − 12γ z 6 * L 10 q(v * ) 2 ≡ S(l, t) . where, z 0 is the UV cut-off of the theory. It turns out that numerical computation of the entanglement entropy becomes somewhat difficult in Weyl corrected gravity, even at first order in the Weyl coupling. This is explained in details in the last section of this paper.