Strictly local tensor networks for short-range topological insulators

TL;DR

This paper demonstrates that for free fermions, strictly local tensor networks can accurately describe thermal averages of gapped Hamiltonians with short-range couplings, overcoming previous no-go theorems for topologically nontrivial states.

Contribution

It introduces a method to use strictly local tensor networks for thermal states of free fermions with short-range Hamiltonians, including a truncation-reconstruction scheme.

Findings

Tensor networks can describe thermal averages of gapped free fermion systems.

The scheme applies to any dimensionality without topological obstructions.

Illustrated on the 2D Haldane model for both topological and trivial phases.

Abstract

Despite the success in describing a range of quantum many-body states using tensor networks, there is a no-go theorem that rules out strictly local tensor networks as topologically nontrivial groundstates of gapped parent Hamiltonians with short-range (including exponentially decaying) couplings. In this work, we show that for free fermions, strictly local tensor networks may describe nonzero temperature averages with respect to gapped Hamiltonians with exponentially decaying couplings. Parent Hamiltonians in this sense may be constructed for any dimensionality and without any obstructions due to their topology. Conversely, we also show that thermal averages with respect to gapped, strictly short-range free-fermion Hamiltonians can be calculated by tensor networks whose links decay exponentially with distance. We also describe a truncation-reconstruction scheme for such tensor networks that leads to a controlled approximation of exact averages in terms of a sequence of related thermal averages. We illustrate our scheme on the two-dimensional Haldane honeycomb model considering both topological and nontopological phases.

Equations

$$⟨\hat{X}(\hat{a},{\hat{a}}^{†})⟩=N\int (\prod _{j=1}^n\mathrm{d}{\theta }_j\mathrm{d}{\bar{\theta }}_j)\tilde{X}(\theta ,\bar{\theta }){\mathrm{e}}^{-{\bar{\mathbf{\theta }}}^{\top }{\mathsf{h}}_T\mathbf{\theta }},$$

$$N^{-1}=\int (\prod _{j=1}^n\mathrm{d}{\theta }_j\mathrm{d}{\bar{\theta }}_j){\mathrm{e}}^{-{\bar{\mathbf{\theta }}}^{\top }{\mathsf{h}}_T\mathbf{\theta }}=\det (-{\mathsf{h}}_T)$$

$${\mathsf{h}}_T=\coth \left(\frac{\beta }{2}\mathsf{H}\right).$$

$$|f(\mathsf{M}{)}_{ij}|\le C{\chi }^{-R_{ij}/R^{\ast }}.$$

$$⟨\hat{X}(\hat{a},{\hat{a}}^{†})⟩=N\int (\prod _{j=1}^n\mathrm{d}{\theta }_j\mathrm{d}{\bar{\theta }}_j)\tilde{X}(\theta ,\bar{\theta }){\mathrm{e}}^{-{\bar{\mathbf{\theta }}}^{\top }\mathsf{g}\mathbf{\theta }},$$

$$N^{-1}=\int (\prod _{j=1}^n\mathrm{d}{\theta }_j\mathrm{d}{\bar{\theta }}_j){\mathrm{e}}^{-{\bar{\mathbf{\theta }}}^{\top }\mathsf{g}\mathbf{\theta }}=\det (-\mathsf{g}),$$

$${\beta }_{\text{par}}{\mathsf{H}}_{\text{par}}=2\text{arccoth}(\mathsf{g}).$$

$${\mathsf{h}}^{\mathrm{t},\mathrm{r}}=\alpha ({\mathsf{h}}^{\mathrm{t}}-\gamma 1).$$

$$⟨\hat{X}⟩\approx \frac{\int ({\prod }_{j=1}^n\mathrm{d}{\theta }_j\mathrm{d}{\bar{\theta }}_j)\tilde{X}(\mathbf{\theta },\bar{\mathbf{\theta }}){\mathrm{e}}^{-{\bar{\mathbf{\theta }}}^{\top }{\mathsf{h}}^{\mathrm{t},\mathrm{r}}\mathbf{\theta }}}{\int ({\prod }_{j=1}^n\mathrm{d}{\theta }_j\mathrm{d}{\bar{\theta }}_j){\mathrm{e}}^{-{\bar{\mathbf{\theta }}}^{\top }{\mathsf{h}}^{\mathrm{t},\mathrm{r}}\mathbf{\theta }}},$$

$${\eta }_{\pm }^{\mathrm{t},\mathrm{r}}=\coth \left(\frac{\beta {ϵ}_{\pm }}{2}\right),$$

$$S^{\mathrm{t},\mathrm{r}}=-\sum _k[p_k^{\mathrm{t},\mathrm{r}}\ln p_k^{\mathrm{t},\mathrm{r}}+(1-p_k^{\mathrm{t},\mathrm{r}})\ln (1-p_k^{\mathrm{t},\mathrm{r}})],$$

$$⟨{\hat{a}}_j^{†}{\hat{a}}_l⟩=\frac{1}{2}⟨{\delta }_{jl}+{\bar{\theta }}_j{\theta }_l⟩.$$

$$\frac{1}{2}⟨{\delta }_{jl}+{\bar{\theta }}_j{\theta }_l⟩=\frac{1}{2}({\delta }_{jl}-({\mathsf{h}}_T^{-1}{)}_{lj}).$$

$$C_{lj}=⟨{\hat{a}}_j^{†}{\hat{a}}_l⟩=\frac{1}{Z}\text{Tr}\left({\hat{a}}_j^{†}{\hat{a}}_l{e}^{-\beta \hat{H}}\right).$$

$$\mathsf{C}={\mathsf{U}}^{†}\mathsf{DU},$$

$${\mathsf{D}}_{kq}=⟨{\hat{a}}_q^{†}{\hat{a}}_k⟩={\delta }_{kq}⟨{\hat{a}}_k^{†}{\hat{a}}_k⟩$$

$$⟨{\hat{a}}_k^{†}{\hat{a}}_k⟩=\frac{1}{\exp (\beta {ϵ}_k)+1}=\frac{1}{2}\left[1-\tanh \left(\frac{\beta {ϵ}_k}{2}\right)\right],$$

$$\mathsf{C}=\frac{1}{2}\left[\mathsf{I}-{\mathsf{U}}^{†}\tanh \left(\frac{\beta \mathsf{E}}{2}\right)\mathsf{U}\right]=\frac{1}{2}\left(\mathsf{I}-{\mathsf{h}}_T^{-1}\right),$$

$$⟨{\hat{a}}_j^{†}{\hat{a}}_l⟩=\frac{1}{2}\left({\delta }_{jl}-({\mathsf{h}}_T^{-1}{)}_{lj}\right)$$

$$\Gamma_{ij}=i\left\langle\left(\begin{array}[]{c}\hat{c}_{i}^{1}\\ \hat{c}_{i}^{2}\end{array}\right)\left(\begin{array}[]{cc}\hat{c}_{j}^{1}&\hat{c}_{j}^{2}\end{array}\right)\right\rangle-i\mathsf{I}_{2}\delta_{ij}.$$

$$\!\!\!\!-iW^{\dagger}\Gamma_{ij}W\!=\left(\begin{array}[]{cc}\delta_{ij}-2\langle\hat{a}_{j}^{\dagger}\hat{a}_{i}\rangle\\ &2\langle\hat{a}_{i}^{\dagger}\hat{a}_{j}\rangle-\delta_{ij}\end{array}\right).$$

$$-iW^{\dagger}\mathsf{\Gamma}W=\!=\!\left(\begin{array}[]{cc}\tanh\left(\frac{\beta\mathsf{H}}{2}\right)\\ &-\tanh\left(\frac{\beta\mathsf{H}^{T}}{2}\right)\end{array}\right).$$

$$⟨X(\hat{a},{\hat{a}}^{†})⟩=N{P}_0\int \prod _j{T}_j^{(X_j)}\prod _{k<l}{G}_{kl},$$

$$G_{kl}=\exp [-{\bar{\theta }}_{k_l}{\theta }_{l_k}{g}_{kl}]\exp [-{\bar{\theta }}_{l_k}{\theta }_{k_l}{g}_{lk}].$$

$$T_j^{(1)}=\left(\sum _{J\in j}d{\theta }_J\right)\left(\sum _{J\in j}d{\bar{\theta }}_J\right)-g_{jj}.$$

$$T_j^{({\hat{a}}_j)}=-\frac{1}{\sqrt{2}}\sum _{J\in j}d{\bar{\theta }}_J,T_j^{({\hat{a}}_j^{†})}=\frac{1}{\sqrt{2}}\sum _{J\in j}d{\theta }_J,$$

$$T_j^{({\hat{a}}_j^{†}{\hat{a}}_j)}=\frac{1}{2}\left[\left(\sum _{J\in j}d{\theta }_J\right)\left(\sum _{J\in j}d{\bar{\theta }}_J\right)-(g_{jj}-1)\right],$$

$$E_k(\phi ,I)=\inf \left\{\underset{x\in I}{\max }|\phi (x)-p(x)|:p\in P_k\right\},$$

$$x(u)=\frac{{\eta }_2^2-{\eta }_1^2}{2}u+\frac{{\eta }_1^2+{\eta }_2^2}{2}$$

$$E_k(\tilde{f}[x(u)],[-1,1])\le \frac{2M(\chi )}{{\chi }^k(\chi -1)},$$

Full text

I Introduction

It is highly nontrivial to obtain the eigenstates in general quantum many-body problems, even for the ground state. Among the possible approaches, tensor networks TensorNetwork1 ; TensorNetwork2 ; ground have established themselves as promising generalizations of matrix-product-state methods White:DMRG ; Schollwoeck:2011 ; TEBD to systems beyond one spatial dimension. Conceptual and methodological aspects of tensor networks are rapidly developing, with scope including not only condensed matter systems, but also fields such as quantum chemistry Nakatani:2013 , holographic dualitiesholographic , etc. The efficiency and versatility of tensor networks is due to their capability of encoding correlations and entanglement in quantum states such that they effectively capture the relevant Hilbert space based on the area law of entanglement entropyarea-law . Thus tensor networks can be viewed as a natural language for quantum problems where area laws are relevant Orus:2014 .

Tensor networks are parameterized by tensors on discrete sites, and links between them describing the contraction of indices of tensors at the link ends. For a tensor network description to be efficient, the degrees of freedom (DOF) in the tensor variables should not grow exponentially with system size. This, in particular, rules out tensor networks with long-range link structure, since each additional link to a tensor implies an additional tensor index. For gapped systems, the area law, as well as the exponential decay of correlations Hastings:2004 ; Kitaev:2006 suggest some local structure for ground states. One may thus wonder whether a description is possible using local tensor networks, where links between sites are unimportant, or even absent, beyond a certain finite range.

Such local tensor networks, may however be at odds with topologically nontrivial states, as these display characteristic features largely independent of local details. Indeed, even though a range of time-reversal symmetric topological states (e.g. string-net condensates stringnet ) are known to admit local tensor network representations, obstacles exist Wahl:2013 ; Dubail:2015 in the local tensor network description of chiral topological states (e.g. quantum Hall states or chiral $p$-wave superconductors), and more generally for any topologically nontrivial free-fermion state beyond one dimension Read:2017 . Recent no-go theorems Wahl:2013 ; Dubail:2015 ; Read:2017 in fact show that local tensor networks for such topological states can only be the ground states of parent Hamiltonians parent ; Wahl:2013 ; Yang:2015 that are either gapless or have power-law-decaying couplings. (Useful, albeit nonlocal, tensor networks for such states with gapped short-range parent Hamiltonians Iblisdir:2007 ; Feiguin:2008 ; Changlani:2009 ; Motruk:2016 , as well as numerical studies using strictly local tensor networks that account for spurious long-range correlations Chen:2018 , do however exist.)

An alternative approach to the tensor network description of topological phases was suggested in Ref. Beri:2011, . Instead of constructing ground states, it was shown that there exist simple local tensor network representations of groundstate expectation values for a range of topological states, including chiral and strongly correlated cases.

The tensor networks of Ref. Beri:2011, are “exponentially local”: the importance of their links decays exponentially with distance. This is distinct from the “strictly local” tensor networks in the no-go theorems Wahl:2013 ; Dubail:2015 ; Read:2017 mentioned earlier; there links are entirely absent beyond a certain finite distance. Such strictly local tensor networks also naturally arise in applying the approach of Ref. Beri:2011, in practice, since keeping tensor sizes from scaling exponentially with system size requires truncation.

Motivated by these observations, in this work we ask where the tensor networks for expectation values Beri:2011 lead us if we insist that the networks therein be strictly local. Focusing on free fermions as in Refs. Wahl:2013 ; Dubail:2015 ; Read:2017 , we consider two complementary perspectives:

(i) In the spirit of the no-go theorems Wahl:2013 ; Dubail:2015 ; Read:2017 , we take a tensor network as the starting point, and ask about the parent Hamiltonian. We show that a parent Hamiltonian may indeed be defined; with respect to this the tensor network describes nonzero temperature averages. Remarkably, the parent Hamiltonians may be both gapped and have exponentially short-range couplings with no obstructions imposed by their topology.

(ii) Conversely, in the spirit of Ref. Beri:2011, , we take a gapped strictly short-range Hamiltonian $\hat{H}$ as a starting point. We show that, under a suitable truncation-reconstruction procedure, an exact exponentially local tensor network for groundstate averages is approximated by a strictly local network for nonzero temperature averages with respect to a parent Hamiltonian. This process is expected to be exponentially controlled: the elimination of exponentially small couplings should only have small quantitative effects. Reducing the degree of truncation (i.e., increasing the range of permitted network links) corresponds to reducing the fictitious temperature and approaching $\hat{H}$ with the parent Hamiltonian. Throughout this process, the parent Hamiltonian remains gapped; for gapped systems, this truncation-reconstruction procedure thus preserves topological features and connects smoothly to the ground state without significantly affecting physical properties. The same considerations also apply for nonzero temperature averages with respect to $\hat{H}$.

Our results are based on a generalization of the approach of Ref. Beri:2011, to describe the expectation values not just in the ground state but at general nonzero temperatures. We therefore start the presentation, in Sec. II, by explaining the basis of constructing tensor networks for thermal averages with respect to gapped free fermion systems. In Sec. III, we show how suitable mappings between Hamiltonians and such tensor networks relate strictly local to exponentially local objects, including strictly local tensor networks to gapped, exponentially short-range, parent Hamiltonians for thermal averages [(i) above]. In Sec. IV, we discuss how parent Hamiltonians for thermal averages can be reconstructed after truncating an exponentially local tensor network [(ii) above]. We numerically illustrate this in Sec. V taking the two-dimensional (2D) Haldane model Haldane:1988 as an example. We close in Sec. VI with a summary and discussion. The Appendices contain further details on the Grassmann expression underlying our tensor network averages (App. A), topological aspects of thermal states (App. B), and the outline of the proof of the theorem invoked in Sec. III establishing the strictly local to exponentially local correspondence (App. D).

II Grassmann tensor networks and free fermion thermal averages

We consider a free fermion system with quadratic Hamiltonian $\hat{H}=\sum_{i,j}\hat{a}_{i}^{\dagger}H_{ij}\hat{a}_{j}$, where $\hat{a}_{i}^{({\dagger})}$ are fermionic annihilation (creation) operators on a set of $n$ lattice sites labelled by $i=1,\ldots n$. We denote the eigenstates of the single-particle Hamiltonian matrix $\mathsf{H}$ (with components $H_{ij}$ in the basis of lattice sites) by $|\psi_{k}\rangle$ and the corresponding energies (measured from the chemical potential) by $\epsilon_{k}$. The many-body ground state $|\Psi\rangle$ is the Slater determinant of all occupied ($\epsilon_{k}\leq 0$) eigenstates. We will be considering insulators: particle number conserving gapped systems such that $\epsilon_{k}\neq 0$ for all system sizes, i.e., the chemical potential lies within the gap of the system. (Our considerations can be also straightforwardly extended to gapped superconductors.) It will be useful to introduce the single-particle spectral projector $\mathsf{P}=\sum_{\epsilon_{k}<0}|\psi_{k}\rangle\langle\psi_{k}|$ onto the single-particle states occupied in the ground state. The matrix representation of this projector in the lattice basis describes two point groundstate correlations of the lattice operators, $P_{ij}=\langle\Psi|\hat{a}_{j}^{\dagger}\hat{a}_{i}|\Psi\rangle$. It will also be useful to introduce the so-called flat band Hamiltonian $\mathsf{h}=1-2\mathsf{P}=\text{sgn}\,{\mathsf{H}}$ Kitaev:2006 , a matrix obtained from the single-particle Hamiltonian $\mathsf{H}$ by flattening its spectrum to $-1$ ($+1$) for occupied (unoccupied) states, but keeping intact its gapped nature and eigenvectors, and hence its topological features. The tensor networks of Ref. Beri:2011 for groundstate expectation values are based on $\mathsf{h}$.

For nonzero temperatures $\beta^{-1}$, the system can be described by the normalized density matrix $\hat{\rho}=\frac{1}{Z}\mathrm{e}^{-\beta\hat{H}}=\frac{1}{Z}\mathrm{e}^{-\sum_{k}\beta\epsilon_{k}\hat{a}_{k}^{\dagger}\hat{a}_{k}}$ with partition function $Z=\prod_{k}(1+\mathrm{e}^{-\beta\epsilon_{k}})$ and annihilation (creation) operators $\hat{a}_{k}^{(\dagger)}$ of eigenstates $|\psi_{k}\rangle$. Using the same Grassmann integral approach Bravyi:2005 as in Ref. Beri:2011, , we next show that the expressions for nonzero temperature thermal expectation values have the same structure, and thus relation to Grassmann tensor networks Gu:2010 , as those in Ref. Beri:2011, for groundstate averages.

The thermal average $\langle\hat{X}(\hat{a},\hat{a}^{\dagger})\rangle\!=\!\mathrm{Tr}[\hat{X}(\hat{a},\hat{a}^{\dagger})\hat{\rho}]$ of an operator $\hat{X}(\hat{a},\hat{a}^{\dagger})$ can be conveniently expressed as an integral over a set of Grassmann variables. Following Bravyi Bravyi:2005 and Ref. Beri:2011, , the Grassmann mapping substitutes $\frac{\theta_{j}}{\sqrt{2}}$ for $\hat{a}_{j}$, $\frac{\bar{\theta}_{j}}{\sqrt{2}}$ for $\hat{a}_{j}^{\dagger}$, and $\frac{\bar{\theta}_{j}\theta_{j}+1}{2}$ for $\hat{a}_{j}^{\dagger}\hat{a}_{j}$, where $\theta_{j}$ and $\bar{\theta}_{j}$ are Grassmann variables. (We assume that operators on the same site appear in the order $\hat{a}_{j}^{\dagger}\hat{a}_{j}$, without loss of generality.) This maps $\hat{X}(\hat{a},\hat{a}^{\dagger})$ to some function $\tilde{X}(\theta,\bar{\theta})$. As shown in Appendix A,

[TABLE]

where

[TABLE]

is an $X$ independent normalization factor and

[TABLE]

Importantly, Eqs. (1) and (2) are almost identical to the Grassmann integrals that form the starting point of the analysis in Ref. Beri:2011, ; the only difference is the appearance of $\mathsf{h}_{T}$ instead of $\mathsf{h}$. [In the zero temperature limit, $\beta\rightarrow\infty$, $\epsilon_{k}\neq 0$ and $\coth(\pm\infty)=\pm 1$ ensures $\mathsf{h}_{T}\rightarrow\mathsf{h}$.] It is based on these expressions that Ref. Beri:2011, develops a Grassmann tensor network Gu:2010 description and the same approach thus applies to the nonzero temperature case that we consider here. The key step involves Beri:2011 the Grassmann Leibniz rule Zinn-Justin:2002 to change the integral over Grassmann variables on sites to variables on links, with a nonzero matrix element $(\mathsf{h}_{T})_{jl}$ corresponding to there being a link between sites $j$ and $l$. For our considerations below, the transformation into a link-variable structure is not necessary, and thus it will be Eq. (1) that will form the basis of our subsequent analysis. We will, however, use that it is $\mathsf{h}_{T}$ that determines the link structure of the tensor network we could convert to, and thus we shall henceforth refer to $\mathsf{h}_{T}$ as the link matrix. (For completeness, in Appendix C we briefly summarize the explicit relation between our expressions and the Grassmann tensor network contractions of Ref. Gu:2010, .)

III Strictly local tensor networks and their parent Hamiltonians

Having identified $\mathsf{h}_{T}$ as the key object encoding the tensor network connectivity properties, we can now turn to discussing how the locality properties of corresponding Hamiltonian – tensor network pairs are related. The basis of this relation is the following theorem on “strictly short-range” matrices, i.e., matrices $\mathsf{M}$ with indices $i$, $j$ labeling the sites of a finite lattice, and for which there is a range $R^{*}>0$ much smaller than the linear size of the lattice such that $M_{ij}=0$ if $R_{ij}>R^{*}$, where $R_{ij}$ denotes the distance between sites $i$ and $j$:

Theorem 1

Let $\mathsf{M}$ be a Hermitian strictly short-range matrix with range $R^{*}$ and eigenvalues within $D\!=\!D_{-}\!\cup\!D_{+}$ where $D_{-}\!=\![-\eta_{2},-\eta_{1}]$, $D_{+}\!=\![\eta_{1},\eta_{2}]$ with $0<\eta_{1}<\eta_{2}$. Let $f$ be such that on $D_{-}$ ($D_{+}$) it is the restriction of a function $f_{1}$ ($f_{2}$) analytic in a neighborhood of $D_{-}$ ($D_{+}$) in the complex plane. Then, there exist $C>0$ and $\chi>1$ such that

[TABLE]

This result is based on the work of Benzi et al. Benzi:2012 , who proved a variant of it [for $f_{1}$ ($f_{2}$) analytic on the left (right) half of the complex plane] building on approximation theorems of Chui and Hasson Chui:1983 and Bernstein Bernstein . Our formulation of the theorem differs from that of Ref. Benzi:2012, only in qualifying the domain of analyticity of $f_{1,2}$ and is thus proven following nearly the same steps. We outline the proof in Appendix D.

The utility of the Theorem 1 for our purposes lies in taking $\mathsf{M}$ as the strictly short-range object that we wish to relate to its counterpart. Taking $\mathsf{M}=\mathsf{H}$, the gapped strictly short-range single-particle Hamiltonian, the exponentially local nature of the tensor network for groundstate averages Beri:2011 follows if we take $f(x)=\text{sgn}(x)$, that is $f(\mathsf{H})=\mathsf{h}$, the flat band Hamiltonian, which, as we noted in Sec. II, is the groundstate link matrix. (The exponential decay of $\mathsf{h}$ was also noted in Ref. Kitaev:2006, in a form of an approximation theorem; see also Ref. Ringel:2011, .)

If we again take $\mathsf{M}=\mathsf{H}$ but now we use $f(x)=\coth(\beta x/2)$, that is $f(\mathsf{H})=\mathsf{h}_{T}$, we find that the thermal averages remain described by exponentially local tensor networks even for nonzero temperatures.

A particularly interesting implication of Theorem 1 arises when one applies it in reverse: take a tensor network defined via its link matrix $\mathsf{g}$. Requiring the spectrum of $\mathsf{g}$ to lie within the range of $\coth$ ensures that

[TABLE]

with

[TABLE]

corresponds to the thermal average with respect to a parent Hamiltonian $\mathsf{H}_{\text{par}}$ defined by

[TABLE]

We note that this procedure defines $\mathsf{H}_{\text{par}}$, and thus also the temperature $\beta_{\text{par}}^{-1}$, only up to an overall energy scale. However, this already determines the decay properties of the couplings and the presence or absence of a spectral gap in $\mathsf{H}_{\text{par}}$, which are the features we are interested in.

Suppose now that $\mathsf{g}$ defines a strictly local tensor network, and that the spectrum of $\mathsf{g}$ is contained within the union of finite intervals $[-g_{\text{max}},-g_{\text{min}}]\cup[g_{\text{min}},g_{\text{max}}]$, $1\!<\!g_{\text{min}},g_{\text{max}}\!<\!\infty$ such that the number of eigenvalues contained in each of the intervals increases linearly with the number of sites $n$. (That is, viewed as a single-particle Hamiltonian, $\mathsf{g}$ is a strictly short-range finite bandwidth insulator with spectrum in the analytic domain of arccoth.) Then, by Theorem 1 and Eq. (7), $\mathsf{H}_{\text{par}}$ will be an insulator with exponentially decaying couplings. Furthermore, the eigenvectors of $\mathsf{g}$ and $\mathsf{H}_{\text{par}}$ are the same; in particular the eigenvectors of $\mathsf{g}$ corresponding to eigenvalues $G_{j}$ of a given sign are eigenvectors of $\mathsf{H}_{\text{par}}$ with eigenvalues of the same sign. Therefore, if the $G_{j}<0$ eigenvectors of $\mathsf{g}$ have topological features, including being those of a chiral topological state with nonzero Chern number, then $\mathsf{H}_{\text{par}}$ will define a ground state, and thermal density matrix, with the same topological features. (The density matrix topology follows from the same eigenvectors rhotop1 ; rhotop2 , as we review in Appendix B.) We thus find that thermal expectation values of gapped parent Hamiltonians with exponentially decaying couplings can be represented by strictly local tensor networks, without any restrictions imposed by topology.

It is useful to compare this finding to the no-go theorems of Refs. Dubail:2015, ; Read:2017, : these rule out the existence of gapped parent Hamiltonians with short-range (including exponentially decaying) couplings for a large class of strictly local tensor networks encoding topologically nontrivial ground states. Our construction, in contrast, works with different settings: firstly we consider expectation values instead of a direct representation of states. Secondly, at least for topologically nontrivial insulators $\mathsf{H}_{\text{par}}$, the strictly local tensor networks in our approach always correspond to thermal averages at strictly positive temperature. This is not a mere artefact of using Eq. (7) to define $\mathsf{H}_{\text{par}}$: for any $\mathsf{H}_{\text{par}}$, zero temperature averages correspond to $\mathsf{g}=\text{sgn}\,\mathsf{H}_{\text{par}}$, which cannot be both topologically nontrivial and strictly short-range Chen:2013 ; Read:2017 .

IV Truncation and Reconstruction

Starting from an exponentially local tensor network, a strictly local tensor network naturally arises if one wishes to perform a numerical computation: In a practical calculation, to have a numerically tractable number of DOF, one has to truncate some of the longer links, be the tensor network constructed from the flat band Hamiltonian $\mathsf{h}$ in the zero temperature case Beri:2011 , or from $\mathsf{h}_{T}$ at nonzero temperatures. We shall refer to the resulting truncated link matrix as $\mathsf{h}^{\rm t}$. If the goal is simply to obtain a controlled approximation for the purposes of numerical results, one may obtain $\mathsf{h}^{\rm t}$ by directly discarding links beyond a certain cutoff length; the error will be controlled exponentially due to the exponentially local nature of the exact network. However, more care should be taken if one wants to associate a valid free-fermion state to the network after truncation, because merely truncating $\mathsf{h}$ or $\mathsf{h}_{T}$ may place some of the eigenvalues $\eta^{\rm t}_{k}$ of $\mathsf{h}^{\rm t}$ in the interval $(-1,1)$, thus placing the spectrum of $\mathsf{h}^{\rm t}$ outside the domain of arccoth [cf. Eq. (7)] thereby precluding the correspondence to a parent Hamiltonian, and thus parent density matrix $\hat{\rho}=Z^{-1}\exp\left[-\beta\sum_{i,j}\hat{a}_{i}^{\dagger}(\mathsf{H}_{\text{par}})_{ij}\hat{a}_{j}\right]$.

Whether $|\eta^{\rm t}_{k}|<1$ arise depends both on the degree of truncation and on the temperature. For a thermal system at high temperature (small $\beta$) the eigenvalues $\eta_{k}=\coth({\beta}\epsilon_{k}/2)$ of $\mathsf{h}_{T}$ are far from $[-1,1]$ if the $\epsilon_{k}$ spectrum is bounded. One then expects $|\eta^{\rm t}_{k}|>1$ to continue to hold under any moderate truncation. However for the zero temperature system, represented by the flat-band Hamiltonian $\mathsf{h}$ with eigenvalues $\pm 1$, getting $|\eta^{\rm t}_{k}|<1$ for at least some eigenvalues is inevitable however small the truncation is.

To address the cases where $|\eta^{\rm t}_{k}|<1$, we modify the truncated link matrix $\mathsf{h}^{\rm t}$ into a “reconstructed” link matrix $\mathsf{h}^{\rm t,r}$ by a shift and rescaling

[TABLE]

Here $\alpha$ and $\gamma$ are chosen to ensure that the eigenvalues of this truncated and reconstructed link matrix, $\eta^{\rm t,r}_{k}=\alpha(\eta^{\rm t}_{k}-\gamma)$, satisfy $|\eta^{\rm t,r}_{k}|\geq 1$ (see below). As $\mathsf{h}^{\rm t,r}$ has a spectrum in the domain of arccoth, it now corresponds to an effective parent Hamiltonian $\beta^{\rm t,r}{{\mathsf{H}^{\rm t,r}}}$ with eigenvalues $\beta^{\rm t,r}\epsilon^{\rm t,r}_{k}=2\,\mathrm{arccoth}\eta^{\rm t,r}_{k}$. After this reconstruction we obtain the approximate Grassmann formula

[TABLE]

which replaces $\mathsf{h}_{T}$ by $\mathsf{h}^{\rm t,r}$ in Eqs. (1) and (2).

To specify $\alpha$ and $\gamma$, we impose two constraints, one each on $\eta^{\rm t,r}_{+}$ and $\eta^{\rm t,r}_{-}$, the minimum positive and the maximum negative eigenvalue of $\mathsf{h}^{\rm t,r}$, respectively. These, via arccoth, correspond to the largest (in magnitude) positive and negative eigenenergies, $\epsilon^{\rm t,r}_{+}$ and $\epsilon^{\rm t,r}_{-}$, respectively, of $\mathsf{H}^{\rm t,r}$. To set $\eta^{\rm t,r}_{\pm}$, we choose $\beta^{\rm t,r}\epsilon^{\rm t,r}_{\pm}=\beta\epsilon_{\pm}$, that is,

[TABLE]

where $\epsilon_{+}$ and $\epsilon_{-}$ are maximum and minimum eigenenergies of $\mathsf{H}$, respectively, and $\beta$ is the physical inverse temperature. For zero temperature, as $\beta\epsilon_{\pm}\rightarrow\pm\infty$, the constraints become $\eta^{\rm t,r}_{\pm}=\pm 1$. We note that the corresponding infinite thermal eigenvalues $\beta^{\rm t,r}\epsilon^{\rm t,r}_{\pm}$ only indicate the corresponding states of the approximate nonzero temperature problem always being occupied or empty, but not a divergent approximation for the total energy $\langle\hat{H}\rangle$, which should be calculated via Eq. (9).

We also note that instead of Eq. (10), one could also place the $\eta^{\rm t,r}_{k}$ spectrum into the domain of arccoth via a requirement on $\min(|\eta^{\rm t,r}_{-}|,|\eta^{\rm t,r}_{+}|)$ only. As this is just one constraint, the parameters $\alpha$, $\gamma$ may be chosen such that an additional requirement is satisfied, e.g., that the exact average particle number is maintained under the truncation reconstruction procedure: $\langle\hat{N}\rangle\!=\!\sum_{k=1}^{n}\frac{1}{\mathrm{e}^{\beta\epsilon_{k}}+1}$ be equal to $\langle\hat{N}^{\rm t,r}\rangle\!=\!\sum_{k=1}^{n}\frac{1}{\mathrm{e}^{2\mathrm{arccoth}\eta^{\rm t,r}_{k}}+1}$.

With the truncation-reconstruction procedure described above, a remaining question is whether the parent Hamiltonian $\mathsf{H}^{\rm t,r}$ is gapped. The gap of $\mathsf{H}^{\rm t,r}$ closes if the spectrum of $\mathsf{h}^{\rm t,r}$ approaches $\pm\infty$. Provided the spectrum of $\mathsf{h}_{T}$ is bounded (as is the case for insulators), no gap closing of $\mathsf{H}^{\rm t,r}$ will occur as the changes due to truncation (and thus the corresponding adjustments in reconstruction) are exponentially controlled. We expect the gapless case to arise only near topological phase transitions, where the gap of the physical Hamiltonian $\mathsf{H}$ vanishes.

V Numerical results for the Haldane model

We now turn to studying our truncation reconstruction procedure on a concrete model. The 2D Haldane model on honeycomb lattice is a minimal, two-band, system permitting nontrivial topological phases Haldane:1988 . Its Hamiltonian involves on-site energy offsets $\pm M$ for A and B sublattices, nearest neighbour hoppings with uniform amplitude $t_{1}$, and next-nearest neighbour hoppings with coefficients $t_{2}\mathrm{e}^{\pm\mathrm{i}\phi}$ where the positive (negative) sign corresponds to the hoppings forming a closed triangle in one hexagon with clockwise (counterclockwise) direction. This model is known to yield Chern numbers $C=\pm 1$ when $|M|<3\sqrt{3}|t_{2}\sin\phi|$ (cf. insets in Fig. 1). We utilize this model to demonstrate the truncation and reconstruction approach described above for gapped systems, including its applicability to topological phases. For comparison, we also show results for a gapless system arising at a topological phase transition.

In what follows we focus on the approximation of groundstate expectation values, i.e., the truncation starting from $\mathsf{h}$, where reconstruction is always necessary as discussed before. For simplicity, we choose a lattice of rhombus shape with periodic boundary conditions.

We can obtain the groundstate expectation value of any operator using Eq. (1) with $\mathsf{h}_{T}$ replaced by the flat-band link matrix $\mathsf{h}$. As $\mathsf{h}$ is only exponentially but not strictly local, we wish to truncate it into $\mathsf{h}^{\rm t}_{m}$ with various truncation lengths $L_{m}$ by setting all entries $h_{ij}\!=\!0$ if $R_{ij}\!>\!L_{m}$. (We choose $L_{m}$ as the $m$-th nearest neighbor distance.) Then we reconstruct to $\mathsf{h}^{\rm t,r}_{m}$, which is not only gapped around $(-1,1)$ by construction, but, as expected and verified numerically, also has finite bandwidth if $\mathsf{H}$ is gapped, thus corresponding to a gapped $\mathsf{H}^{\rm t,r}$.

For the simplicity of reconstruction, we consider the system at half filling; here $\mathsf{H}$ is particle-hole symmetric. The truncation procedure is chosen such that the same particle-hole symmetry holds for $\mathsf{h}^{\text{t}}$; this implies $\gamma\!=\!0$. At half filling, this two-band model has average particle number $N_{\mathrm{tot}}$ for a system with $N_{\mathrm{tot}}=N\times N$ unit cells, where $N$ is the period along one edge of the rhombus; the preservation of spectral symmetry during truncation ensures that $\langle\hat{N}^{\rm t,r}\rangle=N_{\mathrm{tot}}$ is maintained. We find that for gapped systems, the finite correlation length renders the $N$ dependence in our numerics negligible beyond a certain system size; the results in Figs. 1 and 2 are for $N=40$ which is already in this $N$ independent regime.

We consider a number of measures to assess the effects of truncation. Since the spread and the location of the $\eta^{\rm t}_{k}$ spectrum determines the properties of the subsequent reconstruction, in Fig. 1 (left panel) we show the smallest and largest positive $\eta^{\rm t}_{k}$ as a function of truncation length $L_{m}$. (The behavior of the negative eigenvalues follows by spectral symmetry.) The rate of $\eta^{\rm t}_{k}$ deviating from $\pm 1$ as $L_{m}$ is reduced depends on how far the system is from the topological phase transition. The closer to the transition, the wider the $\eta^{\rm t}_{k}$ spectrum becomes, as truncations result in larger errors due to $\mathsf{h}$ becoming less local. Due to $\gamma\!=\!0$, and spectral symmetry, the reconstruction amounts to $\eta^{\rm t,r}_{k}=\alpha\eta^{\rm t}_{k}$ with $\alpha=1/\eta^{\rm t}_{+}$. The $\eta^{\rm t,r}_{k}$ spectrum is thus a further factor of $1/\eta^{\rm t}_{+}$ wider than that of $\eta^{\rm t}_{k}$.

Next we turn to the entropy corresponding to $\beta^{\rm t,r}\mathsf{H}^{\rm t,r}$, which provides a useful single-number measure of how thermal the state corresponding to $\beta^{\rm t,r}\mathsf{H}^{\rm t,r}$ is and thus how far it is from the ground state of the physical Hamiltonian. We have

[TABLE]

where $p^{\rm t,r}_{k}=[\exp(\beta^{\rm t,r}\epsilon^{\rm t,r}_{k})+1]^{-1}$ and the summation is over all eigenstates. The results for $S^{\rm t,r}$ (per unit cell) are shown in Fig. 1 (right panel). Similarly to the spectrum plots, the truncation effect is smaller for cases deep in the gapped (topological or nontopological) phases but larger near a phase transition. It is noticeable that $S^{\rm t,r}$ grows from [math] with decreasing $L_{m}$ faster than the spectra do in Fig. 1, because the entropy formula increases dramatically from zero with $\eta^{\rm t,r}_{k}$ moving away from $\pm 1$. In the extreme case at the topological phase transition, e.g., with $M\!=\!0.5$, $\phi\!=\!\pi/6$, $S^{\rm t,r}$ approaches the maximum for such a two-site unit cell, i.e., $2\ln 2\!\approx\!1.386$.

Finally, to obtain a spatially resolved measure of the effects of truncation, we turn to the one-body density matrix $\langle\hat{a}_{i}^{\dagger}\hat{a}_{j}\rangle$, focusing on the gapped case where our approach is expected to work. Before truncation, as we are considering the ground state, $\langle\hat{a}_{i}^{\dagger}\hat{a}_{j}\rangle=P_{ji}\!=\!h_{ji}/2$ ($i\!\neq\!j$). Calculating the same after truncation and reconstruction using Eq. (9) leads to the approximate result $\langle\hat{a}_{i}^{\dagger}\hat{a}_{j}\rangle^{\rm t,r}\!\equiv\!P^{\rm t,r}_{ji}=\![(\mathsf{h}^{\rm t,r})^{-1}]_{ji}/2$ ($i\!\neq\!j$). Given that $\gamma\!=\!0$ in the particular case we are studying, the change in the off-diagonal correlations after reconstruction is a simple scaling: $P^{\rm t,r}_{i\neq j}\!=\!\alpha^{-1}[(\mathsf{h}^{\rm t})^{-1}]_{i\neq j}/2$. In Fig. 2, we show the approximate correlations $P^{\rm t,r}_{ij}$ and the corresponding errors, for the topologically nontrivial case of Fig. 1 with $M=0.5$, $\phi=\pi/2$, and Chern number $C=1$. (The topologically trivial insulator at $\phi=0$ has similar behavior.) For both $P^{\rm t,r}_{ij}$ and the errors, the absolute values shown are isotropic, i.e., the spatial dependence is only through $R_{ij}$. As the $L_{m}$ dependence of the errors shown illustrates, the approximation converges uniformly and exponentially upon increasing the truncation length; the short-distance correlations are in the exponentially converging regime already for moderate $L_{m}$. The accurate recovery of the one-body density matrix confirms that the Chern number $C$ is preserved throughout the truncation-reconstruction: owing to the exponential decay of correlations, $C$ is encoded already in a few, short-range elements $P_{ij}$ Kitaev:2006 ; Ringel:2011 (a conclusion that remains true for thermal states as characterized in App. B and references therein); this also allows the calculation of the Chern number via local Chern markers Bianco:2011 ; Caio:2019 . The preservation of $C$ is also as expected from the absence of gap closing during our procedure.

VI Summary and Discussion

We have formulated the calculation of thermal expectation values in free fermion systems in a form equivalent to Grassmann tensor networks, thereby generalizing the results of Ref. Beri:2011, to nonzero temperatures. The consideration of nonzero temperature averages allowed us to demonstrate how strictly local tensor networks can describe gapped, exponentially short-range parent Hamiltonians. Conversely, a tensor network for thermal averages was shown to be exponentially local provided we work with a gapped, strictly short-range Hamiltonian. Topology presents no obstructions to these results; they only require the corresponding spectrum to lie in the analytic domain of the function mapping between the link matrix and the parent Hamiltonian Eq. (7) or vice versa, Eq. (3). The parent Hamiltonians here, however, are Hamiltonians with respect to which nonzero temperature expectation values are calculated using the tensor networks, thus are understood in a slightly different sense than parent Hamiltonians in earlier works Wahl:2013 ; Dubail:2015 ; Read:2017 where tensor networks encode the corresponding ground states.

We have also described a truncation and reconstruction procedure applicable to the exact, exponentially local, tensor networks above, be those for thermal or groundstate averages. The strictly local tensor networks arising in this process approximate exact averages by averages with respect to a sequence of parent Hamiltonians at fictitious nonzero temperatures; the corresponding thermal states approach their exact form upon increasing the truncation length. We have illustrated this procedure by approximating groundstate averages for the 2D Haldane honeycomb model and found that the truncation-reconstruction works well both for nontopological and topological gapped phases.

While we formulated our analysis in terms of thermal states of Hamiltonians, our approach, at its core, only uses the Gaussian nature of the state of interest and therefore may be extended (by suitable choices for $f$ in Theorem 1) to include Gaussian states arising, e.g., from dissipative dynamics rhotop1 ; rhotop2 . In that context, a trade-off analogous to that observed here, namely between strictly short-range dissipative dynamics, nontrivial topology, and the purity of the achievable states, has been pointed out recently Goldstein:2018 ; the reason there is the same as for tensor-networks (cf. Sec. III): the impossibility for the flat band Hamiltonian $\mathsf{h}$ to be both strictly short-range and topologically nontrivial Chen:2013 ; Read:2017 .

Our results represent the first step towards developing efficient approaches bridging tensor networks, mixed Gaussian states, and topology. Key to our approach is exploiting the advantages offered by focusing on expectation values directly and bypassing representations for states themselves. This initial study based on simple free-fermion Hamiltonian leaves much to be investigated in the future. One of the most interesting questions is to what extent our approach can be generalized to capture and approximate strongly correlated states.

Acknowledgements

This research was supported by the NSFC grant No. 11704072, the EPSRC grants EP/M02444X/1, EP/K030094/1, EP/P034616/1 and EP/P009565/1, a Simons Investigator award, the Royal Society, and the ERC Starting Grant No. 678795 TopInSy.

Appendix A Thermal expectation values via Grassmann integrals

In this Appendix we verify the Grassmann integral expression Eq. (1). The averages on both sides of Eq. (1) can be calculated using Wick’s theorem. The substitution rules $\hat{a}_{j}\rightarrow\frac{\theta_{j}}{\sqrt{2}}$, $\hat{a}_{j}^{\dagger}\rightarrow\frac{\bar{\theta}_{j}}{\sqrt{2}}$, and $\hat{a}_{j}^{\dagger}\hat{a}_{j}\rightarrow\frac{\bar{\theta}_{j}\theta_{j}+1}{2}$ from the operator $\hat{X}(\hat{a},\hat{a}^{\dagger})$ to $\tilde{X}(\theta,\bar{\theta})$ thus imply the following correspondence between two-point averages:

[TABLE]

Eq. (1) holds if calculating the left hand side of Eq. (12) with the thermal density matrix gives the same result as calculating the right hand side via the Grassmann integral expression.

For the right hand side, the Grassmann integral yields

[TABLE]

To evaluate the left hand side, we define

[TABLE]

In terms of the matrix $\mathsf{U}$ that diagonalizes the single-particle Hamiltonian as $\mathsf{H}=\mathsf{U}^{\dagger}\mathsf{E}\mathsf{U}$ with $\mathsf{E}=\text{diag}(\epsilon_{k})$, we have

[TABLE]

with

[TABLE]

in terms of the annihilation (creation) operators $\hat{a}_{k}^{(\dagger)}$ of single-particle eigenstates. Using

[TABLE]

we find

[TABLE]

or

[TABLE]

as required.

Appendix B Topological characterization of free fermion thermal states

In this Appendix, we briefly summarize how the topological features of free fermion thermal states may be characterized, adopting the approach of Ref. rhotop2, to the particle number conserving case we consider here. The topological characterization of Ref. rhotop2, is in terms of the correlation matrix $\Gamma_{ij}^{pq}=(i/2)\langle[\hat{c}_{i}^{p},\hat{c}_{j}^{q}]\rangle$, where $i,j=1,\ldots,n$, $p,q=1,2$ with Majorana fermions $\hat{c}_{j}^{1}=\hat{a}_{j}+\hat{a}_{j}^{\dagger}$, $\hat{c}_{j}^{2}=(-i)(\hat{a}_{j}-\hat{a}_{j}^{\dagger})$. Collecting the $p,q$ components of $\Gamma_{ij}^{pq}$ in a matrix,

[TABLE]

For a particle number conserving system, one can relate the correlation matrix to the correlator $\langle\hat{a}_{j}^{\dagger}\hat{a}_{l}\rangle$ using $W=\frac{1}{\sqrt{2}}\left(\begin{smallmatrix}1&1\\ -i&i\end{smallmatrix}\right)$ as

[TABLE]

For a free-fermion thermal state at temperature $\beta^{-1}$ and with single-particle Hamiltonian $\mathsf{H}$ this takes the form

[TABLE]

Viewing $\mathsf{\Gamma}$ as a single-particle Hamiltonian, it is gapped around zero energy if and only if $\mathsf{H}$ is an insulator. For such gapped $\mathsf{\Gamma}$, the topological properties (determined by invariants whose concrete form depends on the presence or absence of time-reversal and sublattice symmetries as well as the dimensionality of the system) are set by the “ground state” (negative “energy”) eigenvectors of $\mathsf{\Gamma}$ rhotop2 . Recognising the right-hand side as the Bogoliubov-de-Gennes Hamiltonian corresponding to $\tanh(\beta\mathsf{H}/2)$, we note that, up to a particle-hole redundancy, these eigenvectors have the same topology as the groundstate eigenvectors of $\tanh(\beta\mathsf{H}/2)$. These, in turn, are the same as those of $\mathsf{H}$ due to $\tanh(x)$ being a monotonically increasing odd function of $x$. Thus, the topology of the groundstate eigenvectors of $\mathsf{H}$ also determines that of the thermal density matrix.

Appendix C Relation to Grassmann tensor networks

In this Appendix, we briefly review the relation between our expressions and contractions of the Grassmann tensor networks introduced in Ref. Gu:2010, . The structure of such network contractions, specialized for the averages we consider, is illustrated in Fig. 3: links $ij$ between sites $i$ and $j$ are associated with link tensors $G_{ij}$; sites are associated with site tensors $T_{j}^{(X_{j})}$ (the superscript indicates that the concrete form of these depends on the action $X_{j}$ of operator $X$ on site $j$). The Grassmann variables $\theta_{j_{i}}$, $\bar{\theta}_{j_{i}}$ of the network are associated with the legs $j_{i}$ of sites $j$ (with $j_{i}$ denoting the end of link $ij$ running into site $j$). The link tensors $G_{ij}$ are Grassmann even expressions of $\theta_{i_{j}},\bar{\theta}_{i_{j}},\theta_{j_{i}},\bar{\theta}_{j_{i}}$; the site tensors $T_{j}^{(X_{j})}$ are expressions involving Grassmann differentials $\{d\theta_{J},d\bar{\theta}_{J}|J\in j\}$ with $J\in j$ denoting leg $J$ of site $j$. The Grassmann parity of $T_{j}^{(X_{j})}$ equals the fermion parity of $X_{j}$.

The tensor network contractions are of the form

[TABLE]

where $P_{0}$ denotes a projection on the scalar part of the resultant integral and where we separated a normalization factor as in Eq. (1). As described in Ref. Gu:2010, , such contractions may be efficiently approximated using the Grassmann tensor entanglement renormalization group; for works applying this scheme to concrete problems see e.g., Refs. Gu:2013a, , Gu:2013b, .

The conversion of the expression Eq. (1) or Eq. (5) is achieved using a combination of the Grassmann Leibniz rule Zinn-Justin:2002 and inserting suitable scalar part projections; this transforms the initial integral from one over Grassmann variables on the sites of the lattice into one over an increased number of Grassmann variables on the legs, as in Fig. 3. This leads to link tensors [in terms of Eq. (5), for concreteness]

[TABLE]

The site tensors where $X$ does not act ($X_{j}=1$) are

[TABLE]

For $X_{j}=\hat{a}_{j}$, $X_{j}=\hat{a}_{j}^{\dagger}$ , and $X_{j}=\hat{a}_{j}^{\dagger}\hat{a}_{j}$ we have

[TABLE]

and

[TABLE]

respectively. We emphasize, however, that for our purposes the conversion to the structure in Fig. 3 is not required; we directly work with the integrals of the form Eq. (5) and the resultant matrix expressions.

Appendix D Proof of Theorem 1

In this Appendix we briefly outline the proof of Theorem 1. This is based on the proof of Theorem 8.12 of Benzi et al. Benzi:2012 and its matrix generalizations discussed in the same reference.

A preliminary step is to show that a function $\tilde{f}$ that takes real values on $[\eta_{1}^{2},\eta_{2}^{2}]$, and is analytic in a neighbourhood of $[\eta_{1}^{2},\eta_{2}^{2}]$ in the complex plane, can be approximated by polynomials to exponential accuracy. Based on this, one can show that $f$ in Theorem 1 can also be approximated by polynomials to exponential accuracy. Then, the matrix decay bound in Theorem 1 follows via the spectral decomposition of $\mathsf{M}$. Our proof only differs from that of Benzi et al. Benzi:2012 in the preliminary step, and only in the choice of the ellipse $\mathcal{E}_{\chi}$ below, which accommodates the slightly weaker analyticity requirements on $\tilde{f}$ (and thus on $f$) compared to those in Ref. Benzi:2012, .

To make our formulation more precise, for integer $k\geq 0$, we define $E_{k}(\varphi,I)$ as the error for the $k$-th best approximation of a function $\varphi$ on a finite set $I$ of finite intervals:

[TABLE]

where $P_{k}$ is the set of all polynomials with real coefficients of degree no greater than $k$.

Consider now an ellipse $\mathcal{E}_{\chi}$ with foci at $\eta_{1}^{2}$ and $\eta_{2}^{2}$, and denote its semiaxes by $a$ and $b$ ($a\!>\!b$). With the constraint $a^{2}\!-\!b^{2}\!=\!c^{2}$ where the fixed semifocal length $c\!=\!\frac{\eta_{2}^{2}\!-\!\eta_{1}^{2}}{2}$, the ellipse is parametrised by $\chi\!=\!\frac{a\!+\!b}{c}$, which determines $a\!=\!\frac{\eta_{2}^{2}-\eta_{1}^{2}}{4}(\chi\!+\!\chi^{-1})$ and $b\!=\!\frac{\eta_{2}^{2}\!-\!\eta_{1}^{2}}{4}(\chi\!-\!\chi^{-1})$. The infimum of $\chi$ is $1$, when $\mathcal{E}_{\chi}$ reduces to the segment $[\eta_{1}^{2},\eta_{2}^{2}]$ on the real axis. The supremum $\bar{\chi}$ is set by the analytic properties of $\tilde{f}$; it corresponds to $\mathcal{E}_{\chi}$ first touching a singularity of $\tilde{f}$ as it is gradually blown up with increasing $\chi$. By the coordinate transformation

[TABLE]

we map $u\in[-1,1]$ to $x\in[\eta_{1}^{2},\eta_{2}^{2}]$ and the ellipse $\mathcal{E}^{0}_{\chi}$ with vertices $\pm(\chi\!+\!\chi^{-1})/2$ and $\pm\mathrm{i}(\chi\!-\!\chi^{-1})/2$ to $\mathcal{E}_{\chi}$. For $1<\chi<\bar{\chi}$, the function $\tilde{f}[x(u)]$ is analytic on $\mathcal{E}^{0}_{\chi}$ and takes real values on $[-1,1]$. One can thus invoke Bernstein’s theorem Bernstein (see also Theorem 8.7 in Ref. Benzi:2012, ):

[TABLE]

where $M(\chi)=\max_{z\in\mathcal{E}_{\chi}}|f(z)|$. Because $x(u)$ is a linear transformation, we also have

[TABLE]

Now taking $\chi$ such that $\pm\sqrt{\mathcal{E}_{\chi}}$ are in the interiors of the neighborhoods in Theorem 1, the proof regarding the polynomial approximation of $f$ on $D$ is identical to that of Theorem 8.12 in Ref. Benzi:2012, : this uses suitable combinations of $f_{1(2)}$ and $\sqrt{\cdot}$ for $\tilde{f}$ to show that there exists $K>0$ such that

[TABLE]

The final step is to use Eq. (32) to prove Theorem 1. As the strictly short-range Hermitian matrix $\mathsf{M}$ has range $R^{*}$, the matrix power $\mathsf{M}^{k}$ has range $kR^{*}$, and so does the polynomial function $p_{k}(\mathsf{M})$ with $p_{k}(x)$ the $k$-th best approximation of $f(x)$ on $D$. Now consider a pair of sites $i\neq j$, and take $k$ to be the largest integer such that $(k+1)R^{*}\geq R_{ij}>kR^{*}$. The matrix elements $[p_{k}(\mathsf{M})]_{ij}$ are zero. Therefore,

[TABLE]

where we use the matrix spectral norm, $||\mathsf{A}||_{2}=\max_{||x||=1}\sqrt{\langle x|\mathsf{A}^{\dagger}\mathsf{A}|x\rangle}$ or, equivalently, the maximum singular value of $\mathsf{A}$. (The inequality $|A_{ij}|\leq||\mathsf{A}||_{2}$ follows, e.g., from the Cauchy-Schwartz inequality applied to $A_{ij}=\langle e_{i}|\mathsf{A}e_{j}\rangle$ with basis vectors $|e_{i,j}\rangle.$) As $\mathsf{M}$ is Hermitian, we can diagonalise it using a unitary matrix $\mathsf{U}$, i.e., $\mathsf{M}=\mathsf{U}^{\dagger}\Lambda\mathsf{U}$ where $\mathsf{\Lambda}$ is a diagonal matrix containing the eigenvalues $\lambda_{j}$. We also have $f(\mathsf{M})=\mathsf{U}^{\dagger}f(\Lambda)\mathsf{U}$ and $p_{k}(\mathsf{M})=\mathsf{U}^{\dagger}p_{k}(\Lambda)\mathsf{U}$. As $f(\mathsf{M})-p_{k}(\mathsf{M})$ is also Hermitian, its singular values are the absolute values of its eigenvalues, $|f(\lambda_{j})-p_{k}(\lambda_{j})|$. Furthermore,

[TABLE]

where the first inequality holds because $\lambda_{j}\in D$, and the second because $p_{k}$ is the $k$-th best polynomial approximation. Defining $\tilde{K}=K\chi^{1/2}$ and $\tilde{\chi}=\chi^{1/2}$ we find

[TABLE]

which, upon renaming variables, proves Theorem 1.