Dynamics of entanglement in expanding quantum fields
Jürgen Berges Stefan Floerchinger Raju Venugopalan
hep-th hep-ph nucl-th

Device-independent Certification of

One-shot Distillable Entanglement
    [
    
======================================================================





§ 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 <cit.>. Nowadays entanglement is seen as a key resource for quantum computation and quantum communication devices <cit.>. Entanglement in quantum field theory was investigated initially mainly with a view on black hole physics <cit.> and more recently in the context of holography <cit.>. 

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 ρ ,

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,

    S_A  = - Tr{ρ_A lnρ_A }.

(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 <cit.>. 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 <cit.>. The formalism has even been extended to discuss nonequilibrium dynamics <cit.>. Further, a replica method has also been developed to compute the relative entanglement entropy  between two density matrices <cit.>. 

Novel methods have also been developed for free field theories, in particular for massive (or massless) scalars and fermions <cit.>. 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. <cit.>, 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. <cit.> and develops it further. (See also <cit.> 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 <cit.>. 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 <cit.>.

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 <cit.> or the closely related Unruh effect <cit.> 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 <cit.>. 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 <cit.> and underlies, for example, the PYTHIA event generator <cit.>. 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. <cit.>, 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 <cit.>. 

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. <ref>, this interval can be described by a reduced density matrix corresponding to a trace over the complement region B as shown in eq. (<ref>).


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πτ) <cit.>. This effect is governed by an interplay of entanglement and the expansion dynamics. Our results in this direction are summarized in a recent letter <cit.>. 

We note that in the context of high energy hadron and heavy ion collisions, entanglement and entropy related ideas were discussed some time ago <cit.> as well as more recently <cit.>. 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 <cit.> and  another concerns experiments with ultracold atomic quantum gases where entanglement can actually be investigated directly <cit.>.

This paper is organized as follows. In section <ref>, we discuss Gaussian density matrices in quantum field theory and their entanglement properties. We start in section <ref> 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 <ref>, we discuss general mixed Gaussian states in terms of a variant of the Glauber-P representation <cit.> of the density matrix. A key formula here is an expression for the most general Gaussian density matrix in eq. (<ref>). In subsection <ref>, 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 <ref>. Subsequently, in subsection <ref>, 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. (<ref>) and (<ref>). 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 <ref> 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. (<ref>).

In section <ref> we illustrate the general formalism of section <ref> 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 <ref>. This is used in subsection <ref> for the computation of field-field and conjugate momentum correlation functions, and in subsection <ref> for the computation of the entaglement entropy.

In section <ref>, 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 <ref>, we investigate the Schwinger model and in subsections <ref> and <ref> 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 <ref>. 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 <ref> and provide a simple intuitive picture of how the apparent thermal character of these excitations arise. Our conclusions are stated in section <ref>. 

Appendix <ref> generalizes our formal discussion of entanglement entropies of Gaussian density matrices to relative entanglement entropies  and appendix <ref> provides a brief review of the  symmetries, quantum anomalies and bosonization of the Schwinger model for quantum electrodynamics in 1+1 dimensions.







§ 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

    Φ_m | ϕ⟩ = ϕ_m | ϕ⟩

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.[For instance, in the absence of any other components, we have ϕ_1 = ϕ and ϕ_2 = ϕ^*.] 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 <ref>. 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[For an introductory discussion about most general Gaussian states, see for example <cit.> (section 2.4) for the case of real bosonic fields.]     

    ⟨ϕ| ψ⟩  =     exp[ -1/2ϕ_m^* h_mnϕ_n + i/2ϕ_m^* λ_m + i/2λ^*_m ϕ_m +1/2ϕ^*_m  κ_m + 1/2κ_m^* ϕ_m ] 
      =     exp[ -1/2ϕ^†  h ϕ + i/2ϕ^†λ + i/2λ^†ϕ + 1/2ϕ^†κ + 1/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 (<ref>) being quadratic. 
Correspondingly, we have  

    ⟨ψ | ϕ⟩ = exp[ -1/2ϕ^†  h^†ϕ - i/2ϕ^†λ - i/2λ^†ϕ + 1/2ϕ^†κ + 1/2κ^†ϕ].

The norm is given by a functional integral over ϕ,

    ⟨ψ | ψ⟩  =     ∫ D ϕexp[ -1/2ϕ^† (h+h^†) ϕ + ϕ^†κ + κ^†ϕ] 
      =     exp[ -1/2ln (h+h^†) + 2 κ^† (h+h^†)^-1κ] ,

where we use ∫ D ϕ = ∏_n ∫ Dϕ_n with Dϕ_n= d Reϕ_n d Imϕ_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 by h̅, λ̅ and κ̅ is given by the functional integral

    ⟨ψ̅| ψ⟩ = ∫ D ϕ exp[ -1/2ϕ^† (h+h̅^†) ϕ + i/2ϕ^† (λ-λ̅) + i/2 (λ-λ̅)^†ϕ + 1/2ϕ^† (κ+κ̅) + 1/2 (κ+κ̅)^†ϕ].

Using two sets of eigenstates of (<ref>), whose eigenvalues we denote by ϕ_+ and ϕ_-, the pure-state density matrix can be written as

    ρ_λ,κ[ϕ_+, ϕ_-]   =  ⟨ϕ_+ | ψ⟩⟨ψ | ϕ_- ⟩
        =  exp[ -1/2ϕ_+^†  h ϕ_+    -1/2ϕ_-^†  h^†ϕ_- + i/2 (ϕ_+^† - ϕ_-^†) λ + i/2λ^† (ϕ_+ - ϕ_-) +1/2 (ϕ_+^† + ϕ_-^†) κ + 1/2κ^† (ϕ_+ + ϕ_-)  ].

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

    π_m = - i δ/δϕ_m  .

The representation (<ref>) implies the canonical commutation relations

    [ϕ_m, π_m] = i δ_mn,        [ϕ_m,ϕ_n] = [π_m,π_n]=0,

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[Deviations from Gaussianity can be treated in this formalism perturbatively but shall not be discussed any further here.]. We  can write such a mixed state density matrix in the form,

    ρ[ϕ_+, ϕ_-] = ∫ D λ Dκ  P[λ,κ]  ρ_λ,κ[ϕ_+, ϕ_-],

where ρ_λ,κ[ϕ_+, ϕ_-] is the pure state density matrix in (<ref>) (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. (<ref>) is closely related (although not identical) to the Glauber-P representation of a density matrix <cit.>.



In the following, for simplicity, we will take P[λ, κ] to also be of the Gaussian form,

    P[λ,κ] =  exp[ - [ λ^† - j^†, κ^† ]Σ^-1[ λ - j;     κ ] + const],

with a hermitian operator Σ that we take to be of the form[More general Gaussian forms for P[λ,κ] are possible but will not be needed for our construction.]

    Σ=Σ^†=[     Σ_a   i Σ_b; - i Σ_b     Σ_a ],       with  Σ_a=Σ_a^†,  Σ_b= Σ_b^† .

For this to be a properly normalizable probability distribution, the eigenvalues of Σ should be positive. One may also introduce the linearly transformed parameter fields,

    μ_m = i/√(2) (λ_m - j_m) + 1/√(2)κ_m,     ν_m = -i/√(2) (λ_m - j_m)+ 1/√(2)κ_m,

in terms of which the exponent in (<ref>) becomes diagonal,

    P[μ,ν] =  exp[ - μ^† (Σ_a-Σ_b)^-1μ - ν^† (Σ_a+Σ_b)^-1ν + const].

Substituting (<ref>) in (<ref>), 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 ϕ_+→ϕ_+-ϕ̅, ϕ_-→ϕ_- - ϕ̅ gives 

    ρ[ϕ_+, ϕ_-] = exp [    - 1/2[ ϕ_+^† - ϕ̅^†, ϕ_-^† - ϕ̅^† ][           h - Σ_a + Σ_b; - Σ_a - Σ_b         h^† ][ ϕ_+ - ϕ̅; ϕ_- - ϕ̅ ]
        + i/2 (ϕ_+^† - ϕ_-^†)j + i/2 j^† (ϕ_+ - ϕ_-)  ].

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 (<ref>) for suitable Gaussian P[λ,κ].

The trace of (<ref>) is given by

    Tr{ρ} = ∫ Dϕ ρ[ϕ,ϕ] = exp[ -1/2ln det(h+h^† - 2Σ_a) ] ,

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^† whereby

    ϕ̃_m = (δ_mn-P_mn) ϕ_n^+ = (δ_mn-P_mn) ϕ_n^-

are the fields we want to trace out, and

    ϕ̂_m^± = P_mnϕ_n^±

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

    ρ_R[ϕ̂_+, ϕ̂_-] = ∫ D ϕ̃ ρ[P ϕ̂^++(1-P)ϕ̃, P ϕ̂^-+(1-P)ϕ̃] .

(The projectors have been inserted here for clarity.) For the Gaussian density matrix (<ref>), one can formally perform the functional integral over ϕ̃ and obtain yet again a Gaussian form for the reduced density matrix,

    ρ_R[ϕ̂_+, ϕ̂_-] = exp [    - 1/2[ (ϕ̂_+^†  - ϕ̅^†) P, (ϕ̂_-^† - ϕ̅^† ) P ][           h -d_(++)   - Σ_a+Σ_b -d_(+-); - Σ_a - Σ_b -d_(-+)         h^† - d_(–) ][ P(ϕ_+ - ϕ̅); P(ϕ_- - ϕ̅) ]
        + i/2 (ϕ_+^† - ϕ_-^†)P j + i/2 j^† P (ϕ_+ - ϕ_-)  ].

We have used here the abbreviations

    d_(++) =      (h-Σ_a+Σ_b) (1-P) [(1-P)(h+h^†-2Σ_a)(1-P)]^-1(1-P) (h-Σ_a-Σ_b), 
    
    d_(+-) =      (h-Σ_a+Σ_b) (1-P) [(1-P)(h+h^†-2Σ_a)(1-P)]^-1(1-P) (h^† -Σ_a + Σ_b), 
    
    d_(-+) =      (h^† -Σ_a - Σ_b) (1-P) [(1-P)(h+h^†-2Σ_a)(1-P)]^-1(1-P) (h -Σ_a-Σ_b) , 
    
    d_(–) =      (h^†-Σ_a-Σ_b) (1-P) [(1-P)(h+h^†-2Σ_a)(1-P)]^-1(1-P) (h^† -Σ_a + Σ_b) .

The operator inversions used here are to be done within the subspace that is integrated out. Note that the reduced density matrix (<ref>) is of the same general structure as the original density matrix (<ref>); in particular, it still satisfies the hermiticity property (<ref>). However if one starts with the density matrix of a pure state where Σ_a=Σ_b=0, the reduced density matrix (<ref>) contains in general terms that mix ϕ_+ and ϕ_- (the terms d_(+-) and d_(-+) in (<ref>)) 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 (<ref>) one has for the field ϕ_m and the canonical conjugate momentum field π_m=-i δ/δϕ_m,

    ⟨ϕ_m ⟩ = ∫ D ϕ ϕ_m  ρ[ϕ, ϕ]/∫ Dϕ ρ[ϕ, ϕ] = ϕ̅_m ,       ⟨π_m ⟩ =  ∫ D ϕ ( -i δρ[ϕ_+, ϕ_-] / δϕ_+m)_ϕ_+=ϕ_-=ϕ/∫ Dϕ ρ[ϕ, ϕ] = j_m^*.

For the connected correlation functions ⟨ A B ⟩_c = ⟨ AB ⟩ - ⟨ A ⟩⟨ B ⟩ one finds

    ⟨ϕ_m ϕ_n^* ⟩_c =     [(h+h^† - 2 Σ_a)^-1]_mn, 
    ⟨π^*_m π_n ⟩_c =     [h - (h-Σ_a+Σ_b) (h+h^†-2Σ_a)^-1 (h-Σ_a-Σ_b) ]_mn
    
    =      [h^† - (h^†-Σ_a-Σ_b) (h+h^† - 2 Σ_a)^-1 (h^† - Σ_a+Σ_b) ]_mn
    
    
    =      [Σ_a - Σ_b + (h -Σ_a+Σ_b) (h+h^†-2Σ_a)^-1 (h^†-Σ_a+Σ_b)]_mn
    
    =      [Σ_a + Σ_b + (h^†-Σ_a-Σ_b)(h+h^†-2Σ_a)^-1 (h-Σ_a-Σ_b)]_mn , 
    ⟨ϕ_m π_n ⟩_c =      i [ (h+h^†-2Σ_a)^-1 (h -Σ_a - Σ_b) ]_mn, 
    ⟨π_n ϕ_m ⟩_c =      -i [ (h+h^†-2Σ_a)^-1 (h^† -Σ_a + Σ_b) ]_mn, 
    ⟨ϕ^*_m π^*_n ⟩_c =      i [ (h -Σ_a+Σ_b) (h+h^†-2Σ_a)^-1]_mn, 
    ⟨π^*_n ϕ^*_m ⟩_c =     - i [ (h^† -Σ_a - Σ_b) (h+h^†-2Σ_a)^-1]_mn .

These are compatible with the canonical commutation relations,

    ⟨ϕ_m π_n - π_n ϕ_m ⟩_c= i δ_mn, 
       ⟨ϕ^*_m π^*_n - π^*_n ϕ^*_m ⟩_c= i δ_mn .


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

    a_mn≡[(h+h^†-2Σ_a)^-1Σ_a]_mn =    ∑_k ⟨ϕ_m ϕ_k^* ⟩_c ⟨π^*_k π_n ⟩_c - 1/4δ_mn
        - 1/4∑_k,l,r⟨ϕ_m ϕ_k^* ⟩_c ⟨ϕ_k^* π_l^* + π^*_l ϕ_k^*⟩_c [⟨ϕϕ^* ⟩_c ]^-1_lr⟨ϕ_r π_n + π_n ϕ_r ⟩_c , 
    
    
    b_mn≡[(h+h^†-2Σ_a)^-1Σ_b]_mn =     i/4⟨ϕ_m π_n +  π_n ϕ_m ⟩_c  - i/4⟨ϕ_m ϕ_k^* ⟩_c ⟨ϕ_k^* π_l^* + π_l^* ϕ_k^* ⟩_c [ ⟨ϕϕ^* ⟩_c ]^-1_ln.

They follow with some algebra from (<ref>). 



 §.§ Entropy and entanglement entropy



In this subsection we will determine the Rényi entropy for the general Gaussian density matrix (<ref>) and extract from it the von Neumann entropy. The Rényi entropy is defined by (we assume standard normalization tr{ρ}=1),

    S_N(ρ) = 1/1-Nlntr{ρ^N },

and the von Neumann entropy follows from this <cit.> as a limit

    S(ρ) = - tr{ρ ln ρ}= lim_N→ 1 S_N(ρ).

For the purposes of the calculation that follows, it is clear[This can be seen directly by writing out the expressions for tr{ρ^N} in the functional integral formalism or more formally by noting that ϕ̅ and j can be changed by unitary transformations.] that one can drop ϕ̅ and j, as they do not enter tr{ρ^N }.  Using otherwise the general expression for ρ in (<ref>) leads after a straight-forward exercise in Gaussian integration to <cit.>

    Tr{ρ^N } = exp[ - 1/2Trln( M_N  ) ],

which contains the N dimensional cyclic matrix (with operator valued entries)

    M_N = (1+2a) 1_N - (a+b) Z_N - (a-b)Z_N^T  .

We have used here the abbreviations introduced in (<ref>). 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

    M_N = A_N(a,b)  A_N^T (a^T,-b^T) ,

with

    A_N(a,b) = ( √(14+a + b^2) + b + 1/2) 1_N - ( √(14+a + b^2) + b - 1/2) Z_N.

The determinant of the matrix A_N(a,b) is found to be 

    detA_N(a,b) = ( √(14+a + b^2) + b + 1/2)^N - ( √(14+a + b^2) + b - 1/2)^N.


Combining terms leads to a compact expression for the Rényi entropy,

    S_N(ρ) =  1/2(N-1){   Trln(  ( √(14+a+b^2)+b+1/2)^N - ( √(14+a+b^2)+b-1/2)^N )  
        +Trln(  ( √(14+a+b^2)-b+1/2)^N - ( √(14+a+b^2)-b-1/2)^N ) }.

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 <cit.>.

From the result above, one can directly obtain an expression for the von Neumann entropy,

    S =    1/2Tr{(√(1/4+a+b^2)+b+1/2) ln(√(14+a+b^2)+b+1/2) }
        - 1/2Tr{(√(1/4+a+b^2)+b-1/2) ln(√(14+a+b^2)+b-1/2) }
        + 1/2Tr{(√(14+a+b^2)-b+1/2) ln(√(14+a+b^2)-b+1/2) }
        - 1/2Tr{(√(1/4+a+b^2)-b-1/2) ln(√(14+a+b^2)-b-1/2) } .

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 (<ref>). The above expression simplifies for b=0 to

    S =    Tr{(√(14+a)+1/2) ln(√(14+a)+1/2) - (√(14+a)-1/2) ln(√(14+a)-1/2) } ,

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

    Δ^S_ϕϕ(x⃗-y⃗) =     ⟨ϕ(t,x⃗) ϕ(t,y⃗) ⟩ =  ∫_p⃗e^i p⃗ (x⃗- y⃗)/√(p⃗^2+M^2)[ 1/2 + n(p⃗) ], 
    Δ^S_ππ(x⃗-y⃗) =     ⟨π(t,x⃗) π(t,y⃗)  ⟩ = ∫_p⃗ e^i p⃗ (x⃗- y⃗)√(p⃗^2+M^2)[ 1/2 + n(p⃗) ],

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 (<ref>) is diagonal in momentum space. One finds

    S = ∫_p⃗{(n(p⃗)+1) ln(n(p⃗)+1) - n(p⃗) ln(n(p⃗))  } ,

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 <ref>.





 §.§ 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,

    ϕ_m → U_mnϕ_n ,          ϕ_m^* →ϕ_n^* (U^†)_nm , 
    π_m^* →  U_mnπ_n^* ,           π_m →π_n (U^†)_nm .

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 (<ref>) 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,

    [χ_m, χ^*_n] =  Ω_mn ,

where, symbolically,

    Ω = Ω^† = [     0   i 1; - i 1     0 ].

The transformations

    χ_m → S_mnχ_n,       χ_m^* →χ_n^* (S^†)_nm ,

such that 

    S Ω S^† = Ω ,

are compatible with the canonical commutation relations. This defines a symplectic transformation. Written in terms of the Lie algebra, S=exp[iθ^A J^A], the condition (<ref>) becomes

    Ω J^A = (J^A)^†Ω .

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

    χ_m^* = R_mnχ_n ,       χ_n = R^-1_nmχ_m^* = χ^*_m (R^†)_mn .

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 (<ref>). More specifically, one has

    R → U^* R (U^†) .

The matrix R transforms also by the symplectic transformations (<ref>),

    R → S^* R (S^-1) .

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 (<ref>) and the generator J^A transform under the unitary block diagonal transformations (<ref>) as

    S → U S U^†,        J^A → U J^A U^†.


We need to show that the symplectic transformations (<ref>) 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

    J^A = - (J^A)^* = Ω (J^A)^†Ω = - Ω (J^A)^T Ω .

There is now a representation of the Lie algebra specified by (<ref>) in terms of the operators (<cit.>, see also <cit.>)

    X^A = 1/2χ^*_m Ω_mn (J^A) _nkχ_k ,

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 (<ref>) 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 (<ref>) is invariant under the block diagonal unitary transformations (<ref>) 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

    Δ_mn = 1/2⟨χ_m χ_n^* + χ_n^* χ_m ⟩_c = 
    [                   ⟨ϕ_m ϕ_n^* ⟩_c         12⟨ϕ_m π_n + π_n ϕ_m ⟩_c; 12⟨π_m^* ϕ_n^* + ϕ_n^* π_m^* ⟩_c                ⟨π_m^* π_n ⟩_c  , ],

which is a key ingredient in our discussion of Gaussian states, transforms under symplectic transformations as

    Δ→ S Δ S^† .

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 <cit.> for a discussion) however that there  must exist a symplectic transformation that brings Δ to diagonal form,

    Δ = diag(λ_1,λ_2,…, λ_1,λ_2,…),

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 

    ΔΩ→ S Δ S^†Ω = S ΔΩ S^-1 ,

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 (<ref>) of the symmetric correlation matrix (<ref>) results in a very simple form for the quantities in (<ref>):

    a_ij = (λ_j^2 - 1/4) δ_ij ,        b_ij = 0 .

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 (<ref>) can be directly expressed in terms of the symplectic eigenvalues as 

    S = ∑_j {(λ_j+1/2) ln(λ_j+1/2) - (λ_j-1/2) ln(λ_j-1/2)  }.

Moreover, the symplectic eigenvalues λ_m follow as pairs of conventional eigenvalues ±λ_j of the combination ΔΩ. Following Sorkin (<cit.>, see also <cit.>), a further simplification can be obtained by considering the matrix 

    D=ΔΩ + 1/21 ,

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 (<ref>) simply as 

    S = ∑_j {ω_j^+ ln(ω_j^+) + ω_j^- ln(-ω_j^-) } ,

where the sum is over pairs of eigenvalues.  More simply,

    S = ∑_m ω_m ln(|ω_m|) = 1/2Tr{ D ln(D^2) }.

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,

    D_mn = 
    [       121 - i/2⟨ϕ_m π_n + π_n ϕ_m ⟩_c                      i ⟨ϕ_m ϕ_n^* ⟩_c;                     -i ⟨π_m^* π_n ⟩_c 121+ i2⟨π_m^* ϕ_n^* + ϕ_n^* π_m^* ⟩_c ]=
    [    -i ⟨ϕ_m π_n ⟩_c   i ⟨ϕ_m ϕ_n^* ⟩_c;  -i ⟨π_m^* π_n ⟩_c i ⟨π_m^* ϕ_n^* ⟩_c ].

Hence, the entropy associated with a Gaussian density matrix is fully determined from the set of connected correlation functions of (<ref>) evaluated in the domain of interest. It is understood that the operator trace in (<ref>) 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 <ref>. (See ref. <cit.> 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 (<ref>)) 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 (<ref>)) commute, i.e. [Δ^(ρ)Ω, Δ^(σ)Ω] = 0. One can then write the relative entropy as

    S(ρ | σ) = ∑_m ω_m^(ρ)( ln | ω_m^(σ) | - ln | ω_m^(ρ) | ),

where the sum goes over all pairs of simultaneous eigenvalues (ω_m^(ρ), ω_m^(σ)) of D^(ρ)=Δ^(ρ)Ω+1/21 and D^(σ)=Δ^(σ)Ω+1/21.





§ 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 <cit.>. 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 <ref>.

The approach we developed in section <ref> 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 <ref>, 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 a=(h+h^†-2Σ_a)^-1Σ_a defined in (<ref>) becomes

    a(x,y) =    ∫_-L/2^L/2 dz  Δ^S_ϕϕ(x-z)  Δ^S_ππ(z-y) - 1/4δ(x-y) 
    
    =    ∫_p,q{√(q^2+M^2)/4√(p^2+M^2)sin( 1/2(p-q)L )/1/2(p-q) - 2π/4δ(p-q)}  e^ipx-iqy.

Note that this is a nondiagonal matrix in momentum space for finite L. Only in the limit L→∞ does one obtain sin(1/2(p-q) L) / (1/2(p-q)) → (2π)δ(p-q) 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,

    φ(x) = φ^(s)(x) + φ^(a)(x) ,     φ^(s)(x)=φ(x)+φ(-x)/2 ,     φ^(a)(x)=φ(x)-φ(-x)/2 .

The symmetric part can be expanded into a Fourier series

    φ^(s)(x) = 1/L∑_n=-∞
     n  even^∞φ_n   e^inπ x/L,       φ_n = ∫_-L/2^L/2 dx  φ^(s)(x)   e^-inπ x/L   (n  even) .

In a similar fashion, one can expand the anti-symmetric part

    φ^(a)(x) = 1/L∑_n=-∞
     n  odd^∞φ_n   e^inπ x/L,       φ_n = ∫_-L/2^L/2 dx  φ^(a)(x)   e^-inπ x/L   (n  odd).

We can summarize this as

    φ(x) = 1/L∑_n=-∞^∞φ_n   e^inπ x/L,       φ_n = ∫_-L/2^L/2 dx  φ(x)  1/2[e^-inπ x/L + (-1)^n e^inπ x/L] .

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. (<ref>) becomes formally (with p=nπ/L and φ̂(q) = φ_n),

    φ(x) = 2 ∫dp/2π e^ipxφ̂(p) ,       φ̂(p) = ∫ dx   e^-ipx 1/2[φ(x) + (-1)^pL/πφ(-x) ]   .

It is useful to relate these expressions to the standard momentum space representation defined as usual by 

    φ(x) = ∫dq/2π e^ipxφ̃(p),       φ̃(p) = ∫ dx   e^-ipxφ(x).

One has for finite interval length L

    φ_n = ∫dp/2π sin(pL2-nπ2) [ 1/p-nπ/L + 1/p+nπ/L] φ̃(p),

and for very large L formally

    φ̂(p) =1/2[ φ̃(p) + (-1)^pL/πφ̃(-p) ].

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

    Δ_ϕϕ^S(x-y) = Δ_ϕϕ^(ss)(x,y) + Δ_ϕϕ^(aa)(x,y),       Δ_ππ^S(x-y) = Δ_ππ^(ss)(x,y) + Δ_ππ^(aa)(x,y) ,

with

    Δ_ϕϕ^(ss)(x,y) =    ⟨ϕ^(s)(x) ϕ^(s)(y)⟩ = ⟨12[ϕ(x) + ϕ(-x) ] 12[ ϕ(y)+ ϕ(-y)] ⟩ , 
    Δ_ϕϕ^(aa)(x,y) =    ⟨ϕ^(a)(x) ϕ^(a)(y)⟩ = ⟨12[ϕ(x) - ϕ(-x) ] 12[ ϕ(y) - ϕ(-y)] ⟩ ,

and similar for the canonical momentum field correlator. The cross terms like Δ_ϕϕ^(sa) vanish by parity symmetry. Further, Δ_ϕϕ^(ss)(x,y) is symmetric with respect to x→ -x as well as y → - y while Δ_ϕϕ^(aa)(x,y) is anti-symmetric with respect to parity. 

For the combined operator in (<ref>), it is natural to decompose a(x,y)=a^(ss)(x,y) + a^(aa)(x,y) with

    a^(ss)(x,y) =    ∫_-L/2^L/2 dz  Δ_ϕϕ^(ss)(x,z) Δ_ππ^(ss)(z,y) - 1/8δ(x-y) - 1/8δ(x+y) ,
    
    a^(aa)(x,y) =    ∫_-L/2^L/2 dz  Δ_ϕϕ^(aa)(x,z) Δ_ππ^(aa)(z,y) - 1/8δ(x-y) + 1/8δ(x+y) .

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 (<ref>), using (<ref>) (for n(p⃗)=0) one finds that 

    1/L[Δ^S_ϕϕ]_m(-n) =    1/L∫dp/2π sin(pL2-mπ2) sin(pL2-nπ2)
    [ 1/p-mπ/L + 1/p+mπ/L]
    [ 1/p-nπ/L + 1/p+nπ/L] 
    1/2√(p^2+M^2), 
    1/L[Δ^S_ππ]_m(-n) =    1/L∫dp/2π sin(pL2-mπ2) sin(pL2-nπ2)
    [ 1/p-mπ/L + 1/p+mπ/L]
    [ 1/p-nπ/L + 1/p+nπ/L] 
    √(p^2+M^2)/2 .

For a state with nonvanishing occupation number, one has to insert factors (1+2 n(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

    1/L[Δ^(ss)_ϕϕ]_m(-n) =    (-1)^m-n/2/8π L∫_-∞^∞ dp  [ 1/p-mπ/L + 1/p+mπ/L]
    [ 1/p-nπ/L + 1/p+nπ/L] 
    1/√(p^2+M^2)[1 - cos(pL)  ], 
    1/L[Δ^(aa)_ϕϕ]_m(-n) =    (-1)^m-n/2/8π L∫_-∞^∞ dp  [ 1/p-mπ/L + 1/p+mπ/L]
    [ 1/p-nπ/L + 1/p+nπ/L] 
    1/√(p^2+M^2)[1 + cos(pL)  ],

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 

    [1∓cos(pL)] = 12(1 ∓ e^ipL) + 12(1 ∓ e^-ipL) .

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) .

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 (<ref>) 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 (<ref>) and in its counterpart for Δ_ππ. 

The contribution from the branch cuts to the field-field correlator, as determined by the right hand side of (<ref>), is given by the integral

    (-1)^m-n/2L/π∫_ML^∞ dy y^2/[y^2+(mπ)^2][y^2+(nπ)^2]√(y^2-(ML)^2)(± e^-y-1) .

We shall first discuss the result without the exponential in the last bracket. This should be a good approximation in particular for (ML)≫ 1. For the opposite limit (ML)≪ 1, we will add the contribution from the exponential term as well. Performing the integral gives for m^2≠ n^2,

    -(-1)^m-n/2L/2π^3(m^2-n^2)[ |mπ|/√(m^2π^2 + (ML)^2)ln( √(m^2π^2+(ML)^2)+|mπ|√((mπ)^2+(ML)^2)-|mπ|) -
    |nπ|/√(n^2π^2 + (ML)^2)ln( √(n^2π^2+(ML)^2)+|nπ|√((nπ)^2+(ML)^2)-|nπ|)
     ],

while for m^2=n^2 the result is

    -(-1)^m-n/2L/4π^2 |m|((mπ)^2+(ML)^2)^3/2[ 2√(m^2π^2 (m^2π^2 + (ML)^2)) + (ML)^2 ln( √(m^2π^2+(ML)^2)+|mπ|√(m^2π^2+(ML)^2)-|mπ|) ].

Finally, for m=n=0 one obtains

    (-1)^m-n/2L/π(M L)^2.

Adding up these expressions gives the field-field correlation function in the massive case for ML≫ 1. 

In the massless case M=0, one can also perform all integrals and we obtain the field-field correlators,

    1/L[Δ^(ss)_ϕϕ]_m(-n) =     L ×(-1)^m-n/2/π^3[ c̃ĩ(mπ) - c̃ĩ(nπ) /m^2-n^2]       (m^2 ≠ n^2, m,n  even) 
    Si(mπ)/2π^2 m       (m^2=n^2≠ 0, m, n  even) 
    1/4π[ 3-2 γ - 2 ln(k L) ]      (m=n=0)
    
          0      (otherwise)
    1/L[Δ^(aa)_ϕϕ]_m(-n) =     L ×(-1)^m-n/2/π^3[ c̃ĩ(mπ) - c̃ĩ(nπ) /m^2-n^2]       (m^2 ≠ n^2, m,n  odd) 
    Si(mπ)/2π^2 m - 1/π^3 m^2       (m=n, m, n  odd) 
    
          -Si(mπ)/2π^2 m + 1/π^3 m^2       (m=-n, m, n  odd) 
    
          0      (otherwise)

For the zero mode with m=n=0, we have introduced an infrared momentum regulator k. Analytic continuation from nonvanishing m=n suggests k L = 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 by

    c̃ĩ(x) =    Ci(x) - ln(x) - γ =  ∫_0^x cos(t)-1/t =  
    ∑_k=1^∞(-x^2)^k/2k (2k)! ,
    Si(x) =    ∫_0^x sin(t)/t .


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) .

The contribution from the branch cut of the square root gives

    (-1)^m-n/2/π L∫_ML^∞ dy y^2/[y^2+(mπ)^2][y^2+(nπ)^2]√(y^2-(ML)^2)(1 ∓ e^-y).

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 (ML)≫ 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,

    (-1)^m-n/2/π Lln( 2Λ/e M).

We subtract now this contribution from the integral and obtain for the remainder when m^2≠ n^2,

    (-1)^m-n/2/π L   ( 1 - 1/2(m^2-n^2)π^2[ |mπ| √(m^2 π^2 + (ML)^2) ln( √(m^2π^2 + (ML)^2)+|mπ|√(m^2π^2+(ML)^2)-|mπ|) .. 
       . .                   - |nπ| √(n^2 π^2 + (ML)^2) ln( √(n^2π^2 + (ML)^2)+|nπ|√(n^2π^2+(ML)^2)-|nπ|)
     ] ),

and for m^2=n^2,

    (-1)^m-n/2/π L( 1/2 - 2m^2 π^2 + (ML)^2/4 |mπ| √(m^2 π^2 + (ML)^2)ln( √(m^2π^2 + (ML)^2)+|mπ|√(m^2π^2+(ML)^2)-|mπ|)  ).

By adding these expressions together, one obtains the field-field correlator in the massive case (ML)≫ 1. 

One can also perform the integrals for the canonical momentum correlation function in the massless limit to obtain,

    1/L[Δ^(ss)_ππ]_m(-n) =    1/L×(-1)^m-n/2/π[ m^2 c̃ĩ(mπ) - n^2 c̃ĩ(nπ)/m^2-n^2 + γ+ ln(Λ L) ]       (m^2≠ n^2, m, n  even)
    1/π[ c̃ĩ(mπ) + 12 mπ Si(mπ) + γ + ln(Λ L) ]       (m^2= n^2 ≠ 0, m, n  even)
    1/π[ γ + ln(Λ L)  ]       (m=n=0) 
    
    0       (otherwise)
    
    1/L[Δ^(aa)_ππ]_m(-n)=    1/L×(-1)^m+n/2/π[ m^2 c̃ĩ(mπ) - n^2 c̃ĩ(nπ)/m^2-n^2 + γ+ ln(Λ L) ]       (m^2≠ n^2, m, n  odd)
    1/π[ c̃ĩ(mπ)  + 12 mπ Si(mπ) -1 +  γ + ln(Λ L) ]       (m= n, m, n  odd)
    
    -1/π[ c̃ĩ(mπ)  + 12 mπ Si(mπ) -1 + γ + ln(Λ L) ]       (m= -n, m, n  odd)
    
    0       (otherwise)

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

    1/La_m(-n) =  ∑_l=-∞^∞{1/L[Δ^(ss)_ϕϕ]_m(-l)1/L[Δ^(ss)_ππ]_l(-n) + 1/L[Δ^(aa)_ϕϕ]_m(-l)1/L[Δ^(aa)_ππ]_l(-n)} - 1/41_m(-n) ,

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 (<ref>). Alternatively, one can determine the spectrum of eigenvalues of the matrix D in (<ref>) and use (<ref>). 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. <ref> 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 <cit.>

    S=c/3ln (Λ L)+const,

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. <cit.>.



§ 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

    ℒ =     - ψ̅γ^μ (∂_μ - i e A_μ) ψ - m ψ̅ψ - 14F_μν F^μν
    
     =      i ψ_L^* (∂_0 + ∂_1-ieA_0 - i e A_1) ψ_L + i ψ_R^* (∂_0 - ∂_1-ie A_0 + ie A_1) ψ_R 
        + i m ψ_L^* ψ_R - i m ψ_R^* ψ_L - 14F_μν F^μν.

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 <cit.>. (For completeness, we include a discussion of this feature of the model in Appendix <ref>.) The model can be reformulated as a scalar field theory with Lagrangian (see eq. (<ref>))

    ℒ = - 1/2∂_μϕ∂^μϕ - 1/2e^2/πϕ^2 - m e exp(γ)/2π^3/2cos(2√(π)ϕ + θ) .
 
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 <ref>, 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, (<ref>) 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. <cit.>. 




 §.§ General coordinates and background evolution



If one considers the free scalar field theory with mass M in arbitrary curved coordinates, one has the action

    S = ∫ d^2 x √(g){ -1/2 g^μν∂_μϕ∂_νϕ - 1/2 M^2 ϕ^2 }.

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 <cit.>

    ∂_τ^2 ϕ̅+ 1/τ∂_τϕ̅+ M^2 ϕ̅=0 .

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 <cit.>. In the limit of vanishing mass M→ 0, the corresponding two independent solutions to the equation of motion (<ref>) 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τ).


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/2M^2 ϕ̅^2 and “pressure” p̅=1/2(∂_τϕ̅)^2-1/2M^2 ϕ̅^2. One may check easily that the ideal fluid equation of motion,

    ∂_τϵ̅+ ϵ̅+ p̅/τ=0 ,

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 state p̅ = ϵ̅ for pure radiation in 1+1 dimensions.  More generally, depending on the initial conditions for ϕ̅ and ∂_τϕ̅, one can have initial conditions between p̅=-ϵ̅ and p̅=ϵ̅. For the solution in terms of J_0(Mτ), one initially has p̅ = - ϵ̅ with the result oscillating between this value and p̅=ϵ̅. 

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

    ℒ_ξ g_μν = ξ^α∂_α g_μν + ∂_μξ^α g_αν + ∂_νξ^α g_μα = 0 .

However, there is a conformal Killing vector field of this form which is a solution of 

    ℒ_ξ g_μν - λ g_μν= ξ^α∂_α g_μν + ∂_μξ^α g_αν + ∂_νξ^α g_μα- λ g_μν = 0 ,

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

    ϕ(τ, η) = ϕ̅(τ) +  φ(τ, η).

The fluctuation field has the equation of motion 

    ∂_τ^2 φ(τ, η) + 1/τ∂_τφ(τ, η) + (M^2-1/τ^2∂^2/∂η^2) φ(τ, η)=0 .

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

    φ(τ,η) = ∫d k/2π{ a(k) f(τ, |k|) e^ikη + a^†(k)   f^*(τ, |k|) e^-ikη},

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 

    ∂_τ^2 f(τ, k) + 1/τ∂_τ f(τ, k) + (M^2+k^2/τ^2) f(τ, k)=0.

The inner product of these mode functions provides the normalization condition 

    -i f(τ, k) τ∂_τ f^*(τ, k) + i f^*(τ, k) τ∂_τ f(τ, k) = 1 .

The canonical conjugate momentum field is 

    π_ϕ(τ,η) = τ∂_τϕ(τ,η) ,

and the canonical equal-time commutation relations [ϕ(τ,η), π_ϕ(τ,η^')] = i δ(η-η^') imply and are implied by the commutator [a(k), a^†(k^')] = 2πδ(k-k^'). 

Note that (<ref>) and (<ref>) do not fix the mode functions uniquely. The different possibile functions are related by the Bogoliubov transformations,

    f̅(τ, k) =    α(k) f(τ, k) + β(k) f^*(τ, k),         f̅^*(τ, k) = α^*(k) f^*(τ, k) + β^*(k) f(τ, k), 
    
    f(τ, k) =    α^*(k) f̅(τ, k) - β(k) f̅^*(τ, k),        f^*(τ, k) = α(k) f̅^*(τ, k) - β^*(k) f̅(τ, k),

with |α(k)|^2-|β(k)|^2=1. The corresponding creation and annihilation operators are related by

    a̅(k) =    α^*(k)a(k) - β^*(k) a^†(k),       a̅^†(k) = - β(k) a(k) + α(k) a^†(k),  
    
    a(k) =    α(k) a̅(k) + β(k) a̅^†(k),          a^†(k) = β^*(k) a̅(k) + α^*(k) a̅^†(k) .

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 a̅(k)| Ω̅⟩ = 0, in the latter case. Note however that the vacuum |Ω⟩ might contain excitations with respect to the operator a̅(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, (<ref>) has the two independent solutions J_i k(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

    f(τ, k) =  √(π)/2e^kπ/2 H^(2)_ik(Mτ),        f^*(τ, k) =  √(π)/2e^kπ/2 H^(1)_-ik(Mτ).

The set of mode functions f(τ,k) in (<ref>) is distinguished by being a superposition of positive frequency modes with respect to time t of standard Minkowski spacetime <cit.>. 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 (<ref>) becomes

    f(τ, k) → e^-ik ln (τ M/2)ie^kπ/2Γ(i k)/2√(π)  - e^ik ln (τ M/2)√(π) e^kπ/2 [cosh(k π)-1]/2Γ(1+i k)    (Mτ→ 0) .

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

    f̅(τ, k) = √(π)/√(2sinh(π k)) J_-ik(Mτ) ,       f̅^*(τ, k) = √(π)/√(2sinh(π k)) J_ik(Mτ) .

In the limit of vanishing mass Mτ→ 0, the mode function (<ref>) becomes

    f̅(τ, k) → e^-ik ln (τ M/2)√(π)/Γ(1-i k)√(2 sinhkπ)  = e^-i k ln (τ)-iθ (k, M)1/√(2k)    (Mτ→ 0) .

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) + (Γ(1-i k)) .

Note, in particular, that the factor multiplying k diverges in the formal limit M→ 0.

The Bogoliubov coefficients that connect the mode functions (<ref>) and (<ref>) are

    α(k) = √(e^π k/2 sinh(π k)),       β(k) = √(e^-π k/2 sinh(π k)) .



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, 

    ⟨ a^†(k) a(k^') ⟩_c     = n(k)  2π δ(k-k^') , 
    ⟨ a(k) a(k^') ⟩_c     = u(k)    2π δ(k+k^') , 
    ⟨ a^†(k) a^†(k^') ⟩_c     = u^*(k)   2π δ(k+k^') .

For example, the vacuum state with n(k)=u(k) = u^*(k) = 0 results in the correlation function

    ⟨ϕ(τ_1, η_1) ϕ(τ_2, η_2) ⟩_c 
     = ∫dk/2π  f(τ_1, k) f^*(τ_2, k)    e^ik(η_1-η_2) .

It is instructive to compare (<ref>) with the corresponding Minkowski space expression

    ⟨ϕ(x_1) ϕ(x_2) ⟩_c = 1/2∫_-∞^∞d^2 k/(2π)^2  e^-i k^0(x_1^0-x_2^0) + i k^1(x_1^1-x_2^1) (2π)  δ(k^2+M^2).

Employing the integral representation

    f(τ, k) e^ik η = i/2√(π)∫_-∞^∞ d z   e^-i M cosh(z) x^0 - i M sinh(z) x^1- i k z ,

where x^0=τcosh(η) and x^1=τsinh(η) are standard Minkowski coordinates, one obtains 

    ⟨ϕ(x_1) ϕ(x_2) ⟩_c  = 1/4πRe∫_-∞^∞ dz   e^-i M cosh(z) (x_1^0-x_2^0) - i M sinh(z) (x_1^1-x_2^1) .

This agrees with (<ref>) after substituting k^0=κcosh(z), k^1=-κsinh(z), thereby confirming that the mode functions in (<ref>) 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 

    ⟨ϕ(τ_1, k) ϕ^*(τ_2, k^') ⟩_c =     2πδ(k-k^')   f(τ_1, k) f^*(τ_2, k), 
    ⟨π^*(τ_1, k) π(τ_2, k^') ⟩_c =     2πδ(k-k^')  ḟ(τ_1, k) ḟ^*(τ_2, k), 
    ⟨ϕ(τ_1, k) π(τ_2, k^') ⟩_c =     2πδ(k-k^')    f(τ_1, k) ḟ^*(τ_2, k), 
    ⟨π^*(τ_1, k) ϕ^*(τ_2, k^') ⟩_c =     2πδ(k-k^')  ḟ(τ_1, k) f^*(τ_2, k).

where we have used the abbreviation ḟ(τ, k) = τ∂_τ f(τ,k). 
One can directly verify that with the above correlators, the matrix D in (<ref>) 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  (<ref>) 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,

    ⟨a̅^†(k) a̅(k^') ⟩_c     = n̅(k)  2π δ(k-k^') =  |β(k)|^2  2π δ(k-k^'), 
    ⟨a̅(k) a̅(k^') ⟩_c     = u̅(k)    2π δ(k+k^') = -α^*(k) β^*(k)    2π δ(k+k^'), 
    ⟨a̅^†(k) a̅^†(k^') ⟩_c     = u̅^*(k)   2π δ(k+k^') = -α(k) β(k)   2π δ(k+k^') .

In this alternative basis, the correlation functions do not look like those of an empty state but rather of one with occupation number n̅(k) = |β(k)|^2. From (<ref>) one obtains,

    n̅(k)= |β(k)|^2 = 1/e^2π k-1 .

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 (<ref>)
corresponds to a thermal spectrum with the time dependent temperature 

    T=1/(2πτ) .

Such a thermal interpretation is not possible for a nonvanishing mass M, but the fact that the quasiparticles defined by the mode functions f̅(τ, k) have a nonvanishing occupation number remains true. Of course, a thermal state would have vanishing entries for the other correlators appearing in (<ref>). The “off-diagonal occupation function” is

    u̅(k) = u̅^*(k)= -1/2 sinh(π k) .


One may use the relations in (<ref>) to express the correlation functions in  (<ref>) in the alternative basis,

    ⟨ϕ(τ_1, k) ϕ^*(τ_2, k^') ⟩_c =     2πδ(k-k^')  {f̅(τ_1, k) f̅^*(τ_2, k) [1+n̅(k)] + f̅^*(τ_1, k) f̅(τ_2, k) n̅(k) . 
                   .+ f̅(τ_1, k) f̅(τ_2, k) u̅(k) + f̅^*(τ_1, k) f̅^*(τ_2, k) u̅^*(k)   }, 
    ⟨π^*(τ_1, k) π(τ_2, k^') ⟩_c =     2πδ(k-k^')  {ḟ̅̇(τ_1, k) ḟ̅̇^*(τ_2, k) [1+n̅(k)]  + ḟ̅̇^*(τ_1, k) ḟ̅̇(τ_2, k) n̅(k) . 
                   .+ ḟ̅̇(τ_1, k) ḟ̅̇(τ_2, k) u̅(k) + ḟ̅̇^*(τ_1, k) ḟ̅̇^*(τ_2, k) u̅^*(k) }, 
    ⟨ϕ(τ_1, k) π(τ_2, k^') ⟩_c =     2πδ(k-k^')   {f̅(τ_1, k) ḟ̅̇^*(τ_2, k) [1+n̅(k)] + f̅^*(τ_1, k) ḟ̅̇(τ_2, k) n̅(k) .
                   . f̅(τ_1, k) ḟ̅̇(τ_2, k) u̅(k) + f̅^*(τ_1, k) ḟ̅̇^*(τ_2, k) u̅^*(k) }, 
    ⟨π^*(τ_1, k) ϕ^*(τ_2, k^') ⟩_c =     2πδ(k-k^')  {ḟ̅̇(τ_1, k) f̅^*(τ_2, k) [1+n̅(k)] + ḟ̅̇^*(τ_1, k) f̅(τ_2, k) n̅(k) . 
                   . + ḟ̅̇(τ_1, k) f̅(τ_2, k) u̅(k) + ḟ̅̇^*(τ_1, k) f̅^*(τ_2, k) u̅^*(k)  }.

From these relations, one obtains the equal time correlation functions by setting τ_1=τ_2=τ. Moreover, in the limit Mτ≪ 1, one can use (<ref>) for the mode functions f̅(τ, k). This gives the correlators 

    ⟨ϕ(τ, k) ϕ^*(τ, k^') ⟩_c =     2πδ(k-k^') 1/|k|{[12+n̅(k) ] + cos[ 2 k ln(τ) + 2 θ(k, M) ]  u̅(k)   }, 
    ⟨π^*(τ, k) π(τ, k^') ⟩_c =     2πδ(k-k^') |k| {[12+n̅(k) ] + cos[ 2 k ln(τ) + 2 θ(k, M) ]  u̅(k)   }, 
    ⟨ϕ(τ, k) π(τ, k^') ⟩_c =     2πδ(k-k^') {i2 - sin[ 2 k ln(τ) + 2 θ(k, M) ]  u̅(k)   }, 
    ⟨π^*(τ, k) ϕ^*(τ, k^') ⟩_c =     2πδ(k-k^') { -i2 + sin[ 2 k ln(τ) + 2 θ(k, M) ]  u̅(k)   }.

Note that the “off-diagonal occupation functions” u̅(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[2 k 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 <ref> 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 anti-symmetric components. The calculation proceeds by determining correlation functions as in eq. (<ref>) at fixed Bjorken time τ and in the finite interval. As previously, one can use eq. (<ref>) to determine the entanglement entropy. 





The correlation functions on the finite rapidity interval can be most conveniently obtained from the momentum space representations (<ref>) or (<ref>), together with the analog of the relation (<ref>) 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 τ)

    ϕ_n = ∫dk/2πsin(k Δη2-nπ2) [ 1/k-n π/Δη +  1/k+n π/Δη] ϕ(k) .
 

When one calculates correlators  in the discrete basis, such as  ⟨ϕ_n ϕ_-m⟩_c, from  (<ref>) using the kernel (<ref>), one observes that the terms proportional to the off-diagonal occupation number u̅(k) in (<ref>) 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 (<ref>), 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 numbers n̅(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 (<ref>). 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 entropy of an interval of length L is given by <cit.>

    S(T,l) = c/3ln( 1/π T ϵsinh(π L T) ) + const ,

with central charge c=1 in the present situation and where ϵ is a small length serving as ultraviolet regulator. Setting now L=τΔη and T=1/(2πτ) leads to the result

    S(τ, Δη) = c/3ln( 2τ sinh(Δη/2)/ϵ)+const.


Alternatively, one can obtain the result in (<ref>) for the expanding system from the Minkowski space  computations in our framework, as discussed in section <ref>. 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. <ref>. Using (<ref>), this leads to (<ref>) as well. This reasoning has been used for the derivation of the results first obtained in <cit.>.

The additive constant in eq. (<ref>) 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 <cit.>. In this case, the derivative with respect to Δη of the entanglement entropy gives

    ∂/∂ΔηS(τ, Δη) =  ∂ S/∂lnΔ z∂lnΔ z/∂Δη
        =  c_E(2Mτsinh(Δη/2))  (Δη/2)/2 ,

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)→ x K_1(2x)/4, which decays exponentially. 

We can employ the general expression for c_E(x) <cit.> for massive free bosons to compute dS/dΔη for the massless Schwinger model. The result is displayed in Fig. <ref> 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. <cit.>. The result is also shown in Fig. <ref> (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 (<ref>) 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 <cit.>.




 §.§ 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 (<ref>), when projected to the finite interval with the kernel (<ref>), 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,

    ρ = 1/Z e^-K,        Z = Tr  e^-K,

where K is the so-called modular or entanglement Hamiltonian defined as 

    K = ∫_Σ dΣ^μβ^ν T_μν  .

Equations (<ref>) and (<ref>) 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 (<ref>). Note also that β^μ=(1/T, 0) (in coordinates τ, η) is actually a conformal Killing vector according to (<ref>). 

This result can also be obtained as a limit of a more general result, as discussed already in Ref. <cit.>. A conformal field theory in the vacuum state in a region with a boundary formed by the intersection of two light cones (see figure <ref>), has a reduced density matrix of the form (<ref>) on any hypersurface in that region that has the same boundary. The modular Hamiltonian is a local expression given by <cit.> (see also <cit.>)

    K = ∫_Σ dΣ^μξ^ν T_μν  .

Here T_μν is again the energy-momentum tensor of excitations and ξ^ν is a vector field that can be written as 

    ξ^μ = 2π/(q-p)^2[ (q-x)^μ (x-p)(q-p) + (x-p)^μ (q-x)(q-p) - (q-p)^μ (x-p)(q-x) ] ,

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 (<ref>) is again of the same form as a density matrix of a thermal state if one identifies ξ^μ = 1/Tu^μ 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 <ref>, we have

    S(ρ^1 | ρ) = Δ⟨ K ⟩ = β^μ P_μ  ,

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 <ref> - <ref> and had discussed previously in <cit.>, 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. (<ref>). At very early Bjorken time τ, the interaction term ∼ m e 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.




§ RELATIVE ENTROPY



Extending the discussion of section <ref>, 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 σ) }.

The relative entropy measures in an information theoretic sense to which extend one can distinguish the distribution ρ from the distribution σ, see <cit.> 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

    σ= 1/Z e^-β H,

with β=1/T and Hamiltonian H. The partition function can be written as Z=e^-β F with free energy F=E-TS, dF=-SdT-pdV. The relative entropy becomes

    S(ρ | σ) =    Trρ ln ρ + β( ⟨ H ⟩_ρ - F_σ) 
    
    =     - S_ρ + S_σ + β( ⟨  H ⟩_ρ - ⟨ H ⟩_σ),

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 (<ref>). (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 (<ref>). Note that for equal energy, ⟨  H ⟩_ρ = ⟨ H ⟩_σ, the relative entropy (<ref>) 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 (<ref>) as a limit of the Rényi relative entropies <cit.> (the latter defined in <cit.>)

    S_N(ρ | σ) = 1/1-NlnTr{ρ σ^N-1}/Tr{ρ^N } = 1/1-NlnTr{ρ σ^N-1}- S_N(ρ).

With these definitions one has <cit.>

    S(ρ | σ) = lim_N→ 1 S_N(ρ | σ).

When both ρ and σ are Gaussian density matrices, they are fully characterized by their respective expectation values

    ⟨ϕ_m ⟩_ρ =    Tr{ρ ϕ_m } = ϕ̅_m^(ρ),        ⟨π_m ⟩_ρ = Tr{ρ π_m } = j_m^(ρ), 
    ⟨ϕ_m ⟩_σ =    Tr{σ ϕ_m } = ϕ̅_m^(σ),        ⟨π_m ⟩_σ = Tr{σ π_m } = j_m^(σ),

as well as connected two-field correlation functions, specified for ρ in terms of h^(ρ), h^(ρ)†, Σ_a^(ρ) and Σ_b^(ρ) and similar for σ in terms of h^(σ), h^(σ)†, Σ_a^(σ) and Σ_b^(σ) (see discussion in section <ref>).

For simplicity we will assume in the following that the field expectation values with respect to ρ and σ agree, i. e. ϕ̅_m^(ρ) = ϕ̅_m^(σ) and j_m^(ρ) = j_m^(σ). We make no such assumption about the connected correlation functions, however, and in fact keep them fully general. 

Using similar manipulations as in section <ref> we obtain

    1/1-N   ln Tr{ρ σ^N-1} = 
       1/2(N-1){Trln ((√(14+a^(ρ)+b^(ρ)2)+b^(ρ)+1/2 )(√(14+a^(σ)+b^(σ)2)+b^(σ)+1/2 )^N-1
                      - 
      (√(14+a^(ρ)+b^(ρ)2)+b^(ρ)-1/2 )(√(14+a^(σ)+b^(σ)2)+b^(σ)-1/2 )^N-1 )
                +Trln ((√(14+a^(ρ)+b^(ρ)2)-b^(ρ)+1/2 )(√(14+a^(σ)+b^(σ)2)-b^(σ)+1/2 )^N-1
                      - 
      (√(14+a^(ρ)+b^(ρ)2)-b^(ρ)-1/2 )(√(14+a^(σ)+b^(σ)2)-b^(σ)-1/2 )^N-1 )},

with

    a^(ρ) =     (h^(ρ) + h^(ρ)† - 2 Σ_a^(ρ) )^-1Σ_a^(ρ), 
    
    b^(ρ) =     (h^(ρ) + h^(ρ)† - 2 Σ_a^(ρ) )^-1Σ_b^(ρ), 
    
    a^(σ) =     (h^(σ) + h^(σ)† - 2 Σ_a^(σ) )^-1Σ_a^(σ), 
    
    b^(σ) =     (h^(σ) + h^(σ)† - 2 Σ_a^(σ) )^-1Σ_b^(σ).

By subtracting (<ref>) from (<ref>) 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/2Tr   {(√(14+a^(ρ)+b^(ρ)2)+b^(ρ)+12)  [ln (√(14+a^(σ)+b^(σ)2)+b^(σ)+12 )
                                      - ln (√(14+a^(ρ)+b^(ρ)2)+b^(ρ)+12 ) ]
        - (√(14+a^(ρ)+b^(ρ)2)+b^(ρ)-12)  [ln (√(14+a^(σ)+b^(σ)2)+b^(σ)-12 )
                                      - ln (√(14+a^(ρ)+b^(ρ)2)+b^(ρ)-12 ) ]
        +(√(14+a^(ρ)+b^(ρ)2)-b^(ρ)+12)  [ln (√(14+a^(σ)+b^(σ)2)-b^(σ)+12 )
                                      - ln (√(14+a^(ρ)+b^(ρ)2)-b^(ρ)+12 ) ]
        - (√(14+a^(ρ)+b^(ρ)2)-b^(ρ)-12)  [ln (√(14+a^(σ)+b^(σ)2)-b^(σ)-12 )
                                      - ln (√(14+a^(ρ)+b^(ρ)2)-b^(ρ)-12 ) ]}.

In the special case of b^(σ)=b^(ρ)=0 this simplifies to

    S(ρ| σ) = Tr   {(√(14+a^(ρ))+12)  [ln (√(14+a^(σ))+12 )  - ln (√(14+a^(ρ))+12 ) ]
        - (√(14+a^(ρ))-12)  [ln (√(14+a^(σ))-12 ) - ln (√(14+a^(ρ))-12 ) ]}.





§ SYMMETRIES, ANOMALIES AND BOSONIZATION



In this appendix we recall how the Schwinger model corresponding to vector-like QED in two dimensions can be bosonized <cit.>. One proceeds via a discussion of gauge symmetries and associated quantum anomalies. 



 §.§ Symmetries, conservation laws and anomalies

Consider the Lagrangian

    ℒ = - ψ̅( γ^μ∂_μ - i e A_μγ^μ - i e B_μγ^μγ_5 ) ψ.

The vector gauge symmetry is as usual

    ψ→ e^iαψ,     ψ̅→ψ̅e^-iα,      A_μ→ A_μ + 1/e∂_μα,      B_μ→ B_μ,

or for left and right handed fermions

    ψ_L→ e^iαψ_L,     ψ_R→ e^iαψ_R,     ψ_L^* →ψ_L^* e^-iα,     ψ_R^* →ψ_R^* e^-iα .

One can define the regularization such that the vector gauge symmetry above is not anomalous. The associated Noether vector current is in our notation

    J^μ = i ψ̅γ^μψ,

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∂_μβ ,

or for left and right handed fermions

    ψ_L→ e^iβψ_L,     ψ_R→ e^-iβψ_R,      ψ_L^* →ψ_L^* e^-iβ,     ψ_R^* →ψ_R^* e^iβ,

is anomalous. The axial vector current

    J_5^μ = i ψ̅γ^μγ_5 ψ

is not conserved. There is an associated anomaly (see e. g. ref. <cit.>, eq. (5.41))

    ∂_μ⟨ J^μ_5 ⟩ = - 1/2πϵ^μν e F_μν - 1/π e ∂_μ B^μ .

The Schwinger functional defined by

    e^iW= ∫ Dψ̅D ψ  e^i ∫ d^2 x ℒ

changes by an infinitesimal axial transformation δβ according to

    δ W= - ∫ d^2 x  δβ(x) ∂_μ⟨ j_5^μ(x) ⟩ =  ∫ d^2x {δβ(x) e/2πϵ^μν F_μν(x) + δβe/π∂_μ B^μ}.




 §.§ Bosonization

Consider now the partition function

    e^iW[J] = ∫ Dϕ Dψ̅D ψ D A   e^iS+i∫_x J^μ A_μ

where S=∫_x ℒ 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∂^νξ ϵ_νμ,

with two functions ζ and ξ. The Lagrangian is

    ℒ = -1/2∂_μϕ∂^μϕ - ψ̅(γ^μ∂_μ - i ∂_μζγ^μ - i ∂_μξγ^μγ_5 ) ψ - 14F_μν F^μν.

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πξ∂_μ∂^μξ.

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

    e^iW[J] = ∫ Dϕ D A   e^i ∫_x { - 1/2∂_μϕ∂^μϕ -1/4F_μνF^μν - e/2πξϵ^μν F_μν - e/2πξ∂_μ∂^μξ + J^μ A_μ} .

Use now A_μ=1/e∂^νξϵ_νμ which implies ϵ^μνe F_μν= -2 ∂_μ∂^μξ, perform a shift of variables ϕ+1/√(π)ξ→ϕ and use ∂^μξ = -e ϵ^μν A_ν to obtain

    e^iW[J] = ∫ Dϕ DA   e^i ∫_x {- 1/2∂_μϕ∂^μϕ -1/4F_μν F^μν -e/√(π)∂_μϕ ϵ^μνA_ν+ J^μ A_μ}.

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

    ∂_μ F^μν 
    = e/√(π)∂_μϕϵ^μν,

and one can formally solve this

    F^μν = e ( 1/√(π)ϕ + θ/2π) ϵ^μν,

where θ appears as an integration constant and is the vacuum angle. Using this in the partition function gives (for J=0)

    e^iW = ∫ Dϕ  e^i ∫_x {- 1/2∂_μϕ∂^μϕ - 1/2e^2/πϕ^2 + e^2/8π^2θ^2}.

This is now the partition function of a massive scalar particle with mass m=e/√(π). The vacuum angle is here just an irrelevant overall factor and therefore drops out.

For the massive Schwinger model the standard bosonization is more involved because the mass term transforms non-trivially under the axial gauge transformations. Order by order in a perturbative series in m one finds that the theory is equivalent to the bosonic theory <cit.>,

    ℒ = - 1/2∂_μϕ∂^μϕ - 1/2e^2/πϕ^2 - m e exp(γ)/2π^3/2cos(2√(π)ϕ + θ).
 
Here, γ is the Euler constant but note that the factor in front of the cosine term could be absorbed into a rescaling of the fermion mass m. This is now an interacting theory, and it depends on the vacuum angle θ.


This work is part of and supported by the DFG Collaborative Research Centre “SFB 1225 (ISOQUANT)”. R. V.'s research is supported by the U. S. Department of Energy Office of Science, Office of Nuclear Physics, under contracts No. DE-SC0012704. R. V. would like to thank ITP Heidelberg and the Alexander von Humboldt Foundation for support, and ITP Heidelberg for their kind hospitality.




Einstein:1935rr 
  A. Einstein, B. Podolsky and N. Rosen,
  “Can quantum mechanical description of physical reality be considered complete?,”
  Phys. Rev.  47, 777 (1935).
  
  
  

NielsenChuang
  M. A. Nielsen, I. L. Chuang, “Quantum Computation and Quantum Information”, (Cambridge University Press, Cambridge, 2000).


Bombelli:1986rw 
  L. Bombelli, R. K. Koul, J. Lee and R. D. Sorkin,
  “A Quantum Source of Entropy for Black Holes,”
  Phys. Rev. D 34, 373 (1986).
  
  
  


Srednicki:1993im 
  M. Srednicki,
  “Entropy and area,”
  Phys. Rev. Lett.  71, 666 (1993).
  
  
  
  


Callan:1994py 
  C. G. Callan, Jr. and F. Wilczek,
  “On geometric entropy,”
  Phys. Lett. B 333, 55 (1994).
  
  
  
  


Ryu:2006bv 
  S. Ryu and T. Takayanagi,
  “Holographic derivation of entanglement entropy from AdS/CFT,”
  Phys. Rev. Lett.  96, 181602 (2006); 
  
  
  
  


  
  “Aspects of Holographic Entanglement Entropy,”
  JHEP 0608, 045 (2006).
  
  
  
  


Holzhey:1994we
  C. Holzhey, F. Larsen and F. Wilczek,
  “Geometric and renormalized entropy in conformal field theory,”
  Nucl. Phys. B 424, 443 (1994).
  
  
  
  


Vidal:2002rm 
  G. Vidal, J. I. Latorre, E. Rico and A. Kitaev,
  “Entanglement in quantum critical phenomena,”
  Phys. Rev. Lett.  90, 227902 (2003).
  
  
  
  


Korepin:2004zz 
  V. E. Korepin,
  “Universality of Entropy Scaling in One Dimensional Gapless Models,”
  Phys. Rev. Lett.  92, 096402 (2004).
  
  
  


Calabrese:2004eu
  P. Calabrese and J. L. Cardy,
  “Entanglement entropy and quantum field theory,”
  J. Stat. Mech.  0406, P06002  (2004).
  
  
  
  


Calabrese:2009qy 
  P. Calabrese and J. Cardy,
  “Entanglement entropy and conformal field theory,”
  J. Phys. A 42, 504005 (2009).
  
  
  
  

Calabrese:2005in 
  P. Calabrese and J. L. Cardy,
  “Evolution of entanglement entropy in one-dimensional systems,”
  J. Stat. Mech.  0504, P04010 (2005).
  
  
  
  


Calabrese:2007mtj 
  P. Calabrese and J. Cardy,
  “Entanglement and correlation functions following a local quench: a conformal field theory approach,”
  J. Stat. Mech.  0710, no. 10, P10004 (2007).
  
  
  
  


Calabrese:2016xau 
  P. Calabrese and J. Cardy,
  “Quantum quenches in 1 + 1 dimensional conformal field theories,”
  J. Stat. Mech.  1606, no. 6, 064003 (2016).
  
  
  
  


Lashkari:2014yva
  N. Lashkari,
  “Relative Entropies in Conformal Field Theory,”
  Phys. Rev. Lett.  113, 051602 (2014).
  
  
  
  


Lashkari:2015dia 
  N. Lashkari,
  “Modular Hamiltonian for Excited States in Conformal Field Theory,”
  Phys. Rev. Lett.  117, 041601 (2016).
  
  
  
  


Ruggiero:2016khg
  P. Ruggiero and P. Calabrese,
  “Relative Entanglement Entropies in 1+1-dimensional conformal field theories,”
  JHEP 1702, 039 (2017).
  
  
  
  


Casini:2005zv 
  H. Casini and M. Huerta,
  
  J. Stat. Mech.  0512, P12012 (2005).
  
  
  
  


Casini:2007bt 
  H. Casini and M. Huerta,
  “Analytic results on the geometric entropy for free fields,”
  J. Stat. Mech.  0801, P01012 (2008).
  
 
  
  


Casini:2009sr 
  H. Casini and M. Huerta,
  “Entanglement entropy in free quantum field theory,”
  J. Phys. A 42, 504007 (2009).
  
  
  
  
 
 
Koksma:2010zi 
  J. F. Koksma, T. Prokopec and M. G. Schmidt,
  “Entropy and Correlators in Quantum Field Theory,”
  Annals Phys.  325, 1277 (2010).
  
  
  
  


Koksma:2011fx 
  J. F. Koksma, T. Prokopec and M. G. Schmidt,
  “Decoherence and Dynamical Entropy Generation in Quantum Field Theory,”
  Phys. Lett. B 707, 315 (2012).
  
  
  
  
 

Sorkin:2012sn 
  R. D. Sorkin,
  “Expressing entropy globally in terms of (4D) field-correlations,”
  J. Phys. Conf. Ser.  484, 012004 (2014).
  
  [arXiv:1205.2953 [hep-th]].
  
  
  

Saravani:2013nwa 
  M. Saravani, R. D. Sorkin and Y. K. Yazdi,
  “Spacetime entanglement entropy in 1 + 1 dimensions,”
  Class. Quant. Grav.  31, no. 21, 214006 (2014).
  
  
  
  


Birrell:1982ix 
  N. D. Birrell and P. C. W. Davies,
  “Quantum Fields in Curved Space,” (Cambridge University Press, Cambridge, 1984).
  
  
  


Mukhanov:2007zz 
  V. Mukhanov and S. Winitzki,
  “Introduction to quantum effects in gravity,” (Cambridge University Press, Cambridge, 2007).
  
  


Gelis:2006yv 
  F. Gelis and R. Venugopalan,
  “Particle production in field theories coupled to strong external sources,”
  Nucl. Phys. A 776, 135 (2006).
  
  
  
  


Gelis:2006cr 
  F. Gelis and R. Venugopalan,
  “Particle production in field theories coupled to strong external sources. II. Generating functions,”
  Nucl. Phys. A 779, 177 (2006).
  
  
  
  


Hawking:1974sw 
  S. W. Hawking,
  “Particle Creation by Black Holes,”
  Commun. Math. Phys.  43, 199 (1975); 
  Erratum: [Commun. Math. Phys.  46, 206 (1976)].
  
  
  


Unruh:1976db 
W. G. Unruh,
“Notes on black hole evaporation,”
Phys. Rev. D 14, 870 (1976).




Becattini:1995if 
  F. Becattini,
  “A Thermodynamical approach to hadron production in e+ e- collisions,”
  Z. Phys. C 69, no. 3, 485 (1996).

Castorina:2007eb 
  P. Castorina, D. Kharzeev and H. Satz,
  “Thermal Hadronization and Hawking-Unruh Radiation in QCD,”
  Eur. Phys. J. C 52, 187 (2007).

Andronic:2008ev 
  A. Andronic, F. Beutler, P. Braun-Munzinger, K. Redlich and J. Stachel,
  “Thermal description of hadron production in e+e- collisions revisited,”
  Phys. Lett. B 675, 312 (2009).


Becattini:2009fv 
  F. Becattini and R. Fries,
  “The QCD confinement transition: Hadron formation,” in “Relativistic heavy ion physics”
  Landolt-Bornstein 23, 208 (2010).
  
  
  
  
  
Becattini:2010sk 
  F. Becattini, P. Castorina, A. Milov and H. Satz,
  “A Comparative analysis of statistical hadron production,”
  Eur. Phys. J. C 66, 377 (2010).

Andersson:1983ia 
B. Andersson, G. Gustafson, G. Ingelman and T. Sjöstrand,
“Parton Fragmentation and String Dynamics,”
Phys. Rept.  97, 31 (1983).


Andersson:1998tv 
  B. Andersson,
  “The Lund model,”
  Camb. Monogr. Part. Phys. Nucl. Phys. Cosmol.  7, 1 (1997).
  
  

Sjostrand:2006za 
  T. Sjostrand, S. Mrenna and P. Z. Skands,
  “PYTHIA 6.4 Physics and Manual,”
  JHEP 0605, 026 (2006).
  
  
  
  


Sjostrand:2007gs 
  T. Sjostrand, S. Mrenna and P. Z. Skands,
  “A Brief Introduction to PYTHIA 8.1,”
  Comput. Phys. Commun.  178, 852 (2008);
  
  
  
  



Fischer:2016zzs 
  N. Fischer and T. Sjöstrand,
  “Thermodynamical String Fragmentation,”
  JHEP 1701, 140 (2017).
  
  
  
  


Bialas:1999zg 
  A. Bialas,
  “Fluctuations of string tension and transverse mass distribution,”
  Phys. Lett. B 466, 301 (1999).
  
  
  
  
  

Berges:2017zws 
  J. Berges, S. Floerchinger and R. Venugopalan,
  “Thermal excitation spectrum from entanglement in an expanding quantum string,”
  arXiv:1707.05338 [hep-ph].
  


Elze:1994hj 
  H. T. Elze,
  “Entropy, quantum decoherence and pointer states in scalar 'parton' fields,”
  Phys. Lett. B 369, 295 (1996).
  
  
  
  
  

Elze:1994qa 
  H. T. Elze,
  “Quantum decoherence, entropy and thermalization in strong interactions at high-energy. 1: Noisy and dissipative vacuum effects in toy models,”
  Nucl. Phys. B 436, 213 (1995).
  
  
  
  


Muller:2011ra 
  B. Müller and A. Schäfer,
  “Entropy Creation in Relativistic Heavy Ion Collisions,”
  Int. J. Mod. Phys. E 20, 2235 (2011).
  
  
  
  

Kovner:2015hga 
  A. Kovner and M. Lublinsky,
  “Entanglement entropy and entropy production in the Color Glass Condensate framework,”
  Phys. Rev. D 92, no. 3, 034016 (2015).


Akkelin:2016rhm 
  S. V. Akkelin,
  “Hot origin of the Little Bang,”
  Eur. Phys. J. A 53, no. 12, 232 (2017).
  
  
  
  

Kharzeev:2017qzs 
  D. E. Kharzeev and E. M. Levin,
  “Deep inelastic scattering as a probe of entanglement,”
  Phys. Rev. D 95, no. 11, 114008 (2017).

Martens:2017cvj 
  J. C. Martens, J. P. Ralston and J. D. T. Takaki,
  “Quantum tomography for collider physics: Illustrations with lepton pair production,”
  arXiv:1707.01638 [hep-ph].

Shuryak:2017phz 
  E. Shuryak and I. Zahed,
  “Regimes of the Pomeron and its Intrinsic Entropy,”
  arXiv:1707.01885 [hep-ph].


Muller:2017vnp 
  B. Müller and A. Schäfer,
  “Why does the thermal model for hadron production in heavy ion collisions work?,”
  arXiv:1712.03567 [nucl-th].
  


Baker:2017wtt 
  O. K. Baker and D. E. Kharzeev,
  “Thermal radiation and entanglement in proton-proton collisions at the LHC,”
  arXiv:1712.04558 [hep-ph].
  

FisherInformation
H. Strobel, W. Muessel, D. Linnemann, T. Zibold, D. B. Hume, L. Pezzè, A. Smerzi and M. K. Oberthaler,
“Fisher Information and entanglement of non-Gaussian spin states,”
Science 345, 424 (2014).

Susceptibilities
P. Hauke, M. Heyl, L. Tagliacozzo, and P. Zoller,
“Measuring multipartite entanglement via dynamic susceptibilities,”
Nature Physics 12, 778 (2016).

MeasuringEntanglementEntropy
R. Islam, R. Ma, P. M. Preiss, M. E. Tai, A. Lukin, M. Rispoli, and M. Greiner,
“Measuring entanglement entropy in a quantum many-body system,”
Nature 528, 77 (2015).

EPRsteering
P. Kunkel, M. Prüfer, H. Strobel, D. Linnemann, A. Frölian, T. Gasenzer, M. Gärttner and M. K. Oberthaler, 
“Spatially distributed multipartite entanglement enables Einstein-Podolsky-Rosen steering of atomic clouds,”
arXiv:1708.02407.


Glauber:1963tx 
  R. J. Glauber,
  “Coherent and incoherent states of the radiation field,”
  Phys. Rev.  131, 2766 (1963).
  
  
  


Mandel:1995seg 
  L. Mandel and E. Wolf,
  “Optical Coherence and Quantum Optics,”
  (Cambridge University Press, Cambridge, 1995).
  


Berges:2015kfa 
  J. Berges,
  “Nonequilibrium Quantum Fields: From Cold Atoms to Cosmology,”
  arXiv:1503.02907 [hep-ph].
  
  

Holevo1999
  A. S. Holevo, M. Sohma and O. Hirota,
  “Capacity of quantum Gaussian channels,”
  Phys. Rev. A 59, 1820 (1999).


Mukunda:1988cb 
  N. Mukunda, R. Simon and G. Sudarshan,
  “Gaussian pure states in quantum mechanics and the symplectic group,”
  Phys. Rev. A 37, 3028 (1988).
  
  
  

ARL2014
  G. Adesso, S. Ragy and A. R. Lee,
  “Continuous Variable Quantum Information: Gaussian States and Beyond,”
  Open Syst. Inf. Dyn. 21, 1440001 (2014).
  
  

Simon1999
R. Simon, S. Chaturvedi and V. Srinivasan,
  “Congruences and canonical forms for a positive matrix: Application to the Schweinler-Wigner extremum principle,”
  J. Math. Phys. 40, 3632 (1999).
  


Vedral:2002zz
  V. Vedral,
  “The role of relative entropy in quantum information theory,”
  Rev. Mod. Phys.  74, 197 (2002).
  
  
  


Coleman:1975pw 
  S. R. Coleman, R. Jackiw and L. Susskind,
  “Charge Shielding and Quark Confinement in the Massive Schwinger Model,”
  Annals Phys.  93, 267 (1975).
  
  
  

Goykhman:2015sga 
  M. Goykhman,
  “Entanglement entropy in 't Hooft model,”
  Phys. Rev. D 92, 025048 (2015).


Loshaj:2011jx 
  F. Loshaj and D. E. Kharzeev,
  “LPM effect as the origin of the jet fragmentation scaling in heavy ion collisions,”
  Int. J. Mod. Phys. E 21, 1250088 (2012).
  
  
  
  


Hebenstreit:2013baa 
  F. Hebenstreit, J. Berges and D. Gelfand,
  “Real-time dynamics of string breaking,”
  Phys. Rev. Lett.  111, 201601 (2013).
  
  
  
  

 

Casini:2011kv 
  H. Casini, M. Huerta and R. C. Myers,
  “Towards a derivation of holographic entanglement entropy,”
  JHEP 1105, 036 (2011).
  
  
  
  
  

Arias:2016nip 
  R. Arias, D. Blanco, H. Casini and M. Huerta,
  “Local temperatures and local terms in modular Hamiltonians,”
  Phys. Rev. D 95, no. 6, 065005 (2017).
  
  
  
  


Candelas:1978gf 
P. Candelas and J. S. Dowker,
“Field Theories On Conformally Related Space-times: Some Global Considerations,”
Phys. Rev. D 19, 2902 (1979).





Anomalies
F. Bastianelli and P. van Nieuwenhuizen, 
“Path Integrals and Anomalies in Curved Space,” 
(Cambridge University Press, Cambridge, 2006).


