Dynamics of entanglement in expanding quantum fields

We develop a novel real-time approach to computing the entanglement between spatial regions for Gaussian states in quantum field theory. The entanglement entropy is characterized in terms of local correlation functions on space-like Cauchy hypersurfaces. The framework is applied to explore an expanding light cone geometry in the particular case of the Schwinger model for quantum electrodynamics in 1+1 space-time dimensions. We observe that the entanglement entropy becomes extensive in rapidity at early times and that the corresponding local reduced density matrix is a thermal density matrix for excitations around a coherent field with a time dependent temperature. Since the Schwinger model successfully describes many features of multiparticle production in $e^+ e^-$ collisions, our results provide an attractive explanation in this framework for the apparent thermal nature of multiparticle production even in the absence of significant final state scattering.


Introduction
Entanglement is one of the most puzzling features of quantum theory. It is the subject of intense debates about the foundations of quantum mechanics at least since the influential work of Einstein, Podolsky and Rosen [1]. Nowadays entanglement is seen as a key resource for quantum computation and quantum communication devices [2]. Entanglement in quantum field theory was investigated initially mainly with a view on black hole physics [3][4][5] and more recently in the context of holography [6].
For globally pure states, a good measure of the entanglement between a region A and its complement region B is the entanglement entropy. If one considers the reduced density matrix for A that follows from tracing over the Hilbert space associated with region B, ρ A = Tr B ρ , (1.1) this reduced density matrix is of mixed state form as a result of the entanglement between A and B. This can be quantified by the entanglement entropy defined as the von Neumann entropy associated with ρ A , (1.2) (More general Rényi entanglement entropies will be discussed later in the text.) A detailed understanding of entanglement entropy exists for 1+1 dimensional conformal field theory [7][8][9][10][11]. Technically, a replica trick in the Euclidean formulation of the theory leads to a partition function on an n-sheeted Riemann surface which can be evaluated. One can, as a result, compute the entanglement entropy not only of vacuum states but also of finite temperature states [9,10]. The formalism has even been extended to discuss nonequilibrium dynamics [12][13][14]. Further, a replica method has also been developed to compute the relative entanglement entropy between two density matrices [15][16][17]. Novel methods have also been developed for free field theories, in particular for massive (or massless) scalars and fermions [18][19][20]. These methods are based on the Euclidean formulation of these theories and also employ the replica trick. Analytical insights are supplemented by detailed numerical investigations using lattice techniques. In Ref. [20], alternative real time methods are mentioned but not fully developed.
In this work, we will be interested in the real-time dynamics of entanglement resulting from the rapid expansion of a system. We will concentrate on Gaussian states and on bosonic theories, a setup that will allow us to discuss the dynamics of an expanding string. The latter, as we shall soon discuss, will be a significant focus of this work. To follow the dynamics in real time, we will use a Schrödinger functional formulation for the corresponding density matrix. In order to determine the entanglement entropy of general Gaussian states (equilibrium or nonequilibrium, pure or mixed), it will be convenient to start with a somewhat abstract but fully general discussion of the corresponding mathematics. This discussion follows an approach that goes back to Refs. [3,4] and develops it further. (See also [21][22][23][24] for more recent work in this direction.) We will arrive at results that express the von Neumann and Rényi entanglement entropies (as well as relative entanglement entropies) directly in terms of traces over combinations of two-point correlation functions within the region A.
Quantum field theory when applied to the temporal evolution of systems, as for example in early universe cosmology or in the presence of time dependent background fields, reveals surprising features [25,26]. For example, excitations can be described with different mode functions; these (and the accompanying creation and annihilation operators corresponding to inequivalent vacuum states) are related by nontrivial Bogoliubov transformations. In a static Minkowski space problem, one basis is preferred by having mode functions with positive frequencies; this is not the case in time dependent situations with fewer symmetries. It is possible that an incoming vacuum state has nonvanishing particle number with respect to mode functions that are distinguished in having positive frequency states at asymptotically late times. In this case, the time dependence of the background field or of expanding geometry lead effectively to multiparticle production [25][26][27][28].
Entanglement plays an important role in a deeper understanding of such intriguing phenomena. Effective particle production happens typically in terms of entangled Einstein-Podolsky-Rosen pairs. When they separate in space, they contribute to entanglement between spatial regions. This is particularly relevant in the presence of horizons, for example close to a black hole or in the closely related Rindler wedge of spacetime, because causality dictates that observers there cannot recover the full information about a quantum state. As a result, Hawking radiation [29] or the closely related Unruh effect [30] lead to a thermal spectrum of particles. While Hawking and Unruh radiation concern idealized static situations governed by an event horizon, similar horizon phenomena can occur in time dependent situations. For example, observations at a given space-time point are, by reasons of causality, only sensitive to the interior of their past light cone. If the quantum field theoretic state is specified on some Cauchy surface in the past of an observer, only regions on this Cauchy surface inside the light cone are of relevance to her. This intersection of the light cone with the Cauchy surface constitutes a kind of particle horizon similar to the cosmological light horizon. The interior of this particle horizon typically fills a large volume; in this case, the entanglement entropy, if it scales with the area of its boundary, has a negligible effect. The situation can be very different in an expanding situation and can therefore lead to interesting unanticipated consequences. Further, this phenomenon can help to explain certain puzzling experimental observations in high energy collision experiments that have thus far defied explanation.
A long standing puzzle in elementary electron-positron collisions is that experimental results for particle multiplicities are well described by a thermal model corresponding to Boltzmann weighted distributions with a certain temperature T [31][32][33][34][35]. This is surprising because the theoretical picture we have about these collisions makes thermalization by multiple collisions unlikely. A popular theoretical model for soft QCD processes has been developed in Lund [36,37] and underlies, for example, the PYTHIA event generator [38,39]. It is based on expanding QCD strings from which hadrons and resonances are produced by tunneling processes via the Schwinger mechanism. In the standard implementation of the Lund model, as noted in ref. [40], the thermal-like features seen in experimental data are hard to understand. One attempt to cure this problem by allowing for fluctuations of the string tension, is discussed in [41].
In a quantum field theoretic description, different regions in a QCD string are in fact entangled. If one considers a sub-region, for example the region A in fig. 1, this interval can be described by a reduced density matrix corresponding to a trace over the complement region B as shown in eq. (1.1). Entanglement between different regions ensures that the reduced density matrix ρ A is of mixed state form. One may now ask oneself whether this reduced density matrix could resemble locally a thermal state and whether this could explain the close-to-equilibrium distribution of hadron ratios as found experimentally. We will explore this question within the Schwinger model of 1+1-dimensional QED, which is the framework underlying the above mentioned phenomenological Lund model. Remarkably, we find that at very early times, the density matrix for excitations around a coherent field describing the string in a rapidity interval is indeed of thermal form, with a temperature that decays with proper time, T = /(2πτ ) [42]. This effect is governed by an interplay of entanglement and the expansion dynamics. Our results in this direction are summarized in a recent letter [42].
We note that in the context of high energy hadron and heavy ion collisions, entanglement and entropy related ideas were discussed some time ago [43][44][45][46][47] as well as more recently [48][49][50][51][52]. Besides high energy collider physics, the dynamics of entanglement for quantum field theoretic models is also of relevance in other contexts. One prominent example is of course cosmology [25,26] and another concerns experiments with ultracold atomic quantum gases where entanglement can actually be investigated directly [53][54][55][56].
This paper is organized as follows. In section 2, we discuss Gaussian density matrices in quantum field theory and their entanglement properties. We start in section 2.1 with a description of pure Gaussian states described at fixed time or on some appropriate Cauchy surface in the Schrödinger functional representation. We use a notation with abstract indices which has the advantage that it can be applied to many concrete physics situations. In section 2.2, we discuss general mixed Gaussian states in terms of a variant of the Glauber-P representation [57,58] of the density matrix. A key formula here is an expression for the most general Gaussian density matrix in eq. (2.14). In subsection 2.3, we discuss the density matrix that emerges when traces over parts of the configuration space are performed. A key result here is that the reduced density matrix remains of Gaussian form with new covariance matrices that can be related to the original covariance matrix by certain projection operators. It is clear that Gaussian density matrices can be fully characterized by expectation values and connected correlation functions of fields and their canonically conjugate momenta and we discuss corresponding relations in subsection 2.4. Subsequently, in subsection 2.5, we derive expressions for the Rényi and von Neumann entropies of general Gaussian density matrices using the replica method at an intermediate step. The resulting expressions are directly governed by connected two-point correlation functions as displayed in eq. (2.33) and (2.34). Further simplifications of these expressions can be obtained by the use of canonical transformations, in particular transformations that employ the symplectic structure of fields and their canonical conjugate momenta. We discuss the corresponding mathematics in subsection 2.6 and emphasize the powerful insight provided by Williamson's theorem into the resulting form of the density matrix. A simple representation of the von Neumann entanglement entropy emerges in eq. (2.59).
In section 3 we illustrate the general formalism of section 2 by applying it to the derivation of the entanglement entropy of an interval of length L in 1 + 1 dimensional Minkowski space governed by a free scalar field theory. We consider the corresponding eigenvalue problem in a discrete Fourier representation and discuss the nontrivial role of boundary terms in subsection 3.1. This is used in subsection 3.2 for the computation of field-field and conjugate momentum correlation functions, and in subsection 3.3 for the computation of the entaglement entropy.
In section 4, we compute the entanglement entropy and observe that it generates a locally thermal state for relativistic particles in the rapidly expanding environment of a high energy electron-positron collision. In subsection 4.1, we investigate the Schwinger model and in subsections 4.2 and 4.3 the background field and the dynamics of perturbations of the expanding string solution that capture essential aspects of the collision dynamics. The entanglement entropy of an expanding string is computed in subsection 4.4. We describe the excitations around an expanding coherent field in terms of local density matrices that are of thermal form with a time dependent temperature T = /(2πτ ) in subsection 4.5 and provide a simple intuitive picture of how the apparent thermal character of these excitations arise. Our conclusions are stated in section 5.
Appendix A generalizes our formal discussion of entanglement entropies of Gaussian density matrices to relative entanglement entropies and appendix B provides a brief review of the symmetries, quantum anomalies and bosonization of the Schwinger model for quantum electrodynamics in 1+1 dimensions.
2 Entropies and entanglement of Gaussian states

Gaussian pure states
We shall consider bosonic quantum field operators Φ m in the Schrödinger picture and employ a set of eigenstates with eigenvalues φ m . For a compact notation, the index m can label both continuous degrees of freedom such as position or momentum as well as discrete ones such as spin or internal quantum numbers. In particular, for complex fields, we will also take their complex conjugate to be part of the "Nambu field" φ m . 1 Nevertheless, it will often be convenient to also talk about the complex conjugate field φ * m although it is not independent of φ m . If the quantum system is in a pure state |ψ , then φ|ψ represents its Schrödinger wave functional (to be understood as a functional of the field φ). We will be particularly interested in expanding systems, where at early times the expansion rate can be much larger than typical interaction or scattering rates of field excitations. In this case, the quantum dynamics is typically well described in terms of Gaussian wave functionals, as will be discussed in more detail in section 4.1. The most general form of a Gaussian Schrödinger wave functional representing a pure state for a complex bosonic field theory can be written as 2 In the second equation, we have used a condensed notation with φ denotes the collection of fields, and h, λ and κ are complex quantities that parametrize the state. Their physical meaning will be discussed below. Gaussianity refers to the maximum power of the field appearing in the exponential (2.2) being quadratic. Correspondingly, we have The norm is given by a functional integral over φ, where we use Dφ = n Dφ n with Dφ n = dReφ n dImφ n /π for complex fields and Dφ n = dφ n / √ π for real fields. Canonical normalization corresponds to this expression being unity. More generally, the scalar product between a state functional ψ[φ] specified by h, λ and κ and another oneψ[φ] specified byh,λ and κ is given by the functional integral Using two sets of eigenstates of (2.1), whose eigenvalues we denote by φ + and φ − , the pure-state density matrix can be written as (2.6) The condition for a pure state density matrix Tr{ρ 2 λ,κ }/Tr{ρ λ,κ } 2 = 1 is satisfied as it should be. We will also work with the canonically conjugate momentum field π m . In the Schrödinger representation employed, it is given by a functional derivative operator The representation (2.7) implies the canonical commutation relations and correspondingly for the complex conjugate fields.

Gaussian density matrices
The Gaussian density matrix we discussed was for a pure state. In the following, we will make the transition to mixed states in terms of the density matrix, but restrict ourselves to situations where the density matrix remains of the Gaussian form 3 . We can write such a mixed state density matrix in the form, where ρ λ,κ [φ + , φ − ] is the pure state density matrix in (2.6) (dependent on the fields λ and κ) and P [λ, κ] is a quasi-probability distribution. When positive, P [λ, κ] can be seen as the probability distribution for statistical noise in the parameter fields λ and κ. More generally, however, P [λ, κ] need not be positive semi-definite. (The density operator ρ should of course be hermitian and positive semi-definite.) Note that eq. (2.9) is closely related (although not identical) to the Glauber-P representation of a density matrix [57,58]. In the following, for simplicity, we will take P [λ, κ] to also be of the Gaussian form, with a hermitian operator Σ that we take to be of the form 4 For this to be a properly normalizable probability distribution, the eigenvalues of Σ should be positive. One may also introduce the linearly transformed parameter fields, 12) in terms of which the exponent in (2.10) becomes diagonal, Substituting (2.10) in (2.9), one obtains straightforwardly that the mixed state density matrix ρ[φ + , φ − ] is also Gaussian. In the limit where Σ → 0, the functional P [λ] approaches a delta distribution functional and one recovers the pure state we started with. Performing the functional integral over λ and an overall shift of fields (2.14) Interestingly, this is the most general Gaussian density matrix that satisfies the hermiticity property As a consequence, any Gaussian density matrix can be written in the form (2.9) for suitable Gaussian P [λ, κ]. The trace of (2.14) is given by 16) and in the following we will assume canonical normalization, tr{ρ} = 1, or divide by the appropriate power of the above expression.

Projections and reduced density matrix
For the computation of the entanglement entropy and related quantities, we shall require reduced density matrices that result from performing partial traces over some of the degrees of freedom. Towards this end, one may introduce a projection operator P = P † wherebỹ More general Gaussian forms for P [λ, κ] are possible but will not be needed for our construction.
are the fields we want to trace out, andφ ± m = P mn φ ± n (2.18) are those we wish to retain. For example, P might be a projector in position space which equals unity in some interval and is zero in the complement region. The reduced density matrix is formally given by The projectors have been inserted here for clarity.) For the Gaussian density matrix (2.14), one can formally perform the functional integral overφ and obtain yet again a Gaussian form for the reduced density matrix, (2.20) We have used here the abbreviations (2.21) The operator inversions used here are to be done within the subspace that is integrated out. Note that the reduced density matrix (2.20) is of the same general structure as the original density matrix (2.14); in particular, it still satisfies the hermiticity property (2.15). However if one starts with the density matrix of a pure state where Σ a = Σ b = 0, the reduced density matrix (2.20) contains in general terms that mix φ + and φ − (the terms d (+−) and d (−+) in (2.21)) and describes therefore a mixed state, as expected. Because the reduced density matrix is again Gaussian, it is completely determined by the field expectation values and the two-field correlation functions. As we will discuss, these can be computed directly on the domain of interest without further reference to projection operators.

Correlation functions
A Gaussian density matrix can be completely characterized in terms of field expectation values and correlation functions of two fields. More specifically, for (2.9) one has for the field φ m and the canonical conjugate momentum field π m = −iδ/δφ m , For the connected correlation functions AB c = AB − A B one finds (2.23) These are compatible with the canonical commutation relations, (2.24) Note that the matrices h, h † , Σ a and Σ b are fixed in terms of the connected correlation functions of φ, φ * and the momenta π, π * . For later convenience, we note the relations (2.25) They follow with some algebra from (2.23).

Entropy and entanglement entropy
In this subsection we will determine the Rényi entropy for the general Gaussian density matrix (2.14) and extract from it the von Neumann entropy. The Rényi entropy is defined by (we assume standard normalization tr{ρ} = 1), 26) and the von Neumann entropy follows from this [7, 10] as a limit For the purposes of the calculation that follows, it is clear 5 that one can dropφ and j, as they do not enter tr{ρ N }. Using otherwise the general expression for ρ in (2.14) leads after a straight-forward exercise in Gaussian integration to [19] Tr{ρ which contains the N dimensional cyclic matrix (with operator valued entries) We have used here the abbreviations introduced in (2.25). We define further 1 N to be the N -dimensional unit matrix and Z N is the N -dimensional cyclic matrix (Z N ) mn = δ (m+1)n . Here m, n are in the range 1, . . . , N and the index m = N + 1 is to be identified with the index m = 1. One can write Combining terms leads to a compact expression for the Rényi entropy, We have rederived here in completely general terms, and in the the functional integral representation, a result previously obtained only for a special case in the operator formalism [60]. From the result above, one can directly obtain an expression for the von Neumann entropy, Note that this is positive semi-definite because a − b and a + b are positive semi-definite. Note also, that a and b can be expressed in terms of correlation functions of fields and canonical momenta using (2.25). The above expression simplifies for b = 0 to S =Tr and it vanishes as expected for b = a = 0.
As a first example and check of this formalism, consider a free real massive scalar field in n-dimensional infinite Minkowski space. The correlation functions in this case can be written as where n( p) are occupation numbers. The symmetric mixed correlation function vanishes: φπ + πφ /2 = 0 and therefore b = 0. Here one can evaluate the (full) entropy easily, because a = (h + h † ) −1 Σ a as given in which is the standard result for free bosonic fields. As expected, the entropy vanishes for n( p) → 0. We note that similar considerations also hold for relative entropies involving more than one density matrix. This concept is discussed further later and the relative entropy for free scalar fields is discussed in more detail in appendix A.

Symplectic transformations, Williamson's theorem and entanglement entropy
The above expressions can be further simplified with the help of canonical transformations. One considers unitary changes of the field basis, They can be used to diagonalize hermitian operators such as h + h † as is the case in going from position to momentum space. These transformations have unitary representations as transformations of the Schrödinger functionals. This is clear as the scalar product (2.5) remains unchanged by unitary transformations of the field basis due to Dφ = D(U φ).
In addition to this, there is a larger class of transformations, which transform fields and momenta into each other. Consider the combined field Their canonical commutation relation defines a symplectic metric, where, symbolically, The transformations Indeed one can confirm that this relation defines a Lie algebra.
Recall that φ and π contain also the corresponding complex conjugate fields so that there are relations One may assume without loss of generality that there is a field basis where all fields are real such that there R mn = δ mn . Of course, the matrix R changes under the unitary, block diagonal transformations (2.38). More specifically, one has The matrix R transforms also by the symplectic transformations (2.42), Notice that in the field basis where φ m is real, there is no change in R for real symplectic transformations. Hence S * mp = S mp , as expected. In a field basis that differs from this by a unitary transformation, the symplectic transformation has different form (and is not necessarily real). Specifically, the symplectic transformation matrix in (2.42) and the generator J A transform under the unitary block diagonal transformations (2.38) as We need to show that the symplectic transformations (2.42) have unitary representations as transformations of the states of the field theory. This is best done in the field basis where φ m are real fields and R = 1. The symplectic transformations are then real and one has There is now a representation of the Lie algebra specified by (2.44) in terms of the operators ( [61], see also [62]) acting on a Schrödinger functional. Indeed, one can confirm that they have the same commutation relations as the generators J A . Moreover, one has (X A ) † = X A in the sense of the bilinear form (2.5) so that the symplectic transformations indeed have unitary representations. This is an important result because it allows one to use the symplectic transformations to simplify calculations, for example of the entropy. Because symplectic transformations have unitary representations, they do not change the entropy by construction. Finally, we note that (2.50) is invariant under the block diagonal unitary transformations (2.38) and can therefore be used in any field basis. It is also clear that the corresponding unitary transformation maps Gaussian states to Gaussian states.
In particular, the hermitian and positive covariance matrix corresponding to the symmetrized correlation function which is a key ingredient in our discussion of Gaussian states, transforms under symplectic transformations as This is not a similarity transformation because S † = S −1 . In other words, the eigenvalues of ∆ are not invariant under symplectic transformations. Williamson's theorem states (see [63] for a discussion) however that there must exist a symplectic transformation that brings ∆ to diagonal form, with real and positive λ j > 0. These latter are the symplectic eigenvalues of the symmetrized covariance matrix. This is realized by considering the combination ∆Ω. One can show that it transforms as which indeed satisfies the properties of a similarity transformation. The eigenvalues of this combination, ±λ j , are directly related to the symplectic eigenvalues. It is therefore convenient to determine the eigenvalues of ∆Ω and to thereby relate observables such as the entanglement entropy to the Williamson form. We first note that the Williamson form expression (2.53) of the symmetric correlation matrix (2.51) results in a very simple form for the quantities in (2.25): Since the Heisenberg uncertainty relation tells us that a ij has to be positive-definite, this indicates that λ j ≥ 1/2. The symplectic transformations we have discussed and the resulting use of Williamson's theorem leads to a very convenient form for entropies. The entropy in (2.35) can be directly expressed in terms of the symplectic eigenvalues as Moreover, the symplectic eigenvalues λ m follow as pairs of conventional eigenvalues ±λ j of the combination ∆Ω. Following Sorkin ([23], see also [24]), a further simplification can be obtained by considering the matrix which has the eigenvalues ω + j = 1/2 + λ j and ω − j = 1/2 − λ j . As noted, the uncertainty relation gives us λ j ≥ 1/2; therefore, ω + j ≥ 1 and ω − j ≤ 0. One can then write (2.56) simply as where the sum is over pairs of eigenvalues. More simply, In the last expression, the sum is now over all the eigenvalues of D. Each negative eigenvalue ω m < 0 is paired with a positive one 1 − ω m . A pure state without entropy has ω m ∈ {0, 1}. Finally, we note that the matrix D in symbolic form can be expressed as, Hence, the entropy associated with a Gaussian density matrix is fully determined from the set of connected correlation functions of (2.60) evaluated in the domain of interest. It is understood that the operator trace in (2.59) is also restricted to this domain. Thus far, we have concentrated on the Rényi and von Neumann entropies of a single Gaussian density matrix. It is also possible to determine relative entropies between two Gaussian density matrices ρ and σ in a similar way and we discuss this in appendix A. (See ref. [64] for a general exposition on the concept of relative entropy.) Let us remark here that Williamson's theorem is not as useful for the determination of the relative entropy of two Gaussian density matrices as it is for the determination of the entropy of a single one. This is because it is not guaranteed that there is a basis in which the covariance matrices ∆ (ρ) and ∆ (σ) (defined in (2.51)) simultaneously assume their Williamson diagonal form. However, this should be the case when the matrices ∆ (ρ) Ω and ∆ (σ) Ω (which transform under symplectic transformations as similarity transforms, see (2.54)) commute, i.e. [∆ (ρ) Ω, ∆ (σ) Ω] = 0. One can then write the relative entropy as where the sum goes over all pairs of simultaneous eigenvalues (ω

Entanglement entropy in Minkowski space
In this section, we will illustrate our general formalism by applying it to the derivation of the entanglement entropy of an interval of length L in 1+1 dimensional Minkowski space governed by a free scalar field theory. This is a well studied problem with alternative approaches and results are available from both numerical and analytical methods [20]. The conventional numerical method to deal with this problem is to discretize the entire theory on a spatial lattice. The interval L corresponds to a finite subset of lattice sites and they are entangled with the lattice sites in the complement region. Numerically, one can study the continuum theory in the limit of finer and finer lattices. The entanglement entropy itself is ultraviolet (UV) divergent and a regulator is provided by the lattice spacing. However universal quantities such as derivatives of the entanglement entropy with respect to the interval length L can be extracted as well. Alternatively, one may consider relative entropies that are finite in the continuum limit as discussed in appendix A. The approach we developed in section 2 has the advantage that it depends only on the correlation functions inside the interval considered. Information about entanglement with degrees of freedom outside this interval is taken into account by appropriate boundary conditions in a nontrivial way. While it is rather straightforward to treat a field theory on a finite interval with periodic boundary conditions, which corresponds in fact to an isolated system, it is more involved to properly treat the theory on an interval that is not fully isolated but entangled with a complement region. To illustrate the subtle differences properly is a major focus for the following discussion.

Eigenvalue problem and boundary conditions
We will consider an interval (−L/2, L/2) in Minkowski space with one spatial dimension. A free scalar field φ will be governed by a Gaussian reduced density matrix on this interval. Moreover, the corresponding matrix entries, namely the functions h, h † , Σ a and Σ b introduced in section 2, will be such that the correlation functions have the same form as in infinite space; they are just restricted to the interval. The technical difficulties arise now from the fact that products of these functions involve integrals over the interval (−L/2, L/2), only. For example, the quantity Note that this is a nondiagonal matrix in momentum space for finite L. Only in the limit L → ∞ does one obtain sin( 1 and a becomes diagonal in momentum space (and zero). The challenge is now to find the eigenvalues of the matrix a(x, y) on the interval (−L/2, L/2).
To solve the eigenvalue problem, we will use a discrete basis involving Fourier expansion on the interval (−L/2, L/2). In doing so, we will not assume periodic boundary conditions. We first divide the relevant function (or field) into a symmetric and an anti-symmetric part, The symmetric part can be expanded into a Fourier series In a similar fashion, one can expand the anti-symmetric part We can summarize this as For ϕ(x) ∈ R, one has ϕ n = ϕ * −n . Note that this type of Fourier expansion does not assume periodic boundary conditions for ϕ(x).
In the limit of large interval length L → ∞, eq. (3.5) becomes formally (with p = nπ/L andφ(q) = ϕ n ), It is useful to relate these expressions to the standard momentum space representation defined as usual by One has for finite interval length L ϕ n = dp 2π sin( pL and for very large L formallyφ One recovers the standard Fourier transformed field 2φ(p) →φ(p) only if the second, strongly alternating, term above is neglected.

Field-field and conjugate momentum correlation functions
We shall now determine the matrix expressions for the correlation functions ∆ S φφ (x − y) = φ(x)φ(y) and ∆ S ππ (x − y) = π(x)π(y) . We first note that the correlation functions in position space can be decomposed as and similar for the canonical momentum field correlator. The cross terms like ∆ φφ (x, y) is anti-symmetric with respect to parity.
For the combined operator in (3.1), it is natural to decompose a(x, y) = a (ss) (x, y) + a (aa) (x, y) with (3.12) The products of symmetric and anti-symmetric operators vanish under parity symmetry. With respect to the discrete indices m and n, one infers that a mn is only nonzero when both indices are even or when both are odd, but that there can be no cross-terms.
For the momentum space representation (2.36), using (3.8) (for n( p) = 0) one finds that (3.13) For a state with nonvanishing occupation number, one has to insert factors (1 + 2n( p)) on the right hand side of these expressions. For the parity even-even (ss) (m, n even) and odd-odd (aa) (m, n odd) components one has and likewise for ∆ ππ , where the square root appears in the numerator.
In computing these integrals, note first that there are no poles on the real p-axis; we can therefore pull the contour slightly below the real axis. In addition, we can write The integral involving the first bracket on the r. h. s. can be closed above the real p-axis while the second can be closed below. The integral that is closed above picks up a contribution from poles on the real p-axis for m 2 = n 2 as well as a contribution from the branch cut of the square root. The integral that is closed below has contributions from only the branch cut. The contribution from the integral over the poles is simply δ m−n,0 ± δ m+n,0 4 ( mπ L ) 2 + M 2 .

(3.16)
This result is in fact the ground state correlator one would have obtained by quantization of the scalar field theory directly on the interval (−L/2, L/2) with periodic boundary conditions for parity-even modes and anti-periodic boundary conditions for parity-odd modes. It corresponds to a pure state without entanglement. Indeed, keeping only (3.16) together with the corresponding (poles only) approximation for the momentum-momentum correlation function, would lead to a vanishing entanglement entropy. We therefore see concretely from this example that the nontrivial contribution to the correlation function taking entanglement properly into account arises actually from the branch cut contribution to the integral in (3.14) and in its counterpart for ∆ ππ . The contribution from the branch cuts to the field-field correlator, as determined by the right hand side of (3.14), is given by the integral We shall first discuss the result without the exponential in the last bracket. This should be a good approximation in particular for (M L) 1. For the opposite limit (M L) 1, we will add the contribution from the exponential term as well. Performing the integral gives for m 2 = n 2 , Finally, for m = n = 0 one obtains (3.20) Adding up these expressions gives the field-field correlation function in the massive case for M L 1. In the massless case M = 0, one can also perform all integrals and we obtain the field-field correlators, For the zero mode with m = n = 0, we have introduced an infrared momentum regulator k. Analytic continuation from nonvanishing m = n suggests kL = e 1 2 −γ , where γ is Euler's constant; we will make this choice in the following. In the above formula, we used the modified cosine integral and sine integral functions given byc

(3.22)
Turning now to the canonical momentum field correlator ∆ ππ , following steps similar to the above, we find the contribution from poles to be (δ m−n,0 ± δ m+n,0 ) 1 4 ( mπ L ) 2 + M 2 . Note that this expression has a logarithmic UV divergence. We will follow the same strategy as for the field-field correlation function and shall discuss first the integral without the exponential term. This is a good approximation in the massive or long distance limit (M L) 1. (Its contribution will be added in the massless limit.) To deal with the UV divergence, we first consider the case m = n = 0 where one obtains with a UV momentum cutoff Λ and assuming Λ M , We subtract now this contribution from the integral and obtain for the remainder when m 2 = n 2 , (3.26) and for m 2 = n 2 , By adding these expressions together, one obtains the field-field correlator in the massive case (M L) 1. One can also perform the integrals for the canonical momentum correlation function in the massless limit to obtain,  Here again Λ is the UV momentum cutoff that we introduced to regularize the momentum integral. We assume ΛL m, n.

Entanglement entropy
After laying the groundwork above, one may now proceed to determine numerically the eigenvalues of the matrix in an approximation that consists in truncating the infinite sums to a finite range. This can be used to determine the entanglement entropy in terms of (2.35). Alternatively, one can determine the spectrum of eigenvalues of the matrix D in (2.60) and use (2.59). Typically N max = O(100) terms are enough. The result for the entanglement entropy in the massless case obtained thus is shown as a function of ΛL in Fig.  2 for N max = 300. We have verified that the result depends only very weakly on the precise value of N max as long as N max ΛL. Note that this is the condition under which the expressions for the correlators shown above are valid. Our result shows clearly that the entanglement entropy for a single massless bosonic field is well described by the anticipated result [7-10] with central charge c = 1, as expected.
The above formalism can also be used for a massive scalar field although the treatment is numerically more involved. The reason is that now the problem has three length scales given by the interval length L, the (inverse of) the UV cutoff Λ −1 and the inverse of mass M −1 . An additional UV momentum scale is set by N max /L. One must be careful to interpolate observables into the physical regime without leaving the range of validity of approximations made in due course. For numerical investigations, it may in this case be advantageous to use a lattice discretization as discussed for example in Ref. [20].

Entanglement entropy for expanding systems
In relativistic dynamics, the separation between a causally connected spacetime region and its acausal complement leads to a natural description of the finite observable region in terms of a reduced density operator. Moreover, the presence of an effective "lightcone" for the propagation of signals, the spread of entanglement, and the generation of correlations is also the subject of topical studies in nonrelativistic systems realized with ultracold quantum gases. In these latter cases, the dynamics is often well described in terms of generalized coordinates that reflect the presence of a lightcone. The expressions we will derive in this section will be useful in describing the generation of entanglement entropy in both these cases.
We will compute here the entanglement entropy and obtain its description in terms of a thermal density matrix in the bosonized massless Schwinger model in 1+1-dimensions. As noted previously, this is a popular model whose dynamics is observed to describe key features of particle production in elementary electron-positron collisions. After motivating the massless Schwinger model, we we will discuss its bosonized version in terms of massive scalar fields. This will allow us to employ the machinery we developed in the previous sections to compute the time evolution of the entanglement entropy. We will also provide the corresponding results for fermionic fields towards the end of this section.

The Schwinger model
The Schwinger model describes vector-like QED in 1 + 1 space-time dimensions. For a single massive fermion, the Lagrangian is The parameters of the model are the fermion mass m and the U(1) charge e. Note that the latter has dimensions of mass in two spacetime dimensions. The parameter e is related to the string tension, as one can see from the following consideration. The energy per unit length of the electric field between charges e and −e is given by E 2 /2. Moreover, from the Gauss law in 1+1 dimensions one has E 2 = e 2 such that the string tension is σ = e 2 /2. As is well known, the Schwinger model in two dimensions can be bosonized exactly [65]. (For completeness, we include a discussion of this feature of the model in Appendix B.) The model can be reformulated as a scalar field theory with Lagrangian (see eq. (B.20)) Here γ is the Euler constant and θ is the vacuum angle. Note that QED in two dimensions has topologically non-trivial vacuum structure because the gauge group U(1) and the "boundary of Euclidean space at infinity" have the same topology S 1 . As is clear from the explicit bosonization procedure in appendix B, the scalar bosons are quadratic in the original field and correspond to fermion-antifermion bound states φ ∼ψψ. Hence the fermions themselves are not part of the physical spectrum -the theory displays (geometric) confinement.
For the general case of a nonvanishing fermion mass m, (4.2) is still a nontrivial, interacting theory that cannot be solved easily. In particular, one expects a nontrivial renormalization of the effective potential as well as propagator terms. Also additional bound states that are of quadratic or higher order in φ could arise. The situation simplifies substantially in the strong interaction limit e 2 m 2 where one may set m = 0. The massless Schwinger model becomes a free scalar field theory after bosonization with a scalar boson of mass M = e/ √ π. Entanglement in the Schwinger model, as well as in the 't Hooft model, has been investigated for static situations in ref. [66].

General coordinates and background evolution
If one considers the free scalar field theory with mass M in arbitrary curved coordinates, one has the action For the string that forms between a highly energetic quark-antiquark pair produced in electron-positron collisions, it is natural to consider a boost invariant expansion formulated in Bjorken coordinates with proper time τ = (x 0 ) 2 − (x 1 ) 2 and rapidity η = arctanh(x 1 /x 0 ). The metric is g µν = diag(−1, τ 2 ) and √ g = τ . The non-vanishing Christoffel symbols are Γ τ ηη = τ , Γ η τ η = Γ η ητ = 1 τ . We will be interested in a situation where the external charges at the endpoints of a string generate a field expectation valueφ = φ , where the physics dictates that this background field is boost invariant in the rapidity variable. The free-field equation of motion in the Schwinger model for the "Bjorken expansion" of the η-independent field expectation valueφ(τ ) is then given by [67] The two independent solutions areφ(τ ) ∼ J 0 (M τ ) andφ(τ ) ∼ Y 0 (M τ ). While the former is regular for τ → 0, the latter diverges. Both oscillate for late times τ 1/M , which can be understood as multiple string breaking [68]. In the limit of vanishing mass M → 0, the corresponding two independent solutions to the equation of motion (4.4) areφ(τ ) ∼ const andφ(τ ) ∼ ln(τ ). To fulfill the Gauss law, the electric field E = eφ/ √ π must approach the U(1) charge of the external quarks, E → ±e for τ → 0 + . Therefore we obtain the solutionsφ (τ ) = ± √ πJ 0 (M τ ). (4.5) The energy-momentum tensor has the form T µν =¯ u µ u ν +p(u µ u ν + g µν ) with the "fluid velocity" u µ = (1, 0), energy density¯ = 1 2 (∂ τφ ) 2 + 1 2 M 2φ2 and "pressure"p = 1 2 (∂ τφ ) 2 − 1 2 M 2φ2 . One may check easily that the ideal fluid equation of motion, is satisfied. However here energy density and pressure are not related by a fixed equation of state as in local thermal equilibrium but they are both determined dynamically byφ. The massless case, M = 0 is an exception which satisfies the equation of statep =¯ for pure radiation in 1+1 dimensions. More generally, depending on the initial conditions forφ and ∂ τφ , one can have initial conditions betweenp = −¯ and p =¯ . For the solution in terms of J 0 (M τ ), one initially hasp = −¯ with the result oscillating between this value andp =¯ . Note that the Bjorken metric g µν = diag(−1, τ 2 ) does not have a Killing vector field pointing in τ direction; there is no solution ξ µ = (f (τ ), 0) of the equation However, there is a conformal Killing vector field of this form which is a solution of given by ξ µ = (c τ, 0) with λ = 2c with some constant c. These observations allow us to infer that no local thermal equilibrium state is conceivable which has Bjorken boost symmetry except for the case of a conformal field theory. This insight will play a role in our discussion later on the emergence of a local thermal state and the dynamics generating entanglement entropy.

Dynamics of perturbations
In the following, we consider perturbations or fluctuations around the background solutionφ. We begin by writing the field as φ(τ, η) =φ(τ ) + ϕ(τ, η). (4.9) The fluctuation field has the equation of motion We will now discuss the quantization of the field ϕ. Time dependent quantization problems of this type are best solved in terms of convenient mode functions. One writes the quantized field as where a(k) and a † (k) are annihilation and creation operators.
The mode functions depend only on the magnitude of the wave vector |k| and are solutions to the differential equation The inner product of these mode functions provides the normalization condition The canonical conjugate momentum field is with |α(k)| 2 − |β(k)| 2 = 1. The corresponding creation and annihilation operators are related bȳ The different sets of mode functions correspond to the vacuum states |Ω and |Ω respectively. These vacuum states contain no excitations in the sense that they are annihilated by the operators a(k)|Ω = 0 in the formed case and byā(k)|Ω = 0, in the latter case. Note however that the vacuum |Ω might contain excitations with respect to the operatorā(k) and the vacuum |Ω with respect to a(k). However, the Bogoliubov transformations that connect the different choices are unitary transformations and therefore do not change entropy. Typically, the vacuum with respect to one set of mode functions corresponds to a squeezed state with respect to other sets of mode functions. In particular, (4.12) has the two independent solutions J ik (M τ ) and Y ik (M τ ) or, equivalently, H (2) ik (M τ ) and its complex conjugate H (1) −ik (M τ ). The normalized mode functions corresponding to the latter choice are The set of mode functions f (τ, k) in (4.17) is distinguished by being a superposition of positive frequency modes with respect to time t of standard Minkowski spacetime [25]. This means that the standard Minkowski vacuum will be also a vacuum with respect to these mode functions in Bjorken coordinates. In the limit of vanishing mass M τ → 0, the mode function (4.17) becomes One observes that it contains both positive and negative frequency contributions with respect to the logarithm ln(τ ) of Bjorken time.
An alternative choice of mode functions is In the limit of vanishing mass M τ → 0, the mode function (4.19) becomes In this case, we see that it has only positive frequency contributions with respect to ln(τ ). The phase in the last equation is given by θ(k, M ) = k ln(M/2) + arg(Γ(1 − ik)). Note, in particular, that the factor multiplying k diverges in the formal limit M → 0. The Bogoliubov coefficients that connect the mode functions (4.17) and (4.19) are .
The Gaussian density matrix or the vacuum state for this problem can be specified in terms of field expectation values and the connected correlation functions, For example, the vacuum state with n(k) = u(k) = u * (k) = 0 results in the correlation function It is instructive to compare (4.24) with the corresponding Minkowski space expression Employing the integral representation where x 0 = τ cosh(η) and x 1 = τ sinh(η) are standard Minkowski coordinates, one obtains This agrees with (4.25) after substituting k 0 = κ cosh(z), k 1 = −κ sinh(z), thereby confirming that the mode functions in (4.17) are indeed the ones corresponding to the standard Minkowski space vacuum. Identifying π(τ, k) = τ ∂ τ φ * (τ, k), a complete set of connected correlation functions for the vacuum with n(k) = u(k) = u * (k) = 0 in momentum space is given by where we have used the abbreviationḟ (τ, k) = τ ∂ τ f (τ, k). One can directly verify that with the above correlators, the matrix D in (2.60) has pairs of eigenvalues {0, 1} such that the entropy associated with the entire expanding string is zero. In fact, one does not need the precise form of f (τ, k) to show this; the normalization condition in (4.13) alone is sufficient. Now with respect to the alternative set of mode functions, with only positive frequency solutions, the state with n(k) = u(k) = u * (k) = 0 has the set of correlation functions, In this alternative basis, the correlation functions do not look like those of an empty state but rather of one with occupation numbern(k) = |β(k)| 2 . From (4.22) one obtains, Recalling that the single particle energy in the expanding situation is E = k 2 /τ 2 + M 2 , for massless bosons, this distribution appearing in the "diagonal" elements of the correlation matrix (4.29) corresponds to a thermal spectrum with the time dependent temperature Such a thermal interpretation is not possible for a nonvanishing mass M , but the fact that the quasiparticles defined by the mode functionsf (τ, k) have a nonvanishing occupation number remains true. Of course, a thermal state would have vanishing entries for the other correlators appearing in (4.29). The "off-diagonal occupation function" isū . (4.32) One may use the relations in (4.29) to express the correlation functions in (4.28) in the alternative basis,  Note that the "off-diagonal occupation functions"ū(k) that appear here are always multiplied with sine or cosine functions that are strongly oscillating with k in the limit M τ → 0. An interpretation of this term in position (rapidity) space is obtained by Fourier transforming the correlators above back to position (rapidity) space and examining the structure, for instance, of φ(τ, η)φ(τ, η ) . The strongly oscillating terms ∼ cos[2k ln(M τ /2)] correspond then to pronounced structures in the rapidity difference η − η of the spatial correlators at separation |η − η | ∼ 2| ln(M τ /2)|. This rapidity separation becomes large in the limit M τ → 0.

Entanglement entropy of an expanding string
We shall now investigate the entanglement entropy of an interval with length ∆η in rapidity at some fixed Bjorken time τ . Following the discussion in section 3 for the static Minkowski space case, one has to use a discrete Fourier expansion scheme for the finite interval (−∆η/2, ∆η/2). As discussed there, to ensure boundary conditions are not restricted, it is most convenient to split the fields into symmetric and antisymmetric components. The calculation proceeds by determining correlation functions as in eq. fixed Bjorken time τ and in the finite interval. As previously, one can use eq. (2.59) to determine the entanglement entropy. The correlation functions on the finite rapidity interval can be most conveniently obtained from the momentum space representations (4.28) or (4.33), together with the analog of the relation (3.8) expressing the discrete field basis in terms of standard Fourier modes with an appropriate integral kernel. In our case here, this becomes (at fixed time τ ) (4.35) When one calculates correlators in the discrete basis, such as φ n φ −m c , from (4.34) using the kernel (4.35), one observes that the terms proportional to the off-diagonal occupation numberū(k) in (4.34) contain a term that oscillates very fast with k in the limit M τ → 0. This is because the combination 2k ln(τ ) + 2θ(k, M ), with the phase θ(k, M ) in (4.21), is strongly dependent on k with infinite derivatives contributing in the limit M τ → 0. Hence these terms effectively do not contribute to the correlators such as φ n φ −m for a finite length interval ∆η in the limit M τ → 0.
One is then left with correlators that describe vacuum fluctuations and the occupation numbersn(k) corresponding to a thermal distribution at very early times τ . This shows that the two limits M τ → 0 and ∆η → ∞ do not commute. If one considers M τ → 0 for finite ∆η, one finds the entanglement entropy of a thermal state with the time dependent temperature given by (4.31). This holds even if one then takes the limit ∆η → ∞. In contrast, if one considers an infinite interval ∆η → ∞ for finite M τ , one finds a pure state with vanishing entanglement entropy. This holds also if one considers subsequently the limit of vanishing mass M τ → 0.
Remarkably, one finds that at the very early times M τ → 0, this entanglement entropy in a finite rapidity interval ∆η is equivalent to that of a 1+1 dimensional conformal field theory at finite temperature when T = 1/(2πτ ). In such a 1+1 dimensional conformal field theory at temperature T , the entanglement  Alternatively, one can obtain the result in (4.37) for the expanding system from the Minkowski space computations in our framework, as discussed in section 3. In obtaining this result, we first make use of the fact that the coherent field does not contribute to the connected correlation functions in the covariance matrix and can be dropped from the computation of the entanglement entropy. Furthermore, the entanglement entropy of an interval is unchanged by unitary evolution as long as the boundaries are kept fixed. This allows one to reduce the problem to finding the entanglement entropy of an interval in Minkowski space at a constant time with length ∆z = L = 2τ sinh (∆η/2) at a fixed time t = τ cosh(∆η/2); this corresponds to the dotted red line in Fig. 3. Using (3.30), this leads to (4.37) as well. This reasoning has been used for the derivation of the results first obtained in [42].
The additive constant in eq. (4.37) is not universal but the derivatives of S with respect to τ and ∆η are universal. This implies τ ∂S/∂τ = c/3 and ∂S/∂∆η = (c/6) coth(∆η/2). For a large rapidity interval ∆η 1, one has S = (c/6)[∆η + 2 ln(τ )] + const. This shows the existence of a time-independent piece of the entanglement entropy that is extensive in rapidity and a ∆η-independent piece that grows logarithmically with the proper time.
At later times, for the nonconformal case of free massive scalars, the universal part of the entanglement entropy behaves as in the conformal case for M ∆z 1 and decays for M ∆z 1 [20]. In this case, the derivative with respect to ∆η of the entanglement entropy gives  where c E (M ∆z) = ∆z∂S/∂∆z is the entanglement entropy c-function for a massive scalar field. Taking the conformal limit, one obtains, as anticipated, that c E (0) = c/3 = 1/3. For large values of the argument, this function has the form c E (x) → xK 1 (2x)/4, which decays exponentially.
We can employ the general expression for c E (x) [20] for massive free bosons to compute dS/d∆η for the massless Schwinger model. The result is displayed in Fig. 4 for the different values of M τ that are shown in the caption (dashed lines). For short times τ M 1, one observes a significant entanglement over rapidity intervals ∆η = O(1). At intermediate values of τ M and ∆η, one observes that ∂S/∂∆η approaches a plateau at 1/6 as a function of ∆η at early times. One also sees that it decays both for very large ∆η and for later times τ . The plateau is governed by the conformal limit M τ → 0 which is shown by the solid black line in the figure.
Since the entanglement entropy c-function in the conformal limit is identical for real massless scalar bosons and for massless Dirac fermions, this indicates that the conformal limit is consistently described in the Schwinger model with or without bosonization. One may also determine the entanglement entropy of free massive Dirac fermions using similar manipulations as described above and using the corresponding c-function given in Ref. [20]. The result is also shown in Fig. 4 (dotted lines). The approach to the universal plateau at dS/d∆η = 1/6 for M τ → 0 is even faster for the free fermion case.
A more intuitive (but also more heuristic) description of why the vacuum state |Ω in the conformal limit looks thermal in a finite rapidity interval ∆η is as follows. Although a pure state, the state describing the expanding string contains entangled pairs of quasiparticles with opposite momentum in the Bogoliubov basis (4.19) which consists of modes with positive frequencies with respect to ln(τ ) for M → 0. This particular basis is special because only there can one interpret quasiparticles (in the classical quasiparticle limit) as moving in space on well defined trajectories. For a rapidity interval ∆η, this implies that quasiparticles constantly come in via the left and right boundaries. They are entangled with other quasiparticles moving in the opposite direction but that is not seen locally. Because these quasiparticles have a thermal spectrum, local observables will effectively look thermal. A related argument was employed previously to understand the time evolution of entanglement entropy after a quantum quench [12,14].

Local density matrix of an expanding string
One can go even further and make stronger statements regarding the thermal character of entanglement entropy in a finite interval of the expanding string. Note that the correlators in (4.34), when projected to the finite interval with the kernel (4.35), are at M τ → 0 exactly those of the 1+1-dimensional conformal field theory in thermal equilibrium if the temperature is identified to be T = 1/(2πτ ). Because Gaussian density matrices are fully specified by one-point expectation values and two-point correlation functions, the density matrix obtained by first taking the limit M τ → 0 and then ∆η → ∞ has the thermal form, where K is the so-called modular or entanglement Hamiltonian defined as Equations (4.39) and (4.41) specify the density matrix of a locally thermal state where β µ = u µ /T is the so-called inverse temperature vector given by the fluid velocity u µ (pointing in τ -direction) and the temperature is T = 1/(2πτ ). Moreover, T µν is the energy-momentum tensor of excitations ϕ(τ, η) around the coherent fieldφ(τ ) according to (4.9). Note also that β µ = (1/T, 0) (in coordinates τ, η) is actually a conformal Killing vector according to (4.8). This result can also be obtained as a limit of a more general result, as discussed already in Ref. [42]. A conformal field theory in the vacuum state in a region with a boundary formed by the intersection of two light cones (see figure 3), has a reduced density matrix of the form (4.39) on any hypersurface in that region that has the same boundary. The modular Hamiltonian is a local expression given by [69,70] (see also [71]) Here T µν is again the energy-momentum tensor of excitations and ξ ν is a vector field that can be written as where x is the space-time position on the hypersuface we consider, q is the end-point of the future light cone and p the starting point of the past light cone that form the boundary of the region. Note that (4.41) is again of the same form as a density matrix of a thermal state if one identifies ξ µ = 1 T u µ the vector of (inverse of) temperature T and fluid velocity u µ . The vector ξ µ vanishes on the boundary of the region enclosed by the two light cones corresponding formally to an infinite temperature.
Consider now a situation where the two enclosing light cones intersect on a constant τ hypersurface in Bjorken coordinates with a rapidity difference ∆η. If we concentrate on the point in the middle of this rapidity interval, the vector ξ µ points in τ -direction and has length 2πτ (cosh(∆η/2) − 1)/ sinh(∆η/2). The associated temperature approaches precisely T = 1/(2πτ ) for ∆η → ∞.
The state ρ we are considering here is vacuum-like but may have a nonvanishing coherent field and therefore, a corresponding nonzero energy. Further, excitations that can be formed from local unitary (entropy preserving) transformations lead to a modified density matrix ρ 1 . Using the notion of a relative entropy to ρ as explained in appendix A, we have where P µ = Σ dΣ ν T µν is the four-momentum associated with the perturbation. This is precisely the characteristic of a thermal state. Within a large but finite rapidity interval, the relative entropy of the state with a small perturbation compared to the coherent field state of the expanding string has the same value as in a thermal state with temperature T = 1/(2πτ ). In the expanding geometry of interest, quantum fluctuations of this kind are therefore as likely as in a thermal state. Such fluctuations should therefore be observed with a distribution corresponding to that of a grand canonical ensemble.

Conclusions
We developed in this paper a powerful formalism exploiting Gaussian density matrices to examine different entanglement entropies that arise in a wide variety of equilibrium and nonequilibrium problems in quantum field theory. The most important results of this exercise are expressions for the entanglement entropy in terms of two-point correlation functions that can be directly evaluated within the finite spacetime interval of interest.
In particular, we revisited the computation of the entanglement entropy of an interval of length L in static, two-dimensional Minkowski space and discussed how it can be computed within our formalism. We find that it can be fully determined by correlation functions within the interval. Our results clearly reveal the importance of boundary conditions in the emergence of the entanglement entropy in finite spacetime intervals.
We then applied the insights gain from our general treatment to study entanglement in the expanding string formed between a highly energetic quark-antiquark pair after an electron-positron collision. In order to do so, we employed the Schwinger model of quantum electrodynamics, which is the basis of much of the phenomenology describing multiparticle production in these collisions. We exploit the fact that the Schwinger model can be bosonized and is (in the limit of vanishing fermion mass m) equivalent to a free bosonic theory. The expanding string is best described in Bjorken coordinates of proper time τ and rapidity η and corresponds then to a coherent field solution of the Klein-Gordon equation in an expanding geometry. The entanglement properties of this state are fixed in terms of correlation functions for excitations around the coherent field. We discuss them in terms of appropriate mode functions. We observe that different sets of mode functions are possible and that they are related by Bogoliubov transformations.
We find that the state corresponding to the standard Minkowski space vacuum appears at early time τ as an occupied squeezed state with mode functions that have positive frequency with respect to the logarithm of Bjorken time ln(τ ). Moreover, as we showed very explicitly in sections 4.3 -4.5 and had discussed previously in [42], this state appears in any finite rapidity interval as a thermal state governed by a τ -dependent temperature T = /(2πτ ). This is a rather remarkable finding, since it implies that excitations around the coherent field solution of the expanding sting are in fact thermal at very early times! Although our model for a QCD string is not fully realistic, we believe that the above described mechanism provides a compelling candidate for a deeper understanding of the approximately thermal distributions of hadron ratios found in many collider experiments even in contexts where scattering effects contributing to thermalization of particles in the final state are likely small.
A crucial question beyond the framework investigated here is what happens when interactions are taken into account. For the bosonized massless Schwinger model, they arise from a nonvanishing fermion mass m, leading to a Sine-Gordon type interaction term as in eq. (4.2). At very early Bjorken time τ , the interaction term ∼ me is expected to be irrelevant but could start to play a role at intermediate times and -for large enough m -even before the boson mass term ∼ e 2 /π. We expect that further correlations, and therefore entanglement, can be built up by interactions between excitations and we plan to investigate their consequences for the entanglement dynamics in future work.
While we have concentrated here on dynamics in 1+1 dimensions, similar entanglement dynamics may be at play in a Bjorken-type expanding geometry with additional transverse dimensions. If this is the case, quantum entanglement may be an important feature of early time dynamics in heavy-ion collisions. We plan to explore this exciting possibility in future work.
Similar considerations also apply to nonrelativistic condensed matter systems such as ultra-cold quantum gases whose "horizon" is set by the speed of sound. Entanglement in such many-body systems can be efficiently explored with table-top experiments enabling direct confrontation of theoretical predictions with experimental measurements.

A Relative entropy
Extending the discussion of section 2.5, we now consider two density matrices ρ and σ. They may be reduced density matrices originating from the trace over some part of the Hilbert space. We will be interested in the quantum relative entropy or Kullback-Leibler divergence of ρ with respect to σ, defined by S(ρ|σ) = Tr{ρ (ln ρ − ln σ)}. (A.1) The relative entropy measures in an information theoretic sense to which extend one can distinguish the distribution ρ from the distribution σ, see [64] for a review. In particular, for S(ρ|σ) = 0, the two distributions cannot be distinguished and agree, ρ = σ (up to a set of measure zero). For any other choice ρ = σ, the relative entropy is positive, S(ρ|σ) > 0. An interesting special case is when σ is the thermal density matrix with β = 1/T and Hamiltonian H. The partition function can be written as Z = e −βF with free energy F = E − T S, dF = −SdT − pdV . The relative entropy becomes where · ρ and · σ denote expectation values with respect to the density matrices ρ and σ, respectively, while S ρ and S σ denote the corresponding von Neumann entropies. Possible UV divergent contributions to the entanglement entropy are independent of the state and cancel between the first and second term in the second line of (A.3). (This is actually a general statement independent of the specific choice for σ made here.) Moreover, also possibly UV divergent contributions to the expectation values of energy, e. g. from the zero-point fluctuations of various modes cancel on the right hand side of (A.3). Note that for equal energy, H ρ = H σ , the relative entropy (A.3) equals the the difference of entropies. One may also consider a situation where ρ can be obtained from σ by a unitary operation. In that case S ρ = S σ and S(ρ|σ) = β( H ρ − H σ ), which is reminiscent of the definition of temperature starting from the microcanonical ensemble, dS = βdE. These considerations can be extended away from the case where σ is thermal in terms of the modular Hamiltonian K = − ln σ, S(ρ|σ) = K ρ − K σ ≥ 0.
Similar as for the von Neumann entropy we will perform the calculation in the replica formalism and obtain the relative entropy in (A.1) as a limit of the Rényi relative entropies [15] (the latter defined in [17]) With these definitions one has [15] S(ρ|σ) = lim (see discussion in section 2.4). For simplicity we will assume in the following that the field expectation values with respect to ρ and σ agree, i. e.φ (A.8) By subtracting (2.33) from (A.7) one obtains directly an expression for the relative Rényi entropy. In particular, we obtain for the quantum relative entropy by taking the limit N → 1, the expression S(ρ|σ) = 1 2 Tr

B Symmetries, anomalies and bosonization
In this appendix we recall how the Schwinger model corresponding to vector-like QED in two dimensions can be bosonized [65]. One proceeds via a discussion of gauge symmetries and associated quantum anomalies.

B.1 Symmetries, conservation laws and anomalies
Consider the Lagrangian The vector gauge symmetry is as usual or for left and right handed fermions One can define the regularization such that the vector gauge symmetry above is not anomalous. The associated Noether vector current is in our notation and it is classically and quantum theoretically conserved, ∂ µ J µ = 0. In contrast, the axial gauge transformation ψ → e iβγ 5 ψ,ψ →ψe iβγ 5 , A µ → A µ , B µ → B µ + 1 e ∂ µ β, (B. 5) or for left and right handed fermions is anomalous. The axial vector current J µ 5 = iψγ µ γ 5 ψ (B.7) is not conserved. There is an associated anomaly (see e. g. ref. [72], eq. (5.41)) where S = x L is essentially the microscopic action of the mass-less Schwinger model and we have added an irrelevant constant in the form of an integral over a free bosonic field φ. Decompose the gauge field as A µ = 1 e ∂ µ ζ + 1 e ∂ ν ξ νµ , (B.12) with two functions ζ and ξ. The Lagrangian is The idea is now to perform now a vector gauge transformation with α = −ζ (which is simple to do) and an axial vector gauge transformation with β = −ξ. One needs to do this in infinitesimal steps taking into account the presence of the axial vector potential B µ = ∂ µ (ξ + β). The integration of the anomaly gives −ξ 0 dβ e 2π µν F µν + e π ∂ µ ∂ µ (ξ + β) = − e 2π ξ µν F µν − e 2π ξ∂ µ ∂ µ ξ. (B.14) The fermionic functional integral becomes now a free one which contributes only an irrelevant factor to the partition function. Taking the anomaly into account leads therefore to This constitutes the bosonized form of the Schwinger model. In last step one can integrate out the gauge field by performing the Gaussian integral or, equivalently, solving the corresponding field equation.
In two dimensions, the gauge field A µ is not dynamical. One may chose the axial gauge A 1 = 0 and the field strength is F 10 = −F 01 = E 1 = ∂ 1 A 0 . (There is no magnetic field in two dimensions.) The equation of motion is in the bosonized theory and one can formally solve this