Universal terms of the entanglement entropy in a static closed universe

Subdominant contributions to the entanglement entropy of quantum fields include logarithmic corrections to the area law characterized by universal coefficients that are independent of the ultraviolet regulator and capture detailed information on the geometry around the entangling surface. We determine two universal coefficients of the entanglement entropy for a massive scalar field in a static closed universe $\mathbb{R} \times \mathbb{S}^3$ perturbatively and verify the results numerically. The first coefficient describes a well known generic correction to the area law independent of the geometry of the entangling surface and background. The second coefficient describes a curvature-dependent universal term with a nontrivial dependence on the intrinsic and extrinsic geometries of the entangling surface and curvature of the background. The numerical calculations confirm the analytical results to a high accuracy. The first and second universal coefficients are determined numerically with a relative error with respect to the analytical values of the orders $10^{-4}$ and $10^{-2}$, respectively.


I. INTRODUCTION
Correlations of vacuum fluctuations of a quantum field in a curved background carry information on the geometry of the background. A key example of the interrelation between the entanglement of field fluctuations and geometry is the celebrated area law [1,2], which states that the dominant contribution to the entanglement entropy of the vacuum for any finite region A of space is proportional to the area A of its boundary Σ = ∂A, where ε is an ultraviolet regulator and d is the dimensionality of spacetime. This geometric entropy was originally proposed as a source of entropy for black holes in [1,2]. Since then, several works have been devoted to clarifying its role in the quantum mechanics of spacetime, from the formulation of the holographic principle [3] and the description of geometric entropies through the AdS/CFT correspondence [4] to the the analysis of emergent properties of spacetime from quantum information-theoretic properties of quantum fields [5][6][7][8][9]. The fact that properties of the background geometry are imprinted in the field correlations allows one to extract information on the geometry from the network of correlations among the entangled spatial subsystems. In this regard, the area law (1) suffers from the drawback * rsoldati@usp.br † yokomizo@fisica.ufmg.br that, in general, it depends on the choice of regularization. However, subdominant terms in the entropy formula (1) include logarithmic divergences that are expect to be regularization independent [10]. These can be used to reliably extract geometric information from the entropy function, and provide a natural tool to study the relation between entanglement and geometry. Universal contributions for the entropy of massive fields were first obtained in [11], where the case of a flat waveguide background geometry was considered. For a free massive scalar field and an even number of spacetime dimensions d, a contribution of the form α 1 (d) m d−2 log(mε)A (2) was identified, with a universal coefficient α 1 (d) that depends only on the dimension d. In four dimension, α 1 (d = 4) = 1/(24π). In [10], a curved spherical waveguide background was considered, leading to the identification of new curvature-dependent universal terms. For an even d ≥ 4 and a spherical entangling surface, a universal term was found of the form where α 2 (d) depends on the Ricci scalar of the intrinsic metric on the entangling surface Σ and on the coupling constant ξ that describes the strength of the interaction of the field with the scalar curvature of the background. This universal term describes how the the intrinsic curvature at (d − 2)-dimensional surfaces Σ manifests itself in the entanglement of the field. In more general geometries, the universal coefficients α i are expected to capture more detailed information on the curvature around the entangling surface. Similarly to what happens for conformal field theories [12], one may expect the universal coefficients to depend on several scalars constructed from the the background curvature tensor and the intrinsic and extrinsic curvatures of the entangling surface, as discussed, for instance, in [10]. A strategy for computing universal terms of massive theories is provided by the perturbative approach introduced in [13,14]. This approach was applied to a scalar field in de Sitter space in [15], and indeed led to the discovery of curvature-dependent universal terms with a more complex dependence on the curvature tensors.
The analytical result for the universal term of a conformal field theory in four dimensions obtained in [12] was verified numerically for a massless scalar field and spherical entangling surfaces in [16], where the numerical approach originally explored for establishing the area law [2] was improved to allow for the determination of subleading corrections. The same universal term was derived by different methods in [17,18]. The analogous curvaturedependent universal terms for massive theories obtained in [15] have not been verified numerically, however, or confirmed by independent alternative derivations.
In this work, we determine universal terms of the entanglement entropy of a massive scalar field for spheres in the Einstein universe S 3 × R, both analytically, by the application of the perturbative approach introduced in [13,14], and numerically, through the application of the real-time approach originally introduced in [2] and discussed in the reviews [19,20]. The background geometry in the vicinity of a spherical entangling surface in the Einstein universe is a perturbed spherical waveguide, so that the results of [10] describe the zeroth-order term in the perturbations series. The lowest nonzero order terms describe new universal terms of the form (3), but with a universal coefficient α 2 (d) that depends nontrivially on several scalars constructed from the intrinsic, extrinsic and background curvatures at the entangling surface.
Being spatially finite, the Einstein universe has a natural infrared cutoff at the scale of its spatial radius, which is convenient for the numerical calculations. In addition, by considering spherical entangling surfaces of distinct radii, the intrinsic and extrinsic curvatures can be varied at the entangling surface, allowing the dependence of the universal coefficients in the distinct curvature terms to be analyzed. As the Einstein universe is static, there is no question as to the choice of the vacuum state, which is unique. Moreover, the discretization required for the numerical calculations can be implemented in a timeindependent manner. These properties single out the Einstein universe as a specially convenient geometric background for the study of curvature-dependent universal terms and, in particular, for a numerical test of the analytical techniques employed for their determination.
In Section II, we describe the relevant features of the theory of a massive neutral scalar field in the Einstein universe and its discretization. In Section III, we briefly review the perturbative approach for the calculation of the universal terms of the entanglement entropy and then apply it to the case of spheres in the Einstein universe. Next, we describe the techniques employed for the numerical calculation of these terms. The numerical results are presented in Section IV. We summarize and discuss our results in Section V.

II. THE MODEL
A. Massive scalar field in the Einstein universe The metric of the Einstein universe M = R × S 3 in spherical coordinates reads where dΩ 2 = dθ 2 + sin 2 θ dα 2 is the metric of the unit 2sphere, the coordinates are defined on the intervals t ∈ R, χ, θ ∈ [0, π] and α ∈ [0, 2π], and R is the constant radius of the 3-spheres describing spatial sections at fixed time. The volume element is The geodesic distance between antipodal points on the 3-spheres of constant time is Rπ. The area of a spherical surface of fixed χ is given by We consider a real massive scalar field Φ(x) on this background with a generic coupling to the scalar curvature. The action for the theory in the continuum is where the last term describes the interaction of the field with the scalar curvature 6/R 2 of the metric (4). The case of ξ = 0 describes the minimally coupled theory, while ξ = 1/6 corresponds to a conformal coupling. After an integration by parts in the angular coordinates θ, α, the Lagrangian assumes the form: where is the Laplace-Beltrami operator on the unit 2-sphere. The momentum associated with the field is the scalar density and the Hamiltonian is obtained as usual through the application of a Legendre transformation, The canonical fields satisfy the usual Poisson brackets, The analogue of the momentum representation in spacetimes with spherical spatial sections is obtained by expanding the field in real spherical harmonics, Integrating on the angular variables and using the orthogonality of the spherical harmonics, we obtain the following representation for the Lagrangian: The system was thus decomposed into a collection of independent fields Φ µ (χ) living on a one-dimensional space, as in [2]. In contrast with [2], however, where the background geometry is that of flat Minkowski spacetime, the one-dimensional space associated with the radial direction is now finite, reflecting the compactness of the spatial sections of the Einstein universe.

B. Discretization of the model
We now discretize the fields Φ µ (χ). The interval χ ∈ [0, π] can be partitioned into a union of N subintervals bounded by the equally spaced points At any fixed time, each χ j defines a 2-sphere with area The subregions bounded by these surfaces provide a decomposition of S 3 into a union of N − 2 thick spherical shells and two 3-balls (at the North and South Poles). The subspace formed by the n first subregions is bounded by a surface of area S n . Denote by ε the geodesic radial Fixing an angle θ, the discretization vertices χj span the range 0 ≤ χ ≤ π, covering a physical distance Rπ = N ε.
Each vertex is at a distance jε − N ε/2 from the equator.
distance between successive boundary surfaces. The maximum distance Rπ then becomes N ε (from the North to the South pole, see Fig. 1), and we have: In addition, the one-dimensional field Φ µ is replaced by a set of N variables that we interpret as living at the center of each subinterval of the decomposition, and its partial derivatives can be approximated by finite differences, We now approximate the Lagrangian (13) by a sum of contributions from each spatial subregion depending on the discretized field. Since the number of independent spatial derivatives is smaller than the number of field variables, terms with spatial derivatives must be handled differently from terms without derivatives in the fields. It is natural to require the discretization to preserve the symmetry under a spatial reflection about the equatorial surface χ = π/2. This can be achieved by discretizing the volume element differently for terms with and without spatial derivatives: By doing so, we discretize the integral over χ into a sum of N contributions for terms without spatial derivatives, and a sum of N −1 contributions for the term including spatial derivatives. In both cases the sum of the discretized volumes gives the total volume of the space: where the integral was replaced by a sum, The discretized Lagrangian then reads: The passage to the Hamiltonian formalism is straightforward. The conjugate momenta are readily obtained, and we find for the discretized Hamiltonian: Applying a canonical transformation that makes the transformed Φ and Π have the same dimensions, the Hamiltonian becomes with a potential term characterized by the potential matrix The discretized Hamiltonian is an explicit sum of decoupled normal modes µ. Each mode describes a set of coupled harmonic oscillators on a finite one-dimensional lattice. The potential matrix V ij , which arises from the spatial derivative terms in the original theory in the continuum, is responsible for coupling neighbouring vertices i and j. This is the source of entanglement between complementary regions of space, since V ij couples degrees of freedom inside any given region of space to degrees of freedom outside it, which become thus correlated.

A. Universal coefficients: perturbative calculation
In the continuum, the entanglement entropy of a quantum field is in general divergent, with a leading divergence proportional to the area A of the surface of the subregion, that is, the entropy satisfies an area law S ∝ A/ε 2 , where ε is an ultraviolet regulator [19]. The proportionality constant depends on the details of the regularization. Corrections to the leading divergence include universal contributions that are independent of the regularization, characterized by the couplings of the theory (see, for in- stance, [10,11]). The numerical constant appearing in a universal term is called a universal coefficient.
We are interested in computing universal contributions S univ to the entanglement entropy of spatial subregions in the Einstein universe. The spatial sections of constant time for the Einstein universe in the metric (4) are all isometric to a 3-sphere S 3 of radius R. We consider entangling surfaces Σ that are 2-spheres S 2 of area A embedded in a constant time slice (see Fig. 2). Any such entangling surface splits the spatial 3-sphere into the union of two complementary subregions, each isometric to a 3-ball.
We restrict to the case of a massive field of mass m, and consider the regime where correlations in the fluctuations of the field in the complementary subsystems are concentrated in the vicinity of the entangling surface. We thus assume that the correlation length m ∼ m −1 of the scalar field is much smaller than the radius of the entangling surface. The entanglement entropy must then be determined by the form of the metric near the entangling surface, which is a perturbation of a spherical waveguide geometry R 2 × S 2 .
In this regime, the perturbative approach introduced in [13,14] can be applied to the problem. This approach allows the entropy to be computed as that of the unperturbed metric together with corrections arising from its perturbations. The technique was applied to spherical entangling surfaces in de Sitter space in [15]. In this subsection, we review the relevant results of [13][14][15] and apply them to the case of spheres in the Einstein universe. The results so obtained will later be confronted with numerical calculations for the regularization of the unperturbed theory discussed in Section II.
We now proceed to the description of the perturbed spherical waveguide geometry around the entangling surface. Let us introduce a new variable r through: The variable r ∈ [−πR/2, πR/2] provides an arc length parametrization along the radial direction with origin at the equator, situated at χ eq = π/2. Following [15], we consider a Wick rotated metric with a temporal coordinate τ = it. The spacetime coordinates are denoted by x µ = (τ, r, θ, α), µ = 1, . . . , 4. Let g µν be the Riemannian metric obtained by Wick rotating the line element (4).
Expanding it around an entangling surface at r e , with r = r e + ∆r, we find The metric near the entangling surface, ∆r/R 1, has the form where the background geometryḡ describes a spherical waveguide of radius R cos(r e /R), and the perturbation h is given by up to second order in ∆r/R. As discussed in detail in [13], the perturbed metric (31) can be expressed in terms of the extrinsic and intrinsic curvatures of the entangling surface Σ. Following their approach, we relabel the coordinates parametrizing the entangling surface as y i = (θ, α) and the transverse coordinates as x a = (τ, ∆r), so that x µ = (x 1 , x 2 , y 1 , y 2 ). We can choose τ = 0 for the spatial section of interest, so that both transverse coordinates vanish at Σ. The unit normals to the entangling surface along the transverse directions have coordinates and the metric is such that g ia = 0 on Σ, as required in the formalism of [13]. We denote the intrinsic metric of the entangling surface by γ = g| Σ . The line element can then be written in the form: where the extrinsic curvature is defined as and R µνρσ is the curvature tensor of the Einstein universe evaluated at the surface, which has nonzero components only in the purely spatial part: (38) The term R acbd x c x d present in the transverse part of the metric in Eq. (36) in fact vanishes in the Einstein universe, but we keep it explicitly in the perturbation of the metric in order to obtain results that are valid in a more general class of spacetimes that include, in particular, spheres in the Einstein universe and de Sitter spacetime. For the extrinsic curvature, we find whereγ ij is the metric of the unit two-sphere. Substituting Eqs. (38) and (39) and is nonvanishing only tangentially to the surface. It includes contributions from the curvature tensor of the background and from the extrinsic curvature of the entangling surface.
As the starting point for the calculation of the entropy, we consider the unperturbed waveguide geometry. In [11], the heat kernel method was used to compute universal mass term contributions to S for a scalar field on a waveguide geometry via the replica trick, for a halfspace decomposition with a flat entangling surface. A universal logarithmic term of the form m 2 log(mε) was identified in (3 + 1)-dimensions. On a spherical waveguide, in addition to the universal term identified in [11], new universal contributions that are sensitive to the curvature of the entangling surface and background were obtained in [10]. We quote their result for the universal terms in where we set the radius of the spherical waveguide equal to R cos(r e /R). The first term, which is independent of the curvature, agrees with the universal contribution found in [11]. A different expression for the entropy (41) was presented in the Appendix B of [15], however, which follows from a different treatment of the heat kernel, where an extra overall factor of (1 − 6ξ) appears. In particular, in a conformally coupled theory this universal contribution would then vanish identically. We take Eq. (41) as our formula for the universal terms on the unperturbed geometry, as it will provide a better fit for the numerical results presented later in the paper.
The contribution δS of first order in the metric perturbation h can be calculated using the perturbative approach introduced in [13,14], as done on a de Sitter background in the Appendix B of [15]. As discussed in [15], such first order contributions have the form: where T µν is the energy-momentum tensor defined as The integrals in Eq. (42) refer to the transverse and longitudinal directions of the entangling surface in the spherical waveguide background. The average value T µν K must be computed in the spherical waveguide background.
Since we are working in the regime where the correlation length of the field is small compared to the radius of the entangling surface, T µν K can be approximated by that in a flat geometry, with small curvature corrections. As this quantity is multiplied by the first order perturbation h, we can safely neglect the curvature corrections and use the value of T µν K in flat space, computed in [21], as done for the case of de Sitter space in [15]. From [22], we can write the relevant components of such averages, expressed in our coordinate system, as where the spectral functions are given by The function K 0 is the modified Bessel function of order zero. In our case, the dimensionality of spacetime is d = 4.
Substituting Eq. (44) and the expression (40) for the metric perturbation in the formula (42) for the entropy, and using the identities we obtain an explicit formula for the contributions of the background and extrinsic curvatures for the variation of the entanglement entropy: A similar expression was obtained in [15] for the equatorial surface at the spatial section of minimum radius in de Sitter space. In that case, there are no contributions from the extrinsic curvature, which vanishes at the equator. We further discuss the relation between the results for de Sitter spacetime and the Einstein universe in the Appendix A.
The expression (45) can be regularized with the introduction of a hard cutoff at µ = 1/ε ≡ δ that eliminates ultraviolet contributions from length scales smaller than δ. The cutoff can then be sent to infinity. The leading divergence comes from the integration of the spectral function c (2) (µ), which gives: Substituting this leading divergence in (45), integrating over the two-sphere and setting d = 4, we obtain the universal contribution We are left with the task of computing the contractions of the curvature tensors, and finally obtain The first term is a contribution that depends on the background curvature tensor at the entangling surface. The second term is the contribution from the extrinsic curvature. It depends explicitly on the radius r e of the entangling surface and vanishes at the equator. The total entropy is the sum of the entropy of the spherical waveguide, given in Eq. (41), and the contribution from the metric perturbation. Our final result for the universal terms in the entanglement entropy of a scalar field for spheres in the Einstein universe is: with universal coefficients The first universal coefficient α 1 is independent of all parameters of the model. It describes a generic subleading logarithmic correction to the area law for the entanglement entropy that is independent of the coupling ξ to the scalar curvature and geometry of the entangling surface and background. It corresponds precisely to the universal term first identified in [11]. The second universal coefficient α 2 can be expressed in terms of the scalar curvature of the background and scalars constructed from the intrinsic and extrinsic curvatures of the entangling surface and that orthogonal to it, as follows.
At the entangling surface, the scalar of curvature (4) R of the perturbed background geometry is given by The tangential components of the curvature tensor are related to the components of the intrinsic and extrinsic curvature of the entangling surface by the Gauss-Codazzi identity [13]: where (2) R ijkl is the intrinsic curvature. An identical relation holds for the transverse surface Σ parametrized at each point of the entangling surface by the transverse coordinates x i . For a metric of the form (36), the extrinsic curvature of Σ vanishes, however, and the Gauss-Codazzi identity reduces to the simpler form: where (2) R abcd is the intrinsic curvature of Σ. Eqs. (51)-(53) allow us to express the formula (47) for δS univ exclusively in terms of scalars of curvature: with K a = γ ij K a ij , and where (2) R, (2) R are the Ricci scalars of the entangling and transverse surfaces, respectively, As the Gauss-Codazzi identity relates distinct curvature terms, the Eq. (54) can be equivalently expressed in terms of other sets of independent contractions of the curvature tensors.
The formula for δS univ in Eq. (54) is valid for any metric of the form (36). In the Einstein universe, , and we recover Eq. (48). The formula can also be directly applied to spheres in de Sitter space, in which case we recover the result of [15], as described in the Appendix A.
In addition, the contribution of the zeroth-order geometry to the universal terms given in Eq. (41) can be expressed in terms of the intrinsic curvature of the spherical entangling surface as The coefficient 1/(24π) corresponds to α 1 , while the second term within the brackets is a contribution to the second universal coefficient α 2 .
In the special cases of minimal and conformal coupling, the second universal coefficient reduces to (59) At the equator, r e = 0, the extrinsic curvature vanishes and the formula further simplifies to which also provides a good approximation near the equator, r e /R 1. In addition to the universal terms, the full entropy of the field includes regularization-dependent terms, both divergent, as the leading contribution to the area law, S ∝ A/ε 2 , and finite, as a whole tower of terms involving products of factors of the form (mε) 2p and (ε/R) 2q , with p, q ∈ N, multiplied by A, that show up in the integration of the spectral functions during the calculation of δS.

B. Entanglement entropy in the regularized theory
The entanglement entropy of spheres in the Einstein universe can also be computed numerically, exploring the discretization of the theory of a scalar field on this background discussed in Section II. The discretization introduces an ultraviolet regulator at the length scale set by the lattice spacing ε. The entropy of the ground state is then expected to include the universal logarithmic contributions described in Eq. (49). An area law term that scales with ε −2 should also be present [1,2,19], in addition to terms that remain finite in the limit of ε → 0.
For sufficiently small ε, terms that diverge in the limit ε → 0 will dominate. In this regime, the universal coefficients α 1 , α 2 can be determined from calculations of the entropy for different masses. In what follows, we will verify this numerically in order to corroborate our analytical results and, more generally, the perturbative approach developed in the works [13][14][15] and employed in our calculations. We describe the approach adopted for the calculation of the entropy in the discretized theory in this section, and discuss its numerical implementation in the next section.
The system of interest is the canonical quantization of the discretization of the scalar field on the Einstein universe introduced in Section II. The basic observables of the model are the canonical pairs {(Φ µj , Π µj )}. The Hamiltonian is given by Eqs. (27) and (29). It describes a set of coupled oscillators labelled by a multi-index a = ( , µ, j) over a one-dimensional lattice with nodes j = 1, . . . , N , with nearest-neighbor interactions described by the off-diagonal components of the potential matrix (29). The Hilbert space of the system is the tensor product where each H µj is the Hilbert space of an individual degree of freedom, i.e., the representation space for the canonical pair (Φ µj , Π µj ). The system naturally decomposes into a set of spatially localized subsystems, each associated with a single node j: Let N = {1, . . . , N } be the set of all nodes and N A ⊂ N a generic subset of nodes. The subsystem A associated with the set of nodes N A is described by the Hilbert space It consists of the representation space for the set of canonical pairs with j ∈ N A . The complement B of the subsystem is defined analogously with N A replaced with its complement N B = N \ N A . We thus obtain a bipartition H = H A ⊗ H B . The subsystem A describes the scalar field restricted to the spatial region formed by the union of the regions associated with nodes in N A , which are thick spherical shells, for i = 1, N , or a 3-ball at the South or North pole, for j = 1 and j = N .
In order to reproduce in the discrete theory the decomposition of S 3 into a union of two glued 3-balls, we consider a subsystem A formed by the first n nodes. From Eqs. (14) and (15), such a subsystem is bounded by an entangling surface of area The reduced density matrix is obtained by taking a partial trace over the last N − n nodes of the full state . The bipartition under consideration has thus the form: We wish to compute the entanglement entropy S of the subsystem A for the ground state of the Hamiltonian (27). Since the distinct angular momentum modes ( , µ) are decoupled, they constitute independent subsystems over the one-dimensional lattice with nodes N . The entropy of the subsystem A is then additive over the angular momentum modes, For a given ( , µ), the Hamiltonian of the associated mode has the form where the quadratic potential is independent of µ. As a result, modes with the same index and distinct index µ contribute equally, and we have We can then focus on the calculation of S µ . The ground state of a Hamiltonian that is quadratic in the canonical variables is a Gaussian state. General techniques for computing the entropy of Gaussian states are known, and can be directly applied to our problem. A general formula for the entanglement entropy was first derived in [23], reformulated in terms of symplectic invariants in [24] (see also the reviews [25][26][27]), and expressed in terms of the complex structure that characterizes Gaussian states in [28] (see the review [29]). The case of discretized free field theories has been considered in several works [19,20], and numerical results have been reported for varied fields and lattices [2,16,[30][31][32]. Of particular relevance to our purposes, an efficient algorithm for the numerical calculation of the vacuum entropy of a discrete free field is provided by the real time formalism reviewed in [19,20]. This technique was first applied to the study of entanglement of quantum fields already in the original works where the area law for the entropy was established [1,2]. Let us consider the application of the real time approach to the Hamiltonian (63).
A Gaussian state is completely characterized by the oneand two-point functions of the configuration variables Φ a and their conjugate momenta Π b , where a = ( , µ, i). The ground state of a quadratic Hamiltonian has vanishing onepoint functions, Φ a = Π a = 0. Hence, it is sufficient to specify its two-point functions. These are gathered in the covariance matrix C, which for vanishing one-point functions has the form: For a Hamiltonian without mixed terms consisting of products of a field and a momentum operator, the twopoint functions of the ground state satisfy the relations so that the off-diagonal blocks of the covariance matrix vanish. The relevant information is then encoded in the symmetric matrices which in our case take the following form: See Appendix B for details. Restricting the indices in the correlators to the subset A, we obtain the covariance matrix C A of the subsystem. Denote its diagonal blocks by X A and P A . All information on observations performed within the subsystem is encoded in the restriction of the covariance matrix to it. This allows the reduced density matrix to be reconstructed from C A . The entanglement entropy can then be computed from it and expressed directly in terms of the two-point functions. The entropy of the subsystem thus calculated has a simple expression in terms of the positive eigenvalues ν i of the matrix √ X A P A [19], The eigenvalues ν i satisfy ν i ≥ 1/2. This is the formula, together with Eq. (64), that we use for the numerical calculations of the entropy. We are thus provided with a convenient shortcut for the calculation of S. Instead of taking the partial trace of the full density matrix and using the formula for the von Neumann entropy, which involves traces in a Hilbert space of infinite dimension, it is sufficient to take the first n × n entries of the diagonal blocks X, P of the covariance matrix, which defines the matrices X A , P A , and compute the eigenvalues of √ X A P A . Such simple manipulations of linear algebra in finite dimensions can be implemented numerically in a straightforward way.

C. Numerical techniques
The discrete model under consideration has three characteristic length scales: the radius of the universe R, the lattice spacing ε and the inverse mass m −1 . The entanglement entropy also depends on the choice of the subsystem, which introduces a new length scale R e , the radius of the entangling surface, defined by A = 4πR 2 e . We need to fix these parameters for each evaluation of the entanglement entropy. This must be done for an interval of masses respecting the approximations involved in the perturbative calculation of the entropy, and in such a manner that the universal contributions are sufficiently large in comparison with the finite terms so that the universal coefficients can be accurately determined from the numerical results. We consider the case of conformal coupling, ξ = 1/6. The size N of the lattice determines the time required for an individual evaluation of the entropy. Let us first fix a choice of N . In the numerical calculations, an initial value of N will be progressively refined in order to give increasingly accurate estimates of the universal coefficients. A choice of N establishes a relation between two characteristic scales, as the radius of the universe and the lattice spacing are related through R = N ε/π, as described by Eq. (16). Hence, only one of these length scales is independent, say ε, for a fixed N .
For all calculations, we choose the equatorial surface as the entangling surface, in which case R e = R. In the lattice, this can be done by taking an even N , and letting the subsystem A be formed by the first N/2 nodes. The universal coefficients could be equivalently determined for any other entangling surface, but this choice allows us to avoid introducing an extra length scale in the problem.
We now observe that the entropy is invariant under a certain scaling of the parameters of the model. The existence of such a symmetry can be expected from the fact that the entropy is adimensional. From Eq. (67), for a transformation of the potential matrix of the form V → λV , the components of the covariance matrix transform as X → λ −1/2 X, P → λ 1/2 P , and the matrix √ X A P A whose eigenvalues determine the entropy remains unchanged. Therefore, the entropy is invariant under a transformation V → λV . Moreover, the factors of 2ε present in X and P as described in Eq. (67) cancel in the product √ X A P A and can be disconsidered for the calculation of the eigenvalues ν i . In the explicit formula (29) for the potential, apart from an overall factor of 1/ε, the potential matrix V depends on the length scales only through the combination εm. Therefore, the entropy is invariant under: We can then fix the lattice spacing arbitrarily. Results for any ε can be obtained through the application of the scaling transformation (69). In short, for any given N , we have: where ε can be fixed arbitrarily, and we consider the case of conformal coupling, ξ = 1/6. From Eqs. (61) and (70), the area of the entangling surface is From Eqs. (49), (60) and (70), the universal coefficients determined analytically in the perturbative approach are given for these parameters by: These are the analytical results that we wish to verify numerically. We will keep the product N ε fixed while refining the lattice to keep the coefficient α 2 constant under the variation of the lattice size N . We choose an arbitrary unit for lengths and express both ε and m −1 numerically in terms of this unspecified unit. The calculation of the entropy for a given mass follows three steps. First, for a given angular momentum mode , the potential matrix is diagonalized and the blocks X and P of the covariance matrix are computed using (67). Next, we take the restrictions X A and P A and compute the eigenvalues of √ X A P A . The entropy of the mode is given by Eq. (68). Finally, we sum over the modes , introducing a cutoff max , in order to obtain the total entropy. The cutoff max is increased until the numerical fit of the universal coefficients, which we will discuss next, stabilizes.
In addition to the universal terms described in Eq. (49), the entropy includes the non-universal area law term that diverges as ε −2 in the limit of ε → 0 and is proportional to the area of the entangling surface. Finite terms proportional to (mε) p A and (ε/R) q A, p, q ∈ 2N, are also expected, as discussed in Section III A. If R and ε are kept fixed, the latter finite terms describe a mass-independent contribution proportional to the area A, while the former finite terms become proportional to m p A. Accordingly, we model the dependence of the entropy on the mass, for fixed R and ε, with the function: (72) The coefficients α r , β s can be fitted for a given set of numerical evaluations of the entropy, {(S(m i ), m i )}. The dependence of the entropy on the unknown coefficients is linear. We estimate them using a multilinear regression based on the least squares method. The maximal power p max of the finite terms can be varied, allowing us to find an optimal choice that minimizes the uncertainties in the numerical estimates of the universal coefficients. If too few finite terms are included in the model, the universal terms absorb contributions of the finite part of the entropy in the fit, which affects the estimation of the universal coefficients, but including too many finite terms can lead to overfitting, making the linear regression more sensitive to numerical noise or systematic errors and thereby increasing the uncertainty in the results.
We are then left with the task of determining a suitable interval of masses to fit the coefficients of the function S(m) given in Eq. (72). In the analytical calculation of the universal coefficients, it is assumed that the correlation length of the field is small in comparison with the radius of the entangling surface, m ∼ m −1 R e = R. In the lattice, this is satisfied if In addition, the lattice must be sufficiently fine so as to provide a reliable approximation of the theory in the continuum. Accordingly, we require the lattice spacing to be small in comparison with the correlation length, An optimal value for the mass that takes into account both inequalities on the same footing is determined by the condition For sufficiently large N , one can expect to find an adequate interval of masses near m −1 = √ N ε that allow for a reliable numerical fit of the coefficients in the entropy function (72).
To determine such a suitable window of masses, we first computed the entropy for a large set of masses M scan = {m i } in an interval including the optimal mass m −1 = √ N ε, and then fitted the curve (72) for subsets {m i } imax i=imin ⊂ M scan , varying the position of the mass window, determined by the initial mass m imin , and the number i max − i min + 1 of masses in the windows. This allowed us to determine an interval of parameters where the fit is stable under the variation of the initial mass and the width of the mass window, for which the estimates of the universal coefficients do not change considerably.
The bulk code was written in FORTRAN 90 using the libraries LAPACK [33] and OpenBLAS [34] as linear algebra solvers, and the package OpenMPI as parallelization framework. For the lattice with N = 1500 sites, the time required for the calculation of the contribution of each angular momentum mode is of order ∼ 11 s/mode on an Intel i3 8100 processor. Our main results were obtained with a cutoff max = 5000, so that the time required for the calculation of the entropy S(m i ) of a single mass m i is approximately 15 h.

IV. NUMERICAL RESULTS
We report results obtained for lattices with N ≥ 1000 sites representing a universe of radius R = 1000/π. Preliminary results indicated that the estimation of the coefficient α 1 = 1/24π can be done accurately in coarser lattices, but not that of the numerically much smaller coefficient α 2 = 1/(360πR 2 ), which is strongly affected by errors in the estimation of α 1 . We observed that the coefficient α 1 must be determined with a relative error roughly at the order of 10 −4 in order that α 2 can be determined at the percent level, which required the lattices  to have at least N ∼ 1000 sites. We considered the case of conformal coupling and set the entangling surface at the equator, as discussed in Section III C. Consider a lattice with N = 1000 sites and lattice spacing ε = 1. From Eq. (71), the universal coefficients are then given by According to Eq. (75), the model (72) should provide a reliable fit of the entropy function S(m) for masses near m −1 = √ 1000 31.6. In order to determine an adequate number max of angular modes in the numerical calculations, we first chose a set of 48 equidistant masses M = {m i } in the interval m −1 ∈ (30, 50) and studied the variation of the universal coefficients under changes of max . We will discuss the choice of the interval of masses in more detail latter.
For a given mass, the contribution of a mode µ to the total entropy is given by Eq. (68), where ν i ≥ 1/2. The function that describes the contribution of each eigenvalue ν i to S µ satisfies lim νi→1/2 s(ν i ) = 0 , so that only for ν i > 1/2 we have nonvanishing contributions. Numerically, however, the function s(ν) is not well behaved at ν = 1/2; in addition, numerical noise can produce eigenvalues that are numerically less than 1/2. We first removed such contributions from the calculation of the entropy by introducing a cutoff tol and restricting the sum (68) to include only terms associated with eigenvalues such that ν i − 1/2 > tol. We decreased the value of the cutoff tol until no change was observed in the computed entropies S µ . We verified that this can be attained with tol = 10 −35 . We also checked that the total entropy S obtained by summing over the angular momentum modes as described in Eq. (64) remained unchanged under further decrease of tol for max = 10 4 . Next we computed the entropy S(m i ; max ) for several choices of max and estimated the universal coefficients α 1 , α 2 by fitting the curve (72) to the numerical data for each choice of max . This was done for small values of p max ∈ N. We observed that for max ∼ 5000, the relative variation ∆α 2 /α 2 in the estimated coefficient with the inclusion of higher modes reached the percent level. The results for max = 5000 and 10000 are compared in Table I. The relative variation of the coefficient α 1 is of the order ∼ 10 −6 for p max = 2, and rapidly increases with the inclusion of more finite terms in the fit, reaching ∼ 10 −2 for p max = 8. The relative variation of the coefficient α 2 is at the percent level for p max = 2, 4, and increases for larger p max , reaching approximately 70% for p max = 8. The computation time is proportional to max . In order to be able to considerably refine the lattice and increase the number of masses in the fit with the available computational resources, we set max = 5000 for the remaining computations. As a result, our estimation of the universal coefficient α 2 will be affected by systematic errors due to the cutoff in max which we estimate to be at the percent level for p max ≤ 4, and can only provide an order of magnitude estimate for larger values of p max . The uncertainties given in Table I are the statistical uncertainties in the multilinear regression used in the estimation of the universal coefficients. In these fits, the systematic errors are of the same order of magnitude as the statistical errors in the estimation of the coefficients (except for α 1 , p max = 2, when ∆α 1 /α (1) 1 is even an order of magnitude smaller than the statistical error).
Comparison with the analytical values (76) shows that a fit with a single finite term, p max = 2, is inconsistent with the theoretical predictions for N = 1000. Moreover, for p max = 4, the estimated values of α 1 and α 2 are within 4σ and 2σ from the theoretical predictions, including only statistical errors in the uncertainties. Adding a second finite term thus improves the accuracy of the fit, but the further inclusion of additional finite terms is not advantageous due to the increasing uncertainty in the estimation of the coefficients for fits with more variables, as discussed before. We fix then p max = 4 for our best estimates of the universal coefficients.
Another source of systematic errors is the lattice approximation of the continuum theory. These should decrease by refining the lattice. In finer lattices, the finite terms become less pronounced in comparison with the divergent contributions, and the accuracy of the estimation is expected to improve. With the same set of masses M, we computed the universal coefficients in lattices with a variable number of sites N = 1000, 1250, 1500 at fixed R = N ε/π = 1000/π. The results are described in Table II. We found that the coefficient α 1 does indeed gradually approach the analytical value given in Eq. (76). The coefficient α 2 reaches the analytical value, within its uncertainty, for N = 1250 and N = 1500. The estimation of the coefficients is thus stable under refinement of the lattice and approaches the analytical value for both universal coefficients. From the discussion in Section III C, for the lattice with N = 1500 sites and R = 1000, we expect the fit to be stable in some window of masses near m −1 25.8. In order to verify this, we selected a set M scan = {m i } of 256 masses ranging from m −1 = 2.45 to m −1 = 50.57, equally spaced in the axis m −1 , and fitted the universal coefficients for subsets {m i } imax i=imin ⊂ M scan , varying the position of the mass window, determined by the initial mass m imin , and the number of masses W = i max − i min + 1. For each mass window, the distribution of inverse masses m −1 i is centered at an inverse mass m −1 med , where m med is the median of the mass window.
In the first panel of Fig. 3, we plot the mean squared error χ 2 /W for fits with different numbers of masses W , for mass windows centered at m −1 med in the axis m −1 . The quantity χ 2 /W is used to evaluate the quality of the fit. We see that it stabilizes for larger m med , and reaches a common order of magnitude for all displayed W roughly around m −1 med ∼ 25. At m −1 med ∼ 30, the qualities of all displayed fits are comparable, except for the shorter window W = 20, which has larger fluctuations in χ 2 /W . At this m −1 med , the fits involve masses such that m −1 20. The quantity χ 2 /W stabilizes at larger m −1 med for fits involving a larger number of masses W . As the inverse masses are equidistant, a fit with larger W involves a larger width ∆m −1 of inverse masses than a fit with smaller W . The fits with W = 111 were performed on intervals of width ∆m −1 20. Keeping such a width ∆m −1 20 fixed, we verified that, for distinct numbers of points W ranging from 21 to 111, the quantity χ 2 /W stabilizes roughly at the same m −1 med for all W . Hence, the tendency observed in Fig. 3 of fits with larger W to stabilize at larger m −1 med can be attributed to the wider mass window employed in such fits.
We conclude that fits with W > 40 and masses satisfying m −1 20 over an interval of width ∆m −1 20 provide estimates of the universal coefficients with a stable quality. In the middle and right panels of Fig. 3, we show the estimated values of the universal coefficients α 1 and α 2 for mass windows with W = 110 centered at masses m med . We see that the fits are stable for m −1 med 25 for p max = 4 and that the values of the estimated universal coefficients agree with the analytical values.
In order to obtain our best estimate of the universal coefficients, we performed a fit using 110 masses such that m −1 > 30 for the lattice with N = 1500 sites, setting p max = 4. We obtained: errors. These are slightly lower than those for the estimates displayed in Table II for N = 1500, due to the increased number of masses used in the fit. The error in the estimation of α 2 is at the order of the estimated systematic error due to the cutoff max in the sum over angular momentum modes. Further increasing the number of masses in the fit might decrease the statistical error, but in order that the fit provides a more accurate estimation of the coefficients it would be necessary to also decrease the systematic errors, by increasing the number of angular momentum modes in the calculation of the entropy and perhaps further refining the lattice. For the purpose of numerically testing the predictions of the perturbative approach for the calculation of the universal coefficients, we consider that our best estimates already provide strong evidence for the correctness of the entropy formula (49) with universal coefficients given by Eq. (50), obtained with the application in Section III of the perturbative approach developed in [13,14]. The first universal coefficient α 1 was determined with a relative error of the order 10 −4 , and the second, curvature-dependent universal coefficient α 2 , which to the best of our knowledge has not been obtained numerically before, was determined up to a relative error at the percent level. The fitted entropy function is plotted against the numerical data in Fig. 4.

V. DISCUSSION
We have determined two logarithmic universal terms of the entanglement entropy for a massive scalar field in four spacetime dimensions, both analytically and numerically. The universal terms are characterized by numerical coefficients α 1 and α 2 , We considered the case of spherical entangling surfaces in the Einstein universe R × S 3 . The universal coefficients were first computed perturbatively following [13][14][15] and expressed in terms of geometric properties of the background and entangling surface. They were then estimated numerically with the application of the real-time approach [19,20] to a discretization of the theory. We observed a close agreement between the analytical and numerical results.
For the perturbative determination of the universal coefficients, we considered a general class of spherically symmetric spacetimes allowing for the variation of the intrinsic, extrinsic and background geometry at the spherical entangling surface. The geometry around the entangling surface is a perturbed spherical waveguide characterized by the background curvature R µνρσ and the extrinsic curvature K a ij evaluated at the entangling surface (Eq. (36)). Restricting to the case of de Sitter space, we recovered the universal terms described in [15]. We focused on the case of the Einstein universe in order to confront the general perturbative formula for the universal coefficients with numerical calculations.
The first universal coefficient α 1 is independent of all parameters of the model, α 1 = 1/(24π), describing a universal term of the entropy characteristic of massive scalar field theories in general geometries and for any coupling ξ to the scalar curvature. It corresponds to the universal term first obtained in [11]. We determined it numerically with a relative error of the order of 10 −4 . Such a high accuracy in the numerical estimation of α 1 is required in order that the second universal coefficient α 2 can also be estimated from the numerical data.
The coefficient α 2 includes a contribution from the unperturbed spherical waveguide geometry that depends on the coupling constant ξ to the background scalar curvature and on the intrinsic curvature of Σ, described by Eq. (58), and contributions from the metric perturbations that include terms proportional to the scalar curvatures (4) R, (2) R, (2) R of the background, entangling surface Σ and the surface Σ orthogonal to it, respectively, as well as terms proportional to the contractions K ail K ail , K a K a of the extrinsic curvature at Σ, as described in Eq. (54). The numerical value of α 2 at the equator of the Einstein universe is identical to that for an entangling surface of the same radius at the equator of de Sitter space, computed in [15]. For an entangling surface of generic radius, the universal coefficients in the Einstein universe are given by Eqs. (54) and (58), with the required contractions of curvature tensors given explicitly in Eq. (57). The universal coefficient α 2 was computed numerically at the equator of the Einstein universe. Our best estimate agrees with the analytical value up to a relative error of ∼ 3.5%. The close agreement between the numerical and analytical results for both universal coefficients provides a stringent numerical test of the perturbative approach to the calculation of the entanglement entropy of massive fields. For the numerical determination of the universal coefficients, the real-time approach was applied to a lattice model describing a discretization of the scalar field along the radial direction, in a straightforward adaptation of the approach introduced in the original numerical verification of the area law [2] to the case of the curved background of the Einstein universe and a massive field. The entropy was then computed for a sufficiently large set of masses {m i } and the universal coefficients were obtained by fitting a model for the entropy curve including the universal terms (78) to the numerical data S(m i ). Two sources of systematic errors are introduced in this approach: a required cutoff max in the set of angular momentum modes of the field and the lattice regularization itself. These can be reduced by increasing the cutoff max and decreasing the lattice spacing ε, at the cost of increasing the CPU time required for the calculation. We have progressively improved the numerical calculations in this manner until the estimation of both universal coefficients became reasonably stable. The mass window {m i } was chosen so that the approximations involved in the perturbative calculations were valid. A lattice with N = 1500 sites, a momentum angular cutoff max ∼ 5000 and at least ∼ 40 masses were required for the determination of the coefficient α 2 at the percent level. Our best estimate was obtained with a set of 110 masses and the total CPU time was approximately ∼ 1680 h.
The numerical approach explored here for the case of a massive scalar field in the Einstein universe can also be applied to the estimation of universal coefficients of the entanglement entropy in other spherically symmetric static geometries and to the case of massless theories. In particular, it should be possible to determine universal terms for a massless scalar field in the Einstein universe with a similar procedure. In the massive case, the entropy curve could also be further studied by extending the range of masses to regimes where the approximations assumed in the perturbative approach do not hold. The perturbative calculation of universal terms of the entropy for a spherical entangling surface in a perturbed spherical waveguide was discussed in Section III A. We considered the case of the Einstein universe in detail, but a general formula for the universal terms was first derived for any metric of the form (36). It includes a contribution from the unperturbed spherical waveguide, given by Eq. (41), and a contribution from the metric perturbation, expressed entirely in terms of scalars of curvature in Eq. (54). Here we apply this formula to the case of de Sitter space and show how the results obtained in [15] are reproduced from these formulas.
The metric of de Sitter space after a Wick rotation of the temporal variable is that of a four-sphere. It can be written in the form: where R is the radius of the four-sphere. Denote by y i the angular variables parametrizing the two-sphere with metric dΩ 2 . Introducing the new variables: x 1 = r cos τ , x 2 = r sin τ , and expanding around the equator (r = 0), we obtain the perturbed metric, ds 2 = δ ab dx a dx b + γ ij dy i dy j where γ ij is the metric of the two-sphere of radius R. The first line in Eq. (A1) corresponds to the unperturbed spherical waveguide, and the remaining terms describe the metric perturbation h. The curvature tensor is given in the full metric by: The extrinsic curvature vanishes at the equator, K a ij = 0 , at x a = 0 .
The perturbed metric (A1) is written in terms of the curvature tensor evaluated at the equator as: which has the form (36) with K a ij = 0. This allows us to apply the formula (54) to compute the contribution δS univ of the metric perturbations to the universal terms of the entanglement entropy. The contribution from the unperturbed spherical waveguide geometry is the same as before, and given by Eq. (41).
In de Sitter space, the relevant scalars of curvature are: (4) R = 12 R 2 , (2) R = (2) R = 2 R 2 , and we find δS univ = 1 360πR 2 , which corresponds to Eq. (B.7) of [15]. The result is numerically identical to that at the equator of the Einstein universe, but the curvature contributions from the metric perturbations to the result are distinct. At the equator of de Sitter space, the extrinsic curvature vanishes, and the transverse surface has a nonzero intrinsic scalar curvature (2) R. In the Einstein universe, there are contributions from the extrinsic curvature, while the curvature of the transverse surface vanishes. The contributions from the background are also distinct. In both cases, the entangling surface is a sphere of radius R, and the contribution from its intrinsic curvature (2) R is the same.