Universality in fast quantum quenches

We expand on the investigation of the universal scaling properties in the early time behaviour of fast but smooth quantum quenches in a general $d$-dimensional conformal field theory deformed by a relevant operator of dimension $\Delta$ with a time-dependent coupling. The quench consists of changing the coupling from an initial constant value $\lambda_1$ by an amount of the order of $\delta \lambda$ to some other final value $\lambda_2$, over a time scale $\delta t$. In the fast quench limit where $\delta t$ is smaller than all other length scales in the problem, $ \delta t \ll \lambda_1^{1/(\Delta-d)}, \lambda_2^{1/(\Delta-d)}, \delta \lambda^{1/(\Delta-d)}$, the energy (density) injected into the system scales as $\delta{\cal E} \sim (\delta \lambda)^2 (\delta t)^{d-2\Delta}$. Similarly, the change in the expectation value of the quenched operator at times earlier than the endpoint of the quench scales as $<{\cal O}_\Delta>\sim \delta \lambda (\delta t)^{d-2\Delta}$, with further logarithmic enhancements in certain cases. While these results were first found in holographic studies, we recently demonstrated that precisely the same scaling appears in fast mass quenches of free scalar and free fermionic field theories. As we describe in detail, the universal scaling refers to renormalized quantities, in which the UV divergent pieces are consistently renormalized away by subtracting counterterms derived with an adiabatic expansion. We argue that this scaling law is a property of the conformal field theory at the UV fixed point, valid for arbitrary relevant deformations and insensitive to the details of the quench protocol. Our results highlight the difference between smooth fast quenches and instantaneous quenches where the Hamiltonian abruptly changes at some time.

, δλ 1/(∆−d) , the energy (density) injected into the system scales as δE ∼ (δλ) 2 (δt) d−2∆ . Similarly, the change in the expectation value of the quenched operator at times earlier than the endpoint of the quench scales as O ∆ ∼ δλ (δt) d−2∆ , with further logarithmic enhancements in certain cases. While these results were first found in holographic studies, we recently demonstrated that precisely the same scaling appears in fast mass quenches of free scalar and free fermionic field theories. As we describe in detail, the universal scaling refers to renormalized quantities, in which the UV divergent pieces are consistently renormalized away by subtracting counterterms derived with an adiabatic expansion. We argue that this scaling law is a property of the conformal field theory at the UV fixed point, valid for arbitrary relevant deformations and insensitive to the details of the quench protocol. Our results highlight the difference between smooth fast quenches and instantaneous quenches where the Hamiltonian abruptly changes at some time.

Introduction
In recent years, there has been a great deal of interest in studying quantum quenches [1], i.e., studying of the quantum evolution of an isolated system in the presence of a time dependent parameter in the Hamiltonian. Amongst other things, these processes are theoretically interesting as probes of two related issues: thermalization and critical points. Considering the first of these, suppose we start with a system in its ground state. If a parameter in the Hamiltonian, e.g., an external field, undergoes a rapid change, the system driven to some highly excited state but one would expect that after sufficient time the system will approach a steady state which resembles a thermal state. The question then is to understand the sense in which the final pure state is close to a thermal state, and to understand the approach to such a state. Similar questions can be studied in a thermal quench, where the initial state is a thermal state. Of course, these questions lie at the heart of the foundations of statistical mechanics and they are typically difficult to investigate, especially when the system is strongly coupled. Recent experiments with cold atom systems and heavy ion collisions are beginning to yield valuable experimental insights into such processes, which pose both greater motivation and interesting challenges for the theoretical community. A second class of interesting quenches are those which cross a critical point. That is, suppose the time dependent parameter passes through a value which would correspond to a critical point in equilibrium. One would then expect that the subsequent evolution of the system will carry universal signatures of the critical point. An early example of such a signal is Kibble-Zurek scaling [2,3]. Suppose one starts in a gapped phase of the system, with the quench rate slow compared to the scale set by the gap. Initially the evolution of the system would be adiabatic. However, as the parameter approaches the critical point, the instantaneous gap vanishes and adiabaticity is lost, producing an excited state. Kibble [2], and subsequently Zurek [3], argued that the density of defects at late times scales as a universal power of the quench rate with the exponent determined by the equilibrium and near-equilibrium critical exponents. In recent years, this argument has been extended to quantum phase transitions and the same arguments have been shown to lead to scaling of other one point functions and correlation functions [4]. The arguments which lead to Kibble-Zurek scaling are based on rather drastic approximations; nevertheless, there are several model systems where such scaling appears to hold. There is no theoretical framework analogous to the renormalization group which justifies such scaling, and strongly coupled systems remain beyond the reach of current theoretical tools. At the other extreme, Cardy, Calabrese and Sotiriadis [5,6] derived a set of exact universal results for instantaneous quenches in two-dimensional field theories from a gapped phase to a critical point, us-ing powerful methods of boundary conformal field theory. Yet another set of scaling relations can be derived from time dependent perturbation theory when the amplitude of an instantaneous quench to a critical point is small [7].
In the past few years, the AdS/CFT correspondence has been used to study both quantum and thermal quenches in strongly coupled quantum field theories which possess a gravity dual. In this approach, the couplings in the field theory are related to boundary conditions for the metric and other fields in the dual gravity theory. Therefore studying a quench process reduces to solving of a set of partial differential equations with specified initial conditions and time dependent boundary conditions -a problem which is much easier to tackle than the original quantum problem in a strongly coupled field theory. The dual description of thermalization becomes the collapse of an incoming shell leading to the formation of a black hole horizon [8][9][10]. One of the interesting results which emerged from these studies is that few body correlation functions thermalize rapidly -a phenomenon which is indeed observed in heavy ion collisions. For quenches across critical points, holographic studies point towards a mechanism for emergence of scaling solutions in the critical region [11] and has led to novel dynamical phases [12]. Further, progress has been made towards observing Kibble-Zurek scaling of defect densities in symmetry breaking phase transitions [13].
Recently, holographic studies also revealed a new set of scaling relations in the early time behaviour of fast but smooth quenches in a critical theory deformed by a relevant operator O ∆ with conformal dimension ∆ [14,15]. The quenches in question involve introducing a time dependent coupling λ(t) for the latter operator. If the coupling varies by an amount δλ in a time δt, a fast quench means δλ (δt) d−∆ 1 . (1.1) In this fast regime, studying quenches where the relevant coupling goes from being zero initially to δλ at late times, it was found that the change of the holographically renormalized energy density δE scales as Similarly, the peak of the renormalized expectation value of the quenched operator was found to scale as O ∆ ren ∼ δλ (δt) d−2∆ , (1.3) consistent with certain Ward identities. These same results also hold for reverse quenches where the relevant coupling goes from δλ at early times to zero at late times. For ∆ > d/2, this implies that δE and O ∆ grow with the quench rate, i.e., as δt shrinks. In fact, the growth in O ∆ is enhanced by a logarithmic factor for even d and integer ∆ and for odd d and half-integer ∆. Implicitly, eqs. (1.2) and (1.3) indicate that for ∆ > d/2, these quantities diverge in the limit of an infinitely fast quench, i.e., with δt → 0. Hence these results seem to be at odds with known results for instantaneous quenches, e.g., [5] - [7]. In these works, a parameter in the Hamiltonian is taken to change instantaneously from one constant value to another value at some time t 0 , and the dynamics is treated in the sudden approximation. This means that in the Schroedinger picture, the state at t = t 0 is treated as an initial condition for standard evolution by the new time independent Hamiltonian. Naïvely, one may think that such an instantaneous quench should correspond to the δt → 0 limit of a smooth quench but this is clearly not the case since in the setup just described, the renormalized expectation values are certainly not divergent.
Of course, the holographic studies [14,15] were implicitly considering strongly coupled quantum field theories whereas the work on instantaneous quenches typically considered free (or weakly coupled) field theories, e.g., [5,6], except in two space-time dimensions. Hence, one possibility is that the new divergences appearing as δt → 0 are only a feature of the special class of strongly coupled theories which have gravity duals. However, we recently showed that this is not the case [16]. In fact, precisely the same scaling as in eqs. (1.2) and (1.3) was found to be exhibited in mass quenches of free field quantum theories. Further, we argued that this behaviour is rather generic. In the present paper, we provide more details of the calculations presented in [16] and report several new results. We also provide a new argument that the universal scaling in the early time response shown in eqs. (1.2) and (1.3) holds for a quench from any constant value of the relevant coupling λ 1 to any other value λ 2 as long as the time scale δt is small compared to all other physical length scales in the problem, , δλ 1/(∆−d) (1.4) In the following, we first consider free bosonic and fermionic field theories in arbitrary dimensions, with a time-dependent mass which evolves smoothly in some time interval δt. We consider a variety of different protocols, i.e., profiles for m(t), which allow us to solve the problem exactly for arbitrary δt. Hence, we are able to calculate O ∆ for finite δt and then examine the result in the fast regime where δλ (δt) d−∆ 1. We find that the (renormalized) expectation value indeed obeys the universal scaling law (1.3), originally found in the holographic models studied in [14,15]. Our analysis clearly exposes the difference between fast but smooth quenches arising in the limit δt → 0, and instantaneous quenches where one works with the sudden approximation. Since we are considering a quantum field theory, the quench rate 1/δt and the quench amplitude m (the initial mass) are not the only scales in the problem. There is, in addition, the UV (momentum space) cutoff Λ. Implicitly, our fast quench limit involves a quench rate which is large compared to the initial mass but still small compared to cutoff, i.e., Λ 1/δt m . (1.5) Although the quench rate never appears in the discussion of instantaneous quenches, they can be considered as having 1/δt ∼ Λ. However, local quantities like O ∆ receive contributions from all scales, and are therefore sensitive to whether or not the quench rate is comparable to the cutoff scale. Indeed we show explicitly that the correlation functions of individual momentum modes for fast and smooth quenches reduce to those in the instantaneous quenches only when the quench rate is large compared to the momenta -hence matching local quantities would require rates comparable to the cutoff scale. The regime of our interest is quite distinct from the latter and arguably more physical. Nevertheless, we expect that for certain quantities, e.g., correlation functions at finite distances, the results for both types of quenches will agree when the distance is large compared to δt since one expects that only small momenta contribute to the result. Similarly, one might expect that the results of smooth fast quench should agree with those of instantaneous quench at late times, t δt. For free fields we will find that the late time behavior indeed agrees for d = 3. However, in higher dimensions the late time results for smooth and instantaneous quenches differ [17]. While we trace the technical origin of the difference, we do not have a good physical reason why this should be the case.
A key ingredient in our work 1 is the renormalization of the underlying quantum field theory. The bare quantity O ∆ is, of course, UV divergent and we need to add suitable counterterms to extract physical quantities at resolutions much coarser than the cutoff scale. Our problem is quite similar in spirit to quantum field theories in curved space-times, e.g., see [18][19][20]. In that case, the required counterterms involve operators made out of quantum fields, as well as curvature tensors of the background space-time. Further, in this context, diffeomorphism invariance provides an important guide restricting the form of the counterterms, which may appear. In the present case with a time dependent mass, we find that we need to add counterterms which involve time derivatives of the mass function, in higher dimensions (d ≥ 6) where stronger divergences appear. Further, the underlying theory is invariant under coordinate transformations if we treat the mass as a background scalar field. Hence diffeomorphism invariance is again a useful guide in restricting the form of the required counterterms.
However, we are still left with the problem of determining the precise coefficients of the counterterms which render the renormalized observables finite. We find that these coefficients can be determined by examining the quenches in an adiabatic limit. 2 That is, the counterterm coefficients determined for an adiabatic quench still remove all of the UV divergences in O ∆ for fast (but smooth) quenches. We argue that this result can be anticipated as follows: in renormalizing the theory, we are always considering quench rates 1/δt which are much smaller than the UV cutoff Λ. In this situation, we expect the high momentum modes, near the cutoff scale, do not care if the quench rate is large or small compared to the mass. Hence any UV divergences should be the same in fast quenches with 1/(δt) m and in adiabatic quenches where 1/δt m. Of course, the latter adiabatic limit is relatively straightforward to analyze since one is performing an expansion in derivatives with respect to time.
To conclude the introduction, we outline our key results and provide their locations throughout the paper.
1. We show in detail that the adiabatic expansion provides the correct counterterms which renormalize one point functions for free bosonic and fermionic theories for time-dependent masses. The bosonic case is discussed in section 2.1 while the fermionic case is contained in section 3.
the UV cutoff scale. The comparison between the fast smooth quenches and the instantaneous quenches will be discussed in much greater detail in [17]. 5. In section 2.7, we compare the late-time response (i.e., t δt) of a smooth fast quench with that of an abrupt quench for free bosonic field theory. In particular, we explicitly show that for d = 3, that the response is independent of δt at late times and leads to a logarithmic growth of φ 2 with time, in exact agreement with the abrupt quench result. For d = 5, we show once again that at late times, the δt → 0 limit is smooth. 6. In section 4, we argue that the universal scaling discussed in this paper is a property of any quantum field theory whose UV limit is a conformal field theory, e.g., a conformal field theory in any number of dimensions deformed by a relevant operator. For quenches which take the system from any nonzero value of the corresponding coupling to some other value of the coupling this universal scaling holds for early time response -so long as the time scale of quench is the smallest physical scale in the problem, as in eq. (1.4). The scaling is purely a property of the UV conformal field theory.

Quenching a free scalar field
We start by analyzing mass quenches for the simple case of a free scalar field φ in d spacetime dimensions, i.e., d − 1 spatial dimensions. In particular, we focus on varying the mass with the following profile: Hence the mass goes smoothly from the value m 2 (A−B) in the infinite past to m 2 (A+B) in the infinite future but the transition occurs essentially in a time period of duration δt centered around t = 0. 3 While much of our discussion does not depend on specific values of A and B, we will begin with a discussion of the case A = −B = 1/2, with which the theory is massive with mass m 2 in the past and becomes massless in the future. 4 As we will show that the scaling behaviour in eqs. (1.2) and (1.3) is recovered with this particular choice.
In section 2.3, we also examine quenches with the mass profile m 2 (t) = m 2 sech 2 (t/δt) (2.2) where the mass vanishes both the infinite past and the infinite future. We again find that the renormalized expectation values show the same scaling as in eqs. (1.2) and (1.3). Finally in section 2.4 we show that the analysis easily extends to general A and B and we again find the same scaling as long as the coefficients satisfy in accord with eq. (1.1), and also The particular protocols or mass profiles in eqs. (2.1) and (2.2) were chosen because they allow us to completely solve the corresponding quantum field theory. That is, the mode functions for the scalar field can be written in closed form, as we will show below. In fact, for the profile (2.1), we can use results first derived in studying quantum fields propagating in curved spacetimes [18,19]. Specifically, in that case, scalar field was examined in an expanding flat Freedman-Robertson-Walker cosmology, which corresponds to a conformally flat geometry described by metric (2.5) A (minimally coupled) free massive scalar field φ, with a constant mass m, propagating in this cosmological background obeys the equation of motion where denotes the ordinary flat space d'Alembertian. That is, the scalar field equation in this curved geometry is identical to that of a scalar field in flat space but with a time-varying mass m 2 (t) = m 2 a 2 (t). Further, it was noted in [18,19], that with a 2 (t) = (A + B tanh t/δt), i.e., with the mass profile (2.1), the corresponding mode functions are given in terms of hypergeometric functions. Hence we may use these results but now interpret the theory as a scalar field undergoing a mass quench. It is also important to mention that with these closed form solutions, we are able to study the behaviour of the theory for arbitrary quenches rates 1/δt and hence we can take the limit δt → 0 to approach an instantaneous quench.
Let us begin with analyzing the theory with the mass profile (2.1). We start by decomposing our field in mode functions As a boundary condition, we will choose the u k to be the in-modes which behave as plane waves in the infinite past. Similarly there will be a corresponding set of outmodes which become plane waves in the infinite future. The operators a k above are then defined to annihilate the in-vacuum, i.e., a k |in, 0 = 0. Exact solutions for these in-modes are [18,19] where 2 F 1 is the usual hypergeometric function and 9) ω ± = (ω out ± ω in )/2 .

Regularization and Renormalization
The quantities we are interested in involve a sum over all modes and are typically UV divergent and need to be renormalized by adding suitable counterterms. In this subsection we show how this can be done. The discussion is valid for generic m(t)in fact we will find the counterterm in terms of the function m(t) and its derivatives. However it is useful to begin the discussion with the mass profile (2.1).
First focus on the case where A = −B = 1/2, in which case we have ω in = k 2 + m 2 and ω out = | k|. Now we adopt the perspective presented in the holographic analysis of [14,15] in the following. In particular, we think of the scalar field theory as a CFT deformed by the operator O ∆ = φ 2 , with conformal dimension ∆ = d − 2. Further, the quenches are made by varying the corresponding coupling in time, i.e., λ(t) = m 2 (t). Our first calculation will be to determine the expectation value of O ∆ = φ 2 , which is straightforward given the mode decomposition above Of course, this expectation value (2.10) contains UV divergences associated with the integration of k = | k| → ∞. The standard approach to deal with these UV divergences is to add suitable counterterms involving the time-dependent mass to the effective action, as in the holographic renormalization of [14]. We turn to the determination of the counterterms in section 2.1.3. However, as described in [16], it is straightforward to find the counterterms which render the expectation value (2.10) finite. Hence let us write the renormalized expectation value as . (2.11) where f ct (k, m(t)) designates the counterterm contribution and Ω d−2 denotes the angular volume of a unit (d-2)-dimensional sphere, i.e., As a first attempt to evaluate f ct (k, m(t)), we might naïvely think that the counterterm contributions needed to regulate φ 2 are those related to the divergences in the constant mass case. That is, with a constant mass, we can identify the UV divergences by expanding With the simple substitution m 2 → m 2 (t), we might then conjecture that eq. (2.11) becomes finite with where we would only include the terms proportional k n with n ≥ −1. As we will see below, this conjecture is only correct for d ≤ 5. For higher spacetime dimensions (i.e., d ≥ 6 in the scalar case), new counterterms involving time derivatives of the mass are allowed by dimensional counting. For example, in d = 6, the term proportional to m(t) 4 /k is associated with a logarithmic divergence in φ 2 . However, by dimensional analysis, f ct (k, m(t)) could also contain a term of the form ∂ 2 t m(t) 2 /k, which might cancel a new logarithmic divergence proportional to ∂ 2 t m(t) 2 in d = 6. Of course, in the case of a constant mass (2.13), no such divergence appears but in the present case of a mass quench, a new UV divergence of this form will be found. As we go to higher and higher dimensions, the set of dimensionally allowed terms involving time derivatives of the mass quickly grows and in fact, the corresponding divergences (generically) do appear, as we will see below. However, let us note that the same dimensional arguments would have identified a potential contribution of the form ∂ t m(t) 2 /k in d = 5 but no corresponding divergence is found. Hence this makes evident that these terms are subject to constraints beyond simple dimensional analysis. In particular, we will show that this single-derivative contribution can be ruled out by diffeomorphism invariance.
Finally, let us comment that in holographic calculations [14,15], these kind of terms naturally appear since couplings are not just constants but boundary values of spacetime-dependent bulk fields. Holographic renormalization then requires introducing counterterms in the gravitational action constructed out of derivatives of the boundary values.

Regulating the theory using an adiabatic expansion
An elegant way to find the necessary counterterm contributions is to look at the divergences appearing in eq. (2.10) for an adiabatic quench, i.e., an infinitely slow quench. In that way, one can organize all contributions with an adiabatic expansion and exactly find the divergent pieces. The discussion below is for a general function m(t).
The adiabatic expansion is an expansion in time derivatives, more precisely in powers of ∂ n t m/m n+1 1. These ratios are, of course, small if the time variation of the mass is infinitely slow. In a generic quantum mechanical system, this expansion is achieved by expanding the state as a linear superposition of instantaneous eigenstates and solving the resulting differential equations for the coefficients in a derivative expansion. For a free field theory, the procedure is easier -one can obtain mode solutions of the equations of motion for each momentum mode, 15) in a WKB type approximation. That is, we wish to find solutions of this equation which are of the form (2.16) Demanding that this ansatz solves eq. (2.15) requires that Ω k satisfies The adiabatic expansion is then obtained in eq. (2.17) by expanding the solution as where Ω (n) k is n th order in time derivatives. We can now substitute this expansion into eq. (2.17) and solve it order by order. The first two orders are trivial, yielding where the latter yields Ω (1) k = 0. The next two orders produce k Ω (1) which are solved by and Ω Again substituting these results into eq. (2.17), we find the next order equation (2.23) As we will see, it is enough to expand up to this order to get all the necessary counterterm contributions to regulate present theories up to d = 9. Now, we want to extract the large-k behaviour of and so we will need to expand 1/Ω k for large k, as well as in time-derivatives. Using 25) where each line in the last expression corresponds to a particular order in time derivatives, e.g., the first line is zeroth order; the second line, second order; etcetera. The ellipsis at the end of each line indicates terms that are higher order in 1/k, i.e., 1/k 8 and higher. Multiplying by k d−2 , those are all the divergent terms in spacetime dimensions less or equal to d = 9. We can see that the first line corresponds to the terms discussed in eq. (2.13). But this is only the zeroth-order adiabatic approximation and there are additional divergent terms at higher orders in the expansion in time derivatives. 5 Of course, the results match those reported in [16], where we found For a fixed spacetime dimension d, we would only keep the terms up to the power k −1 and drop any terms with more negative powers of k. Again, the contributions explicitly written above are sufficient to regulate theories up to and including d = 9. The above discussion applies for a general space-time dimensions, however, we should distinguish between odd and even dimensions. For odd d, all of the powers of k appearing in eq. (2.26) are even (or zero) and f ct essentially subtracts a series of power-law divergences Λ n , where Λ is the UV cutoff scale. When d is even, the powers of k are now odd, and similar power-law divergences are appearing for the positive powers of k. However, apart from these divergences, we may also find a logarithmic divergence when eq. (2.26) contains a 1/k term. If we considered this term alone, the k integral in eq. (2.11) is divergent both in the UV and in the IR. Hence, we also need to introduce a lower bound µ i for each such integral, which then yields log(Λ/µ i ). Hence we see this amounts to introducing an extra renormalization scale in defining the renormalized expectation value (2.11) for even d. The appearance of these new scales reflects certain scheme-dependent ambiguities in defining the renormalized theory, and in particular, as observed in previous holographic studies [14], new ambiguities can arise with time-dependent couplings. Of course, the potential Λ divergences are all eliminated in eq. (2.11) and we take the limit Λ → ∞ in evaluating the renormalized expectation value. Hence in the final result, an infrared scale must replace the UV cutoff in the logarithmic dependence on the renormalization scale, e.g., log(δt µ i )see sections 2.2 and 2.2.2 for further discussion.
With the subscript on µ i , we are emphasizing that in principle one can introduce a separate renormalization scale for each such integral corresponding to a separate coun-terterm. For example, with d = 6 in eq. (2.26), there can be a separate renormalization scale associated with the integrals proportional to m 4 (t) and ∂ 2 t m 2 (t), since they correspond to contributions coming from distinct counterterms -see section 2.1.3 for further discussion. However, in our explicit calculations in the following, we will set all of these scales to be equal, i.e., µ i = µ. The effect in the computation is to divide the integral in the expectation value (2.11) into two parts. The first, from k = 0 to k = µ does not include the 1/k contribution in f ct while in the second, from k = µ to k = ∞, we use the full expression for f ct including the 1/k term. Now we claim that the large-k terms appearing in the adiabatic expansion provide the correct counterterm contributions to regulate φ 2 for general quenches. This claim may seem surprising since the adiabatic expansion should be only valid for slow quenches. However, one can easily verify numerically that with eq. (2.26), the renormalized expectation value (2.11) is finite, e.g., with the mass profile (2.1) even outside of the adiabatic regime. The point is that we are considering a quench rate 1/δt which is always slow compared to the UV cutoff scale, though it may be fast compared to m, e.g., as described by eq. (1.5). The condition for validity of the adiabatic expansion isω k ω 2 k . For this condition to hold for all k, we must have m δt 1. However, for high momenta k m, this condition still holds as long as k δt (m/k) 2 , which is always satisfied for sufficiently large k. Hence the fact that we are interested in studying fast quenches where m δt 1 does not matter for the very high momentum modes, whose contributions are producing the UV divergences. This explains why the adiabatic expansion provides a consistent and convenient framework to find the divergent pieces of the expectation value. In fact, the counterterms are universal and, in particular, independent of the rate at which the mass varies.

Explicit verification for tanh profile
We now show explicitly that eq. (2.26) provides the correct counterterm contributions for general quenches, we return to the tanh profile in eq. (2.1) with A = −B = 1/2. Recall that in this case, the bare expectation value is given in eq. (2.10) where the details of hypergeometric functions appear in eqs. (2.8) and (2.9). Now we proceed to expand these hypergeometric functions for large momentum. In the series representation, the hypergeometric function is defined as 27) where (x) n ≡ x(x + 1) · · · (x + n − 1) and (x) 0 = 1. Further a, b, c, z are given in eq. (2.8). In particular, we recall that the argument z is given by z = 1 2 (1 + tanh(t/δt) ) , (2.28) This variable is the same regardless the choice of A and B. Now we can expand ω in and ω − for large k and see how the first few terms of this series behave. Then we have to take the absolute value squared to get the counterterms for the expectation value of eq (2.10). By checking the behaviour of the series, it can be verified that each successive term begins with a lower power of k. Hence in order to get all the divergent terms up to d = 9, it is sufficient to work with only the first five terms in eq. (2.27).
Focusing again on the mass profile (2.1) with A = −B = 1/2 and expanding these terms for large k, we find At first sight this expression does not look similar to eq. (2.26), but we will now show that they are both actually the same. To start with, we should notice that we can write m 2 (t) as a function of z as m 2 (t) = m 2 (1 − z). Then, for instance the k d−5 term in eq. (2.29) is just −m 2 (t)/2, matching the corresponding term in eq. (2.26). The same happens with all the terms that are independent of the value of δt; i.e., they give m 4 (t) and m 6 (t), as they should. The appearance of terms which are inversely proportional to δt reflects the appearance of time-derivatives in those terms. In fact, we can use trigonometric identities to express derivatives of the mass in terms of powers of the same mass function. This is because derivatives of the hyperbolic tangent are formed by terms proportional to tanh and sech. For instance, the first derivative of m 2 (t) gives ∂ t m 2 (t) = − 1 example, ∂ 2 t m 2 (t) = 4m 2 δt 2 z(1 − 3z + 2z 2 ) and up to an extra minus sign, this expression matches exactly the last three terms appearing in the k d−7 term of eq. (2.29). In the same way, we can translate all the terms of eq. (2.29) to match the universal form that we found in eq. (2.26).
We emphasize that the above calculations are valid for any value of δt and hence this verifies that a single set of counterterms can be chosen to regulate φ 2 independent of the quench rate. In particular, the same counterterms should be valid in the limit δt → 0. We have also performed the same calculations with expanding the hypergeometric function for the reverse quench and found the same counterterms, now as functions of the new m 2 (t) = m 2 2 (1 + tanh(t/δt)). We will also see that for a pulsed quench (2.2), as studied in section 2.3, the same counterterm contributions (2.26) again regulate the expectation value for any value of δt. Hence all of these examples provide a verification of our claim that studying an adiabatic quench is sufficient to determine the correct counterterm contributions to regulate φ 2 for general quenches.

Counterterms in the path integral
Up to this point, we have been interested in finding the necessary contributions which render eq. (2.10) finite and allow us to calculate the renormalized expectation value in eq. (2.11). However, we may also be interested in computing other observables, e.g., the expectation value of the energy-stress tensor -see section 2.2.4. Of course, expectation values of other operators will again generally be UV-divergent and also need regularization. The point we would like to emphasize in this section is that all such divergences should be eliminated by a common set of counterterms regulating the effective action or partition function. Once we have the regulated partition function, we can find the renormalized expectation value of the operators of interest by taking functional derivatives with respect to the appropriate sources. Suppose the path integral is regulated by a UV cutoff Λ, then we have which includes the free field action and the counterterm action 6 where R is the Ricci scalar of the metric g µν and s ij are finite numbers. Of course, for a fixed dimension d, we only retain the terms above with positive powers of Λ and in cases, where the naïve power is zero, it should be replaced by a logarithmic divergence log(Λ/µ) -as discussed in the previous section. Now the expectation values of the 'mass operator' and the stress tensor are given by (2.34) Again we have explicitly shown all of the possible counterterms in eq. (2.32) which would be needed to regulate these two expectation values up to d = 9. Now several comments are in order: First we have introduced a background curved space metric in the partition function (2.30), even though we are evaluating the final expectation values in flat space. This is, of course, because the metric serves as the source of the stress tensor as in eq. (2.34). Further, in this vein, we have included counterterms linear in Ricci scalar in eq. (2.32) since even though these terms vanish in flat space, their variation still contributes to regulating the expectation value of the stress tensor in eq. (2.34). Of course, these terms are not needed to evaluate φ 2 ren in eq. (2.33). We have ignored terms involving higher powers of the Ricci Scalar since they do not contribute to the two one-point functions in eqs. (2.33) and (2.34). Further we have dropped any total derivative terms in the counterterm action, as well as terms that can be related to those appearing in eq. (2.32) by using integration by parts and the identity ∇ µ R µν = 1 2 ∇ ν R. As a result, we were able to eliminate any counterterms linear in the Ricci tensor.
Implicitly above, we are treating the mass-squared as a background scalar field which is a function of all of the spacetime coordinates, i.e., m 2 = m 2 (x µ ). For example, this assumption is evident in eq. (2.34) where the variation yields the expectation value of the local operator operator φ 2 (x µ ). Now if the path integral (2.30) is performed with a covariant action, the counterterm action, as well as the entire partition function, will be diffeomorphism invariant, as assumed with the presentation in eq. (2.32). In particular, the derivatives of the mass only appear there as powers of the covariant d'Alembertian operator. 7 Of course, in applying this counterterm action to study (global) mass quenches, we only consider the mass to be a function of time but the structure revealed here readily explains why all of the counterterm contributions in eq. (2.26) have an even number of time derivatives. We might also comment that in the curved background geometry we have and hence these derivative terms also contribute nontrivially to regulating the stress tensor in eq. (2.34). Let us observe that there are four terms at order k d−9 in eq. (2.26) but only three corresponding counterterms at order Λ d−8 in eq. (2.32). Hence the four counterterm contributions are not all independent. In fact, it is straightforward to show that for a time-dependent mass, the variation of the counterterm with s 41 m 4 m 2 is proportional to , which has precisely the ratio of coefficients with which these two terms appear in eq. (2.26). In fact, by carefully comparing eqs. (2.26) and (2.32), we can identify the coefficients: In principle, the adiabatic expansion in the last subsection could also be used to find the remaining coefficients in eq. (2.32), which would be needed to regulate the expectation value of the stress tensor (2.34) -in dimensions up to d = 9. However, as we will explain in section 2.2.4, we can avoid this calculation, at least in evaluating the expectation value of the energy density (the tt component of the stress energy tensor). The latter can be related to φ 2 ren using a diffeomorphism Ward identity. We will explicitly apply this approach for d = 5 in section 2.2.4 and for d = 3 in section 2.2.5.

Response to the mass quench
In this subsection we calculate renormalized quantities which measure the response to a mass quench of the form (2.1) with A = −B = 1/2.

Numerical results
Given eq. (2.11) for the renormalized expectation value and eq. (2.26) for the necessary counterterm contributions, we are in position to compute φ 2 ren for spacetime dimensions from d = 3 to 9. We first perform this computation numerically. The evolution of the resulting expectation value is shown in figs. 1, 2 and 3 for different values of the quench rate δt. In these plots, the expectation value for an 'adiabatic' quench is subtracted, where the latter actually corresponds to δt = 10. We have verified that the expectation value is essentially independent of δt for larger values. Further, as discussed in section 2.1.1, regulating the expectation value in even dimensions requires the introduction of additional renormalization scales. In the plots presented here, we have set all of these to one, i.e., µ i = 1. Further we have also set m = 1 in the mass profile (2.1). Renormalized expectation values φ 2 ren as a function of time t/δt, for d = 3 and 4. In each plot, the different curves correspond to different quench rates: δt = 1/1, 1/2, · · · , 1/10 where the curves exhibiting higher peaks (in absolute value) correspond to smaller values of δt. Note that the expectation value is multiplied by the numerical constant Further, at each time, the expectation value for an 'adiabatic' quench is subtracted. ren as a function of time t/δt, for d = 5, 6 and 7. In each plot, the different curves correspond to different quench rates: δt = 1/1, 1/2, · · · , 1/10 where the curves exhibiting higher peaks (in absolute value) correspond to smaller values of δt. Note that the expectation value is multiplied by the numerical constant Further, at each time, the expectation value for an 'adiabatic' quench is subtracted.
We can see in figs. 1, 2 and 3 that the peaks in the expectation value grow (in absolute value) as δt becomes smaller, and that this growth becomes even faster when the spacetime dimension is increased. To quantify the growth more precisely, fig. 4 shows φ 2 ren (t = 0) over a broad range of δt, going from δt −1 = 1 to δt −1 = 200, in a log-log plot for d = 3 to 9. Furthermore, for each value of d, the linear fits were made to the curve and the results indicate that the expectation value scales as φ 2 ren ∼ δt 4−d for small δt. 8 For the special case of d = 4, where this formula seems to indicate no scaling, we found that there is actually a logarithmic scaling. Both of these facts match 8 Recall that we have set m = 1 and hence δt 1 should be interpreted as m δt 1, in agreement with the fast quench condition (2.3). the scaling found in holographic analysis [14,15]. In particular, the quenched operator is φ 2 with ∆ = d − 2 and hence the exponent in eq. (1.3) becomes d − 2∆ = 4 − d, precisely the scaling found with the linear fits. Further given that ∆ is an integer, the holographic results also suggest that there should be an extra logarithmic enhancement for even dimensions [15], i.e., φ 2 ren ∝ δt 4−d log(δt). The logarithmic scaling found for d = 4 certainly agrees with this expected enhancement, although there was no evidence of such an enhancement in d = 6 or 8. In section 2.2.2, we will see this occurs simply because for the particular tanh profile, the logarithmic contribution simply vanishes at t = 0. Fig. 5 shows similar plots of φ 2 ren (t = δt/2) over a broad range of δt for d = 6 and 8. There the fit with the extra logarithmic enhancement is clearly preferred over the linear fit. 9 Hence we have found that effectively the mass quenches of a free scalar theory quench reproduces precisely the same early time scaling that was discovered with a holographic analysis [14,15]. Renormalized expectation values φ 2 ren as a function of time t/δt, for d = 8 and 9. In each plot, the different curves correspond to different quench rates: δt = 1/20, 1/21, · · · , 1/30 where the curves exhibiting higher peaks (in absolute value) correspond to smaller values of δt. As in the previous figure, the expectation value is multiplied by the numerical constant σ s = 2(2π) d−1 Ω d−2 . Further, at each time, the expectation value for an 'adiabatic' quench is subtracted.
Note that the holographic result is even valid for d = 3, where there is no divergence but a linear relation to δt. We leave the detailed analysis of this particular case after we discuss the analytical results in section 2.2.2. d=3 The slope of the linear fit in each case is shown in the brackets beside the labels. The results support the power law scaling φ 2 ren ∼ δt 4−d .

Analytical leading contributions: d ≥ 5
The numerical results above revealed that fast mass quenches in the free scalar theory have the same early time scaling (1.3) as in the holographic quenches [14,15]. However, looking at the curves of figs. 2 and 3, the entire time profile of the expectation value seems to take a relatively simple and possibly universal form. In particular, for odd spacetime dimensions, one can easily verify that the response takes a form similar to a certain time-derivative of the mass profile. In particular, it seems that φ 2 ren ∝ ∂ d−4 t m 2 (t), where the power of the time derivative in this ansatz was chosen as it matches the power-law scalings already discussed. In this section, we will verify this universal form by developing an expansion of the hypergeometric functions which allows us to extract the leading behaviour of the expectation value in the limit in which δt → 0. In fact, we will show that this leading behaviour is in perfect agreement with the numerical response presented in previous subsection. In the case of even d, we perform a similar expansion to again extract the leading universal response for small δt and we will find an enhancement by a logarithmic factor.
We show in a blue solid curve the best fit by a function f (δt) = δt −α (a log δt + b), where we get α = 1.9995 for d = 6 and α = 4.0097 for d = 8. The purple curve is the best fit for a function f (δt) = aδt 4−d . The plots clearly show that there is an extra logarithmic divergence in expectation values. The results support the scaling φ 2 ren ∝ δt 4−d log(δt) for even d.
We first define dimensionless parameters. The relevant physical variables in the quenches here are the initial mass m, the momentum k and the quench rate δt. With those, we define Now we want to expand the hypergeometric function for small κ and fixed q. We will need to expand the hypergeometric series in eq. (2.27) to second order in κ, which gives where the notation ( ) n is as defined below eq. (2.27). Also note that given our definition (2.38), terms with higher powers of κ will contain extra factors of δt and so in the limit of δt → 0, these contributions will be subleading, giving a slower scaling with δt. From eq. (2.39), we see here that each term in the infinite series has an order κ 2 contribution. Indeed the contribution proportional to κ 2 is an infinite series in powers of z. However, we are only interested in computing | 2 F 1 | 2 and then integrating over all momenta.
Remarkably it turns out that for a given d these integrals which multiply factors of z p vanish for all p ≥ p d where p d is an integer which depends on d. Therefore we can calculate the nonvanishing contributions to φ 2 explicitly, with only the first few terms. We also need to regulate the expectation value after making this expansion. So in the same way as before, we expand the series for large q and subtract the divergent contributions. Of course, this procedure produces leading order expansion in κ 2 of the counterterm contributions that were found in eq. (2.26). We are able now to compute the leading contribution in κ 2 to the expectation value of φ 2 . This gives, for odd (2.40) Note that to get this universal result, we need to use the same relations that were used in computing the counterterms in order to relate z with m 2 (t). As we are considering the mass profile m 2 (t) = m 2 2 (1 − tanh t/δt), eq. (2.40) supports the early time scaling that was found numerically above. Here, it emerges from the leading term in an analytical expansion when δt → 0. A nice way to visualize this behaviour is to replot the numerical results as δt d−4 φ 2 ren and compare the curves with the leading order contribution (2.40). This is shown for d = 5 and 7 in fig. 6. As we see, the numerical curves collapse down onto the leading analytical profile as δt gets smaller and smaller. Further, the plots demonstrate that that the numerical curves converge to the leading behaviour (2.40) more quickly in higher dimensions, as might be expected since the power law scaling is more pronounced. Now let's turn to the case of even dimensions where the situation is more subtle. First, we have the IR regulator µ which we use to produce the dimensionless variable ν = µδt along with κ and q, as in eq. (2.38). Now we follow the same procedure as before: expanding to leading order in κ 2 and further expanding for large q to find the counterterm contributions. The difference in this case is that in evaluating φ 2 ren , the integration over the momentum is divided into two regions, as described in subsection 2.1.1, and this is where the ν dependence will appear. In fact, in a manner similar to that found above, we find that the entire ν contribution is encoded in the first few terms of the expansion of hypergeometric functions and after some manipulation, those terms simplify to yield where we already wrote the expectation value in terms of dimensionful µ and the dots indicate terms independent of µ. However, let us note that we will see that the latter contributions include terms that still scale as δt 4−d . Further let us re-iterate the comment in footnote 9 that for the above behaviour contribution to become dominant, we need µδt 1 as well as m δt 1 to be in the fast quench regime. Hence eq. (2.42) reveals a further logarithmic enhancement of the leading response over the power law scaling in eq. (1.3). Rather for even d, we find where we have set µ = 1 above. In fact, this logarithmic enhancement is exactly the kind of behaviour found in the holographic studies [14,15]. If we present the numerical response as δt 2 φ d−4 ren , as is shown in fig. 7(a) for d = 6, the peaks in the curves continue to grow as δt becomes smaller and smaller. This growth reflects the additional logarithmic factor appearing above in eq. (2.43). . Panel (a) shows δt 2 φ 2 ren . As we reduce δt from δt = 1/10 to δt = 1/1000, the peaks in the response continue to grow. In particular, the curves correspond to δt = {1/10, 1/20, 1/50, 1/100, 1/1000}. Panel (b) shows δt 2 φ 2 ren / log δt. Here as δt decreases, the curves converge to the analytic expression (red dashed line). In this case, the amplitude of the left peak increases monotonically as δt shrinks and the various curves correspond to δt = {1/10, 1/20, 1/30, · · · , 1/100, 1/1000}.
In fig. 7, we show instead δt 2 φ d−4 ren / log δt for d = 6. There is also red dashed line that corresponds to the leading order expression derived analytically and as expected, the numerical response collapses down onto this analytical profile as δt decreases. However, there is still a part of this analytic response that we need to describe. Basically, the hypergeometric function will give us a structure like where φ 1 (t) is given by eq. (2.42). Unfortunately, φ 2 (t) cannot be expressed as neatly as in the case of φ 1 (t), possibly indicating that the form of this contribution is not universal. In fact, all of the terms in the expansion (2.39) of the hypergeometric functions contribute to this profile. The result for even dimensions d ≥ 4 can be written We have written the double sum in terms of a limit because we found that we can approximate the entire expression for φ 2 (t) well with the expression above where h is kept finite but large. In particular, the analytical profile shown in fig. 7(b) corresponds to eq. (2.44) evaluated with δt = 10 −3 and taking h = 25 in eq. (2.45), as well as m = µ = 1. Again, as shown in fig. 7(b), there is essentially exact agreement between the numerical solution and this analytic profile. Note also that even for δt = 10 −3 , log(δt) ∼ −6.9 and so both terms in eq. (2.44) contribute significantly to the expectation value, i.e., one must go to much smaller values of δt before φ 2 (t) can be neglected. Finally, let us turn to the question of why we did not see the logarithmic enhancement in the original numerical results, i.e., in figs. 1, 2, 3 and 4. Recall that in those plots, we were examining φ 2 ren (t = 0) as a function of the quench times δt. The key point here is that we choose to evaluate the response at time t = 0. Here we might note that in fig. 7(a), all of the curves go through the same point at precisely t = 0, i.e., the entire scaling has been removed by multiplying by δt d−4 at this time. This effect arises because we are studying the specific mass profile m 2 (t) = m 2 (1 − tanh(t/δt))/2. In this case, any even number derivatives of this profile precisely vanishes at t = 0. Hence we were simply unlucky in our choice of the time at which to sample the response. As shown in figs. 5 and 7(a), the logarithmic enhancement can be seen in the numerical results when we examine the response at any other value of t.

Analytical leading contributions: Low dimensional spacetimes
There are a number of reasons to treat d = 3 and 4 separately. First, eq. (2.40) does not make sense when d = 3 since the latter would give a negative number of time derivatives in this formula. Moreover, for both d = 3 and 4, all terms in the hypergeometric series expansion (2.39) contribute. Finally, our expansion in powers of κ has some problems in the IR related to simultaneously taking the limits κ, q → 0.
Let us illustrate the latter problem with d = 3. In this case, the counterterm contribution (2.26) reduces to f ct (k, m(t)) = 1 and hence in terms of dimensionless variables, eq. (2.11) can be written as However, if we now first expand the integrand in powers of κ and then consider the limit q → 0, we find an extra divergent term: −κ 2 /(2q 2 ). Of course, this ill-behaved term arises because we are expanding q/ω in δt = 1/ 1 + κ 2 /q 2 for both κ and q around zero. For general dimensions, this term becomes q d−2 /ω in δt = q d−3 / 1 + κ 2 /q 2 and the order κ 2 term becomes − κ 2 2 q d−5 . Therefore a similar logarithmic divergence appears for d = 4 but no extra divergence appears at order κ 2 for d ≥ 5. Furthermore, we observe that we do not encounter any IR divergence coming from the same expansion for the reverse quench (i.e., A = B = 1/2) in any d. In the latter quenches, we have simply ω in δt = q. Finally, we note that no such IR divergence appeared in the numerical calculations for either d = 3 or 4. Therefore we conclude that this is not a physical divergence of our system. Rather it is a spurious problem generated by our expansion in powers of κ.
Hence we remove this divergence by simply subtracting the spurious term as an extra counterterm contribution, which yields for d = 3, 47) where σ s was defined in eq. (2.37). Above, the second expression is just the sum that appears when expanding the hypergeometric function and the third one is the result of summing all the terms in the sum. We might also mention that for the reverse quench in d = 3, we find φ 2 ren = π 4σs m 2 δt log 1−tanh t/δt 2 + O(δt 3 ), which is just the negative of the above result.
Let us also say that subtracting that extra counterterm has its effect on the final expression for the expectation value. In fact, by carefully comparing the full numerical integration with the analytic answer we found that they are shifted by a factor √ m 2 . We write m in this way to emphasize that this extra term is non-analytic in m 2 , so in fact what we are finding is that where we have substituted σ s = 4π for d = 3 using eq. (2.37). We can recognize, though, that this extra term is due to the κ-expansion because, for instance, it does not appear in the reverse quench where ω 2 in = k 2 + m 2 and there is no problem in taking both limits. This difference is illustrated for both types of quenches in figs. 8(b) and 8(c). In section 2.7, we will see that this constant term simply corresponds to the renormalized expectation value for a fixed mass. The red curves correspond to the leading order analytic expression (2.48) while the blue curves are the full numerical solution. By comparing panels (a) and (b), we can see that the difference between the two solutions is roughly of order O(δt). We can observe that apart from having an extra minus sign difference, the reverse quench in panel (c) starts from zero without needing to be shifted by the factor of m. Now if the leading term in eq. (2.47) is evaluated in the middle of the mass quench, we have φ 2 ren (t = 0) = log 2 16 m 2 δt. Hence as observed in fig. 4, this result is linear in δt and so it actually approaches zero in the limit δt → 0. This behaviour should be contrasted with the growing response found in higher dimensions, e.g., as shown in eq. (2.41). In fact, the same diminishing response will be found in d = 3 when the expectation value is evaluated for any finite value of t/δt. However, this scaling is deceptive as it may lead one to expect that the quench has a vanishing effect in the limit δt → 0. Considering eq. (2.47) but in the limit t/δt 1 instead, we find which is independent of the quench rate! We will analyze the late time behaviour of the quench in greater detail in section 2.7. However, the above expression (2.49) clearly indicates that the apparent scaling behaviour shown in fig. 5 and eq. (2.47) does not give an accurate characterization of the overall effect of the mass quenches in three dimensions. Pushing our numerical results to longer times, we were able to go as far as t/δt ∼ 18. At these 'late' times, we found that the response is indeed linear and independent of δt. For instance, a linear fit to the numerical results in fig. (8) certainly respects the analytic limit. 11 However, in section 2.7, we will show that for very late times, where m 2 t 2 1, the growth is no longer linear but rather logarithmic.
This result also highlights another key difference between eq. (2.47) and the leading behaviour (2.40) found in higher dimensions. In higher dimensions, the time profile of the leading analytic term approaches zero exponentially fast (with a 'tanh' mass profile) for t/δt 1, while in eq. (2.47), the corresponding time profile grows without bound at large times.
The situation for d = 4 is quite similar, but now the κ expansion generates an extra logarithmic divergence in the calculation of the response, as already commented above. However, the same discussion as in the case of d = 3 still applies. Being in an even number of dimensions, the leading order response has two components as in eq. (2.44) and so for d = 4, we have Here we find φ 1 (t) = m 2 4 (1 + tanh(t/δt)) and φ 2 (t) is given by eq. (2.45) with d = 4. We also note that for the reverse quench in d = 4, the only difference is that we find φ 1 (t) = − m 2 4 (1+tanh(t/δt)). Evaluated at t = 0, the leading contribution for small δt is just φ 2 ren = m 2 4 log(µδt). Hence as the numerical results in fig. 4 showed, the leading contribution in d = 4 scales logarithmically when δt → 0. Again, this logarithmic scaling was also agrees with the behaviour found in holographic quenches [14,15]. We might also comment that for large times, i.e., t/δt 1, both φ 1 (t) and φ 2 (t) approach a constant in eq. (2.50).

The stress-energy tensor
There is an elegant and independent consistency check of our results involving the energy density. In particular, we can consider the diffeomorphism Ward identity [14] ∂ where E is the (renormalized) energy density. In the case of a constant mass, this identity simply expresses the conservation of energy for the system, i.e., the RHS vanishes identically. But in the case of a time-dependent coupling, eq. (2.51) determines the work done by the quench. Of course, again in our case, λ = m 2 (t) and O ∆ = φ 2 , so with our previous analysis, we already have all of the information needed to compute the RHS of the identity. The independent consistency check will then consist of evaluating the time derivative of the energy density, i.e., computing the LHS of eq. (2.51) directly. The energy density, defined by the T tt component of the stress-energy tensor, is given by where the index i is summed over the spatial dimensions. Given our mode expansion (2.7), this expression results the following expectation value, Now it is straightforward to check analytically that taking the time derivative of the above expression and simplifying the result with the equations of motion for the scalar field, yields exactly 1 2 ∂ t m 2 (t) φ 2 , as required by eq. (2.51). We can also verify this agreement numerically. For simplicity, we will set d = 5 and in this case, we know that all counterterms come from the zeroth order terms in the adiabatic expansion, which can be extracted from the constant mass expectation value -see discussion below. In this case, eq. (2.53) reduces to (2.54) Hence we know the necessary counterterms in d = 5 will to regulate the expectation value of the energy density by subtracting off these first three terms, with m 2 replaced by m(t) 2 . With this subtraction, we can evaluate the finite part of eq. (2.53) to get fig.  (9(a)). By numerically differentiating it with respect to time we should get exactly the RHS of eq. (2.51) and that is indeed the result, as shown in fig. 9(b). As a further check of our analysis, we can verify that the counterterm contributions have the expected form. That is, even though we are finding them separately and independently in eqs. (2.26) and (2.54), they should actually come from the same counterterm action, as discussed in section 2.1.3. In the present case of d = 5, the action in eq. (2.32) reduces to five terms S ct (m 2 , g µν , Λ) = − d d x √ −g s 00 Λ 5 + s 10 m 2 Λ 3 + s 20 m 4 Λ + R s 50 Λ 3 + s 51 m 2 Λ (2.55) The counterterm contributions to φ 2 and E are then determined from this action by eqs. (2.33) and (2.34), respectively. It is clear that the terms involving the Ricci scalar do not contribute to φ 2 when the latter is evaluated in flat space. Similarly, the variation of the s 50 term to E , coming from the variation with respect to the metric, vanishes in flat space. Finally, the variation of the s 51 term yields a contribution of the form: 12 T µν ∼ (∂ µ ∂ ν − η µν ) m 2 . However, since the mass only depends on time, one finds that this particular contribution vanishes for the energy density, E = T tt . Hence, in fact, only the first three counterterms in eq. (2.55) will contribute in the present case. That is, we should find Now if we integrate eq. (2.54) up to a momentum k max and compare to the analogous result in eq. (2.13), we find (2.57) Hence we find that the coefficients of the cubic and linear divergences match between the two expectation values, as desired . Further, we can supplement the list of coefficients in eq. (2.36) with s 00 = −1/(d σ s ), after generalizing eq. (2.54) to d dimensions.

The energy density in three dimensions
As in the case of section 2.2.3, it is interesting to repeat the above analysis but focusing on the d = 3 case separately. In this case, the scaling found for φ 2 ren in section 2.2 is proportional to δt. Hence on the RHS of the identity (2.51), this is multiplied by ∂ t m 2 which gives a factor of 1/δt and so one would find that ∂ t E ren does not scale at all with δt. Since the quench essentially takes place over an interval δt, this would then reproduce the naïve scaling δ E ren ∼ m 2 δt as suggested by eq. (1.2), i.e., no work is done in the limit δt → 0. However, we will show below that this is not really the case and rather we find that ∂ t E ren scales as 1/δt and that δ E ren ∼ m 3 -see figures 10 and 11. Note that the latter result indicates that the work done is not analytic in the mass coupling, i.e., δ E ren ∼ (m 2 ) 3/2 .
Let us start by computing the expectation value of the energy density for a constant mass. In this case and with d = 3, eq. (2.53) yields (2.58) The first two divergent contributions would be removed by the counterterm contributions and hence the renormalized expectation value of the energy density would be E ren = −m 3 /(12π) in the case of a constant mass. Again it is notable that this result is not analytic in the mass coupling. However, we can easily extract the counterterm contributions to regulate the expectation value in the case of a time-varying mass from eq. (2.58) to find f d=3 ct = k 2 + m 2 (t)/2. Subtracting these terms in the integral in eq. (2.53) with d = 3 then yields the renormalized expectation value of the energy density. Then we computed this expectation value numerically for different values of δt ranging from δt = 1/10 to δt = 1/100, as shown in fig. (10). We observe that the energy density grows from its corresponding value at minus infinity -as we set m = 1, this means 4π E ren (t 0) = −1/3 -to a certain constant value at late times. In particular, as δt becomes smaller, the latter constant seems to be independent of δt. Hence, from this figure, we can see that the naïve power counting does not work, because as described above, it suggests that the change in energy density would be proportional to δt. Further, we can also compute the time derivative of this profile and compare it with the RHS of the Ward identity (2.51), using our previous results for the expectation value of φ 2 . Again, we get perfect agreement, as shown in fig. 11. There we also see that ∂ t E scales as 1/δt.
How can we understand this scaling?
The key point is that the change in the expectation value of φ 2 has a scaling proportional to δt but the full expectation value does not start from zero. Recall from eq. (2.48) that the expectation for d = 3 is given Now, as δt → 0, the second term and all the subleading will go to zero and then φ 2 ren −m/(4π). So if we integrate the Ward identity (2.51) in this limit, we find Further at very early times (i.e., t 0), the energy density will match that found in the case of a constant mass. Hence given the results in eq. (2.58), we have E t=−∞ = −m 3 /(12π). Hence for late times (and small δt), we should find the energy density to approach E t=∞ = m 3 /(24π), which is exactly what is shown for the long time behaviour in fig. 10.
To close this section, we reiterate that eq. (1.2) suggests the scaling of the energy should be δ E ren ∼ m 2 δt for d = 3. This scaling was not realized here in eq. (2.60) but this result depended on the fact that φ 2 = −m/(4π) in the past, i.e., at the start of the quench. On the other hand, if we considered a 'reverse' quench, where the mass starts at zero and rises to some finite m, this initial expectation value would vanish and hence the expected scaling would be fulfilled. That is, zero work is done by the reverse quench in the limit δt → 0.

Universal scaling of higher spin currents
It is known that free scalar field theory has an infinite set of higher spin conserved currents j i 1 ···is [21,22]. Apart from being conserved, these currents are symmetric in their indices and, in the case of massless theory, traceless. It is interesting, then, to analyze how these currents behave in the present quenches. In particular, we will be interested in determining how the higher spin currents scale in the fast quench limit.
Higher spin currents for a massless complex scalar field are given by [22] (2.61) where indices i 1 , · · · , i s should be symmetrized above. In case of a complex scalar field, the even spin currents are symmetric under the interchange φ ↔ φ * , while odd spin currents are antisymmetric. In our calculations, we are dealing with real fields and so the odd spin currents trivially vanish. Hence we will only consider the even spin currents.
Let us start by revisiting the spin-2 current, i.e., the stress-energy tensor of the conformally coupled scalar. Hence we can obtain this current by varying the scalar field action with respect to the metric, We note that the equation of motion, φ = m 2 φ, was used to simplify the above expression. Further, we can verify that if we set m 2 = 0, the above result reproduces the s = 2 current in eq. (2.61), up to an overall numerical factor. It will be convenient for the following to split the current into two parts: the minimally coupled current (obtained by setting ξ = 0) and the remaining contribution coming from the conformal coupling term proportional to R in eq. (2.63). Then we have Of course, we have restored the terms involving φ using the equations of motion in j (2)conf ab . The reason for doing so is that it makes apparent that j (2)conf ab is a total derivative, i.e., j Then, in our case (where φ 2 only depends on time), we find that the a = b = t component of this part vanishes and we are only left with the minimally coupled current. Therefore the energy density calculated with the full stress tensor (2.64) agrees with that found with the minimal stress tensor (2.66), as was done in the previous sections. Of course, for a constant mass, the spin-two current (2.64) is conserved. However, if we allow for a time-varying mass, the divergence of this current yields Of course, we have reproduced the diffeomorphism Ward identity (2.51), from which we can determine the energy which the quench injects into the system if we are given the expectation value φ * φ . The reason for revisiting this result for the spin-2 current is that we will now apply the analogous analysis with the spin-4 current and we will find the scaling of this higher spin current in the limit of fast quenches. Further, we will use this approach to argue for the scaling of all of the higher even spin currents. First we must build the spin-4 current for the massive theory as follows: Take eq. (2.61) and explicitly symmetrize the indices. Then introduce all the necessary trace terms with the necessary coefficients to ensure that the result is traceless in the massless case. The next step is to generalize this current for a massive field. Here, we take the divergence of the massless expression and add all the necessary terms proportional to the mass to ensure that the divergence vanishes upon evaluation on the massive equation of motion. This procedure is explicitly carried out for the spin-4 current in Appendix A. The final result is j (4) where again the indices in last term should be symmetrized. Now we are interested in obtaining the analogous Ward identity for the spin-4 current. In particular, we make the mass time-dependent and evaluate the time-derivative of the j tttt component. Note that the part proportional to the conformally coupled spin-2 current will vanish and hence we find To determine the scaling of this spin-4 'charge density' in the limit of fast quenches, we can use the scaling of the energy density E ∼ m 4 /δt d−4 to find: Hence in the fast quench limit, the spin-4 charge diverges with precisely the same power of δt as the spin-2 charge and the spin-0 charge (i.e., φ 2 ), while an extra power of m 2 appears to make up the necessary dimension of the new operator. Extending the construction of the spin-4 current, described above, to obtain higher spin currents in the massive theory is straightforward, though tedious. We expect that the massive terms for the spin-s current can decomposed, as in the spin-4 case, in terms of the spin-(s-2) current and a total derivative term. Then, it is easy to see that for a time-varying mass, we will get a hierarchy of generalized Ward identities, (2.74) Now integrating these identities will similarly yield a hierarchy of scalings for the final currents in the fast quench limit, i.e., j . Hence the scaling of all of the higher spin currents would be determined by that originally found from how φ 2 scales. Then in general we should find that (2.75) Of course it would be interesting to explicitly construct the currents in the massive theory and derive these scalings for the higher spin currents. However, our expectation is that after a quench, all of currents that will scale with precisely the same power of δt. In particular then, for d ≥ 4, all of these currents will diverge as δt → 0.

CFT to CFT quenches
This subsection is devoted to study the response of the scalar field under a quench whose mass profile is asymptotically zero at both infinite past and future. We smoothly turn on the mass up to some m 2 and then go back to the critical point. The whole process is again characterized by a time length δt. We may proceed analytically if we choose the following mass function m 2 (t) = m 2 cosh 2 (t/δt) . (2.76) Analysing this system is interesting because it provides a further check of our previous analysis. In particular, we should expect to have the same scaling behaviour for the renormalized expectation values in the limit of fast quenches. Moreover, the counterterms should be the same as in the previous case with the only difference that we should change the mass function (and its derivatives) to the new profile. Even though this is expected, it is not at all trivial : rather it provides a good confirmation of our results. Lastly, we will return to such CFT-to-CFT quenches later in section 4 to give a general argument that should be valid for arbitrary CFTs and hence the present section provides an explicit example of these processes. Finally, as in the case of the tanh profile, it is straightforward to extend the present analysis of these pulse-like quenches to include a constant mass, i.e., We will explicitly analyze quenches with this profile in section 2.4. However, our intuition suggests that universal scaling in eqs. (1.2) and (1.3) should still hold if we satisfy both m 2 δt 2 1 and m 2 0 δt 2 1. Again, this emphasizes that what is important is that the theory has a UV fixed point, i.e., , the UV description of the theory is a CFT. The IR details become unimportant in the fast quench limit, i.e., when 1/δt dominates all of the IR scales.
As with the tanh profile (2.1), we can exactly solve this problem by decomposing the scalar field into momentum modes, as in eq. (2.7). This modes will satisfy the Klein-Gordon equation with mass given in eq. (2.76), This equation can be written in hypergeometric form by expressing it in terms of variable y = cosh 2 (t/δt), We are interested in those solutions that behave purely as positive frequency waves in the infinite past, so we need to fix initial conditions so that u where ω k is just ω k = k because in the infinite past we are in the massless theory. 13 Then, the complete solution for the modes in terms of k and y is given by Now, as we did in the previous case, we integrate over momentum modes in evaluating the expectation value of φ 2 and this integral is UV divergent. To get the finite renormalized expectation value, we must subtract the appropriate counterterm contributions. In section 2.1.1, we used an adiabatic expansion to obtain the counterterms supposing only that the mass depends on time. Hence we can expect the counterterm contributions in eq. (2.26) will regulate φ 2 for any mass profile. Hence, we use the same expression here and only change the profile of m 2 (t) to the pulsed one (2.76). In this way, we obtain which is UV-finite, as we will see below. A more nontrivial result is that these expectation values should yield the same leading order behaviour, as derived in section 2.2.2, where the results were expressed in terms of derivatives of the mass profile. In fact, we found that this same universal behaviour indeed emerges for the pulsed profile and so eq. (2.40) also gives the correct result in this example. In particular, fig. 12 13 The way to take this limit is to use identities that relate hypergeometric functions of argument z with a linear combination of hypergeometric functions of argument 1/z -see, for instance, [23]. In our case, as t → −∞, 1 − y → −∞ and then, such identities are useful.
shows the renormalized expectation value of φ 2 for odd dimensions d = 5, 7, 9. As we did in the original quenches, here, we divide by the expected scaling m 2 /δt d−4 and plot eq. (2.82) for different time intervals δt. We see that the curves rapidly converge to the analytic expression given in eq. (2.40) as δt goes to zero. This clearly shows that both the expected scaling in eq. (1.2) and the leading analytical behaviour in eq. (2.40) are valid in the present example of a pulsed quench. The solid curves correspond to δt = 1, 1/2, · · · , 1/10 with δt decreasing as they converge to the analytical leading expression (2.40), plotted with dashed red curve. This leading term has φ 2 (d) For even d, we expect the scaling to be enhanced by a logarithmic factor, as discussed in section 2.2.2. In the case of the previous case with the tanh profiles, we could not see this enhancement in our numerical results [16] because the leading order term vanishes at t = 0. However, in the present case, the even derivatives of the mass are not zero at t = 0 and hence, we should be able to see the expected behaviour even at zero time. This can be seen exactly in fig. 13, where the fits of the curves support the scaling φ 2 ren ∼ m 2 /δt d−4 log δt. In contrast, for odd d, the corresponding derivatives of pulse profile (2.76) vanish at t = 0. However, we can instead evaluate φ 2 ren (t = −δt/2) to reveal the same scaling applies in odd d, as shown in fig. (14). Of course, this scaling was already confirmed above by matching the leading analytic behaviour.

Universal scaling for arbitrary initial and final mass
In this section, we would like to show that the universal scaling in eq. (1.2) is not exclusive to quenches which involve a critical theory at the initial and/or final times, but are also found for arbitrary initial and final mass under certain assumptions. Basically what we need is 1/δt to be the only relevant scale of the problem. So long the initial and final mass (and their difference) are much smaller than 1/δt, we will find the same scaling.
For the tanh profile (2.1), this scaling can be explicitly seen by extending the analysis of the renormalized expectation value (2.11) to general initial and final masses, i.e., general A and B in eq. (2.1). Hence we have where (2.85) and the counterterm contributions f ct (k, m(t)) are given by eq. (2.26). Now let us redefine the integration variable in eq. (2.83). We definẽ and hence With this choice, eq. (2.83) starts to look like the expectation value for a quench from an initial mass-squared (δm 2 ) to the massless case. In fact, the absolute value of the hypergeometric function in the integrand will look exactly like that. We have to take care about the rest of the integral. Applying the change of variables (2.86), eq. (2.83) becomes Further as in section 2.2.2, we introduce a dimensionless momentumq =kδt, which then yields (2.90) In the limit of m f δt 1, the expectation value becomes , (2.91) up to contributions suppressed by m 2 f δt 2 . Hence we have reproduced exactly the expression computing the renormalized expectation value for a quench starting at (δm 2 ) and ending at zero mass. Then, as shown previously, in the case where δ(m 2 )δt 2 1, the expectation value of φ 2 scales as δm 2 δt 4−d .
Hence to obtain the universal scaling (1.2) in quenches with arbitrary masses, we need to satisfy two conditions: It is easy to check that these two conditions are equivalent to those in eqs. (2.3) and (2.4), i.e., (m 2 i − m 2 f )δt 2 1 and (m 2 i + m 2 f )δt 2 1. Finally, we will comment on the case of the pulsed quench around any arbitrary mass. If our mass profile becomes m(t) 2 = m 2 0 + m 2 cosh 2 t/δt , then it is easy to verify that the only change in eq. (2.79) is to add a term proportional to m 2 0 ending up with In analogy to eq. (2.86), we definek 2 = k 2 + m 2 0 , so that the equation becomes the same but with k →k. Then the solution for the modes will be the same with the only difference that in eq. (2.81), k is replaced withk (in a and b). To obtain the expectation value, we will have to integrate over all momenta. In a way completely analogous to the previous case, we can perform a change of variables to integrate ink and in the limit of m 2 0 δt 2 1, we will get exactly the same integral as in section 2.3. Hence it is expected that the same scaling will appear. In conclusion, for the pulsed quench, the expectation value of φ 2 will scale as m 2 δt 4−d provided that m 2 δt 2 1 and m 2 0 δt 2 1.

Comparison to linear response
The results of sections 2.1 and 2.4 for the tanh profile are leading order in the dimensionless variable (mδt) 2 . Therefore they should agree with a linear response calculation.
In this section, we compute 0|φ 2 |0 in in linear response theory for the quench starting from a CFT and ending with a massive theory and show that the result is in exact agreement with the expansion of the exact answer to O(m 2 ), for each of the k modes individually. This agreement should hold for the other kinds of protocols as well, such as the pulse profile in section 2.3. The linear response result for the expectation value 0|φ 2 |0 in is given by the expression where the retarded correlator is given by The correlation functions are to be evaluated in the initial theory, which is the massless free field theory. The right hand side can be computed exactly leading to (2.97) We will express the right hand side of eq. (2.98) as a power series expansion in η = exp(2t/δt) . (2.98) In eq. (2.97), we write m 2 (t ) = m 2 2 (1 + tanh(t /δt)) = m 2 η 1+η . Then expanding this expression as a series in η and performing the intergral over η , we obtain Let us now consider the O(m 2 ) contribution to 0|φ 2 (x, t)|0 in from the exact answer. This is given by We need to expand the hypergeometric function to O(m 2 ) and express the answer as a power series expansion in η. It turns out that Substituting eq. (2.102) into eq. (2.100), it is easily seen that the O(m 2 ) contribution to the exact answer matches the answer from linear response theory (2.99).

Comparison with instantaneous quenches
The results of the previous sections appear to be at odds with the well studied examples of instantaneous (or abrupt) quenches in field theories, in particular [5,6]. The behavior of e.g., eq. (1.3) suggests that for ∆ > d/2, the expectation value of the operator O and hence the rate of energy production diverges in the limit δt → 0. In contrast, the results of instantaneous quenches indicate that there is a smooth limit. In this section we resolve this apparent discrepancy.
The main point is that the fast quench limit, considered here, involves a quench rate, i.e., 1/δt, which is fast compared to the scale set by the relevant coupling, but slow compared to the UV cutoff. This is implicit in the above since we are working with renormalized quantities, where in fact the UV cutoff has been sent to infinity. However the abrupt quenches which are considered in the literature involve an instantaneous change of the Hamiltonian at some time, e.g. t = 0. The wave function evolves from early times according to one time independent Hamiltonian H in up to time t = 0. The resulting wavefunction at t = 0 then acts as an initial condition for evolution with a different time independent Hamiltonian H out . This process can be considered as a limit of a smooth time dependent Hamiltonian provided the scale of variation is infintely fast compared to all scales in the problem. In a field theory, this means that 1/δt is large compared to all momentum scales including the UV cutoff scale Λ. This is clearly not the limit considered in our work.
To make this point explicit, we will now compute the two-point correlation function in position space in the free bosonic field theory with the time dependent mass given by eq. (2.1). We are interested in this at late times. For this purpose, it is convenient to work in terms of the "out" modes, with the various frequencies defined in eq. (2.9). These modes have the usual planewave behaviour at late times, t δt, (2.105) The in and out sets of modes (u k and v k , respectively) are related by a Bogoliubov transformation The Bogoliubov coefficients have been evaluated in [18], The correlation function of the field is then given by Using (2.107) one finds (2.109) Consider now the limit 2.110) in which the quantities appearing in eq. (2.109) become Let us now compute the correlation function (2.108) taking both the limit ( Note that δt has disappeared from the result. In fact this reproduces the result for an instantaneous quench from a mass m 2 in = m 2 (A − B) to a mass m 2 out = m 2 (A + B), e.g., see eq. (8) of [6].
In this paper, we have concentrated on local quantities like φ 2 or the energy density. These involve integrals over all momenta all the way to the cutoff, and clearly the limit (2.110) is not appropriate for large UV momenta in these integrals. In our analysis, we have worked with renormalized quantities which, as we noted above, implicitly involves taking the UV cutoff Λ much larger than 1/δt. This is why our limit of fast quenches is physically different from the instantaneous quenches, studies elsewhere, where the quench rate is necessarily fast compared to Λ. In fact in the continuum limit, it is unphysical to consider such an instantaneous quench. It would be interesting to investigate these issues in a theory with finite cutoff. In such a theory one would expect that the scaling discussed in this paper should hold in a protocol where Λ −1 δt m −1 . On the other hand, when δt is the same order as Λ −1 , the answers should approach those for an instantaneous quench.
Nevertheless, for distances | x − x | δt, only momenta much less than δt −1 should be making a substantial contribution to the correlation functions. For such quantities, the condition (2.110) is effectively satisfied and so by the above analysis, one should expect only small differences between a fast smooth quench and an instantaneous quench at late times, i.e., when eq. (2.112) is also satisfied. Details of this comparison will be discussed in [17] -see also discussion in the following section and in section 5.

Late time behaviour
In section 2.2.3, we observed some interesting late time behaviour for the expectation value φ 2 (x, t) in three dimensions, i.e., at late times, the expectation value is independent of δt. This may lead us to suspect that this late time behavior agrees with the results of an instantaneous quench. In this section, we will show that in a suitable regime this is indeed true for d = 3, but one finds that the same agreement does not generally occur in higher dimensions [17].
As we noted in section 2.2.3, the numerical analysis only allowed us to evaluate the expectation values out to times of order t ∼ 10 δt. Hence given the values of m and δt that we were using, we were always in a regime where m 2 t 2 1. We therefore first compare the result obtained in section 2.2.3 with the result of an instantaneous quench in this regime. We will show that with the results already obtained we can reproduce our previous results in this limit but also go beyond them and evaluate the proper long time behaviour of the scalar field for m 2 t 2 1. The starting point will be to consider the correlator for instantaneous quench, eq. (2.113) and evaluate this expression at coincident points in space and time, i.e., x = x , t = t . This gives, Focusing on the quench to the critical point, i.e., A = −B = 1/2 in eq. (2.1), for which ω 2 out = k 2 and ω 2 in = k 2 + m 2 , we find Of course, this expectation value is divergent in the UV, so it must be regulated as described in section 2.1.1. While in general this is a somewhat involved procedure, we begin here by considering d = 3 in which case there is a single mass-independent UV divergence -see eqs. (2.13) and (2.14). Hence the difference between the quenched expectation value and that for a fixed mass m will produce a finite result. 14 That is, we subtract from eq. (2.115) and then evaluate the finite difference In these expressions, we have substituted σ s = 4π for d = 3 using eq. (2.37). Above in the second line, we also introduced the dimensionless momentum p = k t.
The first thing to verify is that we recover our previous results for d = 3 in the regime where m 2 t 2 1 -see discussion in section 2.2.3. In this limit, we can drop the m 2 t 2 appearing in the denominator of the integrand to find This is exactly the same result we found in eq. (2.49), showing a linear growth in the expectation value of φ 2 with a slope that is independent of δt. From these results, we can also identify the constant displacement in eq. (2.48) and in fig. 8 as the renormalized expectation value for a constant mass, i.e., φ 2 f ixed,ren = −m/(4π). However, given eq. (2.117), we can go further and analyze the behaviour of the expectation value for any value of m 2 t 2 . In particular, this expression can be integrated exactly for any m 2 t 2 and evaluated in terms of generalized hypergeometric functions, (2.119) Fig. 15 shows a plot of this expectation value as a function of mt. From the figure, we observe that the linear growth (2.118) of the expectation value is only valid for mt 1. After that, the expectation value continues to grow but in a slower rate. In fact, one can take the limit mt → ∞ in eq. (2.119) to find Hence we see that for very late times, i.e., mt 1, the expectation value continues to grow but only logarithmically. In any event, if we look into infinite future time, the expectation value is divergent. At first sight, this unbounded growth may seem counterintuitive since, for example, it may seem that the physical work done by the quench will also diverge. However, if we recall that eq. (2.51), the time rate of change of the energy density is given by the product of this expectation value with the derivative of the mass coupling. For the mass profile (2.1), the latter decays exponentially in time and hence the corresponding integral for the energy density remains finite and well-defined, despite the logarithmic growth of the expectation value (2.120).
Let us now examine the question: why does the long time answer for smooth fast quench as defined in this paper agree with the instantaneous quench result for mt 1. We need to consider the validity of assumptions which were implicit in the above discussion. Our starting point was eq. (2.114) which was found by taking the limit of coincident points in eq. (2.113). However, the latter correlator was simplified by assuming late times as in eq. (2.112) but also small δt as in eq. (2.110). While the late time assumption is certainly valid here, it is not clear that the second assumption should hold. In particular, one expects that for sufficiently large momenta that the inequalities in eq. (2.110) will be violated. However, if we examine the form of the integrand in eq. (2.117), we see that it decays as roughly 1/p 2 for large (dimensionless) momentum. Hence we can expect that the dominant contributions to the integral come from small and finite values of p. Further given that p = kt, we will certainly satisfy kδt 1 in the late time limit and hence eq. (2.110) will be satisfied. For example, one can make a simple estimate of the error introduced in ignoring the very high momenta as follows: Certainly, eq. (2.110) is satisfied for k ∼ m and hence the integrand in eq. (2.117) is accurate for dimensionless momenta at least up to p = mt. Then an upper bound on the error in our result is given by the integral from p = mt to ∞ but removing the factor of sin 2 p. The final result of this integration is a fixed constant, i.e., approximately 0.07 m. Hence at large mt, this upper bound on the error is small compared to the results given in eqs. (2.118) and (2.120). In fact, given that the numerical fits to the constant term were also good, this suggests that even this approximation is a gross over-estimate of the error.
In fact we can check the validity of our approximation by comparing the full integrand of eq. (2.108) in the limit of late times, i.e., by using eq. (2.105) for the out-modes, with the approximate eq. (2.115). Let's recall that eq. (2.108) is not assuming any relation between the energies in the system and δt, while eq. (2.115) assumes ωδt 1 for every ω. Essentially, we want to compare the integrands of where α k and β k are given by eq. (2.107), and (2.122) for d = 3. Fig. 16(a) plots the integrands in these two expressions as a function of k and in fact, there is no visible difference. The figure uses mδt = 10 −3 and mt = 10 but similar results hold for different values of these parameters. It is clear that the integrand decays rapidly, i.e., in fig. 16(a), it has become negligibly small around mk ∼ 5. Therefore the approximation kδt 1 is effectively satisfied since even though we are integrating over all momenta in the expressions above, the main contribution comes from very low momenta. The latter is explicitly verified in fig. 16(b) which shows φ 2 after both the smooth and the instantaneous quench. As we show before, the instantaneous expression can be integrated analytically and the final result is given by the right-hand side of eq. (2.119) plus φ 2 f ixed,ren = −m/(4π). This result is shown in the figure with the solid blue curve. The purple points correspond to integrating numerically eq. (2.121). We see good agreement between both expectation values. In fact, if we compute the relative difference between them at late times, we see that it is of order 10 −6 and hence we verify that both approaches give the same result at late times.  However, the same agreement does not hold in higher dimensions, as we will discuss in detail in [17]. However, let us present the late-time limit of the smooth quench in d = 5 here. Recall that the desired expectation value is given by eq. (2.121) with d = 5. Although this expression is quite complicated, we can integrate it numerically for different values of mt and follow the evolution of the expectation value at late times, as shown in fig. 17 for mδt = 10 −1 . 15 In the figure, the blue dots are obtained by evaluating the absolute value of the hypergeometric, as we did in section 2.2. However, that analysis only allowed us to go relative short times, in units of 1/m. The evaluation of eq. (2.121) is shown in purple dots for late times and we can see a nice continuity between the two approaches, showing its consistency. The exponential fit of the purple dots also shows that the expectation value is decaying as expected due to the exponential nature of the mass profile. Note that even though the decay is exponential, it does not decay to zero, but to a finite value. One can perform the analysis for different δt's and see that in the limit of δt → 0, that constant approaches to φ 2 (t → ∞) 0.168 m 3 /σ s . It would be interesting to have an analytical understanding of this asymptotic value and also to generalize these results to higher dimensions.

Quenching a free fermionic field
Another way to test our universal scaling formulae in eqs. (1.2) and (1.3) is to quench an operator with a different conformal dimension. In this section, we will be quenching the mass of a free Dirac fermion ψ in d dimensional spacetimes. Then, our operator of interest will be O ∆ = ψ ψ , whose conformal dimension is ∆ = d − 1 and the corresponding coupling is the mass λ(t) = m(t). It is interesting that in this case we should expect divergences to appear as O ∆ ren ∼ δλ/δt 2∆−d = δλ/δt d−2 and hence, even for low dimensional spacetimes with d = 2, 3 we should be able to find divergent behaviours. The calculations are analogous to those for the scalar field. In partiuclar, the situation can be related to one of fermions in curved space-times, where analytic solutions are known for specific mass profiles. Then one can compute the expectation values and find numerical solutions. We will also be able to find analytical leading order solutions in the fast quench limit when δt → 0. The main conclusion is that the scaling relations, (1.2) and (1.3), which were originally discovered by the holographic analysis also hold in this case. Further, the universal power-law scaling is enhanced by a logarithmic factor in the case of even d.
We will be following the conventions and notation of [20], where the problem of fermions in flat FRW backgrounds is discussed. In this case, the equations of motion for a Dirac field Ψ are not directly those that we are interested in, i.e., the Dirac equation with a time-dependant mass. However, it is possible to do a confomal rescaling of the fields, Ψ = C(t) 1−d 2 ψ, where C(t) 2 is the expansion factor, and then, one finds that ψ satisfies the Dirac equation of motion, (iγ µ ∂ µ − m C(t)) ψ = 0, (3.1) where we will define our time-dependent mass as m(t) = m C(t). The exactly solvable model here requires C(t) = A + B tanh t/δt [20], and so m(t) = m (A + B tanh t/δt) , (3.2) in contrast to the scalar case (2.1), where it was the mass squared that had the tanh profile. Solutions to eq. (3.1) are given by where j denotes spatial coordinates and φ k is a scalar field that satisfies For simplicity in the last expression we are not writing the time dependence on C or the fields any more. Now, the full solution for fermionic field ψ can be written in terms of the in modes as (3.6) and the sum over the spinor index λ runs up to 2 (d−3)/2 if d is odd. Here u(0, λ) and v(0, λ) are constant basis spinors, the ω's here and below are defined as in eq. (2.9) and m in = m(t = −∞). Further, a in and b in are operators that annihilate the in-vacuum. It can be shown that this solution reproduces the corresponding solutions for flat space at infinite past and infinite future -see [20]. For the tanh mass profile (3.2), there exist analytic solutions for the scalar field φ k that are of the form where 2 F 1 is the usual hypergeometric function. Note that this solution is similar but not equal to that appearing for the scalar field modes (2.8). Further, we will again focus on quenches to the critical point, where A = −B = 1/2. Now we are interesting in finding the time evolution of the mass operatorψψ through the quench. This is given by that after some algebra it turns to where φ k is actually φ in(+) k and σ f is a numerical factor that depends on the spacetime dimension as As in the case of scalar fields, this expectation value is in general UV divergent so we need to regulate the result by subtracting the appropriate counterterm contributions ψ ψ ren ≡ σ −1 f dk(ψ div (k) − f ct (m(t), k)) . (3.11) These counterterm contributions can again be found as for the scalar field in section 2.1.1. In this way, we find that are all the necessary terms needed to regulate theories up to d = 7. Note that again contributions with time derivatives of the mass profile appear, now from d = 4 onwards. Further, the first line of eq. (3.12) corresponds to the counterterms that would appear in order to regulate the expectation value for a constant mass. Given this finite expectation value (3.11), we are able to evaluate it numerically for all dimensions and different values of the quench rate δt. The results are shown in fig. 18. Note that in these plots, we are subtracting the expectation value in the adiabatic case, for which we are using δt = 10. We verified that the adiabatic expectation value is independent of δt as long as δt is large enough.
As in the scalar case, we should distinguish between odd and even spacetime dimensions. In the case of even d, we also get logarithmic divergences (apart from the usual power-law divergences), that need extra renormalization scales in order to avoid infinities near k = 0. Much in parallel to the scalar case, this will generate a logarithmic enhancement of the scaling behaviour in the expectation value.
We can appreciate how the expectation values grow as we decrease δt in different dimensions in fig. 18. Note that in contrast with the scalar case, now we can see large growth appears as low as d = 2. In order to quantify the precise nature of this growth, we compute the expectation value at a fixed time t = 0 for a larger range of δt and plot it in a log-log scale in fig. 19. Even though the choice t = 0 appears not to be appropriate to find the expected logarithmic enhancement, the figure and the linear fits there support completely the expected power-law scaling ψ ψ ren ∼ m/δt d−2 . Again, if we do the same exercise but at a slightly shifted time, we find that in even dimensions there is a logarithmic enhancement of the divergences. This behaviour will be supported soon by analytical results in computing the expectation value. For now, we can only say that there is a logarithmic growth in d = 2 that is in agreement with previous holographic results.
As in the case of scalars, one can recognize certain relationship between the expectation values and time-derivatives of the mass by looking at the plots of fig. 18. What we will show next is that if we compute the leading contribution in the limit of δt → 0, we'll find precisely those mass derivatives. The different curves correspond to δt = 1/1, 1/2, · · · , 1/10. The curves are so that higher peaks (in absolute value) correspond to smaller δt. Note also that we are plotting the expectation value multiplied by the numerical constant σ f that depends on the spacetime dimension. Also note that we are subtracting at each time the expectation value in the adiabatic case, for which we are using δt = 10. In even spacetime dimension d, the plots corresponds to having the renormalization scale set to k 0 = 1. The procedure is exactly the same as in the scalar case. We define dimensionless variables q = kδt and κ = mδt and then use the hypergeometric series expansion to get the leading terms in a κ-expansion. To get the counterterms to that order we can also expand for large q. The difference in this case is that the expectation value given in eq. (3.9) also requires to compute the time-derivative of the hypergeometric function and in general, this can be an involved task. However, we should notice that the only time dependence in the hypergeometric function is in the last argument, i.e., z = (1 + tanh t/δt)/2. The rest of the coefficients do not depend on time. Then, where we use the usual trigonometric identities to express the time-derivative of z(t) as a function of z(t) itself. With this in mind, we can expand our hypergeometric series and note again that only the few first terms are needed in order to get the leading order κ behaviour. For odd d ≥ 3, we obtain ψ ψ ren = (−1) 14) which correctly gives the expected scaling behaviour mδt 2−d . Fig. 20 shows how the numerical solutions approximate this leading order analytic term as δt → 0. As the curves approach the leading order analytic solution (3.14) shown with the dashed red line, δt gets smaller, with solid (numerical) curves going from δt = 1 to δt = 1/10. For d = 3 we also included the curves with δt = 1/50 and 1/100.
For even d the situation is again a little bit different, since we have an extra logarithmic term. Then we can define where k 0 is the renormalization scale. Then, the universal term yields, for d ≥ 4, (3.16) In contrast, ψ 2 is much more complicated and we expect that it is not universal. For the present tanh quenches, ψ 2 can be written as (3.17) where again z(t) = 0.5 + 0.5 tanh(t/δt). As in the case of the scalar field, these results support the holographic scaling where power-law growth is enhanced by a logarithmic factor in even d. We can appreciate these additional logarithmic factors by looking at figs. 21(a) and 21(b). There we divide out by the expected power-law scaling and we still see that the expectation value is growing as we decrease δt. Finally, by using both eqs. 3.16 and 3.17, in figs. 21(c) and 21(d), we can see that the numerical solutions approach the analytical leading term for sufficiently small δt's. In panels (a) and (b), we only divide by the expected power-law scaling. As we reduce δt from δt = 1/10 to δt = 1/100 for d = 4 and from δt = 1 to δt = 1/10 for d = 6, we see that the expectation value still grows, indicating the presence of an extra logarithmic factor. If we take the latter into account and divide by it as well, we find panels (c) and (d), where we see that the curves now converge towards the analytic expression (red dashed line).

Quenches in general interacting theories
Both the results in [16] and in this paper show that different observables in free field theories after a smooth fast quench obey the same universal scaling relations as in quenches in holographic theories, as shown in [14,15]. As the holographic CFT's are implicitly strongly coupled, we seem to have found the same scaling at two ends of the spectrum of possible interacting quantum field theories. Hence we should expect that the same result holds for a large variety of quenches in a wide range of interacting theories. In this section, we give arguments that the universal scaling in eq. (1.2) appears quite generally for fast quenches. The crucial assumption will be that the interacting theory which is being quenched approaches a UV fixed point, i.e., its UV properties can be described by an appropriate CFT.
To motivate the general argument, we begin by considering quenches with a pulse profile in a CFT, as presented in [16]. In this case, we can use conformal perturbation theory. The starting point is a CFT which is deformed as follows where O ∆ is a relevant operator with dimension ∆ < d. We assume the profile for the corresponding coupling λ(t) has the form where δλ is the maximum coupling value and h(y) is some smooth function that goes from 0 to 1 and back to 0 (at least roughly) in the interval y = 0 to 1. Then, our coupling (4.2) has the form of a pulse in the time interval t ∈ [0, δt] with a maximum δλ. Note that this form essentially matches that of the profile (2.76) that we analyzed in section 2.3. There the mass profile was a pulse that goes from the critical point (i.e., the massless theory) to the same critical point after passing through some maximum mass at t = 0. Clearly this condition is not strictly necessary to obtain the universal scaling (1.2), as we have shown in quenches with a tanh profile for both the scalar and fermion masses yield the same scaling. However, the above framework will help to formulate our general argument. Basically, since our theory is critical at both infinite past and infinite future (and anywhere outside the interval t ∈ [0, δt]), we can calculate the expectation value of our operator using conformal perturbation theory, which yields where all expectation values are evaluated in the critical (conformal) field theory. Now, the first term in the RHS vanishes because O ∆ is a relevant operator, so its expectation value O ∆ (0) CFT must vanish. The second term is the linear response term where the retarded correlator is given by The next-to-leading term is given by three-point correlator, 4.5) As in the free field cases analyzed in this paper, the expression in eq. (4.3) is usually UV divergent so we need to regulate it by adding counterterms in order to get a finite expectation value. We will assume that renormalization can be done without problems to have a finite value that only depends on two renormalized parameters, δt and δλ. This assumption relies on the fact that, as mentioned in [25,26], we do not expect that the quench protocol would lead to any unconventional RG flows.
It is natural to expect that these counterterms are precisely given by the adiabatic expansion, as we have seen explicitly for the free field theory. Again, the reason is that the UV contributions to these quantities are insensitive to the quench rate so long as the rate is slow compared to the UV cutoff scale. Our protocol is chosen such that the rate is fast compared to the scale of the relevant coupling, as in eq. (1.1), but always slow compared to the UV cutoff, i.e., Λδt 1. For free field theories, it is easy to perform the adiabatic expansion since all we had to do is solve the wave equation in a WKB expansion, as described in section 2.1.1. For interacting theories, this is no longer the case and we have to use the standard procedure in quantum mechanics starting with an expansion of the wave functional in terms of instantaneous eigenstates of the (time-dependent) Hamiltonian. Now, as all correlators in eq. (4.3) are CFT correlators, they should be independent of the parameters δλ and δt. So, basically, δt will set the scale for the integrals and dimensional analysis will fix the form of all the possible terms in the expectation value. This means that O ∆ (0) ren = a 1 δλ δt d−2∆ + a 2 δλ 2 δt 2d−3∆ + · · · , (4.6) where the constants a n are finite numbers by assumption. Then, we can see that the first term, the linear response, is responsible of producing the universal scaling found, i.e., O ∆ ∼ δλ/δt 2∆−d . However, we still have an infinite set of nonlinear contributions and so the next step is to show that the these become negligible once we take the limit of fast quenches (1.1). For that, it will be easier to define a dimensionless effective coupling, g ≡ δλδt d−∆ , so that eq. (4.6) becomes simply That is, conformal perturbation theory (4.3) has expressed the expectation value in terms of a series expansion in terms of the dimensionless coupling. The quenches we are considering correspond to keeping δλ fixed, while taking δt → 0. This means that we are taking our effective coupling small, i.e., g → 0, since we are quenching a relevant operator with ∆ < d. Hence with these protocols, the expansion (4.7) is a very effective perturbation expansion and the leading behaviour is determined by just the first term, Of course, as we already noted, this term gives the desired scaling that was found in our previous calculations. We should note that any pulsed quench will then give the desired scaling, independent, for instance, of the underlying CFT. Hence the present argument encompasses both the holographic CFTs of [14,15] and the massless free fields studied here in previous sections. In our discussion of free field quenches, we found that the universal scaling behavior (1.2) is valid for profiles which are lot more general than the pulses considered above. Indeed, we now argue that the same scaling applies for general profiles subject to certain constraints and for general field theories subject to the assumption that the UV properties are described by a conformal fixed point. 16 That is, we regard our original theory as emerging from an RG flow away from some perturbed CFT in the UV with the action where λ 0 is the coupling constant for some relevant operator O ∆ (x). Now consider a quench where the profile of the coupling λ(t) only varies in the time interval t ∈ [0, δt]. At early times, λ(t) will simply be fixed at λ 0 while after the quench it will take another constant value λ 1 . For example, consider a coupling which interpolates between constant values λ 0 and λ 1 We leave the details of the function F (y) unspecified other than F (y ≤ 0) = 0 and F (y ≥ 1) = 1 and the maximum is finite with F max ≥ 1. Further this profile may dip below zero by some finite amount and so we specify the minimum as F min ≤ 0. Implicitly, we are also assuming that the profile is smooth. Now we will work in the regime where We will calculate the expectation value of the operator at some time t which is earlier than (or soon after) t = δt. Now motivated by the conformal perturbation expansion in eq. (4.3), we evaluate the change in O(t) relative to the initial theory (4.8) by 16 In many respects, the following argument closely resembles the holographic analysis in [15]. expanding in δλ, i.e., where G R,λ 0 denotes the retarded Green's function for the deformed CFT in eq. (4.8) and similarly K λ 0 denotes the analogous three-point correlator (4.5) in this deformed theory. Of course, the first term in this expansion corresponds to the linear response.
In writing the explicit range for the time integrals in eq. (4.10), we have used the fact that the function F (y) vanishes for y ≤ 0. Since all of the correlators in the above expansion are retarded, i.e., only have support within the past light-cone, the spatial integrals are also limited to a range of order t ≤ δt. That is, the integrals in eq. (4.10) only receive nonvanishing contributions from correlators where the operators are separated by a proper distance of less than O(δt). Now the fast quench regime defined by eq. (4.9) implies that these separations are all small compared to the inverse mass scales of the quenched theory. Hence the correlators will basically be the same as the CFT correlators, in eq. (4.3) and up to small corrections, the integrals again all scale with the power of δt determined by dimensional analysis. Therefore the change in the expectation value takes a general scaling form, 11) with none of the IR scales defining the deformed theory appearing in the problem. Again, the leading behaviour is determined by the linear response, i.e., the term linear in the dimensionless coupling, and hence the change in the expectation value has the desired scaling, δλ δt d−2∆ . Further, the diffeomorphism Ward identity (2.51) still applies in the present context. Hence energy is only injected into the system while ∂ t λ(t) is non-vanishing, i.e., only in the interval 0 < t < δt. But this is precisely the interval in which the change in the expectation value was evaluated in eq. (4.11) above. Now if we further assume that ∆ > d/2, then this change will be large compared to the unquenched expectation value in the fast regime. Hence integrating the right hand side of eq. (2.51) will lead to the expected scaling of the energy, as given in eq. (1.3).
Hence we have argued that the universal scaling in eqs. (1.2) and (1.3) will emerge for a broad variety of quenches in a wide class of interacting field theories. Of course, the above framework could be made even more elaborate, e.g., by introducing further deformations in the initial theory (4.8). Again, the first essential ingredient in our argument was that the interacting field theory under study can be considered to emerge in the infrared from an RG flow away from a conformal fixed point in the UV. Further, we are considering fast quenches where the quench rate 1/δt is much larger than any of the IR mass scales defining the initial theory or appearing in the quench protocol. The upshot of this is that when δt is the smallest physical length scale in the problem, the early time response is entirely governed by the conformal field theory at the UV fixed point, which explains its universality. In particular, the scaling behavior is independent of the details of the protocol so long as eq. (4.9) is obeyed.
We can extend this discussion to make explicit the independence of the scaling behavior from the initial and final mass scales appearing in the quenches of the free field theory in section 2.4. As we have seen above, the leading result is given by the linear response. Hence we consider the linear response answer for φ 2 in free scalar field with a mass profile similar to one considered there, i.e., a profile which interpolates between m 2 i and m 2 where the function F (y) rises from zero around y = 0 and quickly settles to 1 soon after y = 1. The intitial and final masses, as well as δm are small compared to the quench rate The change in the expectation value is given by a generalization of eq. (2.97) where we have used the fact that the function F (y) vanishes for y < 0.
To estimate this, consider for example a function F (x) which is piecewise constant Then for any t ≤ δt, eq. (4.14) becomes where q = kδt. Clearly m f has dropped out of this expression. Furthermore to the leading order in the limit (4.13), the integral in (4.15) becomes independent of m i as well. This leading answer is the same as in the case of a quench from a CFT. Note that if we use the expression (4.10) for times much longer than δt, we will need to address issues of infrared divergences associated with conformal perturbation theory for constant deformations [28]. For the question we are addressing here, we do not need to do this. For a recent discussion of our scaling result in a theory with an infrared regulator, see [29].

Conclusions
In this paper, we have expanded on the results of fast but smooth quantum quenches that we previously presented in [16], and extended the results to more general quenches. We have given details of our calculations in free field theories, where both numerically and analytically we obtain the same scaling relations as in previous holographic studies of the same kind of quenches [14,15]. This universal behavior in the early time response was found in a variety of quench protocols which interpolate between arbitrary constant masses so long as the quench rate 1/δt is large compared to all other physical mass scales in the problem. In section 4, we provided a general argument that the universal scaling in eqs. (1.2) and (1.3) will appear in fast quenches of any quantum field theory which flows from a conformal fixed point in the UV, i.e., for any theory that can be described as a CFT deformed by some relevant operator(s). The scaling is purely a property of the UV conformal field theory, which emerges at early times as long as the duration of the quench is short compared to all other physical length scales in the problem, as in eq. (1.4).
A key ingredient in our work, and in the corresponding holographic studies [14,15], is the renormalization of the underlying quantum field theory. Bare quantities, such as the expectation value O ∆ , are UV divergent and counterterms are needed in order to extract physically meaningful quantities. The problem here for quenches is quite similar in spirit to quantum field theories in curved space-times, e.g., [18][19][20]. In that case, the required counterterms involve operators made out of quantum fields, as well as curvature tensors of the background space-time. As discussed in section 2.1.3, for a global quench with a time-dependent mass, we need to add counterterms involving time derivatives of the mass function. In fact, to properly renormalize the expectation value of the stress tensor, we should also consider the theory in a curved background and include additional counterterms involving curvatures -even if we are only considering these expectation values in a flat space background.
However, we are still left with the problem of determining the precise coefficients of the counterterms which render the renormalized observables finite. We argued that these coefficients can be determined by examining the quenches in an adiabatic limit and in section 2.1.1, we demonstrated explicitly how to construct the necessary counterterms order by order in the expansion for slow quenches. Moreover, this procedure does not depend on any specific mass profile and so, the resulting counterterms should be universal. We verified that claim by correctly regulating quenches with a variety of different mass profiles using the same counterterms. Of course, it may appear surprising that an adiabatic expansion, which is an expansion in time derivatives, yields the correct counterterms for a fast quench. We argued that the physical reason behind this is that high momentum modes won't see whether the quench is fast or slow, so long as the quench rate is smaller than the cutoff scale Λ. In our cases we managed to take that cutoff to infinity while renormalizing the physical quantities, so we could expect that the counterterms would be the same in both slow and fast quenches. It would be interesting to test these assumptions in interacting field theories. Quenches in the large-N vector model, for instance, have been studied previously in the literature [27]. This would be a good place to make explicit calculations and verify whether our intuition holds even when we have interacting theories.

Renormalized observables
As emphasized above, our considerations refer to the renormalized quantities which require 'removing' various UV divergences in our calculations. While this is, of course, the standard approach in quantum field theory, one may still ask how our renormalized observables would be related to measurements made in a physical experiment, where implicitly there is a finite UV cutoff? As a simple analogy, let us consider a quench which consists of suddenly applying external pressures to a crystal. The phonons in the crystal would provide the analog of our quantum fields, i.e., at least in a certain regime, they would have a QFT description. The quench will 'excite' the final state of crystal in two ways. Naturally, the quench will generate phonon excitations in the crystal but the external pressure may also deform the crystal structure in the final configuration e.g., modifying the dispersion relation for the phonons. The work done in deforming the crystal structure would then be the analog of the changes in the divergent 'zeropoint' energy that appears in the bare expectation value E and which is subtracted by introducing mass-dependent counterterms to produce the renormalized energy density. Similarly, the energy available in the phonon excitations would correspond to the final E ren . The latter is the energy that can be accessed and manipulated by probing the system with local operators. Let us add for the analogy of the crystal quench becomes more precise if we also insist that the quench time δt is larger than the lattice spacing, which provides the UV cutoff scale.
While the above analogy should make clearer the role of bare and renormalized quantities in a physical system with a finite cutoff, one may still ask what quantities would appear in experimental measurements. Answering this question becomes even more complicated for even dimensions, where in section 2.1.1 we found that logarithmic divergences introduced various renormalization ambiguities. Such ambiguities were also discussed in the holographic context in [14,15]. The resolution there is that various fiducial experimental measurements would be made to fix these ambiguities. For example, examining eq. (2.32) for d = 4, we find that there will be two such logarithmic terms. 17 However, with some thought, we can see that the associated renormalization scales can be fixed by first measuring φ 2 ren and E ren at some fixed finite mass. Implicitly in the previous discussion but more generally, we can work with quantities which are free of UV divergences by comparing expectation values at different times or in different quenches. For example, as discussed in section 2.7, the difference φ 2 quench − φ 2 f ixed appearing in eq. (2.119) is completely finite for quenches in d = 3. Similarly, one can produce UV finite quantities by comparing the results for different quench protocols or by quenching different initial states with the same quench protocol, as in briefly discussed in appendix B. Of course, another family of UV finite observables would be correlators measured with finite separations, e.g., as in eq. (2.108). Further, one may be able to find evidence of universal scaling in the early time response with a strategic choice of the positions in the correlator.
Of course, it would also be interesting to analyze cases where the cutoff remains finite, e.g., in some lattice model. Though the analysis would be more complex in such a case, we expect that our universal scaling properties should emerge in the regime where the scales are properly distinguished, i.e., in a regime where Λ 1/δt m. In fact, one might expect that as 1/δt approaches the cutoff scale, one would recover the results of an instantaneous quench.

Comparison to instantaneous quenches
Finally, we should comment on the relation between our smooth quenches and the instantaneous (or abrupt) quenches that are usually studied in the literature [5][6][7]. Some preliminary discussion of the comparison between these two classes of protocols was given in section 2.6 -see also section 2.7 -and a more detailed discussion will appear in [17]. Here, the universal scaling in eqs. (1.2) and (1.3) suggests that divergences will appear as δt → 0 (whenever ∆ ≥ d/2). This would seem to contradict instantaneous quench results. However, as we already discussed in [16], these two types of quenches are different: while the present quenches evolve smoothly in a time-dependent scheme, the instantaneous quench approach can be thought as the evolution of a farfrom-equilibrium initial state evolving under a fixed, time-independent, Hamiltonian. The scalings discussed in this paper hold for renormalized quantities and as emphasized above, the renormalization procedure demands that the quench rate is slow compared to the UV cutoff scale. On the other hand, instantaneous quenches necessarily involve quench rates which are fast compared to all scales, including the UV cutoff. Indeed we have explicitly shown that for free field theories, the momentum space correlator agrees with that for an instantaneous quench only when the momenta are small compared to the quench rate 1/δt. Local quantities, like the one-point function of the mass operator or the energy density, involve an integral over all momenta and so this condition does not hold.
However, this constraint above may still hold effectively if the contributions to the integral at high momentum are suppressed for other reasons. One case where the latter might apply is at late times after the quench. The intuition behind is that at late times, we expect only low energies (or momenta) contribute and hence the observables for fast smooth quenches and for instantaneous quenches may agree at late times. Section 2.7 presents some preliminary evidence for this conclusion. In the free bosonic theory for d = 3, we found that at sufficiently late times, the one point function becomes independent of δt and grows logarithmically in time. Further, this late time growth exactly agrees with the result from an instantaneous quench. For d = 5, we showed that the late time result for a smooth fast again becomes independent of δt as δt → 0. However, as we will discuss in [17], this answer only roughly agrees with the expectation value at late times after an instantaneous quench. More generally, the expectation values generated by the two different protocols fails to agree even at late times in higher dimensions [17] and hence the precise agreement in d = 3 is quite exceptional. However, the disagreement found more generally should come as no surprise since the expectation value φ 2 involves a momentum integral up to the cutoff scale where, as we already argued, agreement should not be expected. However, we should add that a detailed analysis reveals agreement for the late time correlators at finite spatial separations which are large compared to δt. Again a full discussion of these issue will be presented in [17].

Higher spin currents
One interesting feature of the free field theories studied here is that they contain an infinite family of conserved higher spin currents. In section 2.2.6, we began a study of the response of the higher spin currents in fast smooth quenches. In particular, we discussed the construction of the higher spin currents in the case of massive free fields -see also appendix A. This construction naturally leads directly to a generalization of the diffeomorphism Ward identity (2.51) for the higher spin currents. In general, there is a hierarchy of generalized Ward identities (2.74), which can be used to understand how the 'work' done in varying the mass parameter changes the various higher spin 'charge densities.' In the fast quench regime (1.1), the latter yields a simple universal scaling property (2.75) for these higher spin densities, i.e., j (s) t···t ∼ (m 2 ) s 2 +1 /δt d−4 . Of course, for spin-2 and spin-0, these are scalings of the energy density and the mass operator, in accord with eqs. (1.2) and (1.3). We explicitly carried out this construction and demonstrated the corresponding scaling for the spin-4 current. Hence it will be interesting to explicitly construct all of the massive higher spin currents and to explicitly find the corresponding generalized Ward identities. Of course, it would also be interesting to develop a better intuition for the physical meaning of this hierarchy of Ward identities and the resulting universal scaling.
where the traces coefficients have been chosen to make the current traceless and the whole expression should be symmetrized in all four indexes. When we take the divergence, it is straightforward to verify that this current is conserved: as all terms have ∂ a a φ or its conjugate, which vanishes because of the equations of motion in the massless case. To generalize eq. (A.1) to the massive case, we will need to add terms proportional to the mass squared, so that all these terms now are cancelled but upon evaluation in the massive equation of motion ∂ a a φ − m 2 φ = 0. So, for instance, the first term in the RHS of eq. (A.2) should be cancelled with one of the form + m 2 48 ∂ b φ∂ cd φ * . Now, all the possible m 2 terms we can add to the current are of the form where A, B, C are constants to be determined (as always, the most general form of the current should be symmetrized). By taking the divergence we obtain, ∂ a j (4)m 2 abcd = 3(A − 2B) (∂ b φ∂ cd φ * + ∂ b φ * ∂ cd φ) + 3A (φ * ∂ bcd φ + φ∂ bcd φ * ) + 3(A − B + 2C) (η bc ∂ a φ∂ ad φ * + η bc ∂ a φ * ∂ ad φ) + (A.3) +3A (η bc φ * ∂ a ad φ + η bc φ∂ a ad φ * ) − 3B (η bc ∂ a a φ * ∂ d φ + η bc ∂ a a φ∂ d φ * ) and so, we are left with a 3 × 3 system to solve for A, B, C in order to have j (4) abcd + j and then we have the full generalized 4-spin conserved current for the case of a massive scalar field.
We can do an analogue procedure in any spacetime dimension and we get j (4) abcd = j (4)m=0 abcd + m 2 η ab 4(d/2 + 2)!(d/2)! d 2 j (2) 3 + m 2 j 0 , (A.5) where j 3 = η cd ∂ e φ * ∂ e φ, (A.8) The crucial fact for our discussion in section 2.2.6 is that we can write the current as the sum of the minimally coupled and the conformally coupled spin-2 current and then it is direct to evaluate the generalized Ward identity.
Finally, we should say that this procedure is, in principle, easily generalized to any higher spin current. However, the procedure becomes tedious as the number of terms in the massless current grows quickly with the spin and so does the number of possible terms that should be canceled with mass terms.

B Scaling of excited states in the scalar quench
It is interesting to also analyse the behaviour of excited states under a quench. This gives an extra observable to evaluate and it may be particularly useful in case one wants to make explicit contact with experiment. If we take the case of even dimensions, for instance, extra regulator ambiguities appear in the problem, as discussed in section 2.1.1. So, if someone is performing an experiment, before looking at the scalings and so on, one should establish a way to fix these ambiguities. Interestingly, after being fixed, one should be able to compare different states using the same protocol, so excited states become useful observables to evaluate the behaviour of the system. We can think of different possible excited states such as giving the system some excitations of in-modes or even think about a coherent state of in-modes. In any case, one interesting example is to compute φ 2 n ≡ in, 0|a n k φ 2 a †n k |in, 0 , for any momentum k and any number of excitations n.
However, we already know the exact solution to φ under the quenches we are considering, i.e., see eqs. (2.7) and (2.8). So, in order to evaluate the excited expectation values we only need to use repeatedly the property of commutation of the a k modes -see eq. (2.7). Explicitly, what we need to find are expressions of the type in, 0|a k · · · a k a k a † k a † k · · · a † k |in, 0 and in, 0|a k · · · a k a † k a k a † k · · · a † k |in, 0 . After some algebra we get where φ 2 0 means the vacuum expectation value and the hypergeometric function is evaluated at the same arguments as in the main body of this article but at a fixed momentum k. Now one can ask whether the difference between the excited states expectation value and the vacuum also scales as δt → 0. However, it is quite direct to show that as δt goes to zero with fixed momentum k, the hypergeometric function goes to 1, and so the difference φ 2 n − φ 2 0 would go to some constant depending on k and n but would not scale with some power of δt, as found for the vacuum expectation value.