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ényi 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 δt\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta t$$\end{document} can be either monotonic or oscillating. In the fermionic case the situation is more 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 a e-mail: dwfalves@gmail.com b e-mail: gcamilo@iip.ufrn.br 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 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 [3,4], where mode solutions for some specific choices of mass profiles 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 [5,6] to study the behaviour of 1-point functions following the quench (see also [7][8][9] for related work).
Apart from local quantities such as one-point functions, there are non-local quantities whose behaviour during a quan-tum quench is also of interest given that they can probe the dynamics at different distance scales. These include the entanglement entropy (EE) S = − Tr(ρ log ρ), (2) where ρ is the reduced density matrix of a given subsystem, or the closely related one-parameter family of Rényi 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 [10], numerical simulations [11], and holography [12,13]. The study of scaling properties of the EE as a function of the quench rate δt was recently initiated in [14] 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. 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 the logic of [15] (see also [16]), which discussed the entanglement produced between momentum modes due to the expansion of the universe for a specific choice of FLRW scale factor. We will consider the same choice, which in our case translates to massincreasing or mass-decreasing quenches, as well as another class of quench profiles that recovers the initial mass at late times. For previous work on momentum space entanglement in quantum field theory we refer the reader to [17][18][19]; see also [20,21] for similar work with spin chains.
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 [22] using a system of Bose gases in one dimension. However, later in [23] 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., [24] 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. This reinforces the relevance of considering momentum-space entanglement as an interesting probe of thermalization of quantum systems. Unlike its real-space counterpart, the entanglement between momentum modes is not tied to any particular spatial sub-region of the system, so one should expect it to capture different physics characterized by non-local correlations in real space. It is worth to recall that the most accepted explanation for the spreading of real-space entanglement after a quench relies on the quasi-particle picture introduced in [10], where it is assumed that pairs of quasi-particles with opposite momenta created within the correlation length are entangled, while those far apart from each other are not. As these particles travel through the system, real-space entanglement happens to grow ballistically. This intuition is able to correctly reproduce many of the qualitative features of the entanglement dynamics, and has been used to derive new results [25]. From our calculations we can make this concrete by computing exactly how entangled the different particles produced are as a function of the parameters controlling the quench, such as its amplitude and speed. Hence, understanding the dynamics of momentum-space entanglement may also shed light into our understanding of how real-space entanglement grows. The free field example is chosen for convenience, since it is amenable to exact analytical results while still allowing for interesting phenomena such as the approach to the GGE, but we hope that our study can be a useful benchmark in future studies of thermalization in more complicated interacting models.
The paper is organized as follows. In Sect. 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 Sect. 3 we calculate the total amount of momentum space entanglement produced by different quenches at late times by following the same logic used in [15] for the curved space picture. Section 4 presents the generalization for a fermionic field, while Sect. 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, Sect. 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 [5,6]. 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 [26] 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 ϕ(t, x) ≡ a(t) − φ(t, x) the action can be rewritten in terms of the flat spacetime metric as (there is no Ricci scalar contribution since R[η] = 0) where 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 [26] 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 two representative cases where the mass profile m(t) = m a(t) allows for exact solutions, but before particularizing to specific choices of m(t) let us first sketch the general strategy (see next section for explicit expressions in the cases of interest).
The Klein-Gordon equation arising from (10) is 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 ones behaving as positivefrequency ∼ 1 √ 2ω in e −i ω in t at early times t → −∞ are referred to as "pre-quench" modes or simply "in" modes χ in k (t) (borrowing terminology from the curved space description), while the ones that behave as positive-frequency ∼ 1 √ 2ω out e −i ω out t at late times t → +∞ are called "postquench" or "out" modes χ out k (t). 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.
Translational invariance of the model in the spatial directions ensures that this transformation is of the following blockdiagonal form with coefficients α k , β k to be determined. Accordingly, it follows from (13) and (14) that the operators a k , a † k in both bases relate as 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.
The quantization of φ then proceeds by postulating the standard creation/annihilation algebra for the operators a k , a † k (either "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 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 oper- (15)) 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 observer, the full process of quenching the scalar field mass from m in to m out has produced a bath of |β k | 2 particles carrying momentum k.
It is interesting to note that after using the mode expansion (13) the momentum space Hamiltonian describing the mass quench above reads where χ k here can be either the "in" or "out" modes, is the vacuum energy contribution, and ω k (t) 2 ≡ k 2 + m(t) 2 . We see 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, where the modes are simply (20) 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. This shall be done in Sect. 3, but first let us illustrate the results above explicitly for the two cases of interest.

Tanh profile
The first quench profile of interest is which smoothly interpolates (during a time scale of roughly δt) from an initial "pre-quench" value m 2 in at t = −∞ to another final "post-quench" value m 2 out at t = +∞. This is reminiscent of the massive scalar field model studied in [3,4], defined on a FLRW background with scale factor a(t) 2 where the cosmological parameters A and B are related to the quench data 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 [27,28]. For m 2 out = m 2 in (or B = 0) we recover the static case of no quench.
The in and out mode solutions for this system were originally obtained in [29] and are given by is the hypergeometric function and we have introduced the shorthand notation and By working out well-known algebraic properties of the hypergeometric functions that convert the argument z into 1 − z (here z = 1 2 (1 + tanh t δt )), the Bogoliubov coefficients (14) for the present model have been obtained in [4], namely In particular, their absolute values are easily checked to take the simple forms An interesting limiting case is that of abrupt quenches (δt → 0), where the Tanh profile (21) becomes the step function m(t) 2 = m 2 in θ(−t) + m 2 out θ(t) and the in and out modes are simple plane waves with frequency ω in and ω out , respectively. In this case the coefficients above simplify to

Sech profile
The second quench of interest is the Gaussian-like profile which has the same initial and final values of mass m 2 0 and was studied previously in [6]. The in and out solution read [30] χ in/out where we have defined The Bogoliubov coefficients relating in and out solutions are given by [30] We will use these expressions to compute the entanglement and Rényi entropies below.

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 [15]. 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 negative-momentum 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 (20), 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 (14) only mixes opposite momenta k and −k, each fixed k piece of (33) 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 (15) to (34), i.e., (notice from (25) which is solved by The remaining coefficient c 0 is fixed by requiring (34) to have unit norm, namely Therefore we see that the in-vacuum (33), 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. We focus on a single pair of modes (k, −k), since opposite-momenta are the only case allowed by the Bogoliubov mode mixing and the result for multiple pairs of 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 (36) and (37) 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 (26) 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 5 The proper way to construct it is by tracing out all momentum modes in the initial vacuum (33) except k, which is a trivial operation giving an overall factor of unity for each mode and yielding (38) at the end. 6 The subscript B stands for "boson", to be contrasted with the fermionic case in the next section.
closed-form expression for the Renyi entropies In particular, the limit q → 1 gives the entanglement entropy This agrees exactly with the expression originally obtained in [15] 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 (39) would yield the same Renyi and entanglement entropies as long as the full density matrix ρ is a pure state, since in this case the Schmidt form (38) guarantees the two reduced density matrices ρ k and ρ −k to have the same eigenvalues |c n | 2 .
The entropies (43) and (44) quantify the total amount of entanglement produced at late times by the quench. They depend on the quench rate δt, on the mode k, and on the masses m 2 in , m 2 out or m 2 0 for the Tanh or Sech profile, respectively. But the dependence on these four physical parameters is of a very peculiar kind: only combined inside the single parameter γ B defined in (42). In this sense, γ B contains all the information concerning the late-time entanglement between opposite-momentum modes. Its explicit form for the Tanh and Sech profiles introduced in Sect. 2 is readily found from (26) and (31), namely with μ given by (30). 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 [15]. 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, k with q = 2, 3, 4, 5 (from top yellow to bottom purple) produced by the tanh mass quench for a fixed quench rate δt = 0.1. The black dashed curve shows the corresponding particle production rate |β k | 2 discussed in (19). In the plot we have set m 2 in = 1 and m 2 out = 4, but the shape of the plots remains unchanged for other values. The magnitude of the entropies is controlled by the mass difference δm 2 = m 2 out − m 2 in , namely it grows as δm 2 is increased (the precise proportionality law is not known though) In any case, having obtained γ B = γ B (S) should be enough to express the quench parameters in terms of it using the definition (42).
The physical interest, however, is on the k-and δtdependence of the entropies themselves, which we now analyze in detail.

Tanh profile
We start with the Tanh profile. Since γ B is symmetric under ω in ↔ ω out one can choose to focus on quenches that increase the mass, i.e., m 2 out > m 2 in . 7 An interesting special case is that of m 2 in = 0, where the pre-quench 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 first few integer Rényi entropies (q = 2, 3, 4, 5).
For a fixed quench rate δt, all the entropies decrease monotonically with k ≡ |k| as shown in Fig. 1. This shows that more entanglement is produced between IR (low k) modes than between UV (high k) modes. The magnitude of the entropies is proportional to the mass difference δm 2 = m 2 out − m 2 in between the initial and the final states. The maximal value corresponds to quenches that start from the massless boson CFT, although this case is subtle since there is formally a divergent zero mode contribution ∼ log 1/k at k = 0. In practice we can simply ignore this fact since in this case there is no (k, −k) splitting of modes at all to begin with. 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 (21)). with q = 2, 3, 4, 5 (from top yellow to bottom purple) as functions of the quench rate δt. The black dashed curve shows the corresponding particle production rate |β k | 2 discussed in (19). For the plot we fix a single mode k = 0.5 and m 2 in = 0.5, m 2 out = 2, but the behavior is qualitatively the same for other values Let us recall that the usual upper bound S E E ≤ 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 for a single mode k. It 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.

Sech profile
For the case of the Sech profile, we plot the results on Figs. 3 and 4. By fixing δt and m 0 , we see that all Rényi and entanglement entropies decrease monotonically with k. That is, light degrees of freedom are more entangled than heavy ones, similarly to the Tanh profile case. However,there is an interesting difference when fixing k and exploring the δt dependence, as shown in Fig. 4. The entanglement does not decrease monotonically as in the Tanh profile case, but rather oscillates with an amplitude that decreases as δt increases. The origin of this oscillatory behavior can be seen in (45). Also, notice that again there is a divergence as k approaches zero that we should not worry about as we commented above for the Tanh profile.

Entanglement per particle
It is also interesting to consider the ratio of the entanglement production by the number of particles as a function of δt, for k with q = 2, 3, 4, 5 (from top yellow to bottom purple) produced by the sech mass quench for a fixed quench rate δt = 0.1. The black dashed curve shows the corresponding particle production rate |β k | 2 discussed in (19). In the plot we have set m 0 = 2, but the shape of the plots remains unchanged for other values fixed k, and as a function of k for fixed δt. This is presented in the Figs. 5 and 6 below for both the Tanh and Sech profiles. It is important to note that this ratio is not constant, besides our earlier remarks that both quantities behave qualitatively similar as a function of the various quench parameters.

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 time-dependent 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 (48) 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 (48) have been worked out in [4] 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 with q = 2, 3, 4, 5 (from top yellow to bottom purple) as functions of the quench rate δt. The black dashed curve shows the corresponding particle production rate |β k | 2 discussed in (19). The mass is fixed to be which is morally the same as (21) 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ÿ  (22)).
Notice the strong parallel with (12). 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 (12), but now with an imaginary contribution to the time-dependent oscillator mass. The solutions with "in" and "out" boundary condi-tions 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 (24) and we have introduced 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 [4], 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 (59) 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 [16] 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 (62) 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 (68) 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 (68). 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 (66) to simplify it. This is the fermionic analog of γ B in the bosonic case (see (42)). 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 (70) is just a two-dimensional diagonal matrix (compare with the infinite sum over occupation numbers in (39)). 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 [16] 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 late-time entanglement properties between opposite momentum modes. The resemblance between the formulas (72), (73) and (43), (44) 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 7 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 8 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 ≥ 0 (this is the left-right entropy studied, e.g., in [31]). 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 [23] 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 [32] 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 (38), 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 (75), namely Tr ρ GGE k N out Now it is straightforward to check (see (19) 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 (26)) we get where γ B is defined by (42). 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 . This gives support to the results of [33], where thermalization to a GGE was argued to hold at the level of two-point correlators.
Indeed, if we calculate the Renyi entropy associated to the reduced GGE state (78) we obtain which is in perfect agreement with our result (43) after identifying the parameters according to (80). This shows that the GGE prediction for the late-time dynamics of free bosons as proposed in [32] 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 (71).

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 is inspired by [15,16], which studied the entanglement production for quantum fields in an expanding universe, and can be translated to the problem of a quench through 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 Rényi 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, for both the Tanh and the Sech mass profiles, we saw that at any fixed quench rate δt more entanglement is produced between light modes (the ones carrying small momentum k) than between heavy ones, as expected. The entanglement production is monotonically decreasing with k and its magnitude grows with the magnitude of the quench, characterized by the absolute value of the difference δm 2 = m 2 out − m 2 in of the initial and final mass in Tanh case and by m 2 0 (the maximum value of the mass during the quench) in the Sech case. The picture is qualitatively similar to the particle production rate given by n k = |β k | 2 .
The dependence on δt for a given mode k is more interesting: for both profiles it is true that more entanglement is produced for faster quenches (small δt) with respect to slow ones, while adiabatic quenches (δt → ∞) do not produce entanglement at all. In particular, this means that, among all the quenches reaching some m 2 out = 0, maximal entropy production is achieved for the one that started from a CFT (i.e., m 2 in = 0). However, for the Tanh profile the entanglement is found to decrease monotonically as a function of δt, while for the Sech profile the decrease is non-monotonic, being given by damped oscillations modulated by k.
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 the Tanh profile for 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 function 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 Rényi 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.
We believe these results shall hold (at least qualitatively) for generic mass quenches in free field theories, that is, regardless of the precise form of the quench profile m(t). For instance, we expect the same behavior as in the Tanh case for any other quench profile that goes from the same m in to the same m out within a finite scale δt. This expectation is supported even by comparing the qualitative results for the Tanh and the Sech profiles (which do not share the same initial and final masses), namely the fact that the entanglement production is more significant for light modes and faster quenches. Of course we have no formal proof of that, but it would be very surprising if by simply changing m(t) the k or δt dependence of the entropies happened to be qualitatively different.
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. Another interesting generalization would be a geometric prescription for computing momentum-space entanglement in AdS/CFT, by generalizing the Ryu-Takayanagi prescription [34,35]) which is appropriate to deal only with real-space entanglement entropy. This was already speculated on [36] but the problem remains open. In fact, as discussed in [36], one of the main difficulties for putting forward a holographic proposal on momentum-space entanglement is the fact that this quantity is not well studied on the field theory side (specially for interacting theories). Hence, we hope our work can contribute to this issue.