Analytical results for the entanglement Hamiltonian of a free-fermion chain

TL;DR

This paper derives explicit analytical expressions for the entanglement Hamiltonian of a free-fermion chain, revealing persistent deviations from conformal predictions even at large system sizes.

Contribution

It provides the first explicit large-$N$ expressions for the entanglement Hamiltonian in free-fermion chains at arbitrary filling.

Findings

Analytical expressions match numerical results

Deviations from conformal field theory persist at large $N$

Explicit formulas relate to a commuting operator

Abstract

We study the ground-state entanglement Hamiltonian for an interval of $N$ sites in a free-fermion chain with arbitrary filling. By relating it to a commuting operator, we find explicit expressions for its matrix elements in the large-$N$ limit. The results agree with numerical calculations and show that deviations from the conformal prediction persist even for large systems.

Equations

$$\mathcal{H}=2\pi L{\int }_0^{L}dx\frac{x}{L}\left(1-\frac{x}{L}\right)T_{00}$$

$$\hat{\mathcal{H}}=-\frac{1}{2}\sum _{n}t(c_n^{†}{c}_{n+1}+c_{n+1}^{†}{c}_n)+\sum _{n}d{c}_n^{†}{c}_n$$

$$C_{m,n}={\int }_{-q_F}^{q_F}\frac{dq}{2\pi }{e}^{-iq(m-n)}=\frac{\sin q_F(m-n)}{\pi (m-n)}.$$

$$\mathcal{H}=\sum _{i,j=1}^N{H}_{i,j}{c}_i^{†}{c}_j$$

$$H_{i,j}=\sum_{k=1}^{N}\,\phi_{k}(i)\;\varepsilon_{k}\;\phi_{k}(j)\$$

$${\epsilon }_k=\ln \frac{1-{\zeta }_k}{{\zeta }_k}.$$

$$T=\left(\begin{array}[]{ccccc}d_{1}&t_{1}&&&\\ t_{1}&d_{2}&t_{2}&&\\ &t_{2}&d_{3}&t_{3}&\\ &&\ddots&\ddots&\\ &&&t_{N-1}&d_{N}\end{array}\right)$$

$$t_i=\frac{i}{N}\left(1-\frac{i}{N}\right),d_i=-2\cos q_F\frac{2i-1}{2N}\left(1-\frac{2i-1}{2N}\right).$$

$${\epsilon }_k(q_F)=-{\epsilon }_{N+1-k}(\pi -q_F),{\lambda }_k(q_F)=-{\lambda }_{N+1-k}(\pi -q_F)$$

$${\varphi }_k(q_F,i)=(-1{)}^i{\varphi }_{N+1-k}(\pi -q_F,i).$$

$$\mathcal{T}=\sum _{i,j=1}^N{T}_{i,j}{c}_i^{†}{c}_j$$

$$H_{i,j}=(1-(-1{)}^{i+j})\sum _{k=1}^{N/2}{\epsilon }_k{\varphi }_k(i){\varphi }_k(j)$$

$$h_{i,j}=-\underset{N\to \infty }{\lim }\frac{H_{i,j}}{N}.$$

$$-\underset{N\to \infty }{\lim }\frac{{\epsilon }_k}{N}=\eta (\kappa ),\kappa =\frac{k}{N}\mathrm{fixed}.$$

$$\eta (\kappa )={\int }_0^B\sqrt{\frac{B-\xi }{\xi (1-{\xi }^2)}}\mathrm{d}\xi$$

$$\kappa (B)=\frac{1}{\pi }{\int }_B^1\sqrt{\frac{\xi -B}{\xi (1-{\xi }^2)}}\mathrm{d}\xi$$

$$\underset{N\to \infty }{\lim }{\lambda }_k=\lambda (\kappa )=\frac{B}{2},\kappa =\frac{k}{N}\mathrm{fixed}.$$

$$\frac{1}{\sqrt{1-{\xi }^2}}=\sum _{m=0}^{\infty }{\alpha }_m{\xi }^{2m},{\alpha }_m=\frac{1}{\sqrt{\pi }}\frac{\Gamma (m+1/2)}{\Gamma (m+1)}$$

$$\eta =\sum _{m=0}^{\infty }{\alpha }_m{I}_m.$$

$$I_m={\int }_0^B\sqrt{B-\xi }{\xi }^{2m-1/2}\mathrm{d}\xi ={\beta }_m{\left(\frac{B}{2}\right)}^{2m+1},{\beta }_m=\sqrt{\pi }{2}^{2m}\frac{\Gamma (2m+1/2)}{\Gamma (2m+2)}.$$

$$\eta =\sum _{m=0}^{\infty }{\alpha }_m{\beta }_m{\lambda }^{2m+1}.$$

$$h=\sum _{m=0}^{\infty }{\alpha }_m{\beta }_m{T}^{2m+1}.$$

$$(T^{2m+1}{)}_{i,i+r}=\sum _{l_1,\dots ,{\text{l}}_{2m}}{T}_{i,l_1}{T}_{l_1,l_2}\dots T_{l_{2m},i+r}.$$

$$|l_j-l_{j-1}|=1,j=1,\dots ,2m+1,l_0=i,l_{2m+1}=i+r$$

$$(T^{2m+1}{)}_{i,i+2p+1}=(\genfrac{}{}{0ex}{}{2m+1}{m-p})z^{2m+1},z=t_{i+p}=\frac{i+p}{N}\left(1-\frac{i+p}{N}\right)$$

$$h_{i,i+2p+1}=\pi \frac{(2p-1)!!(4p-1)!!}{2^{p}p!(2p+1)!}{z}^{2p+1}{}_3{F}_2\left(p+\frac{1}{4},p+\frac{1}{2},p+\frac{3}{4};p+1,2p+2;(4z{)}^2\right)\text{vspace}0.2cm$$

$$h_{i,i+1}=\pi x(1-x){}_3{F}_2\left(\frac{1}{4},\frac{1}{2},\frac{3}{4};1,2;{\left[4x(1-x)\right]}^2\right)$$

$$H_{i,j}=\sum _{k=1}^{\bar{N}}{\epsilon }_k{\varphi }_k(i){\varphi }_k(j)-(-1{)}^{i+j}\sum _{k=1}^{N-\bar{N}}{\hat{\epsilon }}_k{\hat{\varphi }}_k(i){\hat{\varphi }}_k(j)$$

$$h=h_{-}-J{\hat{h}}_{-}J,T=T_{+}-J{\hat{T}}_{+}J.$$

$$\eta (\kappa )={\int }_A^B\sqrt{\frac{B-\xi }{(\xi -A)(1-{\xi }^2)}}\mathrm{d}\xi ,\lambda (\kappa )=\frac{B-A}{2},A=\cos q_F.$$

Full text

Dedicated to John Cardy on the occasion of his 70th birthday.

1 Introduction

In entanglement studies [1, 2], one divides a system into two parts and asks how these are coupled in the quantum state under consideration. This information is contained in the reduced density matrix for one of the subsystems which can be written in the form $\rho=\exp(-\mathcal{H})/Z$ with an operator $\mathcal{H}$ now commonly called the entanglement Hamiltonian [3]. Its eigenvalue spectrum determines the degree of the entanglement and thus is of primary importance. However, also its structure is interesting since it shows universal features which depend on the way in which the system is divided.

For non-critical integrable chains divided into two halves, the relation of $\rho$ to corner transfer matrices leads to an operator $\mathcal{H}$ with terms which increase linearly away from the boundary between the subsystems, see [4]. This holds for spin chains [5] as well as for coupled oscillators [6] and is also found in the continuum [7]. For critical quantum chains, conformal field theory gives results also for other cases, as discussed most recently by Cardy and Tonni [8]. Thus, for a subsystem of length $L$ in an infinite chain, one finds

[TABLE]

where $T_{00}$ is the energy density in the physical Hamiltonian [9, 10, 11, 8]. The parabolic weight factor becomes linear near the ends and can be viewed as the combined effect of the two boundaries.

On the lattice, such a form was found in numerical studies of a fermionic hopping model equivalent to an XX chain [12] and in an XXZ chain for the value $\Delta=1/2$, when the ground state is particularly simple [13]. However, the nearest-neighbour terms in $\mathcal{H}$ were not exactly parabolic and there were couplings to more distant neighbours, albeit of smaller magnitude. Due to the small subsystems, no clear answer could be given for the asymptotic behaviour.

In this paper, we look once more at the free-fermion model and settle the problem for this system. Starting from the insight obtained by numerical calculations, we develop an analytical treatment in the large-$N$ limit. It is based on a paper by Slepian who studied the matrix problem underlying the determination of $\mathcal{H}$ already in the seventies in the context of time- and band-limited signals [14, 15]. The strategy is to express $\mathcal{H}$ in terms of a closely related tridiagonal matrix $T$ and powers of it. For a half-filled system, this leads to explicit formulae for all the hopping matrix elements in terms of generalized hypergeometric functions. For other fillings, the coefficients in the resulting series are more complicated and the summations have to be done numerically. In this way we confirm that the deviations mentioned above do persist in the large-$N$ limit. In particular, the hopping to more distant neighbours for half filling shows a parabolic law raised to the power $r$ near the boundaries, where $r$ is the range. In the bulk, there is an additional enhancement and the maximal amplitude decreases with $r$ as $r^{-3}$. The analytical expressions are checked against numerical results obtained directly for large subsystem sizes. This requires extreme accuracy in the diagonalization of the correlation matrix from which $\mathcal{H}$ follows (roughly $N$ digits for $N$ sites) and such calculations have been done so far only in one very recent study where the continuum version of the present problem was addressed [16].

In the following Section 2 we formulate the problem and give the basic expressions. Then, in Section 3, we present numerical results for the single-particle spectrum and the matrix elements of $\mathcal{H}$ for half filling. In Section 4 we recall some asymptotic results of Slepian and present the application to our problem with the resulting formulae for a half-filled system. In Section 5, the more complicated case of arbitrary filling is discussed and in Section 6 we conclude with a discussion. Results for the extremal eigenvalues and some mathematical details are collected in two appendices.

2 Setting and basic formulae

We consider the ground state of an infinite chain with the Hamiltonian

[TABLE]

where we take the hopping $t=1$ and use the site energy $d$ to adjust the ground-state filling. With the single-particle energies $\omega_{q}=-\cos q+d$, the Fermi momentum is given by $\cos q_{F}=d$. All entanglement properties are determined by the correlation matrix $C_{m,n}=\langle c^{{\dagger}}_{m}c_{n}\rangle$ in the ground state

[TABLE]

In the following, we always consider a subsystem of $N$ consecutive sites denoted by $i$ or $j$. Since the ground state is a Slater determinant, the entanglement Hamiltonian $\mathcal{H}$ has again free-fermion form [17]

[TABLE]

and the matrix $H$ is given in terms of the eigenfunctions $\phi_{k}(i)$ and eigenvalues $\zeta_{k}$ of the correlation matrix $C_{i,j}$ restricted to the subsystem

[TABLE]

where

[TABLE]

An important (and intriguing) feature is that the restricted correlation matrix $C$ and thus the matrix $H$ commute with a tridiagonal matrix $T$ of the form

[TABLE]

with elements

[TABLE]

This matrix can already be found in Slepian’s paper [14] (with diagonal terms differing from (8) by an additive constant) and was later rediscovered in [18]. Its eigenvalues will be denoted by $\lambda_{k}$ and they will play a central role in the analytical treatment of the problem.

The matrices $C$, $H$ and $T$ have common eigenvectors and a common particle-hole symmetry. If one changes $q_{F}$ into $\pi-q_{F}$, they transform into $J(1-C)J,-JHJ$ and $-JTJ$, respectively, where $J$ is the diagonal matrix $J_{i,j}=(-1)^{i}\delta_{i,j}$. For the eigenvalues, when ordered according to magnitude, this means

[TABLE]

and the corresponding eigenvectors satisfy

[TABLE]

These relations will be used later in the construction of $H$. Pictorially speaking, the $\varepsilon_{k}$ are mostly positive for small filling (since then most of the $\zeta_{k}$ which are occupation numbers lie near zero) and mostly negative for large filling. At half filling, the spectrum is symmetric.

Finally, if one writes in analogy to (4)

[TABLE]

the operator $\mathcal{T}$ has, up to a factor of $-1/2$, exactly the same form as the physical Hamiltonian (2) but with nearest-neighbour hopping and site energies which are both multiplied by essentially the same parabolic factor increasing from the boundaries of the subsystem towards the middle. Thus it has the structure of the entanglement Hamiltonian (1) as found in conformal field theory. The actual entanglement Hamiltonian $\mathcal{H}$ will turn out to be similar but not identical.

3 Half filling: Numerical results

We present here the basic properties of $H$ and $T$ as obtained from numerical diagonalization for half filling, $q_{F}=\pi/2$. It was noted already in [18] that in this case the matrices $H$ and $\pi NT$ have the same low-lying eigenvalues. This is shown in Fig. 1 on the left where $\varepsilon_{k}$ and $-\pi N\lambda_{k}$ are plotted for a subsystem of $N=40$ sites. (The minus sign was left out in [18].) One sees that the small eigenvalues are practically indistinguishable on the chosen scale. Only for the large ones there is a small difference. Note that the largest $\varepsilon_{k}$ are about $\pm 66$ which means that the corresponding $\zeta_{k}$ differ from 0 or 1 by only $e^{-66}\simeq 10^{-29}$. Thus one needs at least 30 digits in the diagonalization of $C$ in order to find these $\varepsilon_{k}$ reliably. The actual calculation with Mathematica used 40 digits.

The dependence of the maximal values $\varepsilon_{max}$ and $\pi N\lambda_{max}$ on the subsystem size $N$ is shown on the right of Fig. 1. Basically, both quantities increase linearly. However, the slope is not the same and therefore their difference also increases. This scaling with $N$ is also found analytically and will be discussed further in Section 4. The figure shows the excellent agreement between the numerics and the asymptotic results already for small $N$. The scaling is not restricted to the maximal eigenvalues and suggests to work with the intensive quantities $H/N$ and $T$.

In Fig. 2 the leading matrix elements in $-H/N$, calculated from (5), are displayed. Shown are the amplitudes for hopping to first, third and fifth neighbours as functions of the position, again for $N=40$ sites. Hopping to second, fourth, etc. neighbours and diagonal terms are strictly zero due to (9),(10). The figure resembles the one in [12] for $N=16$, but the scale factors to make the longer-range hopping visible are smaller here, i.e. these terms are relatively larger. Shown for comparison is also the nearest-neighbour hopping in $\pi T$ with its strictly parabolic form. Since the two matrices have the same eigenfunctions, all differences in their matrix elements are tied directly to the differences in the spectra. Although these were seen to be relatively small, they occur for those eigenvalues which enter with the largest weight into the spectral representation (5) and thus produce significant effects. In the next section, we will derive analytical formulae for the shown quantities.

4 Half filling: Analytical treatment

The half-filled case is particularly simple due to the invariance of $H$ under a particle-hole transformation. Indeed, using the symmetry properties (9) and (10) of the spectrum and the eigenvectors, one can rewrite the entanglement Hamiltonian as

[TABLE]

where the sum involves only the negative branch of the spectrum. The matrix elements $T_{i,j}$ have an analogous expression, with $\varepsilon_{k}$ replaced by the positive $\lambda_{k}$. As remarked in the previous section, the operator $H$ has extensive character. It is thus useful to define its density which, in the thermodynamic limit, is given by

[TABLE]

The main idea of the following calculation is to relate the spectra $-\varepsilon_{k}/N$ and $\lambda_{k}$ for $N\to\infty$ in terms of a series expansion. From Fig. 1 it is obvious, that to leading order the relation is linear with a coefficient $\pi$. However, the behaviour for larger eigenvalues hints at the presence of higher order terms. To obtain the series analytically, we shall use the results of [14] where the following uniform asymptotic limit was found

[TABLE]

The corresponding $\zeta_{k}$ are exponentially close to one, $1-\zeta_{k}\sim\mathrm{e}^{-N\eta(\kappa)}$, which is valid up to the centre of the spectrum, $0<\kappa<1/2$, in the scaling limit $k,N\gg 1$ with the ratio $\kappa=k/N$ kept fixed. The scaled eigenvalues are then given by the integral**111see Eqn. (133) in [14]. The quantity $\theta_{k}$ there is $N^{2}(\lambda_{k}/2+1/4\cos{q_{F}})$ due to the different definition of the tridiagonal matrix.00footnotetext: see Eqs. (43),(47) and (59) in [14]. Our $\zeta_{k}$ is called $\lambda_{k}$ there, our $\eta(\kappa)$ corresponds to $L_{3}$ and our $q_{F}$ to the quantity $2\pi W$.

[TABLE]

where the dependence on the spectral parameter $\kappa$ is implicit in the parameter $B=B(\kappa)$. The inverse of this function is given by a similar integral

[TABLE]

where the limits of integration are changed. In particular, the maximal eigenvalue $\kappa=0$ corresponds to $B=1$, where (15) yields $\eta(0)=2\,\mathrm{artanh\,}(1/\sqrt{2})\simeq 1.7627$. The subleading terms in the asymptotics of the extremal eigenvalues can also be determined with the results given in Appendix A. On the other hand, the centre of the spectrum, $\kappa=1/2$, corresponds to $B=0$ and thus trivially $\eta(1/2)=0$. In fact, these low-lying eigenvalues obey a different type of asymptotics, scaling as $1/\ln N$ [14, 18].

In a similar fashion, the eigenvalues $\lambda_{k}$ of $T$ can be analyzed and, in the same uniform asymptotic limit, one finds**111see Eqn. (133) in [14]. The quantity $\theta_{k}$ there is $N^{2}(\lambda_{k}/2+1/4\cos{q_{F}})$ due to the different definition of the tridiagonal matrix.

[TABLE]

Thus the scaling limit $\lambda(\kappa)$ of the eigenvalues is, up to a factor of two, simply given by the parameter $B$. Therefore, instead of expressing $\eta(\kappa)$ and $\lambda(\kappa)$ in terms of the spectral parameter $\kappa$, one can find a direct relation $\eta(\lambda)$ by expanding the integral (15) in terms of $B$. Using the Taylor series

[TABLE]

one rewrites (15) as an infinite sum of integrals

[TABLE]

These integrals can be evaluated explicitly and yield

[TABLE]

Hence, using Eq. (17), one has the series expansion

[TABLE]

As expected, the coefficient of the linear ($m=0$) term is $\alpha_{0}\beta_{0}=\pi$.

Our next goal is to rewrite the relation between the eigenvalues in terms of the corresponding matrices. This can be done by substituting (21) into (12) and (13) and taking the sum over the eigenvalues in the proper scaling limit. Since only odd powers of $\lambda$ appear in (21), the symmetry used in (12) holds for all terms and the relation can immediately be lifted to the matrix level

[TABLE]

The above exact series representation of $h$ already explains the two important qualitative features found numerically in the previous section, namely the deviation from the parabolic profile in the nearest-neighbour hopping, and the appearance of longer-range hopping. Indeed, the $n$-th power of the matrix $T$ generates hopping terms up to the $n$-th neighbour.

To evaluate the matrix elements $h_{i,i+r}$ in the scaling limit, one needs a closed-form expression for the corresponding matrix elements of the powers of $T$

[TABLE]

The summation indices run over the values $l_{j}=1,\dots,N$ for each $j=1,\dots,2m$. However, due to the tridiagonal structure of $T$ with the diagonal being identically zero for half-filling, the only non-vanishing terms in the sum occur for indices that satisfy the constraints

[TABLE]

where we defined extra indices $l_{0}$ and $l_{2m+1}$ that are kept fixed. These constraints define a simple random walk, with the position of the walker given by $l_{j}$, and with the requirement that the walk starts at site $i$ and arrives at $i+r$ in exactly $2m+1$ steps. Since the number of steps is always odd, so must be the distance covered $r=2p+1$, and thus the only non-vanishing matrix elements are $h_{i,i+2p+1}$, as already clear from Eq. (12).

The expression in Eq. (23) further simplifies if one takes a proper continuum limit $i,N\to\infty$ with $i/N$ kept fixed. Indeed, considering the definition of the matrix elements in (8), all the factors $t_{l_{j}}$ with $j=0,\dots,2m$ appearing in (23) can be chosen to be equal. Even though the range of sites visited during the random walk grows as $|i-l_{j}|\leq m-p$, for any finite and fixed values of $m$ and $p$ one has $t_{l_{0}}\simeq t_{l_{1}}\simeq\dots\simeq t_{l_{2m}}$ up to corrections of the order $1/N$ that vanish in the continuum limit. In choosing the proper scaling variable one should, however, take into account the reflection symmetry on the lattice, i.e. $h_{i,i+2p+1}$ should be invariant under $i\to N-2p-i$. The continuum limit of (23) satisfying this constraint then reads

[TABLE]

where the combinatorial factor is just the number of random walks satisfying (24) with $r=2p+1$ and the quantity $i+p$ in $z$ is essentially the midpoint between initial and final site in the hopping.

Finally, the result (25) must be inserted into (22). Interestingly, the infinite sum can be carried out and yields the following closed-form expression

[TABLE]

where $(2p-1)!!=1\cdot 3\cdot 5...(2p-1)$ and ${}_{3}F_{2}$ is a generalized hypergeometric function. Its definition as well as the derivation of the result (26) are summarized in Appendix B. In particular, for the nearest-neighbour hopping one has

[TABLE]

with $x=i/N$. One thus finds that the parabolic profile is multiplied with a generalized hypergeometric function which increases monotonously from $1$ at $x=0$ up to its maximum of approximately $1.076$ in the middle of the interval, $x=1/2$. Hence, despite the continuum limit involved, the lattice result for the nearest-neighbour hopping is different from the CFT prediction, with a relative deviation of approximately $8\%$ around the maximum.

The analytical result (26) for the scaled entanglement Hamiltonian can be tested against direct lattice calculations as in Section 3 and for increasing values of $N$. This is shown in Fig. 3 for the nearest ($p=0$) and third-neighbour ($p=1$) hopping. While for $p=0$ the agreement is excellent already for relatively small $N$, the case $p=1$ already requires to take $N=100$ to obtain a good convergence. This can be attributed to the much smaller value of the matrix element. The finite-size corrections that have been neglected in $(\ref{tcont})$ are then more pronounced. Nevertheless, the data shows a perfect convergence towards the analytical result.

This trend continues also for larger $p$ values. The fifth-neighbour hopping is shown on the left of Fig. 4, and here the largest interval size had to be taken as $N=200$ in order to reach a convincing agreement. Note that the convergence is always the slowest around the peaks. In fact, it is interesting to have a look at the decay of these maxima as a function of the distance $r=2p+1$ in the hopping, i.e. evaluating (26) at the argument $4z=1$. This is shown on the right of Fig. 4 on a logarithmic scale. One can see a clear power-law decay with a power $3$ which was obtained by fitting the data. Formally, it results from an increase of ${}_{3}F_{2}$ with $r$ which is overcompensated by a decrease of the prefactor. While the drop from $r=1$ to $r=3$ is rather big, the following maxima decrease less rapidly. Note that if one plotted the maxima for a fixed value of $N$, the decay would be faster due to the slower convergence for higher $p$.

As to the hopping profiles, the curves are more and more concentrated in the middle of the system as $r$ becomes larger. Near the boundaries, this simply reflects the factor $z^{2p+1}$ in (26) which originates from the lowest power of $T$ leading to the corresponding matrix element. In the centre, the corrections contained in ${}_{3}F_{2}$ enhance the values and lead to a sharpening of the profile.

5 Arbitrary filling

We proceed to the treatment of the general case, where the filling is given by $q_{F}/\pi$ with arbitrary $0<q_{F}<\pi$. We first show the spectra for several values of $q_{F}$ in Fig. 5. The symmetry stated in (9) is clearly visible, even though the largest eigenvalues fall outside the figure. The quantities $-\pi N\lambda_{k}$ are divided here by an additional factor $\sin q_{F}$ which is motivated by the analytics below and seen to produce good agreement with the low-lying $\varepsilon_{k}$. However, the higher $\varepsilon_{k}$ now lie either fully above or fully below the $\lambda$-values which hints at even powers in a series expansion.

The idea of the analytical treatment is the same as before, but the formulas and technical details are more involved. We start again by using the symmetries (9) and (10) to write

[TABLE]

where $\bar{N}\simeq Nq_{F}/\pi$ is the number of $\varepsilon_{k}$ eigenvalues in the negative branch and the hat denotes quantities corresponding to $\pi-q_{F}$. An analogous formula holds for $T_{i,j}$, with the sums taken over the positive branches of $\lambda_{k}$ and $\hat{\lambda}_{k}$ eigenvalues. Introducing the density (13) as in the half-filled case, one has the matrix representations

[TABLE]

The goal is again to relate the above operators in terms of a series expansion. However, in contrast to the half-filled case, one has now $h_{-}\neq\hat{h}_{-}$ and $T_{+}\neq\hat{T}_{+}$ and thus the two contributions have to be treated separately at first. Nevertheless, it is easy to see that the two pieces in the decomposition (29) are orthogonal, $h_{-}(J\hat{h}_{-}J)=0$ and similarly for $T_{+}$, which will be used later.

As before, the first step is to consider the uniform asymptotic limit of the eigenvalues. Instead of (15) and (17), one has now

[TABLE]

Analogously, the general form of the spectral parameter is given by

[TABLE]

where the maximal eigenvalue, $\kappa=0$, corresponds again to $B=1$, whereas setting $B=A$ yields $\kappa=q_{F}/\pi$, i.e. the endpoint of the branches with $\eta(q_{F}/\pi)=\lambda(q_{F}/\pi)=0$. For $B\simeq A$, (30) can be evaluated by setting $\xi^{2}=A^{2}$ in the integrand with the result

[TABLE]

which was already used in plotting the numerical data. Note that the factor $\sin q_{F}$ can be seen as the Fermi velocity in (2) for general filling.

To obtain the function $\eta(\lambda)$ in general, an expansion around $\xi=A$ is, however, not useful, since the resulting series do not converge for all $A$. Instead, we insert again the Taylor series (18) into the integral in (30). This gives the infinite sum (19) with the integrals

[TABLE]

where the coefficients read

[TABLE]

Hence, $\eta$ is now given in terms of a double series in powers of both parameters $A$ and $B$. Using the expression for $\lambda$ in (30), this can be written as

[TABLE]

For half filling, where $A=0$ and only the term $n=0$ contributes, one recovers the previous formula (21). Furthermore, the inversion of the filling, $q_{F}\to\pi-q_{F}$, corresponds to $A\to-A$, and thus an analogous expression with the appropriate sign change holds between $\hat{\eta}$ and $\hat{\lambda}$. The fact that $\eta$ vanishes for $\lambda=0$ is not directly visible but hidden in the coefficients $\beta_{m,n}$, of which $\beta_{m,2m+1}$ is negative.

Lifting (35) to the matrix level requires some extra considerations. Although $\eta(\lambda)$ (respectively $\hat{\eta}(\hat{\lambda})$) can straighforwardly be written as a relation between the matrices $h_{-}$ and $T_{+}$ (respectively $\hat{h}_{-}$ and $\hat{T}_{+}$), this does not immediately yield a relation between $h$ and $T$ as given in Eq. (29). Incorporating the matrix factors $J$ is no problem, but all contributions involving the $T_{+}$ have to combine correctly. Now, expanding in powers in (35), one obtains terms of the form $T_{+}^{k}A^{l}$ where the sum $k+l=2m+1$ is always odd. In the expression for inverted filling, the terms with odd powers of $A$ (and thus even powers of $T_{+}$) will enter with changed sign. Hence, in the sum for $h_{-}-J\hat{h}_{-}J$ one has only the combinations

[TABLE]

where the orthogonality of the two pieces was used. Therefore, the series can indeed be rewritten in terms of $T$ only and one has

[TABLE]

In the final step, we generalize the calculation for the matrix elements to powers of the shifted matrix $\tilde{T}=T+A/2$

[TABLE]

The power $s$ can now take both even and odd values and the summation indices $l_{j}$ must satisfy the following constraints

[TABLE]

Contrary to the half-filled case, now $l_{j}-l_{j-1}=0$ is also allowed since the diagonal elements of $\tilde{T}$ are nonvanishing. Therefore, Eq. (39) defines a random walk, where at each step $j$ the walker is allowed to pause, such that the final site $i+r$ is reached from $i$ in exactly $s$ steps. Denoting by $\ell$ the number of stops, one has in the continuum limit

[TABLE]

where the first binomial factor gives the number of walks reaching the destination site at a distance $r$ in $s-\ell$ jumps, while the second one gives the number of ways the $\ell$ stops can be chosen within $s$ steps. Clearly, the number of jumps $s-\ell$ must have the same parity as the distance $r$, which restricts the possible terms in the sum. Finally, the two factors at the end of (40) are just the diagonal and offdiagonal matrix elements of $\tilde{T}$ in the continuum limit, raised to the power $\ell$ and $s-\ell$, respectively. As for the half-filled case, the scaling variable $z$ must be chosen such that the lattice reflection symmetry $i\to N+1-r-i$ is satisfied for $h_{i,i+r}$, which yields

[TABLE]

One thus arrives at the final expression

[TABLE]

where (40) with $s=2m+1-n$ has to be inserted for the matrix elements. Since this expression is significantly more complicated than the one for the half-filled case, we refrain from a further analysis. Instead, we simply evaluate (42) numerically, using a sufficiently large cutoff $M$ in the infinite sum.

The scaling functions obtained in this way are compared to the lattice results in Fig. 6 for the diagonal ($r=0$, left) and nearest-neighbour hopping terms ($r=1$, right). In both cases we have considered an interval of length $N=40$ and the cutoff was set to $M=30$, yielding an excellent agreement for all the different fillings shown. We have also checked that decreasing the cutoff to $M=10$ has no visible effects on the scaling functions, thus the convergence of (42) is very fast for the dominant matrix elements of $h$. As the diagonal entries $T_{i,i}$ in (8), the matrix elements $h_{i,i}$ are negative for $q_{F}<\pi/2$ and we plotted their absolute value.

It is interesting to have a look also at the longer-range hopping terms, shown in Fig. 7 for $r=2$ and $r=3$, which have again much smaller amplitude. As for half-filling, this requires to take larger intervals for a good agreement with the continuum limit result. We thus used $N=100$ in Fig. 7 and increased also the cutoff to $M=50$, leading again to a perfect match. Note that setting $M=30$ results in tiny changes around the peaks only, barely visible on the scale of the figure. Qualitatively, the hopping for $r=3$ resembles the one for half-filling in Fig. 4. However, the hopping for $r=2$ (which is zero at $q_{F}=\pi/2$) has a double-peak structure.

Returning to the two dominant matrix elements, the question is again, how they compare with a parabolic law. Near the boundaries, i.e. for small $z$, this can be answered very simply. The leading contribution to the hopping term is obtained in (40) if the walker makes only one move, $\ell=s-1$. This gives

[TABLE]

and with the identity

[TABLE]

one finds

[TABLE]

Similarly, the diagonal term is obtained if the walker does not move at all, $\ell=s$, which leads to

[TABLE]

Hence, up to a rescaling by the Fermi velocity, both matrix elements are the same as in $\pi T$. This corresponds to lifting (32) to the matrix level and means that, close to the boundaries, only the small $\eta$ determine the profile.

For the maxima, $z=1/4$, the expression (40) also simplifies. Then the walker is not allowed to pause and only the $\ell=0$ term remains. As a result, one obtains for $r=0$

[TABLE]

while for $r=1$

[TABLE]

These expressions are shown in Fig. 8 for a number of filling parameters $A$. Whereas $h^{\mathrm{max}}_{i,i}$ is antisymmetric in $A$ and vanishes for half filling, $A=0$, the hopping element $h^{\mathrm{max}}_{i,i+1}$ is symmetric and reduces to the expression (27) at $A=0$, since then only the term $n=0$ survives in (48). Both amplitudes diverge as $A\to 1$.

A simple approximation for the two quantities is suggested by looking at the representation (28). If one assumes that the two pieces scale roughly as the corresponding extremal eigenvalues $\gamma N$ and $\hat{\gamma}N$, one can write

[TABLE]

With the explicit formula (52) in Appendix A, one can determine the constants from the lowest terms in $A$ and then evaluate the right hand sides. This gives $c_{0}\simeq 0.532$ and $c_{1}\simeq 0.240$ and leads to the solid lines in the Figures. One sees, that they fit the asymptotic result very well up to $A\sim 0.9$.

Comparing now to the maxima of (45) and (46), which are shown as dotted lines, one sees that the diagonal elements always lie below the parabolic law, while the hopping lies above it for small $A$ (as found for $A=0$ in section 4) and then also below it beyond $A\simeq 0.47$ ($q_{F}\lesssim\pi/3$ or $q_{F}\gtrsim 2\pi/3$) because $1/\sin q_{F}$ diverges more strongly.

6 Discussion

We have determined the entanglement Hamiltonian for an interval in a free-fermion chain. While this can be done with normal numerics for up to $N=20$ sites, larger systems require special routines and precisely for this case we were able to provide an analytical treatment. This was done by writing the single-particle eigenvalues $\varepsilon_{k}$ as a power series in the eigenvalues $\lambda_{k}$ of the matrix $T$, lifting the relation to the matrix level and calculating the elements of $T^{n}$ by counting paths.

As a result, we obtained a rather complete picture of $\mathcal{H}$. For half filling, there is the well-known structure with hopping only between different sublattices, which is largest in the middle of the system. In a contour plot of the matrix $H$, this longer-range hopping shows up as a ridge along the anti-diagonal as found in [16] for a slightly non-critical system. For general filling, there are also diagonal terms and hopping elements within the same sublattice. The latter show a double-peak structure, which results from negative contributions in higher powers of $T$. A general filling corresponds to a critical system with Fermi velocity $v_{F}=\sin{q_{F}}<1$, and one might have expected that a division of $H$ by $v_{F}$ is sufficient to obtain the results from those for half filling. However, while this is true for the low-lying part of the spectrum and the linear part of the hopping profile, the situation is not so simple in general.

The deviations from the conformal prediction should be seen as a lattice effect. The longer-range hopping is not really contained in a conformally invariant system. In [16], such terms were lumped together in order to obtain a local continuum expression and a rather good agreement with a parabolic law was found. It seems, however, that this is connected with the slight off-criticality treated there. In our case, all hopping elements for half filling have the same sign, and combining them would make the deviation from a parabola even larger.

We would like to close with some further remarks. Firstly, the determination of $\mathcal{H}$ differs in a characteristic way from the usual entanglement calculations. For the entanglement entropy, one needs essentially only the low-lying $\varepsilon_{k}$, the higher ones give exponentially small contributions. For $\mathcal{H}$, on the other hand, the high-lying eigenvalues give the largest contributions, and in this sense one is sampling different parts of the spectrum in the two cases. This was pointed out also in [16]. However, the high eigenvectors are oscillator functions located in the centre and do not contribute to the linear variation of the nearest-neighbour hopping near the boundaries. As noted already in section 5, this behaviour comes from the low-lying eigenvalues.

This allows to understand, why the recent approach to interpret the parabolic weight factor in $\mathcal{H}$ as a local temperature and to calculate the entanglement entropy from the thermal equilibrium expression gives the correct asymptotic result [10, 16, 19]. The logarithm $\ln L$ comes from integrating $1/x$, which is the inverse of the linear slope and thus connected with the low-lying part of the spectrum. Closely related is the so-called entanglement contour introduced by Chen and Vidal [20, 21, 22]. In this case, the term $1/x$ can be linked to the envelope $[x(1-x)]^{-1/2}$ of the low-lying eigenfunctions $\phi_{k}$ which enter in the quantity with the largest weight.

Secondly, the commuting matrix $T$ which played a key role in our analysis, is a bit of a mystery, because in contrast to $H$ it has exactly the conformal form. Even more astonishingly, the same situation is found for a finite ring of $M$ sites. Then such a $T$ also exists and has matrix elements which are trigonometric functions [23]

[TABLE]

where, for an odd number of $(2m+1)$ particles, the Fermi momentum is defined as $q_{F}=(2m+1)\pi/M$. This is again the conformal form [10, 8]. Of course, the diagonal terms are only determined up to a global constant, but still the feature is intriguing and one wonders if there is a deeper reason.

Thirdly, the case of an interval at the end of a half-infinite system is closely related to the situation considered here. It can be obtained from a subsystem of $2N+1$ sites in an infinite chain by restricting oneself to the odd eigenfunctions $\phi_{k}$ and the corresponding eigenvalues [15]. This means that the spacing and the maximal eigenvalues are about twice as large as for an interval of length $N$ in the infinite chain. Therefore, with simple numerics, one can only reach a size of $N=10$ here. But a commuting matrix $T_{s}$ exists again and is one-half of the matrix $T$ for $2N+1$ sites so that all steps in the asymptotic treatment can be repeated.

It is a pleasure to thank Ming-Chiang Chung for help and important remarks in the initial stage of this work. We also thank John Cardy for correspondence and Erik Tonni for discussions. V. E. acknowledges funding from the Austrian Science Fund (FWF) through Lise Meitner Project No. M1854-N36.

Appendix A: Extremal eigenvalues

From the results in [14] one can also find asymptotic expressions for the upper and lower ends of the spectra. The $\zeta_{k}$ closest to one gives**222see Eqn. (58) in [14] for $k=0$ and $2\pi W=q_{F}$.

[TABLE]

where

[TABLE]

The largest eigenvalue follows from the relation (9) and is

[TABLE]

where $\hat{\gamma}$ and $\hat{c}$ are obtained from $\gamma$ and $c$ via $q_{F}\to\pi-q_{F}$ which replaces the sines with cosines. Note that $\gamma$, $\hat{\gamma}$ can vary between zero and infinity. For half filling, $q_{F}=\pi/2$, the result is explicitly

[TABLE]

For $N=40$, for example, both subleading terms are approximately 2 while the leading term is 70. The logarithmic correction does not appear in the situation treated in Section 4.

Analogous formulae exist for the eigenvalues of $T$. They can be obtained also by taking a continuum limit and converting the matrix equation into the differential equation for the harmonic oscillator. The result is

[TABLE]

Here, there is no logarithmic correction and the constants cannot exceed one. The total width of the spectrum also equals one to leading order. For half filling, one has

[TABLE]

The expressions (54) and (56) are compared with the numerical data in Fig. 1.

Appendix B: Generalized hypergeometric functions

The function ${}_{3}F_{2}$ encountered in the main text is a special case of the functions ${}_{p}F_{q}$ which generalize the hypergeometric series to more parameters in the coefficients. They are defined as [24]

[TABLE]

where $(a)_{n}=a(a+1)(a+2)...(a+n-1)$ is the Pochhammer symbol with the property

[TABLE]

In particular, $(a)_{1}=a$ and $(1)_{n}=n!$. The series for ${}_{p}F_{q}$ converge absolutely for $|z|<1$ and also for $z=1$ if the sum over the parameters in the denominator is larger than the one in the numerator. For $p=2,q=1$, they reduce to the normal hypergeometric functions.

To obtain (26), one proceeds as follows. The explicit expression for $h_{i,i+2p+1}$ reads

[TABLE]

Since the summation starts at $m=p$, one shifts the index setting $m=p+n$ with $n\geq 0$. This gives

[TABLE]

The second gamma function in the numerator can be rewritten with the doubling formula

[TABLE]

After that, all gamma functions can be expressed by Pochhammer symbols using (58) and one finds

[TABLE]

where the sum is seen to be a ${}_{3}F_{2}$ function and the prefactor

[TABLE]

can be brought into the form given in (26) via the relation $\Gamma(p+1/2)=\sqrt{\pi}\;2^{-p}(2p-1)!!$.

Incidentally, also the relation (21) can be written in terms of a ${}_{3}F_{2}$ function. In this case, the formula is

[TABLE]

The results for arbitrary filling involve sums over generalized hypergeometric functions. For example, (47) is obtained by setting $r=0,z=1/4$ in (40) and changing $n\to 2n+1$, since only odd $n$ are allowed, giving

[TABLE]

For fixed $n$, the sum starts only at $m=n$. Thus, shifting $m\to m+n$ gives as coefficient of $A^{2n+1}$

[TABLE]

and inserting the $\alpha$ and $\beta$, which consist of Gamma functions, one recognizes it as a ${}_{3}F_{2}$ function times some factors. Similarly, the analogue of (64) obtained from (35) is an infinite sum over ${}_{3}F_{2}$ functions.