An Inverse Mass Expansion for the Mutual Information in Free Scalar QFT at Finite Temperature

We study the entanglement entropy and the mutual information in coupled harmonic systems at finite temperature. Interestingly, we find that the mutual information does not vanish at infinite temperature, but it rather reaches a specific finite value, which can be attributed to classical correlations solely. We further obtain high and low temperature expansions for both quantities. Then, we extend the analysis performed in the seminal paper by Srednicki [1] for free real scalar field theories in Minkowski space-time in 3 + 1 dimensions at a thermal state. We find that the mutual information obeys an area law, similar to that obeyed by the entanglement entropy at vanishing temperature. The coefficient of this area law does not vanish at infinite temperature. Then, we calculate this coefficient perturbatively in an 1/μ expansion, where μ is the mass of the scalar field. Finally, we study the high and low temperature behaviour of the area law term. 1 ar X iv :1 90 7. 08 50 8v 1 [ he pth ] 1 9 Ju l 2 01 9


Introduction
Quantum entanglement is a fundamental phenomenon without classical analogue, which plays an important role in quantum physics. When a complex quantum system lies in an entangled state, even if this is a pure state, there is no answer to the question "what is the state that describes the subsystem A?". However, the latter can be described by a density matrix, the reduced density matrix, which can be derived from the density matrix of the composite system, via the tracing out of the degrees of freedom of the subsystem A C , which is complementary to A, (1.1) Under the assumption that the quantum composite system lies in a pure state, at the limit the subsystems A and A C become disentangled, the reduced density matrix ρ A corresponds to a pure state, and, thus, the question "what is the state of the subsystem A?" acquires an answer. Therefore, it is natural to claim that entanglement is encoded in the spectrum of the reduced density matrix, and the most natural measure of entanglement is Shannon entropy applied to the latter, i.e., which is the so called entanglement entropy. Entanglement is a property that depends on the specific separation of the composite system to the pair of complementary subsystems A and A C . Naturally, one would require from a measure of entanglement to obey the property which can indeed be shown to hold, when the composite system lies in a pure state. In a seminal paper [1], Srednicki showed that entanglement entropy has a particularly interesting property in free massless scalar quantum field theory: Assuming that the system lies at its ground state, and separating the degrees of freedom to two subsystems, one containing the degrees of freedom inside a given sphere of radius R, and another being its complementary, it was shown that entanglement entropy is proportional to the area of the sphere. This property is somehow expected from the physics of entanglement: As already mentioned, entanglement characterizes the separation of the composite system to two subsystems and not the subsystems themselves. Thus, the entanglement entropy cannot depend on the properties of any of the two subsystems (such as the volume of subsystem A), but on those of their only common feature, i.e. their boundary. Nevertheless, this finding is highly intriguing, since it resembles the area law of the black hole entropy. This similarity motivates the investigation of whether the black hole entropy can be attributed completely or partially to quantum entanglement. Bombelli et.al. [2] also motivated by the similarity to black hole physics calculated the entanglement entropy for a scalar field in the background of a Schwarzschild black hole resulting in similar conclusions.
However, the entanglement entropy is a good measure for entanglement, or more generally of correlations between the subsystems, only when the composite system lies in a pure state. If this is not the case, the entanglement entropy will inherit contributions that originate from the classical entropy of the composite system, and, thus, they do not characterize the entanglement between the two subsystems. In general, when the composite system lies in a mixed state, (1. 4) In field theory, the above argument implies that when the composite system lies in a thermal state, the entanglement entropy will have contributions originating from the thermal entropy of the composite system, and, thus, will be proportional to the volume of the subsystem. Entanglement in field theory at finite temperature has been studied mainly in the context of two-dimensional conformal field theory [3][4][5][6] with the use of the replica trick [7,8]. Much fewer works focus on gapped systems [9] or to higher dimensional theories [10][11][12]. In more recent years, entanglement in thermal states has also been studied through the holographic duality. The issue has been posted in the original works that established the Ryu-Takayanagi conjecture [13,14], where the problem of the non-symmetry property of the entanglement entropy is resolved by the existence of more than one minimal surfaces, due to the presence of the black hole, which are homologous to complementary boundary region. This study has been extended in several works (see e.g. [15]). Most of these focus mainly on the geometry of the BTZ black hole [16][17][18][19], which is also relevant to two-dimensional CFTs, as this is the only black hole geometry where minimal surfaces can be expressed analytically. Entanglement in harmonic lattice systems at finite temperature has been studied in [20]. However, there is not much attention to the study of entanglement in field theory at finite temperature via the techniques originally used in [1]. When the composite system lies in a mixed state, a better measure of the correlation between the two subsystems is the mutual information, which has the symmetric property by construction. Following the arguments above, the mutual information should characterize the separation of the composite system to two subsystems and, thus, in field theory it should depend only on the properties of the entangling surface, even at mixed, e.g. thermal, states. It has been shown that in lattice spin systems the mutual information obeys an area law bound [21]. It was recently shown [22] that appropriate generalization of the techniques of [1] can be used to calculate the mutual information in free scalar field theory at finite temperature and indeed it is proportional to the area of the entangling surface.
In [23], the authors developed a perturbative approach in order to study the area law of the entanglement entropy in scalar field theory at its ground state, analytically, bypassing the numerical part of the original calculation in [1]. In this paper, we extend this method, in order to calculate perturbatively the entanglement entropy and the mutual information in free scalar field theory at a thermal state. In section 2, we study the system of two harmonically coupled oscillators at finite temperature. In section 3 we generalize to a coupled harmonic system with an arbitrary number of degrees of freedom at a thermal state. In section 4 we develop the hopping expansion for chains of coupled oscillators, i.e. systems where only neighbouring oscillators are coupled. In section 5 we use the results of the previous sections, in order to study the entanglement entropy and the mutual information in free scalar field theory in 3+1 dimensions. Finally, in section 6 we discuss our results. There are also several appendices containing more details of the related algebra.

A Pair of Coupled Harmonic Oscillators
In order to study entanglement entropy and mutual information in free scalar field theory at finite temperature, we first study systems of coupled harmonic oscillators with a finite number of degrees of freedom. The simplest such system, which is the subject of this section, is a system of two coupled harmonic oscillators at finite temperature. The analysis closely follows the original treatment presented in [1], in the sense that it is performed in coordinate representation and presents several technical similarities. A short account of this analysis recently appeared in [22].

A Single Harmonic Oscillator at Finite Temperature
First, we would like to recall some formulas related to the problem of a single harmonic oscillator at finite temperature in coordinate representation [24], which will be useful in the following. Without loss of generality, we consider the mass of the harmonic oscillator to be equal to one, i.e. the Hamiltonian of the system is In coordinate representation, the energy eigenstates and the corresponding eigenvalues of the harmonic oscillator are ψ n (x) = 1 √ 2 n n! 4 ω π e − ωx 2 2 H n √ ωx , E n = ω n + 1 2 , where H n is the Hermite polynomial of order n. The equation (2.2) trivially implies that the density matrix describing a quantum harmonic oscillator at finite temperature T is given by As a consequence of Mehler's formula, where we defined the quantities a and b as Finally, it is a matter of simple algebra to show that the thermal entropy of the single quantum harmonic oscillator at temperature T equals (2.7) Expanding the above equation at high temperatures yields whereas expanding it at low temperature yields

Two Coupled Harmonic Oscillators
Now, let us consider a system of two coupled oscillators at finite temperature. The oscillator described by coordinate x and canonical momentum p is constituting the subsystem A, whereas the other oscillator, which obviously coincides with subsystem A C , is described by coordinate x C and canonical momentum p C . All oscillator masses are taken equal to one. The Hamiltonian of the system is When the Hamiltonian is written in terms of the canonical coordinates, it assumes the form where ω ± are the eigenfrequencies of the normal modes, namely, ω + = √ k 0 and ω − = √ k 0 + 2k 1 . The Hamiltonian (2.12) describes two decoupled oscillators, corresponding to the two normal modes of the system. It follows that the density matrix that describes the composite system at finite temperature can be trivially written as the tensor product of the thermal density matrix (2.5), for each of the two normal modes, (2.14) In order to find the reduced density matrix of the subsystem A, this density matrix has to be expressed in terms of the original coordinates x and x C prior to tracing out the A C degrees of freedom, We proceed to trace out the degree of freedom of the subsystem A C , integrating out x C . After some simple algebra we find Similarly to the ground state case analysis [1], one can show that the functions 19) are the eigenfunctions of the reduced density matrix. The respective eigenvalues are This can be expressed in terms of the physical quantities of the problem, i.e. the eigenfrequancies of the normal modes and the temperature, Then, it is straightforward to calculate the entanglement entropy, which equals As we argued in the introduction, the entanglement entropy is not a very good measure for the quantum entanglement when the overall system lies at a mixed state, like the scenario under consideration. In general, it contains a contribution from the thermal entropy of the overall system. Indeed, the entanglement entropy does not vanish at the limit k 1 → 0 as one would expect from a good measure of quantum entanglement. It rather tends to the thermal entropy of a single oscillator with eigenfrequency √ k 0 at temperature T . In the case of the two coupled oscillators that we study here, it holds that S A C = S A , due to the symmetry of the system. Therefore, the mutual information is given by, where S A is given by (2.23) and S th is obviously given by the sum of two versions of equation (2.7), one for each normal mode.

Similarity to a Single Harmonic Oscillator
One may observe that the reduced density matrix (2.16) is identical to the thermal density matrix of a single harmonic oscillator (2.5), after some appropriate identifications. More specifically, there is no experiment that can be performed to the one of the two coupled oscillators at finite temperature T that can distinguish it from a single effective harmonic oscillator with eigenfrequency equal to at an effective temperature equal to The latter is always higher than the physical temperature T . This identification obeys some obvious consistency checks. For example, at the limit k 1 → 0, the two oscillators become decoupled, each having eigenfrequency equal to √ k 0 . It follows that at this limit, the system is separable, i.e. ρ = ρ 1 ⊗ ρ 2 , and, thus, the reduced density matrix should be identical to ρ 1 , i.e. the thermal density matrix of a single harmonic oscillator with eigenfrequency √ k 0 at temperature T . Indeed, expanding ω eff and T eff around k 1 = 0 yields Similarly, at the limit T → 0, one finds the following where (2.31) Therefore, we recover correctly the ground state result [1]. At low temperatures the corrections to the zero-temperature values of ω eff and T eff are exponentially suppressed and tend to reduce the eigenfrequency of the effective oscillator, whereas they tend to increase its temperature. This expansion is an asymptotic expansion and not a usual Taylor series. This is due to the fact that the involved functions are not analytic at T = 0. The results are expressed at first order in the exponentials e − ω ± T , but one has to be careful with this kind of expansion; for example, depending on the values of ω ± , the second order term in the exponential of ω + may be a more significant contribution that the first order term in the exponential of ω − .
In a similar manner at high temperatures we find This implies that at high temperatures, the eigenfrequency of the effective oscillator tends to a finite given value, whereas the effective temperature is dominated by the physical temperature of the composite system. A very interesting question that can be posted is whether the fact that the subsystem A can be described by an effective thermal reduced density matrix can be attributed to the eigenstate thermalization hypothesis [25]. Naturally, this should not be expected, since the system under consideration is integrable.
When we consider either a thermal state or the ground state for the overall system, its density matrix is time independent. This implies that the same holds for the reduced density matrix, which describes the considered subsystem. However, the subsystem is an open system, and, thus, a time-independent state, has to be a non-trivial state that describes a system in equilibrium with its environment (not necessarily thermal).
This behaviour becomes clearer in the case of many harmonic oscillators that we are about to study in next section. There, we will analyse a system of N coupled oscillators, considering as subsystem A an arbitrary subset comprising of n oscillators. Although we are not going to discuss on the similarity of the reduced density matrix to the density matrix of a harmonic system of n oscillators at an appropriate state, the entanglement entropy is identical to the sum of the thermal entropies of n effective oscillators, each lying at a different temperature. This is consistent with the picture of a harmonic system with n degrees of freedom, where each normal mode has been heated to a different temperature. Since, the normal modes of a harmonic system do not interact, this is an equilibrium, time-independent state, which nevertheless is not thermal. It follows that the reduced system is not thermalized; actually, it is as far as possible from a thermalized state, as imposed by its integrability.
In the case of the two coupled oscillators, the considered subsystem contains a single degree of freedom, and thus, such a state is a thermal one. Thus, the fact that the reduced density matrix appears to be thermal is not a consequence of thermalization, but rather a technical coincidence due to the specific selection of the state of the overall system and the number of the degrees of freedom.

High and Low Temperature Expansions
At temperatures much higher than the system eigenfrequencies, the entanglement entropy and mutual information have asymptotic expansions of the form and respectively. Notice that the coefficients of the high temperature expansion of the mutual information do vanish when the oscillators are decoupled, i.e. k 1 → 0, as expected. Furthermore, the coefficient of the 1/T 2 term in the mutual information vanishes, which turns out to be a more general feature, as we will show in next section. Finally, it is evident that the mutual information does not vanish at infinite temperature, but rather it tends to the value It is well known that in qubit systems, the mutual information vanishes at infinite temperature. It is natural to wonder what is the underlying reason for this difference between qubits and oscillators. The answer to this seeming inconsistency is related to the dimensionality of the Hilbert space of our problem. In all qubit systems, the related Hilbert spaces are finite dimensional. Trivially, at the infinite temperature limit, the density matrices of the composite system tends to This is a separable density matrix, implying trivially that It follows that the entanglement entropies tend to whereas the thermal entropy tends to lim T →∞ The above imply that the mutual information vanishes at infinite temperature, However, in our case the corresponding Hilbert spaces are infinite dimensional and the above arguments cannot be applied equally well. Both entanglement entropies S A and S A C diverge at infinite temperature as ln T . This divergence is cancelled in the mutual information, via the same mechanism that enforces the mutual information to vanish in qubit systems; however, there is a finite remnant.
In general, the mutual information measures both classical and quantum correlations. So, another natural question concerns the origin of this mutual information remnant at infinite temperature. The mutual information I ∞ coincides with the mutual information that one can calculate via a classical analysis, as shown in the appendix A (see also [20]). Therefore, this infinite temperature remnant should be attributed solely to classical correlations. As intuitively expected, at infinite temperature the classical fluctuations completely dominate and yield the quantum fluctuations irrelevant. This is also in line to the fact that another measure of quantum entanglement, the entanglement negativity, also vanishes at infinite temperature. Actually, the negativity vanishes above a finite critical temperature, as shown in appendix B, a phenomenon widely known as sudden death of entanglement. However, this does not necessarily imply that there is really such a finite temperature phase transition in the system of coupled oscillators. The absence of negativity is not a proof of lack of entanglement in infinite dimensional Hilbert spaces, as in finite dimensional ones [26,27]. This issue requires further investigation.
At low temperatures, the entanglement entropy tends to the zero temperature result, plus exponentially suppressed corrections Similarly, the mutual information is equal to As shown in figure 1, where the mutual information is plotted as a function of the temperature, the mutual information may be a monotonous function of the temperature or not. This depends on the relevant magnitude of the couplings k 0 and k 1 , which determines the sign of the coefficient of the 1/T 4 term in the high temperature expansion.

Entanglement Entropy and Mutual Information
Building on the results of Section 2, we proceed to study a system of N coupled harmonic oscillators. In this analysis, the subsystem A C coincides with any subset of n oscillators. Without loss of generality, all oscillators are considered having unit mass. The Hamiltonian is given by The matrix K is symmetric and all its eigenvalues are positive, since the above Hamiltonian should describe an oscillatory system around a stable equilibrium. Writing down the Hamiltonian in terms of the normal coordinates y i , which are related to the initial coordinates x i via an orthogonal transformation O, yields where ω i are the frequencies of the normal modes. In other words, the orthogonal transformation O diagonalizes the matrix K, or These matrices can be related to the eigenfrequencies of the system as Since the oscillators corresponding to the normal modes are decoupled, the density matrix of the overall system can be written as the tensor product of the density matrices corresponding to each of the normal modes, We express the density matrix in terms of the original x coordinates, using the orthogonal transformation O, In the following, we use the block form notation We will also write any symmetric matrix M in block form, using the notation The indices A, B and C will always indicate the corresponding blocks of such matrices. Then, the density matrix ρ (x, x ) can be expressed as, We proceed to trace out the first n degrees of freedom to find the reduced density matrix for the remaining N − n ones. Simple algebra with Gaussian integrals yields Similarly to the ground state case [1], one may find the spectrum of the reduced density matrix, via the explicit construction of its eigenfunctions. It reads where the quantities ξ i are given by and β Di are the eigenvalues of the matrix γ −1 β. It follows that the entanglement entropy is given by Notice that this formula is identical to the formula that would provide the thermal entropy of n independent oscillators, each with eigenfrequency 1 − β 2 Di and at temperature − 1 − β 2 Di / ln ξ i . As a consistency check, let us consider the special case where the two subsystems are decoupled, i.e. K B = 0. It holds that In this case, it is straightforward that Therefore the eigenvalues β Di of the matrix γ −1 β can be expressed in terms of the eigenvalues of the matrix K C , i.e. the eigenfrequencies ω i of the decoupled subsystem A. Notice that the eigenfrequencies, as well as the thermal entropy of the subsystem A are well defined in this limit, since the two subsystems are decoupled. The eigenvalues β Di read It follows that Comparing equations (2.7) and (3.17), we conclude that in the K B = 0 case, the entanglement entropy is simply equal to the thermal entropy of the subsystem A. This is expected, since at this limit, the composite system density matrix is separable. This also implies that the mutual information vanishes at this limit.

High and Low Temperature Expansions
A high temperature expansion of the above result can be performed. The details are included in the appendix C. The high temperature expansions of the entanglement entropy and mutual information are and respectively. Interestingly, the coefficient of 1/T 2 in the high temperature expansion of the mutual information vanishes for any system. It is trivial to show that in the case of the two oscillators, where the matrices of the above formula are simply numbers, namely, K A = K C = k 0 + k 1 and K B = −k 1 , the above formulae reproduce the expansions (2.35) and (2.36). Furthermore, in the case where the two subsystems are decoupled, i.e. the matrix K B vanishes, the above formula implies that the first terms in the expansion of the mutual information are vanishing, as expected. At low temperatures, the situation is a little less transparent. As in the case of the two oscillators, the involved functions are not analytical at T = 0. Nevertheless, we may obtain an asymptotic expansion, approximating the hyperbolic functions with exponentials. It turns out that the matrix γ −1 β, whose eigenvalues determine the entanglement entropy is given in this expansion by The details of this calculation are included in the appendix D. It is not possible to obtain a generic expression for the low temperature expansion of entanglement entropy or the mutual information in this limit. However, the equation (3.27) implies that at low temperatures the corrections to the zero temperature result are exponentially suppressed as exp (−ω i /T ), where ω i are the eigenfrequencies of the overall system. In the case of the two oscillators, it can be shown that the above formula correctly reproduces the results (2.43) and (2.44).

Chains of Oscillators
In this section, we consider systems of coupled oscillators, with the specific property that only adjacent degrees of freedom are coupled. In other words, we consider a coupling matrix K of the form We will refer to such systems as "chains of oscillators". This class of harmonic systems, apart from their own interest, will be essential in the study of the free scalar quantum field theory in next section.

A Hopping Expansion
Assuming that the diagonal elements of the matrix K are much larger than the offdiagonal ones, one may follow the approach of a hopping expansion, in the spirit of [23], in order to calculate the entanglement entropy and the mutual information for this class of systems perturbatively. One may define and perform an expansion in ε (or equivalently in l/k).
In the following, we will adopt a particular notation for the matrix elements of all the involved matrices. The subscript denotes the line of the element, when it lies on top of the main diagonal and its column, when it lies below the main diagonal. The superscript denotes the diagonal (i.e. superscript 0 implies that the element lies in the main diagonal, superscript 1 implies that it lies in the first superdiagonal, superscript −1 implies that it lies in the first subdiagonal and so on). In other words . Obviously for symmetric matrices M it holds that M j i = M −j i and we will not post the results for both. Finally, the second superscript, which will appear into parentheses, denotes the order of the element in the ε expansion.
Furthermore, for simplicity we define the functions which will appear throughout the calculations of this section. Expanding the matrix γ −1 β in ε, one can show that the zeroth and first order terms are given by and and all other elements are vanishing. The second order result is given by a little more complicated expressions. We provide here only its diagonal part, which is crucial for the following There is a special contribution in the very first element, which originates from the term of the γ and β matrices. This is going to play an important role in the following. More details are provided in the appendix E.1.
The eigenvalues of the matrix γ −1 β have to be perturbatively calculated in the ε expansion. The problem is more difficult than the zero temperature problem [23]; In this case, the elements of the matrix γ −1 β obey an hierarchy in both its directions, i.e. the leading contribution to the element (γ −1 β) ij is of order i + j. This hierarchy is inherited to the eigenvalues, setting their perturbative calculation a simple task. However, in the case of finite temperature, the thermal contributions have changed this structure; The leading contribution to the element (γ −1 β) ij is of order |i − j|. It follows that a more systematic approach is required.
In order to obtain the expressions (4.9), (4.10) and (4.11), we only needed to demand that the diagonal elements of the matrix K are larger than the non-diagonal ones. However, this does not suffice for the perturbative specification of the eigenvalues of the matrix γ −1 β. In order to clarify this, we post a simple, indicative example: Assume the Hamiltonian where the diagonal elements are much larger than the off-diagonal ones. In order to calculate its eigenvalues perturbatively, naively one would consider the diagonal part of this Hamiltonian as an exactly solvable unperturbed Hamiltonian and the off-diagonal elements as a perturbation. However, this is not necessarily a good approach. This is evident in this two by two example, since the problem is simple enough to find its answer analytically, Following this approach is equivalent to Taylor expanding the above eigenvalues with respect to the parameter g. However, this expansion does not converge whenever g > In this case, one should perform a Taylor expansion in h 1 − h 2 , which implies that another setup for the perturbative calculation of the eigenvalues should have been considered. The unperturbed Hamiltonian should be considered proportional to the identity matrix. Then, there are two perturbations: one that consists of the non-diagonal part of the Hamiltonian and a manifestly smaller one, which is diagonal and proportional to the difference of the two diagonal elements. Now the unperturbed problem is degenerate and the basic eigenvectors are determined by the large perturbation.
Thus, the appropriate structure of the perturbation theory depends on the ratio of the diagonal elements to the difference of the diagonal ones. The assumption we have made for the matrix K does not determine this ratio. It follows that there are two distinct approaches in determining the eigenvalues of γ −1 β, which we will call "non-degenerate" and "degenerate" perturbation theory. They are presented in appendices E.2 and E.3, respectively.
When the diagonal elements have differences of the same order of magnitude as themselves, the non-degenerate perturbation theory applies and it yields 14) The unique second order contribution to (γ −1 β) 11 has affected a single eigenvalue at this order. This is similar to the zero temperature case; however, the other eigenvalues do not vanish. These eigenvalues imply that the entanglement entropy for subsystem A equals The first two lines of the the above expression contain the contributions from the generic eigenvalues. The rest originates from the special eigenvalue β D1 . The entanglement entropy S A C has a similar structure. The contributions to the entanglement entropy from all the generic eigenvalues are identical to those of the thermal entropy, and, thus, at this order in l/k, the mutual information receives contributions only from the two special eigenvalues, one from each subsystem. It is equal to (4.16) Expanding for high temperatures the above result yields which coincides with the l/k expansion of the high temperature formula for the generic oscillatory system (3.26).
In the case the differences of the diagonal elements are smaller than the nondiagonal ones, one should apply degenerate perturbation theory. We will focus on a subclass of this kind of problems that emerges from the discretization of 1 + 1 dimensional field theory, namely the case where the matrix K is of the form It is a matter of algebra (see appendix E.3) to show that the matrix γ −1 β can be perturbatively calculated as where The above imply that the eigenvalues at zeroth order are As expected, they are all equal, and, thus, they do not determine the eigenvectors. At first order the matrix γ −1 β is proportional to the matrix δ i+1,j + δ i,j+1 . The determination of its eigenvectors is a simple problem. The normalized eigenvectors v j are v j i = 2 N + 1 sin ijπ N + 1 (4.25) and the eigenvalues of the matrix γ −1 β at first order equal Now we may apply degenerate perturbation theory to determine the eigenvalues at second order. They equal It is a matter of algebra to show that (4.28) The above eigenvalues imply that the entanglement entropy equals Interestingly enough, a similar cancellation between the contributions from all eigenvalues, but two, one from each subsystem, occurs in the calculation of mutual information in this case too. One can show that at this order The above formula may look quite dissimilar to the formula (4.16) that we found in the case of the non-degenerate perturbation theory. However, it is exactly the smooth limit of the latter as k i → k and l i → l, i.e. (4.31) The non-degenerate and degenerate perturbation theories resulted in different results for the entanglement entropy, but in the same result for the mutual information. This hints that the mutual information is determined by an underlying matrix object, which has the same double hierarchy as the matrix γ −1 β at zero temperature, and, thus, at this order in the l/k expansion it has only two non-vanishing elements. This is not unexpected, since the symmetry property of the mutual information enforces the latter to depend only on the entangling surface (in this case the point that separates the two subsystems) and not the subsystems. Whether the two approaches provide different results at higher orders is an issue that requires further investigation. At leading order, the difference of the two approaches, is restricted to the thermal contributions to the entanglement entropy, thus, irrelevant to our interests.
The formula (4.31) also has a high temperature expansion of the form which coincides with the l/k-expansion of the high temperature formula (3.26).

Low Temperature Expansion
In the previous section, we managed to find an l/k expansion for the mutual information in the case of a chain of oscillators. Although there is an ambiguity at the process of the perturbative calculation of the eigenvalues of the matrix γ −1 β, as long as the mutual information is considered, this ambiguity disappears, at least at this order in the perturbation theory.
We also showed that the expressions agree with the expected form for the high temperature expansion of the mutual information. However, as we will see in the next subsection with the study of two indicative example chains of oscillators, at low temperatures, the expressions we obtained with the l/k expansion fail to approximate successfully the actual mutual information. The underlying reason for this is the fact that at low temperatures, most eigenvalues tend to zero (at least at this order in the perturbation theory). As a result, the expansive formulae for the calculation of the contribution of an eigenvalue to the entanglement entropy are not correct, since they reach a singular point. The case of low temperatures should be dealt separately, making the appropriate adaptations of the relevant formulae. This is performed in the appendix F. It turns out that the low temperature expansion of the mutual information is given by n is the non-vanishing eigenvalue of the matrix γ −1 β at zero temperature, which at this order in the l/k expansion reads and k (2) i is the second order correction of the eigenvalues of the matrix K in a nondegenerate perturbation theory approach, namely The first line of the equation (4.33) is trivially twice the zero temperature entanglement entropy. The second line is the thermal correction to the mutual information at low temperatures, which clearly is exponentially suppressed.

Two Characteristic Examples
Let us now consider two characteristic example chains of oscillators. The one is a chain, whose couplings matrix is of the form In an obvious way, this is a chain, where the non-degenerate perturbation theory is appropriate for the determination of the eigenvalues of the matrix γ −1 β. We compare the l/k expansion (4.16), its high temperature expansion (4.17) and its low temperature expansion (4.33) to numerical results. The numerical calculation of entanglement entropy and the mutual information is performed with the use of Wolfram's Mathematica. The comparison of the numerical and analytic results for various values of k are shown in figure 2. In all cases l is considered equal to −1.
Furthermore, in all cases we assume N = 60 and n = 30. It is evident that the perturbative formulae approximate the numerical results successfully, especially for large values of the parameter k.
The second chain of oscillators that we consider has a couplings matrix of the form Obviously, this is the basic example where the degenerate perturbation theory applies. This case is also very interesting, as it can be obtained from the discritization of the degrees of freedom of 1 + 1 dimensional free massive scalar field theory.
In this case one can obtain another analytic formula. Whenever, the couplings matrix is of the form of a chain of oscillators, i.e. only neighbouring oscillators are coupled, the high temperature expansion formula (3.26) assumes a simple form, as the block K B contains only one non-vanishing element, which is equal to l n , namely In the case of the chain (4.37), it is possible to calculate exactly the above expression, since the eigenvectors of the block K A are known (see e.g. appendix E.3), Therefore, in this case we also have an expression for the high temperature expansion of the mutual information, which is exact in l/k. As in the previous example, the analytic formulae are compared with numerical calculations for various values of k. All examples have l = −1, N = 60 and n = 30. The perturbation theory is in good agreement with the numerical results, whenever the parameter k is large. Notice that there is an interesting change in the behaviour numerical ε expansion-high T ε expansion-low T ε expansion  of the mutual information as k gets lower. There is a critical value of k, where the dependence of the mutual information on the temperature ceases being monotonous. This is exactly the value where the coefficient of 1/T 4 in the exact high temperature expansion (4.38) vanishes. This critical k, for large values of n and N tends exponentially fast to the value k = −5/2l. As k further reduces, another more dramatic change occurs. The mutual information at infinite temperature becomes larger than that at zero temperature.

Discritizing the Degrees of Freedom in a Spherical Lattice
In this section, we extend the results of sections 3 and 4 to quantum field theory. We restrict our attention to the case of a free real scalar field in 3+1 dimensions. The analysis closely follows that of [1]. The Hamiltonian equals where Y lm are the real spherical harmonics namely, which form an orthonormal basis on the sphere S 2 . The moments ϕ lm (x) and π lm (x) obey the canonical commutation relations [ϕ lm (x) , π lm (x)] = iδ (x − x ) δ ll δ mm . The Hamiltonian expressed in terms of ϕ lm (x) and π lm (x) assumes the form (5.5) Had we descritized the radial coordinate appropriately, we would have resulted in an expression of the Hamiltonian containing countably infinite, canonically commuting variables, thus a Hamiltonian that can be dealt with the techniques of section 3.
In order to achieve that, we introduce a lattice of spherical shells with radii x i = ia with i ∈ N and 1 ≤ i ≤ N . The radial distance between consequent spherical shells sets a UV cutoff equal to 1/a to our system, while the overall size of the lattice sets an IR cutoff equal to 1/(N a). The Hamiltonian of the discretized system can be obtained from equation (5.5) substituting, x → ja, ϕ lm (ja) → ϕ lm,j , π lm (ja) → π lm,j a , which results in (5.7) Different l and m indices do not mix and moreover m does not enter directly in the Hamiltonian, thus, the problem can be split to infinite independent sectors identified by index l, each being composed by 2l + 1 identical subsectors. Thus, the overall entanglement entropy can be calculated as the series where S l (N, n) is the entanglement entropy corresponding to the Hamiltonian The latter contains a finite number of degrees of freedom, thus, S l (N, n) can be calculated using the results of section 3. For large l, the matrix describing the N oscillators becomes almost diagonal and as a result for large l the system is almost disentangled. As a consequence, it can be shown that the series (5.8) is converging [1].
It follows that the mutual information can also be calculated as the series where I l (N, n) is the mutual information corresponding to the Hamiltonian (5.9).

The Large R Expansion
We intend to study the dependence of the entanglement entropy and more importantly the mutual information, as a function of the size of the entangling sphere. For this purpose, we assume that the entangling sphere lies in the middle between the n-th and (n + 1)-th site of the spherical lattice. It follows that if we define then the radius of the entangling sphere will be R = n R a. (5.12) In the following we will study the expansion of the entanglement entropy and the mutual information for large radii R of the entangling sphere, or equivalently for large n R . The series (5.8) or (5.10) cannot be summed directly. Instead we will approximate them using the Euler-MacLaurin formula, closely following [23]. This reads where the coefficients B k are the Bernoulli numbers defined so that B 1 = 1/2. Using this formula, we may approximate the series (5.10) with the integral I ∞ 0 d (2 + 1) I (N, n, ( + 1)). (5.14) We are interested in the behaviour of this integral for large R. This behaviour cannot be isolated trivially, since n R appears in the integrand in the form of the fraction ( +1)/n 2 R and takes arbitrarily large values within the integration range. This can be bypassed performing the change of variables ( + 1)/n 2 R = y. Then the integral formula (5.14) assumes the form which can be expanded for large n R . The term that is proportional to the highest power of n R that appears in this expansion is the one which is proportional to n 2 R , i.e. the "area law" term. When the size of the entangling sphere is sufficiently large, the mutual information is dominated by the area law term, in agreement with [22]. This term receives contributions only from the integral term of the Euler-MacLaurin formula (5.13). Therefore, the large R behaviour of the mutual information is determined by equation (5.15).

The Hopping Expansion for the Area Law Term
The Hamiltonian (5.9) describes a system of coupled oscillators with couplings matrix, which can be approximated as for the purpose of the determination of the leading "area law" term in the large R expansion. Trivially, the Hamiltonian (5.9) describes a chain of oscillators and we may use the results of section 4. Substituting the hopping expansion of the mutual information for a chain of oscillators (4.16) with the couplings (5.16) to the integral formula (5.15) and expanding for large n R yields

(5.17)
This formula has the high temperature expansion which unlike the general formula for coupled oscillators contains an 1/T 2 term. This seeming contradiction is due to the fact that we have integrated contributions from arbitrary high angular momenta . The high temperature expansion (3.26) holds for temperatures higher than the eigenvalues of the matrix K. However, when one considers arbitrarily high angular momenta, these eigenvalues become arbitrarily large. This would be resolved had one introduced a physical cutoff to the angular momenta. We will return to this at the next subsection.
As we have seen in section 4, the 1/µ expansion fails at low temperatures. In the same section, we obtained the appropriate low temperature expansion for the mutual information (4.33). Substituting this low temperature expansion into the Euler MacLaurin formula yields . (5.19) The first term, I T =0 , is the zero temperature mutual information, which is simply twice the zero temperature entanglement entropy. Perturbative expressions for this term in the l/k expansion may be obtained from [23]. Unlike the general case, the integral in the above formula cannot be performed analytically. However, its behaviour is dominated by the exponential part. The exponent, i.e. the function f (y) = 2 + µ 2 a 2 + y T 1 + 3 4y (2 + µ 2 a 2 + y) (5.20) has only one minimum in (0, ∞), which lies at at this order in l/k. Therefore, a saddle point approximation may be performed. The value of the function f and its second derivative at the minimum, at this order equal It is then a matter of algebra to show that Figure 4 shows the dependence of the coefficient of the "area law" term of the mutual information on the temperature, for various values of the mass parameter. For each mass, the first order result in the l/k (5.17), as well as the high temperature (5.18) and low temperature (5.23) expansions are displayed. The analytic formulas are compared with a numerical calculation of the mutual information, which is based on the direct numerical specification of the eigenvalues of the matrix γ −1 β and is performed with the use of Wolfram's Mathematica. We used the third order result for the entanglement entropy at zero temperature, derived in [23], in order to approximate the I T =0 term in the low temperature formula (5.23). It is evident that the analytic formulae that we obtained in this section are in good agreement to the numerical results, especially for large masses.

Dependence on the Regularization
As explained in [23], the regularization scheme that we use in this section is quite peculiar. The radial and angular excitations of the field are treated differently; while numerical high T of ε expansion ε expansion at low T ε expansion there is a UV cutoff equal to 1/a for the radial ones, the angular ones are integrated up to infinite scale. One can enforce a more uniform regularization introducing a cutoff at the angular momenta of the form l max = cR/a. The appropriate selection for c in 3+1 dimensions, so that the density of the degrees of freedom at the region of the entangling surface is homogeneous, is c = 2 √ π. Then, the results of the previous subsection serve as an upper bound for the area law term. It has to be noted that had one desired to generalize these results to an arbitrary number of dimensions, they would have found that the integral without the angular momentum cutoff diverges at 4 + 1 and higher dimensions; this upper bound exists only in 2 + 1 and 3 + 1 dimensions. Obviously, the introduction of this cutoff yields the coefficient of the area law term of the mutual information finite at all dimensions. Returning to 3 + 1 dimensions, such a regularization yields This formula has the high temperature expansion (5.25) This is exactly what should be expected from the general high temperature formula (3.26). The 1/T 2 term is vanishing, whereas the 1/T 4 contains only the leading term in the 1/µ expansion (the last term of equation (3.26)), which is equal to 1/(1440a 4 T 4 ) from each angular momentum sector. As we have cutoff the angular momenta at l max = cR/a c(n R + 1/2), at leading order in n R there are c 2 n 2 R such sectors, which is consistent with our result.
The low temperature behaviour is determined by the low angular momenta. Naturally, the introduction of the angular momenta cutoff does not alter the procedure of deriving the low temperature expansion of the mutual information, as long as c > 3/2. For these values of c the formula (5.23) provides a good approximation of the mutual information at low temperature. Figure 5 shows the dependence of the coefficient of the "area law" term of the mutual information on the temperature, with the use of an angular momentum cutoff l max = 2 √ πR/a, for various values of the mass parameter. The first order expansion, as well as the low and high temperature expansions are compared to numerical calculations performed with the use of Wolfram's Mathematica. numerical high T of ε expansion ε expansion at low T ε expansion

Discussion
In a seminal paper [1], Srednicki calculated the entanglement entropy in massless scalar field theory at its ground state when the entangling surface is a sphere. It turns out that the entanglement entropy is proportional to the area of the sphere and not its volume, resembling the well-known property of the black hole entropy. This behaviour continues to hold at massive scalar field theory, where perturbative methods have been applied to calculate the entanglement entropy for a spherical entangling surface [23].
When the mass of the field is very large, the area law can be understood as a result of the locality. In such cases only correlations between nearest neighbours are important, therefore the entanglement entropy should be expected to be proportional to the number of neighbouring degrees of freedom that have been separated by the entangling surface. These are obviously proportional to the area of the entangling surface. However, the area law holds in the massless case, too. The underlying cause of this behaviour is the symmetric property of the entanglement entropy. Whenever the composite system lies in a pure state it holds that S A = S A C . Therefore, a volume term cannot appear as it should be proportional to the volume of the interior and simultaneously to the volume of the exterior of the sphere. Naturally, the entanglement entropy has to depend on the geometric characteristics of the only common feature that the interior and exterior of the sphere share, i.e. the sphere itself.
In this work we study free scalar field theory at a thermal state, generalizing the perturbative methods of [23]. It turns out that the entanglement entropy contains volume terms, which are inherited from the thermal entropy of the overall system. The presence of such terms should not be considered surprising, since the symmetry property of the entanglement entropy does not hold whenever the composite system lies in a mixed state. The entanglement entropy is not a good measure of quantum entanglement in such cases; a better measure of the correlations between a subsystem and its complement is the mutual information. This, by definition obeys the symmetry property, and, thus, it should be expected that in field theory it obeys an area law, even at finite temperature. Indeed, our perturbative calculations, as well as the numerical calculations that we performed, verify this intuitive prediction.
The coefficient of the area law term of the mutual information exposes an interesting behaviour as a function of the temperature. This coefficient reduces as the temperature increases; this is expected as the thermal effects tend to wash out the quantum correlations between the considered subsystems. However, as the temperature tends to infinity, the coefficient does not vanish, but it rather tends to a given finite value. This is a property of any harmonic oscillatory system. It turns out that the asymptotic value of the mutual information at infinite temperature is identical to the mutual information of the equivalent classical system of coupled oscillators at finite temperature.
Following the approach of [23], we found a perturbative expression for the area law coefficient, expanding in the inverse mass of the scalar field. The calculation is performed in the lowest order. It is in good agreement with the numerical calculations, especially for large values for the field mass. The calculation, although significantly more complicated than the zero temperature one, can be directly performed at higher orders, improving the accuracy of the analytic results.
Similarly to the zero temperature case, due to the particular discretization of the field degrees of freedom in radial shells, the expansion continues to work even at the massless field limit in 3 + 1 dimensions. This is due to the fact that the angular momentum effectively acts as a mass term for the corresponding moments of the field. However, it fails in 1 + 1 dimensions at the massless limit.
The original calculation of Srednicki implements a peculiar regularization. Although a lattice of spherical shells is used, introducing a UV cutoff at the radial field excitations, the angular momenta are integrated up to infinity. This scheme provides a finite result only at 2 + 1 and 3 + 1 dimensions. One may apply a more uniform scheme, introducing an angular momentum cutoff so that a similar UV cutoff applies at the angular degrees of freedom on the entangling surface. Such a regularization scheme exposes the fact that the area law term is regularization scheme dependent. Furthermore, similarly to the zero temperature case, the Srednicki regularization in 2 + 1 and 3 + 1 provides an upper bound for the coefficient of the area law term. In higher dimensions there is no such bound, however, the introduction of this more uniform regularization leads to a finite result for the area law coefficient.
Finally, another interesting property concerns the high temperature expansion of the mutual information in any harmonic oscillatory system. This expansion naturally contains even powers of 1/T . However, it turns out that the first term, namely the 1/T 2 term, always vanishes.

A The Classical Mutual Information for a Pair of Coupled Oscillators
In order to understand the nature of the remnant of the mutual information at infinite temperature, we present the classical analysis [20]. First we consider a single harmonic oscillator with eigenfrequency ω. Without loss of generality we assume that the mass of the oscillator is equal to one. In the classical limit, the probability of finding the particle at position x is inverse proportional to the magnitude of the velocity.
It follows from energy conservation, 1 2 v 2 + 1 2 ω 2 x 2 = E, that when the system has energy E, the above probability distribution assumes the form Now we turn on the temperature, introducing a canonical ensemble of harmonic oscillators. As a consequence of the fact that the period of the motion is independent of the energy, the phase space volume per energy is constant. It follows that the appropriately normalized probability distribution for the energies is This implies that the spatial probability distribution at finite temperature T is where the lower bound of the integration was taken equal to 1 2 ω 2 x 2 , since at least that much energy is required is order to reach the position x.
Let us now consider the system of two coupled oscillators of section 2, which is described by the Hamiltonian (2.10). As usual, one may introduce the canonical coordinates (2.11), which allow the re-expression of the Hamiltonian in the form (2.12), which describes two decoupled oscillators, one for each mode. Therefore, The probability distribution of the position of the first of the two coupled oscillators can be calculated integrating out the position of the second one. Simple algebra yields . We remind the reader that this is not the first time we meet this frequency. It is identical to the limiting value at infinite temperature (2.34) of the eigenfrequency of the effective single oscillator (2.25) that reproduces the reduced density matrix at the appropriate effective temperature (2.26).
It is now straightforward to find the classical version of the "entanglement" entropy, i.e. the Shannon entropy of the classical probability distribution p (x 1 ; T ), and the thermal entropy It follows that the classical mutual information is equal to This does not depend on the temperature and is equal to the asymptotic value of the quantum mutual information at infinite temperature (2.37). It follows that the quantum mutual information at infinite temperature should be attributed to classical correlations. One can trivially show that in a similar manner the classical mutual information coincides with the infinite temperature limit of the quantum mutual information in the case of an arbitrary number of coupled harmonic oscillators [20].

B Entanglement Negativity in Systems of Coupled Oscillators
In section 2, we showed that there is a finite remnant of mutual information at infinite temperature, unlike the usual behaviour in qubit systems. This remnant can be attributed to classical correlations, as we showed in appendix A. A verification check is the specification of entanglement negativity. This is defined as the opposite of the sum of the negative eigenvalues of the partially transposed density matrix, ρ T A , i.e. if λ i are the eigenvalues of ρ T A , then the negativity N will be equal to The entanglement negativity is a measure of quantum entanglement 1 . Although a non-vanishing negativity implies the presence of quantum entanglement, the opposite does not hold, when the subsystems have sufficiently high-dimensional Hilbert spaces [28]. Obviously, this is the case for harmonic oscillators, since the corresponding Hilbert spaces are infinite dimensional. Thus, finding vanishing negativity at infinite temperature is not a proof of the classical origin of the mutual information, but it is consistent with such an interpretation. In qubit systems, typically negativity vanishes at a given finite temperature and it remains vanishing at temperatures higher than that. We will show that this also holds in harmonic oscillatory systems. The techniques of section 3 can be easily generalized for the calculation of entanglement negativity.
The density matrix of a system of N oscillators in a thermal state reads (see equation (3.8)), We calculate the entanglement negativity between the first n (system A) and the last N − n (system A C ) oscillators. As usually, we decompose x as Taking the partial transpose ρ T A is equivalent to the interchange of x C and x C , which is also equivalent to the interchange of x and x . It is easy to show that after this action the density matrix assumes the form The spectrum of the density matrix is given by where the quantities ξ i are related to the eigenvalues λ i of the matrix γ −1 β as Let us consider first the case of two coupled harmonic oscillators. In this case the elements of the matrices γ and β in the expressions (B.6) are not blocks but single elements. These matrices equal The eigenvalues of the matrix γ −1 β are (B.10) Clearly, one of those, namely λ 2 , is negative at zero temperature, since whereas they are both positive at infinite temperature since Both eigenvalues are monotonous functions of the temperature, therefore there is a specific finite critical temperature T neg , defined as the single solution of the equation where λ 2 vanishes. At temperatures higher than this critical temperature, the negativity vanishes. Figure 6 shows the dependence of T neg on the ratio ω − /ω + . Appro- Figure 6 -The critical temperature T neg , as function of the ratio ω − /ω + priate expansions can be used to show that the critical temperature for large values of the ratio ω − /ω + is approximately equal to where c is the solution of the equation tanh 1 2c = 2c, namely c 0.41678. It is a matter of simple algebra to show that below the critical temperature T neg , the entanglement negativity equals Figure 7 shows the dependence of the entanglement negativity on the temperature. Figure 7 -The eigenvalues of the partially transposed density matrix (left) and the entanglement negativity (right), as functions of the temperature. For these plots it is assumed that ω − /ω + = 3/2, which implies that lim T →0 In the case of a system of N coupled oscillators, the eigenvalues λ i are determined by the equation 16) or equivalently by (B.17) The eigenvalues λ i can be re-expressed as where Λ i are the eigenvalues of the matrix Since the matrix a+b tends to the zero matrix at infinite temperature, it follows that all eigenvalues Λ i tend to infinity, or equivalently all eigenvalues λ i tend to one. This implies that the negativity vanishes at infinite temperature. Actually, since all λ i 's tend to one, it follows that they all become positive at a finite critical temperature, similarly to the two oscillators case.
One the contrary at zero temperature, the b matrix vanishes and the a matrix tends to the matrix Ω = √ K. Therefore, the eigenvalues λ i are determined by the equation 20) or equivalently by These eigenvalues come in min(n, N − n) pairs in view of Sylvester's determinant identity. There are always negative eigenvalues, therefore the system exhibits quantum entanglement. This is obviously expected, since the system lies in a pure state and has non-vanishing entanglement entropy.

C The High Temperature Expansion for Coupled Oscillators
In this appendix, we obtain the high temperature expansion for the entanglement entropy and the mutual information for systems of coupled harmonic oscillators. For this purpose, we first need to expand the matrices a, b and a + b, which are defined in equation (3.4), for high temperatures. It is simple to show that In the following, we will need the A, B and C blocks of the matrices K 2 and K 3 , in order to substitute them into formulae (C.1), (C.2) and (C.3). These are given in terms of the corresponding blocks of the matrix K by and We need to specify the high temperature expansion of the eigenvalues of the matrix γ −1 β. We recall that the matrices γ and β are defined as γ = a C − d/2 and . As a direct consequence of the equation (C.3), we have Then, definingK C ≡ K C − K T B (K A ) −1 K B and using the notation we find Adopting a similar notation for the high temperature expansions of the matrices β and γ, their definitions (3.13) and (3.14) yield and γ (1) = 1 12 The calculation of the high temperature expansion of the matrix γ −1 β, is facilitated by the use of the iterative formulae The specification of the high temperature expansion of the eigenvalues of the matrix γ −1 β is now a straightforward perturbation theory problem. The zeroth order result is obviously 1 and the eigenvectors are arbitrary. Let |v i be the eigenvectors of the matrixK C , i.e.K We expand the eigenvalues of the matrix γ −1 β as As a direct consequence of the equation (C.28), we have The specification of the next corrections to the eigenvalues is a problem identical to the usual perturbation theory in quantum mechanics. The role of the unperturbed Hamiltonian is played by −(γ −1 β) (1) and there are two perturbations, one which is of first order in the expansive parameter 1/T 2 , namely −(γ −1 β) (2) , and a second order one, namely −(γ −1 β) (3) . Therefore, D gets contributions from both perturbations Given the expansion (C.32), the corresponding quantities ξ i and the contribution of each eigenvalue to the entanglement entropy are (1) Di respectively. Notice that although odd powers of T are absent in the expansion of β Di , they appear in ξ i due to the presence of 1 − β 2 Di in the definition of ξ i . We expand the entanglement entropy as We recall the definition of the mutual information I A : The formula (2.8) implies that in the case of N coupled oscillators the thermal entropy has a high temperature expansion of the form It follows that the logarithmic terms cancel and the mutual information has a high temperature expansion of the form At zeroth order we find In an obvious manner, S Then the zeroth order contribution to the mutual information is The two last forms for I (0) , although they are expressed as determinants of matrices of different dimensions, they are equal and they are connected through the Sylvester's determinant formula. Similarly,

S
(1) A C = 1 24 TrK A and thus, Tr The two terms that are written as a sum, simplify if we write the double sum term as the symmetrized sum, (C.47) The latter implies Tr Using the definition ofK C and expressing K 2 C in terms of (K 2 ) C , using formula (C.7), yields Finally, the above equation implies that Putting everything together, the high temperature expansions of the entanglement entropy and the mutual information are given by the equations (3.25) and (3.26), respectively.

D The Low Temperature Expansion for Coupled Oscillators
At zero temperature, the matrices a and b (3.4) are not analytic functions of the temperature. Acquiring a low temperature expansion of the entanglement entropy or the mutual information is not as straightforward as the respective high temperature expansion presented in appendix C. In an obvious manner, at exactly T = 0, a = √ K and b = 0, resulting in the usual results for the ground state of the system, presented in [1]. Beyond that, we may obtain an asymptotic expansion, approximating the hyperbolic functions as a series of exponentials. More specifically, a = Ω I + 2  where the superscript in parentheses indicates the power ofΩ that appears in each term. Using the same notation for the matrices γ, β, γ −1 and γ −1 β, it is easy to show that γ −1 = γ −1 (0) γ (0) − γ (1) + γ (1) 2 − γ (2) + . . . γ −1 (0) , (D.7) thus, at leading order one recovers the zero temperature result As a result, we obtain (D.14) After some algebra, we obtain and (D. 16) which will become handy later. In order to obtain the expansions of the matrices a and b in ε, we first need to find the corresponding expansion of the powers of the matrix K. The latter equals, (1) . (E.8) Therefore, writing Since K (0) is diagonal, it is trivial to find its powers. Therefore it is a matter of simple algebra to show that at zeroth order (E.14) At first order K N (1) ij = K N 1(1) i δ i+1,j + K N 1(1) j δ i,j+1 , (E. 15) where Finally, at second order Throughout this appendix, we will use the shorthand notation .
(E. 28) In a similar manner, one can obtain the expansion for the matrix b. The formulae are identical upon the substitution of the function f 1 with the function f 2 . We proceed to calculate the matrix γ −1 β. We define two approaches for the perturbative calculation of the eigenvalues of the matrix γ −1 β and consequently the entanglement entropy and the mutual information. In this subsection, we consider the case where the off-diagonal elements of the matrix K are smaller than the differences of the diagonal ones. We refer to this approach as "non-degenerate" perturbation theory.
In this case one can consider (γ −1 β) (0) as an unperturbed, exactly solvable problem, and treat (γ −1 β) (1) and (γ −1 β) (2) as perturbative corrections. Since (γ −1 β) Now we may proceed with perturbation theory to determine the eigenvalues at second order. They equal β There are three contributions to the above formula. The first one is trivial and comes from the part of (γ −1 β) (2) that is proportional to the identity matrix. It equals The quantities ξ i are where oscillators, since we managed to perturbatively calculate these eigenvalues, we can perform this task. As we have already encountered in section D, the functions that we need to expand are not analytic at T = 0. However, they can be expanded in a series of exponentials and we expect to find that all deviations from the T = 0 result should be exponentially suppressed. First we expand the eigenvalues, which have been calculated up to second order in the l/k expansion in the previous appendix (see equations (E.39), (E.40) and (E.41)). The generic eigenvalues, i.e. all eigenvalues except β Dn and β Dn+1 , can be expanded as . . , i = n, n + 1. (F.1) The two special ones are a little different. Since they do not vanish at zero temperature they can be expanded around this value to yield