Momentum-space entanglement after smooth quenches

We compute the total amount of entanglement produced between momentum modes at late times after a smooth mass quench in free bosonic and fermionic quantum field theories. The entanglement and R\'enyi entropies are obtained in closed form as a function of the parameters characterizing the quench protocol. For bosons, we show that the entanglement production is more significant for light modes and for fast quenches. In particular, infinitely slow or adiabatic quenches do not produce any entanglement. Depending on the quench profile, the decrease as a function of the quench rate $\delta t$ can be either monotonic or oscillating. In the fermionic case the situation is subtle and there is a critical value for the quench amplitude above which this behavior is changed and the entropies become peaked at intermediate values of momentum and of the quench rate. We also show that the results agree with the predictions of a Generalized Gibbs Ensemble and obtain explicitly its parameters in terms of the quench data.


Introduction
A quantum quench is one of the simplest protocols to put a quantum system away from equilibrium. The typical setup is to prepare a state at some initial time slice (e.g., the ground state of a Hamiltonian H 0 ) and suddenly let it evolve in time with a different Hamiltonian H 1 which acts on this state in a non-trivial way. The problem can be equivalently formulated as that of a time-dependent Hamiltonian H(t) that changes from H 0 to H 1 as one of its parameters change in time. This has the advantage of allowing smooth transitions that happen within a finite time scale δt rather than instantaneously. In any case, one is usually interested in the dynamics of various physical observables such as correlation functions and entanglement measures during the process.
The study of quenches provides a window to explore a number of important questions concerning the non-equilibrium quantum dynamics, such as the mechanism underlying thermalization (or not) of isolated systems [1]. It has attracted increasingly more attention due to recent developments in cold atom physics that made possible to experimentally probe the real-time dynamics following a quantum quench [2]. A particular situation of interest corresponds to quenches that cross a critical point, in which case the time evolution of observables may inherit universal properties associated with the critical point. One of the first observations made in this direction was the so-called Kibble-Zurek scaling [3], namely, if the quench rate δt is much larger than any other scale available (i.e., slow quenches) a universal scaling behaviour appears for observables like correlation functions and the entanglement entropy [4].
A special class of quenches that admits many exact results is that of mass quenches in free quantum field theories. The simplest example is that of a free massless boson with action where m(t) 2 is a smooth function that asymptotes to constant values m 2 in and m 2 out at t = −∞ and t = +∞, respectively, changing in time during a scale of roughly δt. This is to be seen as a simple representative of the class of generic quenches of scaling dimension ∆ operators O ∆ (x) in a conformal field theory (CFT) (here the massless scalar CFT, with O ∆ = φ 2 and ∆ = d−2 2 ). Smooth mass quenches of this kind have been recently studied in the literature since certain mass profiles m(t) are amenable to analytical solution for any quench rate δt. The main reason behind that is the observation that the problem can be equivalently understood as the one of a standard (constant mass) quantum scalar field placed in an cosmological background -more precisely the Friedmann-Lemâitre-Robertson-Walker (FLRW) spacetime. The latter has been studied back in the 70's by Bernard and Duncan [5,6], where the mode solutions of the model were obtained and its quantization was performed, so the results together with the intuition from quantum field theory in curved spacetime can be adapted to understand the mass quenches above. This approach has been used in [7,8] to study the behaviour of 1-point functions following the quench. In particular, for fast (but smooth) quenches, i.e., the ones where δt is small compared to all the other physical scales in the theory (but still parametrically larger than the UV cutoff), they found a new universal scaling behaviour for the early time response of renormalized 1-point functions [7]. This scaling had been previously found by the same authors in holographic (hence strongly coupled) theories [9] and was then argued to be a general feature of relativistic quantum field theories. The transition from slow to fast mass quenches as a function of δt was analyzed in [10] for free scalar and fermion theories, showing that there is a smooth crossover to the Kibble-Zurek scaling as the quench rate becomes slow. For a recent review on scaling laws in quantum quenches, see [4].
Apart from local quantities such as one-point functions, there are non-local quantities whose behaviour during a quantum quench is also of interest given that they can probe the dynamics at different distance scales. These include for instance the entanglement entropy (EE) between two subsystems, S = − Tr(ρ log ρ) , where ρ is the reduced density matrix of a subsystem, or the closely related one-parameter family of Renyi entropies with 0 < q < ∞ (q = 1), which includes the EE as the limiting case q → 1. Their time evolution has been thoroughly studied in the literature for the case of spatial subregions (i.e., where the system is split in such a way that the degrees of freedom live either inside or outside a given spatial region) in a variety of models using CFT techniques [11], numerical simulations [12], and holography [13], [14]. The study of scaling properties of the EE as a function of the quench rate δt was recently initiated in [15] for the harmonic chain, where the authors found consistence with the fast quench scaling behaviour of correlation functions and with Kibble-Zurek scaling in the appropriate regimes.
In the present paper, we focus instead on momentum-space entanglement (for previous work on momentum space entanglement in quantum field theory we refer the reader to [16,17,18]). That is, we divide the Fock space in terms of positive and negative-momentum modes of the quantum field and calculate the entanglement production due to the quench between a single mode and all the others (though we shall see that only modes carrying opposite momenta are actually entangled). We investigate mass quenches in both free scalar and free fermion theories. The work follows closely the logic of [19] (see also [20]), which discussed the entanglement produced between momentum modes due to the expansion of the universe. There, in order to have access to exact mode solutions to the equations of motion, the authors had to resort to a unphysical cosmological scale factor which asymptotes to constant values at early and late times (whereas scale factors tipically found in cosmology are of power law type). Fortunately, as we shall see, this weakness becomes actually an advantage from the point of view of quenches, where well-defined asymptotic states are required before and after the quench. In particular, we will show that, starting from the ground state of the pre-quench Hamiltonian at early times (which is a product state of positive and negative-momentum modes according to a pre-quench observer), field modes of momentum k and −k become entangled at late times as a result of the mass quench. This entanglement can be quantified by the Renyi and entanglement entropies.
An important aspect is that, unlike the majority of interacting quantum systems which are known to thermalize after a quench in the sense of approaching a thermal (or Gibbs) ensemble at late times, integrable systems such as the free field models studied here do not reach thermal equilibrium in this usual sense. This is because they possess an infinite number of conserved charges that conspire to constrain the dynamics, a fact that was first observed experimentally in [21] using a system of Bose gases in one dimension. However, later in [22] it was proposed that integrable systems do thermalize in a broader sense to a new kind of equilibrium state called the Generalized Gibbs Ensemble (GGE) (see, e.g., [23] for a review and validity checks for a variety of 1D systems). Having that in mind, we will also show that our results for the entanglement production precisely match the GGE prediction and, moreover, that the parameters characterizing this ensemble can be expressed in closed form in terms of the parameters defining the quench protocol.
The paper is organized as follows. In Section 2, we review the exact mode solutions and quantization of a free massive scalar field in FLRW spacetime and show how a conformal rescaling provides the solution to a mass quench in flat spacetime. In Section 3 we calculate the total amount of momentum space entanglement produced by the quench at late times by following the same logic used in [19] for the curved space picture. Section 4 presents the generalization for a fermionic field, while Section 5 compares the results with the prediction using a generalized Gibbs ensemble and shows an explicit expression for it in terms of quench parameters. Finally, Section 6 contains the closing remarks.

Smooth quantum quenches from FLRW fields
We begin by reviewing the connection between quantum quenches and fields in an expanding spacetime initially raised in [7,8]. A simple way to think of it is to start with the action for a real scalar field with mass m conformally coupled to a curved background metric g µν (x) in any number d of spacetime dimensions, Here, R = R[g] is the Ricci scalar curvature associated with g µν . This gives rise to the equation of motion For the special choice of parameter it is well-known [24] that the action in the massless case (m = 0) is invariant under the Weyl rescaling At the level of the equation of motion this property is nicely summarized by the following operator identity which shows that when acting on a scalar field that transforms simultaneously as shown in (7) the equation is left invariant. Such a symmetry is particularly useful for the purposes of quantization of the field when the curved background g µν (x) is a conformally flat spacetime, i.e., g µν (x) = Ω(x) 2 η µν , since one can use a Weyl rescaling to map to the Minkowski metric η µν and then resort to intuition from canonical field quantization in flat spacetime.
The conformally flat background we shall be interested in is the FLRW spacetime with vanishing spatial curvature. In terms of the conformal time t the metric can be written as g µν (t) = a(t) 2 η µν , i.e, which describes a spatially homogeneous and isotropic geometry that expands according to the scale factor a(t). Unlike the massless case discussed above, the massive scalar field action (4) does not remain invariant under Weyl rescalings. In fact, by introducing the new scalar field φ according to x) the action can be rewritten in terms of the flat spacetime metric as (there is no Ricci scalar contribution since This is nothing but a quantum quench m(t) in the mass of a scalar field in flat spacetime. Hence we see that the problem of mass quenches in flat space can be recast as the "dual" or equivalent problem of a (constant mass) conformally coupled quantum scalar field ϕ living in a curved FLRW background. 1 Of course such a relation is not a peculiarity of mass quenches. It is a straightforward exercise to check that a general interaction term ∼ λϕ n (n ≥ 2) in the action (4) in FLRW gets mapped under the same Weyl rescaling to a quench ∼ λ(t)φ n of the transformed scalar field that lives in flat space, with λ(t) ≡ λ a(t) n+(2−n)d/2 . We are dealing here with the special case of n = 2 only for practical reasons, since the equation of motion is linear.
In flat spacetime, the canonical procedure for quantizing fields [24] begins with first finding the positive-frequency normal modes u k (t, x) ∼ e i (k·x−ωt) that solve the Klein-Gordon equation in order to express the quantum field φ in the basis formed by u k and its complex conjugate. Such a meaningful classification into positive-frequency modes, however, is only possible because flat spacetime admits a timelike Killing vector field K ≡ ∂ t whose corresponding conservation law guarantees well-defined "energy" eigenvalues at any time, i.e., i ∂ t u k = ω u k with ω ≥ 0. In FLRW spacetime this timelike isometry K is clearly not present in general since the metric (9) depends explicitly on time, so the very first step of the quantization procedure seems to fail once we go beyond the flat space scenario. This is a general feature of curved spacetimes that can be boiled down to the diffeomorphism invariance of general relativity, namely the inexistence of a preferred time coordinate to which a meaningful notion of energy can be associated.
Fortunately, this problem can be worked around at least when the conformal factor a(t) 2 in (9) asymptotically approaches constant values as t goes to ±∞. In this case a timelike Killing vector emerges asymptotically in the past and future infinity and one can still make sense of the quantization procedure (i.e., asymptotic positive-frequency modes, the vacuum state, particle excitations, and so on). This particular subset of FLRW spacetimes will be the one of interest in the following. Besides this technical reason, incidentally it turns out to be a useful toy model from the point of view of mass quenches, where well-defined initial and final equilibrium states before and after the quench are required. It is important to stress though that it has little relevance from the perspective of cosmology, where none of the relevant scale factors a(t) happens to saturate to constant values.
We will be interested in the exactly solvable massive scalar field model of [5,6] defined on a FLRW background with scale factor a(t) 2 = A+B tanh t δt , which flows between two asymptotic values A±B at t = ±∞. Equivalently, as discussed above, after a Weyl rescaling the model corresponds to another scalar field in flat space subject to the mass quench that smoothly interpolates (during a time scale of roughly δt) from an initial "pre-quench" value m 2 in at t = −∞ to a final "post-quench" value m 2 out at t = +∞. Here the connection to the original curved spacetime parameters is given simply by The extreme limit of δt → ∞ corresponds to an infinitely slow or adiabatic quench, while δt → 0 would correspond to a step function or instantaneous quench profile of the type discussed in [25] (but the comparison here is subtle and must be done carefully, as discussed in [26]). For m 2 out = m 2 in (or B = 0) we recover the static case of no quench (or no cosmological expansion). The Klein-Gordon equation to be solved is −∂ 2 t +∂ 2 x +m(t) 2 φ = 0. On the grounds of translation invariance in the spatial directions one can seek for normal modes of the Fourier form u k (t, x) = (2π) −(d−1)/2 e i k·x χ k (t), where k ∈ R d−1 is the spatial momentum. 2 The mode functions χ k (t) are easily shown to satisfy a simple harmonic oscillator equation with a time-dependent fundamental frequency, namely Different solutions can be found depending on the boundary conditions imposed. The solutions behaving as positive-frequency ∼ e −i ω in t at early times t → −∞, which we shall refer to as "prequench" modes or simply "in" modes (borrowing terminology from the curved space description), were originally obtained in [27] and are given by is the hypergeometric function and we have introduced the shorthand notation 2 Sometimes for a better handling of UV divergences it is convenient to restrict k to a finite range such as a torus S d−1 of length L (i.e., taking k i ∈ [0, L] with periodic boundary conditions at the two extrema). We will be implicitly doing this in the following when writing for instance Kronecker deltas δ kk instead of Dirac ones, which serves well our physical purposes here without unnecessary technical complications of dealing with infinities.
One can also find another similar set of mode functions that behave as positive-frequency ∼ e −i ωoutt at late times t → +∞, the "post-quench" or simply "out" modes, namely Both constitute equally good basis functions in terms of which the field φ can be expressed, i.e., so they must be related by a Bogoliubov transformation. The special character of our exactly solvable model ensures that this transformation is of the following block-diagonal form [6] χ in with coefficients α k , β k to be determined. Accordingly, it follows from (18) and (19) that the operators a k , a † k in both bases relate as By working out algebraic properties of the hypergeometric functions, the Bogoliubov coefficients for the present model have been obtained in [6], namely From this it is easy to check the following result that will be useful later. The quantization of φ then proceeds by postulating the standard creation/annihilation algebra for the operators a k , a † k (either the "in" or "out" ones), Annihilation by the operator a in k defines the vacuum state |0 in interpreted as the state with no particles as seen by an observer at t = −∞, or before the quench, i.e., Acting on it with a † in k we can construct one-particle (as well as multi-particle) states of the "in" type carrying momenta k and energy ω in . Similarly, annihilation by a out k defines another vacuum state |0 out that now has no particles according to an observer at t = +∞ or after the quench, from which we obtain particle states of the "out" type carrying momenta k and energy ω out through the action of a † out k . The crucial fact here with no analog in the constant mass case is that even the notion of particles in the present case is time-dependent, not universal. From the equivalent point of view of quantum fields in FLRW spacetime this is just the well-accepted statement that the definition of particles is observer-dependent in a curved space. Intuitively, this property is already manifested above in the fact that particle excitations experienced by the two different asymptotic observers carry different energies ω in = ω out , or by the fact that |0 out is not annihilated by a in k (and vice-versa). Indeed, a no-particle state for one observer can look like a complicated particle bath as told by the other. For instance, the number of particle excitations carrying momentum k counted by an asymptotic observer at past infinity using the number operator N in (20)) it is non-zero in the "out" vacuum, This phenomenon is often referred in the general relativity literature as particle production by the gravitational field. In the quench picture it states that, from the point of view of the initial configuration, the full process of quenching the scalar field mass from m 2 in to m 2 out has produced |β k | 2 particles with momentum k, with |β k | 2 given in (22).
By using the mode expansion (18) one can easily show that the momentum space Hamiltonian describing the mass quench above reads where χ k here can be either the "in" or "out" modes, We note that, in the presence of the quench m(t), in addition to the usual particle number piece ∼ a † k a k the Hamiltonian contains also interaction terms between opposite momentum modes k and −k. Such terms disappear when m(t) = const, in which case χ k (t) = 1 √ 2ω k e −i ω k t and H reduces to its standard form H = d d−1 k ω k a † k a k + E 0 . Equation (27) suggests that it might be interesting to split the Fock space into k and −k modes and calculate the entanglement properties of this bipartite system, which shall be done in the next section.

Entanglement production by the quench
The goal of the present section is to study the total amount of entanglement produced between scalar field modes by the mass quenching process introduced above. The discussion follows closely the one of [19], originally done in the context of quantum field theory in curved spacetime in order to estimate the entanglement production due to a cosmological expansion. In what follows, we will be working in the Heisenberg picture, so states are not supposed to change in time.
Let us first give a qualitative picture of what is going on here before delving into the calculations. We have seen in last section that there are two equally good bases in which one can express the field, one adapted to early and the other to late time observers. In the cosmological problem, both observers use the same time coordinate t (defined by the metric (9)) since its tangent vector provides a timelike Killing vector adapted to each of them. Given this time coordinate, both observers can then define a local notion of particle excitation and vacuum, which is valid only at either the asymptotic past or future. That is why, even by working in the Heisenberg picture, we are able to talk about particle and entanglement production after the expansion of the universe.
The same reasoning translates directly to the problem of quantum quenches, which does not involve curved spacetimes at all (it is defined in flat spacetime) but we know to be equivalent to a FLRW field. Now, both the pre-quench and post-quench observers use the same time coordinate t associated to some inertial reference system. But if an experimenter preparing the state at early times is to have a meaningful notion of particles, he or she must use the in-modes, while an observer who will analyze the system at late times, after the quench is finished, naturally picks the out-modes. As a result, we will now show that an initial product state (with respect to positive and negativemomentum bipartitioning of the Fock space) is seen by the post-quench observer as entangled. This entanglement production is obviously tied to our time-dependent Hamiltonian (even though the states do not evolve in time), since the in and out solutions are obtained from it.
We now provide the details for the argument above. As suggested by (27), opposite momentum modes k and −k provide a natural splitting of the Fock space which can be used to explore entanglement properties of the model. For simplicity it is assumed that the initial state is the vacuum |0 in of the Hamiltonian with mass m 2 in , which can be thought of as the state 3 having zero excitation number in any of the momentum modes (as told by the pre-quench observer). 4 Since the Bogoliubov transformation (19) only mixes opposite momenta k and −k, each fixed k piece of (28) admits the following Schmidt decomposition in terms of the out basis where the coefficients are all real and n out labels the excitation number according to a post-quench observer. An explicit expression for the c n can be obtained by applying (20) to (29), i.e., (notice from (21) that α −k = α k and β −k = β k ) which is solved by The remaining coefficient c 0 is fixed by requiring (29) to have unit norm, namely Therefore we see that the in-vacuum (28), which was perceived at initial times as a simple product (i.e., unentangled) state of opposite-momentum modes, is seen by a post-quench observer as a highly entangled state made of infinitely many particle excitations of the out type, with the Schmidt coefficients determined by the Bogoliubov ones α k , β k .
We can now proceed to quantify how much entanglement has been produced by the quench. For simplicity we focus on a single pair of modes (k, −k) with a specific k, since opposite-momenta are the only case allowed by the Bogoliubov mode mixing and the result for multiple modes is easily obtained from this. The density matrix for our bipartite system of opposite momentum modes is therefore 5 This will be the state at any time and the physical properties of pre-or post-quench stages are entirely encoded in which basis of the Hilbert space we decide to use. In order to obtain the entanglement produced between the modes at the end of the process, we need to trace out with respect to the out-basis the negative-momentum modes to obtain the reduced state describing the positive-momentum ones. Notice that had we traced out with respect to the in-basis we would have obtained the unentangled state ρ k = |0 in k 0 in k |. The Rényi entropies (3) can be used to quantify the amount of entanglement between the modes k and −k. They are constructed from the reduced state ρ k as From (31) and (32) it is seen that |c n | 2 depends on the index n essentially as a power law, i.e., where we have introduced the parameter 6 and used (22) in the second equality to put it into a compact form. The summation over n thus becomes just a geometric series which is easily carried out to yield the following closed-form expression for the Renyi entropies In particular, taking the limit q → 1 gives another important entanglement measure -the entanglement entropy S k = − Tr {k} ρ k log ρ k , namely This agrees exactly with the expression originally obtained in [19] in the context of a scalar field in a cosmological setup.
It is important to notice that tracing out the positive momentum modes instead of negative ones in (34) would yield the same Renyi and entanglement entropies as long as the full density matrix 5 The proper way to construct it is by tracing out all momentum modes in the initial vacuum (28) except k, which is a trivial operation giving an overall factor of unity for each mode and yielding (33) at the end. 6 The subscript B stands for "boson", to be contrasted with the fermionic case in the next Section.
ρ is a pure state, since in this case the Schmidt form (33) guarantees that the two reduced density matrices ρ k and ρ −k have the same eigenvalues |c n | 2 . The entropies (38) and (39) quantify the total amount of entanglement produced at late times by the quench. They depend on the initial and final values of the masses m 2 in , m 2 out (and of course on the selected mode k) via the frequencies ω ± , and additionally on the quench rate δt, but both in a very peculiar way: all combined inside the parameter γ B shown in (37). So in this sense γ B contains all the information concerning the late-time entanglement between opposite-momentum modes.
Both S (q) k and S k are monotonically increasing functions of γ B , so they can be inverted to yield γ B (S), as noticed in [19]. This is rather interesting, meaning that the quench protocol m(t) can in principle be fully reconstructed only from entanglement data. The expression for the EE is too complicated to invert and find γ B (S) analytically, but for the Renyi's the situation is simple enough so that this can be done, namely, all we have to do is solve for γ B the q-th order equation with s q ≡ e (1−q)S (q) k . The second Renyi entropy provides the simplest example, In any case, having obtained γ B = γ B (S) should be enough to express the quench parameters in terms of it using the definition (37). The physical interest, however, is on the k and δt-dependence of the entropies instead of the parameter γ B . Since γ B is symmetric under ω in ↔ ω out one can choose to focus on quenches that decrease the mass, i.e., m 2 in > m 2 out . 7 An interesting special case is that of m 2 out = 0, where the postquench Hamiltonian is that of a conformal field theory (the free massless boson). In the following we shall limit our numerical analysis to the entanglement entropy and the second Renyi entropy, although we have checked that the conclusions hold also for the other Renyi entropies with q = 2.
For a fixed rate δt, the entropies decrease monotonically with k ≡ |k| as shown in Figure 1 for many different combinations of the two masses. This shows that more entanglement is produced between IR (low k) modes than between UV (high k) modes. Also, for each given k the entropies are bigger the larger is the difference δ(m 2 ) ≡ m 2 in − m 2 out between the initial and final masses, the maximal value corresponding to quenches that end in the massless boson CFT (top curves of each color). Actually when the post-quench theory is a CFT there is formally a divergent zero mode contribution ∼ log 1/k at k = 0, but we can simply ignore this fact since in this case there is no (k, −k) splitting of modes at all to begin with. Let us recall that the usual upper bound S EE ≤ log dim(H) for the EE is infinite here since the Hilbert space for the reduced state ρ k is infinite-dimensional, so there is nothing contradictory about the fact illustrated in the plot that the EE is not limited from above. Figure 2 illustrates the dependence of the entropies on the time scale δt. In (a) we show the result for a single mode k, while (b) corresponds to the total entanglement production by the quench obtained by integrating over all k ∈ (0, ∞). Again it can be seen that more entanglement is produced the farther away the initial and final states are from each other (a property that holds at any rate δt). It also becomes clear that faster quenches produce more mode entanglement than slower ones. In particular, as δt grows all curves approach zero asymptotically, indicating that infinitely slow or adiabatic quenches (those with δt → ∞) do not create any entanglement between field modes. 7 When m 2 in = m 2 out , γ B vanishes and there is no entropy production, which is trivial since in this case there is no quenching at all (see (12)

Fermionic case
A completely analogous construction holds for the case of a mass quench of free Dirac fermions. At the level of equation of motion, the problem amounts to solving the Dirac equation with a timedependent mass, subject to appropriate boundary conditions that allow the identification of positive-frequency modes and asymptotic observers to proceed with the quantization of the model. Here, γ µ are the usual Dirac matrices in d-dimensional flat spacetime, defined by the Clifford algebra {γ µ , γ ν } = 2η µν . Just as in the bosonic case, the same equation appears when one analyzes a free fermion with constant mass m placed on a curved FLRW spacetime. Specifically, the Dirac equation for such a fermion Ψ minimally coupled to gravity reads whereγ µ are curved space Dirac matrices (satisfying {γ µ ,γ ν } = 2g µν ) and Γ µ the spinorial affine connection associated to the FLRW metric g µν . By performing the conformal rescaling together with g µν = a(t) 2 η µν we get equation (42) for the transformed fermion ψ that lives in flat spacetime, where the quench profile is related to the cosmological scale factor by Analytical mode solutions to (42) have been worked out in [6] for an exactly solvable FLRW model that we now briefly review with appropriate adaptations in order to translate the results for our case of interest, the problem of mass quenches. The model is characterized by the mass profile which is morally the same as (12) in the bosonic case (but notice that here the tanh profile is for m(t), while in that case it was for m(t) 2 ). 8 Translational symmetry in the spacelike directions allows the ansatz where the function φ k (t) can be checked to satisfÿ Notice the strong parallel with (14). One subtlety, however, is that φ k here has d components and this is a matrix equation. Writing where v(0, λ), u(0, λ) are constant zero-momentum eigenspinors of γ 0 with opposite eigenvalues ±i , it follows that the scalar functions φ (±) k also satisfy a harmonic oscillator equation similar to (14), but now with an imaginary contribution to the time-dependent oscillator mass. The solutions with "in" and "out" boundary conditions at t = ±∞, i.e., behaving as positive-frequency for pre-quench and post-quench observers respectively, are the following where ω in , ω out are the same as previously defined in (16) and we have introduced 8 Equivalently, from the curved spacetime point of view the bosonic mass profile (12) corresponds to a FLRW metric with scale factor a(t) = A + B tanh t δt while the fermionic one (46) to a(t) = A + B tanh t δt (with the cosmological parameters A, B related to the ratios min m , mout m as in (13)).
The general solution for the fermion ψ can then be written as with the sum over the spinor index running up to λ max = 2 d/2−1 for d even and 2 (d−3)/2 for d odd, and the curved space spinor mode solutions Of course there are analogous expressions for the out modes as well.
A special property of this exactly solvable model, as in the bosonic case, is that the Bogoliubov transformation connecting pre-quench and post-quench modes is diagonal in momentum space. In terms of the functions φ (±) k it takes the simple form that allows us to relate the corresponding creation and annihilation operators as Here, andū out (k, λ) ≡ i u † out (k, λ) γ 0 refers to the polarization λ and momentum k spinor The Bogoliubov coefficients above are known analytically [6], namely With the definitions above we can now proceed as in the previous section and calculate the entanglement production between opposite momentum modes of the fermionic field due the mass quench. We focus on d = 1 + 1 dimensions, where the sum over λ in (53) runs over a single value and considerably simplifies the analysis (also, since the spatial momentum k has only one component it can be denoted simply by k). A similar analysis has been carried out in [20] where the idea was to quantify the entanglement production for a fermionic system due to the cosmological expansion of the FLRW spacetime.
The pre-quench vacuum is defined by and the post-quench one |0 out can be defined in a similar way. The Fock space can again be split into opposite momentum modes ±k, e.g. |0 in = ⊗ k |0 in k 0 in −k , and since the Bogoliubov transformation (56) between in and out creation/annihilation operators only mixes k and −k modes one can repeat the steps done previously and express the initial vacuum in terms of the out basis. The result takes the form with Here, |1 out k , 1 out −k ≡ a out † k,λ b out † −k,λ |0 out denotes the state containing a particle with momentum k and an antiparticle with momentum −k as told by the late time observer. Therefore, we see that the initial vacuum is populated by particle-antiparticle pairs of the out type carrying opposite momenta. This should be contrasted with (62) for bosons, in which case there was an infinite tower of multiparticle excitations thanks to the absence of Pauli's principle in that case.
The density matrix corresponding to the in vacuum state is simply ρ = |0 in 0 in | with |0 in expressed in terms of the out basis by (62). By focusing on a particular pair (k, −k) of momentum modes and tracing out all the antiparticles, we get the very simple reduced state where we have introduced the parameter and used (60) to simplify it. This is the fermionic analog of γ B in the bosonic case (see (37)). At this point the technical difficulties of solving the Dirac equation with a time-dependent mass faced above really pay off when it comes to calculating the entanglement properties: the reduced state (64) is just a two-dimensional diagonal matrix (compare with the infinite sum over occupation numbers in (34)). The Rényi entropies of order q are immediately found to be and taking the limit q → 1 we obtain the entanglement entropy which agrees with the expression originally obtained in [20] in the framework of an expanding spacetime. Notice that the result for the entropies would have been the same had we traced out the particles instead of the antiparticles in the beginning, since the total state is pure. The total amount of entanglement produced by the quench at late times is measured by these entropies and depends on the masses m in , m out , the mode k, and the time scale δt. Similarly to the free boson case, it is fair to state that the parameter γ F encodes all the the late-time entanglement properties between opposite momentum modes. The resemblance between the formulas (66),(67) and (38),(39) for the fermionic and bosonic entropy production is striking. In fact, the expressions for the Renyi entropies can be converted into minus one another under γ F ↔ −γ B . Notice that this simple relation between the bosonic and fermionic Rényi entropies does not commute with the limit q → 1, i.e., it is not shared by the expressions for the EE.
The entropies are again monotonic functions of γ F and can be inverted to determine γ F (S), i.e., the information concerning mode entanglement is in principle enough to determine all the quench parameters. An important difference with respect to the bosonic result, however, is that the EE in the present case is limited from above by log 2 ≈ 0.7, since the reduced state for particles with momentum k lives in a two-dimensional Hilbert space. In particular, even when one of the masses vanishes (i.e., when the quench crosses a critical point) this upper bound prevents the occurrence of the divergent zero mode contribution that takes place for bosons. Another important difference is that now the expressions are not symmetric under ω in ↔ ω out , so the behavior for m in > m out is distinct from m in < m out (or δm > 0 and < 0, respectively) and must be analyzed separately, as we shall see. Figure 3 shows the k-dependence of the EE and the second Renyi entropy for a fixed quench rate δt. The lefthand side figure corresponds to quenches that decrease the mass, while the other two to increasing-mass quenches. For the former the entropies are always monotonically decreasing functions of k, showing that entanglement production is more noticeable between IR modes. However, for the case of increasing-mass quenches this is only true up to a critical value of δm (see (b) and (c)). Above this critical value of δm this monotonic behavior is broken as shown in the right figure.
Interestingly, this indicates that for large enough final masses the maximal entanglement production is not achieved at IR modes but rather at an intermediate momentum value k = k max . Figure 4 illustrates the dependence of the entropies on the time scale δt. Cases (a), (b), (c) show the result for a single mode k in a quench with δm < 0 (the former) and δm > 0 (the latter two), respectively, while (d), (e), (f) show the total entanglement produced in each of the two situations, obtained after integration over all k. We see that in mass-decreasing quenches the entanglement production is always bigger the faster the quench is done. For a single mode k, this is still true for mass-increasing quenches up to a limiting value δm * > 0 ((a) and (b)) but fails to be true for deformations stronger than this (case (c)), where maximum production is then achieved at an intermediate value of δt. This unusual feature disappears when one integrates over all modes k, as shown in (d). Notice also that again adiabatic quenches do not produce any mode entanglement, since all curves asymptote to zero as δt grows.

Connection with the Generalized Gibbs Ensemble
In this section, we will show that our above-mentioned results for the Renyi and entanglement entropies agree with the predictions of a Generalized Gibbs Ensemble, showing that this steady state correctly describes the late time dynamics after smooth quenches even from the point of view of entanglement measures (in addition to the well-understood case of correlation functions).
Recall that while most quantum systems are known to thermalize in the usual sense of approaching a Gibbs ensemble at late times, integrable systems such as our free field models thermalize in the more general sense of a Generalized Gibbs Ensemble (GGE). The density matrix for a GGE state is given by [22] ρ where {I j } denotes a full set of conserved charges, Z = Tr e − j λ j I j is the partition function, and {λ j } are Lagrange multipliers, fixed by demanding that the GGE value of the conserved charges coincides with their initial value, The precise statement about the equilibration of integrable systems, therefore, is that expectation values of local operators in the model converge to their GGE ensemble values at late times or, equivalently, that the reduced state of any finite subsystem approaches the corresponding reduced density matrix of GGE. The question remains of which integrals of motion should appear in the GGE for a given model. For free boson theories in 1 + 1 dimensions and Gaussian initial states (such as the vacuum of quadratic Hamiltonians), it was shown in [28] that a good choice for these operators is the particle number N k = a † k a k at each momentum mode (the j then becomes dk). In our case of interest (33), where we deal with a single pair of modes (k, −k) in the out-basis, after tracing over all the other irrelevant momenta to get a bunch of unit factors the GGE is simply where the partition function is .
After tracing out the negative-momentum modes one is left with the reduced GGE density matrix describing the positive ones. The lagrange multipliers λ k are defined by the condition (69), namely Tr (ρ GGE k N out k ) = 0 in |N out k |0 in ≡ n k , which gives Now it is straightforward to check (see (26) for a similar calculation) that the occupation number n k is determined by the Bogoliubov coefficient β k as n k = |β k | 2 , and using the algebraic relation |β k | 2 + 1 = |α k | 2 (see (22)) we get where γ B is defined by (37). This is a remarkable result which is worth emphasizing: the Lagrange multipliers that characterize the GGE are completely determined by the quench parameters m 2 in , m 2 out , δt and the momentum mode k via the parameter γ B . To the authors's knowledge, such a complete construction of the GGE state was unknown in the literature for the case of smooth quenches.
Indeed, if we calculate the Renyi entropy associated to the reduced GGE state (72) we obtain which is in perfect agreement with our result (38) after identifying the parameters according to (74). This shows that the GGE prediction for the late-time dynamics of free bosons as proposed in [28] indeed works for our model (as told by entanglement measures). A similar construction should hold also for the fermionic model, in which case the Lagrange multipliers characterizing the GGE would be fully determined by the fermionic parameter γ F of (65).

Conclusions
In this work, we have calculated the late-time entanglement in momentum space produced after smooth mass quenches for free scalar (in d dimensions) and fermion theories (in 1+1 dimensions). The strategy was to follow [19,20], which studied the entanglement production for quantum fields in an expanding universe, and showing how this can be immediately translated to the problem of a quench by using a Weyl rescaling of the field. The initial state was taken to be the vacuum as defined by a pre-quench observer. As the quench is performed, particle excitations carrying all possible momentum modes are produced as told by a post-quench observer. We then calculated the entanglement and Renyi entropies between a single quantum field mode k and its opposite mode carrying momentum −k, which is the only non-trivial case in our exactly solvable model (i.e., no entanglement is produced between other pairs (k, k ) of modes other than this). The results have simple analytical formulas and are expressed in terms of a single parameter (γ B for bosons and γ F for fermions) that encodes all the late-time entanglement properties.
For (1 + 1)-dimensional bosons, we have shown that our results match the predictions of a Generalized Gibbs Ensemble where the conserved charges are taken to be the mode number operators as defined by the post-quench observer, namely N out k = a † out k a out k (for all k). We were able to precisely calculate the Lagrange multipliers that characterize the GGE in terms of the parameters m in , m out , δt defining the quench, namely λ k = − log γ B . Hence, having fully specified the GGE state, the long-time behavior of any local observable after the quench can now be calculated (not only the entanglement and Renyi entropies presented here).
In the bosonic case, our results show that either for a fixed mode k or integrating over all of them, more entanglement is produced for faster quenches (small δt) with respect to slow ones, while adiabatic quenches (δt → ∞) do not produce entanglement at all. Also, the amount of entanglement grows with the intensity of the quench, characterized by the magnitude δ(m 2 ) = m 2 out − m 2 in of the difference between the initial and final state masses, the result being symmetric under m 2 in ↔ m 2 out . In particular, this means that, among all the quenches ending up with some m 2 out = 0, maximal entropy production is achieved for the one that started from a CFT (i.e., m 2 in = 0). We also showed that, for any fixed quench rate δt, more entanglement is produced between light or IR modes (the ones carrying small momentum k) than between heavy ones, as expected.
In the fermionic case the results are more subtle. First we noted that the sign of δm = m out − m in becomes important, that is, the result here is not invariant under m in ↔ m out and, as a consequence, whether the quench increases or decreases the mass has a significant impact on the final result. For a fixed rate δt, and for a quench decreasing the mass, we find a very similar result to that of bosons. For quenches that increase the mass, however, the results remain similar to that of bosons only up to a certain critical value δm * , above which the entanglement ceases to be a monotonic funtion of |k|. Instead, as we increase the mass beyond this critical value, we observe a peak at some intermediate value k max that grows as δm grows. In other words, for quenches that increase too much the mass it fails to be true the statement that more entanglement is produced between light modes. Also, for a particular mode k, the behaviour of the entanglement and Renyi entropies as a function of δt has the same features discussed in the last paragraph. Namely, for quenches that decrease the mass the result is very similar to the bosonic one, while for quenches that increase the mass this holds true only up to a critical value δm * above which the entropies cease to decrease monotonically with δt, becoming instead peaked at an intermediate value. If we integrate over all modes k to get the total entanglement production the entropies remain monotonically decreasing, although for δm > δm * we can still notice a small bump at an intermediate value of δt.
As future prospects, one would like to perform similar calculations for weakly coupled quantum field theories in order to access how the presence of interactions modify the results. However, subtleties are expected to appear in this context when factorizing the Hilbert space owing to mode mixing, which should be clarified at first place. It would be also interesting to compare the entanglement production with that of instantaneous quenches instead of the smooth ones considered here, since (as noted in [26]) there are subtleties in naively taking the limit δt → 0 here that need to be dealt with properly. Another interesting generalization would be evaluating momentum space entanglement in strongly coupled theories via the AdS/CFT correspondence by generalizing the Ryu-Takayanagi prescription [29] (see also [30]) that is appropriate to deal only with real-space entanglement entropy.