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 a 1/μ expansion, where μ is the mass of the scalar field. Finally, we study the high and low temperature behaviour of the area law term.

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, 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 a natural measure of entanglement is the von Neumann entropy of 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 JHEP02(2020)091 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, 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 on higher dimensional theories [10][11][12]. In more recent years, entanglement at 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], which connects the entanglement entropy in the boundary CFT to the area of minimal surfaces in the bulk. When thermal states are considered, the non-symmetry of the entanglement entropy corresponds to the existence of more than one minimal surfaces, due to the presence of the black hole, which are homologous to complementary boundary regions. This study has been extended in several works (see e.g. [15]). Most of these focus 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. It follows that 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 JHEP02(2020)091 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 formulae 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 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,

JHEP02(2020)091
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 H = 1 2 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 two thermal density matrices of the form of (2.5), one for each of the two normal modes, (2.14)

JHEP02(2020)091
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

JHEP02(2020)091
In the introduction, we argued that 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 contributions originating 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. 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 ω eff = α (2.25) 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

JHEP02(2020)091
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, but it is 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 of the considered subsystem. However, the subsystem is an open system, and, thus, a time-independent state, has to be a 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.

JHEP02(2020)091
Since, the normal modes of a harmonic system do not interact, this is an equilibrium, timeindependent 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. when k 1 → 0, as expected. Furthermore, the coefficient of the 1/T 2 term in the mutual information vanishes, which is 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 (2.37) 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 any qubit system, the related Hilbert spaces are finite dimensional. Trivially, at the infinite temperature limit, the density matrix of the composite system tends to This is a separable density matrix, implying trivially that

JHEP02(2020)091
It follows that the entanglement entropies tend to lim T →∞ whereas the thermal entropy tends to lim T →∞ The above imply that the mutual information at infinite temperature vanishes, lim T →∞ 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, a natural question concerns the origin of this mutual information remnant at infinite temperature. The mutual information I ∞ coincides with the mutual information that one calculates 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.
Discerning the classical and quantum contributions to the mutual information requires the introduction of other entanglement measures. A widely used one is the quantum discord Q [26][27][28]. In this approach, the mutual information is written as where C is the difference between the entropy of the subsystem A, S A to the conditional entropy S A|A C , maximized over all possible measurement bases of A C . This is a natural definition since C at the classical limit tends to the mutual information. The calculation of the quantum discord is a highly complicated task (it is actually an NP-complete problem), due to the problem of the specification of the basis that maximizes C. Typically, these measures are applied to qubit systems, which do not have a classical equivalent system. Unlike these systems, in our case, the classical equivalent is well-defined and the equivalent classical thermal state is also well-defined. As we commented above, the mutual information of the classical system does not depend on the temperature. Therefore, a natural definition for the classical and quantum parts of the mutual information for the coupled harmonic oscillators is (2.44) The above are directly extendable to systems of an arbitrary number of coupled harmonic oscillators and free field theory, which we are going to study in next sections.
JHEP02(2020)091 Attributing the infinite temperature remnant of the mutual information to classical correlations solely 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 [29,30]. 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 of the mutual information.
In view of the discussion above, this dependence of the mutual information on the temperature is the equivalent to the quantum "freezing" of the degrees of freedom in the context of entanglement.

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 are related to the eigenfrequencies of the system as Since the normal modes are decoupled, the density matrix of the overall system can be written as the tensor product of the thermal density matrices corresponding to each of the normal modes,

JHEP02(2020)091
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

JHEP02(2020)091
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 when K B = 0, 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. JHEP02(2020)091

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 the mutual information are 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 the 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.45) and (2.46).

JHEP02(2020)091 4 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 couplings matrix K of the form In the following, we will refer to such systems as "chains of oscillators". Apart from their own interest, this class of harmonic systems 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 off-diagonal 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. We define and then, we perform an expansion in ε (or equivalently in l/k).
In the following, we adopt a particular notation for the elements of all the involved matrices. The subscript denotes the line of the element, when it lies on top of the main diagonal, whereas it denotes its column, when it lies below the main diagonal. The superscript denotes the diagonal (i.e. the superscript 0 implies that the element lies in the main diagonal, the superscript 1 implies that it lies in the first superdiagonal, the superscript −1 implies that it lies in the first subdiagonal and so on). In other words M i,j ≡ M j−i min(i,j) . 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 term 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 JHEP02(2020)091 whereas all other matrix 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 definitions of the γ and β matrices (3.13) and (3.14). This is going to play an important role in what follows. 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 that 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 JHEP02(2020)091 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 perturbative approach is equivalent to Taylor expanding the above eigenvalues with respect to the parameter g. However, this expansion does not converge whenever 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. The formulae (3.16) and (3.17) imply that the contribution of a single eigenvalue of the matrix γ −1 β to the entanglement entropy, is equal to JHEP02(2020)091 and, thus, The first two lines of 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.17) 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 non-diagonal 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 20) JHEP02(2020)091 where (4.24) The above imply that the eigenvalues at zeroth order are Now we may apply degenerate perturbation theory to determine the eigenvalues at second order. They equal β j(2) It is a matter of algebra to show that The above eigenvalues imply that the entanglement entropy equals (4.30)

JHEP02(2020)091
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.17) 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.32) 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.32) 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 perturbative formulae for the calculation of the contribution of an eigenvalue JHEP02(2020)091 to the entanglement entropy are not correct, since they reach a singular point. Namely, the contribution to the entanglement entropy from an eigenvalue of the matrix γ −1 β of the form β Di in general is given by equation (4.15). However, as β (0) Di → 0, the quantity ξ β (0) Di also tends to zero. It follows that the series (4.15) fails being a good approximation and it has to be substituted by S Di ε 2 . Although there is no problem to the perturbative calculation of the eigenvalues of the matrix γ −1 β, this technicality enforces to deal with the case of low temperatures (or equivalently small eigenvalues) 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 where β Dn 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.34) 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.17), its high temperature expansion (4.18) and its low temperature expansion (4.34) to numerical results. The numerical calculation of entanglement entropy and the mutual information is performed via the numerical diagonalization of the matrix γ −1 β and then the substitution of its eigenvalues to the formulae (3.16) and (3.17). This task is performed using Wolfram's Mathematica. The comparison of the numerical and analytic results for various values of k is 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.38), 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.  parameter k is large. Notice that there is an interesting change in the behaviour 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 the 1/T 4 term in the exact high temperature expansion (4.39) 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

JHEP02(2020)091
We define, 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 Had we descritized the radial coordinate appropriately, we would have resulted in an expression of the Hamiltonian containing a countable number of canonically commuting variables, i.e. 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 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 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 the radius of the entangling sphere is In the following we study the expansion of the entanglement entropy and the mutual information for large radii R of the entangling sphere, i.e. 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 where we defined I (N, n, ( + 1)) = I (N, n) , (5.14) taking advantage of the fact that appears in I (N, n) only in the form of the product ( + 1). We are interested in the behaviour of this integral for large R. This behaviour JHEP02(2020)091 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.13) 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 this term, in agreement with [22]. The "area law" term receives contributions only from the integral term of the Euler-MacLaurin formula (5.12). 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  This formula has the high temperature expansion which unlike the general formula for coupled oscillators contains a 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.

JHEP02(2020)091
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.34). Substituting this low temperature expansion into the Euler MacLaurin formula yields .
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 factor of the integrand. The exponent, i.e. the function f (y) = 2 + µ 2 a 2 + y aT 1 + 3 4y (2 + µ 2 a 2 + y) (5.20) has only one minimum in (0, ∞), which lies at y min = 3/2, 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 equal f (y min ) = 2 + µ 2 a 2 /(aT ) and f (y min ) = √ 2/ aT 3 (2 + µ 2 a 2 ) , respectively. It is then a matter of algebra to show that (5.21) 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.21) expansions are displayed. The analytic formulae are compared with a numerical calculation, performed as in section 4.3. For these numerical calculations N is taken to be equal to 60, similarly to past calculations (e.g. [1]). The linear part of the curve is stable for much smaller values of N , as shown in figure 5. Further increasing the value of N does not alter the accuracy of the results significantly for the purpose of this work. The mutual information is always dominated by an area law term, as shown in [22]. The coefficient of this area law term is determined by scanning n from the value 10 to the value 50. 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.21). It is evident that the analytic formulae that we obtained in this section JHEP02(2020)091 numerical high T of ε expansion ε expansion at low T ε expansion are in good agreement to the numerical results, especially for large values of the scalar field mass.

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 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 the angular momentum 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 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.21) provides a good approximation of the mutual information at low temperatures.  Figure 6. The area law term coefficient of the mutual information as function of the temperature with an angular momentum cutoff l max = 2 √ πR/a. The dashed lines are the low and high temperature expansions of the mutual information, whereas the dotted lines are the asymptotic values for T → ∞. Figure 6 shows the dependence of the coefficient of the dominant "area law" term of the mutual information on the temperature, with 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 with the same parameters as in the previous subsection. As in the previous subsection, we used the third order result for the entanglement JHEP02(2020)091 entropy at zero temperature from [23], in order to approximate the I T =0 term in the low temperature formula (5.21). For large values of the scalar field mass, the analytic formulae that we obtained in this section are in good agreement to the numerical results.

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 entangling 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, even at finite temperature, it behaves similarly to the entanglement entropy at zero temperature. Indeed, our perturbative calculations, as well as the numerical calculations that we performed, verify this intuitive prediction; the mutual information is dominated by an "area law" term.
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 JHEP02(2020)091 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 by 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 JHEP02(2020)091 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. In a similar manner, one can trivially show that 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 consistency 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 nonvanishing negativity implies the presence of quantum entanglement, the opposite does not hold, when the subsystems have sufficiently high-dimensional Hilbert spaces [31]. 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.

JHEP02(2020)091
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 in section 3, 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 partially transposed density matrix is given by where the quantities ξ i are related to the eigenvalues λ i of the matrix γ −1 β as First, let us consider 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 Clearly, one of those, namely λ 2 , is negative at zero temperature, since lim T →0 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 7 shows the dependence of T neg on the ratio ω − /ω + . Appropriate 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 to the equation tanh 1 2c = 2c, which is approximately equal to c 0.41678.
It is a matter of simple algebra to show that below the critical temperature T neg , the entanglement negativity equals T neg 2ω + 6ω + 10ω + T neg Figure 8. 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 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 and they are continuous functions of the temperature, it follows that they all become positive at a finite critical temperature, similarly to the two oscillators case.
On 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 det

JHEP02(2020)091
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; at this limit the system lies at its ground state, which is a pure, entangled 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), at infinite temperature. 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

JHEP02(2020)091
and 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

JHEP02(2020)091
which yield 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 ma-trixK 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 D gets contributions from both perturbations

JHEP02(2020)091
Given the expansion (C.32), the corresponding quantities ξ i and the contribution of each eigenvalue to the entanglement entropy are (1) Di and 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 : 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 JHEP02(2020)091 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

JHEP02(2020)091
Using the definition ofK C and expressing K 2 C in terms of K 2 C , using formula (C.7), yields S (2) Tr K 2 C − 1 720 Tr K T B (K A ) −1 K B 2 + 1 2880 Tr K T B K B . (C.49) 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, defined in equation (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 will-known 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 As a result, we obtain At next to leading order it holds where we used the following shorthand notation (D.14) After some algebra, we obtain and (D. 16) It is straightforward to substitute the above into (D.11) and show that It is not possible to obtain analytic expressions for the eigenvalues of γ −1 β in the low temperature expansion. However, the above formula implies that the corrections to the zero temperature result are exponentially suppressed. 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 Finally, at second order K N (2) ij = K N 0(2) i δ ij + K N 2(2) i δ i+2,j + K N 2(2) j δ i,j+2 , (E. 17) where 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 Similarly to the previous steps of this calculation, we expand γ −1 β as Although γ −1 and β are symmetric, this is not the case for γ −1 β. At zeroth order we get  Now we may proceed with perturbation theory to determine the eigenvalues at second order. They equal β j(2) 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 . (E.79)