Entanglement Entropy of Eigenstates of Quadratic Fermionic Hamiltonians

TL;DR

This paper develops methods to compute the average entanglement entropy of eigenstates in quadratic fermionic systems, revealing how entanglement scales with subsystem size and system size, and showing eigenstate thermalization in certain limits.

Contribution

It provides exact bounds for the average entanglement entropy of eigenstates in translationally invariant quadratic fermionic Hamiltonians, highlighting deviations from typical pure states.

Findings

Average entanglement entropy is less than maximal for finite fraction subsystems.

Eigenstates exhibit maximal entanglement when the subsystem is a vanishing fraction of the system.

Results demonstrate eigenstate thermalization in the limit of small subsystems.

Abstract

In a seminal paper [D. N. Page, Phys. Rev. Lett. 71, 1291 (1993)], Page proved that the average entanglement entropy of subsystems of random pure states is $S_{\rm ave}\simeq\ln{\cal D}_{\rm A} - (1/2) {\cal D}_{\rm A}^2/{\cal D}$ for $1\ll{\cal D}_{\rm A}\leq\sqrt{\cal D}$, where ${\cal D}_{\rm A}$ and ${\cal D}$ are the Hilbert space dimensions of the subsystem and the system, respectively. Hence, typical pure states are (nearly) maximally entangled. We develop tools to compute the average entanglement entropy $\langle S\rangle$ of all eigenstates of quadratic fermionic Hamiltonians. In particular, we derive exact bounds for the most general translationally invariant models $\ln{\cal D}_{\rm A} - (\ln{\cal D}_{\rm A})^2/\ln{\cal D} \leq \langle S \rangle \leq \ln{\cal D}_{\rm A} - [1/(2\ln2)] (\ln{\cal D}_{\rm A})^2/\ln{\cal D}$. Consequently we prove that: (i) if the subsystem size is a finite fraction of the system size then $\langle S\rangle<\ln{\cal D}_{\rm A}$ in the thermodynamic limit, i.e., the average over eigenstates of the Hamiltonian departs from the result for typical pure states, and (ii) in the limit in which the subsystem size is a vanishing fraction of the system size, the average entanglement entropy is maximal, i.e., typical eigenstates of such Hamiltonians exhibit eigenstate thermalization.

Figures

Equations

$$\mathrm{i}J$$

$$S_m=-\mathrm{Tr}\left\{\left(\frac{11+[\mathrm{i}J{]}_A}{2}\right)\ln \left(\frac{11+[\mathrm{i}J{]}_A}{2}\right)\right\}.$$

$$S_m=V_{\mathrm{A}}\ln 2-\sum _{n=1}^{\infty }\frac{\mathrm{Tr}[\mathrm{i}J{]}_{\mathrm{A}}^{2n}}{4n(2n-1)},$$

$$0\le \mathrm{Tr}[\mathrm{i}J{]}_{\mathrm{A}}^{2(m+1)}\le \mathrm{Tr}[\mathrm{i}J{]}_{\mathrm{A}}^{2m}\le 2V_{\mathrm{A}}$$

$$⟨S⟩=V_{\mathrm{A}}\ln 2-\sum _{n=1}^{\infty }\frac{⟨\mathrm{Tr}[\mathrm{i}J{]}_{\mathrm{A}}^{2n}⟩}{4n(2n-1)},$$

$$V_{\mathrm{A}}\ln 2-\frac{⟨\mathrm{Tr}[\mathrm{i}J{]}_{\mathrm{A}}^2⟩}{2}\ln 2\le ⟨S⟩\le V_{\mathrm{A}}\ln 2-\frac{⟨\mathrm{Tr}[\mathrm{i}J{]}_{\mathrm{A}}^2⟩}{4}.$$

$$\langle\mathrm{Tr}[{\rm i}J]_{\rm A}^{2}\rangle=2\sum_{l=1}^{V}\sum_{i,j\in A}\!\!\begin{array}[]{l}\\ \left(|\alpha_{il}^{*}\alpha_{jl}-\beta_{il}^{*}\beta_{jl}|^{2}\right.\\ \,\,\,\left.+|\alpha_{il}\beta_{jl}-\beta_{il}\alpha_{jl}|^{2}\right).\end{array}$$

$$⟨\mathrm{Tr}[iJ{]}_{\mathrm{A}}^2⟩$$

$$=\frac{2V_{\mathrm{A}}^2}{V}-\sum _{l=1}^V\frac{8|v_{{\vec{k}}_l}{|}^2|u_{{\vec{k}}_l}{|}^2}{V^2}\prod _{\eta =1}^d\frac{{\sin }^2(L_{\eta }[k_l{]}_{\eta })}{{\sin }^2([k_l{]}_{\eta })},$$

$$S_1^{+}$$

$$S_1^{-}$$

$$⟨\mathrm{Tr}[\mathrm{i}J{]}_{\mathrm{A}}^2⟩=\frac{2}{L^2}\sum _{i,j=1}^{L_{\mathrm{A}}}\sum _{k,k^{\prime }=1}^L⟨N_k{N}_{k^{\prime }}⟩e^{\mathrm{i}\frac{2\pi }{L}(k-k^{\prime })(i-j)}=\frac{2L_{\mathrm{A}}^2}{L}$$

$$S_2^{+}$$

$$S_2^{-}$$

$${\Sigma }_S^2=\frac{⟨S^2⟩-⟨S{⟩}^2}{(L_{\mathrm{A}}\ln 2{)}^2}=\frac{1}{(L_{\mathrm{A}}\ln 2{)}^2}\sum _{m,n=1}^{\infty }{F}_{m,n},$$

$$F_{m,n}=\frac{⟨\mathrm{Tr}[\mathrm{i}J{]}_{\mathrm{A}}^{2m}\mathrm{Tr}[\mathrm{i}J{]}_{\mathrm{A}}^{2n}⟩-⟨\mathrm{Tr}[\mathrm{i}J{]}_{\mathrm{A}}^{2m}⟩⟨\mathrm{Tr}[\mathrm{i}J{]}_{\mathrm{A}}^{2n}⟩}{4m(2m-1)4n(2n-1)}.$$

$${\Sigma }_S^2=\frac{1}{L^2}\frac{1}{(\ln 2{)}^2}\left(\frac{L_{\mathrm{A}}}{3}+\frac{1}{6L_{\mathrm{A}}}\right).$$

$${\sigma }^2=\frac{1}{2L_{\mathrm{A}}}⟨\sum _j{\lambda }_j^2⟩=\frac{1}{2L_{\mathrm{A}}}⟨\mathrm{Tr}[\mathrm{i}J{]}_{\mathrm{A}}^2⟩=\frac{L_{\mathrm{A}}}{L}.$$

$$\bar{\mathcal{S}}\equiv \frac{⟨S⟩}{L_{\mathrm{A}}\ln 2}=1-\frac{1}{2\ln 2}\frac{L_{\mathrm{A}}}{L}+\mathcal{O}\left[{\left(\frac{L_{\mathrm{A}}}{L}\right)}^2\right].$$

$$\frac{1}{2}\sum _b{\left(11+\mathrm{i}J\right)}^a{}_b{\hat{\xi }}^b|m⟩=0,$$

$$[\mathrm{i}J{]}_{\mathrm{A}}{\vec{w}}_A=P_{A}J(w_1,\cdots ,w_{2V_A},0,\cdots ,0{)}^{⊺},$$

$$⟨\mathrm{Tr}[\mathrm{i}J{]}_A^{2n}⟩=2\sum _{i_1,\cdots ,i_{2n}=1}^{L_A}⟨j(i_1-i_2)\cdots j(i_{2n}-i_1)⟩,$$

$$\text{contraction}j(d_1)j(d_2)\cdots j(d_{2n-1})j\text{contraction}j(d_1)j(d_2)\cdots jj(d_1)j(d_2)\cdots j(d_{2n-1})j(d_{2n})=c_n\frac{\delta ({\sum }_{i=1}^{2n}{d}_i)}{L^{2n-1}},$$

$$\begin{array}{rl}⟨j(d_1)\cdots j(d_{2n})⟩& =\sum (\text{all possible contractions})\\ & =\text{contraction}jjjj\cdots \text{contraction}jjjj+\cdots +\text{contraction}jj\cdots jj\text{contraction}jj\cdots jjj\cdots jj,\end{array}$$

$$\begin{array}{rl}& ⟨j(d_1)j(d_2)j(d_3)j(d_4)⟩=\text{contraction}jjjj\text{contraction}jjjj+\text{contraction}[2ex]jjjj\text{contraction}jjjjjjj+\text{contraction}[2ex]jjj\text{contraction}jjjjjjjj+\text{contraction}jjjj\text{contraction}jjjjjjj\\ & =c_1^2\frac{\delta (d_1+d_2)\delta (d_3+d_4)}{L^2}+c_1^2\frac{\delta (d_1+d_4)\delta (d_2+d_3)}{L^2}\\ & +c_1^2\frac{\delta (d_1+d_3)\delta (d_2+d_4)}{L^2}+c_2\frac{\delta ({\sum }_{i=1}^4{d}_i)}{L^3}.\end{array}$$

$$⟨\mathrm{Tr}[\mathrm{i}J{]}_A^4⟩=\frac{8L_A^3+L_A}{3L^2}-\frac{2L_A^4}{L^3},$$

$$\begin{array}{rl}& ⟨j(d_1)\cdots j(d_{2n})⟩=\text{contraction}jjjj\cdots \text{contraction}jjjj+\cdots +\text{contraction}jj\cdots jj\text{contraction}jj\cdots jjj\cdots jj\\ & =c_1^n\frac{\delta (d_1+d_2)\cdots \delta (d_{2n_1}+d_{2n})}{L^n}+\cdots +c_n\frac{\delta ({\sum }_{i=1}^{2n}{d}_i)}{L^{2n-1}}.\end{array}$$

$$\mathcal{C}=\prod _{i=1}^l{c}_{n_i}\frac{\delta (D_i)}{L^{2n_i-1}}.$$

$$\sum _{x_1,\cdots ,x_{2n}}^{L_A}\delta (D_1)\cdots \delta (D_l)=\mathcal{O}(L_A^{2n-l+1}),$$

$$\sum _{x_1,\cdots ,x_{2n}}^{L_A}⟨j(d_1)\cdots j(d_{2n})⟩=\frac{\mathcal{O}(L_A^{n+1})}{L^n}+\cdots +\frac{\mathcal{O}(L_A^{2n})}{L^{2n-1}},$$

$$F_{m,n}=\frac{⟨\mathrm{Tr}[\mathrm{i}J{]}_A^{2m}\mathrm{Tr}[\mathrm{i}J{]}_A^{2n}⟩-⟨\mathrm{Tr}[\mathrm{i}J{]}_A^{2m}⟩⟨\mathrm{Tr}[\mathrm{i}J{]}_A^{2n}⟩}{4m(2m-1)4n(2n-1)}.$$

$$\begin{array}{rl}& ⟨j\cdots j|j\cdots jjj\cdots jj⟩\\ & =\text{contraction}j\cdots j|j\text{contraction}[2ex]j\cdots j|j\cdots jj\cdots j|j\cdots j\text{contraction}jjjj\cdots \text{contraction}jjjj+\cdots +\text{contraction}j\cdots j|j\cdots jjj\cdots j\text{contraction}j\cdots j|j\cdots jj\text{contraction}j\cdots j|j\text{contraction}j\cdots j|j\cdots jjj\cdots jjj\cdots j|j\cdots jjj\cdots jj\\ & =\frac{\delta (D_1)\cdots \delta (D_{n+m})}{L^{m+n}}+\cdots +\frac{\delta ({\sum }_{i=1}^{n+m}{d}_i)}{L^{2(n+m)-1}}.\end{array}$$

Full text

Introduction. The concept of entanglement is a cornerstone in modern quantum physics. Different measures of entanglement have been extensively used to probe the structure of pure quantum states eisert_cramer_10 , and they have started to be measured in experiments with ultracold atoms in optical lattices islam_ma_15 ; kaufman_tai_16 . Here, we are interested in the bipartite entanglement entropy (referred to as the entanglement entropy) in fermionic lattice systems. In such systems, an upper bound for the entanglement entropy of a subsystem A (smaller or equal than its complement) is $S_{\rm max}=\ln{\cal D}_{\rm A}$, where ${\cal D}$ and ${\cal D}_{\rm A}$ (${\cal D}_{\rm A}\leq\sqrt{\cal D}$) are the dimensions of the Hilbert space of the system and of the subsystem (see Fig. 1 for an example for spinless fermions). Note that $\ln{\cal D}_{\rm A}\propto V_{\rm A}$, where $V_{\rm A}$ is the number of sites in A, i.e., this upper bound scales with the “volume” of A. (When A is larger than its complement, the Hilbert space of the complement is the one that determines $S$.) Almost 24 years ago, motivated by the puzzle of information in black hole radiation page93a , Page proved page93 that typical (with respect to the Haar measure) pure states nearly saturate that bound lubkin78 ; lloyd_pagels_88 ; foong_kanno_94 ; sanchezruiz95 ; sen96 . Their reduced density matrices are thermal at infinite temperature goldstein_lebowitz_06 ; popescu_short_06 ; tasaki_98 .

In stark contrast with typical pure states, ground states and low-lying excited states of local Hamiltonians are known to exhibit an area-law entanglement eisert_cramer_10 . Namely, their entanglement entropy scales with the area of the boundary of the subsystem. On the other hand, most eigenstates of local Hamiltonians at nonzero energy densities above the ground state are expected to have a volume-law entanglement entropy (with the exception of many-body localized systems nandkishore_huse_review_15 ; altman_vosk_review_15 ). Within the eigenstate thermalization hypothesis (ETH) deutsch_91 ; srednicki_94 ; rigol_dunjko_08 , one expects volume-law entanglement in all eigenstates (excluding those at the edges of the spectrum) of quantum chaotic Hamiltonians santos_polkovnikov_12 ; deutsch13 ; garrison15 ; beugeling15 ; dalessio_kafri_16 , with those in the center of the spectrum exhibiting maximal entanglement dalessio_kafri_16 .

Thanks to the availability of powerful analytical and computational tools to study ground states, many remarkable results have been obtained for the entanglement entropy of such states srednicki_93 ; osterloh_amico_2002 ; osborne_nielsen_02 ; vidal_latorre_03 ; calabrese_cardy_04 ; hastings_07 . On the other hand, for excited states there is a wide gap between what is expected and what has been shown. For interacting Hamiltonians, computational studies are severely limited by finite-size effects so it is difficult to know what happens to the entanglement entropy with increasing the subsystem size. This question was recently addressed for quadratic storms14 ; lai15 ; nandy_sen_16 and nonquadratic but integrable alba15 Hamiltonians, for which one can study much larger lattices, revealing that randomly generated eigenstates are generally maximally entangled in the limit in which the size of the subsystem is a vanishing fraction of the size of the system (in short, a vanishing subsystem fraction).

In this Letter we prove that, for a nonvanishing subsystem fraction, the average entanglement entropy of all many-body fermionic eigenstates of translationally invariant quadratic Hamiltonians departs from Page’s result. Still, the average entanglement entropy exhibits a volume law scaling. In the limit of a vanishing subsystem fraction, we prove that the overwhelming majority of eigenstates are maximally entangled. Our proof stands on calculations of spectrum averages of eigenstate entanglement entropies, which are based on the insight that such averages can be obtained as traces over even powers of a correlation matrix, without the need of calculating its eigenvalues.

Entanglement entropy of energy eigenstates. We study the most general quadratic Hamiltonian of spinless fermions: $\hat{H}=-\sum_{i,j=1}^{V}(\Delta_{ij}\hat{f}_{i}^{\dagger}\hat{f}_{j}^{\dagger}+\Delta_{ij}^{*}\hat{f}_{j}\hat{f}_{i}+t_{ij}\hat{f}_{i}^{\dagger}\hat{f}_{j})$, where $\Delta_{ij}=-\Delta_{ji}$ and $t_{ij}=t_{ji}^{*}$, and $\hat{f}_{i}$ is the fermionic annihilation operator at site $i$. A Bogoliubov transformation $\hat{f}_{i}=\sum_{l=1}^{V}(\alpha_{il}\hat{c}_{l}+\beta_{il}\hat{c}_{l}^{\dagger})$ rotates the Hamiltonian so that it commutes with the quasiparticle number operator $\hat{N}_{l}=2\,\hat{c}^{\dagger}_{l}\,\hat{c}^{\phantom{\dagger}}_{l}-1$. Hence, the many-body energy eigenkets $|m\rangle$ satisfy $\hat{N}_{l}|m\rangle=N_{l}|m\rangle$ with $N_{l}=\pm 1$, and we adopt the binary representation $m=1+\sum_{l=1}^{V}\frac{1+N_{l}}{2}\,2^{l-1}$ ($m$ runs from $1$ to ${\cal D}=2^{V}$, $V$ is the number of lattice sites).

Correlations of a state $|m\rangle$ are encoded in $V\times V$ one-body correlation matrices. They form a $2V\times 2V$ matrix $J$, which is a linear complex structure suppmat

[TABLE]

Since the many-body eigenstates $\{|m\rangle\}$ are Gaussian states, the matrix ${\rm i}J$ fully characterizes them peschel_03 ; peschel_eisler_09 ; dierckx_fannes_08 ; hackl_bianchi_17 . Correlations of a subsystem $A$ containing $V_{\rm A}$ sites are encoded in the restricted complex structure $[{\rm i}J]_{\rm A}$, the $2V_{\rm A}\times 2V_{\rm A}$ matrix obtained by restricting the matrix ${\rm i}J$ to the entries with $i,j\in A$. The entanglement entropy of subsystem $A$ in the eigenstate $|m\rangle$ can be computed as peschel_03 ; peschel_eisler_09

[TABLE]

Expanding Eq. (4) in powers of $[{\rm i}J]_{\rm A}$, about $[{\rm i}J]_{\rm A}=0$, allows one to compute the entanglement entropy without calculating the eigenvalues of $[{\rm i}J]_{\rm A}$

[TABLE]

where we use the compressed notation ${\rm Tr}[{\rm i}J]_{\rm A}^{2n}\equiv{\rm Tr}\{[{\rm i}J]_{\rm A}^{2n}\}$. Since the restricted complex structure satisfies $[{\rm i}J]_{\rm A}^{2}\leq 1\!\!1$ suppmat , one has

[TABLE]

and the series in Eq. (5) is convergent.

Equation (5) allows one to compute the average over the ensemble of all eigenstates $\{|m\rangle\}$ as

[TABLE]

where we define $\langle O\rangle\equiv{\cal D}^{-1}\sum_{m=1}^{\cal D}O_{m}$.

A remarkable property of the series in Eq. (7) is that every higher-order term lowers the average entanglement entropy. Hence, any truncation gives an upper bound. Using the inequality in Eq. (6), one can also produce lower bounds for the average entanglement entropy. To obtain the first-order lower and upper bounds, $S^{-}_{1}\leq\langle S\rangle\leq S^{+}_{1}$, we only need to compute $\langle\mathrm{Tr}[{\rm i}J]_{\rm A}^{2}\rangle$ since: (i) truncating the series in Eq. (7) after the first term results in $S^{+}_{1}$, and (ii) substituting all averages of higher-order traces by $\langle\mathrm{Tr}[{\rm i}J]_{\rm A}^{2}\rangle$ results in $S^{-}_{1}$. This gives

[TABLE]

For a given eigenstate $|m\rangle$ of our Hamiltonian, $[{\rm i}J]_{\rm A}$ is linear in the quantum numbers $N_{l}$. In fact, $\langle m|\hat{f}_{i}^{\dagger}\hat{f}_{j}-\hat{f}_{j}\hat{f}_{i}^{\dagger}|m\rangle=\sum_{l}N_{l}(\alpha^{*}_{il}\alpha_{jl}-\beta_{il}^{*}\beta_{jl})$ and $\langle m|\hat{f}_{i}^{\dagger}\hat{f}_{j}^{\dagger}-\hat{f}_{j}^{\dagger}\hat{f}_{i}^{\dagger}|m\rangle=\sum_{l}N_{l}(\alpha_{il}^{*}\beta_{jl}^{*}-\beta_{il}^{*}\alpha_{jl}^{*})$. The average $\langle{\rm Tr}[{\rm i}J]_{\rm A}^{2n}\rangle$ can therefore be computed from the binomial correlation function $\langle N_{l_{1}}\cdots N_{l_{2n}}\rangle$. In particular, to compute $\langle\mathrm{Tr}[{\rm i}J]_{\rm A}^{2}\rangle$, we use that $\langle N_{l}N_{l^{\prime}}\rangle=\delta_{ll^{\prime}}$ to get

[TABLE]

Whenever $\langle\mathrm{Tr}[{\rm i}J]_{\rm A}^{2}\rangle/V_{\rm A}$ does not vanish in the thermodynamic limit, $\langle S\rangle/V_{\rm A}<\ln{2}$.

Bounds for translationally invariant Hamiltonians. The Bogoliubov coefficients for a translationally invariant system in $d$ dimensions are: $\alpha_{il}=e^{{\rm i}\vec{k}_{l}\cdot\vec{x}_{i}}u_{\vec{k}_{l}}/\sqrt{V}$ and $\beta_{il}=e^{-{\rm i}\vec{k}_{l}\cdot\vec{x}_{i}}v_{\vec{k}_{l}}/\sqrt{V}$, with $u_{\vec{k}}=u_{-\vec{k}}$, $v_{\vec{k}}=-v_{-\vec{k}}$, and $|u_{\vec{k}}|^{2}+|v_{\vec{k}}|^{2}=1$. Substituting these in Eq. (9) leads to

[TABLE]

where the space sum runs within an $d$-dimensional hypercube $A$ with side lengths $L_{\eta}$ containing $V_{\rm A}=\prod_{\eta=1}^{d}L_{\eta}$ sites. One can bound $\langle\mathrm{Tr}[{\rm i}J]^{2}_{A}\rangle$ from below using $|u_{\vec{k}_{l}}|^{2}|v_{\vec{k}_{l}}|^{2}\leq 1/4$. As $V\rightarrow\infty$, one can substitute $\sum_{l}\to\frac{V}{(2\pi)^{d}}(\int^{\pi}_{-\pi}\!\!dk)^{d}$. Since $\int^{\pi}_{-\pi}\!\!dk\sin^{2}(L_{\eta}k)/\sin^{2}(k)=2\pi L_{\eta}$ then $\langle\mathrm{Tr}[{\rm i}J]^{2}_{A}\rangle\geq 2V_{A}^{2}/V-2V_{A}/V$. In the thermodynamic limit, we get that $\langle\mathrm{Tr}[{\rm i}J]^{2}_{A}\rangle=2V_{A}^{2}/V$. The corresponding universal first-order bounds are:

[TABLE]

Note that: (i) $S^{+}_{1}$ and $S^{-}_{1}$ fulfill a volume law as they are proportional to $V_{\rm A}$, and (ii) for any nonvanishing subsystem fraction, $\lim_{V\to\infty}V_{\rm A}/V\neq 0$, $S^{+}_{1}<V_{\rm A}\ln 2$ in the thermodynamic limit, i.e., the average departs from the result for typical states in the Hilbert space. If the subsystem fraction vanishes in the thermodynamic limit, $\lim_{V\to\infty}V_{\rm A}/V=0$, the lower and the upper bounds coincide and $\lim_{V_{\rm A}/V\to 0}S^{-}_{1}=\lim_{V_{\rm A}/V\to 0}S^{+}_{1}=V_{\rm A}\ln 2$. In this limit, the average entanglement entropy is maximal, i.e., typical eigenstates of the Hamiltonian have a typical (à la Page page93 ) entanglement entropy.

Entanglement entropy bounds for free fermions. We now apply our construction to free fermions on a translationally invariant chain with $L$ lattice sites, described by the Hamiltonian $\hat{H}=-\sum_{i,j=1}^{L}t_{j-i}\hat{f}_{i}^{\dagger}\hat{f}_{j}$. In this case, the Bogoliubov coefficients are $u_{k}=1$ and $v_{k}=0$, so that the eigenstates are plane waves. This allows us to obtain closed form expressions for finite systems. We denote the linear subsystem size as $L_{\rm A}$. Figure 1 shows the distribution of $S_{m}$ for all eigenstates $|m\rangle$ in a lattice with $L=36$ sites, as well as the corresponding average $\langle S\rangle$. It is remarkable that when $L_{\rm A}$ departs from 1, the entanglement entropy of the eigenstates with the most weight departs from $S_{\rm max}=L_{\rm A}\ln 2$ storms14 ; nandy_sen_16 .

Using $\langle N_{k}N_{k^{\prime}}\rangle=\delta_{kk^{\prime}}$, we can explicitly compute

[TABLE]

for finite systems. The first-order bounds for $\langle S\rangle$, for finite systems, are then given by Eq. (11) upon replacing $V\to L$ and $V_{\rm A}\to L_{\rm A}$.

It is straightforward to calculate bounds beyond the first order. A general procedure to compute averages of traces of $[{\rm i}J]_{\rm A}^{2n}$ is presented in Ref. suppmat . The main insight from our analysis is that the term $\langle{\rm Tr}[{\rm i}J]_{\rm A}^{2n}\rangle/L_{\rm A}$ is a polynomial that, when $L\to\infty$, only contains powers from $(L_{\rm A}/L)^{n}$ to $(L_{\rm A}/L)^{2n-1}$. For the second order upper ($S^{+}_{2}$) and lower ($S^{-}_{2}$) bounds suppmat , one gets

[TABLE]

In order to obtain Eq. (13), we neglected finite-size corrections of order ${\cal O}(1/L)$ and higher.

In Fig. 2(a), we compare the first- and second-order bounds with the average $\langle S\rangle$ computed on a lattice with $L=36$ sites. The bounds can be seen to be very close to the numerically computed average. At $L_{\rm A}/L=1/2$, where the relative deviation is largest, we get that $0.52<\langle S\rangle/[(L/2)\ln 2]<0.59$. In Fig. S1 of Ref. suppmat , we extrapolate numerical results for $\langle S\rangle$ to the limit $L\to\infty$. Finite-size effects are found to be exponentially small in $L$ iyer15 , and we obtain $\lim_{L\to\infty}\langle S\rangle/[(L/2)\ln 2]=0.5378(1)$.

Entanglement entropy variance for free fermions. In order to understand whether the average of the entanglement entropy over all eigenstates is representative of the entanglement entropy of typical eigenstates, we calculate the variance

[TABLE]

where

[TABLE]

The computation of $F_{m,n}$ is, in general, a daunting task. However, by using a summation technique to compute higher-order traces suppmat , we are able to extract key properties of $\Sigma_{S}$. In particular, we are able to prove that $\Sigma_{S}$ vanishes with increasing the system size as $\Sigma_{S}\sim 1/\sqrt{L}$ or faster suppmat . Furthermore, in the limit of vanishing subsystem fraction (fixed $L_{\rm A}$ for $L\rightarrow\infty$), we obtain the lowest order term in $L$ to be

[TABLE]

Numerical results for $\Sigma_{S}$ in this limit, reported in Fig. 2(b), confirm the accuracy of this prediction. Numerical results for $L_{\rm A}/L=1/2$, reported in Fig. 2(c), confirm that $\Sigma_{S}\sim 1/\sqrt{L}$ for a nonvanishing subsystem fraction. The vanishing of the variance proves that the average and typical entanglement entropies are identical.

Eigenvalue distribution for free fermions. Our results for the average entanglement entropy allow us to unveil some remarkable properties of the eigenvalues $\lambda_{j}$ of $[{\rm i}J]_{\rm A}$ in energy eigenstates. They satisfy $|\lambda_{j}|\leq 1$ suppmat . It is also straightforward to prove that the average of the sum of eigenvalues vanishes: $\langle\sum_{j}\lambda_{j}\rangle=\langle{\rm Tr}[{\rm i}J]_{\rm A}\rangle=L_{\rm A}(1-2\langle N\rangle/L)=0$, where $\langle N\rangle=L/2$ is the average number of particles. On the other hand, the average of the variance of the eigenvalues of $[{\rm i}J]_{\rm A}$ can be calculated using Eq. (12), yielding

[TABLE]

This allows us to conclude that $\sigma^{2}$ vanishes if $\lim_{L\rightarrow\infty}L_{\rm A}/L=0$ (implying $\langle S\rangle$ is maximal), and cannot vanish if $\lim_{L\rightarrow\infty}L_{\rm A}/L\neq 0$. In Fig. 3, we report results of numerical calculations of the distribution of eigenvalues of $[{\rm i}J]_{\rm A}$ (for small values of $L_{\rm A}/L$) in a large ensemble of randomly chosen eigenstates. This distribution can be seen to have a universal form that only depends on the ratio $L_{\rm A}/L$, and whose width is $\sqrt{L_{\rm A}/L}$.

The variance of the distribution of eigenvalues of $[{\rm i}J]_{\rm A}$ ($\sigma^{2}=L_{\rm A}/L$) is important as it determines how the average entanglement entropy reaches the maximal value in the thermodynamic as $L_{\rm A}/L\rightarrow 0$. The lowest order correction to $\langle S\rangle=L_{\rm A}\ln 2$ in terms of $L_{\rm A}/L$ can be read from Eqs. (13), in which the upper and lower bounds coincide up to ${\cal O}[(L_{\rm A}/L)^{2}]$,

[TABLE]

A comparison of the latter expression to numerical results, reported in the inset in Fig. 3, yields an almost perfect agreement for $L_{\rm A}/L\lesssim 0.05$.

For a vanishingly small subsystem fraction, the facts that (i) the average entanglement entropy is maximal, and (ii) the distribution of eigenvalues of $[{\rm i}J]_{\rm A}$ is universal (see Fig. 3), hints that a random ensemble may explain those results. We construct such an ensemble, the Toeplitz Gaussian ensemble (TGE). In the TGE, the entries of $[{\rm i}J]_{\rm A}$ are replaced by random complex numbers whose absolute value is that of a normally distributed variable with zero mean and variance $1/L$, and whose phase is uniformly distributed between 0 and $2\pi$. As shown in Fig. 3, the corresponding eigenvalue distribution is nearly indistinguishable from the numerical calculation over $10^{6}$ random eigenstates. (See also Fig. S2 of Ref. suppmat , which shows that taking the limit $L\to\infty$ first, followed by $L_{\rm A}\to\infty$, results in two distributions that are identical.) This shows that, in the limit of vanishing subsystem fraction, no specific information beyond the symmetries of the model appears to be encoded in the restricted complex structure $[{\rm i}J]_{\rm A}$ of typical eigenstates.

Discussion. Our work introduces a novel methodology that enables the rigorous study of the entanglement entropy of excited eigenstates of quadratic Hamiltonians. The derivation of exact bounds for the average entanglement entropy of translationally invariant quadratic Hamiltonians reveals a fundamental difference between the results for vanishing and nonvanishing subsystem fractions, which is not captured by the analysis of random pure states in the Hilbert space page93 . This highlights the difference in information content between typical eigenstates of physical Hamiltonians, such as those considered here, and typical states in the Hilbert space. The fact that, for a vanishing subsystem fraction, typical eigenstates are maximally entangled constitutes a proof that typical eigenstates satisfy ETH for local observables.

We note that Eqs. (3)–(9) also apply to quadratic models of much current interest such as those appearing in disordered khemani_nandkishore_15 and periodically driven (Floquet) bukov_dalessio_15 systems. While our study focuses on the von Neumann entanglement entropy, the upper bounds derived remain valid for higher-order Renyi entropies, which are bounded from above by the von Neumann entanglement entropy. This is of particular relevance for current experiments with ultracold atoms on optical lattices islam_ma_15 ; kaufman_tai_16 , which can now measure the second Renyi entropy.

Acknowledgments. L.V. and M.R. acknowledge support from the Office of Naval Research, Grant No. N00014-14-1-0540. E.B. acknowledges support from the NSF Grant PHY-1404204. L.H. is supported by a Frymoyer fellowship. This research was supported in part by the Perimeter Institute for Theoretical Physics. The computations were done at the Institute for CyberScience at Penn State.

S1 Linear complex structure

Equation (1) in the main text introduces the linear complex structure $J$ as a way to parametrize eigenstates of quadratic Hamiltonians. $J$ is best understood as a linear map on the vector space $W_{\mathbb{C}}=\mathbb{C}^{2V}$ spanned by $V$ creation operators $\hat{f}_{i}^{\dagger}$ and $V$ annihilation operators $\hat{f}_{i}$. Its eigenvalues are given by $\pm{\rm i}$ with equal multiplicity leading to a decomposition $W_{\mathbb{C}}=W^{+}\oplus W^{-}$ into a direct sum of eigenspaces. The two eigenspaces are related by complex conjugation and define new sets of creation and annihilation operators. This means that every vector in $W^{-}$ represents a linear combination of a new set of annihilation operators, while vectors in $W^{+}$ represent linear combinations of the conjugate creation operators. We can use $J$ to write the projector $\frac{1}{2}\left(1\!\!1+{\rm i}J\right)$ that maps $W_{\mathbb{C}}$ onto $W^{-}$. If we define the vector $\hat{\xi}^{a}=(\hat{f}^{\dagger}_{i},\hat{f}_{j})$ consisting of the $V$ creation and $V$ annihilation operators, we can write the condition

[TABLE]

that defines $|m\rangle$ as the state that is annihilated by all operators in $W^{-}$.

Linear complex structures are used in mathematics to represent the imaginary unit ${\rm i}$ on an even dimensional real vector space. In particular, we have $J^{2}=-1\!\!1$. Complex structures are used in bosonic and fermionic systems to define Gaussian states. They are ideal for the study of entanglement in subsystems because their restriction to a subsystem encodes all relevant correlations.

To compute the entanglement entropy of a subsystem of size $V_{A}$, we need to restrict $J$ to the $2V_{A}\times 2V_{A}$ matrix $[J]_{A}$ associated to the degrees of freedom $\hat{f}_{i}^{\dagger}$ and $\hat{f}_{j}$ in $A$. Let us prove that the eigenvalues $\lambda_{j}$ of the restricted matrix $[{\rm i}J]_{A}$ are in the interval $[-1,1]$. The action of $[{\rm i}J]_{A}$ onto a column vector $\vec{w}_{A}=(w_{1},\cdots,w_{2V_{A}})^{\intercal}$ is

[TABLE]

where $P_{A}(w_{1},\cdots,w_{V})^{\intercal}=(w_{1},\cdots,w_{2V_{A}})^{\intercal}$ is the orthogonal projection (with respect to the standard inner product) onto the first $V_{A}$ components. This implies $\lVert[{\rm i}J]_{A}\vec{w}_{A}\rVert\leq\lVert\vec{w}\rVert$ because: (i) ${\rm i}J$ is norm-preserving being a symmetric matrix with eigenvalues $\pm 1$, and (ii) the orthogonal projection $P_{A}$ cannot increase the norm of a vector. This implies that the eigenvalues $\lambda_{j}$ of $[{\rm i}J]_{A}$ satisfy $-1\leq\lambda_{j}\leq 1$.

S2 Traces of even powers of the complex structure for free fermions

We now focus on free fermions in a one-dimensional periodic lattice with $L$ sites and replace $V\to L$ and $V_{A}\to L_{A}$. In order to compute the higher-order traces

[TABLE]

with $j(d_{i})=\frac{1}{L}\sum^{L}_{k=1}e^{{\rm i}\frac{2\pi kd_{i}}{L}}N_{k}$, we develop a systematic method of computing higher order correlation functions $\langle j(d_{1})\cdots j(d_{2n})\rangle$. Technically, these are Fourier transformed correlation functions of the binomial distribution. They can be computed by adopting a strategy analogous to the one used for computing correlation functions of Gaussian distributions:

The building blocks of correlation functions are the so-called $2n$-contractions given by

[TABLE]

where $\delta(D)=1$ if $D=0\pmod{L}$ and zero otherwise, i.e., $D=\sum^{2n}_{i=1}d_{i}$ is restricted to be an integer multiple of $L$. The prefactors $c_{n}$ can be computed systematically as $c_{n}=L^{2n-1}\langle j(1)^{2n-1}j(1-2n)\rangle$ by using the explicit form $j(d_{i})=\frac{1}{L}\sum^{L}_{k=1}e^{{\rm i}\frac{2\pi kd_{i}}{L}}N_{k}$. For the results reported in the main text, it is sufficient to compute $c_{1}=1$ and $c_{2}=-2$. 2. 2.

Once the $2n$-contractions are known, one can compute a general correlation function as

[TABLE]

where each contraction consists of a product of different pairings, quadruplings, etc., of the $2n$ $j$’s. This is a generalized Wick theorem where we have not only $2$-contractions, but also higher-order $2n$-contractions.

To illustrate this method, let us apply it to the second order correction containing $\langle\mathrm{Tr}[{\rm i}J]_{A}^{4}\rangle$.

S2.1 Second order term

We compute the $4$-point correlation function $\langle j(d_{1})j(d_{2})j(d_{3})j(d_{4})\rangle$. We find three $2$-contractions and one $4$-contraction:

[TABLE]

Plugging this expression into our sum for the trace gives

[TABLE]

where the first term comes from the sum over the $2$-contractions, while the second term comes from the $4$-contraction for which $\delta(\sum^{4}_{i=1}d_{i})=1$ due to $\sum^{4}_{i=1}d_{i}=(x_{1}-x_{2})+\cdots+(x_{4}-x_{1})=0$. Note that, in Eq. (S5), we assume $L_{A}\leq L/2$ to avoid the complication that $\delta(d_{i_{1}}+d_{i_{2}})$ is also nonzero for $d_{i_{1}}+d_{i_{2}}=\pm L$. The term $L_{A}/(3L^{2})$ represents a finite size correction of order $1/L^{2}$. Finite size corrections are expected to appear also at all higher orders, while they are absent at order $n=1$.

Using the result in Eq. (S5) and neglecting finite-size corrections of order ${\cal O}(1/L)$ and higher, we arrive to Eqs. (11) in the main text. For the lower bound, we replace all higher traces by $\langle\mathrm{Tr}[{\rm i}J]_{A}^{4}\rangle$.

S2.2 Higher order terms

To compute higher order traces up to order $n$, we need to determine all prefactors up to $c_{n}$. Here we investigate the scaling in powers of $L_{A}/L$, which can appear in $\langle\mathrm{Tr}[{\rm i}J]_{A}^{2n}\rangle/L_{A}$. A general $2n$-contraction is schematically given by

[TABLE]

To compute $\langle\mathrm{Tr}[{\rm i}J]_{A}^{2n}\rangle/L_{A}$, we sum over $x_{i}$ with $1\leq i\leq 2n$, with $\sum^{2n}_{i=1}d_{i}$ automatically ensured to vanish, which implies that there is one redundant delta in each addend of Eq. (S6).

Let us consider a specific contraction $\mathcal{C}$ of the $2n$ correlation function $\langle j(d_{1})\cdots j(d_{2n})\rangle$ that consists of $l$ $2n_{i}$-contractions with $1\leq i\leq l$. This contraction will give $l$ delta functions $\delta(D_{i})$ with $D_{i}$ being specific sums of the $d_{i}$s according to the chosen contraction:

[TABLE]

This product of deltas appearing in the sum (S3) gives

[TABLE]

where we use the fact that each of the $l$ deltas reduces the dimension of the sum by one, except the one that is redundant. Therefore, the sum gives a polynomial of degree $2n-l+1$. Applying this result to all contractions, we find

[TABLE]

where ${\cal O}(L_{A}^{n+1})$ refers to a polynomial of degree $n+1$. Therefore we conclude that $\lim_{L\to\infty}\langle\mathrm{Tr}[{\rm i}J]_{A}^{2n}\rangle/L_{A}$ is a polynomial function containing only powers from $(L_{A}/L)^{n}$ up to $(L_{A}/L)^{2n-1}$. This result is instrumental in the analysis of the behavior of the entropy for vanishing subsystem fraction $L_{A}/L\to 0$. It implies that it is sufficient to compute the entropy function up to order $n$ to get the dominating terms up to order $(L_{A}/L)^{n}$.

S2.3 Finite-size analysis of the average entanglement entropy for $L_{\text{A}}=L/2$

Finite sums of even powers of the traces in Eq. (S3) provide bounds for the spectrum average of the entanglement entropy $\langle S\rangle$. In Fig. 2(a) of the main text, we compare the bounds in the thermodynamic limit with values of $\langle S\rangle$ in a finite system. Figure S1(a) shows the finite-size scaling of $\langle S\rangle$ for $L_{A}=L/2$. Remarkably, when subtracted by the $L\to\infty$ value ${\cal S}^{\infty}=\lim_{L\to\infty}\langle S\rangle/[(L/2)\ln 2]=0.5378(1)$, we see that the finite-size values approach ${\cal S}^{\infty}$ exponentially fast. This is expected to be the case only in the grand-canonical ensemble in the presence of translational invariance iyer15 .

S3 Scaling of the entanglement entropy variance

Equation (12) in the main text introduces the entanglement entropy variance $\Sigma_{S}$, which contains the term

[TABLE]

Both terms in the numerator of $F_{m,n}$ can be computed as sums over contractions as discussed in Sec. S2. While general contractions appear in the first term, the second term contains no contractions between $\mathrm{Tr}[{\rm i}J]_{A}^{2m}$ and $\mathrm{Tr}[{\rm i}J]_{A}^{2n}$. We can evaluate the difference by computing all contractions of the first term that contain at least one contraction crossing from the first to the second trace, as all contractions without such a crossing will be canceled by the second term. Let us define the symbol $\langle j(d_{1})\cdots j(d_{2m})|j(d_{2m+1})\cdots j(d_{2(m+n)})\rangle$ as the sum over all contractions, for which at least one contraction crosses the separator indicated by the symbol $|$.

The contractions are schematically given by

[TABLE]

In order to find the scaling in $L$, we need to sum over $x_{i}$ coming from the trace expressions, and determine how many deltas are redundant. We find

[TABLE]

where we only have a single redundant delta per term as in the previous case. Even though the two traces imply $\sum^{2m}_{i=1}d_{i}=0$ as well as $\sum^{2(m+n)}_{i=2m+1}d_{i}=0$, which should lead to some contractions with two redundant deltas, these terms are specifically excluded from the sum because they do not contain any crossing contractions. We therefore find that

[TABLE]

is a polynomial function containing only powers $(L_{A}/L)^{m+n}$ up to $(L_{A}/L)^{2(m+n)-1}$. Moreover, in the limit $L\to\infty$, we have

[TABLE]

where $L_{A}^{-1}f_{m,n}(L_{A}/L)$ scales as $1/L$. Finding the coefficient analytically requires the summation of the full series. The coefficient can be computed numerically from a finite-size scaling analysis.

In order to compute the variance scaling for a vanishing fraction $L_{A}/L\to 0$, it is sufficient to compute the first term in Eq. (S14). In fact, our analysis ensures that higher-order terms appear as higher powers in $L_{A}/L$. Truncating the series in Eq. (S14) at $m=n=1$ provides an expansion that is correct up to order $m+n-1=1$ in $L_{A}/L$. At this order, we find

[TABLE]

which is Eq. (14) in the main text, evaluated in Fig. 2(b).

S4 Finite-size analysis of the Toeplitz Gaussian Ensemble

Figure 3 in the main text shows the distribution of eigenvalues ${\cal P}_{\lambda}$ of the restricted complex structure $[iJ]_{A}$. The numerical results for eigenstates are compared to the distribution of eigenvalues ${\cal P}_{\lambda}^{\rm(TGE)}$ of the Toeplitz Gaussian ensemble (TGE). The agreement between the curves is remarkable. Here we quantify their difference. We compute

[TABLE]

where $\delta\lambda=1/100$ is the width used when discretizing ${\cal P}_{\lambda}$.

We first observe that at fixed $L_{A}$, $\delta(L\to\infty)$ extrapolates to a small but finite value, see the inset of Fig. S2 for $L_{A}=20$. We use the fitting function $\delta(L;L_{A})=\Delta(L_{A})+c/L$ to obtain the coefficient $\Delta(L_{A})$. In the second step, see the main panel of Fig. S2, we plot $\Delta(L_{A})$ as a function of $L_{A}^{-1}$. The latter function extrapolates to zero when $L_{A}\to\infty$. In this limit, the eigenvalue distribution of the TGE therefore becomes identical to the eigenvalue distribution of the restricted complex structure $[{\rm i}J]_{A}$.