Quantum Quench in Non-relativistic Fermionic Field Theory: Harmonic traps and 2d String Theory

We investigate a class of exactly solvable quantum quench protocols with a finite quench rate in systems of one dimensional non-relativistic fermions in external harmonic oscillator or inverted harmonic oscillator potentials, with time dependent masses and frequencies. These hamiltonians arise, respectively, in harmonic traps, and the $c=1$ Matrix Model description of two dimensional string theory with time dependent string coupling. We show how the dynamics is determined by a single function of time which satisfies a generalized Ermakov-Pinney equation. The quench protocols we consider asymptote to constant masses and frequencies at early times, and cross or approach a gapless potential. In a right side up harmonic oscillator potential we determine the scaling behavior of the one point function and the entanglement entropy of a subregion by obtaining analytic approximations to the exact answers. The results are consistent with Kibble-Zurek scaling for slow quenches and with perturbation calculations for fast quenches. For cis-critical quench protocols the entanglement entropy oscillates at late times around its initial value. For end-critical protocols the entanglement entropy monotonically goes to zero inversely with time, reflecting the spread of fermions over the entire line. For the inverted harmonic oscillator potential, the dual collective field description is a scalar field in a time dependent metric and dilaton background.


Introduction
A common way to study non-equilibrium properties of quantum field theories is to subject them to a quantum quench, i.e. introduce an explicit time dependence to parameters which appear in the lagrangian. Among other things, this is interesting for several reasons. One motivation is to study equilibration and possible thermalization of these systems. Suppose the time dependent parameters approach constant values in the far past and future, and the system is initially in the ground state. The quench then excites the system. At late times, when the parameter again becomes a constant (which is generally different from the initial value), one would like to know the nature of the excited state, and if it is approximately described by a thermal state in some appropriate sense. A second motivation -which is one of our main interests -is to study dynamics in critical phase transitions [1]. Suppose the initial hamiltonian is gapped, while the quench protocol crosses or approaches a critical point where the gap vanishes. On general grounds one expects that various observables would obey universal behavior.
An early example of such a universal behavior is Kibble Zurek scaling for global quenches [2,3] (where the parameters depend only on time). This holds in many systems when the time scale of the quench δt is large compared to the inverse of the initial energy gap E g . In this case, the initial time evolution is adiabatic. However since the instantaneous gap is descreasing with time, adiabaticity breaks down at some time called the Kibble Zurek time t KZ . This is typically determined by the Landau criterion, where E gap (t) is the instantaneous energy gap. This equation then determines t KZ in terms of δt. According to the assumptions of Kibble and Zurek, the system soon enters a diabatic regime, and the instantaneous correlation length at the Kibble Zurek time is the only length scale in the problem. In the following we will follow standard nomenclature to distinguish several classes of quench protocols. The first protocol is called a trans-critical protocol (TCP). Here the system begins in a gapped phase and the coupling varies monotonically across a critical value, and approaches a final value which also corresponds to a gapped phase. The second is called a cis-critical protocol (CCP) where the time dependence is not monotonic. Here the system starts from a gapped phase, approaches a critical point and reverts back to a constant value which also corresponds to a gapped phase. The third protocol is called a end-critical protocol (ECP). Here the system begins in a gapped phase and monotonically approaches a critical point at infinitely late time. In TCP or CCP, the response at early times then scales as appropriate power of the correlation length, leading to a scaling as some universal power of the quench time scale. For example the one point function of an operator O will scale as O ∼ ξ −∆ KZ (1.2) where ∆ is the dimension of O. For ECP, the appropriate scaling variable is the energy scale at the Kibble Zurek time. At any given time the response will be adiabatic for sufficiently large δt, while for a small enough δt there will be a Kibble Zurek regime [4].
Even though these assumptions appear to be drastic, such a scaling -together with an accompanying mechanism for defect formation in symmetry breaking transitions -appears to hold for many systems. Kibble Zurek scaling has been studied in a variety of solvable models and in holographic setups [5]. The latter provide some insight into the origins of universality. The best known results involve one point functions (e.g. defect density) and correlation functions. However similar scaling holds for the entanglement entropy of a subregion in some model 1 + 1 dimensional systems.
At the other extreme is instantaneous quench where a sudden change of a parameter causes the system to go from a gapped phase to a critical point abruptly. In this case, universal results are known for correlators and entanglement entropies of 1 + 1 dimensional systems [6,7]. Of particular interest is the spread of entanglement with time [6] -this kind of spread has been conjectured to hold for higher dimensional systems [8] and there has been evidence for this in holographic calculations [9] as well as in free field theories.
More recently it has been found that in a relativistic theory there is an intermediate regime where a different universal scaling holds [10]-a result which was first found in holographic calculations [11,12] and later found to hold quite generally. Consider a relativistic quantum field theory in d dimensional space-time which is obtained by the RG flow from a UV fixed point. The action can be then written as Here S CF T stands for the conformal field theory action at the UV fixed point and ∆ denotes the conformal dimension of the operator O ∆ ( x, t) in this CFT. The time dependent coupling λ(t) goes from a constant value λ 0 in the infinite past some other value λ 1 in the distant future, and the time dependence is in some time interval of size δt. Then this regime is defined by where λ ± denote the largest and smallest value of the coupling and δλ is the excursion of the coupling during the quench process. In this regime the one point function soon after the quench is over scales as This is a result in any relativistic field theory, and follows from two basic properties [10]. The first is causality. The second property is that the causal Green's functions of a massive theory become those of the UV conformal theory for space-time separations which are small compared to the inverse mass gap. Once these properties hold, it turns out that the dimensionless parameter which controls time dependent perturbation theory is the combination of the coupling with an appropriate power of δt, and all other scales go away. This combination is small in the fast quench limit and the result (1.5) follows from the lowest order perturbation theory. This regime of scaling has been investigated explicitly in free field theories with time dependent masses and in conformal field theories with relevant and marginal deformations [13]. In continuum free theories there appears to be a smooth transition between Kibble-Zurek and Fast scaling regimes [14], while in lattice theories this connects to the abrupt quench regime at quench rates at the scale of the lattice spacing [15]. Apart from one point functions, the whole range of scaling behavior is visible in quantities like the entanglement entropy [16] as well as circuit complexity [17] 1 . In many situations, particularly in experimental setups, one is interested in non-relativistic systems. Our ultimate goal is to investigate whether there are universal scaling laws which hold in non-relativistic systems. While Kibble Zurek scaling is expected to hold, the status of fast quench scaling is unclear. In specific models where non-relativistic Lifshitz type dispersion relations appear, e.g. the anisotropic critical points of the Kitaev model one indeed finds fast quench scaling with appropriate scaling dimensions [15]. More generally, Lieb Robinson bounds [19] for lattice non-relativistic systems may provide the necessary ingredient. Indeed in recent work in lattice models with dynamical exponent z = 1 it has been found that the spread of entanglement following a sudden quench indeed has an effective finite velocity [20]. However, such a finite speed has been also observed in non-relativistic systems which do not obey Lieb-Robinson bound [21].
In this work we study the issue of scaling in a specific solvable system : a system of N mutually non-interacting non-relativistic fermions in a harmonic or inverted harmonic potential with a time dependent frequency and a time dependent mass. Using the results of [22] will show how the problem of quantum quench with some smooth quench profile in such systems can be solved analytically once one can solve a nonlinear equation (Ermakov-Pinney (EP) equation). The solutions of this equation can be in turn determined in terms of the solutions of the classical equation of motion of a single particle in the same harmonic potential.
Indeed harmonic traps are of considerable interest in experimental cold atom physics : quantum quench experiments often involve release of particles from harmonic traps.
Our interest in the inverted harmonic oscillator potential on the other hand stems from its connection to two dimensional string theory [23]. As is well known, the double scaled limit of the singlet sector of the quantum mechanics of a single hermitian matrix reduces to a set of fermions in an inverted harmonic oscillator potential. The string coupling appears as the mass of the fermion. Thus two dimensional string theory with a time dependent coupling reduces to the problem of fermions with time dependent mass in such a potential 2 . String theory with time dependent string couplings have been studied extensively in the context of AdS/CFT to investigate thermalization via black hole formation [24]. In a different context these have been used as models of AdS cosmology [25][26][27][28][29], but the outcome has been rather inconclusive. Here we hope to obtain exact results in a simplified situation.
In this paper we will set up the formalism necessary to solve both the harmonic and inverted harmonic potential problems. We present detailed results for the problem in harmonic trap : the problem of two dimensional string theory will appear in a future publication [30].
We will solve the quantum mechanical time evolution of such a system for interesting time dependent frequencies of the CCP and ECP type and calculate the early time response of one point functions as well as entanglement entropies for a sub-region for arbitary quench rates to find the scaling behavior in various regimes. We will also explore the late time behavior of the entanglement entropy. We find Kibble-Zurek scaling for slow quenches, while for fast quenches we show that the result scales in a way which is consistent with time dependent perturbation theory. At late times the entanglement entropy in a CCP oscillates with an amplitude which appears to remain constant in time. This reflects the lack of thermalization of the system. For the ECP the entanglement entropy monotonically goes to zero as a power law in time, reflecting the fact that the particles can now spread all over space.
Such solvable systems have played a major role in providing insight into scaling properties of quantum quench in continuum relativistic theories and in spin systems which can be reduced to lattice versions of relativistic fermions [10,14,15]. As we will see, our example may not be the appropriate setup to explore a possible universal scaling at fast rates. Neverthless, we hope that these exact solutions will provide some insight into the general problem.
Abrupt quantum quench in a system of free non-relativistic fermions which arise from Matrix Quantum Mechanics with various potentials has been investigated in several papers [33][34][35][36]38]. In particular [33] has extensively studied the problem in terms of the dynamics of the Wigner phase space density, investigated approach to a generalized Gibbs ensemble and discovered interesting dynamical phase transitions. The papers [34][35][36] deal with the fermion problem directly in the presence of various kinds of abrupt quenches. Other aspects of the dynamics in such fermion systems (e.g. shock wave formation) have been studied in [37].
The paper [38] considers the dynamics of the Wigner phase space density as well as a system of bosons and fermions using methods similar to us, in particular the EP equation. The paper [39] considers slow smooth quenches for bosons also using the EP equation. The EP equation has also been used to study entanglement dynamics following an abrupt quench in a harmonic chain in [40].
Our work is complementary to these papers. We are interested in studying scaling of various quantities as functions of the quench rate. We have been able to find exact analytic solutions to several smooth quench protocols which we use for this purpose.
correspondence [31,32]. Introducing a time dependent mass for such fermions naively corresponds to a time dependent coupling of the Yang-Mills theory. However this breaks supersymmetries : the truncation of matrix models and therefore fermions do not hold any more.
In section 2 we set up the second quantized fermion field theory and show how this can be solved exactly for ±x 2 potentials in terms of a function ρ(t) which satisfies generalized Ermakov-Pinney equation and show how to obtain its solutions. In section 3 we quantize these theories in the Heinsenberg picture "in" state and show how observables can be expressed entirely in terms of ρ(t). In section 4 we provide exact solutions for some CCP and ECP quench protocols for the harmonic problem. Sections 5 -7 contain our results for the one point function of the quenched operator and the entanglement entropy for these protocols and their scaling as functions of the quench rate. Section 8 deals with comments about the behavior of the phase space density.
2 Fermion field theory Consider a system of N non-relativistic fermions in 1+ 1 dimensions with a hamiltonian given by where m(t), ν(t) are real smooth functions. The Schrodinger picture fermion field operators above satisfy the usual anti-commutation relations The condition that the total number of fermions is N then leads to the constraint The plus sign in (2.1) is the hamiltonian of particles with a time dependent mass in a harmonic trap with a time dependent frequency. The minus sign with ν = 1 is the hamitonian of the singlet sector of the double scaled single hermitian matrix quantum mechanics which is dual to two dimensional string theory with a time dependent string coupling g s (t) = m(t). The Heisenberg picture equation of motion is the Schrodinger equation In the following we will set = 1.

The general solution
In terms of a new time variable τ dτ = dt m(t) (2.5) we can transfer the time dependence of the mass to the frequency term and (2.4) becomes A solution to the equation (2.6) can be obtained in terms of the solution of the Schrodinger equation with a constant mass and a constant frequency as follows [22]. First define a new field Φ(x, τ ) by Secondly, make a change of variables and the function ρ(τ ) satisfies a generalization of the Ermakov-Pinney equation [41] Here the positive sign refers to the right-side up harmonic oscillator while the negative sign refers to the inverted harmonic oscillator of relevance to the hermitian matrix model. The latter case will be discussed in detail in [30].
In the adiabatic approximation the function ρ(τ ) is simply 1 √ ω(τ ) . A departure from this value denotes a departure from adiabaticity and describes the exact response. Furthermore the most general solution of (2.12) is given by where A, B, C are constants and f (τ ), g(τ ) are two linearly independent solutions of the classical equation of motion of a single particle moving in a harmonic (inverted harmonic) potential with the same time dependent frequency ω(τ ) where W r(f, g) = f ∂ τ g − g∂ τ f is the wronskian of the two solutions. By the equations of motion this is a constant in time and can be therefore evaluated at any time.
The problem of fermions with a time dependent mass in a harmonic (or inverted harmonic) potential with a time dependent frequency can be therefore reduced to a problem with a constant mass and a constant frequency. The only equation one needs to solve is the classical equation (2.14). As we will see below, many quantities of physical interest can be expressed entirely in terms of the function ρ(t) = ρ(τ ).
A general solution of the equation (2.10) has the form Φ n (y, T ) = N n e −iλ(n)T φ n (y) (2.16) where φ n denote a complete orthonormal set of eigenfunctions of the hamiltonian given by the right hand side of (2.10) with eigenvalue f (n). Then the above discussion implies that a general solution of the equation (2.6) may be written as The orthonormality conditions for the eigenfunctions φ n (y) then imply the orthonormality conditions for the solution (2.17).
The form of the solution (2.17) reveals another physical meaning for the function ρ(τ ). In the wavefunctions of a harmonic oscillator at fixed frequency ω, one can rescale out the frequency by x → √ ωx and τ → ωτ . The normalization of the wavefunction also involves ω 1/4 , as would be required by the rescaling of x. In our problem ρ(τ ) almost plays the role of such a time dependent rescaling. The term which spoils this is the phase factor which involves ∂ τ ρ(τ ). This is of course consistent with the fact that in lowest order of adiabatic approximation ρ .
The function ρ(τ ) is given by (2.13). The independent solutions f (τ ), g(τ ) and the constants A, B, C have to be chosen so that the solutions (2.17) satisfy the correct initial condition.

Solution in terms of Phase Space Density
It will be useful to think in terms of the Wigner phase space density operator The condition (2.2) then becomes where ⋆ denotes the Moyal star product. As shown in [42] the fermion field theory can be expressed as a path integral in terms of these variables with a co-adjoint orbit action.
Formulating the theory in terms of u(p, q, t) is particularly useful in the classical limit → 0, N → ∞ with N held fixed. In this limit the Moyal product reduces to an ordinary product.
In this limit the operator u(p, q, t) satisfies the equation If we make the change of variables provided (2.13) holds. This transformation is in fact a canonical transformation. Therefore the condition (2.20) that u(p, q, t) describes N fermions transforms into the condition dP dQ u(P, Q, T ) = N (2.25) The equation (2.24) is the equation satisfied by the phase space density operator for a system of fermions which is in an external harmonic (or inverted harmonic) potential with unit mass and unit frequency. Therefore once we know the solution for this latter case, we can find a solution of the time dependent case in terms of a solution of the equation (2.12).

Quantization and the "in" state
The quantization of the fermionic field theory proceeds in a standard fashion. Given a complete set of modes {ψ n (x, τ )} which solve the equations of motion the Heisenberg picture field operators may be expressed as where the oscillators satisfy the standard anti-commutation relations Different choices of modes determine different inequivalent quantizations related by Bogoliubov transformations. We will be interested in profiles of m(t), ν(t) such that they approach constant values m in and ν in as t → −∞, and their time derivatives approach zero. Furthermore we will have choices of m(t) such that when t → −∞, one also has τ → −∞. In particular our choices of The corresponding solution ψ n in (3.1) must then have the property that this is positive frequency in the far past, We will consider the Heisenberg picture state which is the "in" ground state,

Observables
The observables we will be interested in are the expectation value of the quenched operator and the entanglement entropy. We will now show that both these quantities can be expressed in terms of the corresponding quantities in the time independent problem and the function ρ(τ ). In the following we will consider the expectation value of the operator This is the operator which comes multiplied by the time dependent coupling ω 2 (τ ) once the theory is expressed in the time variable τ . In the spirit of response theory, the expectation value then measures the response of the system to the external driving. O(τ ) of our problem can be expressed simply in terms of the expectation value of the quenched operator in an auxiliary problem of a harmonic oscillator with unit mass and frequency, using (2.17) Using (2.16) and (2.17) this becomes, after a change of variables, where we have used the fact that the integral on the right hand side is the expectation value of the potential energy of a single harmonic oscillator with unit frequency in the state with quantum number n, and used the standard result. For fermionic systems, the entanglement entropy of a subregion A has an expansion in terms of cumulants of the particle number distribution [43][44][45]. In the leading order of large N the dominant term is the variance of the expectation value of the particle number in A, where the operator N A is given by where the integral is over the region A. This simplifies for the "in" state. Using the mode expansion (3.1) and the state defined in (3.5) it may be easily shown that This quantity can be also expressed entirely in terms of the expectation value of the phase density operator as follows (3.13) Expressing the above expectation values in terms of the mode functions one has Using (2.17) it then follows that the entanglement entropy can be expressed in terms of the entanglement entropy of a rescaled region in the ground state of the theory with a constant mass and frequency. If the subregion A is defined by a ≤ x ≤ b then the rescaled subegion is defined by

Results for fermions in Harmonic Oscillator Potential
For the right side up harmonic oscillator, the two independent solutions of the equation (2.14) may be therefore chosen to be such that To ensure that ρ(τ ) is real we then need to choose Therefore for this solution we have This yields the final form of the solution where H n (x) denotes the n-th order Hermite polynomial. This solution approaches the normalized solutions of the Schrodinger equation with a frequency ω in as τ → −∞. The oscillators in (3.1) with these modes are in the "in" oscillators.
We now provide exactly solvable quench protocols for fermions with a fixed mass m in a harmonic oscillator potential with time dependent frequencies. The two times t and τ are then related by τ = t/m.

Cis-Critical Protocol
The first protocol is a cis-critical-protocol (CCP). As described in the introduction in such a protocol the system starts from a gapped phase, approaches a critical point and then turns back to another constant value. In this work we choose a protocol where the initial and the final values are the same. More specifically we choose This corresponds to a trap which is smoothly removed for a finite interval of time and then re-introduced. The solution to the equation (2.14) which behaves as e −iω 0 τ is then given by where we defined The key function ρ(τ ) is then given by (4.3) with f (τ ) given by (4.6).

End Critical Protocol (ECP)
Another solvable quench protocol is the end critical protocol where the initial theory is a harmonic oscillator with a frequency ω 0 which monotonically descreases smoothly to a vanishing frequency at infinitely late times. This corresponds to a smooth release from a harmonic trap. Consider the slightly more general protocol with the real constants a, b chosen such that a > b to ensure reality of ω(τ ). Then the "in" solution of the equation (2.14) is given by where we defined The end critical protocol we consider has a = −b = 1 2 . The function ρ(τ ) which determines the time dependence of the observables considered above is shown in Figure 1 for both these types of protocol.
At early times ρ(−∞) = 1 √ 2ω in . For ECP ρ(τ ) monotonically increases and behaves as ρ(τ ) ∼ τ at large τ . For CCP ρ(τ ) initially increases and then starts oscillating. At late times these oscillations are around a mean value which is roughly the initial value 1 √ 2ω with an amplitude which remains constant in time and with a frequency approximately given by ω 0 .

The response and scaling : CCP
In this section we present the results of the expectation value of the quenched operator O = dxx 2 ψ † ψ at early times for CCP (equation (4.5)) and investigate their scaling behavior in various regimes. The details of the analytic approximations which lead to these results are given in Appendix A.

Slow Quench Regime
In the slow quench regime ω 0 δt ≫ 1 we can use the asymptotic form of gamma functions to obtain ρ(τ = 0). The leading expression for the one point function O at τ = 0 is, using (3.8), This result is consistent with Kibble-Zurek scaling. The Landau criterion with the instantaneous frequency given by (4.5) leads to which defines the Kibble-Zurek time τ KZ . We expect a scaling behavior only when τ KZ ≪ δt.
In this regime (5.3) leads to The condition τ KZ ≪ δt then becomes consistent with the slow quench condition ω 0 δt ≫ 1. This leads to the instantaneous frequency at the Kibble-Zurek time, According to the Kibble-Zurek argument ρ(τ ) in the middle of the quench (which is τ = 0) is roughly equal to its value at τ = τ KZ . Since the system is approximately adiabatic at τ = τ KZ this is in turn roughly equal to ρ adia (τ KZ ), the value of ρ for the fermions in a harmonic oscillator potential with a constant frequency ω KZ . From (2.12) this is simply which is in agreement with the result from the exact solution (5.2) upto a numerical factor.

Fast Quench Regime
We now consider the regime ω 0 δt ≪ 1. While we have the exact answer anyway, we are able to approximate the answer by suitable expansions and obtain analytic expressions when we have in addition ω 0 τ ≪ 1 . The latter are useful to make a comparison with perturbation calculations. First consider the response at a time τ which is in the range This is the response at early times. At late times, (see equations (A.19) to (A.28)) These results should also follow from usual time dependent perturbation theory. Let us discuss this for a general perturbation δω(τ ) 2 from the initial value. The leading term in the perturbation expansion is 5.12) where ω 0 denotes the expectation value in the ground state of the theory at τ → −∞ which is the harmonic oscillator with a constant frequency ω 0 and The Green's function which appears in the linear response can be calculated. The result is Now consider evaluating the response at a time which is of the order of τ ∼ +δt. In the fast quench regime ω 0 δt ≪ 1. The limits of the integral in (5.14) can be replaced by (−δt, δt). Suppose the form of δω(τ ) 2 is δω(τ ) 2 = δω 2 0 f (τ /δt) where δt is the time scale of the quench and f (x) some smooth function. Using the above form of the Green's function the linear response becomes in the fast quench regime In the protocol we are using δω 2 0 = ω 2 0 and f (τ /δt) = sech 2 (τ /dt). We therefore reproduce the scaling in (5.9).

The exact response
The exact response for CCP is shown in Figure 2. This shows that the analytic approximatations in the fast and slow regime agree very well with the exact answer, and the transition between the two regimes is rather sharp.

The response and scaling : ECP
The investigation of the scaling behavior for the ECP case (4.8) follows along lines similar to CCP.

Slow Quench Regime
In the slow quench regime (ω 0 δt ≫ 1) one expects a Kibble Zurek scaling. For the protocol (4.8) the Landau criterion determining the Kibble-Zurek time t KZ becomes 1 For ω 0 δt ≫ 1 the solution can appear only at late times. This yields τ KZ ∼ δt log(ω 0 δt) (6.2) In this case the instantaneous gap vanishes in the infinite future. This means that in the slow quench regime adiabaticity will fail at late times. The frequency at this time is Therefore the standard Kibble Zurek argument would predict that the response at late times is given by At a time earlier than the Kibble Zurek time i.e. when τ < δt log(ω 0 δt), (6.5) the adiabatic approximation is valid. Therefore if one measures the response at some fixed value of τ /δt = ζ we should have This expectation needs refinement. Using the exact solution we can perform an expansion for ω 0 δt ≫ 1 and for τ ≫ δt log ω 0 δt. We find that the leading term of ρ 2 (τ ) is The additional logarithmic dependence is not easily visible from the naive Kibble-Zurek argument.

Fast Quench Regime
In the fast quench regime one can get an analytic expression (ω 0 δt) 3 τ. (6.8) at late times, i.e. ω 0 δt ≪ ω 0 τ ≪ 1. Details of calculation which leads to (6.8) are summarized in Appendix A.2. The limit δt → 0 is smooth. In this limit the expression (6.8) reduces to the result which is obtained in an abrupt quench where the frequency suddenly changes from ω 0 to zero, In relativistic theories this limit is non-trivial because of UV divergences, as discussed in [10].
Once again the answer should be obtainable by a perturbation expansion in ω 0 δt. Again let δω(τ ) 2 x > 1. (6.10) Then at late times, Thus ω 0 τ ≪ 1 the perturbation expansion gives a good approximation.

The exact response
The above discussion shows that for the ECP it is useful to look at the response for a fixed value of τ /δt = ζ. Our analytic approximations then predict : ω 0 δt ≫ e ζ Figure 3 shows how the exact result compares with the above expectations. Here we plot the quantity ρ 2 /δt = 2 O /(N 2 δt) as a function of ω 0 δt for different values of ζ. For very small ω 0 δt one reproduces the abrupt quench result. For slightly larger ω 0 δt we can see the fast quench correction predicted in (6.8). To investigate the behavior in the fast quench regime, it is useful to subtract the abrupt quench response. The quantity |ρ 2 (τ ) − ρ 2 abrupt (τ )|/δt is plotted in Figure 4. This quantity is close to zero (and slightly negative) for sufficiently small ω 0 δt. For larger ω 0 δt this becomes positive and in a reasonable range of ω 0 δt this is consistent with the (ω 0 δt) 3 term in the fast quench response, equation (6.8) which are shown by solid lines. Note that the cusps in the data appear because the quantity ρ 2 (τ ) − ρ 2 abrupt (τ ) changes sign and we are plotting the absolute value -there is nothing singular here.
For sufficiently large values of ω 0 δt this quantity is proportional to 1/(ω 0 δt) with a proportionality constant which depends on ζ, as expected from an adiabatic response. There is a small window in the intermediate regime where ρ 2 (τ )/δt is roughly constant upto a logarithmic dependence as in (6.7). Figure 3: (Colour online) The response ρ 2 (τ )/δt as a function of ω 0 δt for ECP. The dots are the exact results obtained by using (4.9) for fixed values of ζ = τ /δt = 0, 2,4,6,8,10,12 which are colored from red to blue respectively. The grey dot on each curve corresponds to ω 0 δt = e ζ for that particular ζ. Thus all points in the yellow shaded region are in the adiabatic regime. The points which lie in the blue shaded region have 1 < ω 0 δt < e ζ . For larger values of ζ there is a small window in this regime where ρ 2 (τ )/δt is roughly constant which is the expectation from Kibble Zurek scaling. The slight increase is consistent with the logarithmic term in (6.7). The dark red and dark blue solid lines are the linear fitting (log y = P log x + Q) results of red (τ /δt = 0) and blue dots (τ /δt = 12) when ω 0 δt ≫ e τ /δt (yellow region), respectively. Both the slopes P are approximately −1. The orange, blizzard blue and light blue solid curves in the fast quench regime (ω 0 δt ≪ 1) are the sudden quench result (6.9) for τ /δt = 2, 6, 10, respectively. For ω 0 δt < 1 the data points lie on these solid lines. For ω 0 δt > 1 they continue to lie on the solid lines for a while and then depart from them, reflecting the O(ω 3 0 δt 3 ) terms in (6.8). as a function of ω 0 δt for ECP. The dots are the exact results obtained by using (4.9) for fixed values of ζ = τ /δt = 0, 2,4,6,8,10,12 which are colored from red to blue respectively. The vertical gridline ω 0 δt = 1 is the threshold between fast quench and slow quench. The dashed lines are a set of cubic functions y = ax 3 , where a = 10, 45, 80, 115 from the lowest one to the highest one, respectively to compare with the leading term in (6.8).

Entanglement Entropy
In this section we present the results for the entanglement entropy of a subregion, its scaling at early times and the time dependence at late times. As argued above, the entanglement entropy in a given subregion for a time dependent frequency can be expressed entirely in terms of the entanglement entropy of a scaled subregion for the system at fixed unit frequency, with the scaling factor given by ρ(τ ) (eqn (3.15)). In the following we will examine the behavior of the entanglement entropy for a subregion −a ≤ x ≤ a. We will also be interested in the limit N ≫ 1 so that we can use the expression (3.11).
We will be interested in the entanglement entropy for a subregion size For ECP the function ρ(τ ) monotonically increases with time, so this condition is equivalent to the condition √ ω 0 a ≪ 1 since ρ(−∞) = 1 √ ω 0 -the monotonicity then implies that once we impose (7.1) at the initial time, this will continue to hold for all times. For CCP the function ρ(τ ) oscillates roughly around ρ(−∞) with an amplitude which is roughly constant in time : once we pick a value of a such that this condition is satisfied at some sufficiently large time, this will continue to be satisfied for all times.
The expression for entanglement entropy at large N can be written down using (3.11, 3.12) and (3.15) by using the Christoffel-Darboux formula for orthogonal polynomials This leads to where the notation A P × A P means that the integrals go over the range defined by A P . These expressions simplify in two regimes. First consider the regime Then one gets, using (3.11) where γ E is Euler's constant. The derivation (7.6) is given in Appendix B.
The logarithmic dependence on the subsystem size is characterisic of 1 + 1 dimensional systems. For relativistic systems the scale is provided by a UV cutoff. For free non-relativistic fermions on a line the entanglement entropy is finite with the UV cutoff replaced by N [46]. A similar result holds for fermions in an invererted harmonic oscillator potential [46,47]. For fermions in a harmonic oscillator potential with a constant frequency this logarithmic dependence has been shown in the so called "bulk limit" in [48].
For CCP protocols, ρ(τ ) oscillates and the condition (7.5) continues to hold once it is imposed at early times. However for the ECP ρ(τ ) monotonically increases so that at very late times the condition 1 √ N ≪ a ρ(τ ) will be violated. It turns out, however, that for the regime a ρ(τ ) one can use a different approximation which yields Note that the entanglement entropy is now proportional to a. However the proportionality constant decreases steadily as 1 τ since the function ρ(τ ) ∼ τ at late times. The derivation of (7.8) is given in Appendix B.
Plots of the time dependence of the entanglement entropy in various cases are shown in Figure 5. For the cis-critical protocol, the function ρ(τ ) oscillates after an initial increase, so that the effective size of the interval in the equivalent constant frequency problem also oscillates. This would lead to oscillations in the entanglement entropy as well.
For ECP, however, ρ(τ ) decreases continuously. This means that for a given a the effective value of the interval in the equivalent constant frequency problem keeps decreasing with time. This should also mean that the entanglement entropy keeps desceasing with time. This is basically because as the fermions are released from the trap they simply spread out : both N A and (∆N A ) 2 keep decreasing leading to a loss of entanglement. It follows from (7.8) that at late times the entanglement entropy goes to zero as a power law ∼ 1 τ .

Phase Space Density for Harmonic Oscillator Potential
In this section we present the time evolution of the Wigner Distribution function, also called the phase space density, u(x, p, τ ), under CCP and ECP quench protocols in a right side up harmonic oscillator potential. In the classical limit, which is given by → 0, N → ∞ with N = fixed, u(x, p, τ ) can only take values of 0, 1 since no two fermions can occupy the same position and momentum. A value of 1 corresponds to the presence of one fermion within a phase space volume between q and q + dq and p and p + dp.
Using the canonical transformations in (2.22), we can write the time evolved phase space density as a function of the original coordinates as (e) ω0δt = 10(CCP) (f) ω0δt = 10(ECP) Figure 5: Time evolution of Entanglement Entropy S A (τ ) in various cases. Red dots are exact large N result from (7.4). Blue solid lines are results (7.6) in the regime (7.5). N = 50, a = 1.
Here E f is the fermi level and defines the boundary of the phase space density, i.e. the fermi surface. θ is the Heaviside step function which satisfies the relations Equation (8.1) takes a value of 1 for q, p which satisfy the relation 2E f ≥ ( q ρ ) 2 + (pρ − qρ) 2 . This will produce what we call a phase space 'droplet'. As time evolves, the shape of this 'droplet' will evolve according to the chosen quench protocol. We present the results for the ECP and CCP cases.

ECP case
Here we discuss the time evolution of (8.1) for the ECP case that has a ρ which is given in (4.3) and (4.9).
In Figure 6 we see that the phase space 'droplet' spreads out in the upper right and lower left quadrants. This corresponds to motion along both directions of the infinite line. Since we are quenching to zero potential, we are 'freeing' the fermions from the harmonic trap and they begin to spread over the real line. The rate at which the 'droplet' spreads is related to δt, the timescale of the quench protocol.

CCP case
Here we discuss the time evolution of (8.1) for the CCP case that has a ρ which is given in (4.3) and (4.6).
In Figure 7 we see that the phase space 'droplet' initially spreads out and then begins to rotate in a clockwise fashion. This rotation comes from the oscillatory nature of ρ in the CCP case for τ > 0. We can understand the physical origin of this rotation. We are quenching from a potential of frequency ω 0 , to 0, back to ω 0 over a timescale of δt. The fermions initially just spread along the real line as the potential barrier goes to zero just as in the ECP case. However, when the barrier is restored to its original value, the fermions hit the edge of the restored barrier and then reflect back. This reflection is indicated by the rotation of the stretched droplet in a clockwise fashion. As time evolves the stretched droplet will continue to rotate indefinitely as the electrons keep reflecting off the walls of the potential barrier.

Time evolution of perturbations along fermi surface
In the previous subsection, we demonstrated the time evolution of a phase space 'droplet' under the influence of a right side up harmonic oscillator potential with a time dependent frequency. In this subsection we consider the time evolution of a perturbation of the fermi surface of this 'droplet'. We would like to know how this perturbation evolves in time. To gain a better understanding of what happens in this case, let us first consider the time evolution under a harmonic oscillator potential with a time independent frequency. In figure 8 we plot this evolution. As expected, we find that the perturbation maintains its shape throughout all Figure 6: Time evolution of a contour plot of the Wigner Distribution function in the classical limit for the ECP case. The black region corresponds to u = 1 and the white region corresponds to u = 0. We have taken δt = 1 , ω 0 = 1. The radius of the initial droplet is 2E f = 2 √ 2 and the area, which is conserved in time, is N = 2πE f . of time. This is a consequence of the harmonic oscillator frequency being time independent. As a result, all points of an initial perturbation of the fermi surface will move at the same angular frequency for all subsequent times leaving its shape unaltered. Now consider the case where an initial perturbation of the fermi surface of a phase space 'droplet' evolves under the influence of a right side up harmonic oscillator potential with a time dependent frequency. In particular, we consider the ECP quench protocol. We plot this Figure 7: Time evolution of the Wigner Distribution function in the classical limit for the CCP case. We have taken δt = 1 , ω 0 = 1. The radius of the initial droplet is 2E f = 2 √ 2 and the area, which is conserved in time, is N = 2πE f . evolution in Figure' s 9, 10. In this case we find something quite interesting. We see that the perturbation develops what we call a 'fold'. This is a phenomenon in which a phase space point which is further from the fermi surface moves faster than a phase space point which is closer to the fermi surface. As a result, at some time later than the initial time, the outer most phase space points begin 'folding' over towards the fermi surface.
The feature of an initial perturbation developing a fold is characteristic of a system evolving under the influence of a harmonic oscillator potential with any arbitrary time dependent frequency. One can rewrite the phase space coordinates q, p in terms of polar coordinates r, θ. One can then show that dθ dτ ∝ f (θ, ω(τ ),ω(τ )) and is therefore not constant in time. On the contrary, if dθ dτ = const, then all the phase space points rotate with the same angular frequency. This is exactly the case for the harmonic oscillator potential with a time independent frequency.

Discussion
In this paper we considered quantum quench in a nonrelativistic field theory of fermions in an external harmonic oscillator or an invererted harmonic oscillator potential with time dependent mass and frequency. While the strategy we outlined to obtain exact solutions hold for both these potentials, we gave results for the right side harmonic potential in this paper. Explicit solutions for the inverted oscillator potential, which corresponds to quantum quench in the Matrix Model description of two dimensional string theory, will be presented in a future publication [30]. We examined scaling behavior of observables in the slow and fast quench regime. We found that the slow quench scaling is consistent with Kibble Zurek, and the fast quench scaling is a result of perturbation theory. This system is, however, not suitable to explore if there is a universal fast quench scaling. For the latter we would need to examine a translationally invariant system with an upper bound on the energy spectrum (for example a lattice system) so that a Lieb Robinson bound is possible. We are currently investigating the quench problem in situations like this.

A Approximation of ρ(τ ) 2 in various limits
In this appendix we explicitly derive approximated ρ(τ ) 2 and therefore O in various limits from the exact CCP solution (4.6) and ECP solution (4.9). In appendix A.1 we study the CCP case and in appendix A.2 we study the ECP case.

B Entanglement Entropy
In this appendix we explicitly derive the approximated Entanglement Entropy (7.6) and (7.8).
In appendix B.1 we figure out N A and in appendix B.2, we figure out A P ×A P dxdy|C(x, y)| 2 .
s.t. N A has similar form to A P ×A P dxdy|C(x, y)| 2 in (7.4). One can easily prove that (B.1) and (7.4) are identical.
In the large N limit, the Hermite polynomial shows the following asymptotic behavior e − x 2 2 · H n (x) ∼ 2 n √ π Γ n + 1 2 cos x √ 2n − nπ 2 (B.2) We use this to simplify the integrand on the RHS of (7.4) or (B.1), Now, change the variables of integration by defining and we obtain (B.5)