Towards spacetime entanglement entropy for interacting theories

Entanglement entropy of quantum fields in gravitational settings is a topic of growing importance. This entropy of entanglement is conventionally computed relative to Cauchy hypersurfaces where it is possible via a partial tracing to associate a reduced density matrix to the spacelike region of interest. In recent years Sorkin has proposed an alternative, manifestly covariant, formulation of entropy in terms of the spacetime two-point correlation function. This formulation, developed for a Gaussian scalar field theory, is explicitly spacetime in nature and evades some of the possible non-covariance issues faced by the conventional formulation. In this paper we take the first steps towards extending Sorkin’s entropy to non-Gaussian theories where Wick’s theorem no longer holds and one would expect higher correlators to contribute. We consider quartic perturbations away from the Gaussian case and find that to first order in perturbation theory, the entropy formula derived by Sorkin continues to hold but with the two-point correlators replaced by their perturbation-corrected counterparts. We then show that our results continue to hold for arbitrary perturbations (of both bosonic and fermionic theories). This is a non-trivial and, to our knowledge, novel result. Furthermore we also derive closed-form formulas of the entanglement entropy for arbitrary perturbations at first and second order. Our work also suggests avenues for further extensions to generic interacting theories.


Introduction
Around 1952, Rudolf Peierls found himself "looking for a new expression of some of the basic rules of quantum mechanics, namely the formulation of commutation laws of relativistic field theory" [1]. There already existed consistent formulations of relativistic quantum field theory, developed decades earlier by Heisenberg and Pauli [2,3], as well as more recent interaction picture frameworks developed by Schwinger [4], Tomonaga [5] and others.

JHEP11(2020)114
Despite this, Peierls was worried about the Lorentz invariance of these formulations. He wanted to find an alternative approach in which the Lorentz invariance of the formulation is manifestly evident. An unsatisfactory feature of the existing approaches was that a Hamiltonian and its associated canonical variables were needed, thus tying the formalisms to a choice of time coordinate or frame.
Around this time, Peierls succeeded in deriving the covariant formulation of a relativistic quantum field theory that he had in mind. His formulation, explained in [6], determines a general rule for the commutator at any pair of spacetime points. This formulation uses Green functions of the theory and makes no reference to a Hamiltonian. More precisely, the spacetime commutator (also known as the Peierls bracket or Pauli Jordan function) is the difference between the retarded and advanced Green functions: ∆ = G R − G A . Peierls also held much correspondence about his work with other prominent physicists of his time. When his paper came out, many were interested in it but also expressed skepticism. 1 Nonetheless, the Peierls bracket has stood the test of time and is an important framework for quantum field theory that shows that its time-dependence can be captured in a Lorentz invariant manner. It has also enabled important advances in algebraic quantum field theory [7,8], quantum field theory in curved spacetime [9][10][11], and quantum gravity [12][13][14]. We will return to the Peierls bracket later in the main text. For now we turn to an important topic in the current century and to the work of another physicist adamant about manifest Lorentz invariance. Entanglement entropy has become an increasingly important and useful topic in the past few decades. In this work, we will mainly be interested in it in the context of quantum field theory. The physical content of a quantum field theory is often expressed in terms of the set of its n-point correlations. Under certain circumstances, either intentionally (often the case in condensed matter systems) or unintentionally (often the case in gravitationally interesting systems such as black hole and cosmological spacetimes), one may not have access to the quantum field in the entire spacetime in which it lives. As a result of this limited access, we lose the information that was contained in the n-point correlations involving points both in the accessible and inaccessible regions. This is where entanglement entropy comes into play. Entanglement entropy is a measure of this loss of information.
The prototypical and perhaps most interesting background spacetime in which we would not have access to all the correlation information of a quantum field is that of a black hole. The event horizon of the black hole is a boundary that separates the degrees of freedom we have access to and those that we do not. In fact, the birth of the concept of the entanglement entropy of a quantum field occurred exactly as it was studied by Rafael Sorkin in this context, back in 1983 [15]. In this first work on the topic, a scalar field in a black hole spacetime was considered. The field data on a Cauchy surface was divided into the interior and exterior (of the black hole) components, and the entanglement entropy S = tr ρ ext log ρ −1 ext (1.1) 1 Notably, Pauli who was very much interested in this work was skeptical about it and stated that he thought "it looks like a cemetery for the lorentzinvariant formfactor theory in its present form"! [1] JHEP11(2020)114 was computed. ρ ext is what remains after the interior data is traced out. It was found that the resulting entropy is proportional to the spatial area of the event horizon in units of the cutoff, i.e. A/ 2 in 4d. Without a cutoff the entropy is infinite. Because of this scaling property, entanglement entropy was proposed by Sorkin to be a candidate for the microscopic origin of black hole entropy, which also scales like the area of the event horizon.
The question of the microscopic origin of black hole entropy has been the subject of intense research ever since the semiclassical calculations of Bekenstein and Hawking [16][17][18]. As mentioned, these calculations showed that the black hole entropy is proportional to the spatial area of the event horizon. An understanding of the microscopic degrees of freedom giving rise to this entropy 2 has, however, been elusive. Since Sorkin's proposal, entanglement entropy has been taken seriously as a/the source of the entropy. Numerous studies have already been made on the connection between entanglement entropy and black hole entropy [16][17][18][19], but there is not yet a definitive answer regarding whether or not it is the fundamental description of black hole entropy.
A difficulty in this regard is that the formula is used by defining the density matrix relative to a spatial Cauchy hypersurface Σ. Similarly, the cutoff that renders the entropy finite and quantifies it is defined relative to this hypersurface. Harking back to the worries of Peierls in 1952, such a construction of entanglement entropy lacks Lorentz invariance. Sorkin was also worried about this and in 2012 [20] he derived a covariant definition of spacetime entropy (including but more general than entanglement entropy), using none other than the Peierls bracket. The other major ingredient in this definition is the spacetime correlation function. 3 This definition is reviewed in section 3. This brings us to the topic of the current paper. The definition given in [20] is limited to the Gaussian theory. A question that arises is whether a covariant spacetime definition of entropy can also be found for non-Gaussian and interacting theories. At first sight it might seem like a difficult task. In the Gaussian theory, one has the ease of working with only the two-point function. For a non-Gaussian and/or interacting field theory, one may have to consider the higher n-point functions as well. In this case, there are two possibilities for us to succeed in the generalization we seek: i) we find a natural generalization of the formula in [20] (possibly a(n) (in)finite set of formulas) that includes the contributions from all n-point functions, or ii) we find that not all n-point functions contribute to the entropy.
In the present work, we find that up to first order in perturbation theory possibility (ii) holds: two-point functions suffice to capture the entropy. Beyond first order, we are faced with possibility (i). In particular we show explicitly how all higher-order correlators are needed to fully capture the entropy at second order. We begin by considering a generic non-Gaussian scalar theory with quartic perturbations to the quadratic Gaussian density matrix 2 For example a statistical mechanical understanding of this entropy of the form S ∼ ln(number of microstates).
3 See [21] for a prescription for obtaining the correlation function from the Peierls bracket as well.

JHEP11(2020)114
and we arrive at a formula for the entropy that is essentially the same as the one in [20] up to first order. The difference in our case is that the perturbation-corrected two-point correlation function enters the formula rather than the unperturbed one from the Gaussian theory. We then prove that the same formula captures the first order contribution to the entanglement entropy for any perturbation away from a (bosonic or fermionic) Gaussian theory. This finding may also point towards a deeper understanding of the information content of the (entanglement) entropy of a quantum field. At second order, we find that a part of the entropy is still captured by the same formula as the previous two orders, but the full second order contribution is not captured by it. We show that generically the full second order correction contains contributions of all higher-order correlation functions. Furthermore, we derive closed-form formulas for the entropy for arbitrary perturbations at first and second order. These findings can facilitate several extensions of the current work to more complicated interacting theories, in particular the extension of the formula in [20].

General quantum field theory for a real scalar theory
In this section we briefly discuss aspects of a general quantum field theory with a real scalar field. Since we are interested in a spacetime formulation we will work in the Heisenberg picture, loosely using the language of its axiomatizations such as the Haag-Kastler [22,23] and Wightman axioms [24,25]. We will not attempt to make this presentation rigorous in any way. For further details, see [26][27][28][29][30]. Given a spacetime (M, g) with manifold M and metric g, for any region R ⊂ M we can associate a unital -algebra of observables called A R . These algebras must satisfy certain properties. For example for any subregion R 1 ⊂ R 2 the corresponding algebras are nested A R 1 ⊂ A R 2 . In order to ensure causality, the algebras for any causally disjoint (with respect to the metric g) regions R 1 and R 2 must commute, i.e.
It is also required that if R 1 contains a Cauchy surface Σ of R 2 , then This requirement is essentially about the existence of dynamics. The logic is that observables outside the Cauchy surface Σ are determined by dynamical time evolution of the theory. In other words, operators with support in R 2 outside R 1 can be constructed as (fairly complicated) functions of operators in A R 1 or just A Σ . Other requirements include some form of automorphism under Poincaré transformations in Minkowski spacetime, or a compatibility axiom in more general spacetimes [29]. The set of assignments R → A R is sometimes called a net of local algebras. For a nested family of regions that satisfy ∪ n R n = M, the net of local algebras can be used to construct the full spacetime algebra by ∪ n A Rn = A M . 4

JHEP11(2020)114
In addition to the algebra of observables, we also need a positive functional · : A → C, such that 1 = 1, usually called a functional state. In more conventional language this functional assigns expectation values to any given observable in A, and can therefore be used to the construct correlation functions of the QFT.
It is useful to make the above construction more explicit. 5 We will from now on think of the algebras A R as being generated by a Hermitian scalar field φ(x) and the expectation values from a given vacuum state |0 . We can then generate the Hilbert space (up to other superselection sectors) by the span of states of the form 6 In order to construct this Hilbert space, we do not need the full spacetime algebra A M . Thinking classically, we would expect that operators on a Cauchy surface Σ would be enough to generate the full Hilbert space. But surprisingly according to the Reeh-Schleider theorem [31,32], in quantum field theory much less is needed; any open set R is sufficient to generate a dense subspace of H.

Spacetime entanglement entropy
We now discuss what the entanglement entropy of a state ρ is relative to a spacetime region R. In standard quantum mechanics we have a Hilbert space H associated to a Cauchy surface Σ, where the algebra A Σ acts irreducibly on H, and a "vacuum" state ρ. Since the full spacetime M is the domain of dependence of Σ, we have that A M = A Σ . Therefore A M acts irreducibly on H and we can think of ρ as the global state. Now consider the subregion R ⊂ M. In general, we cannot expect the subalgebra A R to act irreducibly on the Hilbert space H. However, imagine that we can find another Hilbert space H R where it does act irreducibly and a density matrix ρ R in H R such that In such a case we can define the entropy of ρ relative to the spacetime region R as [20] We can think of S(R) as a generalized entropy relative to any spacetime region R. In special cases this can be interpreted as the conventional entanglement entropy; if there exists a Cauchy surface Σ of the spacetime M, such that R ∩ Σ is a Cauchy surface for R, then S(R) corresponds to the standard bipartite entanglement entropy (see figure 1). Even in such cases, S(R) has the advantage over S(R ∩ Σ) of allowing for a covariant spacetime cutoff. The cutoff plays a central role in the definition of the entropy in field theory; without it the entropy would be infinite.
). For such regions the entropy S(R) corresponds to the entropy S(Σ = Σ R ), which captures the bipartite entanglement entropy between Σ R and its complement in Σ.

Quantum Peierls brackets
Let φ(x) be a real scalar field, A an algebra generated by φ(x), H a Hilbert space on which A acts irreducibly, and ρ a density matrix on this Hilbert space. From the associative product on A we can define a general commutator of the generators φ(x) such that [φ(x), φ(y)] ∈ A.
In particular (using the notation φ i for φ(x)) we generally have something of the form Note that not all of these operators are scalars; they can be higher-rank tensors as well.
..in must naturally satisfy a set of constraints, such as causality and the Jacobi identity.
The spacetime commutators (2.6) can be thought of as the commutators of operators in the Heisenberg picture, or as quantum Peierls brackets. 7 Note that products of operators at a single point are often singular. In reality, these products must be defined in a more careful manner. In 1 + 1 dimensional CFTs this is done using point-split regularization and a generalized normal-ordering prescription. More generally, this can be remedied by smearing.

JHEP11(2020)114
Let us define the Wightman correlation tensors for any positive integer n ∈ N as (2.8) In general we expect these to determine the theory fully. Due to hermiticity of φ i , we have the conditions W where is complex conjugation. For the two-point correlation matrix this condition is the standard hermiticiy condition W † (2) = W (2) . Taking the expectation value of the commutator we find that where iΩ ij is the expectation value of the right hand side of (2.6). For a free theory we have that commutators of the Heisenberg operators are c-numbers [33,34] (see also appendix C) and thus W ij The higher-rank correlation tensors are then fixed by W ij (2) through Wick's theorem. However for more general interacting theories the higher-rank correlation tensors will be necessary to specify the theory.
It is now clear that the entanglement entropy of ρ relative to a region R, S(R), can be expressed in terms of W ij (2) for a free theory. But for an interacting theory, we would expect the higher correlators W i 1 ...in (n) to contribute as well (2.12) As a final remark, we note that not all theories are completely specified from the set of n-point functions of local operators. For example, in gauge theories non-local operators can exist such as line and surface operators [35].

Computation using spacetime correlators ∆ and W
In this subsection we review the main results in [20]. We consider a free scalar field theory in a Gaussian state. That the theory is free means that the equations of motion are linear. As we have mentioned already, the spacetime commutator in this case is a c-number.
The two field theory correlators that we need to define the entropy are

JHEP11(2020)114
As we have seen in (2.11), ∆ can be obtained from the anti-symmetric part of W . Other useful relations between these two are where R is real and symmetric. The field theory problem can be divided into a series of calculations each involving a single degree of freedom q and p (see [20]). It suffices to derive the entropy for one degree of freedom and the full entropy becomes a sum of this quantity for all degrees of freedom. The most generic Gaussian density-matrix in a q-basis is where A, B, C are constant parameters and N is a normalization constant that is fixed by the condition trρ = 1. For (3.4) we get N = A/π. We will work with (3.4) in this and the next subsection.
Since the entropy S is dimensionless and invariant under unitary transformations (S(ρ) = S(U ρU † )), we expect it to depend on the symplectic invariant qq pp −(Re qp ) 2 = det R/ det ∆. The strategy will be to obtain ∆ and R in terms of A, B, C and then relate Spec{∆ −1 R} = {iσ, −iσ} to these parameters, thereby getting S(σ). Also notice that σ 2 = det R/ det ∆. In terms of q and p, ∆ and R are We need the correlators qq , qp , pq , and pp to obtain det R/ det ∆ in terms of A and C. B drops out of this expression [20], so the entropy will not depend on it. We carry out the expectation value computations with B set to zero in appendix A.1. The results are Putting the expressions above together we have which gives us the relation we need to switch between σ and the parameters A and C. At this stage, a result from [19] can be used where the entropy is obtained in terms of Using (3.8) we can rewrite (3.9) in terms of σ as

JHEP11(2020)114
From here, further simplifying algebra is used to express (3.10) in terms of the eigenvalues So far we have been considering the contribution of one degree of freedom to the entropy. The final step in our review is to include the contribution of all the degrees of freedom of the scalar field. As mentioned earlier, this can be done by summing over the contributions of each individual degree of freedom. Labelling now the eigenvalues as λ, our final expression for the entropy of a Gaussian scalar field is For some applications of (3.12)-(3.13), where the conventional spatial area law scaling with the UV cutoff is obtained, see [36,37].

Computation using replica trick
In this subsection, we derive the result of the previous subsection using the replica trick with the density matrix (3.4). See also [38]. The replica trick [39] is where n is analytically continued to 1. We must therefore compute tr (ρ n ). The trace is by definition given by (with the notation q n+1 = q 1 ) We can rewrite (3.15) more conveniently as as before, and β = 1 2 A 1 + 2C/A + A + C . In terms of these new parameters the n th power of the normalization constant is N n = A π n = |1 − µ n | n β π n . We next make the change of variables y i ≡ q i − µq i+1 . In terms of the Jacobian matrix where we have substituted in the value of N n in the last step. Finally, we insert (3.17) into (3.14) and take the limit to obtain The entropy expression (3.18) is exactly that of (3.9), and from (3.9) we can follow the same procedure as in section 3.1 that leads to (3.13).

Entropy for a perturbed theory
As a first extension of previous work in the Gaussian theory, we consider a non-Gaussian theory. As in the previous section, we will compute the entropy independently via both the correlators and the replica trick.

Computation using spacetime correlators ∆ and W
In this subsection we conjecture that to first order in perturbation theory, the entropy formula in the Gaussian case stays the same, but with W and ∆ replaced by those of the non-Gaussian theory. We will later return to the justification of this conjecture. 9 With this assumption, we carry out the correlator calculations of section 3.1 to first order in perturbation theory. We consider the density matrix (in a block) which is the most generic (symmetric in q and q ) quartic perturbation to a Gaussian density matrix in a block. By imposing that this density matrix is normalized (to first order in λ i ) we get

JHEP11(2020)114
With this matrix we obtain the following correlators (always to first order in λ i ) where the subscript g refers to the correlators in the Gaussian case. We can express the remaining correlators in terms of q n as and The details of the computations of the above correlators can be found in appendix A.2. The computation of σ 2 (recall from equation (3.8)) reduces to and using (4.3) we get For later comparison we will compute the following quantity to first order in perturbation theory. This is nothing but the quantity appearing in the Gaussian spacetime entanglement entropy formula, computed perturbatively for a non-Gaussian theory. We will investigate whether the entanglement entropy for the non-Gaussian theory can be obtained from the formula (3.9), by replacing µ with µ correlation .

Computation using replica trick
We will now compute the entanglement of the perturbed state ρ qq = q|ρ|q = N e −A/2(q 2 +q 2 )−C/2(q−q ) 2 − λ 1 q 4 +q 4 2 +λ 2 (q 3 q +qq 3 )+λ 3 q 2 q 2 (4.9) JHEP11(2020)114 using the replica trick and perturbation theory. The trace of ρ n can be parameterized in the following way where the quartic perturbation tensor is given by 11) and the quadratic coefficient matrix is given by with periodic convention of the indices i = n+1 = 1 (or in other words, the Kronecker delta is over Z n and indices defined modulo n). In the second line we have used A = β(1 − µ) 2 and C = 2βµ.
It is useful to think of the above as the partition function of an interacting n-particle system and define the following non-interacting correlation functions The Gaussian partition function is given by where in the last expression above we used (4.23). The equation (4.10) can now be expressed in terms of familiar perturbation theory The O(T m ) correction requires a 4m-point · correlation function. In order to compute the 4m-point correlator integrals, it is easier to first compute the integral and then compute the correlators by differentiation . (4.17)

JHEP11(2020)114
Note that this is nothing but Wick's theorem with q α q β = M −1 αβ . We can however reduce this sum significantly as the terms have a lot of internal symmetry. Since M is a symmetric matrix, any permutation of indices within a M −1 factor will give an identical term. This leads to a multiplicity of 2 and thus in total 2 m , if we use a convention that removes this redundancy. Furthermore there is an m!-fold multiplicity since each term is invariant under all permutations of the M −1 factors (or more precisely, their indices). Therefore, we only need to consider a subset of permutations G ⊂ S 2m . (4.18) Here we have defined the quotient of groups where S (i,j) 2 ⊂ S 2m is the subgroup of swaps between the i'th and j'th elements only and S (i,j) p,2 ⊂ S 2m (for odd i = j) are pairwise swaps of the (i, i + 1)'th elements with the (j, j +1)'th elements. All these subgroups are isomorphic to S 2 . The number of inequivalent permutations is thus |G| = (2m)! 2 m m! . For example for a 4-point function (2m = 4), there are 3 different permutations.
As computed in appendix B, the inverse of M can be expressed in different ways This matrix is dense, but for our purposes we only need the diagonal , (4.21) and the next-to-diagonal elements .
Here we have used that M ii = M 11 and M ii+1 = M 12 for any i. The determinant is given by Substituting (4.18) and (4.11) into (4.15), to first order in T we get

JHEP11(2020)114
Since we know the explicit n-dependence of tr (ρ n ), we can find the entropy using the replica trick from equation (3.14) This is the entropy of the non-Gaussian theory up to first order in perturbation theory. We would like to investigate whether there exists a µ replica such that the entanglement entropy can be expressed in the following way We can express this quantity perturbatively as Inserting this into (4.26) and expanding to first order we get (4.28) By setting this equal to equation (4.25) and solving for δµ i we find the following solutions (4.29) We can thus parametrize the entanglement entropy in terms of the parameter (4.30) By comparing this to equation (4.8), we see that In other words µ replica , and thus also the entanglement entropy, is actually a spacetime quantity and can be computed from spacetime correlators. The entanglement entropy formula (3.9) appears to still hold beyond the Gaussian theory (where Wick's theorem holds).

Generalization to arbitrary perturbations and higher orders
Our result above, for the first order contribution to the entanglement entropy, continues to hold for arbitrary perturbations. Namely, the first order (Gaussian or non-Gaussian) perturbation of the entanglement entropy around a Gaussian state (bosonic or fermionic, and involving an arbitrary number of degrees of freedom) is completely captured by the two-point correlation function of the perturbed state. In this section we prove a general theorem stating this, and we discuss some of its implications.

Proof of the general case
We consider a bosonic 10 system with N degrees of freedom, classical phasespace V and separable Hilbert space H V . We can always choose a basis of linear observables 11 ξ = (ξ 1 ,ξ 2 , . . . ,ξ 2N ) ≡ (q 1 ,p 1 , . . . ,q N ,p N ), which satisfy the canonical commutation relations 12 where a, b = 1, . . . , 2N and In the context of field theory, we may refer toq i andp i rather asφ i andπ i , i.e., as field operators and their conjugate momenta. When taking the continuum limit, these operators become operator-valued distributions, as commonly considered in algebraic quantum field theory.

A phase space decomposition V = A ⊕ B into subsystems A and B induces a Hilbert space decomposition
i.e., they only probe the state from the perspective of the subsystem A which could represent a causal diamond in spacetime.
Let us consider a one-parameter family |ψ of pure quantum states in H V , where we require |ψ 0 to be a Gaussian state. This induces a one-parameter family of possibly mixed states by tracing out the degrees of freedom of H B . By construction, ρ 0 is Gaussian and can therefore be written as ρ 0 = e −Ĥ A /Z with Z = tr e −Ĥ A , whereĤ A is known as the "modular Hamiltonian", which is quadratic for Gaussian states. This means there exists a symmetric, real bilinear form h ab with 13 The entanglement entropy is given as a function S = − tr H A (ρ log ρ ), which can be computed at = 0 from the two-point correlation function via (3.12)-(3.13). 14 The first 10 While we focus on the bosonic case, our derivation applies analogously to fermions by treating them in a parallel way, as reviewed in [53]. 11 They are called linear observables, because classical qi : V → R and pi : V → R are linear maps on phasespace. 12 For fermions,ξ a represents Majorana modes that satisfy canonical anticommutation relations, i.e., 14 Compared to our earlier perturbations in terms of λ's, one could think of as a perturbation along a particular directionr in parameter space: (λ1, λ2, λ3) = r. This is essentially writing the parameters in spherical coordinates.

JHEP11(2020)114
law of entanglement entropy [40][41][42] (see a brief derivation in the next subsection) states that we have SinceĤ A is quadratic (for both bosons and fermions), it only probes the two-point correlation function of ρ . It therefore does not matter if |ψ is perturbed in a Gaussian or non-Gaussian way, as the entanglement entropy at linear order will only be sensitive to the two-point correlation function of the state ρ which cannot be distinguished from a Gaussian state. We find 15 Put differently, for any highly non-Gaussian family ρ , we could define a Gaussian oneparameter familyρ , such that tr(ξ a Aξ b A ρ ) = tr(ξ a Aξ b Aρ ). The entanglement entropy would not be able to distinguish between the two at linear order around ρ 0 . This last point ensures that the formula for the linear perturbation of the entanglement entropy will be the same as expanding formulas (3.12)-(3.13) for Gaussian states, which explains the finding of the previous sections. 16 This result can also be interpreted geometrically as indicated in figure 2. The manifold of pure Gaussian states M is a submanifold of the projective Hilbert space P(H V ) consisting of all pure states. At a Gaussian state ψ 0 , we have Due to the fact that these tangent spaces are equipped with a natural Hilbert space inner product enables us to decompose any linear perturbation δψ 0 of a state ψ 0 into two pieces, namely which are orthogonal to each other. One can show that (δψ 0 ) captures the full change of the covariance matrix [47], while (δψ 0 ) ⊥ only feels the change of higher order correlation functions. 17 In this picture, the entanglement entropy as a function of the pure state ψ at ψ 0 only depends on the Gaussian perturbation In summary, we proved the observation of the previous section in full generality (for any non-Gaussian perturbation of a Gaussian state) using the properties of pure and mixed Gaussian states. While we focused on the bosonic case, we commented at the relevant places how we would arrive at the same conclusion for the linear perturbation of the entanglement entropy around fermionic Gaussian states.

Analysis of second order
We would like to compute the higher order perturbations of the entanglement entropy. We define ρ = e −Ĥ , whereĤ is the modular Hamiltonian. We have As discussed in the previous subsection, the first order perturbation is well-known from the first law of entanglement entropy and given by where we used that tr(ρ ) = tr(ρ ) = 0 as the trace of a mixed state is constant and equal to 1. Note that we did not need to worry about the ordering of the inverse of ρ , because the expression was inside a trace where we can use cyclicity. If ρ were not invertible, we mean JHEP11(2020)114 by ρ −1 the Penrose-Moore pseudo-inverse, where we invert on the orthogonal complement of its kernel, without needing to change the resulting equations. As previously discussed, the first order perturbation S 0 away from a Gaussian state (whereĤ 0 is quadratic) will only depend on the change of the two point function of the state ρ .
The second order perturbation can then be computed as We see that this perturbation consists of two pieces. The first one will be the second order change tr(ρ Ĥ ) of the expectation value of the modular Hamiltonian, which will again only depend on the change of the two-point function when perturbing a Gaussian state. The second piece is less straightforward and given by −tr(ρ 2 0 ρ −1 0 ). It is expected that this second piece will probe genuine non-Gaussian properties. Let us take a closer look at this term. Pulling out the derivative from (ρ ) 2 we get 14) The second term above has vanishing contribution to S 0 : The first term yields Since ρ and ρ σ do not commute, unless = σ, we have Note that A σ, is not in general a density matrix, nor even Hermitian. We can however use a polar decomposition to decompose it into a positive Hermitian operator ρ σ, and a unitary operator U σ, . These operators are

JHEP11(2020)114
Here ρ σ, = exp(H σ, ) is a two-parameter family of states where ρ 0,0 = ρ 2 0 is Gaussian, but H σ, is in general some interacting Hamiltonian. Recall that 18 Note that we can write U σ, = exp(ih σ, ), where h σ, is some Hermitian operator that can be expanded as polynomials of ξ a . Thus expanding the exponential of U σ, in the above expression, results in even higher-order correlators. From the higher-order correlators in (5.20) it is evident that non-Gaussian terms will contribute to the second order S 0 term. We can derive very general formulas for the entanglement entropy of a system with a single degree of freedom with a generic perturbation. In particular consider the Gaussian state in (3.4) with the most general perturbation In appendix D we compute S 0 and S 0 for a single degree of freedom with arbitrary perturbations. In particular the second term of (5.13), (S 0 ) 2 = − tr (ρ ) 2 ρ −1 which contains the non-Gaussian contributions, is given by Note that f ∂ ∂J 1 , ∂ ∂J 2 is a differential operator constructed out of the power expansion of the analytic function f (x, y), where the variables are replaced by differential operators. 18 There is a very compact way to compute ρ0 from the 2-point correlation function derived in [46]. For this, it is best to express the two-point function in terms of the linear complex structure J as in [44]. This gives ρ0 = exp (−q abξ aξb )/Z where q = −iω arccoth(iJ) for bosons and q = ig arctanh(iJ) for fermions, where ω is the inverse symplectic form (commutation relations), g the inverse metric (anti-commutation relations) and the equations should be understood as matrix equations. Note that this formula only applies to the Gaussian state ρ0 = ρ0,0.

JHEP11(2020)114 6 Summary, conclusions, and outlook
In this paper we reviewed the main results of [20], where a covariant definition of entropy is proposed for Gaussian theories in terms of the spacetime two-point correlation function.
With the aim of generalizing these ideas to the interacting case, we sketched general properties of interacting theories and discussed expectations for the resulting entropy formula. As a first step towards this goal, we considered generic perturbations away from the Gaussian theory. We found that to first order the same formula holds also for non-Gaussian theories with the correlators replaced by their perturbation-corrected versions. Naively one would expect higher-order correlators to start contributing as we move away from the Gaussian theory, but this is not the case up to first order of perturbation theory. Only from second order on do we see that higher-order correlators contribute and genuine non-Gaussian effects are observed. The formula in (5.20), expresses these genuine non-Gaussian effects partially in terms of spacetime correlators, and thereby extends the formula in [20] to interacting theories (up to second order). In future work, it would be interesting to express this expression entirely and explicitly in terms of spacetime correlators. Such a formulation in terms of spacetime correlators gives us a covariant definition of entropy for interacting theories that is especially useful for the study of entanglement entropy in general curved spacetimes.
Furthermore, we also derived closed-form formulas for the corrections of the entanglement entropy in appendix D. Path-integral methods such as those in [48][49][50] would be a natural setting to generalize these results to the field theory case.
As a final remark, one might also be able to make an interesting connection between our findings and the horizon molecules in [51]; there may be a relation between horizon molecules being links and the entropy up to first order being fully described by two-point functions.

JHEP11(2020)114
the University of Alberta. MZ acknowledges financial support provided by FCT/Portugal through the IF programme, grant IF/00729/2015.

A Computation of p and q correlators
We outline in this appendix the detailed computations of the correlators both for the Gaussian and non-Gaussian theories used in the main text.

A.1 Correlators for Gaussian theory
We start from the density matrix where N g ≡ A π is the normalization constant. With this we obtain Finally, for the pp correlator,

JHEP11(2020)114
With this density matrix we obtain the following correlators (always to first order in λ i ) q n = tr(q n ρ) = N dq q n e −Aq 2 −(λ 1 +2λ 2 +λ 3 )q 4 where the subscript g refers to the correlator in the Gaussian case, , we can further write With this we can write the remaining correlators in a more compact form as follows.

B Inverse of M-matrix
We are interested in computing the inverse of the matrix (4.12) in order to compute correlation functions (4.17). There are several ways to do this. One way is to note that the matrix M can be thought of as the matrix elements of a Hamiltonian for a quantum particle on a ring with n sites and nearest-neighbour hopping where c j and c † j are creation and annihilation operators, respectively. Since this Hamiltonian has translation symmetry, momentum is a good quantum number. We can therefore find the eigenvalues by a discrete Fourier transform over Z n where by G R and G A we will denote the retarded and advanced Green functions. The Peierls bracket of the unperturbed theory is given by [6] {φ 0 (x), φ 0 (y)} = ∆(x, y), (C. 5) where ∆ = G R − G A . Quantizing this theory, we will find c-number commutators of the Heisenberg opeators. Using this Peierls bracket together with (C.4), we find the Peierls bracket of the interacting field to be of the form {φ(x), φ(y)} = ∆(x, y) + λ × (non c-number terms) + O(λ 2 ). (C.6) Here by "non c-number" terms we mean those containing polynomials of φ 0 . Quantizing the interacting theory requires some care as the normal ordering of the non c-number terms is ambiguous and must be chosen in a way that avoids any anomalies. The above illustrates the connection between adding interactions in the Lagrangian and the appearance of non c-number corrections in the interacting field commutator. Generally, theories with c-number commutators are (generalized) free theories [34] and those that have non c-number terms correspond to interacting theories. Actually, if any truncated (connected) Wightman function of order 2n (n > 1) is zero then the theory turns out to a be generalized free theory [33]. Thus interacting theories have non c-number commutators and all their even truncated Wightman functions are non-zero.

D S 0 and S 0 for a single degree of freedom and arbitrary perturbations
Consider a general perturbation for any analytic symmetric function of two variables f (x, y). For the derivation we also need the Gaussian part in operator form, which is given bŷ In the following we will evaluate these expressions for the general perturbed state (D.1) in terms of the perturbation function f (x, y).