On the continuum limit of the entanglement Hamiltonian of a sphere for the free massless scalar field

We study the continuum limit of the entanglement Hamiltonian of a sphere for the massless scalar field in its ground state by employing the lattice model defined through the discretisation of the radial direction. In two and three spatial dimensions and for small values of the total angular momentum, we find numerical results in agreement with the corresponding ones derived from the entanglement Hamiltonian predicted by conformal field theory. When the mass parameter in the lattice model is large enough, the dominant contributions come from the on-site and the nearest-neighbour terms, whose weight functions are straight lines.


Introduction
The reduced density matrix of a subsystem is a central quantity to study in order to understand the entanglement properties of a quantum state for the spatial bipartition. Denoting by A a spatial subregion and byĀ its complement, under the assumption that the Hilbert space of the system factorises as H = H A ⊗ HĀ, the reduced density matrix is ρ A ≡ Tr HĀ ρ, where ρ is the density matrix of the entire system. Considering a system in its ground state, the entanglement entropy S A ≡ − Tr H A (ρ A log ρ A ) measures the bipartite entanglement in this pure state; hence it has been widely explored during the past two decades [1][2][3][4][5][6] (see [7][8][9] for reviews). The reduced density matrix can always be written as ρ A = e −K A /Z A , where the constant Z A = Tr H A (e −K A ) guarantees the normalisation condition Tr H A ρ A = 1 and the operator K A is the entanglement Hamiltonian (also known as modular Hamiltonian).
In certain relativistic quantum field theories (QFTs) and for some particular choices of states and geometric bipartitions, the entanglement Hamiltonian is given by the integral over A of the energy density multiplied by a suitable weight factor. The most important example has been found by Bisognano and Wichmann [10,11]: for a Lorentz invariant quantum field theory in the d + 1 dimensional Minkowski spacetime in its ground state, when A is the half space x > 0, the entanglement Hamiltonian K A is the boost generator in the x-direction.
In conformal field theories (CFTs), this result and the conformal invariance have been employed to obtain other entanglement Hamiltonians. A seminal example is the entanglement Hamiltonian of a sphere B of radius R for a d + 1 dimensional CFT in Minkowski spacetime and in its ground state. It reads [12,13] (see also [14]) where the weight function β(r) is the following parabola (1.2) When d = 1, specific conformal mappings have been constructed to obtain other entanglement Hamiltonians K A in the local form, i.e. written as an the integral over A of the energy density multiplied by the proper weight factor [15][16][17].
It is important to understand how these QFT results can be obtained as the continuum limit of the corresponding entanglement Hamiltonians in many-body quantum systems on the lattice. For free fermionic and bosonic systems, the Gaussian nature of the ground state allows to obtain explicit expressions for the entanglement Hamiltonian of a lattice subsystem for a generic number of spatial dimensions d [8,9,[18][19][20]. These lattice operators are characterised by long-range and inhomogeneous couplings [21][22][23]. In the special case of d = 1, the continuum limit of the entanglement Hamiltonian of a block made by consecutive sites in a chain of free fermions in the ground state has been studied analytically in [24] by employing the results of [22], finding the parabolic weight function (1.2) expected from CFT. In a massless harmonic chain, the corresponding analysis has been performed numerically in [21,25].
The entanglement Hamiltonians of a block of consecutive sites have been studied numerically also for non-critical chains in their ground state, finding also in these cases that the entanglement Hamiltonian matrices contain long-range and inhomogeneous couplings. Far away from criticality, a triangular profile for the weight function has been observed [26], which has been understood through the analytic expressions derived for the entanglement Hamiltonian of the half infinite chain [27,28].
In this manuscript we consider the entanglement Hamiltonian of a sphere B ∈ R d , mainly focussing on the CFT given by the massless scalar field in the d + 1 dimensional Minkowski spacetime and in its ground state. We study the continuum limit that leads to the entanglement Hamiltonian K B of the sphere given by (1.1) and (1.2) specialised to this model. Our analysis is mostly numerical and it is based on the method developed in [21,22,24,25] for d = 1. In the special case of d = 1, we recover the results for the entanglement Hamiltonian of a segment at the beginning of the semi-infinite line with Dirichlet boundary conditions obtained in [25]. In the massive regime, we adapt to the higher dimensional case of the sphere the analysis made in [26] for the entanglement Hamiltonian of the segment in the massive harmonic chain on the infinite line.
The layout of this paper is as follows. In Sec. 2 we introduce the model of the massive scalar field, the lattice regularisations of its Hamiltonian along the radial direction, the CFT prediction for the entanglement Hamiltonian of the sphere (1.1) for this model and the corresponding expressions on the lattice employed in our numerical analysis. In Sec. 3 we focus on the massless case and study numerically the continuum limit of the entanglement Hamiltonian. In Sec. 4 we discuss the regime where the mass parameter in the lattice model is large. Conclusions are drawn in Sec. 5. In the appendices A and B we report further results that clarify and support some discussions in the main text.

Entanglement Hamiltonian of a sphere for the scalar field
In this section we introduce the main expressions employed in this manuscript to study the entanglement Hamiltonian of a sphere in the d + 1 dimensional Minkowski spacetime for the scalar field in its ground state. In Sec. 2.1 we briefly review the Hamiltonian of the scalar field and its lattice regularisation along the radial direction, as done by Srednicki [2] to study the entanglement entropy of a sphere. In Sec. 2.2 we combine the result of this analysis with the expression for the entanglement Hamiltonian of a generic region in harmonic lattices found by Casini and Huerta [9].

Hamiltonian and radial regularisations
The Hamiltonian of the massive real scalar field in the d + 1 dimensional Minkowski spacetime reads where x ∈ R d is the vector identifying the spatial position, ∆ denotes the d-dimensional Laplacian, Φ = Φ(x) is the real scalar field and Π = Π(x) its canonically conjugate momentum field. In order to study the entanglement Hamiltonian K B of a sphere B, it is convenient to employ the hyperspherical polar coordinates of R d with origin in the center of the sphere. These coordinates are given by the radial coordinate r = |x| = x 2 1 + · · · + x 2 d 0 and by the vector Ω = (φ, θ 1 , . . . , θ d−2 ) collecting all the angular coordinates 1 of the d−1 dimensional unit sphere S.
In these hyperspherical polar coordinates, the metric reads ds 2 = − dt 2 + dr 2 + r 2 dΩ 2 and the corresponding volume element to employ in (2.1) is d d x = r d−1 dr dΩ, where dΩ is the volume element of S. The Laplacian operator ∆ reads ∆ = ∂ 2 r + d−1 r ∂ r + 1 r 2∆ , where∆ is the spherical Laplace operator in d − 1 dimensions, whose eigenfunctions are the real spherical harmonics in d dimensions Y l (Ω) of degree l (see e.g. chapter XI in [29] or [30,31]) where l = (l, i) with l 0 and 1 i N d,l labels the linearly independent spherical harmonics of degree l whose total number is which can be obtained from the degeneracy of the SO(d) representations.
1 Denoting by x1 the coordinate along a given vertical axes, the angular coordinates are given by in terms of the Cartesian coordinates of R d , whose ranges are φ ∈ [0, 2π) and θj ∈ [0, π].
The entanglement entropy of a sphere for the massless scalar field has been first studied by Srednicki in [2]. Following his analysis, we decompose the fields in (2.1) as where the sums can be written as l (. . In the special case of d = 1, the scalar field is on the half-line r 0 and the angular coordinates Ω do not occur; hence the decomposition in (2.5) becomes trivial. Instead, when either d = 2 or d = 3, where we have respectively N 2,l = 2 and N 3,l = 2l+1, the decomposition of Φ in (2.5) reads respectively (2.6) By employing the orthonormality condition S Y l 1 (Ω) Y l 2 (Ω) dΩ = δ l 1 ,l 2 for the spherical harmonics, one finds that the partial wave components Φ l (r) and Π l (r) in (2.5) are These fields obey the canonical commutation relations We remark that, for any l, the fields (2.7) satisfy Dirichlet boundary conditions at r = 0. This is a crucial feature in our analysis.
By employing (2.3) and the decompositions (2.5) in the Hamiltonian (2.1), it becomes From (2.9) and the fact that the fields (2.7) corresponding to different l's commute, one concludes that the ground state of H is the direct product of the ground states of H l for different l's [2]. When d = 1, only l = 0 is allowed; hence the coefficient (2.11) vanishes and the sum in (2.9) contains a single term. In this case (2.10) becomes the energy density for the scalar field on the half-line satisfying Dirichlet boundary conditions at the beginning of the half-line [32].
Following the regularisation procedure discussed in [2] (see also [1]), the ultraviolet divergences are regularised by introducing a discretisation of the radial direction (with lattice spacing a) and the infrared ones by confining the system into a finite volume; hence we consider a large but finite number R tot of sites along the discretised radial direction. At the j-th site of the discretised radial direction 2 , the position operatorq l,j and the momentum operatorp l,j are dimensionless, Hermitian and satisfy the canonical commutation relations q l 1 ,i ,p l 2 ,j = i δ l 1 ,l 2 δ i,j .
In the continuum limit a → 0 and R tot → ∞. Since we keep R tot a = R tot finite in this limit, the system in the continuum is enclosed within the finite sphere of radius R tot . The radial position corresponds to r = ja, with 0 r R tot , and the fields Φ l (r) and Π l (r) in the continuum limit (which vanish identically for r > R tot [2]) can be introduced through the position and momentum operators in the standard waŷ (2.12) The lattice regularisation defined above leads to (q l,j+1 +q l,j−1 −2q l,j )/a 2 for the regularisation of ∂ 2 r Φ l . Thus, for the lattice regularisation H l of the operator H l in (2.9) we find wherep l andq l are the vectors whose j-th element is given byp l,j andq l,j respectively. The matrix I tot is the R tot × R tot identity matrix and M l is the R tot × R tot real, symmetric, positive definite and tri-diagonal matrix whose non vanishing entries are where µ d (l) has been defined in (2.11), the addition term 2 in (M l ) j,j comes from the discretisation of Φ l ∂ 2 r Φ l and ω corresponds to the mass parameter in the lattice along the radial direction, which is related to the mass parameter m occurring in the Hamiltonian (2.1) in the continuum limit as ω/a −→ m. Since µ d (l) is quadratic in the parameters d and l and the coefficients of d 2 and l 2 are both positive, the large l regime at given d is similar to the regime given by large d with l fixed.
The ground state correlators q l,iql,j and p l,ipl,j provide the elements of the R tot × R tot correlation matrices Q l and P l respectively, which can be obtained from (2.14) as [1,2,[33][34][35] These correlation matrices are well defined when ω = 0, for any d 1 and l 0. In particular, since translation invariance does not occur in the lattice model along the radial direction, we do not have to deal with the zero mode. We checked numerically that, when µ d (l) = 0 in (2.14) and for large values of R tot , the elements of the correlation matrices (2.15) agree with the analytical expressions for the correlators on the half line with Dirichlet boundary condition at the origin of the half-line given in [36] for ω = 0 and in [25] for a generic ω 0 (see [37] for the corresponding correlators for the harmonic chain on the line, where the zero mode occurs). In the appendix A we compare (2.15) when ω = 0 with the corresponding expressions obtained in the continuum through quantum field theory methods [38,39].

Entanglement Hamiltonian of a sphere
We consider the spatial bipartition of R d given by a sphere B of radius R and its complement in the d + 1 dimensional Minkowski spacetime at any given time slice. When the system is a CFT in its ground state, the entanglement Hamiltonian K B is (1.1), which has been obtained in [12,13] through a conformal transformation of the entanglement Hamiltonian found by Bisognano and Wichmann [10,11].
In the special case of the free massless scalar field, K B is (1.1) with T tt given by the integrand of (2.1) when m = 0. In the resulting K B , we can straightforwardly adapt the steps that have led to write (2.1) as (2.9), finding that with β(r) being defined as the parabola (1.2) and the operator T tt by (2.10) with m = 0. In the regularised model obtained by discretising the radial direction and described in Sec. 2.1), the sphere B has radius R < R tot . In the continuum limit, a → 0 and R → +∞ with Ra = R fixed. The reduced correlation matrices Q l,B and P l,B at a given value of the angular momentum parameter l are the R × R symmetric and positive definite matrices obtained by restricting the corresponding correlation matrices in (2.15) to the sites identifying the sphere B along the radial direction, i.e. (Q l,B ) i,j = (Q l ) i,j and (P l,B ) i,j = (P l ) i,j with 1 i R and 1 j R. by l is given by the following quadratic operator [9] where H l,B is the symmetric, positive definite and block diagonal matrix defined in terms of the reduced correlation matrices Q l,B and P l,B at fixed l as follows The symplectic spectrum of H l,B in (2.19) gives the single particle entanglement energies ε l,k at given l (hence ε 2 l,k are the eigenvalues of V (l) T (l) ), which are related to the symplectic eigenvalues of Q l,B ⊕ P l,B as ε l,k = 2 arccoth(2σ l,k ) = log[(σ l,k + 1/2)/(σ l,k − 1/2)], whose inverse is σ l,k = 1 2 coth(ε l,k /2). Notice that both σ l,k and ε l,k depend also on d through (2.14).
3 Continuum limit of the entanglement Hamiltonian at ω = 0 In this section we study the entanglement Hamiltonian K B of a sphere B ∈ R d for the free massless scalar field by taking the continuum limit of the corresponding entanglement Hamiltonian K B in the regularised model, whose contribution for a given l is (2.18) at ω = 0. The expected result is the CFT expression (1.1) specialised to the free massless scalar field, which can be written as in (2.16), in terms of the operator (2.10) with m = 0. In the resulting CFT expression a non trivial r dependent term occurs whose coefficient (2.11) depends both on d and l. This coefficient vanishes identically when d = 1. In our numerical analysis we follow the procedure employed in [21,22,24,25] for free bosonic and fermionic chains (i.e. for d = 1) at criticality and in their ground states to study the entanglement Hamiltonians of an interval either on the infinite line [21,22,24,25] or at the beginning of the half-line with Dirichlet boundary conditions at the origin [25].
Plugging the matrix (2.19) into (2.18), one finds that the operator K l,B can be written as where the operators H T,l and H V,l are defined through the symmetric matrices T (l) and V (l) as follows These sums can be organised in various ways, as discussed in [25] for d = 1. We find it convenient to write the sums in (3.2) by decomposing the contribution coming from the i-th row of the matrices V (l) and T (l) , as also done in [21,25] in the case of d = 1. This choice leads to Following the analysis made in [22,25] for d = 1, we conjecture the existence of the limits where k 0 and η = ±1 are discrete parameters and Notice that i ± k/2 in (3.5) corresponds to the midpoint between the i-th and the (i ± k)-th site along the discretised radial direction. Since T (l) and V (l) are symmetric matrices, their diagonals labelled by +k and −k in (3.5) contain the same elements.
In our numerical analyses we considered lattices along the radial direction made by R tot ∈ {200, 400, 600, 800} and kept R/R tot 1 fixed (we choose R/R tot = 1/10). This choice is imposed by the fact that the continuum result that we are studying holds in the limit where the volume of the space is infinite. As already remarked in previous studies on entanglement Hamiltonians in free chains [21,22,25], high numerical precision is usually needed. In our case, they are required to evaluate (2.19) because many symplectic eigenvalues of Q l,B ⊕ P l,B are very close to 1/2 and they must be distinguished from 1/2 in order to be employed in (2.20). Higher precision is needed as l or d or ω increases. Thus, for d = 2 and d = 3 we have restricted our numerical analysis to l 10 and ω 5, working up to a precision of 2000 digits for the highest values of these parameters.
In Fig. 1 and Fig. 2 we show some numerical results for the elements in and near the main diagonals of the matrices T (l) and V (l) when ω = 0. From previous analyses for d = 1 [25,26], they are expected to be extensive; hence we consider T i,i+k /R (only the data corresponding to k ∈ {0, 1, 2, 3} are reported). The data points in Fig. 1 and Fig. 2 display good collapse close to the boundary of the sphere B as R increases. For small values of l, this agreement is observed almost all over the range of the spatial index except near the center of the sphere. These data collapses provide some numerical support to the conjecture (3.5) for small values of l. Comparing the left and the right panels in Fig. 1 and Fig. 2, we do not find significant differences between d = 2 and d = 3. As l increases, we expect that larger values of R are needed to observe good collapses of the data points. This lack of convergence for high values of l and our limited capability to treat numerically large systems forces us to focus only on small values of l.   We consider the continuum limit given by a → 0, R → ∞ and R tot → ∞, while Ra ≡ R and R tot a ≡ R tot are kept fixed (hence R/R tot is fixed as well). The parameters R and R tot are the radii respectively of the sphere B and of the entire system, which is a larger concentric sphere.
In the continuum, the radial position is labelled by r = ia with 0 < r < R. This leads to write (3.6) also as which suggests to expand τ l,ηk and ν l,ηk in (3.5) as Taylor series when a → 0. In order to study the continuum limit of the operators in (3.3) and (3.4), first we write these sums , for the operatorsq l,i+ηk andp l,i+ηk in the continuum limit we havê where the Taylor expansions of the fields as a → 0 have been used. Combining (2.12), (3.8) and ( where k max parameterises the number of diagonals included in the sums. The parameter k max plays an important role throughout our numerical analysis. In the continuum limit, also k max is infinite and therefore all the diagonals should be taken into account. However, in order to be consistent with (3.5), where R → ∞ for any finite 0 k k max , in our numerical analysis (where both R and k max are finite) we have to consider k max R. In the appendix B the role of k max is further discussed.
Since a → 0 in the continuum limit, we expand the integrands in (3.10) and (3.11), keeping only the terms that could lead to a non vanishing contribution after the limit. The expansion of (3.10) gives where we have introduced being k (±) max defined as follows In the top panels of Fig. 3 we show numerical results supporting the evidence that the limit (3.13) leads to a well defined finite function (see (3.27)) when k max is large enough, at least for the small values of l explored.
The expansion of (3.11) gives  Fig. 15). The horizontal dashed lines correspond to (2.11), i.e. to l 2 − 1/4 in the left panels and to l(l + 1) in the right panels. In the bottom panels the data corresponding to 1 l 10 are shown by adopting the logarithmic scale for the vertical axes.
where O(a) terms have been neglected and The numerical results displayed in Fig. 4 indicate that the limit (3.17) provides a well defined finite function when k max is large enough (see (3.28)), at least for small the values of l that we have considered.
Since the integrand of the O(1/a) term in (3.16) is proportional to the total derivative ∂ r [ν l,k,η (r) Φ(r) 2 ], the corresponding integral gives the boundary terms [ν l,k,η (r) Φ(r) 2 ]| r=R r=0 .  One of these boundary terms vanishes because ν l,k,η (R) = 0, while the other one does not contribute because of the Dirichlet boundary condition Φ l (0) = 0 imposed at the origin. Thus, (3.16) simplifies to The last term of this expression, whose integrand is ν l,k,η (r) Φ l (r) Φ l (r), can be studied by employing (3.5) and introducing In the bottom panels of Fig. 3 we display numerical results indicating that the limit (3.20) gives a well defined finite function when k max is large enough (see (3.27)).
As for the term in (3.19) whose integrand is ν l,k,η (r) Φ l (r) Φ l (r), since the analytic expressions for ν l,k,η (r) are not known, we approximate ν l,k,η (r) through finite differences, i.e. by replacing this function with [ν l,k,η (r + a) − ν l,k,η (r)]/a. This approximation, combined with (3.5) and (3.19), leads to introduce 1,kmax (i) being defined as follows Similarly, for the term in (3.19) whose integrand is ν l,k,η (r) Φ 2 l (r) we approximate ν l,k,η (r) through finite differences. This leads to define The subindices 1 and 2 in the l.h.s.'s of (3.22)-(3.23) and of (3.24)-(3.25) respectively indicate that the corresponding quantities are related to ν l,k,η (r) and ν l,k,η (r). In the top and bottom panels of Fig. 5 we show some numerical results telling us that the limits in (3.22) and (3.24) respectively give the function that vanishes identically (see (3.29)).
Finally, by employing the expressions (3.13), (3.20), (3.22) and (3.24) as discussed above, we take the limit k max → ∞ in (3.12) and (3.19), finding for the non vanishing contributions the following expression where it is assumed that the weight functions are well defined. Since the CFT expression (1.1) is valid when R tot → ∞, we must consider R R tot in order to compare our numerical results with this CFT formula specialised to the massless scalar field. The main outcomes of our numerical analysis are shown in Fig. 3, Fig. 4 and Fig. 5, where the data reported in the left panels and in the right panels correspond to d = 2 and d = 3 respectively. These results provide some numerical evidence supporting the conjecture that (3.26) provides the CFT prediction (1.1) for the massless scalar field. In Fig. 3 the combinations of diagonals defined in (3.14) and (3.21) are considered. The numerical data shown in this figure lead us to conjecture that where β(r) is the parabola (1.2) restricted to 0 r R, which is independent both of the dimensionality parameter d and of the mode parameter l.
In Fig. 4 we report numerical data points for the combination of diagonals (3.17) which support the following conjecture The horizontal dashed lines in both the panels of Fig. 4 correspond to the coefficient µ d (l) defined in (2.11). We find it worth highlighting that, although the diagonals shown in Fig. 2 for d = 2 and d = 3 seem identical, their combination (3.18) displays the peculiar dependence on d given by (2.11), as shown in Fig. 4. This is a characteristic feature of the fact that we are considering an entanglement Hamiltonian in a Minkowski spacetime with a number of spatial dimensions strictly larger than one; indeed, the term corresponding to (3.28) gives the vanishing function when d = 1. Near the boundary of the sphere, i.e. where i/r ∼ 1 − , the discrepancy between the data points (which are obtained as a ratio of two quantities that are both vanishing at the boundary of the sphere) and (3.28) increases with l. These discrepancies becomes very similar when the logarithmic scale is adopted (see the bottom panels of Fig. 4).
In Fig. 5 we show numerical data for the combinations of diagonals introduced in (3.23) and (3.25). In this case the curves for different sizes R do not collapse and tend to zero as R increases; hence it is natural to conjecture that We emphasise that large values of k max are needed to obtain the numerical results described above. In appendix B we show that the expected CFT results are not obtained when k max is not large enough. This crucial message can be appreciated e.g. for (3.27) and (3.28) by comparing the right panels of Fig. 3 and Fig. 4 with Fig. 14 and Fig. 15 respectively. The conjecture (3.28) naturally leads to ask how k max → ∞ should be taken in the continuum limit. In Fig. 3, Fig. 4 and Fig. 5 the ratio k max /R is kept fixed (and equal to 1/10) as R increases. In the appendix B we report and discuss also numerical results obtained by keeping k max fixed as R increases (see Fig. 14 and Fig. 15). We find that these two different limiting procedures give the same results for the quantities that we are considering. However, this question deserves further investigations. Notice also that all the data reported in Fig. 3, Fig. 4 and Fig. 5 do not display good collapses around the center of the sphere (i.e. where i/R ∼ 0 + ): we expect that larger systems are needed to observe them.
The conjectures (3.27), (3.28) and (3.29) have been formulated for any value of l. However, the data points in Fig. 3, Fig. 4 and Fig. 5 provide numerical support to these conjectures only for the small values of l that we have been able to explore. The validity of (3.27), (3.28) and      We find worth discussing also some numerical results for the entanglement entropy. Since the ground state is the direct product of the ground states corresponding to H l , the entanglement entropy of the sphere S B is obtained by summing the contributions corresponding to all the different values of l [2]. These contributions can be evaluated through the symplectic eigenvalues of the reduced covariance matrix Q l,B ⊕ P l,B at given l or, equivalently, through the single particle entanglement energies at given l [33,34,37,[50][51][52], as briefly anticipated in Sec. 2.2. Since Q l,B ⊕ P l,B depends only on l, the entanglement entropy of the sphere B is computed as follows [2] which has been explored numerically in various studies [2,9,[41][42][43]45].
In Fig. 6 we show the single particle entanglement energies ε l,k /R at given l in terms of k/R when ω = 0, for d = 2 (left panel) and d = 3 (right panel). The summand in S B,l diverges as ε l,k → 0 + ; hence the low-lying part of the single particle entanglement spectrum at fixed l provides the largest contribution to S B,l . Clear differences between ε l,k /R for d = 2 and d = 3 are not visible, despite the fact that Q l,B ⊕ P l,B depends explicitly on d (see (2.14)).
In the left panel of Fig. 7, we show some numerical data for S B,l defined in (3.30) which display a clear dependence on the dimensionality parameter d for this quantity. In this panel the dashed curve corresponds to 1 6 log(R) + const, which is the entanglement entropy of a segment at the beginning of the semi-infinite line [6] and it coincides with S B,0 when d = 3, as expected from the fact that µ 1 (0) = µ 3 (0) = 0. Moreover, the data corresponding to (d, l) = (2, 1) coincides with ones corresponding to (d, l) = (4, 0) because (2.11) gives 3/4 Many interesting results have been obtained for the entanglement entropy of the sphere and other interesting subregions in the continuum limit for generic d through various quantum field theory mehods [6,49,[53][54][55][56][57][58][59]. It would be instructive to explore whether these methods can be employed to study also the corresponding entanglement Hamiltonians.

Entanglement Hamiltonian at ω > 0
In this section we consider the entanglement Hamiltonian of the d dimensional sphere B in the regularised model of [2] introduced in Sec. 2 when ω > 0 and the entire system is in its ground state. We adapt to this case the analysis performed in [26] for the entanglement Hamiltonian of a block made by consecutive sites in the infinite and non-critical harmonic chain. The numerical setup is the one described in Sec. 3 and, since the outcomes for d = 2 and d = 3 are very similar, in the following we report only the results corresponding to d = 3.
In Fig. 8 and Fig. 9 we show the numerical data points for the diagonals of the matrices T (l) and V (l) respectively, evaluated from (2.19) when d = 3 and for ω = 1 (left panels), ω = 3 (central panels) or ω = 5 (right panels). It is instructive to compare these results against the corresponding ones obtained for ω = 0 and displayed in Fig. 1 and Fig. 2. Also for ω > 0, data collapses are observed for small values of l. Presumably, larger systems are needed to observe these collapses also for higher values of l.
The considerations made in [26] for the entanglement Hamiltonian of a block of consecutive sites in the infinite and non-critical harmonic chain can be adapted to this case in a straightforward way, with the crucial difference that only one point separates the subsystem      i,i+1 we can identify a region closed to the boundary of the sphere where the diagonal vanishes. For a given l, the size of this region increases with k at fixed ω and it also increases with ω at a given k. Thus, in the large mass regime only the three diagonals T i,i+1 are non vanishing. Furthermore, when ω l and for the small values of l considered, we find that the analytic results found in [26] (which are based on [27]) can be easily adapted to the case of the sphere. In particular, the dominant matrix elements are well approximated by               where with I(κ) being the complete elliptic integral of the first kind 3 and The analytic results in (4.1) and (4.2) correspond to the straight green dashed lines in Fig. 8 and in Fig. 9 respectively. The massive scalar field in the continuum limit is described 3 The integral representation of the complete elliptic integral of the first kind is by taking ω → 0 and a → 0 while ω/a is kept fixed. When ω → 0, we have that κ → 1 and κ → 0; hence I(κ) → π/2 and therefore 2b(κ) → 2π, which is the value predicted by Bisognano and Wichmann [10,11], as already remarked in [26]. Notice that the analytic expressions in The symplectic spectrum of the reduced covariance matrix Q l,B ⊕ P l,B at a given l provides the corresponding single particle entanglement energies ε l,k at fixed l as discussed in Sec. 2.2. In Fig. 10 we show numerical data for ε l,k (in the case of d = 3) when l is small for some values of ω. By comparing these results with the corresponding ones for ω = 0 displayed in Fig. 6, we observe that ε l,k as function of k/R are well described by a straight line as ω increases. This feature has been already highlighted for d = 1 in the case of the block of consecutive sites in the infinite harmonic chain on the line [26]. Furthermore, the slope of this straight line is given by the same expression found in [26] for the d = 1 case, namely [28] ε l,k = (2k − 1) ε ε ≡ π I(κ) I(κ) (4.6) where κ andκ have been defined in (4.5).
Numerical results for the entanglement entropy S B of the sphere when ω > 0 are shown in the right panel of Fig. 11 (see the right panel of Fig. 7 for S B when ω = 0). They have been obtained through (3.30), with the sum over l restricted to l 500 for ω = 1 and ω = 3, and to l 1000 for ω = 5. The dashed lines correspond to α(d, ω) R d−1 for some fitted values of the constants α(d, ω); hence these data just highlight the expected area law behaviour of the entanglement entropy.
Similarly to the massless case, also when ω > 0 the dependence on d and l is more visible in the quantity S B,l defined in (3.30). In the left panel of Fig. 11 we show S B,l in terms of R for ω = 1 and ω = 3, finding a qualitatively different behaviour with respect to the massless case (see the left panel of Fig. 7). In particular, the horizontal dashed line providing the asymptotic value of S B,l in the left panel of Fig. 11 is given by [8,26] (4.7) We remark that the above numerical analysis does not correspond to the continuum limit of the entanglement Hamiltonian for the massive scalar field, where ω → 0 and R → ∞ while ωR is kept fixed. In this regime we have encountered the same difficulties discussed in [26] for the entanglement Hamiltonians of a block of consecutive sites in non-critical free chains on the infinite line.

Conclusions
We have explored the continuum limit of the entanglement Hamiltonian of a sphere in the d+1 dimensional Minkowski spacetime for a massless scalar field. We have employed a numerical analysis based on the radial lattice discretisation introduced in [2], on the results found in [2,9] and on the procedure already introduced [21][22][23][24][25] to study some entanglement Hamiltonians of a block of consecutive sites in one-dimensional free chains. Our main results for the massless scalar field are the conjectures given by (3.27), (3.28) and (3.29), which are supported by the numerical data reported in Fig. 3, Fig. 4 and Fig. 5 in the cases of d = 2 and d = 3 and for small values of the total angular momentum parameter l. By employing these conjectured results into (3.26), we have obtained the expression (2.17) in the continuum, which can be derived also from the entanglement Hamiltonian of the sphere for a generic CFT in its ground state [12,13] (see (1.1)) specialised to the massless scalar field. An explicit dependence on the dimensionality parameter d and on the total angular momentum parameter l occurs in the term (3.28), which vanishes identically when d = 1. In the special case of d = 1, we recover the results of [25] for the entanglement Hamiltonian of an interval at the beginning of the semi-infinite line with Dirichlet boundary conditions at the origin.
As for the massive regime, we have discussed numerical results obtained for a given ω > 0 (see Sec. 4). At a generic value of ω 0, the matrix (2.19) characterising the quadratic entanglement Hamiltonian contains long-range and inhomogeneous couplings (see Fig. 8 and Fig. 9). However, in the limit of ω l, only the nearest neighbour couplings are non vanishing and the corresponding weight function along the radial direction is well approximated by straight lines whose slopes are independent of d and l and can be determined analytically (see (4.1) and (4.2)), as already observed in [26] also for the entanglement Hamiltonian of a block of consecutive sites in the harmonic chain on the infinite line and in its ground state.
Our analysis can be improved in various directions. All the numerical checks performed in this manuscript correspond to small values of the total angular momentum parameter l; hence an improved numerical analysis is required to explore also the sectors corresponding to higher values of l. For the massless scalar, the existence of the functions τ l,ηk and ν l,ηk introduced in (3.5) is a crucial assumption throughout the derivation of (3.26). It would be interesting to obtain numerical data that support further these conjectures and also to find analytic expressions for these functions, as done in [22] for the entanglement Hamiltonian of a block of consecutive sites in the chain of free fermions on the line. For the massive scalar, it is important to explore through these numerical methods the regime characterised by a given value of ωR when R → ∞, where the entanglement Hamiltonian of the sphere in the continuum is fully non-local [21,60].
It is interesting also to develop methods to write lattice operators including only nearest neighbour couplings that approximate the entanglement Hamiltonians [17,19,22,24,[72][73][74][75] and to understand their relation with the corresponding entanglement Hamiltonian, as done e.g. in [19,22]. A numerical approach to some entanglement Hamiltonains based on a quantum Monte Carlo method has been proposed in [76].

A Correlators in the continuum
In this appendix we report the two-point functions in the continuum computed through quantum field theory methods [38,39] and compare them with the correlators on the lattice (2.15) obtained for different discretizations of the Hamiltonian of the scalar field.
In the spherical coordinates introduced in Sec. 2.1, the two-point function of the free massive scalar field Φ(t, r, Ω) in the d + 1 dimensional Minkowski spacetime reads [39] Φ(t 1 , r 1 , Ω 1 ) Φ(t 2 , r 2 , Ω 2 ) = 1 is the area of the d − 1 dimensional unit sphere and In (A.1) the angle between directions identified by Ω 1 and Ω 2 is denoted by ξ and C q p (x) is the Gegenbauer polynomial of degree p and order q, which can be expressed also through the hypergeometric function as follows (see e.g. section 15.4 in [30]) where (a) b is the Pochhammer symbol.
It is instructive to consider the special case of d = 2, which is the lowest value of d where the sum over l in (A.1) is non trivial. Taking the limit d → 2 in the infinite sum obtained by isolating the l = 0 term in (A.1) and using thatμ = l when d = 2 (from (A.3)), one finds When m = 0 and for r 1 < r 2 (which leads to introduce ρ 12 ≡ r 1 /r 2 < 1), this integral reads which holds for d 2; indeed, for d = 1 and l = 0 we haveμ = −1/2, where this expression is not well defined. From (A.2) we can compute also ∂ t 1 ∂ t 2 G l (t 1 , r 1 ; t 2 , r 2 ) for t 1 = t 2 and the result is  Figure 12: Correlators (2.15) in the massless regime: Q i,i+10 (top panels) and Q i,i+100 (bottom panels) for d = 2 (left panels) and d = 3 (right panels) evaluated from (2.14) with ω = 0 and for systems whose total size is given by R tot ∈ {600, 800}. The dashed lines correspond to (A.7) and the solid ones to (A.11) with m = 0. The magenta data points in the left panels are obtained from (A. 16) with ω = 0.
In Fig. 12 and Fig. 13 we compare the numerical data points obtained for the correlation matrices given by (2.15) and (2.14) with the corresponding expressions in the continuum, which are given by (A.7) and (A.9) for the infinite volume regime (dashed lines) and by ○ Figure 14: Combinations (3.14) and (3.21) when ω = 0 in the case of d = 3, for various k max . These panels should be compared with the right panels of Fig. 3.

B Comments on the role of k max
In this appendix we briefly discuss the role of k max in the summations T Since very similar results are obtained for d = 2 and d = 3, we report only the latter ones. In Fig. 14 we show numerical data for the combinations (3.14) and (3.21), while in Fig. 15 the ones for the combination (3.18) are displayed. The data corresponding to l = 0, l = 2 and l = 5 have been reported in the left, middle and right panels respectively. In each panel, increasing values of the summation parameter k max ∈ {1, 2, 10} are considered.
The aim of these figures is to show that including more diagonals in the summations brings the corresponding numerical curve closer to the expected curve, predicted by CFT. This is in agreement with similar analyses made in previous studies for one dimensional free chains [21][22][23][24][25]. For the systems that we have explored, in all the combinations of diagonals convergence is attained already when k max ∼ 10 for all the values of l considered. Including more diagonals does not lead to visible improvements. In T (l,0) kmax the convergence is faster than the one observed in the other combinations of diagonals; indeed already the main diagonal of T (l) is already close to the expected CFT result (1.2) (see the top panels of Fig. 1). Instead, for V (l,2) kmax and V (l,0) kmax it is evident that k max > 2 is needed to recover the corresponding CFT prediction. In particular, while the data of V (l,2) kmax display a good agreement with the CFT curve already when k max = 2 (see the bottom panels in Fig. 14), the data for V (l,0) kmax in Fig. 15 clearly tell us that higher values of k max must be considered to recover the expected CFT prediction, which is characterised by the horizontal dashed red lines obtained from (2.11).