---
abstract: 'The properties of the reduced density matrix describing an interval of $N$ sites in an infinite chain of free electrons are investigated. A commuting operator is found for arbitrary filling and also for open chains. For a half filled periodic chain it is used to determine the eigenfunctions for the dominant eigenvalues analytically in the continuum limit. Relations to the critical six-vertex model are discussed.'
author:
- |
    Ingo Peschel\
    [Fachbereich Physik, Freie Universität Berlin,]{}\
    [Arnimallee 14, D-14195 Berlin, Germany]{}
title: On the reduced density matrix for a chain of free electrons
---

Reduced density matrices for some portion of a larger system play a central role in the density-matrix renormalization group (DMRG) method [@White92; @White93; @DMRG]. For this reason, they have been investigated for a number of model systems in the last years [@Pescheletal99; @Chung/Peschel00; @Chung/Peschel01; @Cheong/Henley02; @Cheong/Henley03]. It was found that for free electrons or bosons the reduced density matrices have exponential form $exp(-{\cal H})$ and thus look like thermal operators. It was also realized that one can base the determination of $\cal H$ on the one-particle correlation functions in the state one is studying [@Cheong/Henley02; @Peschel03; @Vidal03; @Latorre03]. However, for critical lattice systems, which are particularly interesting in view of DMRG applications, the properties of $\cal H$ have only been studied numerically so far. For free electrons hopping on a chain one can see in this way the typical finite-size effects on the spectra and the concentration of the eigenfunctions near the boundaries. Nevertheless, a more explicit solution, which also makes contact with some results found in field theory [@Gaite03] in calculations of entanglement entropies, would be desirable.\
In the following we present such a solution by looking at the hopping model in a somewhat different way. We show that there exists a relatively simple operator $\cal T$ which commutes with the reduced density matrix for the ground state and thus has the same eigenfunctions. This holds for arbitrary filling and also for a subsystem at the end of a semi-infinite chain. The case of a half-filled periodic chain is then studied in detail. Working with $\cal T$ and taking a proper continuum limit allows to determine the general character of the single-particle eigenstates. For the low-lying states it is found that the eigenvalues of $\cal H$ and $\cal T$ even coincide up to a scale factor. We also discuss the connection with the critical six-vertex model. This allows to derive the spectrum of $\cal H$ by a conformal mapping. It also allows to view the density matrix as a particular transfer matrix and the commutation relation as a special case of similar relations which occur in the treatment of integrable two-dimensional models.\

Basic formulae
==============

We consider a system of free fermions hopping between neighbouring sites of an infinite linear chain. The corresponding Hamiltonian reads =- \_[n]{} (c\_n\^ c\_[n+1]{} + c\_[n+1]{}\^ c\_[n]{}) \[eqn:hop\] where $\hat t$ is the hopping matrix element and the ’hat’ denotes quantities of the total system. The ground state $|0>$ is a Slater determinant corresponding to a certain filling factor $\bar{n}$. As subsystem we take $N$ consecutive sites, labelled by $i,j=1,2,..N$. The reduced density matrix $\rho$ has the general form = [K]{} \[eqn:rho2\] where $\cal K$ is a normalization constant and = \_[i,j]{} H\_[ij]{} c\_i\^ c\_j \[eqn:H1\] The matrix $H_{ij} $ is completely determined by the one-particle correlation function of the total system \_[mn]{} = &lt;0| c\_m\^ c\_n |0&gt; \[eqn:cf\] which in our case is given by $\hat{C}_{mn}=C(m-n)$ where C(m) = \_[-]{}\^ n(q) e\^[iqm]{} \[eqn:cf1\] and $n(q)= 1$ for $|q|< \pi \bar{n}$ and zero otherwise. This leads to C(m)= \[eqn:cf2\] If $C$ is the $N \times N$ submatrix of $\hat{C}$ where the sites are restricted to the subsystem, one has [@Cheong/Henley02; @Peschel03] H = ln \[eqn:HC\] Thus $H$ and $C$ have common eigenfunctions and their eigenvalues $\varepsilon_k$ and $\zeta_k$ are related according to (\[eqn:HC\]). The $\zeta_k$ lie between 0 and 1 while the $\varepsilon_k$ vary between $ -\infty$ and $\infty$. For large $N$, the spectrum of $C$ must approach that of $\hat C$ which is given by the step function $n(q)$. For this reason, most of the $\zeta_k$ lie exponentially close to 0 and 1 [@Latorre03; @Cheong/Henley03] and $C$ is not easy to handle numerically.

Commuting matrices
==================

This difficulty can be circumvented by working with a matrix which commutes with $C$ (and thus $H$), but has a simpler spectrum and also a simpler form. The existence of such a matrix is suggested by the structure of $H_{ij} $ obtained from numerical calculations. One finds that $H_{ij} $ has approximately, but not exactly, tridiagonal form. A similar situation occurs for the transfer matrix of the two-dimensional Ising model, and there a commuting tridiagonal matrix exists [@Suzuki71]. Trying therefore T = (

[ccccc]{} d\_1 & t\_1 & & &\
t\_1 & d\_2 & t\_2 & &\
& t\_2 & d\_3 & t\_3 &\
& & & &\
& & & t\_[N-1]{} & d\_N

) \[eqn:tri\] one finds indeed that $[C,T] = 0$, if the coefficients are chosen as t\_i = , d\_i = - cos(|[n]{}) . \[eqn:coeff1\] The mirror symmetry of the subsystem is reflected directly in these coefficients. Their linear increase near the boundaries is analogous to the situation found for corner transfer matrices [@Baxter82; @Thacker86; @TruongP88; @Davies88]. Physically, the matrix $T$ describes a hopping model with Hamiltonian = -\_[i=1]{}\^[N-1]{} t\_i (c\_i\^ c\_[i+1]{} + c\_[i+1]{}\^ c\_[i]{}) - \_[i=1]{}\^[N]{} d\_i c\_i\^ c\_[i]{} \[eqn:opt\] where both the hopping and the single-site energies increase towards the middle of the subsystem. In the case of half filling ($\bar{n} =1/2$), the diagonal terms are zero and one has only hopping. It should be noted that $T$ is only determined up to a constant factor, and we have chosen a normalization where the largest coefficients are of order one.

A similar result holds if one considers a total system which is semi-infinite and chooses as subsystem the first $N$ sites next to the end. This geometry corresponds to the one in DMRG calculations with open boundaries. In this case, one has to replace $C_{i,j}=C(i-j)$ by the expression $C^s_{i,j} =C(i-j)-C(i+j)$ in all formulae. Then one finds for half filling, that $[C^s,T^s] = 0$ if the $t_i$ are given by t\_i = . \[eqn:coeff2\] The hopping varies again linearly near the interface with the rest of the system and is maximal at the open boundary. Basically, $T^s$ corresponds to the right half of the system described by $T$ which is easy to understand in terms of the geometry. In the following we study $T$ for the case of a half-filled system.

Eigenvectors
============

We now use $T$ to determine the common eigenvectors of $C$, $H$ and $T$. From numerical calculations one sees that for the eigenvalues $\lambda$ of $T$ which are smallest in magnitude the eigenvectors $\varphi$ show rapid oscillations. It is convenient to take these out and also to separate the amplitudes at even and odd sites by writing \_j = (-1)\^[j+1]{} \_[2j-1]{}, \_j = (-1)\^[j+1]{} \_[2j]{} \[eqn:phipsi1\] This leads to the coupled equations -t\_[2j-2]{} \_[j-1]{} + t\_[2j-1]{} \_[j]{} = \_j, t\_[2j-1]{} \_[j]{} - t\_[2j]{} \_[j+1]{} = \_j, \[eqn:phipsi2\] If $N$ is odd, one eigenstate can be found exactly. Then $\lambda=0$ is an eigenvalue and the equations decouple. One finds that $\Psi_j=0$, while $\Phi_j$ follows from the the recursion relation \_[j+1]{} = \_j, \[eqn:rec\] which can be solved explicitly. This “central state” (the $\lambda$ are distributed symmetrically with respect to zero) follows from the structure of the matrix and exists also for $T^s$. It has also been found in the study of Harper’s equation where a similar matrix (with sinusoidal elements) appears [@Hatsugai96]. To see the form of $\Phi$, one can consider large $N$ and take a continuum limit by writing $(2j-1)/N-1/2N= x, t_{2j-1}= t(x),
 \Phi_j= \Phi(x)$. Expanding the quantities, one then finds the differential equation ’(x) = - (x) \[eqn:diff1\] with the solution (x) = = \[eqn:phizero\] which reflects directly the form of the hopping amplitude $t(x)$. This state was studied numerically in [@Cheong/Henley03].\
The same continuum limit can also be used to determine the low-lying eigenstates in general. The equations (\[eqn:phipsi2\]) then become (2t + t’) (x) = (x), -(2t + t’)(x) = (x), \[eqn:diff2\] where $\mu=\lambda N$ absorbs the $1/N$ dependence of the $t_j$. Combining the equations and putting (x) = \[eqn:xi\] one finds for $\chi$ -4 t t = \^2 \[eqn:diff3\] The substitution u = ln(), \[eqn:u\] then leads to the oscillator equation for $\chi(u)$ + \^2 = 0. \[eqn:osc\] The general form of $\Phi$ is therefore (x) = sin \[eqn:sol\] The function $\Psi(x)$ is obtained by replacing the sine with minus cosine according to the relations (\[eqn:diff2\]). Such logarithmic oscillations of the eigenfunctions were also noted in [@Callan94] and had been found earlier for corner transfer matrices [@Davies/Pearce90; @P/Truong91; @P/Wunderling92] where matrices with linearly increasing elements are involved. For the Ising model, the square root appears there, too, and reflects the surface exponent $x_s = 1/2$, while for the Gaussian model there is no such prefactor. To determine the allowed values for $\mu$, one needs proper boundary conditions. From (\[eqn:phipsi2\]) one has $t_1\Psi_1 = \lambda\Phi_1$ and $t_1\Phi_{N/2} = 
 \lambda\Psi_{N/2}$. These lead to the values $\alpha = \pm \pi/4$ for the phase and asymptotically one finds \_k = (2k+1), k= 0,1,2, .... \[eqn:conf\] i.e. equidistant levels scaling as $1/lnN$. For a bosonic field, such a spectrum was also found in [@Callan94]. However, this behaviour can be seen only for very large $N$. For moderate $N$ of the order 10-100 or even 1000, the dispersion relation $\mu_k$ vs. $k$ still shows some curvature [@Chung/Peschel01; @Cheong/Henley03] and the eigenvalues are found to vary as $1/(lnN+b_k)$ rather than as $1/lnN$. The constant $b_0$ for the lowest state, for example, is approximately 2.6. Therefore one has considerable finite-size corrections to the asymptotic law. These features are also known from studies of corner transfer matrices [@P/Truong87; @Davies/Pearce90].

The trigonometric functions $\chi$ vanish near one end of the subsystem and are maximal near the other one. For each eigenvalue, they are related by symmetry such that $\Phi(x) = \pm \Psi(1-x)$.

Direct treatment of C
=====================

With the experience from the treatment of $T$ one can also obtain the eigenvectors directly from $C$. Since $C$ is not a sparse matrix, this leads to an integral equation. Working again with $\Phi$ and $\Psi$, one has instead of (\[eqn:phipsi2\]) the equations \_j [D(2i-2j+1)]{} \_j = | \_i, -\_j [D(2i-2j-1)]{} \_j = | \_i,\
\[eqn:sum1\] where $D(l) = 1/(\pi \;l)$ is the denominator of $C(l)$ and $\bar{\zeta} = \zeta - C(0) = \zeta -1/2$. Combining these, gives a single equation of the form \_j K\_[ij]{} \_j = |[\^2]{} \_i \[eqn:sum2\] Before taking the continuum limit, one has to single out the diagonal terms $K_{ii}$ which are positive, while the others are negative. Since K\_[ii]{} = \_[l=1]{}\^[N/2]{} \[eqn:sum3\] depends only weakly on $i$ and $N$, one approximates it by its limit for large $i$ and $N$, $K_{ii} = K = 1/4$. The equation in the continuum limit then becomes \_[0]{}\^[1]{} dx’ K(x,x’) (x’) = (|[\^2]{}- ) (x) \[eqn:int1\] where K(x,x’) = - . \[eqn:kern\] Using the same substitutions as before, one obtains - \_[-]{}\^ du’ (u’) = (|[\^2]{}- ) (u) \[eqn:int2\] Since one has a difference kernel now, the equation can be solved by Fourier transformation. Writing $\chi$ as a trigonometric function of argument $\mu u$ as in the solution of (\[eqn:osc\]), one finds for $\zeta$ = \[eqn:zeta\] The eigenvalues $\varepsilon$ of $H$ then follow from $\varepsilon = ln((1-\zeta)/\zeta)$ , which gives the simple result = \[eqn:eps\] i.e. up to a factor of $\pi$ they are the same as the eigenvalues of $NT$. For the low-lying states and in the continuum limit, therefore, the matrices $H$ and $T$ are proportional to each other H = N T \[eqn:HT\] This can be checked numerically, and one finds indeed, that $\varepsilon$ and $\lambda$ fulfill this relation with high accuracy. For $N=16$, for example, the relative deviations are of the order of $10^{-3}$ for the lowest states. In general, however, the relation (\[eqn:HT\]) holds only approximately, and the deviations increase for the larger eigenvalues. The close relation of the two matrices can be understood from the structure of $H$ which, to a first approximation, has the same bidiagonal form as $T$ with very similar matrix elements.

Connection to two dimensions
============================

Using the Jordan-Wigner transformation, $\hat{\cal H}$ can be expressed in terms of spin one-half operators and becomes the Hamiltonian of an XX chain = - 2t \_[n]{} (\_n\^x \_[n+1]{}\^x + \_n\^y \_[n+1]{}\^y) \[eqn:xx1\] where $\sigma_n^x,\;\sigma_n^y$ are Pauli matrices. It is known that this operator commutes with the row transfer matrix $\hat{\cal V}$ of the six-vertex model with $\Delta = 0$ [@Sutherland70; @Baxter82]. This model also corresponds to two uncoupled Ising lattices at their critical temperature. Moreover, the ground state of $\hat{\cal H}$ is the state with maximal eigenvalue for the transfer matrix. It follows that one can interpret the reduced density matrix $\rho$ as the partition function of such a system, having infinite size and a cut in it, along which the spins are kept as free variables [@Pescheletal99; @Nishino99]. In previous work on non-critical models, this cut was considered as half-infinite and one was lead directly to standard corner transfer matrices [@Pescheletal99; @Chung/Peschel00]. In the present case, the cut has a finite length $N$.\
This property of $\rho$ can be used to rederive the low-lying spectrum of $\cal H$. For this purpose, one considers the Riemann plane which is obtained by stacking such systems and joining different sides of the cut. If the cut is between $x=-a$ and $x=a$, this is the Riemann plane of the function $ln(z-a)/(z+a)$. For correlations between points on the cut, but in different sheets, $\rho$ plays the role of a transfer matrix, and the asymptotic behaviour of such correlation functions is determined by the largest eigenvalues of $\rho$. Since one is at criticality, one can use the conformal transformation w = ln \[eqn:conf1\] to map the complete Riemann plane to the full $w$ plane. If one further cuts out small circles with radius $b<<a$ around the endpoints $z= \pm a$ of the cut, the map is to a vertical strip $-u_0 < u < u_0$ in the $w$ plane, where $w=u+iv$ and $u_0 = ln(2a/b)$. A vertical correlation function in the strip between points lying $\Delta v =  2\pi$ apart correponds to correlations between equivalent points in consecutive Riemann sheets and involves one power of $\rho$ there. This allows to determine the spectrum of $\rho$ (or $\cal H$) from the correlations in the strip as in [@P/Truong87]. For free boundary conditions along the small circles, this gives the positive $\varepsilon$ as \_k = 2 (x\_s + k) , k= 0,1,2 ... \[eqn:conf2\] where $L = 2 u_0$ is the width of the strip. Inserting the value $x_s= 1/2$ for the surface exponent of the six-vertex model one then has \_k = (2k +1), k= 0,1,2 ... \[eqn:conf3\] If one identifies $2a$ with $N$ and sets $b=1$, this is the continuum result (\[eqn:conf\]), (\[eqn:eps\]) found previously. By choosing $b<1$, one could also obtain an additive correction to $ln N$.

One should mention that a similar mapping was used in [@Holzhey94] to calculate entanglement entropies in the context of black holes. In this case, one needs $\rho$ to calculate the entropy via $S = -tr(\rho\; ln \rho)$. The result was that asymptotically $S$ diverges as $S = (c/3)\;ln N$, where $c$ (equal to one here) is the conformal anomaly. This behaviour was verified recently for various spin chains [@Vidal03; @Latorre03; @Jin03] and can be understood in terms of the logarithmic dependence of $\varepsilon_k$ on the size $N$ [@Callan94]. A discretized spectrum as in (\[eqn:conf3\]) was found previously for finite-size corner transfer matrices. Then $2a$ corresponds to the outer radius of the system and $b$ to the inner one [@P/Truong87].\
The same approach can also be used to discuss $\rho$ for the case of a subsystem at the end of a half-infinite chain. If the end is at $x=0$, one is lead to a half-plane $x > 0$ with a cut from $x=0$ to $x=a$, i.e. from the left edge into the system. This maps to the left half of the previous strip under (\[eqn:conf1\]). Therefore, one has to multiply $\varepsilon_k$ in (\[eqn:conf3\]) by a factor of 2. Such an increase of the single-particle eigenvalues (corresponding to a faster drop of the spectrum of $\rho$ itself) can also be seen in numerical calculations on the lattice and is responsible for the better performance of the DMRG with open boundaries in the present case of a critical system. In the entanglement entropy, this gives a factor of 1/2 which one can interpret as resulting from the reduced number of interfaces between the subsystem and the rest.

Discussion
==========

We have considered the reduced density matrix for a chain of free electrons in its ground state and have found results in two directions. On the one hand, we could derive the structure of the low-lying eigenfunctions which are most important in DMRG calculations and also in general. These have their largest amplitudes near the interface(s) between the subsystem and its surrounding. Away from the interface, they decay with a power law involving the surface exponent $x_s = 1/2$ of the six-vertex model. Thus the influence of the environment propagates far into the interior, as expected for a critical system. The spectrum of $\cal H$ was also obtained and checked by a conformal mapping in the associated two-dimensional model. These results complement previous numerical studies [@Chung/Peschel01; @Cheong/Henley03].\
On a more general level we found the operator $\cal T$, Eqn. (\[eqn:opt\]), which commutes with the reduced density matrix $\rho$. We used it here as a tool for obtaining the eigenfunctions of $\rho$, but it is interesting in itself. Written in terms of Pauli matrices, it has a form analogous to (\[eqn:xx1\]), namely (for half filling) = -2\_[i=1]{}\^[N-1]{} t\_i (\_i\^x \_[i+1]{}\^x + \_i\^y \_[i+1]{}\^y) \[eqn:xx2\] and describes a spatially inhomogeneous XX chain of finite length. The commutation relation $[\rho, {\cal T} ] = 0$ is completely analogous to the relation $[\hat{\cal V},\hat{\cal H}] = 0$ between the row transfer matrix $\hat {\cal V}$ of the six-vertex model and the Hamiltonian $\hat {\cal H}$ of the total chain. It would be interesting, and give further insight, to derive it in a different way. One should mention that reduced density matrices have actually been used in studies of integrable models via vertex operators or via the algebraic Bethe ansatz. In particular, the six-vertex model and the equivalent XXZ spin chain were treated in [@Jimbo96; @Kitanine00]. This lead to exact, but rather complicated expressions for correlation functions along the cut. By contrast, the commuting operator gives a relatively clear physical picture of the problem under consideration.

[99]{}

S. R. White, Phys. Rev. Lett. [**69**]{}, 2863 (1992) S. R. White, Phys. Rev. B [**48**]{}, 10345 (1993) for a review see: I. Peschel, X. Wang, M. Kaulke and K. Hallberg(eds.), [*Density-Matrix Renormalization*]{}, Lecture Notes in Physics Vol. 528, Springer, Berlin (1999) I. Peschel, M. Kaulke and Ö. Legeza, Ann. Physik (Leipzig) [**8**]{}, 153 (1999) M.-C. Chung and I. Peschel, Phys. Rev. B [**62**]{}, 4191 (2000) M.-C. Chung and I. Peschel, Phys. Rev. B [**64**]{}, 064412 (2001) S.-A. Cheong and C. L. Henley, Phys. Rev. B [**69**]{}, 075111 (2004), preprint: cond-mat/0206196 S.-A. Cheong and C. L. Henley, Phys. Rev. B [**69**]{}, 075112 (2004) I. Peschel, J. Phys. A [**36**]{}, L205 (2003) G. Vidal, J. I. Latorre, E. Rico and A. Kitaev, Phys. Rev. Lett. [**90**]{}, 227902 (2003) J. I. Latorre, E. Rico and G. Vidal, Quant. Inf. and Comp. [**4**]{}, 48 (2004), preprint: quant-ph/0304098 for a brief review see J. Gaite, preprint quant-ph/0301120

M. Suzuki, Progr. Theor. Phys. [**46**]{}, 1337 (1971) R. J. Baxter, [*Exactly Solved Models in Statistical Mechanics*]{}, Academic Press, London (1982) H. B. Thacker, Physica D [**18**]{}, 348 (1986) T. T. Truong and I. Peschel, J. Phys. A [**21**]{}, L1029 (1988) B. Davies, Physica A [**154**]{}, 1 (1988) I. Peschel and T. T. Truong, Z. Physik B [**69**]{}, 385 (1987)

Y. Hatsugai, M. Kohmoto and Y.-S. Wu, Phys. Rev B [**53**]{}, 9697 (1996)

C. Callan and F. Wilczek, Phys. Lett. B [**335**]{}, 55 (1994) B. Davies and P. A. Pearce, J. Phys. A [**23**]{}, 1295 (1990) I. Peschel and T. T. Truong, Ann. Physik (Leipzig) [**48**]{}, 185 (1991) I. Peschel and R. Wunderling, Ann. Physik (Leipzig) [**1**]{}, 125 (1992)

B. Sutherland, J. Math. Phys. [**11**]{}, 3183 (1970) T. Nishino, K. Okunishi, Y. Hieida, T. Hikihara and H. Takasaki in Ref. \[3\] Ch. Holzhey, F. Larsen and F. Wilczek, Nucl. Phys. B [**424**]{}, 443 (1994) B.-Q. Jin and V. E. Korepin, J. Stat. Phys. [**116**]{}, 79 (2004), preprint: quant-ph/0304108

M. Jimbo and T. Miwa, J. Phys. A [**29**]{}, 2933 (1996) N. Kitanine, J. M. Maillet and V. Terras, Nucl. Phys. B [**567**]{}, 554 (2000)
