Target space entanglement in quantum mechanics of fermions at finite temperature

We consider the target space entanglement in quantum mechanics of non-interacting fermions at finite temperature. Unlike pure states investigated in arXiv:2105.13726, the (R\'enyi) entanglement entropy for thermal states does not follow a simple bound because all states in the infinite-dimensional Hilbert space are involved. We investigate a general formula of the target space R\'enyi entropy for $N$ fermions at finite temperature, and present numerical results of the entropy in a one-dimensional model. We also argue the large $N$ behaviors with a comparison to the grand canonical ensemble.


Introduction
In quantum field theories, it is interesting to consider the entanglement entropy for subregions in the base space of the QFTs. In particular in the AdS/CFT correspondence, if theories have the gravity duals, the (base space) entanglement entropy is well approximated by area of surfaces in the bulk by the Ryu-Takayanagi formula [2,3]. The formula relates entanglement in QFTs to the geometry in the gravitational theories.
However, we do not have the notion of base space entanglement for an important class of holographic theories, i.e., matrix models. For example, the BFSS model [4], which is believed to have a gravitational dual, is (1 + 0)-dimensional matrix quantum mechanics and does not have the spatial base space. We cannot use the Ryu-Takayanagi formula for theories without base space. On the other hand, it is expected that we generally have a relation between geometry and entanglement in quantum gravity [5,6,7,8,9,10,11,12,13,14]. To discuss such a relation for holographic matrix models, we need another concept of entanglement rather than base space entanglement.
One proposal is the target space entanglement [15,16,17,18,1,19,20,21]. 1 In holographic matrix models, the target spaces of matrix elements can be interpreted as bulk space, e.g. the space where D0-branes moves. Thus, it is natural to consider entanglement associated with partition of target spaces. A problem to define the target space entanglement is that the Hilbert space is generally not tensor-factorized with respect to the target space. Nevertheless, we can associate a subalgebra of operators to any subregion in the target space [15]. We can then define the reduced density matrix associated with the subalgebra, and the entanglement entropy.
Hence, it is interesting to investigate the target space entanglement of matrix models. A simple toy model is one-matrix models. It is known that a two-dimensional string theory is dual to a one-matrix model with a class of potential (see, e.g., [23,24]). The dynamics of onematrix models is reduced to that of non-interacting non-relativistic fermions by diagonalizing the matrix (see also, e.g., [23,24]). Thus, by considering the entanglement of non-relativistic fermions, we can understand the entanglement structure of a matrix model and a string theory [25,26].
In [1], one of the authors developed the target space entanglement of fermions in pure states, especially in the ground states. It is shown that, for pure states with the Slater determinant wave functions, the target space entanglement entropy is bounded as S ≤ N log 2 independently of the details of the subregion and the states, where N is the number of fermions. The derivation cannot be applied to mixed states. It is also shown in [1] that the target space Rényi entropy for the ground state of N free fermions on one-dimensional circle can be computed analytically in the large N limit. It is not sure whether we can perform such an analytical computations for other states.
Thus, it is worth extending the analysis in [1] to mixed states. One important class of mixed states is thermal states. The aim of this paper is to investigate the target space entanglement for non-interacting N fermions at finite temperature. We will see that the target space entanglement entropy for thermal states does not follow a simple formula in [1] which that for pure states does. Because of this difficulty, it is hard to analytically investigate general properties of the target space entanglement for thermal states. Hence, instead of analytical computations, we will develop formulae of the target space Rényi entropy which are useful for numerical computations. Using the formulae, we will numerically study the target space Rényi entropy for thermal states in a simple model.
The paper is organized as follows. In section 2, we summarize the basic setup for quantum mechanics of fermions and the definition of the target space Rényi entropy. In section 3, we develop a formula of the target space Rényi entropy, in particular the 2nd Rényi entropy, for non-interacting N fermions at finite temperature. In section 4, we apply the obtained formula to a concrete model which is free N fermions on a one-dimensional circle, and present the numerical results of the Rényi entropy.
2 Basic setup
We label energy eigenstates of a single particle by integers n, as |n ∈ H. Then, energy eigenstates of non-interacting N fermions are labeled by sets of N distinct integers I = {n 1 , · · · , n N }, where we suppose that n 1 < · · · < n N , as |I = 1 √ N ! σ∈S N (−) σ |n σ(1) ⊗ · · · ⊗ |n σ(N ) = √ N ! |n 1 ∧ · · · ∧ |n N , (2.4) which are normalized as I|J = δ I,J . It is also convenient to introduce the anti-symmetrized operators. Let A 1 , · · · , A N be oper-ators on H. We then define a class of operators on N H as The matrix elements in the basis (2.4) are given by We summarize useful formulae of the wedge products, e.g. the trace and the product of (2.5), in appendix A.

Fermions at finite temperature
Our aim is to develop the target space entanglement for N fermions at finite temperature 1/β. Let |n be energy eigenstates of the Hamiltonian H 1 for a single particle. Then, the states (2.4) are energy eigenstates for non-interacting N fermions with energy E I = E n 1 + · · · E n N . Using this basis, the thermal density matrix for N fermions is given by with Z N (β) = I e −βE I . The density matrix has a basis-independent expression using the wedge product as From the formula (A.17), Z N (β) can be computed from the single-particle partition function Z 1 (β) as . (2.10) The thermal n-th Rényi entropy of ρ N is given by th,(N ) := log tr (ρ n N ) 1 − n . (2.11) This is reduced to the thermal entropy − tr(ρ N log ρ N ) in the limit n → 1. By using the product formula (A.21), we have It means Accordingly, the thermal n-th Rényi entropy (2.11) is written as (2.14) In particular, taking the limit n → 1, the thermal entropy is given by which represents the standard thermal relation S = β(−F + E).

Definition of the target space Rényi entanglement entropy
Here we briefly review the definition of the target space entanglement for quantum mechanics of N fermions (see [1] for details). Let M be the entire position space where fermions live, i.e., the target space of fermions. We take a subregion A in M , and define the Rényi entanglement entropy of A for a given total density matrix ρ N . A difficulty defining the entropy is that the total Hilbert space H N (and also the single particle space H) are not tensor-factorized with respect to the position space in the first quantized picture. However, by adopting an algebraic definition of entanglement, we can define the target space Rényi entanglement entropy.
The key point is decomposing H N into the direct sum of subsectors as H N = N k=0 H k,N where k-th sector H k,N consists of states with k fermions in A and the other (N − k) fermions in the complementĀ. Let Π A , ΠĀ be projections for a single particle onto A,Ā respectively as Then, the projection onto k-th sector, Π k : H → H k,N , is given by where P − is the anti-symmetrizing projection defined in (2.2). For a given density matrix ρ N , the probability resulting in sector H k,N is which satisfies 0 ≤ p k ≤ 1 and N k=0 p k = tr(ρ N k Π k ) = tr(ρ N ) = 1. We also define the projected density matrix on H k,N as which is normalized as tr ρ k = 1. We can compute the reduced density matrix of ρ k by taking the partial trace overĀ, 20) because, in each sector H k,N , the degrees of freedom in A are factorized from those inĀ up to permutations. 2 Then, the (von Neumann) entanglement entropy for subregion A is defined as where S A k (ρ k ) is the entanglement entropy of the projected density matrix on k-th sector H k,N given by Eq. (2.21) consists of two parts. The first term in (2.21) represents a classical part, that is, the Shannon entropy for the probability distribution {p k }. The second term in (2.21), is called a quantum part which is the average of the k-th sector entanglement entropy S A k (ρ k ) for each sector with the probability distribution {p k }. The target space Rényi entropy for subregion A is also defined as (2.25) The limit lim n→1 S (n) agrees with the entanglement entropy (2.21).
In [1], it is shown that, if ρ N is a pure state whose wave function is given by the Slater determinant, the entanglement entropy (2.21) and also the Rényi entropy (2.25) follow simple formulae, and have the following upper bound independently of the Rényi parameter n and choices of the subregion. However, the derivation of the bound (2.26) cannot be applied to general mixed states ρ N . Thus, it is not sure whether the Rényi entropy for thermal states of N fermions follows the bound (2.26). We will discuss that the bound (2.26) does not hold for general thermal states.

Probability for finite temperature systems
We consider the probability p k of finding k particles in region A for thermal state (2.8). By definition (2.18), p k is given by , and we have used ρ N = P − ρ N P − because ρ N is anti-symmetric as (2.8). In the energy eigenstate basis (2.4), (2.27) is computed as It is useful to introduce the weighted overlap matrices Y m,n := e − β 2 Em m|Π A |n e − β 2 En ,Ȳ m,n := e − β 2 Em m|ΠĀ|n e − β 2 En . (2.29) The probability p k is then expressed by the wedge product as where the formula for the trace of wedge products of operators is given in appendix A.

Rényi entropy for N fermions at finite temperature
The target space (Rényi) entanglement entropy is defined by (2.21) and (2.25). Our aim is to compute it for N fermions at temperature 1/β. As we will see in the next subsection, it is easy for a single particle case (N = 1) to compute the entropies (2.21) and (2.25). For multiple fermions, it is also easy to compute the classical part (2.23) (at least numerically) because probabilities {p k } can be computed by (2.30). However, for the following reason, it is difficult to compute the quantum part (2.24) for multiple fermions. When we would like to perform numerical computations, we have to introduce a cutoff making the dimension of Hilbert spaces finite. Let the single particle Hilbert space H be cut off so that the dimension is d max < ∞. Then, the dimension of the Hilbert space H N for N fermions is dmax N which is combinatorially large. Thus, since the reduced density matrix ρ k,A in (2.22) also has combinatorially large size, it is difficult to diagonalize it and to numerically compute S(ρ k ) in the quantum part (2.24). Similarly, it is difficult to compute the Rényi entropy S (n) with general parameters n. Nevertheless, for integer n, we can show that S (n) reduces to a combination of computations of matrices with size d max . More precisely, S (n) with integers n can be written by using the weighted overlap matrices Y,Ȳ like the probability p k in (2.30). For simplicity, we will investigate the case n = 2, that is, the 2nd Rényi entropy The formula for S (2) at finite temperature is given by (3.35) which is in fact written by the weighted overlap matrices (2.29).
In addition, the situation becomes simple for large N . In that case, we may use the grand canonical ensemble, and can compute the entanglement entropy as explained in subsec. 3.5.
3.1 Entanglement entropy for a single particle at finite temperature As a first step, we consider the single particle (N = 1) case. In this case, we have the two sectors k = 0 and k = 1. The projected density matrix on each sector defined in (2.19) is Then, the reduced density matrices on A are We are interested in thermal state ρ = 1 Z 1 e −βH 1 . For the thermal state, the probabilities p k are expressed by the weighted overlap matrices Y,Ȳ in (2.29) as (2.30). For N = 1, they are given by The n-th Rényi entropy (2.25) is also computed as follows. First, note that tr A (ρ n 0,A ) = 1. Next, since ρ n 1 = 1 Thus, from (2.25) and (2.21), the n-th Rényi entropy and the entanglement entropy are given as follows In particular, the 2nd Rényi entropy for N = 1 at finite temperature takes the form

Reduced density matrix
As written in the beginning of this section, we will investigate a formula of the 2nd Rényi entropy (3.1) for multiple fermions at finite temperature. To obtain it, in this subsection, we will consider the reduced density matrix.
We would like to compute the reduced density matrix p k ρ k,A = trĀ(Π k ρ N Π k ) (see (2.19), (2.20)). For this purpose, it is convenient to introduce the projection of e −βH 1 on A,Ā as Y := Π A e −βH 1 Π A ,Ȳ := ΠĀe −βH 1 ΠĀ. (3.9) Since the projection Π k is given by (2.17), Π k ρ N Π k is written as (3.10) Because the products of wedge products of operators can be computed as (A.21), Π k ρ N Π k can be written as wedge products of operators e −βH 1 sandwiched between Π A or ΠĀ, i.e.,Ỹ ,Ȳ and also Π A e −βH 1 ΠĀ, ΠĀe −βH 1 Π A . After taking the trace overĀ of (3.10), e −βH 1 with ΠĀ are connected, and trĀ(Π k ρ N Π k ) can be written by using operators and also several traces with forms tr(Ȳ ∧ · · · ∧Ȳ ). (3.12) We will write these types of traces as Noting also that trĀ(Π k ρ N Π k ) should be an operator on k (Π A H), we can find that the reduced density matrix p k ρ k,A is expanded as where {m} represent sets of k integers as , and C {m} are numerical constants. We will determine the coefficients C {m} below, by computing the matrix elements of the both sides of (3.14) in the position basis.
First, we can compute the matrix elements of the left hand side of (3.14) as follows. By definition (see [1]), we have where y, y represent k-component position vectors in A, and z, z do (N − k)-component vectors inĀ. In the following, we will also use this notation, that is, y denotes coordinates of A, and z does coordinates ofĀ. Since ρ N takes the form (2.8), the matrix elements in the position basis are given by Thus, (3.15) becomes On the other hand, the matrix elements of the right hand side of (3.14) are given by Accordingly, (3.17) should be the same as (3.18). Thus, C {m} should satisfy (3.20) From this equation, we show below that C {m} is determined as where r j is the number of m i satisfying m i = j in {m} = (m 1 , · · · , m k ). For example, for {m} = (1, 1, 2, 3, 3), r 1 = 2, r 2 = 1, r 3 = 2, r 4 = r 5 = 0. The left hand side of (3.20) can be written using a determinant as In the cofactor expansion of the above determinant, there is a term taking the following form: Other terms in the cofactor expansion also becomes the same term after the z-integral, and the number of the terms is (N −k)! (N − i m i )! which is a number of ways to ordering i m i − k elements from z 1 , · · · , z N −k . Thus, in the left hand side of (3.20), there is a term On the other hand, in the right hand side of (3.20), we have the following terms including Y (3.26) , the same type of terms as (3.25) appear in (3.26) as By comparing (3.25) with (3.27), the coefficient C {m} is fixed as To summarize, the reduced density matrix on k-th sector is given by the formula with coefficients (3.28).

Second Rényi entropy with multi fermions
Since we have obtained the reduced density matrices (3.29), it is straightforward to compute the 2nd Rényi entropy Using the expression (3.29) and the product formula (A.21), tr A (p k ρ k,A ) 2 is computed as In this expression, we can replaceỸ ,Ȳ by the weighted overlap matrices Y,Ȳ defined in (2.29) because we have tr(Ȳ n ) = tr(Ȳ n ), tr jȲ = tr jȲ , Therefore, the 2nd Rényi entropy for N fermions at finite temperature is given by We emphasize that (3.35) is written in terms of the weighted overlap matrices Y andȲ whose sizes are independent of the particle number N , and the combinatorial complexity is reduced to the sums over {m}, {n} and σ ∈ S k .

Low temperature approximation
In this section we introduce low-temperature approximation of the 2nd Reǹyi entropy. We do not use the formula (3.35) for low temperature in our numerical computations because if we use it for large β we need to perform highly accurate computations as explained in section 4.2. We therefore use an approximation around the ground states which is valid at low temperature. The density matrix of N fermions at finite temperature is written as At low temperature (large β), the ground states are dominant in the sum. Supposing that the ground states are degenerated with degeneracy d 0 , we represent them by |1 , |2 , · · · , |d 0 . Then, the density matrix becomes in the limit β → ∞. Note that the density matrix even for the ground states is a mixed state as (3.37) when the states are degenerated. At low temperature, it is sufficient to consider the perturbation around the ground state (3.37). We now consider a perturbation up to the first excited states. We suppose that the degeneracy of the ground states with energy E 0 is d 0 and that of the first excited states with energy E 1 is d 1 . 3 Then, the reduced density matrix is approximated as where N is the normalization factor given by We will consider low temperature such as e −β(E 1 −E 0 ) 1, and keep only the linear order of e −β(E 1 −E 0 ) . Then the approximation of the 2nd Rényi entropy is given by where ρ a i ,k,A and ρ b j ,k,A are the k-th sector non-normalized reduced density matrices on A given by Therefore, if we obtain a general formula of tr(ρ a,k,A ρ b,k,A ) for N -fermion states |a , |b , we can compute (3.40). Since the N -fermion states |a , |b are labeled by N distinct one-body states, we represent them by {n 1 , n 2 , · · · , n N } for |a and {m 1 , m 2 , · · · , m N } for |b . Then the N -body wave functions are given by the Slater determinants as where χ k (x) are single-body wave functions normalized as We also introduce the overlap matrices X,X on subregions A,Ā as Matrix elements of non-normalized reduced density matrix ρ a,k,A = Π k |a a| Π k are given by We can rewrite the above expression by using minor determinants as follows. Let F N,k (n) be the set of all subsets of k ordered different elements taken from (n 1 , n 2 , · · · , n N ). For example, when N = 3, k = 2, we have F N =3,k=2 (n) = {(n 1 , n 2 ), (n 1 , n 3 ), (n 2 , n 3 )}. We use I 1 , I 2 to represent elements of F N,k (n), and alsoĪ 1 andĪ 2 for the complements of I 1 and I 2 respectively. In the above example with N = 3, k = 2, if I 1 = (n 1 , n 2 ), the complement is I 1 = (n 3 ). Similarly, we represent the set of subsets of (m 1 , m 2 , · · · , m N ) by F N,k (m). The elements are denoted by I 3 , I 4 , and the complements are byĪ 3 ,Ī 4 . We then define minor determinants det(X I i ,I j ) of k × k submatrix of X associated with sets I i , I j . For example, when I 1 = (n i 1 , · · · , n i k ) and I 4 = (m j 1 , · · · , m j k ), the minor determinant det(X I 4 ,I 1 ) is (3.48) We also define minor determinants det(XĪ i ,Ī j ) associated with setsĪ i ,Ī j . In addition we define the sign sgn(I 1 ) as the sign of permutation from (n 1 , n 2 , · · · , n N ) to I 1 ∪Ī 1 , which is a set combining I 1 andĪ 1 without changing the order. For example, when N = 3 and k = 2, if I 1 = (n 1 , n 3 ), we have I 1 ∪Ī 1 is (n 1 , n 3 , n 2 ) and sgn(I 1 ) is the sign of the permutation from (n 1 , n 3 , n 2 ) to (n 1 , n 2 , n 3 ), i.e., sgn(I 1 ) = −1. We also define the sign for I 2 , I 3 , I 4 in a similar way, and represent the product of sign as Using these minor determinants and sign, (3.47) can be rewritten as Therefore, the 2nd Rényi entropy at low temperature can be computed by (3.40) with applying the formula (3.50) for tr(ρ a i ,A,k ρ b j ,A,k ).

Large N and the grand canonical ensemble
A difficulty of the computations of entropy in the canonical ensemble is due to the condition that the particle number N is fixed. However, for large N , the canonical ensemble with fixed N can be approximated well by the grand canonical ensemble (see, e.g., [27] for the proof). Thus, we may use the grand canonical ensemble for large N , and the computations are easier than the canonical ensemble as explained below. For the grand canonical ensemble, it is convenient to use the second quantized picture (or a field theory description). In the second quantized picture, the entropy is just a base space entanglement. 4 The base space entanglement for non-interacting non-relativistic fermions at zero temperature is investigated, e.g., in [25,28,29,30,31,32,33,34].
The density matrix for the grand canonical ensemble at temperature 1/β with the chemical potential µ is given as where |I are energy eigenstates with energy E I for N -particle states. For non-interacting fermions, the grand potential Ω is simplified as where n runs over all single-body energy eigenstates. Let us work in the second quantized picture. Suppose that c n , c † n are the annihilation and creation operators of the single-body energy eigenstates |n . They satisfy {c m , c n } = 0 and {c m , c † n } = δ m,n . Then, the density matrix (3.51) is written as It is easy to compute the entanglement entropy for this density matrix (see, e.g., [35,36,32]) because it is Gaussian. We introduce the creation operators at x as Then, any (k + )-point correlation functions, G k, (x 1 , · · · , x k ; x 1 , · · · , x ) := tr(c † (x k ) · · · c † (x 1 )c(y 1 ) · · · c(y )ρ), (3.57) satisfy the following Wick's contraction rule: (3.58) For a system satisfying Wick's contraction rule, the reduced density matrix for subregion A can be computed from the restricted two-point function 36,32]. The Rényi entropy is then given by where G n A (x; x ) = A dy G n−1 A (x; y)G A (y; x ) and tr A G n A = A dy G n A (y; y).
For the density matrix (3.54), the two-point function G is given by which is the average number of particles in state |m . Then, we have for any n where Y G is a matrix with elements, Note that X m,n is the overlap matrix on A defined by (3.45). Thus, the Rényi entropy is given by In particular, the von Neumann entropy (n = 1) is given by The Rényi entropy (3.64) is a function of inverse temperature β and chemical potential µ. When we compare it with the result for fixed particle number N , we will fix µ so that the average number of particles is N . The condition is given by

Fermions on circle
We have obtained formulae of the 2nd Rényi entropy for N fermions at finite temperature. We will apply them to a concrete system, which is non-relativistic free fermions on a onedimensional circle.

Setup
We consider non-interacting N fermions on a one-dimensional circle with length L. The energy spectrum of a single particle on the circle is labeled by integers as E n = 1 2m 2πn L 2 (−∞ < n < ∞). We summarize the spectrum for the ground states and the first excited ones for N -particle case in appendix B. In our numerical computations, we cannot deal with the full Hilbert space with infinite dimensions, and have to introduce a cutoff n max so that the labels n are in the range −n max ≤ n ≤ n max . We will choose a sufficiently large n max so that the high temperature behavior of the thermal partition is reproduced in a precise order. With the cutoff, the dimension of the effective Hilbert space of N fermions is 2nmax+1 N , which is combinatorially large. Thus, it is still difficult to directly diagonal the reduced density matrix if we take a sufficiently large cutoff. Our formula (3.35) avoids the obstruction to diagonalizing such combinatorially large size matrices because the right-hand side of (3.35) can be obtained by computing several traces of (2n max + 1) × (2n max + 1) matrices Y,Ȳ .
We define a dimensionless inverse temperatureβ as Then, the Boltzmann factor is given by e −βEn = e −βn 2 . For N = 1, the thermal partition function Z 1 (β) is given by where θ 3 is the Jacobi theta function, θ 3 (q) := ∞ n=−∞ q n 2 . The partition function Z N (β) for N fermions can be computed from Z 1 (β) by using the formula (2.10).
We take a subregion A as an interval with length rL in the circle with the whole length L. The target space entanglement entropy for this subregion can be computed by using the weighted overlap matrix Y m,n ,Ȳ m,n in (2.29). Let |n be single-particle energy eigenstates. We then have Since the overlap matrices Y,Ȳ are given, we can compute the Rényi entropy from the formula (3.35). Note also that Y andȲ cannot be diagonalized simultaneously. This is the reason why it is difficult to compute general Rényi entropies.

Numerical results
We numerically evaluate the 2nd Rényi entropy for N = 1, · · · , 9. 5 We consider high temperature 6β = 0.001, middle oneβ = 0.01 and low oneβ = 1. As we mentioned above, we need to introduce the cutoff n max to perform numerical computations. We take n max so that the partition function at high temperatureβ = 0.001 is computed with good accuracy. 7 We write the partition function with cutoff as Z N (β; n max ). Then, for N = 9, we have Thus, it is reasonable to take n max = 90 for N ≤ 9. We have to take larger n max if we increase N (see subsec. 4.3). At high and middle temperature, we use the formula (3.35) to compute the 2nd Rényi entropy. However, at low temperatureβ = 1, it is not efficient to use the formula. In fact, for largeβ, we need to perform highly accurate computations like the sign problem. To illustrate this problem, let us consider the partition function Z N (β) using the formula (2.10). For large β, the dominant term of GS is the lowest energy for the single particle. Since the right-hand side of (2.10) contains a term Z 1 (β) N , it seems that the right-hand side of (2.10) have a term (e −βE (1) GS ) N . However, we know that the dominant term of GS is the lowest energy for N fermions. Thus, (e −βE (1) GS ) N is canceled 5 The results for large N using the grand canonical ensemble are given in section 4.3. 6 We regardβ as high temperature if β(E 1ES − E GS ) 1, where E GS , E 1ES are energy of the ground and first excited states (see appendix B). Since β(E 1ES −E GS ) ∼βN , we can regardβ = 0.001 as high temperature for N ∼ O(10). 7 At low temperature, we can take small n max because high energy states are suppressed.
out by other terms in (2.10). To precisely see the cancellation in the numerical computation of the right-hand side of (2.10), we need high precision computations because e −βE (N ) GS is very small compared to (e −βE (1) GS ) N for large β. For fermions on the circle, we have βE GS ∼ 8.8 × 10 −27 , which is much less than (e −βE (1) GS ) N = 1, and thus, to see the cancellation, we have to perform numerical computations with 27-digits accuracy if we use the formula (2.10). A situation is worse if we consider the 2nd Rényi thermal entropy which involves Z N (2β), and we need 54-digits accuracy. A similar problem happens, if we use (3.35) for large β. Thus, we instead use the approximation (3.40) at low temperatureβ ≥ 1.
First, we show the N -dependence of the 2nd Rényi entropy in Fig. 1 for half region r = 0.5. We also show the results with other values of r (small region r = 0.1 and large region r = 0.9) in Fig. 16 in appendix C. The plots indicate that the 2nd Rényi entropy for thermal state does not follow the upper bound (2.26) which holds for pure states. The plots also suggest that the 2nd Rényi entropy generally increases with N . However, we will see in subsec. 4.3 that this increase will not continue if we consider the grand canonical ensemble as in subsec. 3.5. Although it is difficult to compute directly the 2nd Rényi entropy for the canonical ensemble with fixed particle number N for large N 10, the equivalence of the canonical and grand canonical ensemble for large N suggests that the 2nd Rényi entropy with fixed temperature β does not grow in linear in N for large N . This result is natural for the following reason. If we increase N with fixed temperature, the system is effectively reduced to the ground states since the energy gap between the ground states and the first excited ones is proportional to N (see appendix. B). Since we know that the Rényi entropy for the ground state is proportional to log N for large N [1], we expect that the 2nd Rényi with fixed temperature is also. We next consider theβ-dependence of the 2nd Rényi entropy. Fig. 3 represents theβdependence for N = 9 and r = 0.5. We also present similar plots with other parameters (N = 8, 9 and r = 0.1, 0.5, 0.9) in Fig. 17 in appendix C. The figures imply that S (2) ∼ r × S (2) th holds qualitatively except for low temperature. It means that most of the (Rényi) entanglement entropy is given by a portion of the thermal (Rényi) entropy in the subregion. We can naively interpret that the deviation from r × S (2) th is related to quantum entanglement between the subregion A andĀ. We will also argue the deviation in subsec. 4.4.  We can also compute the difference S (2) −S (2) cl which can be interpreted as the quantum part of the 2nd Rényi entropy. Fig. 5 is the result for the half region r = 0.5. The results for other regions (r = 0.1, 0.9) are also shown in Fig. 18 in appendix C. These results indicate that the difference increases with increasing temperature and N . Thus, it means that quantum effects are important at high temperature.

Large N results
For large particle numbers N , we expect that the density matrix (2.7) with fixed N can be replaced by that for the grand canonical ensemble (3.51) with an appropriate chemical potential satisfying (3.66). As shown in subsec. 3.5, it is easy to compute the (Rényi) entanglement entropy for the grand canonical ensemble. Fig. 6 shows the plot of the ratio of the 2nd Rényi entropies for fixed N to the grand canonical withβ = 0.01, r = 0.5 (we also show similar plots with other parameters in Fig. 19 in appendix C). It indicates that the grand canonical ensemble is a good approximation of the fixed-N ensemble for sufficiently large N . Thus, we may investigate large N behaviors by using the grand canonical ensemble. In the ensemble, the Rényi entropy is computed by eq. (3.64) with matrix Y G in (3.63). Unlike the case for fixed N , the computational cost is small, and the results for N = 1, · · · , 100 are obtained relatively easily. The result for the 2nd Rényi entropy for the half region r = 0.5 is shown in Fig. 7 (see also Fig. 20 with other r in appendix C). In the computation, we need to introduce a cutoff n max as in the fixed-N case. We take n max = 150 at high temperaturē β = 0.001, 0.01 because the thermal partition function with n max seems to almost converge for n max > 100. At low temperature, smaller n max gives us a good accuracy and we take n max = 75 forβ = 1. In Fig. 7, we also show the asymptotic expression of the Rényi entropy for the ground state explained below.
In [1], the target space entanglement for the ground state of the same system (N free fermions on the circle) was investigated. It was shown that the target space Rényi entropy for a single interval at the ground state takes the following asymptotic form at large N (where N is assumed to be odd numbers): 9 S (n) where Υ n is a constant given by (4.8) In particular, for n = 1 and n = 2, we have Υ 1 ∼ 0.495 and Υ 2 ∼ 0.404. We can also compute the entanglement entropy (n = 1) for the grand canonical ensemble by (3.65).  Figs. 7, 8 indicates that the (Rényi) entanglement entropy in the large N limit with fixed β is reduced to that for the ground state. This is due to the fact that temperature is effectively small in the large N limit because contributions from excited states could be neglected as argued in the previous subsection. If we take a different large N limit such thatβN is fixed, contributions from excited states may remain, and the entropies could grow more than O(log N ). However, it is difficult to perform numerical computations with the limit because, the largerβ is, the larger cutoff n max we have to take.

Entropy inequality
In subsection 4.2, we have seen that the 2nd Rényi entropy S (2) A for interval A is qualitatively the same as a portion of the thermal entropy in the subregion as S (2) A ∼ rS (2) th . It indicates that we have S th . On the other hand, it is known that the von Neumann entropy follows the Araki-Lieb inequality [37]: (4.9) If the total system is a thermal state, we have S (1) th . Thus, the inequalities are written as (4.10) Hence, for von Neumann entropy, we have S (1) th . We would like to check whether a similar inequality holds or not for the 2nd Rényi entropy. For general systems, the Rényi entropies except for n = 1 do not follow such an inequality because they are generally not subadditive. Nevertheless, it is known that the 2nd Rényi for Gaussian states satisfies the strong subadditivity [38]. Thus, we should have S (2) A +S (2) A −S (2) th ≥ 0 for Gaussian states. Here we check whether the target space 2nd Rényi entropy satisfies the following inequalities or not: We show the N -dependence of δ + S (2) for the case where subregion A is the half region r = 0.5 at temperatureβ = 0.001, 0.01, 1 in Fig. 9, and also the case where A is a small interval r = 0.1 in Fig. 10. The figures show that the inequality δ + S (2) ≥ 0 does not hold at high temperature.
th . It shows that δ + S (2) ≥ 0 does not hold at high temperature.
A(r=0.9) |. All of the values are greater than zero, and it indicates that δ − S (2) ≥ 0 holds. β-dependence of the δ + S (2) is shown in Fig. 12 for N = 9, r = 0.5 (and also in Fig. 21 in appendix C for N = 8, r = 0.5). We also draw a continuous curve using the approximation (3.40) although the approximation is not valid for smallβ. They also indicate that δ + S (2) > 0 holds at low temperature while it is violated at high temperature.  Figure 12:β-dependence of δ + S (2) with r = 0.5, that is δ + S (2) = 2S (2) th for N = 9. The purple points are computed by (3.35). The blue curve is the plot of (3.40) although the approximation is not valid for smallβ.
As in the previous subsection, we may use the grand canonical ensemble for large N . We can check the inequalities (4.11) for the grand canonical ensemble. Figs. 13, 14, are the results of δ + S (2) for the subregion r = 0.5, 0.1. Thus, for the grand canonical ensemble, δ + S (2) ≥ 0 holds in contrast to the above Figs. 9, 10 for the fixed-N ensemble. This is consistent with [38] because the grand canonical ensemble is a Gaussian state as the Wick's rule (3.58) holds. Therefore, for small N , the fixed-N ensemble is quite different from the grand canonical ensemble because the behavior of δ + S (2) is different.  Figure 13: N -dependence of δ + S (2) in the grand canonical ensemble with r = 0.5 for N = 1, · · · , 100. All of the points satisfies δ + S (2) ≥ 0 in contrast to Fig. 9. At low temperaturē β = 1 (blue points), the behaviors for odd N are quite different from those for even N . This represents the difference of thermal entropy at low temperature arising from that of the degeneracy of the ground states.

Summary and discussions
We have investigated the target space Rényi entropy for non-interacting N fermions at finite temperature. This is an extension of the work [1] for pure states to mixed states. The 2nd Rényi entropy for thermal states can be computed via the formula (3.35) although it is hard to apply to a large N system because of the combinatorial complexity of the formula. This difficulty is related to the fact that the fixed-N ensemble is not a Gaussian state. This is quite different from the grand canonical ensemble. The grand canonical ensemble for non-interacting fermions is a Gaussian state, and it is relatively easier to compute the entanglement entropy. As shown in [15,16,1], the target space entanglement in the first quantized picture agrees with the conventional base space entanglement in the second quantized picture. However, it may be not convenient to use the second quantized picture for the fixed-N ensemble since the ensemble is not Gaussian because of the fixed-N condition. Our formula (3.35) is a rare example where an explicit formula of the Rényi entropy is available for non-Gaussian states. We applied the formula (3.35) to a concrete system, N free fermions on one-dimensional circle in section 4. We numerically computed the 2nd Rényi entropy for a single interval subregion. The results show that the 2nd Rényi entropy for a single interval is qualitatively the same as a portion of the thermal 2nd Rényi entropy (2.14) in the subregion. We also used the grand canonical ensemble and compared the result with that for the fixed-N ensemble. It seems that they are equivalent for large N . However, although we expect that the fixed-N ensemble is equivalent to the grand canonical one for large N with fixed temperature, it is not sure if the equivalence holds in a different large N limit with N β fixed. In fact, Fig. 19 implies that the agreement of the fixed-N and the grand canonical ensemble is worse at higher temperature. As we discussed, the large N limit with β fixed is somewhat trivial in a sense that the state is reduced to the ground state if the first energy gap grows with N as the model in section 4. An interesting large N limit in this model is N → ∞ with N β fixed. In this limit, it is not sure if the fixed-N ensemble can be replaced by the grand canonical one. For the ground state, the target space Rényi entropy in the large N limit is the same as the base space one for a conformal field theory, i.e., the two-dimensional free compact boson at self dual radius (see [1]). It is interesting to see if thermal states also have an agreement with a conformal field theory.
Properties of the target space (Rényi) entanglement entropy of N fermions are drastically different between thermal states and a class of pure states considered in [1]. In that paper, it is shown that the (Rényi) entanglement entropy for any pure states whose wave functions are the Slater determinants is bounded by N log 2 independently of the states. This represents that only the N states contained in the Slater determinants are involved in the computation of the entropy, and we can effectively regard the dimension of the Hilbert space as 2 N . This fact also makes it easy to numerically compute the entropy. In contrast, for thermal states, all of states in the infinite-dimensional Hilbert space are involved. Therefore, the (Rényi) entanglement entropy for thermal states is not bounded from above. In particular, it diverges in the high temperature limit as the thermal entropy does.
It will be more interesting to investigate the target space entanglement for a one-matrix model dual to a two-dimensional string theory. It is also interesting to investigate a timedependence of the target space entanglement by considering time-dependent states or quantum quench (see [39]). In addition, since the (Rényi) entanglement entropy is not a good measure of entanglement for thermal states (more generally for mixed states), it will be important to consider better measures of entanglement such as the target space negativity or relative entropy.

Acknowledgement
SS thanks support from JSPS KAKENHI Grant Number 21K13927.

A Formulae for wedge products of operators
In section 2.1, we have introduced the wedge product of operators A 1 , · · · , A N as whose matrix elements in the basis (2.4) are In this appendix, we summarize some formulae for the wedge product of operators.
We summarize the trace formula of the wedge product (A.1). The trace is given by We can replace the sum over index I = (n 1 , · · · , n N ) by 1 N ! n 1 ,··· ,n N because the summand is symmetric under the permutation of n 1 , · · · , n N . We thus have This expression can be simplified furthermore as follows. Let us focus on the following part σ ∈S N n 1 ,...,n N (−) σ (A σ(1) ) n 1 ,n σ (1) . . . (A σ(N ) ) n N ,n σ (N ) (A. 6) in (A.5). In the sum over σ ∈ S N , two permutations have the same contributions if they have the same cycle type. Since the conjugacy classes of permutation group S N are classified by the cycle types, the sum over σ can be written as a sum over the conjugacy classes of S N . In addition, the conjugacy classes are in a one-to-one correspondence with the partitions of integer N . We use label λ for the partitions of integer N . The number of elements in a conjugacy class corresponding to a partition λ is given by where r k are defined for partition λ as N = 1 + 1 + · · · + 1 r 1 + 2 + 2 + · · · + 2 r 2 For instance, for N = 3, there are the following three partitions, and r k for each partition is given as follows; 3 = 1 + 1 + 1 → r 1 = 3, r 2 = 0, r 3 = 0, (A.9) = 1 + 2 → r 1 = 1, r 2 = 1, r 3 = 0, (A.10) = 3 → r 1 = 0, r 2 = 0, r 3 = 1. (A.11) These r k denote the numbers of cycles with length k in permutation σ in the conjugacy class corresponding to λ. Thus, the signature (−) σ in the class λ is determined by r k as Furthermore, for each σ in (A.6), if σ is in the class λ, the indices contract as follows where the sum of λ runs over all the partitions of integer N , and the coefficient D λ is given by , (A.15) whose absolute value |D λ | is (A.7) and sign is (A.12). Therefore, we obtain the trace formula of the wedge product (A.1) as (A. 16) If A 1 = · · · = A N (= A), the trace formula (A.16) can be written by using determinant as In a similar way, we can show that the product formula for general N is given by ((A 1 ∧ · · · ∧ A N )(B 1 ∧ · · · ∧ B N )) I,J = 1 N ! σ∈S N (A 1 B σ(1) ) ∧ · · · ∧ (A N B σ(N ) ) I,J . (A.21)

C Other plots
In this appendix, we show some plots similar to those in section 4 with different parameters.  Figure 18: N -dependence of the classical term of the 2nd Rényi entropy and of the entanglement entropy (left), and the difference between them for the 2nd Rényi entropy (right) as similar to Fig.4 and Fig. 5 with r = 0.5. Here we take r = 0.1 (upper) and r = 0.9 (lower). The classical entropies of r = 0.1 and r = 0.9 are the same as seen in the left figures because the probability p k for a region r = 0.1 is the same as p N −k for the complement region r = 0.9.  Figure 21:β-dependence of δ + S (2) for even N . It is the plot of δ + S (2) = 2S (2) A(r=0.5) − S (2) th for N = 8. It is similar to odd N case (Fig. 12).