An Inverse Mass Expansion for Entanglement Entropy in Free Massive Scalar Field Theory

We extend the entanglement entropy calculation performed in the seminal paper by Srednicki for free real massive scalar field theories in 1+1, 2+1 and 3+1 dimensions. We show that the inverse of the scalar field mass can be used as an expansion parameter for a perturbative calculation of the entanglement entropy. We perform the calculation for the ground state of the system and for a spherical entangling surface at third order in this expansion. The calculated entanglement entropy contains a leading area law term, as well as subleading terms that depend on the regularization scheme, as expected. Universal terms are non-perturbative effects in this approach. Interestingly, this perturbative expansion can be used to approximate the coefficient of the area law term, even in the case of a massless scalar field in 2+1 and 3+1 dimensions. The presented method provides the spectrum of the reduced density matrix as an intermediate result, which is an important advantage in comparison to the replica trick approach. Our perturbative expansion underlines the relation between the area law and the locality of the underlying field theory.


Introduction
Quantum entanglement is the physical phenomenon that appears when a composite quantum system lies in a state such that no description of the state of its subsystems is available. In the presence of quantum entanglement, measurements in the entangled subsystems are correlated. The most well known example of an entangled system, the so called EPR paradox [2], requires just two spinors; it was initially conceived as contradictory to causality, and, thus, as an adequate theoretical experiment to question the completeness of the quantum description of nature. However, later on, the corresponding correlations were verified experimentally.
A quantum subsystem A entangled to its environment A C cannot be described by a state; it is rather described by a density matrix ρ A , calculable by tracing out the degrees of freedom of the subsystem A C from the overall density matrix ρ ρ A = Tr A c ρ. (1.1) In the absence of entanglement, there is a state description for the subsystem A, and, thus, this reduced density matrix ρ A corresponds to a pure state; on the contrary, in the case entanglement is present, the reduced density matrix corresponds to a mixed state.
The above indicate that the entanglement is encoded in the spectrum of the reduced density matrix ρ A . It follows that a natural choice for a measure of entanglement is Shannon entropy applied to the spectrum of ρ A , known as Entanglement Entropy, S EE , S EE := −Tr (ρ A ln ρ A ) .
In a seminal paper [1], Srednicki performed a numerical calculation of entanglement entropy for a real free massless scalar field theory at its ground state, considering as subsystem A the degrees of freedom inside a sphere of radius R. The surprising result shows that entanglement entropy is not proportional to the volume of the sphere, but rather to its area. This profound similarity to the black hole entropy [25][26][27], discussed even before Srednicki's calculation [28], became even more intriguing after the development of the holographic dualities [29][30][31] and the Ryu-Takayanagi conjecture [16,17], which interrelates entanglement entropy in the boundary conformal field theory to the geometry of the bulk. The latter may allow the perspective of understanding the black hole entropy as entanglement entropy, and the gravitational interactions as an entropic force associated with quantum entanglement statistics [32][33][34][35].
In this context, the further investigation of the similarities between gravitational and quantum entanglement physics and the development of appropriate tools for their study presents a certain interest. In this work, we extend the original entanglement entropy calculation presented in [1] to massive free scalar field theory and develop a perturbative method for the calculation of entanglement entropy in such systems.
The majority of entanglement entropy calculations in field theory are based on the replica trick [14,[36][37][38][39][40][41]. This technique is based on the calculation of the entanglement Rényi entropies S q for an arbitrary positive integer index q > 1. (The entanglement Rényi entropy S q is defined as S q := 1 1−q ln Trρ q A .) Then, the entanglement entropy is recovered from the analytic continuation of S q as the limit S EE = lim q→1 S q . Although the entanglement Rényi entropies S q in principle contain the whole information of the reduced density matrix spectrum, the process of deriving the latter from the former is complicated. Relevant calculations are usually restricted to the specification of the largest eigenvalue and its degeneracy. Furthermore, holographic calculations cannot provide more information on the reduced density matrix spectrum other than the entanglement entropy, due to the very nature of the Ryu-Takayanagi conjecture [16,17]. An important feature of Srednicki's calculation is the fact that it is not limited to the calculation of entanglement entropy; on the contrary the full spectrum of ρ A is an intermediate result. As we discussed above, quantum entanglement is encoded into the spectrum of ρ A ; the entanglement entropy is just one piece of information. Therefore, although they are old, the methods of [1] present a certain advantage.
The structure of this paper is as follows: in Section 2, we review the derivation of entanglement entropy in systems of coupled harmonic oscillators lying at their ground state and extend the calculation in free scalar field theory including a mass term, closely following [1]. In section 3, we show that the inverse of the scalar field mass can be used as an expansion parameter allowing a perturbative calculation of entanglement entropy and develop the basic formulae of this perturbation theory. In section 4, we perform the perturbative calculation for massive free scalar field theory in 1 + 1, 2 + 1 and 3 + 1 dimensions and show that the leading contribution to the entanglement entropy for large entangling sphere radii obeys an area law; we specify the relevant coefficients and the first subleading corrections and we compare with numerical calculations. In Section 5, we discuss our results. Appendix A contains the details of Srednicki's regularization scheme. Appendix B contains the details of the perturbative calculation of entanglement entropy at second and third order. Finally, appendix C contains the code used for the numerical calculations of entanglement entropy.

Entanglement Entropy of Coupled Oscillators
The first step towards the calculation of entanglement entropy in free scalar field theory is the calculation of the latter in a system of coupled harmonic oscillators at its ground state. This problem has been solved long ago; here, we briefly sketch its solution. More details are provided in [1].
Assume a system of N coupled harmonic oscillators described by the quadratic Hamiltonian where the matrix K is symmetric and has positive eigenvalues, as required for the vacuum stability. Since K has been positively defined, its square root Ω := √ K can be appropriately defined, so that it also has positive eigenvalues.
In the following, without loss of generality, the subsystem A is considered to comprise of N − n oscillators, those described by coordinates x i with i > n. It follows that its complementary subsystem A C comprises of the n oscillators described by coordinates x i with i ≤ n. We may write the matrix Ω in block form as where A is an n × n, C is an (N − n) × (N − n) and B is an n × (N − n) matrix. We define the (N − n) × (N − n) matrices β and γ as, Let λ i , i = n + 1, . . . , N , be the eigenvalues of the matrix γ −1 β. Then, the spectrum of the reduced density matrix ρ A is given by where It follows that the entanglement entropy is given by The ground state of a system of coupled harmonic oscillators is a highly entangled state. The specification of entanglement entropy at the ground state requires a non-trivial, non-perturbative calculation. However, there is a small window allowing a perturbative approach. By definition, the matrix B does not contain any of the diagonal elements of the matrix Ω = √ K. Therefore, the matrix β, and, thus, the eigenvalues of γ −1 β, as well as the entanglement entropy, can be perturbatively calculated, in the case the diagonal elements of the matrix K are much larger than its non-diagonal elements. As the non-diagonal elements of K describe the couplings between the harmonic oscillators, in such an expansion, the zero-th order result is the entanglement entropy in a system of decoupled oscillators at their ground state, i.e. vanishing entanglement entropy.
The entanglement entropy is a valuable measure of entanglement, however, it does not contain the whole information. The latter is contained in the full spectrum of the reduced density matrix ρ A . An important advantage of the approach followed in this work is the fact that it allows the direct calculation of the latter through equation (2.7), as an intermediate step towards the calculation of entanglement entropy.

Entanglement Entropy in Free Scalar Field Theory
In the approach of [1], the degrees of freedom of the scalar field theory are discretized via the introduction of a lattice of spherical shells, and, thus, the introduction of a UV cutoff. Furthermore, an IR cutoff is imposed, putting the system in a spherical box. This inhomogeneous distretization may appear disadvantageous, as it breaks some of the symmetries of the theory; although it preserves rotations, it breaks boosts and translations. However, the consideration of the stationary entangling sphere, which separates the degrees of freedom to two subsystems, has already broken these symmetries. This approach reduces the problem of the calculation of entanglement entropy in field theory to a similar quantum mechanics problem with finite degrees of freedom. Since we are studying free scalar field theory, the latter quantum mechanical system is simply a system of coupled oscillators with a quadratic Hamiltonian at its ground state. More details on this discretazation scheme are provided in appendix A.

+ 1 Dimensions
Let us consider a free real scalar field theory in 3 + 1 dimensions. The Hamiltonian equals Decomposing the field to real spherical harmonics Y m , we find that the corresponding components ϕ m (r) obey canonical commutation relations of the form ϕ m (r) , π m r = iδ r − r δ δ mm , (2.11) where r = | x| is the radial coordinate. The only continuous variable left is the radial coordinate r. We regularize the theory introducing a lattice of N spherical shells with radii r i = ia with i ∈ N and 1 ≤ i ≤ N . The radial distance between consequent spherical shells introduces a UV cutoff 1/a, while the overall size of the lattice imposes an IR cutoff 1/(N a). The introduction of the spherical lattice sets the number of degrees of freedom for each pair ( , m) finite. The discretized Hamiltonian reads (2.12) Different and m indices do not mix and furthermore the index m does not appear explicitly in the Hamiltonian. It follows that the problem can be split to infinite independent sectors, identified by the index , each containing 2 + 1 identical subsectors.
We consider an entangling sphere of radius R = (n + 1/2) a. Then, the entanglement entropy at the ground state is given by (2.13) where S (N, n) is the entanglement entropy corresponding to the ground state of the Hamiltonian 14) The quadratic Hamiltonian (2.14) describes N harmonically coupled oscillators, and, thus, the problem of the calculation of S (N, n) has been reduced to the class of problems solved in section 2.1. For large , the Hamiltonian H becomes almost diagonal. Therefore, for large , the degrees of freedom are almost decoupled, and, thus, the system (2.14) at its ground state is almost disentangled. It can be shown that S (N, n) decreases with fast enough so that the series (2.13) is converging [1,42].

+ 1 Dimensions
In a similar manner, we may study free scalar field theory in 2 + 1 dimensions. The Hamiltonian reads We expand the field to real circular harmonics and then introduce a lattice of circular shells to find the discretized Hamiltonian Different indices do not mix. Therefore, in a similar manner to the problem at 3 + 1 dimensions, the problem can be split to infinite independent sectors, identified by the index . The entanglement entropy at the ground state is given by where S (N, n) is the entanglement entropy corresponding to the ground state of the Hamiltonian The calculation of the latter lies within the class of problems solved in section 2.1.

+ Dimensions
Finally, we consider a free real scalar field theory in 1+1 dimensions. The Hamiltonian reads We may directly apply the same discretization scheme to obtain (2.20)

An Inverse Mass Expansion for Entanglement Entropy
The discretized Hamiltonians (2.14), (2.18) and (2.20) are describing a system of N coupled harmonic oscillators that falls within the class of systems studied in section 2.1. We may thus proceed to calculate the entanglement entropy following the scheme of this section.

An Inverse Mass Expansion
As an indicative example, in 3+1 dimensions, the K matrix describing the interactions between the harmonic oscillators can be directly read from equation (2.14), where i, j = 1, 2, . . . , N . As we have commented in section 2.1, a perturbation theory can be applied when the diagonal elements of the matrix K are much larger than the non-diagonal ones. This criterion clearly is satisfied at the limit of a very large mass µ. A similar approach is followed in [42] focusing in the behaviour of entanglemment entropy for large . It has to be pointed out that the actual expansion parameter is neither m nor , but the diagonal elements of K themselves. In all numbers of dimensions under study, the matrix K is of the form We define the quantities k i and l i so that The parameter ε is the expansion parameter of the perturbation theory that we are about to develop, which is obviously of order 1/µ. The expansion in ε is also a semiclassical expansion; recovering the fundamental constants in the dimensionless expansion parameter ε, the latter assumes the form / (µac). This is in line with the fact that the zeroth order entanglement entropy in this perturbative approach vanishes.
In order to calculate the desired entanglement entropy, we need to calculate the square root Ω of the matrix K, then the matrices β, γ and finally the eigenvalues of γ −1 β, perturbatively in ε. There is one important detail that has to be taken into account in these perturbative calculations. Since the lowest order elements of K are the diagonal ones, this is also going to be the case for its square root Ω. However, the matrix B, being an off-diagonal block of the matrix Ω, does not contain such elements. The lowest order elements of B are the first subleading elements that appear in Ω. As a result, preserving a certain order in perturbation theory requires the calculation of the square root of K at one order higher than the desired order. In the following, we will present the calculation at first non-vanishing order, therefore we will keep two non-vanishing terms in the expansion of Ω.
The square root of the matrix K, with two non-vanishing terms equals (3.25) The blocks A, B and C of the matrix Ω obviously equal It is noteworthy that the above formulae contain only odd powers of ε. Furthermore, the matrix B contains only order ε terms to this order, as it does not contain any diagonal elements of Ω. Had we desired to calculate the eigenvalues of the reduced density matrix with two non-vanishing terms in the ε expansion, we should have calculated Ω at ε 5 order.
From now on, we need to keep only one non-vanishing term in our expressions.
The matrices A −1 and C −1 equal The matrix γ −1 is identical to the matrix C −1 at this order. The matrix β has a single non-vanishing element at this order, namely, Finally, Obviously, the matrix γ −1 β has only one non-vanishing eigenvalue at this order, being equal to its sole non-vanishing element, Thus, the entanglement entropy at first non-vanishing order equals

The Expansion at Higher Orders
Continuing the expansion at higher orders, several patterns appear in the form of the expansions of the related matrices. More specifically, as long as the matrix Ω is considered: • Only odd powers of ε appear in the expansion of Ω.
• The leading term in any element in the k-diagonal is of order ε 2k−1 . Therefore, the matrix Ω calculated with n non-vanishing terms in the perturbation theory contains non-vanishing elements up to the (n − 1)-diagonal.
• Any subleading term in the elements of the matrix Ω is four orders higher than the previous one. Thus, an element in the k-diagonal is written as a series of the form We use the above notation with three indices for the coefficients of the expansion. The subscript denotes the line number if the element lies on the top triangle of the matrix or the column number if it lies in the bottom triangle, the superscript denotes the number of the diagonal, whereas the superscript in parentheses is the order of the term in the ε expansion. The matrices A −1 and C −1 follow the same pattern with an overall increase by 2 to all orders. Namely, The matrix γ is defined as γ = C − β. The expansion for γ −1 β can be acquired using the formula The form of the expansions of Ω, A −1 and C −1 imply that the expansion of the matrix γ −1 β, whose eigenvalues define the spectrum of the reduced density matrix, follows a similar pattern. In this case, the leading order element is the (1, 1) element, which is of order ε 4 . Every offset by a column or a row increases the order of the leading term by 2. Again subleading terms in any element are four orders higher than the previous one, A direct consequence of the above is the fact that the eigenvalues of γ −1 β are naively expected to have a given hierarchy. The largest eigenvalue is of order ε 4 , the second largest is of order ε 8 and so on. The calculation at the next to the leading order is analytically presented in appendix B. It turns out that the second largest eigenvalue vanishes at this order, whereas the largest eigenvalue receives corrections at order ε 8 . At third order the calculation is straightforward but more messy. The result is presented in the appendix only in the appropriate limit for the specification of the "area law" contribution to the entanglement entropy that we will discuss in next section. At this order, the largest eigenvalue receives another correction at order ε 12 , while one more non-vanishing eigenvalue emerges, with a leading contribution at the same order. As a general rule, a new non-vanishing eigenvalue emerges every second order in the perturbation theory. The corrections to the largest eigenvalue at a given order in the expansion have a more important effect to the entanglement entropy than the emergence of new eigenvalues at the same order.

The Leading "Area Law" Term
In section 3 we managed to acquire an expansive formula for entanglement entropy. In order to study the dependence of entanglement entropy on the size of the entangling sphere, we need to expand our results for large entangling sphere radii. We assume that the entangling sphere lies exactly in the middle between the n-th and (n + 1)-th site of the spherical lattice. We define so that R = n R a is the radius of the entangling sphere. In the following, we will study the expansion of entanglement entropy for large n R .

+ 1 Dimensions
In 3+1 dimensions, the entanglement entropy equals the sum of entanglement entropy from all sectors, as shown in equation (2.13). The summation of this series cannot be performed analytically. For this reason, we will use the Euler-Maclaurin formula to approximate the series by the integral The coefficients B k are the Bernoulli numbers defined so that B 1 = 1/2. We would like to expand the above integral for large R. This cannot be done directly, since n R appears in S in the form of the fraction ( + 1)/n 2 R and becomes arbitrarily large within the integration range. This problem may be bypassed performing the change of variables ( + 1)/n 2 R = y to find Now we may expand the integrand for large n R . It is also simple to show that the n 2 R term of entanglement entropy receives contributions only from the integral term of formula (4.42). Therefore, the large R behaviour of entanglement entropy is completely determined by equation (4.44). The matrix elements of K in 3 + 1 dimensions are given by The integral (4.44) can be performed explicitly. At third order in the inverse mass expansion (B.43), we find The leading contribution to entanglement entropy for large entangling sphere radii is proportional to the area of the entangling sphere. This is the celebrated "area law" term calculated at third order in the inverse mass expansion. Using expansive techniques, we managed to acquire an analytic expression for the coefficient connecting the entangling sphere area to entanglement entropy. It is noteworthy to mention that the expansions in the inverse mass and in the size of the entangling sphere are not coinciding; the leading term in the latter expansion, i.e. the area law term, receives corrections at all orders in the inverse mass expansion. In order to verify the validity of our expansion, we compare the perturbative results in the form of formula (4.47) to the numerical calculation of entanglement entropy presented in [1], for various values of the mass parameter and N = 60. The numerical calculation is performed with the use of Wolfram's Mathematica; the code is provided in Appendix C. It is shown in figure 1 that the formula (4.47) is more accurate for large values of the mass parameter, as expected. Furthermore, the entanglement entropy is a decreasing function of the scalar field mass [42]. The divergence of the numerical results from the expansive formula for entangling sphere radii close to N a is an effect induced by the IR cutoff that has been imposed since the theory has been defined in a finite size spherical box.
The numerical calculation requires the introduction of a cutoff in the values of . The convergence of the series (2.13) gets slower as the mass parameter increases. Thus, the perturbative expansion has an additional virtue; it provides a result for entanglement entropy in cases that the numerical calculation is more difficult.
The parameter of expansion in our approach is the ration between the off-diagonal and diagonal elements of the couplings matrix K. This is not exactly the inverse of 1st order 2nd order 3rd order numerical It follows that the perturbative method can be applied even in the massless limit. Of course in such a case, the perturbation series converges much more slowly, nevertheless, it turns out that it does converge to the numerical results, as shown in figure 1. In 3 + 1 dimensions, the coefficient connecting area with entanglement entropy in massless scalar field theory has been calculated numerically in [1] and improved in [43], found approximately equal to 0.295. Setting µ = 0 to the area law derived above, we find S EE 3 + 2 ln 8 32 + 167 + 492 ln 8 36864 + −11 + 2940 ln 8 491520 which is a good approximation to the numerical result and can be further improved continuing the perturbative expansion at higher orders.
In 2 + 1 dimensions, the entanglement entropy equals the sum of all sectors, as shown by equation (2.17). With the use of Euler-Maclaurin formula (4.42), it may be approximated by the integral (4.50) As in 3 + 1 dimensions, in order to find the asymptotic behaviour of this integral for large entangling circles, we perform the change of variables = yn R , Now we may expand the integrand for large n R . In 2 + 1 dimensions, the matrix elements of K are given by . In general the perturbative series converges more slowly than in 3 + 1 dimensions.
In the massless case, we yield
1st order 2nd order 3rd order numerical Especially in the massless case, the perturbative formulae fail completely to capture the logarithmic behaviour of entanglement entropy ( figure 3 top-left). Technically, this happens due to the structure of the couplings matrix K. In all cases this matrix is diagonally dominant, i.e. the sum of the absolute value of all non-diagonal elements does not exceed the diagonal one, in all rows and columns. As a result, the perturbative calculation of its square root and its inverse converges. Only in 1 + 1 dimensions and only in the massless case, the matrix saturates the diagonally dominant criterion. Not unexpectedly, the saturating case, lying between convergence and divergence, leads to a logarithmic dependence on the cutoff scale. However, this logarithmic dependence cannot be evident in a finite number of terms of the perturbation series. We will return to this kind of behaviour in the section 4.2 on the subleading contributions to entanglement entropy.
The area law is the leading contribution to the entanglement entropy for large entangling sphere radii in all number of dimensions. The reason for this fact can be attributed to the locality of the underlying field theory [44][45][46]. The locality is depicted to the fact that the matrix K contains interaction elements only in the subdiagonal and superdiagonal. As a result, no matter what is the size of the entangling sphere (the value of n), there is only one element of K connecting a degree of freedom of subsystem A to a degree of freedom of subsystem A C . This property is inherited to the leading corrections in matrix B, and, thus, to the eigenvalues of the reduced density matrix. Had the theory been non-local, the number of leading contributions to entanglement entropy, would be a complicated function of the entangling sphere radius in general, leading to large divergences from the area law. In a more geometric phrasing, the area law emerges from locality, since the pairs of strongly correlated degrees of freedom (i.e. neighbours) that have been separated by the entangling surface, are proportional to its area.

Beyond the Area Law
The area law term of entanglement entropy is the leading contribution to the entanglement entropy for large radii of the entangling sphere. There are also subleading terms, which can also be calculated in the inverse mass expansion that we developed in section 3.
In 2 + 1 and 1 + 1 dimensions, the subleading terms vanish as a → 0. We will not extend our analysis to these cases. In the case of 3 + 1 dimensions, the first subleading term is a constant. There are two contributions to this term. The first one comes from the integral term (4.44) and can be acquired by appropriate expansion of the integrand. The second contribution comes from the rest of the terms in Euler-Maclaurin formula (4.42). Taking into account that lim →∞ (2 + 1) S EE = 0, the equation (4.42) reads .

(4.59)
Since the parameter appears in S only in the form of the fraction ( + 1) /n 2 R , any action of the derivative operator on S results in a term two orders smaller in the n R expansion. This implies that apart from the S 0 /2 term, we have only one more contribution at n 0 R order, namely the k = 1 term, and specifically the part of latter where the derivative acts on the factor 2 + 1 and not on S . Bearing in mind that B 2 = 1/6, the contribution to the constant term by the discrete part of the Euler-Maclaurin formula is S 0 /3. (4.60) In order to compare the constant subleading term found above to the numerical calculation of entanglement entropy, we performed a linear fit to the outcome of the numerical calculation of the form S EE = c 2 n 2 R + c 0 , for various values of the parameter µ 2 a 2 . The perturbative formulae indeed approximate the numerical results well, as shown in figure 4, apart from the massless limit.  At finite order in the inverse mass expansion, the first subleading term is a constant, even in the massless case. The usual treatment of entanglement entropy in 3 + 1 dimensions in either conformal field theory or in theories with holographic duals through the Ryu-Takayanagi conjecture predicts an expansion for entanglement entropy of the form So, how is the absence of the logarithmic term in our expansion explained? In the case of massive scalar field theory, the answer is the existence of a fundamental scale in the theory, that of the mass, which naturally cutoffs the logarithmic term. As far as the massless limit in 3 + 1 dimensions is concerned, the reason is more complicated and related to the failure to capture the leading entanglement entropy contribution in the same limit in 1 + 1 dimensions. In a similar manner our perturbation theory is unable to capture the constant term in massless 2 + 1 field theory. From a technical point of view, we understand this failure of our perturbation theory as follows: The formulae used in our perturbation theory for the square root of matrix K, as well as the formulae for the inverse of matrices A and C, present some "edge effects" due to the fact that the matrices used in the inverse mass expansion are of finite size. This can be seen in the form of the factors 1 − δ i1 and 1 − δ i,N −1 in formula (B.20) or the factor 1 − δ in in formula (B.29). Such "edge effects" can be treated analytically for arbitrary N and n in our expansion, as long as the order of the expansion is kept lower than the dimension of the matrices. If this is not the case, these "edge effects" will propagate through the matrix and eventually will get reflected at the opposite ends of the matrices that generate them, resulting in spreading all over the matrix elements. This qualitative behaviour implies the following: • The reflections of these "edge effects" will lead to matrices Ω, A −1 , etc, that depend on all the elements of the matrix K. Therefore, at high orders in the perturbation theory, such reflections introduce contributions to the entanglement entropy that depend on the global characteristics of the entangling surface. Such "universal" terms cannot be captured at any finite order in our perturbation series. They are rather non-perturbative effects in this expansion. The logarithmic term in even number of spacetime dimensions [16,17,[47][48][49][50][51], as well as the constant term in odd number of dimensions [48,49] are known to be exactly this kind of universal terms, and, thus, our inability to capture them in the inverse mass expansion should not be considered surprising. Of course such effects are visible in the numerical calculations.
• The terms we capture in our perturbation series cannot sense the global properties of the region defined by the entangling surface. They have the property to depend on the local characteristics of the entangling surface. In a more technical language, this is depicted to the fact that the perturbative expressions for the elements of the matrices Ω, A −1 and C −1 depend on a finite number of the elements of matrix K. This is the reason our method is appropriate to capture the "area law", as well as subleading terms that scale with smaller powers of the entangling sphere radius. Therefore, our method is appropriate to study the dependence of such terms on geometric characteristics of the entangling surface, such as curvature [50], if more general entangling surfaces are considered.
• The introduction of a mass exponentially dumps the propagation of the "edge effects" through the matrix elements [52]. As a result, our expansive calculations accurately converge to the numerical calculations in this case.

Dependence on the Regularization Scheme
Finally, we would like to comment on the dependence of the area law term, as well as the subleading terms of entanglement entropy on the regularization scheme. In our analysis, we have applied a peculiar, inhomogeneous regularization. Namely, we have imposed a cutoff in the radial direction, but not in the angular directions. Thus, the measurables that we have calculated, are those measured by a peculiar observer who has access to radial excitations of the theory up to an energy scale 1/a and to arbitrary high energy azimuthal excitations. We could have applied a different more homogeneous regularization imposing an azimuthal cutoff by constraining the summation series in to a maximum value equal to max . Such a prescription would make our approach more similar to a traditional square lattice regularization. Notice however, that even in the square lattice case, the imposed cutoff is a cutoff to each of the momentum components and not strictly an energy cutoff that would allow direct comparison with formulae like (4.61).
As we discussed above, locality enforces the area law term to depend on the characteristics of the underlying theory in the region of the entangling surface. Therefore, a natural selection for an azimuthal cutoff max , when considering a d-dimensional entangling surface should have the following property: the total number of harmonics with ≤ max should equal the area of the entangling surface divided by a d . In 3 + 1 dimensions this argument implies that a natural selection for the azimuthal cutoff is πR/a, whereas in 2 + 1 dimensions it implies max = πR/a. In all number of dimensions such a cutoff is of the form max = cR/a, where c is a constant. It is not difficult to repeat our analysis including this azimuthal cutoff. The only extra necessary steps are the introduction of a finite upper bound in the definite integral (4.44) and similarly the inclusion of the terms calculated at x = max in the Euler-Maclaurin formula (4.42).
As an indicative example, in 3 + 1 dimensions, the area law term calculated at second order in the inverse mass expansion assumes the form S EE = 3 + 2 ln 4 µ 2 a 2 + 2 (µ 2 a 2 + 2) − 3 + 2 ln 4 µ 2 a 2 + 2 + c 2 (µ 2 a 2 + 2 + c 2 ) 167 + 492 ln 4 µ 2 a 2 + 2 4608(µ 2 a 2 + 2) 3 − 167 + 492 ln 4 µ 2 a 2 + 2 + c 2 (4.62) This equation implies that the coefficient of the area law term depends on the regularization scheme. The coefficients calculated in section 4.1, which correspond to the selection c → ∞, serve as an upper bound for the area law coefficient. In order to investigate whether the inverse mass expansion is still a good approx-imation when an azimuthal cutoff of the form max = cR/a is introduced, the entanglement entropy in 3 + 1 dimensions for µa = 1 and various values of c is numerically computed and compared to the perturbative formulae (4.62) in figure 5. We may conclude the following: • An azimuthal cutoff of the form max = cR/a preserves the dominance of the area law term in entanglement entropy. This is not the case when a more general azimuthal cutoff is chosen (e.g. max = c). The inverse mass expansion is still a good approximation when such a regularization scheme is chosen.
• The area law term, as well as the subleading terms, are strongly affected by the selection of the dependence of the azimuthal cutoff max on the radial cutoff a. This is the expected behaviour comparing with calculations in CFT or holographic calculations via the Ryu-Takayanagi conjecture. The only terms that do not depend on the regularization scheme are the universal terms, which cannot be captured by our perturbation theory.
• The introduction of an azimuthal cutoff would also set the perturbative calculation of the entanglement entropy finite at higher number of dimensions, where the respective integral term diverges as max → ∞.
• Srednicki's calculation, which is equivalent to the specific choice c → ∞, is an upper bound for the area law coefficient. The fact that the integral terms in more than 3 + 1 dimensions diverge, implies that such an upper bound exists only in 2 + 1 and 3 + 1 dimensions.

Discussion
The calculation of entanglement entropy in the ground state of oscillatory systems, which include free scalar field theories, at their ground state is in general a difficult, non-perturbative calculation, since the ground state is highly entangled. We managed to find a perturbative method to calculate it, using as expansive parameter the ratio of the non-diagonal to diagonal elements of the couplings matrix of the system. This parameter in the case of free scalar field theory is being played by the inverse mass of the field. The calculation of entanglement entropy in the inverse mass expansion indicates that the major contribution to entanglement entropy is a term proportional to the area of the entangling surface, i.e. the "area law" term, a well-known fact since [1,28]. The perturbative calculation of the coefficient of this term agrees with the numerical calculation of entanglement entropy, based on the techniques of [1], and provides an analytic method for the specification of such coefficients. Subleading terms in the expansion of entanglement entropy for large entangling sphere radii can also be perturbatively calculated. The inverse mass expansion and the entangling sphere radius expansions can be performed simultaneously, but they are not parallel in any sense. The leading term in the entangling sphere radius expansion, i.e. the area law term, as well as the subleading terms, receive contributions at all orders in the inverse mass expansion.
The area law term, as well as the subleading ones are dependent on the regularization scheme, in line with analogous replica trick calculations. Universal terms that appear in the massless limit and depend on the global characteristics of the entangling surface (logarithmic terms in even dimensions and constant terms in odd dimensions) are non-perturbative contributions in this expansive approach. Furthermore, in this approach, the coefficient of the area law term in 2+1 and 3+1 dimensions has an upper bound, for any regularization scheme. The latter does not exist in higher dimensions.
An interesting feature of the inverse mass expansion is the following: the perturbation parameter is not exactly the inverse mass, but rather the quantity 1/ µ 2 a 2 + 2, where a is the UV cutoff length scale imposed in the radial direction. This fact allows the application of the perturbation series even in the massless field case. Not surprisingly, the perturbation series converges more slowly than in the massive case; however, the values of the first terms strongly suggest that it still converges to the numerical results. In the case of free massless scalar field in 3 + 1 dimensions the inverse mass series for the coefficient of the area law term approaches the value 0.295 found in [1,43].
An important advantage of the presented perturbative method is that it is not limited to the calculation of entanglement entropy, but it provides the full spectrum of the reduced density matrix. The latter, unlike entanglement entropy, contains the full information of the entanglement between the considered subsystems. This is clearly an advantage in comparison to holographic calculations which are constrained to the specification of entanglement entropy due to the nature of the Ryu-Takayanagi conjecture (the latter of course can be applied to strongly coupled system, where it is impossible to apply our perturbative method). As long as calculations based on the replica trick are concerned, they naturally allow the calculation of Rényi entropies S q for all q. Although in principle it is possible to reconstruct the spectrum of the reduced density matrix from the latter, in practise this process is very complicated and usually only the specification of the largest eigenvalue and its degeneracy may be easily achieved.
This perturbative method is an appropriate tool to expose the connection between the "area law" and the locality of the underlying field theory. Locality is encoded into the couplings matrix K as the absence of non-diagonal elements apart from the elements of the superdiagonal and subdiagonal. This in turn results in an hierarchy for the eigenvalues of the reduced density matrix system, leading to the area law. This hierarchy in the spectrum of the reduced density matrix depicts the fact that locality enforces entanglement between the interior and the exterior of the sphere to be dominated by the entanglement between pairs of neighbouring degrees of freedom that are separated by the entangling surface. The latter are clearly proportional to the area and not the volume of the entangling sphere.
It would be interesting to extend the applications of this perturbative expansion to other geometries, e.g. dS or AdS spacetimes, to cases where the overall system does not lie at its ground state (e.g. systems at a thermal state) or to other field theories containing fermionic fields or gauge fields. Furthermore, application of the above techniques for non-spherical entangling surfaces may shed light to the dependence of entanglement entropy on the geometric features of the latter, such as the curvature.

+ 1 Dimensions
We consider a free real scalar field theory in 3 + 1 dimensions. The Hamiltonian reads We define, where r = | x| is the radial coordinate and Y m are the real spherical harmonics, defined as, The real spherical harmonics form an orthonormal basis of harmonic functions on the sphere S 2 . It is easy to show that quantities ϕ m (r) and π m (r) obey canonical commutation relations, ϕ m (r) , π m r = iδ r − r δ δ mm . (A.5) Expanding the field in real spherical harmonics and substituting in (A.1), we acquire an expression of the Hamiltonian in terms of ϕ m (r) and π m (r), The only continuous variable left is the radial coordinate r. We regularize the theory introducing a lattice of N spherical shells with radii r i = ia with i ∈ N and 1 ≤ i ≤ N . The Hamiltonian of the discretized system can be found via the application of the following rules on equation (A.6): . (A.7) The discretized Hamiltonian reads (A.8)

+ 1 Dimensions
We may study free scalar field theory in 2 + 1 dimensions in a similar manner. The Hamiltonian reads We define, where r is the radial coordinate and Y are the real circular harmonics, The functions Y form an orthonormal basis of harmonic functions on the circle S 1 . The quantities ϕ (r) and π (r) obey canonical commutation relations, ϕ (r) , π r = iδ r − r δ . (A.13) We expand the field in real circular harmonics and substitute in (A.9) to find H = 1 2 ∞ 0 dr π 2 (r) + r ∂ ∂r ϕ (r) √ r 2 + 2 r 2 + µ 2 ϕ 2 (r) . (A.14) Using the discretization scheme (A.7), we obtain the discretized Hamiltonian

B The Inverse Mass Expansion at Second and Third Order
It is not difficult to show that there are no corrections to the entanglement entropy at the next to leading order in ε. Therefore, the first corrections appear at third order in the inverse mass expansion. For this purpose it is required that the matrix Ω is calculated with four non-vanishing terms. Following the definitions (3.22), (3.23) and (3.24), we find that the matrix Ω at order ε 5 is given by Trivially, The matrix B has only a finite set of elements not vanishing at this order, namely, We need to acquire the matrices A −1 and C −1 with three non-vanishing terms. They equal .

(B.40)
At third non-vanishing order in the inverse mass expansion, another non-vanishing eigenvalue emerges having a leading contribution of order ε 12 . Considering only the area law term contribution to entanglement entropy is equivalent to approximating all k i and l i with k n R and −1 respectively. At this approximation and without showing more details, the eigenvalues of γ −1 β at third order in the inverse mass expansion equal In section 4, the perturbative formulae for entanglement entropy were compared with a numerical calculation of the latter. The numerical algorithm uses the same regularization scheme as the perturbative expansion, but it calculates the matrices β, γ, as well as the eigenvalues of the matrix γ −1 β, numerically. The numerical calculation was performed with the help of the following code in Wolfram's Mathematica.  The above code applies in the case of 3 + 1 dimensions. For the numerical calculation it is necessary to impose a cutoff max to the values of . This has been selected appropriately large so that the series has converged close enough to the max → ∞ limit (such an appropriate choice of max depends on the value of the mass parameter). Specification of terms beyond the area law term requires much larger values of max , and, thus, running time. The required modifications for the numerical calculation of entanglement entropy in different number of dimensions or the introduction of an angular cutoff that depends on the entangling sphere radius are quite simple.