Euclidean Time Approach to Entanglement Entropy on Lattices and Fuzzy Spaces

In a recent letter, we developed a novel Euclidean time approach to compute R\'{e}nyi entanglement entropy on lattices and fuzzy spaces based on Green's function. The present work is devoted in part to the explicit proof of the Green's matrix function formula which was quoted and used in the previous letter, and on the other part to some applications of this formalism. We focus on scalar theory on 1+1 lattice. We also use the developed approach to go systematically beyond the Gaussian case by considering interacting models, in particular our results confirm earlier expectations concerning the correction to the entanglement at first order. We finally outline how this approach can be used to compute the entanglement entropy on fuzzy spaces for free and interacting scalar theories.


Introduction
Quantum entanglement plays fundamental roles in different areas of physics and future technology, it is generic and ubiquitous, starting from the role it plays in quantum information, being primary resource for quantum computation and information processing, to condensed matter, black holes and high energy physics, see [1][2][3][4], .
Moreover, in view of several recent development and progress made in connection with holography and strongly-couple quantum field theories that have led to some inspirational ideas in defining quantum gravity from quantum many-body entangled states, it seems that quantum entanglement and holography will be two major pieces of our understanding of the fundamental nature of space-time.
Many measures have been developed to quantify the degree of entanglement between different parts of a quantum system [1,4,5]. One of the most important and fruitful measures is of course the entanglement entropy (EE).
In connection with black hole entropy and quantum field theory the finding and observations made in [6,7] led subsequently to a wide interest in EE in the presence of black hole, as one would have to expect some amount of entropy to accompany an event horizon, since it is by definition an information hider.
Several techniques have been devised to compute this entropy depending on the area of research field. Generally those techniques fall into two categories. Real time and Euclidean time formalisms [1].
The real time approach was initially developed in [7,8], during the quiet period of the subject, and later reformulated with deeper insights using the two point correlators of the field variables [9]. This formalism is very suitable for numerical computation once the underlying space is discretized to lattice or for a bosonic or fermionic chain, however this approach seems to be limited to gaussian states or free field theories.
The results of [7,8] trigged attention to the EE in the black hole context and soon an Euclidean time approach was developed in [10,11] based on Euclidean heat kernel and Green's function. This formalism proved to be more powerful in performing analytical computation in QFT and extracting the UV divergence. Moreover it allows one to make connection with string theory [12], extract finite mass dependent but cuttoff independent contribution to the area law [13]. Unlike the real time approach the Green's function approach can be used to go beyond the free case and investigate the dependence of entanglement entropy on the renormalized mass [14].
For 1 + 1 QFT, the tools of CFT were very early employed to investigate the EE [15], ever since a lot of work and substential advances have been made in the understanding of EE in 1+1, see [3,16] and references therein. Now, although the heat kernel and the Green's function approaches offer a powerful tool to compute and investigate the entanglement entropy they seem to be limited to some special cases of continuum geometries. For example when one traces out half of flat d + 1 spacetime Rényi entropy turns out to be governed by a simple geometry, namely a cone with deficit angle 2π(1 − n), for which the Green's function is known or the heat kernel expansion can be made. However, if one traces out a finite region it becomes almost untractable problem to construct the corresponding Green's function or to recognize the resulting n-sheeted geometry.
Motivated by the absence of an Euclidean formalism for the EE on lattices and fuzzy spaces and the limitations of the real time approach to free and gaussian states an Euclidean time approach based on Green's function has recently been proposed in [17]. An alternative formula to compute Rényi EE resulting from tracing out an arbitrary set of intervals was obtained for a one-dimensional chain of coupled H.O's.
By its very formulation the proposed approach furnishes a setting on which a systematic perturbation series for the EE in interacting theories can be obtained.
The fundamental piece upon which this approach rests is the Green's matrix function which was quoted in [17] and the proof of its construction was deferred to the present work. In this paper we show how the Green's function is constructed for an arbitrary set of intervals. We further apply the formalism to a scalar field theory on 1+1 lattice by considering the EE in two important particular cases, one resulting from tracing out half of the space and a second resulting from tracing out an interval of infinitesimal size, in particular an interval made of a single point. For both cases we perform analytical study.
From the technical point of view the case of one-point turned out to be easy to handle analytically within this approach and we could recover explicitly the formalism developed in [9]. However, the case of half-space EE necessitated the finding of an asymptotic inverse for a special Toeplitz matrix with N-dependent symbol. This problem turned out to be a strenuous one due, as far as we are aware, to the lack of any solution for this kind of problem in mathematical literature. Therefore we hope that our modest contribution will pave the way to a better understanding of this class of Toeplitz matrices.
As another important application of the present technique we consider the correction to the EE due to a quartic type interaction for a system made of finite number of coupled oscillators. In particular we derive explicitly the first order correction to the standard EE between two coupled harmonic oscillators due to quartic interaction, λ q 4 i , and show that to the first order in perturbation theory the standard gaussian (free) EE formula stays valid for interacting theories ( non-gaussian) with a shifted parameter due to the perturbation of the ground state. Our result systematically confirm earlier conjecture and expectations that the EE is generally a space-time quantity and captured by the perturbed two-points correlators to the first order in perturbation theory [18] .
In addition to the above quartic interaction, to which we shall refer as of local origin, we consider another natural quartic interaction which can not be seen as arising from a discretization of a non-local term. In contrast to the local interaction this quartic interaction appears to have an IR divergence that can not be normalized away and seems to limit the type of quartic interactions that the entanglement can tolerate .
The present work is organized as follows. In the first section we begin by a subsection in which we review the formalism developed in [17]. The second subsection will be devoted to the detail of the construction of Green's matrix function. In section II we consider a massive scalar theory in 1+1 lattice and consider the EE for two cases, half of the space and one-point, and as a by-product we obtain the standard formula of the EE for two coupled H.Os. A subsection is devoted to the mathematical problem of finding an asymptotic inverse of a particular class of Teoplitz matrices with N-dependent index. In the third section we go beyond the gaussian states by including quartic interaction and derive a general formula for computing the first order correction in perturbation theory. We give explicit first order evaluation of the correction to the EE for 2-coupled H.O and show that for a typical quartic interaction the non-gaussian formula for EE keeps the same form as the standard gaussian one and is captured by the perturbed two-point correlators. We conclude our paper by an outlook on the applications of the developped formalism to free and interacting scalar field theories on fuzzy spaces.
Most of the technical detail of the paper are moved to the appendices section.

The Green's Function Matrix for Arbitrary Set
To make this paper as self-contained as possible we start by reviewing the formalism developed in [17].

The Formalism
Consider a system of N -coupled H.O's with standard quadratic Lagrangian V is a symmetric positive definite N × N matrix. By considering N coupled H.O we cover a wide class of field theories on lattices and fuzzy spaces. This follows from the observation that many physical situations are modeled by a chain of coupled H.O's ( Fermionic and Bosonic), whereas in high energy physics it turns out that for most theories regularized on lattices and fuzzy spaces computing the entanglement entropy reduces to considering independent sectors made of coupled H.Os, and this is independently of the dimension of the spacetime considered and the geometry of the region traced out [1,8,19,20] Let now ρ = |0 >< 0| be the ground state of whole system and consider the reduced density operator ρ A = Tr A ρ resulting from tracing out a subset A of the field variables Q. A could be a union of collection of disjoint subsets . LetĀ denote the complimentary subset.
For pedagogical reasons we shall consider the important special case where A = {q p+1 , q p+2 , .....q N }. It will become clear as we move on that the proof of the final general result is not tied to this particular choice.
Let Q A ≡ (q p+1 , q p+2 , .....q N ) T and QĀ ≡ (q 1 , q 1 , .....q p ) T . The density matrix elements of the ground state can be written using Euclidean path integral representation as The reduced density matrix is then given by, from which it follows that Now let us introduce an n-fold field variable Q with N n components and n-fold potential V as then it is straightforward to write where L E is given by Unlike the continuum flat case we can not interpret the above path integral as the partition function on the n-fold cover manifold (with a conical singularity) of the original manifold. However, it turns out that it is still possible to write TrĀρ n A as a zero temperature partition function with specified Lagrangian.
Let us note that the main obstacle which is preventing us from interpreting (2.6) as a zero temperature partition function is the boundary conditions on the variables QĀ at the cut. This can be surmounted by moving the cut, or the discontinuity, from the integration variables to the potential via a permutation matrix P π which maps (Q 1Ā , Q 1 . Then TrĀρ n A with the correct normalization can be written as it follows immediately that where V(τ ) is a n-fold step-potential given by V and its image V p are N n × N n matrices defined as follows V = I n×n ⊗ V and V p = P π VP T π (2.12) P π is N n × N n permutation matrix, its explicit form is given for the above particular choice of A by P A and PĀ are the projectors on A andĀ resp. More explicitly we have Actually it is easy to see that all permutations matrices for this kind of problems are made of projectors on the available and the unavailable degrees of freedom and zeros.
If we decompose the potential accordingly as A and C are p × p and (N − p) × (N − p) matrices resp. The n-fold potential is given for arbitrary n by The form of this supermatrix is typical in the calculation of the entanglement entropy using Euclidean formalism [10,21], and here it is understood as arising from the action of a special permutation matrix. Now, as mentioned earlier, although the derivation of (2.10) was explicitly made for a particular choice of the subset A, it remains valid for any other choice, one has to only change the permutation matrix. Therefore each tracing operation is characterized by a its permutation matrix P π (A,Ā), and all we need is its action on V or the image potential V p .

Constructing Green's Matrix Function
This section is devoted to the proof of the main formulas which appeared in [17], namely the explicit construction of the Green's function associated with this type of problems.
We start by expressing (2.10) using the heat Kernel matrix associated with the operator A n = − d 2 dτ 2 + V(τ ). Using the eigenvalues and eigenfunctions ofÂ n we can write The eigenfunctions satisfy the following orthogonality and completeness conditions.
Beside these standard conditions the eigenfunctions must also satisfy certain continuity conditions. In view of the fact that the jump in the potential matrix is finite, the eigenfunctions and their first derivative must be continuous at τ = 0. These conditions will be translated into continuity conditions which we must impose on the associated heat kernel and the Green's function ( its Laplace Transform), and they will play a crucial role in constructing the explicit form of the Green's function. Using (2.16) and (2.17) we can write (2.10) as Define now the Laplace transform of the heat kernel K n , G n (τ, τ ′ , E) = ∞ 0 e −Es K(τ, τ ′ , s)ds. The N n × N n matrix function G n (τ, τ ′ , E) is just the Green's function associated with − d 2 dτ 2 + V(τ ) + E satisfying the following differential equation it is then straightforward to show The trace Tr is understood to include integration over τ 1 .
Let us now see how G n (τ, τ ′ ) can be constructed. In view of the fact that V and V p do not generally commute, the construction of the Green's matrix is not straightforward. The construction given here is actually slightly different from the original line of thought that led us to the correct form of this Green's matrix function, but it is more transparent and pedagogical.
1 Equation (2.20) is indeed equivalent to the evaluation of ln Trρ n A for massive scalar theory by first evaluating the derivative of ln Zn with respect to m 2 , with the role of m 2 is here played by E, expressing it in terms of the coincident Green's function Gn(τ, τ, m 2 ) on the n-fold cover space and then integrating with respect to m 2 [4] To construct G n (τ, τ ′ ) we first start by considering the special case where V = V and V p = V p are just real scalars . In this case the problem is reduced to finding the Green's function for one-dimensional step potential for Schrödinger equation (with an appropriate redefinitions of the mass and other parameters). The Green's function for one dimensional step potential has been derived using different techniques [22,23]. Using those results the Green's function is given in the case V = V and V p = V p by, and where we used the standard notation G +− (τ, τ ′ ) = G n (τ, τ ′ ) for τ > 0, τ ′ < 0 ..etc and with Consider now the case where V and V p are two commuting matrices. In this case the problem is practically equivalent to the one-dimensional case, and the Green's function keeps formally the same form without any ambiguity. Of course 1 W and 1 W − +W + would stand for the inverse of the matrices W, and W − + W + , however we do not need to worry about the ordering of the different matrices appearing in equations (2.21) and (2.22).
For the general case where V and V p are non-commuting, equations (2.21) and (2.22) do not work as they stand and the problem seems more challenging.
We note here that we are facing two possibilities. Either we have just to arrange the matrices so that the resulting Green's matrix is still a solution and satisfies all the continuity and boundary conditions, or consider the possibility that extra terms are present and which vanish in the commutative case, or possibly considering more complicated form that involves chronological ordering operations. Fortunately it turns out that all we need is to use the Green's matrix for commutative case and rearrange the matrices multiplications to obtain the correct solution.
We first note that although we were unable to obtain G(τ, τ ′ ) by directly solving the differential equation (2.19), it is possible to write down a series expansion for G(τ, τ ′ ) around the free one V(τ ) = 0 . This is achieved by writing (2.19) as an integral equation and using Laplace transform method to turn it into Riemann-Hilbert problem. Using the technique known for these kind of problems one can show that G(τ, τ ′ ) has the following expansion up to the second order in V(τ ) , for τ, τ ′ < 0, Similar expansion could be obtained for G ±∓ . The detail of the derivation of the above expansion is given in Appendix I. Now it is not difficult to show that the expansion (2.23) and those for other regions coincide with the expansion of the following expression, and where W ± are given by The above matching between the expansion and the forms given by equations (2.24) and (2.25) does not consist a full proof, nevertheless it does almost bring us home. Actually we can forget about the root of (2.24) and (2.25) and limit our effort to showing that they consist the correct solution of (2.19).
Let us remember that the Green function we are seeking for is uniquely defined by the following requirements: • a-It satisfies the following equation • c-All continuity conditions are satisfied.
Requirements (a) and (b) are obviously satisfied. The continuity conditions are not obvious and needed to be listed and checked one by one.
These conditions follow from the continuity conditions that the eigenfunctions and their first derivative must satisfy in view of the fact that the potential has a finite jump at the cut τ = 0.
In Appendix I we list the resulting continuity conditions and show that they are all fulfilled .
It should be noted here that beside these standard continuity conditions there is an extra condition that the solution must satisfy, namely G(τ, τ ′ ) † = G(τ ′ , τ ) . It follows from (2.16) and we shall refer to it as the hermiticity constraint . In Appendix I this condition is shown to be fulfilled by G(τ, τ ′ ) of (2.24) and (2.25). The hermiticity constraint reduces by half the number of continuity conditions we must check, we will have to only consider the continuity with respect to one argument of G(τ, τ ′ ).
This completes the construction of the Green's function. Substituting now G n (τ, τ ) into eqt (2.20) and performing the integration with respect to τ we obtain This our alternative formula which gives Rényi entropy resulting from tracing out an arbitrary subset of H.Os, and it can be used to compute the entanglement entropy via the replica trick.
Some remarks about formula (2.27) are in order. The expression of ln Trρ n A is invariant under rescaling of the potential and vanishes when [V, P π ] = 0 as it should do. Actually it is the nonvanishing of the commutator [V, P π ] which will lead to non zero entanglement entropy, therefore [V, P π ] / V can be considered in a sense a measure of the degree of entanglement. The matrix norm . could be the Hilbert-Schmidt norm or whatever matrix norm one prefers. Actually it turns out that when using the weak norm it is more appropriate to rescale the commutator and use instead √ N [V, P π ] / V as a measure of the degree of entanglement. Finally, we note that in principle one can develop a similar formalism for a Fermionic lattice system by using Euclidean path integral based on Grassman variables, imposing antiperiodic conditions and derive the corresponding Green's matrix and EE formula . This would offer a new tool to investigate the EE for several Fermionic models in connection with critical phenomena and RG, like the Long-Range Kitaev chain and other models.

Massive Scalar Theory on 1-d lattice
As a first application and good warm-up exercise we apply (2.27) to free massive K.G field on 1+1-d Minkowski spacetime. This model in the continuum has been investigated using different techniques. As a by-product of this problem we shall rederive the standard formula of the EE for two coupled H.O's.
The Lagrangian is given by We put this model on a lattice as follow . The x−axis is replaced by a one-dimensional lattice, i.e x i = iǫ where ǫ is the lattice spacing. As an infrared cuttoff we take the size of the lattice to be R = 2N ǫ and impose periodic boundary conditions .
By rescaling the field variable we obtain a dimensionless Lagrangian The potential matrix V is 2N × 2N matrix given by Where µ = mǫ. We note that V is a circulant symmetric real matrix and therefore its elements can be written using its eigenvalues as

The Half-Space Entanglement
Let us consider the entanglement entropy resulting from tracing out the positive axis or Similarly for W − = I n ⊗ (V + E) , and W + = P π W − P T π . In view of the fact that V + E is circulant symmetric matrix, both (V + E) −1 and (V + E) are also circulant symmetric matrices which can be decomposed as A, B, G, F are N × N matrices. Now, putting everything together and using (2.15) it is straightforward to show that where C is a 2n × 2n circulant matrix with elements c j , j = 0, 1. · · · 2n − 1, c 0 = 2, c 2n−2 = c 2 = −1 and zeros otherwise. And Therefore we can use discrete Fourier transformation (DFT) to simultaneously block , perform partial inversion of (W + + W − ) and get where T = G + F cos πj n and we obtain for Rényi entropy The matrices elements of B, F and G follow from (3.4), where ϕ k = µ 2 +E +2(1−cos πk/N ), and l, m run now over the rang l, m = 0, . . . N −1.
Before we move on in studying the large N and small µ behavior of Rényi entropy for this model let us first use (3.9) in the special case N = 1 to obtain the standard EE formula for two coupled H.O's.
In the case N = 1 the potential (3.4) becomes which corresponds to a rescaled potential describing two coupled H.O's with equal frequencies. The case of two coupled H.O's with different frequencies is a little bit more complicated technically and can not be deduced from this model, but it can be worked out exactly [24]. Now, substituting N = 1 in (3.9) and using (3.10), (3.11) and (3.12) we obtain α ± = 1 ± cos πj n , a 0 = µ 2 /2 + 1. The integral inside the sum can be exactly evaluated and after some arrangements and simplifications we end up with Here our α turns out to be related to ξ defined in [8] by α = 1 √ ξ . Using the replick trick and the definition of α we recover the standard formula for EE for two coupled H.O's, Having derived the standard formula for EE for 2-coupled H.O let us now focus on the large N limit .
The main technical obstacle is the inverse of the matrix G + F cos πj n , therefore we shall dedicate a special subsection to this problem .
We first note that our goal is not to find an exact inverse to this matrix for arbitrary N , but we are only interested in the large N limit or an asymptotic inverse.
We start by rewriting equation ( 3.9) as where T N (α) is given in the large N limite by where α ± = 1 ± cos πj n , and c k = a − cos 2πk N .

On the Asymptotic Inverse of a Particular Class of Toeplitz Matrices
The matrix T N (α) defined by (3.19) and (3.20) is a particular Toeplitz matrix, a linear combination of a circulant and anti-circulant matrices. By anti-circulant matrix we refer to the matrixC l . Both C andC can be inverted separatly, however their linear combination has no known inverse. Moreover in the large N limit T N (α) does not admit a Wiener representation with N -independent symbol. Although a great deal is known about the asymptotic behavior of Toeplitz matrix with N -independent symbol under certain conditions, very little, if any, is known about the asymptotic behavior of Toeplitz matrices with N -dependent symbol, as it is the case for T N (α). This kind of matrices is expected to show up within this approach whenever we trace out half of the space, therefore it is necessary to understand the the large N or asymptotic behavior of these type of matrices and their inverses.
Before tackling the problem of finding an asymptotic inverse for T N (α), let us see what the available mathematical literature can say about this inverse.
It can be shown using the machinery of Generalized Locally Toeplitz (GLT) [25,26] that the trace of T N (α) −1 has the following asymptotic limit 2 lim N →∞ where K is the complete elliptic integral of the first kind. The above result can be easily be confirmed numerically. At N = 50 and for several values of a and n we get a perfect coincidence.
We shall not enter into the detail of the derivation of the above asymptotic trace of the inverse, because as it unfortunately turned out the GLT machinary and their * -algebra structure fails to give a good estimation of the full trace which appears in eqt (3.18), they only give weak information on this trace not strong enough to be of any use for our particular problem, namely being just of order one. Nevertheless we shall use the above result about the asymptotic behavior of the trace as a guide in our search for the inverse .
We begin our search for the inverse by noting that it must be a persymmetric matrix as can easily be shown. Analytical and numerical investigation of the structure of both T N (α) and its inverse led us to consider the following ansatz as a candidate for the inverse. where α pk is a function of p and k to be determined later, and k ± are two functions ( a sort of kernels) defined by It is easy to see that M N is persymmetric matrix for arbitrary α pk . Now, the ideal thing is to fix α pk by solving M N T N = I, for finite or atleast large N . Unfortunately this direct approach led to apparently unsolvable equations. However it turns out that a particular value for α pk solves the following necessary condition TrM N T N = N .
Using the form of T N and the ansatz M N it is a straightforward exercise to show that with α pk given by Although this particular choice of α pk only solves the trace condition ( for arbitrary N !) it turns out to be a good candidate for an asymptotic inverse for T N .
Before discussing the validity of M N , with α pk = 1 α + √ c k +α − √ cp , as an asymtotic inverse, some few remarks about this ansatz are in order 3 It is easy to show that M N gives the exact inverse for the two particular cases (α − = 0, α + = 0) and (α + = 0, α + = 0). For large a, the expansion of M N coincides with the neuman series of the inverse of T N , at least for the first few orders we checked explicitly.
Now, what about a close to one, the case of most interest to us?. We have run serval numerical calculations for the following norms (the weak norm or Hilbert-Shmidth).
and ∆ N (a) and as long as a is not very close to 1, a = 1.2 and larger, they are both smaller than 10 −3 for N = 30, and becoming really negigeable as we reach 1.5, larger than this there is almost a perfect matching (for a = 2 both δ and ∆ are of the order 10 −4 ). Of course they become smaller as we move to larger N .
For a close to one, a = 1.001 say, one needs to go to larger N to get small values for δ and ∆, however even for N = 40 and N = 50 the two norms are of the order of 1/N numerically. For a = 1.002, N = 100 , ∆ = 0.06. These numerical results provide strong evidences that both δ N (a) and ∆ N (a) are of the order of 1/N or atleast vanish in the large N limit, as long as a > 1. Therefore we shall make the following conjecture.
Conjecture: Let T N be a Toeplitz matrix given by where α ± = 1 ± cos πj n , and c k = a − cos 2πk N . then the following matrix By asymptotic inverse here we mean that M N fulfil the following condition We shall not try to prove this conjecture at this stage nor do we pretend to know how δ N (a) approaches zero, although our numerical calculation suggests it vanishes like 1/N .
In principle δ N (a) can be computed explicitly, we have the analytical expression of both T N and M N , however upon first inspection it seems to be too complicated to handle directly analytically. On the other hand and interestingly our numerical studies of several Toeplitz matrices similar to T N but with different c k and α ± , showed that our ansatz extends beyond our particular problem. More precisely we have.
If T N is a Toeplitz matrix given by a k e i2πkl/N a k and b k are two functions of k, which need to satisfy certain technical conditions otherwise arbitrary.
Then the following matrix This results, if confirmed by analytic calculation, will offer a general solutions for the asymptotic inverse of a wide class of Toeplitz matrices, with N -dependent symbol, which are not covered by the known results in literature [27,28]. Therefore a more mathematically oriented work is needed to settle this issue.
Now, let us test our conjecture by comparing the trace of M N for N → ∞ with the asymptotic result obtained using GLT machinery, namely equation (3.23).
Using equation (3.24) is straightforward to show In the limit N → ∞ we can turn the above sums over p and k into complex integrations around the unit circle and show This result is proved in Appendix II. The result we have obtained using our ansatz is identical to the asymptotic limit obtained using GLT machinery, equation (3.23). Now, although the above result for the trace does not consist a proof for our conjecture it however provides, along with the numerical results, strong supporting evidences for the its truthfulness. A complete proof of the above conjecture and its use to extract the universal divergence for the Half-space EE in 1+1 will be the subject of a subsequent work.

Infinitismal Intrerval: One point Entanglement
One of the interesting cases which has been discussed extensively in literature is EE resulting from tracing out a single finite interval in 1+1. The universal part of the EE in 1+1 has been derived using different techniques, see [4,5] and references therein. This leading universal divergences is given for a massive scalar theory by L is the length of the interval and ǫ is the UV cutoff of the theroy. The above universal leading divergent term of EE is generally derived by assuming L ≫ ǫ. When the ignored region is of the order of the UV cutoff the subsystem begins to fill most of the universe, thus there is less information to be lost by not measuring outside, and eqt (3.29) suggests that EE becomes of O(1), which would be dependent on the details of the regularization and looses physical meaning . However, with more elaborated tools it was subsequently possible to show the existence of an extra IR divergence (softer) in the massless limit, see for instance [5] .
More precisely we have for m → 0 the following behavior for the entropy difference which suggests an IR divergence given by S(L) = 1 2 ln(− ln(m)) for any set. This IR divergence was given a simple heuristic explanation in terms of the zero mode of the quantum field in the massless limit. As a consequence the mutual information between two sets A, B , I(A, B) ∼ 1 2 ln(− ln(m)) is also IR divergent [5]. According to these results one would expect the entropy of an interval with size comparable to the UV cutoff to be also divergent with a universal coefficient instead of being of the order of one as it was first anticipated, and this is what we are going to show using the approach developed in this paper.
For purely technical reasons we shall slightly modify our lattice regularization and add one point in the middle of the lattice considered previously, the origin; to obtain lattice with 2N +1 point . We impose periodic boundary conditions, the size of the lattice becomes R = (2N + 1)ǫ, with 2N + 1 point.
Consider the entanglement entropy resulting from tracing out the set A made of the point q 0 .
Writing Trρ q 0 using Euclidean path integral , equation (2.6), the corresponding permutation matrix P π can easily be shown to be given by Where 0 2N stands for (2N +1)×(2N +1) zero matrix, P 2N and P 0 are (2N +1)×(2N +1) projectors matrices on the available and unavailable regions defined as follows : 4 . P 0 projects onto the subset A = q 0 and P 2N projects onto the complement made of the 2N points.
The form of the permutation matrix suggests the following decomposition of the potential matrix V Where A and B are N × N matrices, b and b ′ are column vectors and a is the same parameter defined in the previous subsection, a = µ 2 +E 2 + 1. A similar decompositions for V −1 − and W − follow The explicite forms of all these matrices and vectors follow from eqt (3.4) . The next step is to evaluate V + −1 and W + .
Using the explicite form of the permutation matrix it is easy to show that Putting everything together the following product can be readily computed The resulting matrix is block circulant. On the other hand the matrix W − + W + is block circulant as well : This allows us to block diagonalize the following matrices product ( Evaluating the above trace requires determining the inverse of the matrix T given by the above equation.
Note first that this matrix can be rewritten as Where ζ = 1 + e 2πij n . Second W is circulant invertible matrix and both M f and M T f are of rank-one, therefore by recursive use of Miller theorem [29] we can explicitly invert T : Let us note that we so far have not assumed anything about N and all these results are valid for arbitrary N . Now it is a matter of a length but straightforward calculation to compute the trace We note here that it is already possible at this stage to extract the divergent piece of the entanglement entropy in the limit N → ∞, µ → 0 by expanding the integrand in the above equation around E = 0 and using the asymptotic behavior of W N N , W −1 N N and V −3/2 N N (see below). However, it turns out that we can do better than this and obtain an exact compact expression for EE of one-point in 1 + 1.
To that end it is convenient to introduce the following change of variable ω(E) = W N N .W −1 N N . Then eqt (3.45) becomes It is interesting to note to note that equation (3.46) is exactly what one should have expected to obtain from the beginning. The final result we have obtained is only a rediscovery of the formula obtained in [9] for the case of one point entanglement in terms of two-point correlator functions. To bring this out more explicitly we note that the relevant correlation functions of the field variable and its conjugate are given in this case by It therefore follows that C = √ X N N P N N . Consider now the limit N → ∞.
In the continuum limit µ = mǫ → 0 and using the asymptotic form of P 1/2 and P −1/2 we find C = 2 π (− ln µ) 1/2 and we finally obtain ln Trρ n A = (n − 1)( 1 2 ln(− ln mǫ) + ln 2 π ) + ln n (3.50) Therefore our calculation confirms the universal character of the divergence already observed in the finite interval case, however the origin of divergence in our case is slightly different although not fundamentally. Whereas in the finite interval case the divergence shows up in the massless limit, m → 0, in our case it arises as the size of the region shrinks to zero, becoming of the order of the UV cutoff. Nevertheless, the divergences in both cases has the same heuristic interpretation given in [5], it is due to the development of zero modes in the limit µ → 0, when taking ǫ to zero, as can easily be seen from the potential matrix and eigenvalues, eqts (3.3) and (3.4), which in turns leads to a logarithmic divergence in the two-point correlation function, W −1/2 N N | E=0 . This coupling between the IR behavior or the zero mode and the UV cutoff in the massless limit of 1 + 1 scalar theory is a distinctive property that is worth discussing in this connection.
We first note that while the UV divergences of EE in 1 + 1 scalar theory are widely discussed, IR diveregnces are much less studied. A focused numerical discussion of the zero modes and IR diveregences for 1 + 1 scalar theory can be found in [30], see also [31] for more general and related discussion .
Second, it is well known that 1 + 1 massless scalar theory has no normalizable ground state in the full Minkowski space nor in R × S 1 -like the one we considered in our paper by imposing periodic boundary conditions. Therefore if we start with a massive scalar theory in Minkowski space or with its lattice version, impose periodic boundary condition ( R × S 1 ), compute the EE and take the massless limit we should not be surprised to get IR divergences. However what seems to be more interesting is the consequences that this IR divergences has on the UV scaling of EE.
This coupling between the UV and IR cutoffs can be traced back to what may be referred to as the IR nonlocality that accompanies massless field theories in 1 + 1 dim [31]. This IR nonlocality has two related symptoms. The Green function of a massless scalar field grows logarithmically rather than falling off with separation, and the absence of a normalized ground state in both full Minckowski space and R × S 1 .
The concequenses of this coupling between the two scales of the theory can now be used to heuristically understand how the double-logarithmic divergence we obtained in (3.50) emerges as a concequence of the developement of zero mode in the limit ǫ → 0 instead of m → 0.
If we tarce out an interval of order of UV cutoff, one lattice spacing say, in this case one would expect the usual universal logarithmic divergent term, 1 3 ln L ǫ , to become of order one, whereas the double-logarithmic becomes instead a leading term. Of course the this double-logarithmic term blows up not as a consequence of m → 0, which we kept finite, but as we mentioned above it is rather due to the fact that the size of the interval has become of the order of ǫ.
To see how this is related to the IR divergence of the theory we note that for an infintesimal interval and based on dimensional grounds we are left with two scales at our disposal in the theory, the mass m or the IR cutoff, and ǫ or L (but not both L and ǫ, the length of the interval is no longer an independent scale as it is of order of ǫ). Therefore the entropy will be controlled by the dimensionless parameter mǫ = µ ∼ mL, and it is this coupling between m and L or ǫ which brings in the same divergence we see in the massless limit. Technically speaking there is no way to distinguish between taking m to zero or ǫ to zero. The situation is somehow similar to the case of an infinite half line in 1+1, in which the theory only has two scales and they are coupled to give one dimensionless parameter µ = mǫ. The potential, or more precisely the rescaled potential, depends only on this dimensionless parameter and it develops zero modes whether we take m to zero or ǫ to zero. The EE is a function of mǫ, and the logarithmic scaling of the entropy, ln mǫ, can intuitively be related to the IR behavior of the massless theory and give a heuretsic explanation to why the area law in 1 + 1 is logarithmic rather than being just a constant as one would have naively expected, see for instance [31] . Now, if we instead considered an interval made of finite number of points, the size of the interval would be still infinitesimal and we would obtain an EE with the following behavior ln Trρ n A ∼ (n − 1) where L = pǫ, p the number of points making the interval. We shall not try to expand the discussion about the physical relevance of considering the entanglement of an infinitesimal interval, which may loosely speaking be seen as simulating a black hole in its final faith, in one 1+1, because no one would expect our picture of space-time and field theory will continue to make sense when the size of the black hole becomes of the order of the Planck scale say. The relevance of our result is limited to showing the universal character of a softer IR/UV divergence using different regularization scheme and different technique .
Before concluding this section, there is the question of how the current approach can be used to compute EE for a finite interval in 1+1. We shall not try to tackle this problem, however we believe that the case of one point already sheds some light on the technical challenge that one would face. The main mathematical difficulties is finding the inverse of (W − + W + ) −1 . Although partial (block) inversion seems to be always possible, the remaining N × N matrices are generally not easy to invert. In the one-point case it turned out that one has to invert a matrix given by the sum of a circulant matrix and two rank-one matrices, and therefore it was easy to invert the full matrix using Miller theorem. Now, it is not difficult to see that if we considered an interval made of p points say, the resulting matrix would be a sum of circulant matrix and two (may be more) rank-p matrices and here too Miller theorem can be efficiently used, as p will be much smaller than N ( p ≪ N , or L ≪ R ). The addition of new points, p points say, will increase the entropy; for the interval to be of a finite size p should be of the order of 1 ǫ and this is how one would expect the universal divergence ∼ ln p = ln L ǫ to build up.

Interacting Theories, Beyond The Gaussian State
As emphasized in the introduction one of the main motivation for developing Green's function approach for Rényi entropy for models with finite number of degrees of freedom ( lattice or fuzzy models) is to investigate the effects of interaction on the behavior of EE. Indeed free field theories and Gaussian models are generally an approximate idealization of more realistic models for nature. In this section we shall show how the present technique is used to systematically compute the correction due to interaction (quartic type) order by order in perturbation theory. Let us consider λφ 4 type interaction added to our starting lagrangian of equation (2.1). We first note that this quartic interaction can be introduced in two different natural ways. The first is a lattice discretization of φ 4 theory, the lagrangian (3.1), which leads to an interaction term of the form q 4 i . The second possible form is λ(Q T Q) 2 , which can not be interpreted as resulting from the discretization of a local interaction term as we shall show below. The two different interaction terms will be shown to lead to drastically different results.

Local Interaction
Let us start by the following lagrangian This is just a lattice discretization of the 1 + 1 interacting scalar theory .
In view of the fact that the interaction term is invariant under the action of the permutation matrix, a perturbation series for ln Z n (λ) is straightforward to write down and formally similar to the calculation of [14] in the continuum. ln Trρ n A (λ) = ln Z n (λ) − n ln Z 1 (λ) (4.2) A generating partition function Z n (J) can easily be introduced to establish Wick's theorem for this formalism and obtain the following leading order correction to Rényi entropy 5 As an explicit example let us reconsider the case of 2-coupled H.O and compute the first order correction to the entanglement entropy due to an interaction quartic term of the local type.
Considering the correction due to quartic interaction for 2-couple H.O should furnish a good pedagogical exercise for developing a perturbation theory that allows one to compute the effect of interaction order by order; however in view of the limitation of the real time approach to Gaussian models this problem has to our knowledge never been tackled systematically. On the other hand we shall use this interacting model ( non-gaussian) to address the issue of the space-time nature of the EE. This issue was raised by Sorkin in the context of continuum field theory in order to free the entanglement entropy from reference to a density matrix localized to hypersurface [32] . A covariant definition of spacetime entropy ( which includes entanglement entropy) was derived using non other than the Peirls bracket and the related spacetime two-point correlation functions. This covariant definition was limited to gaussian states or free field theories. Related to this definition there is the well known result for bosinic and fermionic lattice systems, where it was shown that reduced density operator can be determined from the properties of the two-point correlation functions for gaussian models [9].
Motivated by the work of Sorkin the authors of [18] tried to generalize the results of [32] and investigated the existence of a covariant spacetime definition of entropy for nongaussian states or interacting theories. Therefore we shall start by briefly reviewing their main results and put them in the context of the model that we will consider below.
Consider first a generic gaussian state ρ(q, q ′ ) for one degree of freedom 6 , The the different correlators can be computed using this density.
using S = lim n→1 ∂ n Trρ n equation (4.5) is recovered. 6 The field theory problem can be divided into a series of calculations each involoving a single degree of freedom [32].
In [18] a conjecture regarding the behavior of the entanglement entropy to the first order in perturbation theory was put forward and justified by explicit calculation for a particular non-gaussian state.
Consider for instance the following non-gaussian state N is a normalization constant, λ is the interaction constant, the perturbation parameter.
We first mention that we shall show exlicitly using the Green's function approach that this non-gaussian state does actually simulate the reduced density operator for 2-coupled H.O with quartic interaction.
Let us conjecture that to the first order in λ the entropy formula in the Gaussian case, eqt(4.5), keeps the same form but with the correlators are now evaluated with respect to non-guassian state ρ λ .
Those correlators were computed in [18] and the modified ξ λ was shown to be given by . On the other hand we can compute the entanglement entropy directly by the replica trick using ρ λ .
This gave to the first order in λ the following entropy It is now easy to show that the above expression of the entropy is none other than the expansion of the following formula to the first order in λ Therefore we reach the conclusion that EE standard formula (4.5) for Gaussian states holds to the first order in perturbation theory beyond the Gaussian theory and EE is still captured by the two-point correlators (the perturbed ones).
Of course the above procedure does not for instance allow one to compute systematically the correction to EE for interacting theory. It only shows that to the first order in perturbation theory the EE is given by the same formula valide for the Gaussian state but with the correlators replaced by their non-gaussian counterparts. Actually this was shown in [18] to be a general future of any perturbation to the first order.
In what follows we shall work out an explicite example, generic enough, of interacting theory with quartic interaction and compute from the first principle the first order correction to EE and compare our results with that of [18].
The lagrangian is given by The gaussian part of this Lagrangian is basically the model we have previously seen, equation (3.13), with matrix potential given by Its free, λ = 0, reduced density matrix ρ(q, q ′ , λ = 0) is of the from given by equation (4.4). For λ = 0 the density matrix is unknown, however within our approach based on Green's function we do not need to construct the reduced density matrix explicitly to compute the correction to EE.
According to (4.3) the correction to Rényi entropy is given to the first order in perturbation theory by In appendix II we show that B ± can be summed using standard summation formula. λ ± can be expressed using ξ as defined in the previous section, equations (3.16) and (3.17).
Using the definition of ξ it is a straightforward exerecise to express the leading order correction to the Reyni entropy δ ln Trρ n A (λ) as where B n = 1+ξ n 1−ξ n , x = B 1 = 1+ξ 1−ξ . As can be easily checked the correction vanishes for n = 1 and for ξ = 0 or x = 1. The vanishing of δ ln Trρ n A (λ) for ξ = 0 reflects the fact that quartic interaction chosen induces no new entanglement and the correction is due to the perturbation of the ground state.
The correction to the entanglement entropy is obtained by the replica trick This leading order correction to the entanglement entropy is very particular and confirms the results we discussed above about the covariant spacetime definition of EE and the validity of Gaussian EE formula beyond the Gaussian states or for interacting theories to the first order [18].
This can be explicitly brought out by noting that δS ent can be obtained from the first order expansion of the following entropy expression Or for ln Trρ n A we have Therefore to the first order in perturbation theory EE is given by the same expression as EE at the free level with the parameter ξ being replaced by its perturbed counter part evaluated using the pertubed two-point correlators.
In addition we can use our result to deduce the corresponding non-guassian reduced matrix ρ λ (q, q ′ ) to the first order by equating ξ λ of equation (4.17) with ξ corr of equation (4.8) and hence fixing α 1 , α 2 , α 3 appearing as arbitrary parameters in (4.7).
The natural question that arises now is what about higher orders? Do the two points correlators continue to capture EE and does the Gaussian standard formula for the entropy continue to be valid beyond the first order of perturbation theory?
According to a general argument presented in [18] higher order correlators start to contribute from the second order on. Therefore the validity of the Guassian formula for interacting theory seems to be a particular future of first order perturbation theory, possibly due to the fact that the tunneling terms of the Green's function, G ±∓ , do not contribute at this order. Nevertheless it would be interesting to push the calculation beyond the first order using our approach. We hope to return to this point in the future.

Non Local Interaction
In this subsection we consider another type of quartic interaction which looks natural to be added to the free lagrangian of a set of coupled H.O. This interaction term is of the form λ 4 (Q T Q) 2 . Using Wick's theorem it is straightforward to show that the first order correction in λ is given by In this correction we can notice that there are two terms, the first comes from what we may refer to as connected Feynman diagrams and can be shown to be finite and formally equal to where F is defined by The second contribution 7 (TrG(τ, τ )) 2 − n(TrG(τ, τ )) 2 ), which can be interpreted as coming from disconnected Feynman diagrams, has unremovable IR divergence coming from the integration over τ . More explicitly we have The first and second term on the right hand side of the above equation give finite contribution. However the last two terms do not cancel each other and lead to an IR unremovable infinity. Indeed unlike what happened in the previous interaction type, q 4 i , or in the Gaussian case , where this kind of IR divergence is cancelled out by the normalization term (−n ln Z), here the two terms scale differently with n, namely we have, This IR divergence may seem surprising or unexpected. Actually it can be shown that the contribution of this kind of disconnected diagrams can be summed to all order . The natural question is what is peculiar about this interaction term?. To try to answer this question we go back to the interaction term (Q T Q) 2 . If we expand this term we easily notice that it is a highly non-local, each oscillator couples to all other oscillators in the chain, and hence can not be understood as arising from the discretization of a local interaction term of scalar theory, in other words it is a long-range interaction and induces couplings between different socillators at the interaction level. Moreover this interaction or this model can not emerge from the discretization of a Lorentz invariant contiunuum theory.
It is important to note here that at the Guaussian level, quadratic lagrangian, longrange interaction or non-local terms in continuum field theories do not seem to have this sort of IR problem. Long-range peotenials or non-local theories only change the scaling of the EE as one would expect, see [33] for an example of non-local field theory and references to earlier works. The absence of this kind of IR divergence at the Guassian level can in fact be seen from our main formula for the EE given by equation (2.27), which is valide for any positive definit potential, whether is a long-range or short -range. However, when it comes to interactions our result strongly suggests that EE does not tolerate non-local interactions, at least not any type of non-local interaction, or interaction which induces higher order coupling between different oscillators. In other words and technically speaking, it is the fact that the chosen interaction induces new entanglement, terms like q 2 i q 2 j , which is the source of the problem upon the following Wick contractions : This type of contraction leads to an IR divergence due to the different scalings of the n-fold Green's matrix function and n = 1 Green function after taking the trace, unlike local interaction where the oscillators degrees of freedom are not coupled by the interaction and no such contraction or disconnected diagrams would be present. It is important here to note that by local interaction we refer only to the type considered in the previous section, i.e q 4 i , and this leads us to the following natural question. Is this IR divergence associated with the non-locality of the quartic interaction, i.e. long-range, or is it a symptom that will accompany any interaction which couples different oscillators, even if it is a short-range (nearest neighbor quartic interaction) 8 Although we do not have a definite general answer for this question, we believe that it is the case for all one dimensional chains with quartic interactions which couple different oscillators. However, for other non-local models such as 2 + 1 fuzzy field theories the story seems a bit different. An outlook on the EE for interacting fuzzy scalar theories is the subject of the next and last section of this paper

Applications to Fuzzy Spaces and Outlook
The main motivation behind the development of this discrete Green's function approach to EE is to investigate EE on fuzzy spaces, in particular the effect of interaction in connection with IR-UV mixing phenomena on EE. Indeed in view of the IR divergence we met when we considered a particular quartic interaction and the phenomena of UV-IR mixing, it is not a priori obvious how the entanglement entropy will respond to a natural quartic interaction on fuzzy or non-commutative underlying geometry. Although this issue is the project of future work and shall not be tackled in the present work, we are going to briefly outline how field theories on fuzzy spaces fit into this framework and the possible challenges to analytically investigate them . To that end we consider scalar theory on 1 + 2 space-time, with space modeled by a fuzzy sphere.
The entanglement entropy on fuzzy sphere was first considered in [19,20] and later by several authors [34,35] .
The starting Lagrangian in [19] was The scalar field φ is an N ×N hermitian matrix with mass parameter µ. The Laplacian L 2 i is the SU (2) Casimir operator given by The L i generate the SU (2) irreducible representation of spin l = N −1 2 . By using the canonical basis where L 3 is diagonal and choosing suitable field variables it was possible to show that the above free scalar theory on a fuzzy sphere into splits naturally into 2(2l) + 1 independent sectors {H m }, m = −(N − 1), · · ··, (N − 1), each sector H m has N − |m| degrees of freedom (N − |m| coupled harmonic oscillator) and described by a Lagrangian L (m) . More explicitly we have and B a = a(N − a) and A a = −a + N +1 2 , a, b = 1 . . . N − |m|. The new canonical variables Q (m) a are introduced after decomposing φ as φ = φ 1 + iφ 2 and defining Φ = φ 1 + φ 2 .
In order to introduce the entanglement entropy for these fuzzy models one needs first to divide the field degrees of freedom into two sets residing in two, or several, disjoint regions. In [19] it was conjectured using heuristic arguments that the field variables represented by the elements of the upper left triangular part of the field matrix Φ + and the right lower triangular part Φ − correspond respectively to the upper and lower fuzzy hemispheres. This conjecture was subsequently proven in [34] using rigourous mathematical argument.
More explicitly for N = 5 we have In terms of the canonical field variables Q m tracing out the upper or the lower hemifuzzy sphere corresponds to taking each individual sector H m and tracing out half of the degrees of freedom in each sector.
In view of the fact that the different sectors do not mix for the free theory, the computation of the Rényi entropy reduces to computing the EE in each sector selectly. Each sector has a Lagrangian of the standard form (2.1). The global ground state operator is given by (5.6) and the reduced one ρ + = T r − ρ is given by where ρ (m) + is obtained by tracing out the degrees of freedom corresponding to the ones in the lower hemi-fuzzy sphere in each sector.
In [19] it was shown numerically and using the Real Time approach that the resulting EE is proportional to the number of freedom residing at the boundary. The result could as well be interpreted as proportional to the circumference of the separating region, a fuzzy circle.
The present Green's function technique can similarly be applied to compute Rényi entropy (at least numerically). Indeed, in view of the canonical splitting of the Hilbert space of the free theory, the Green's function matrix for this model is given by Unlike the case of 1 + 1 lattice the potentials V (m) are not circulant matrices, which poses an additional complication to the analytical investigation of this model, and for the fuzzy disc as well [19]. However, we note that it was shown in [20] that the EE is essentially captured by the near boundary approximation where the potentials are Toeplitz matrices. In addition to this it can be shown that not all sectors, different m, are equally relevant, there are sectors which are weakly entangled and hence give negligible contributions to the entanglement entropy in the large N limit [36], and the relevant sectors are governed by Toeplitz potentials in the near boundary approximation . Putting all this together it is conceivable that in the large N limit one can approximate the potentials for the relevant sectors with circulant matrices which would greatly simplify the problem. Now, as we mentioned this technique is specially devised to go beyond the free gaussian theories and investigate the effects of interaction on the EE on fuzzy spaces. Therefore we conclude this paper by outlining how this could be done.
Let us assume that we have added a quartic interaction term to our starting lagarangian (5.1), In order to apply the Green's function formalism the interaction term should expressed using the Φ or the new canonical variables Q (m) . We shall not go into the detail of this, however we shall mention few important points.
First it is not difficult to see that this interaction term mixes different sectors and introduces new couplings between the degrees of freedom which were not coupled at the free level. Second and related to the first the interaction term once expressed using the n-fold variables is not invariant under the action of the permutation matrix P π within each sector. Therefore the application of Wick's theorem is a bit subtle.
Some natural questions arise here and make the study of the above interaction worth considering. The first is how such interaction would enhance the gaussian entanglement in view of the UV-IR mixing phenomean related to both the interaction and the entanglement entropy in the large N limit. Actually divergences of UV-IR mixing origin had already been observed in [19] at the free level in 2+1 . A second question concerns a possible unremovable IR infinity, for finite N , similar to the one we met in the previous section for the non-local interaction.
As for the second question, i.e. IR divergence, we do not expect such unremovable IR divergence to be present in this model. Although the interaction term once expressed using the above canonical variables, eqt (5.4), seems to behave like the non-local term we considered in the previous section, this behavior is partially apparent. For instance the mixing of different sectors which may cause this infinity is similar to the mixing of different Fourier modes for any local interaction in standard continuum field theory, such a mixing do not show up at the free level. Nevertheless the interaction term (5.9) can be seen to introduce quartic coupling between different oscillators DF within the same sector ( the same m) and therefore only explicite calculation can completely settle this issue. Actually working the case N = 2 would be enough, because as long as this particular IR problem is concerned the size of the fuzzy sphere should play no role. Now, in order to carry out an analytical study of the EE for free and interacting fuzzy field theories, we believe that some mathematical tools should first be developed, in particular understanding the asymptotic behavior of Toeplitz matrices with N -dependent index on which very little if any is known.
Finally, we note that it is conceivable that one may have to take a different path, based on a different decomposition of the field variable, in order to analytically tackle the fuzzy regularization using the Green's function technique. For instance we may reformulate the Green's function developed in the first section by re-expressing the Euclidean path integral formula for the reduced density matrix using a doublet field notation as follows.
We let Φ = Φ + Φ − . With this doublet notation the lagrangian (5.1) takes the form where L ±± and L ±∓ follow from the adjoint actions of the angular momentum operators, viz, ..etc This formally will reduce the problem to a 2-dimensional one with reduced Laplacians and one may hope to obtain a Green's function matrix that makes the problem more tractable analytically. As a by product of the formulation it is possible that one would end up deriving the Laplacian of n-fold cover of a fuzzy space and understanding the constructions of new class of fuzzy spaces . We hope to come to this and other related issues in future work.
6 Appendix I

Series Expansion of Green's Function
This section of the appendix is devoted to proof of the expansion of the Green's matrix function.
Let G(x, x ′ ) be a n × n matrix valued function which satisfies following equation V and V p are two real positive definite n × n matrices, and generally noncommuting. The domain of the solution is R.
To construct G we shall use Weiner-Hopf method to turn the problem into inhomogeneous Hilbert problem and show that it can be obtained order by order around the V (x) = 0 .
We first reformulate (6.1) and convert it into an integral equation.
To that end we introduce the free Green's function G 0 satisfying Using G 0 it is easy to show that G(x, y) satisfies the following integral equation By iteration we obtain the following expansion for G with G n (x, y) given by G n (x, y) = (−1) n dy 1 .. dy n G 0 (x, y 1 )G 0 (y 1 , y 2 )...G 0 (y n−1 , y n )G 0 (y n , y)V (y 1 )...V (y n ) From which the following recurrence relation follows The ordering of all quantities appearing in the above equations and in what follows is crucial, as they are generally noncommuting.
Let us construct the following scaled matrix function Using (6.3) and (6.4) we can write (6.2) as By taking the two sided Laplace Transform of both sides of the above equation and using the fact that G 0 (x, y) depends only on the difference between x and y, the convolution theorem leads to the following recurrence relation between a sequence of inhomogeneous Hilbert problems, where F ± n (s, y) and F 0 (s) are defined by Here we note that in view of the definition of G ± n (x, y) we can infer that F − n (s) is a regular function ( analytic) for Re(s) > 0 and F + n (s) regular for Re(s) < 0. For F 0 (s) we have Now, equation (6.6) is a discontinuity equation , the branch cut being the imaginary axis Re(s) = 0,we can use standard theorems known for this kind of inhomogeneous Hilbert problem, to solve these sequence order by order to obtain F n (s), from which one is computed G n (x, y) by inverting the Laplace transform.
We shall work out explicitly the region y < 0, x < 0 up to the second order. The other regions are either similar or can be deduced simply from this region. The second order turned out to be enough for to guess the general solution.
• First Order We have first to compute F ± 0 (s, y) .
it is then straightforward to show that from which we can compute F ± 1 (s, y) We note that in general one has to add an arbitrary polynomial of s ( matrix elements) to the above solution for F + 1 , however they are excluded by the requirement that lim s→−∞ F + 1 = 0. c 1 and c 2 can be easily computed using contour integration.
Re(s) < 0 For F − 1 we can use the same method or use directly the discontinuity equation to obtain for F + Now, we can compute G 1 (x, y) for x, y < 0 using the inverse Laplace transform and obtain The first term (1) has two contributions To compute (2) and (3) we sum them and use the relation between F + 0 and F − 0 to get y) can then be computed easily by complex contour integration, with x, y < 0 Collecting G 1 and G 2 it is trivial to recover the expansion of G given by equation (2.23).

Continuity Conditions
The continuity of the eigenfunctions and their first derivatives at the cut τ = 0 lead to the following continuity conditions for the Green's matrix function and its first derivatives.
We have just listed half of the set of the resulting continuity equations, namely with respect to one argument( the first or the second), because the others follow from the hermiticity condition which we shall prove lastly.
By exactly similar matrix manipulation we show.
Consider now the continuity of the derivatives. First it is easy to show the following equalities.
We shall work out one conditions explicitly, the others are very similar as can easily be noticed.
Therefore we need to show that where we have used the famous matrix identity to go from the second to the third line.

Appendix II
The aim is to compute the following limit It will be shown that I − is real in the limit N → ∞ and therefore I − and I + give equal contributions.
We start by evaluating I 0 .
Using complex contour integration, it is easy to convert the above contour integral into the following integral where R x = (r + − x)(r − − x)/2x . r + , r − are in addition to z = 0 the branch points of the function √ R z = 2az−z 2 −1
As can be noticed in converting the sum of p into an integral we have added another term proportional to z l−1−N which comes from the reminder when using Euler-Maclaurin, this term can not be neglected within this sum and actually it is this term which survives at the end. Now, using again complex contour integration it is easy to show that which is easily seen to vanish in the limit N → ∞, in view of the fact that x < r − < 1. Therefore we are left only with the second term of equation (7.5). This term after summing over l gives Where we have neglected all terms which vanish in the limit N → ∞.
Where u = e − √ e 2 − 1 is a pole beside x lying inside the unit circle. Putting everything together we find lim N →∞ It is now a matter of trvial exercise to show using Integral Tables that lim N →∞ TrM N (α) N = 1 π 1 √ a + 1 .K 2 a + 1 (7.8) We note that in all the above calculation we assumed α + = α − = 1. This particular case has to be worked out separably although the final result does not depend on the particular values of α + and α − , as long as they sum up to 2.

Appendix III
This appendix is to dedicated to the proof of equation (4.13) of section.
The main task is to compute the diagonal elements of the matrix G, in order to compute the first correction expression : Let us note that all quantities here are evaluated at E = 0.
For (τ < 0) the Green's matrix function is given by : It turns out that it is technically more convenient to write G −− (τ, τ ) as First it is easy to show that With λ ± = √ a ± 1 and similarly for The less straightforward piece is the last term of eqt (8.3). To compute its diagonal elements we first introduce the following orthogonal 2 × 2 matrix and its n-fold version.