Excited state entanglement in homogeneous fermionic chains

We study the Renyi entanglement entropy of an interval in a periodic fermionic chain for a general eigenstate of a free, translational invariant Hamiltonian. In order to analytically compute the entropy we use two technical tools. The first one is used to reduce logarithmically the complexity of the problem and the second one to compute the R\'enyi entropy of the chosen subsystem. We introduce new strategies to perform the computations, derive new expressions for the entropy of these general states and show the perfect agreement of the analytical computations and the numerical outcome. Finally we discuss the physical interpretation of our results and generalise them to compute the entanglement entropy for a fragment of a fermionic ladder.


I. INTRODUCTION
In the last two decades, a great effort has been done to quantify the degree of entanglement between two quantum systems. Indeed, there are different ways to measure quantum correlations [1,2]. One appropriate magnitude which can be used for this purpose is the von Neumann entropy of the reduced state. Suppose we have a bipartite system (i.e. it consists of two subsystems A and B) so the Hilbert space can be expressed as the tensor product H = H A ⊗H B of the Hilbert space of each subsystem. Let ρ be the density matrix which describes the state of the system. The von Neumann entropy of the subsystem A is defined as where ρ A is the reduced density matrix of A, i. e. ρ A = Tr B (ρ) with Tr B denoting the partial trace to the subsystem B. S 1 (A) can be interpreted as a measure of the information which A and B share. By definition S 1 (A) ≥ 0, and the zero value is attained when there are no quantum correlation between the subsystems. Besides, if the total system is in a pure state |ψ , i. e. ρ = |ψ ψ|, then S 1 (B) = S 1 (A). A generalization of the von Neumann entropy is the so-called Rényi entropy defined as In fact when we take the limit α → 1 we recover the expression (1). Both entropies share the same general properties so that the Rényi entropy of the reduced density matrix can also be employed to measure the entanglement.
The von Neumann entropy has been specially studied in extended quantum systems because it is a very suitable magnitude to analyze their universal properties in the neighbors of quantum critical points [3]. In this respect, a well-known result is the von Neumann entropy when A is a single interval of length L in a fermionic, unidimensional chain of N sites, in the ground state of a local and critical Hamiltonian, with periodic boundary conditions. In this case one has [4][5][6] S 1 (A) = c 3 log N π sin πL N + C 1 c 3 log L + C 1 , N → ∞, or for the Rényi entropy Here c is the central charge of the underlying conformal field theory and C α a non universal constant that will be computed latter. However, less attention has been paid to the entropy when the system is in an excited state which can strongly change the behavior of the entropy [7][8][9][10]. In this paper we study the Rényi entropy for certain particular states, not only the ground state, in a unidimensional fermionic chain.
In order to carry out the computation we shall make use of two technical tools. The first one is based on the work of Peschel [11] and allows to reduce the dependence of the complexity of the problem on the size of the system from 2 L to L. Instead of studying the reduced density matrix it is enough to consider the two-point correlation matrix. The second tool can be applied when the correlation matrix has the Toeplitz form, as it is in our case. Then we can find the behavior of the Rényi entropy in the thermodynamic limit using the Fisher-Hartwig conjecture (or rather theorem, as it has been proven in our case [12]).
The paper extends the works of Jin and Korepin [13] and Alba, Fagotti and Calabrese [7] which employ the previously mentioned techniques to compute the entanglement in excited states (see also ref. [9] for a different approach to the subject).
The paper is organized as follows. In the next section we formulate precisely the problem and fix the notation. In section III we discuss the connection between the reduced density matrix and the two-point correlation matrix for cases in which the Wick decomposition applies. In section IV we study the conditions under which the correlation matrix is of the Toeplitz type and the Fisher-Hartwig conjecture holds. Section V is devoted to the analytic evaluation of the Rényi entropy for our general case, while some particular examples and its comparison with the numerical results are discussed in section VI.
Finally, in section VII we present a few conclusions and comments.

II. STATEMENT OF THE PROBLEM
We consider a chain of N identical fermions with a n and a † n , n = 1, . . . , N representing respectively the annihilation and creation operator for the site n. The only non vanishing anticommutation relations are

Now given a different basis
ϕ k (n)a n , k = 1, . . . , N we consider a state where |0 represents the vacuum state in the Fock space, K ⊂ {1, . . . , N }, and is the adjoint of b k . Note that we can carry out a Gram-Schmidt like orthonormalization process and assume, without loss of generality, the canonical anticommutation relations between b k and b † k , in which case the state |Ψ is correctly normalized. In actual applications, we will consider a free, translational invariant Hamiltonian, like that of the Tight Binding Model with periodic boundary conditions where a N +1 ≡ a 1 and the state |Ψ = k∈K b † k |0 is going to be an eigenstate of the Hamiltonian. For the moment, however, we will keep the discussion completely general without making any particular assumptions on the operators b k . Now we decompose the chain into two subsets A = {1, . . . , L} and B = {L+1, . . . , N }.
Adapted to this decomposition we can factor out the Hilbert space H = H A ⊗ H B . The goal is to study the entanglement between both subsystems.
In order to do that we introduce the reduced density matrix ρ A = Tr B (|Ψ Ψ|), that in general does not correspond to a pure state, and compute its Rényi entropy. As it was discussed before the entropy of the subsystem A coincides with that of the subsystem B and provides a measurement for the entanglement between both subsystems.
Once we have obtained the reduced density matrix, we need to compute its eigenvalues in order to evaluate its Rényi entropy. Considering that the dimension of H A is 2 L the computational time grows exponentially with the size of the subsystem.
As we will see in the next section, there is an algorithm due to Peschel [11], that allows to reduce the exponential growth to a potential one.

III. WICK DECOMPOSITION
Wick theorem is a key ingredient in the perturbative expansion of quantum field theory.
It implies that the correlation function of 2J points can be decomposed into the correlation functions of the different possible pairings of the points.
We will say that a state ρ satisfies the Wick decomposition property if the correlation of an odd number of points is zero and for every J we have where d l is any linear combination of creation and annihilation operators, i.e. d l = α l a † l + β l a l , and is the set of permutations that preserve the order in every pair and |σ| is the signature of σ. It is interesting to notice that if the original density matrix of our system ρ satisfies such a property, then, the reduced density matrix ρ A also admits the Wick decomposition, as one can trivially check.
In the following we will further assume that ρ preserves the total fermionic number and therefore Tr(ρa n a m ) = Tr(ρa † n a † m ) = 0. This property is also inherited by the reduced density matrix.
In order to proceed we must solve the inverse problem. That is, given the correlation matrix C nm = Tr(ρa † n a m ) and assuming that the state satisfies the Wick decomposition property, determine the density matrix.
To achieve this goal we diagonalize the correlation matrix that has real eigenvalues in the interval [0, 1]. Denote by {φ 1 , . . . , φ N } the orthonormal basis of eigenvector with respective eigenvalues µ k , k = 1, . . . , N , and by φ k (n)a n the basis of annihilation operators for which the correlation matrix is diagonal. They satisfy, of course, the canonical anticommutation relations. We shall consider three disjoint sets in {1, . . . , N }: E 1 which contains the indices k for which µ k = 0, E 2 that contains the indices such that µ k = 1 and finally E 3 that contains the rest of indices, i.e. those with 0 < µ k < 1.
Consider first the set E 1 . It is clear that if Tr(ρc † p c p ) = 0, i.e. p ∈ E 1 and we denote by P p the orthogonal projector into the image of c p , then we can write ρ = AP p . By the same token if µ q = 1 then ρ = A (1 − P q ). Following these remarks, we introduce the ansatz for the density matrix with h = k∈E 3 k c † k c k and K the normalization constant. It is easy to see that ρ satisfies the Wick decomposition property, the eigenvalues 0 and 1 in the correlation matrix are automatically obtained and, as we will show below, with the adequate choice of h we can account also for the other eigenvalues.
If ρ describes completely the state of the system Tr(ρ) = 1. Using this fact, we can determine the normalization constant Therefore, the normalized ρ is In order to determine the coefficients k we compute the correlation Tr(ρc † k c k ) that, using the decomposition of the density matrix, it is immediately shown to be and, finally, we can relate the eigenvalues of the correlation matrix to those of h.
If, for a moment, we assume that E 1 and E 2 are empty, i. e. 0 < µ i < 1 for any i, then there is a simple relation between h and the correlation matrix in terms of the original creation and annihilation operator basis. In fact, writing M nm a † n a m and denoting by C = (C nm ) and M = (M nm ) the corresponding matrices one has Note that the magic of this procedure is that we have determined the density matrix, of dimension 2 N , through a function of the correlation matrix of dimension N . This fact allows us to go to larger values of N without exhausting the computational capabilities.
Note, however, that this dramatic simplification of the problem relies in the fact that the density matrix is supposed to satisfy the Wick decomposition theorem. A fact that must be checked in every case. One example that satisfies the Wick decomposition theorem is the density matrix associated to the Slater determinant |Ψ in (4). This case corresponds to E 3 = ∅, E 2 = K and E 1 its complementary K c . Therefore, if our state of interest ρ satisfies the Wick decomposition property, according to the previous discussion, also the reduced density matrix ρ A has such a property and hence the previous considerations will be relevant for us.
Once we have determined the density matrix we can compute its Rényi entropy and express it in terms of the eigenvalues of the correlation matrix. In fact, using (6) and (7) we have where in the last equality the product to all the eigenvalues of C (also those that are 0 or 1) can be included without any change in the result. Therefore expresses the Rényi entropy of the density matrix ρ in terms of its correlation. In matrix form we have When α → 1, that is the von Neumann entropy, the previous formulas lead to The discussion has been restricted to our case of interest in which ρ preserves the total fermionic number and, therefore, Tr(ρa † n a † m ) = 0. However a similar result is obtained in the general case, with the particularity that the transformation to the new creation and annihilation operators c † k , c k that diagonalize the correlation matrix, is of Bogoliubov type, and the previous expressions that involve the correlation matrix should be adequately modified [11].

IV. CORRELATION MATRIX OF THE TOEPLITZ TYPE
In order to be able to compute analytically the entropy of the reduced system (in the large N limit) it will be crucial that the correlation matrix is of the Toeplitz type, i.e.
Evidently, this property is a consequence of the translational invariance of the system.
We will see now under which conditions it holds.
Let us consider a state like in (4) ϕ k (n)a n , k = 1, . . . , N satisfy the canonical anticommutation relations. This fact implies that the matrix built with these coefficients ϕ k (i) is unitary (i.e. they form an orthonormal basis). Then the inverse transformation is In the previous section we have shown that the Rényi entropy of a subsystem of the chain can be computed through the correlation matrix of the system.
In particular, for the state |Ψ K , we find that the elements of the correlation matrix are But, if we demand that C nm is a Toeplitz matrix, the ϕ k (n) must satisfy That is, the product ϕ k (n)ϕ * k (m) must be invariant under translations over the lattice of the a-particles. The only functions which constitute an orthonormal basis and satisfy (10) are which are the coefficients of the discrete Fourier transformation of the a-operators. In this case the correlation matrix is of the form When the state of the chain is such that C nm is a Toeplitz matrix then the entropy of a subsystem has some interesting properties: 1. Translational invariance in the coordinate space (a-particles). From (11), it is clear that C nm only depends on the distance between sites. Otherwise we should specify not only the size of the blocks which forms A but also where they begin and end.
2. Translational invariance in the momentum space (b-particles). Let K be a particular configuration of b-excitations and K = K + ∆ which denotes that K has the same excitations which K has but with a constant displacement ∆ on the momentum lattice. The correlation matrix of K is Defining now the unitary matrix U we can rewrite the expression (12) Therefore C results from applying an unitary transformation to C , so that their eigenvalues are the same and the associated entropy is equal.
3. Entropy of the complementary state. Consider a particular configuration K so that K ∪ K c is the configuration in which all sites are excited. The correlation matrix of we can decompose the last sum in two sums of K and K c Therefore, the eigenvalues of C c nm are µ c i = 1 − µ i . Inserting them into (8) we find that the entropies of K c and K are the same.

V. ENTANGLEMENT ENTROPY FOR A SINGLE BLOCK
Suppose that the total correlation matrix of the chain is Toeplitz and a subsystem A which is a single block of L contiguous sites in the coordinate space. In this case, the correlation matrix of the subsystem, which we denote by C L , is a Toeplitz matrix as well.
Using the properties of the Toeplitz matrices, the existence of a generating function and the Fisher-Hartwig conjecture, we can find how entropy scales with the size of A in the thermodynamic limit.
If we denote by v i = 2µ i − 1 the eigenvalues of the matrix V = 2C L − I L (which is a Toeplitz matrix too), the Rényi entropy (8) can be written like It will be convenient to transform the sum in the previous expression into a complex integral on a contour that encloses the eigenvalues of V . To that end we introduce the function and applying the Cauchy's residue theorem we can rewrite the Rényi entropy (13) in terms of a complex integral: The contour C surrounds the poles of the integrand, which are points over the real line  Since we are dealing with a fermionic chain, a configuration K over the b-particles is characterized by the occupied sites. If we consider a large chain it can be approximated by the density of occupied sites, so we can describe its configuration through a periodic function g(θ) that take values in the interval [−1, 1] and is defined by In order to be able to compute the correlation function and compare with numerical data, it will be convenient to take a piecewise constant function with discontinuities at θ 1 , . . . , θ R , Note that the ground state of a free, translational invariant Hamiltonian with short range interactions, corresponds to t r taking the value +1 or −1.
Following Jin and Korepin [13], we use the Fisher-Hartwig conjecture to compute the determinant D L [14,15]. We expressg(θ) = λ − g(θ) in the form R is the number of discontinuities of g.
The Fisher-Hartwig conjecture (which has been proven for our case [12]) states that up to a factor that goes to 1 when L → ∞, the determinant of our Toeplitz matrix, D L , can be approximated by 2π 0 log(ψ(θ))dθ and with G(z) the Barnes G-function. The logarithmic derivative of D L is where the dots stand for terms that vanish in the large L limit. Inserting this expression into (14) we can compute the entropy for large L. Namely and our next task is to compute the coefficients in the expansion.
With respect to the coefficient of the linear term, A α , if we consider that ψ does not depend on θ, we have F [ψ] = ψ. Taking its logarithmic derivative, introducing it in (14) and applying the Cauchy's residue theorem we obtain As for the coefficient of the logarithmic term, B α , we have where the poles and cuts of the r-th term in the sum are represented in the figure. In order to compute the integral and remove the poles we find useful to perform an integration by parts, so that

Now a direct computation of the integral along C gives
and at last So, we write the logarithmic coefficient as an integral that, in general, should be determined numerically. In the next section we will give some examples in which B α is actually computed.
Finally, we must obtain the finite contribution C α . Proceeding like before we have We have to deal with two different integrals. The first one reads where we have performed an integration by parts. The cuts in the complex plain for the above integral are again those of fig. 2. Therefore, after integrating along the contour we get Note that β + r = β − r + 1, so applying the property of the Barnes G-function G(z + 1) = Γ(z)G(z), where Γ is the Euler gamma function, we finally obtain The computation of the second integral goes along similar lines. Define where, again, we find useful to perform the integration by parts. There are two cuts in the region inside the integral contour going from t r−1 to t r and from t s−1 to t s respectively. The integral should be decomposed into three different regions corresponding to the different cuts and computed carefully. But, considering that we only meet harmless logarithmic singularities, we can appropriately group the different terms to obtain Therefore, putting all together we have Note that if our state has g(θ) = ±1, with a finite number of discontinuities between both values, the entropy of A scales with log L because in this case the coefficient of the linear term vanishes. It has been shown in [7] that these configurations correspond to the ground state of a Hamiltonian with local interactions. On the other hand, if there are intervals with g(θ) = ±1 the entropy exhibits both linear and logarithmic contributions.
It has been proposed that these configurations correspond, in general, to the ground state of a non local Hamiltonian [7].

VI. RESULTS
In this section we present the results for the Rényi entanglement entropy of a single interval of length L, when α → 1 (von Neumann entropy) and α = 2. We perform the computation using the analytic results obtained in the previous section for the expansion in L and compare with the numerical results. As we discussed before, thanks to the Wick decomposition property [11], we can easily cover values of L = 10 3 , in the numeric computation (note that ρ A acts on Hilbert space of dimension 2 1000 which is impossible to deal with in any computer).
We consider three different pure states, each one characterized by a set K of bexcitations. In the large N limit they are described by its occupation function. With these three states we can check different aspects of the analytic expressions obtained before.
State 1: half of the b-particle sites, corresponding to the lowest absolute value of the momenta, are occupied, while the others are not. The state can be represented by (0 . . . 01 . . . 11 . . . 10 . . . 0) where 1 represents that the corresponding b-particle is present in the state and 0 that it is not. The occupation function is Note that this is precisely the ground state of the Tight Binding Model with the local periodic Hamiltonian where the constant T > 0. Since it is invariant under translations it is diagonalized by the Fourier transformation so that its eigenstates are In these configurations the correlations between the sites separated by an even distance are null. The reason for this fact is that we have separated our chain into two halves with a different, homogeneous configuration in each of them [20]. If we had divided the system into, for example, three homogeneous subsystems of different size this property would not have been fulfilled. The table I summarizes the numerical values of A α , B α and C α that we get using the expressions found in the previous section for these particular states. We compare them to those obtained by fitting the numerical entropies, computed from the correlation matrix of the subsystem.
The agreement between the outcomes is excellent. Note that when the linear term is nonzero, it provides the best fit to the numerical entropies because it dominates. On the contrary, the adjustment of the finite term is the poorest. From these results we can conclude that the large L expansion (17) works rather well from lengths of L = 10 sites. On the other hand, we must take into account the finite-size effects because the occupation function g(θ) has been defined as an approximation in the thermodynamic  identity, which can be derived analytically with the change of variables x = (t − 1)/(t + 1) that also reveals the relationship with the dilogarithm function.
Concerning the state 2, first note that its linear coefficient does not depend on α.
Actually it is easy to compute to give The ends of the ladder are joined with an inversion as it is indicated by the arrows, forming, therefore, a Moebius strip.
As it is discussed in [7], states whose entropy grows linearly with the size of the interval may correspond to the vacuum of non local Hamiltonian. This can be shown to hold in our case. In fact, consider the simple Hamiltonian that corresponds to a spin ladder as it is depicted in the figure 6 (note that one must perform an inversion when joining the two ends forming a Moebius strip rather than an ordinary ladder). If we take the coupling constants the ground state is precisely the state 2.
A similar result is obtained for state 3. If we consider the Hamiltonian H = n (J 0 a † n a n + J 1 a † n a n+1 + J N/4 a † n a n+N/4 + J N/2 a † n a n+N/2 ) + h.c., with its vacuum is the state 3.
Note that our results follow the rule of thumb according to which the entanglement entropy grows like the number of bonds that are broken when isolating the subsystem.
As it is evident from the figure 6, the number of bonds that we break when isolating a conected subchain depends linearly on the size of the interval, which explains the presence of a linear term in the expansion of the entanglement entropy. The same is true for state 3, while state 1, that is the vacuum of a local Hamiltonian, dos not have this property.

VII. CONCLUSIONS AND FINAL COMMENTS
In this paper we have studied the Rényi entropy of a subsystem of a unidimensional fermionic chain in any, translational invariant state that can be described by a Slater determinant.
Although the computational time, in principle, grows exponentially with the size of the subsystem we can circumvent this limitation and go to large sizes by employing the relationship between the reduced density matrix of the subsystem and the corresponding correlation matrix. In this way, we can compute easily the Rényi entropy through the spectrum of the correlation matrix restricted to the chosen subsystem. It is important to highlight that this method is only valid when the chain is in an state which fulfills the Wick decomposition property. Inside this set of states, there is a kind of configurations (i.e. the factorized momentum states) in which the total correlation matrix is Toeplitz. In this case, if the subsystem is a single interval of contiguous sites in the coordinate space its correlation matrix is Toeplitz as well.
Taking the thermodynamic limit and periodic boundary conditions we have found the behaviour of the Rényi entropy with the length of the interval, using the Fisher-Hartwig conjecture for Toeplitz determinants. We have checked numerically the validity of this expansion for three particular states. The first of them corresponds to the ground state of a local Hamiltonian, in which case we can compare our results with those derived from conformal field theory. We have shown the perfect agreement between both approaches.
The other two states can be viewed as the lowest energy state of non local Hamiltonians, which is illustrated with two simple examples.
A natural extension of our work is the study of the entropy for a subsystem composed of disjoint intervals. There are some recent works where this problem is addressed for the ground state [16][17][18][19]. We would like to consider a general energy eigenstate like in the present paper. In that case, the reduced density matrix of this subsystem is not Toeplitz so we can not apply the Fisher-Hartwig conjecture. Nevertheless, the relationship between this matrix and the correlation function can be exploited not only to study the entropy of a single interval but also in the case of two or more disjoint blocks in the coordinate space.
Thus it can be employed to check numerically the validity of the analytical expressions and approximations which, for the ground state, have been obtained from conformal