Thermal Correction to Entanglement Spectrum for Conformal Field Theories

We calculate the thermal correction to the entanglement spectrum for separating a single interval of two dimensional conformal field theories. Our derivation is a direct extension of the thermal correction to the R\'enyi entropy. Within a low-temperature expansion by including only the first excited state in the thermal density matrix, we approach analytical results of the thermal correction to the entanglement spectrum at both of the small and large interval limit. We find the temperature correction reduces the large eigenvalues in the entanglement spectrum while increases the small eigenvalues in the entanglement spectrum, leading to an overall crossover changing pattern of the entanglement spectrum. Crucially, at low-temperature limit, the thermal corrections are dominated by the first excited state and depend on its scaling dimension $\Delta$ and degeneracy $g$. This opens an avenue to extract universal information of underlying conformal data via the thermal entanglement spectrum. All of these analytical computation is supported from numerical simulations using 1+1 dimensional free fermion. Finally, we extend our calculation to resolve the thermal correction to the symmetry-resolved entanglement spectrum.


Introduction
Quantum entanglement, an inherently quantum mechanical nonlocal phenomenon, plays a central role in modern physics. Besides the general interests in quantum information theory [1,2], nowadays the concepts and anaysis of entanglement has been extended to cover diverse topics, ranging from describing collective behaviors in condensed matters [3,4] to building up spacetime of the universe [5,6,7,8,9,10]. In particular, the study of quantum entanglement brings novel insights into understanding critical phenomena in many-body systems. Especially for (1 + 1)-dimensional [(1 + 1)D] critical systems in equilibration [11,12,13,14,15,16] and out-of-equilibrium [17,18,19,20,21,22,23,24,25], the intrinsic conformal invariance leads to rigid constrains on physical proprieties, including the structure of quantum entanglement. One of the most important measures of quantum entanglement is the entanglement entropy (EE), defined as the von Neumann entropy of the reduced density matrix of the subsystem ρ A : It is analytically predicted that the scaling behavior of the EE in a (1 + 1)D critical ground state is fully determinated by the central charge c of the underlying conformal field theory (CFT) as S = c 3 ln l eff [26,11,12,13,27], where l eff is the effective length of the subsystem. This result leads to great achievement on understanding quantum criticality via entanglement, especially for minimal models that are fully characterized by a single number of the central charge c.
Besides the EE, there are also other measures of quantum entanglement that reflects interesting critical phenomena. Specifically, the entanglement spectrum (ES) has attracted considerable attentions. The ES denotes the eigenvalues of the entanglement Hamiltonian H E that is defined by ρ A = e −H E . Obviously, the EE is just a kind of statistical average of the ES, so that it is naturally to expect that the ES should contain more information [28,29,30,31,32,33]. In this context, it is interesting to consider whether the ES -as a more precise measure of quantum entanglement than the EE -can be determinated by using the conformal symmetry. In (1 + 1)D critical systems, it is found that the ES can be indeed solved through CFT techniques [24], and the conformal data can be resolved by the ES [25]. Moreover, at the continuum limit, the gap of entanglement spectrum closes and the continuous distribution of the ES is found to be dependent only on the central charge [34]. These results demonstrate the powerful constrain of conformal invariance for (1 + 1)D critical systems.
The recent progress on quantum simulations allows to realize critical systems experimentally. These advanced techniques, such as optical lattices [35,36], superconducting circuits [37], and trapped ions [38], are expected to provide a direct access to the entanglement content of many-body states. Although these systems are prepared at extremely low temperature, the thermal corrections are still necessary to be addressed for probing the emergent criticality. In particular, at finite temperature, the ideal pure state becomes a thermal density matrix. The question of whether the observation in ideal (critical) models at zero temperature are stable under temporal perturbation remains open.
In this regard, in the present work we aim to address this question by studying the ES of (1 + 1)D critical systems at finite temperatures. First of all, by using the conformal symmetry, we analytically derive the thermal correction to the ES by cutting a small region in a finite system. Second, these predictions are verified by the numerical simulations at a low temperature T = β −1 ∼ 1/L comparing to the energy gap for a finite total system. These results are not only expected to be valuable for understanding the thermal effect theoretically, but also potentially useful for detecting quantum criticality experimentally.
2 Preliminary: some known results on the entanglement in (1 + 1)D CFT Historically, the attempt to understand quantum entanglement in many-body systems started with the calculation of the EE. By definition, it requires the knowledge of the full spectrum of the reduced density matrix. However, a direct determination of the spectrum of the reduced density matrix is generally hard, which motivates an application of the replica trick to study quantum entanglement in field theory [26,11,12]. It is important to observe that the EE can be treated as a special case of the Rényi entropy S n by considering For positive integer n, the value of Trρ n A can be computed by making n copies (replica) of the theory and sewing them along the entanglement cut along the boundary of A. The analytical continuation of n is then assumed for applying the derivative on it, with taking the replica limit n → 1. Specifically, making n copies of the theory actually forms a non-trivial n-fold manifold, so called replica manifold. To calculate Trρ n A for ground state wavefunction, we further consider an Euclidean path integral representation. It gives where Z n is the partition function on the n-fold replica manifold, and Z 1 is for a single sheet. The determination of Z n is generally hard, however, for two-dimensional conformal field theories it can be calculated by using a conformal mapping of correlation functions from the usual (Euclidean) flat spacetime to the replica manifold. It gives where c is the central charge, c n is a non-universal constant that will be ignored in our further consideration, and the effective length of the subsystem l eff is determinated by the geometry of the manifold with entanglement cut. For example, when cutting a single finite region with length l from an infinite circle (with periodic boundary condition), we have l eff = l a with the lattice constant a (In the rest of this paper, we will simply choice a = 1 as the length scale). Combine with the Eq. (2), we have Meanwhile, if we take n → ∞ in Eq.
This leads to the result of S = −2 ln λ 1 that is observed in critical spin chains.
Although the original proposal of introducing the replica trick is to avoid the calculation of ES, the determination of Trρ n A is actually applicable to give the ES of 2D CFTs [34]. For accessing the spectrum information of the reduced density matrix, we investigate the distribution function P (λ) = i δ(λ − λ i ) for the eigenvalues {λ i }. The trick on evaluating the distribution P (λ) is to consider the following construction with letting z = λ−i . Here the Laurent expansion coefficient is equal to the n-th momentum of the reduced density matrix ρ A This means that knowing Trρ n A for general n can fully determinate the eigenvalues {λ i } of the reduced density matrix, which is equivalent to solving the ES: ξ i = − ln λ i . After some algebra, it gives where I is the modified Bessel function of the first kind. It is obvious that both the value of EE and the distribution of ES are only dependent by the largest eigenvalues of the reduced density matrix, which is given by the central charge of underlying CFT: The thermal correction to the Renyi entropy for a gapped theory was firstly considered by Herzog and Spillane in [39], where they conjectured the correction scales as e −βm and provided numerical evidence for a massive scalar in 1+1 dimensions. Subsequent research conformed to this point both from field theory [40,41,42,43,44,45,46,47,48] and holography [49,50,42,43,48,51,52]. In [41], Cardy and Herzog considered the following low-temperature expansion where |ψ is the first excited state that denoted by the conformal wight ∆, and the total system is chosen to have a finite size L to form a finite gap m = 2π∆/L from the ground state |0 . To the lowest order, it leads to the thermal correction of Trρ n A as where g is the degeneracy of first excited state (the ground state is considered to be unique).
Here the leading term is the zero temperature contribution of Eq. (4) for cutting a single finite region with length l from a ring with circumference L Similar to the dominate contribution, here the subleading finite temperature correction is equivalent to a correlation function on a n-fold replica manifold that can be calculated by a conformal mapping. The result is Next, we will show this thermal expansion, combining with the Eq. (7), leads to the thermal correction of the ES at low-temperature.

The entanglement spectrum in thermal states at low-temperature
In this section, we will provide an analytical derivation of the ES for thermal critical states. As discussed in the previous section, at low-temperature we have The first term in Eq. (14) is the zero-temperature result, and the latter terms come from the finite temperature β. Here, we will use a similar scheme of calculating the zero-temperature ES [34] from Trρ n A to evaluate the thermal correction to ES.
Since the zero-temperature ES is known in previous investigations, here we will only focus on the thermal correction. In the lowest-order expansion, the Eq. (14) can be rewritten to with a temperature-independent parameter This thermal correction does not change the normalization condition of the reduced density matrix Trρ A (β) = 1, since we have lim n→1 i λ n i η i = 0. Moreover, similar to calculating the largest eigenvalue λ 1 in Eq. (6), here we have The corrected largest eigenvalue is immediately obtained As a rough estimation, one can consider that the EE of thermal states is (approximately) determinated by the largest eigenvalue λ 1 (β) as the ground state It then follows At the lowest-order of a small interval limit l L, we find that this is in line with the previously known exact solution [41] However, for higher-order expansions of O l L 4 , Eq. (20) and (21) are not consistent with each other. Physically, this is because that the entanglement content of a thermal state contains more information than the pure ground state, and cannot fully determinated by a single number.
Below we try to solve the thermal correction to the entanglement spectrum for each level i. Here, we consider the distribution function Q(λ) = i δ(λ − λ i )η i , which is expected to be temperature-independent at the lowest-order of e −2π∆β/L . Analog to Eq. (7), we consider the following construction Although each term in the above expansion is known, it is not easy to obtain a compact form with the sine functions in l eff . Hence, we will focus on the small or large interval limits, where the effective length can be simplified to power functions.

Small interval limit
At the small interval limit l L, we have Use the same trick as in the zero-temperature case, this leads to where u = − ln λ 1 ln λ 1 λ . A direct comparison of the distribution function to numerical results on discrete lattice model is hard, instead, we consider the integral of It then follows to give the iterative solution of η i and λ i By using I 0 (0) = 1, I m (0) = 0(m = 1, 2, . . . ), it is easy to get which is consistent with the result in Eq. (17) at the small interval limit l L.

Large interval limit
Since for thermal state the symmetry of S(l) = S(L − l) is broken by mixing the excited states, it is also important to investigate the ES at large interval limit l → L. The expansion of Eq. (16) (for n ≥ 2) in terms of L−l L is given by In this case, the general form of η i could be approximated as where α i is a non-universal constant. The scaling of η i is then able to give the information of underlying CFT: the conformal weight ∆ of first excited state and its degeneracy g. On a lattice model, one can directly diagonalize the reduced density matrix to obtain the sequenced spectrum {λ i } for both zero and finite temperature cases. Then the change can be calculated directly by e −2π∆β/L . With increasing the value of β/L, one should observe the validity of the analytical calculation of η i at lowest-order from numerical results. We will address this in Sec. 4.

Numerical Simulation on the lattice model of free fermions
In this section, to check our analytical results of thermal ES, we present numerical simulations in a (1 + 1)D tight-binding model as a lattice realization of free fermionic CFT where c † i (c i ) is the creation (annihilation) operator on site i, and t = 1 is the hopping amplitude between nearest sites. In this model, we have the lowest conformal wight ∆ = 1/2 and the degeneracy of the first excited state g = 4. Moreover, due to the quadratic nature of the free fermions, the entanglement Hamiltonian can be written in terms of two-point correlators where C ij = c † i c j is the correlation matrix [53, 54, 55].

Probing the CFT prediction on the thermal correction parameter η i
We start with the numerical check of the thermal parameter η i . In Eq. (17), we have obtained the analytical form of the thermal parameter η 1 for the first level of ES, which is exact for all value of the subsystem size l ranging from 1 to L. The dependence of η 1 to the subsystem size is numerically calculated in Fig. 1. As expected, we find that the lattice result of  tends to the low-temperature CFT prediction with increasing the ratio of β/L. For a value of β/L = 2, the numerical results is in good agreement with the analytical form in Eq. (17).
For solving higher levels of ES, we have performed small and large interval expansions to approach the analytical form of general η i in Eq. (27) and (31). They provide different information of the underlying finite-temperature critical theory. In Fig. 2, we numerically calculate η i (i ≥ 2) with l/L ranging from 0 to 0.3. The CFT result predicts a crossover of the thermal correction parameter at a certain high level of ES. This phenomenon is well captured by our lattice simulation. Before the crossover, η i is always negative and approaches CFT prediction with a larger β. Similar trend could be found after the crossover, while η i is positive under this circumstance. Near the crossover region, although the overall tendency still satisfy the CFT prediction, a deviation with the analytical result is apparent. This discrepancy could be traced to our definition of η i . The thermal correction parameter should be, specifically, expressed by

Thermal correction to entanglement spectrum
Here we present a direct calculation of the ES with fixed l, L and β at finite-temperature in Fig. 4.
We select 2 variables to demonstrate our prediction in Figure 4, for l/L = 1/3, where {λ i } are the eigenvalues of the reduced density matrix arranged from the largest one and S(n) is the sum over the largest n eigenvalues.
Additionally, the entanglement spectrum of the ground state in (9) is fixed by previous CFT calculation in [34]. When β becomes large enough, the thermal entanglement spectrum almost conincide with the ground state entanglement spectrum, hence we choose β/L = 0.25, 0.5 and 0.75. It is evident our result could predict the thermal correction to the entanglement spectrum at low temperature limit. The CFT basically assume a continumm limit, while we could see the oscillation of λ i obtained numerically due to lattice effect.
One key observation is the crossover, above which all eigenvalues decrease and vice versa. We speculate this mode might be universal in other situation at low temperature. In addition, this tendency is consistent with the high temperature limit when all eigenstates are occupied by the same filling and the entanglement spectrum becomes uniform.

Thermal correction to the symmetry-resolved entanglement spectrum
In the past few years, the resolution of entanglement with U (1) symmetry has attracted lots of attention from both theoretical and experimental points of view [56,57,58,59,60,61,62,63,64,65,66,67,68,69,70]. Here, we show that our approach to finite-temperature ES has a straightforward extension to the symmetry-resolved entanglement spectrum (SRES).
For systems with an internal U (1) symmetry, e.g. the particle number conservation, one allows to perform a symmetry resolution of the quantum entanglement. Consider the U (1) symmetry is generated by a charge operator Q that satisfies Q = Q A ⊗ Id B + Id A ⊗ Q B for a bipartition of the total system into subsystems A and B, then the commutator [ρ, Q] = 0 implies [ρ A , Q A ] = 0. This leads to a decomposition of the reduced density matrix into a block diagonal form ρ A = ⊕ q ρ A (ρ) that labeled by the eigenvalue q of the charge Q. The symmetry-resolved Rényi entropy is then defined as where p A (q) = Trρ A (q) is the probability that the subsystem A falls into the symmetry sector q, and the symmetry-resolved entanglement entropy is The total EE can be written in terms of where the first term is called configuration entropy, and the second term is refered as fluctuation entropy that reflects the charge fluctuation between the bipartite subsystems A and B. Calculating the partition function Z n (q) from a symmetry resolution in the symmetry sector q is a hard task. It is found that the problem can be transformed to the evaluation of partition function on a replica manifold with generalized Aharonov-Bohm flux: Z n (α) = Tr ρ n A e iαQ A = q Tr ρ n A e iαq [56]. For (1 + 1)D massless bosonic field φ, the additional Aharonov-Bohm phase can be generated by inserting a vertex operator V = e i α 2π φ . This idea leads to an estimation of the charged moments where l eff = L π sin πl L for cutting a single finite region with length l from a ring with circumference L, and {∆ V , ∆ V } is the scaling dimension of the vertex operator. For fermionic theories, the operator V can be interpreted by using bosonization relation ψ ∼ e iφ . This leads to ∆ V = ∆ V = 1 2 α 2π 2 K with the Luttinger parameter K that describes the interactions between fermions, which takes K = 1 for free fermions.
Analog to the case of total ES, for solving the thermal correction to SRES, we consider the low-temperature expansion of the charged partition function as where F n (α) = which has been computed in [71].
At the small interval limit l L, we can simplify the above equation to The charged distribution function Q(λ, α) can be derived accordingly After some algebra, it gives where A(α) = 1 − 3α 2 K π 2 , B(α) = ∆ − 3α 2 K π 2 and u = 2 ln λ 1 B(α) ln λ 1 λ . Analog to the case of solving the total ES, calculating the sum of first few levels allows a self-consistent check.
Here, we have After performing a Fourier transformation from the α to the q-sector, we reach the sum of the thermal correction parameter η i (q) in each symmetry sector as As taking M = 1, this sum gives a single value of λ 1 as which is identical to the result of η 1 = lim n→∞ gF n (α). Moreover, it is easy to check that by taking α = 0, the above calculations are reduced to the case of total ES as discussed in previous section.

Discussion
In this paper, we have derived the thermal correction to the entanglement spectrum for two dimensional CFTs. With the small interval and low-temperature expansion, the thermal correction of order o(e −2π∆β/L ) has been analytically calculated for each eigenvalue. Interestingly, at the large interval limit, the correction displays some scaling behavior depending on the scaling dimension ∆ and the degeneracy g of the first excited state. This allows to extract the universal information of the underlying CFT. Moreover, in all cases, we find a clear crossover changing pattern inside the entanglement spectrum, which encodes how thermal effect reduces the underlying entanglement structure within critical systems. All of the above analytical predictions have been verified in our numerical simulation of a lattice model that realizes a free fermionic CFT.
For zero-temperature cases [34], the conformal symmetry in two dimensional theories is found to be strong enough to determine the entire entanglement spectrum by merely the central charge c of the underlying CFT. As discussed in the main text, with considering the finite-temperature correction, the information of excited states comes into the entanglement spectrum. At low-temperature limit, the correction is dominated by the first excited state and depends on its scaling dimension ∆ and degeneracy g. This leads to the possibility of extracting universal information of underlying CFT via the thermal entanglement spectrum.
We expect more practical applications for the present result on thermal correction of the entanglement spectrum. One possible direction is to investigate the thermalization process of the system from the thermal correction to the entanglement spectrum. Since the expectation value of any operatorÔ defined on the subsystem A depends linearly on ρ A as Ô = Tr(ρ AÔ ), it is then natural to consider an estimation of these expectation values from the spectral information of reduced density matrix. For example, it is straightforward to know that ( Ô (β) − Ô (T = 0) ) ∝ Trρ A (β) − Trρ A (T = 0) ∝ e −2π∆β/L . A careful study would provide more detailed information about thermalizing critical systems from the perspective of quantum entanglement.
We end up with some future directions. Within current work, we focus on two-dimensional case, while similar thermal correction might be able to get for some special cases in higher dimensional spheres S d−1 [45,44]. Secondly, we consider low temperature limit, where the ground state and the first excited states dominates the behaviour of entanglement spectrum. For a generic low-energy excited state, the n-th order Rényi entropy could be computed through a 2n-point correlation function [72]. Additionally, for more general temperature, the Rényi entropy for free fermions or free bosons at S 1 , which involves either Jacobi elliptic theta functions or Riemann-Siegel functions, was computed through the correlation function of twistor operators in high genus surface [49,40,42,73]. It is meaningful to calculate the thermal entanglement spectrum within corresponding cases based on previous work. Thirdly, recent investigation motivated more characterization of entanglement structure for critical systems, such as entanglement negativity [74,75]. It is interesting to know how underlying information emerges in thermal negativity spectrum comparing with the ground state negativity spectrum [76,77]. However, only limited results have been achieved for finite temperature Rényi negativity so far [64,78,79,80]. Finally, it was pointed out that the ground state entanglement spectrum could reproduce more information, such as operator content [81] and Affleck-Ludwig entropy [82], besides the central charge of the underlying CFT. Similar study might unveil more details in the thermal entanglement spectrum, which we leave as future work.
foundation from Westlake University. We thank Westlake University HPC Center for computation support.