Density-of-states of many-body quantum systems from tensor networks

TL;DR

This paper introduces a tensor network-based method to compute the density of states and thermodynamic properties of many-body quantum systems, validated on 1D models with promising extension potential.

Contribution

The authors develop a novel tensor network approach combining Hubbard-Stratonovich transformation and real-time evolution to evaluate spectral and thermodynamic properties.

Findings

Accurately reproduces thermodynamic quantities for 1D models

Method shows good agreement with exact results

Potential for extension to higher-dimensional systems

Abstract

We present a technique to compute the microcanonical thermodynamical properties of a manybody quantum system using tensor networks. The Density Of States (DOS), and more general spectral properties, are evaluated by means of a Hubbard-Stratonovich transformation performed on top of a real-time evolution, which is carried out via numerical methods based on tensor networks. As a consequence, the free energy and thermal averages can be also calculated. We test this approach on the one-dimensional Ising and Fermi-Hubbard models. Using matrix product states, we show that the thermodynamical quantities as a function of temperature are in very good agreement with the exact results. This approach can be extended to higher-dimensional system by properly employing other types of tensor networks.

Figures

Equations

$$\sqrt{\frac{\eta }{\pi }}{e}^{-\eta (H-E{)}^2}=\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }{e}^{-\frac{t^2}{4\eta }+itE}{e}^{-itH}dt,$$

$$D_{\eta }(E)=\frac{1}{\pi }{\int }_0^{+\infty }{e}^{-\frac{t^2}{4\eta }}\mathfrak{Re}\left(e^{itE}\text{mbox}Tr\left[e^{-itH}\right]\right)dt,$$

$${\bar{\Theta }}_{\eta }(E)=\frac{1}{\pi D_{\eta }(E)}{\int }_0^{+\infty }{e}^{-\frac{t^2}{4\eta }}\mathfrak{Re}\left(e^{itE}\text{mbox}Tr\left[\Theta e^{-itH}\right]\right)dt,$$

$$\mathrm{dt}$$

$$\bar{t}$$

$$\text{Complexity of HS-TN}=\mathcal{O}\left({\eta }^{\frac{1}{2}}{L}^2{m}^3\right).$$

Full text

Understanding the properties of strongly correlated quantum many-body systems is one of the main challenges of modern physics, with implications ranging from condensed matter to quantum technologies and high-energy physics. To this aim several different powerful numerical techniques such as Monte Carlo or different variants of numerical renormalisation and variational approaches (see e.g. Anderson (2007); Bulla et al. (2008); Schollwöck (2005)) have been developed in order to explore regimes unaccessible to analytical approaches. Each of those methods has advantages and shortcomings, all together they offer a solid platform for the investigation of strongly correlated phenomena. Despite the tremendous advancements experienced in the last decades by both theoretical and numerical approaches Wang (1994), the complete access to the spectral properties of a many-body quantum system is a long-sighted goal which remained elusive in those cases where accurate Monte Carlo sampling is not possible (due, for example, to sign problems Lyubartsev et al. (1992); Bennett (1976); Langfeld et al. (2012)). Having access to the many-body density of states (DOS) and energy-resolved observable quantities, thus effectively reconstructing its microcanonical ensemble, grants in turn full knowledge of its free-energy, and thus of the whole thermodynamics.

Here we introduce a technique, based on tensor network ansatz states, to compute the microcanonical density of states of a many-body quantum system, and from there to have a direct access to the free-energy of the system. During the last decades Tensor Networks (TNs) methods have proven to be one of the most promising approaches to low-dimensional quantum systems, both in and out of equilibrium Schollwoeck (2011); Orús (2014). The method we describe exploits two key features of TN techniques: the ability to encode efficiently a many-body quantum state and to simulate its (real-time) dynamics. We will show that a smoothened version of the microcanonical density of states can be evaluated using the Hubbard-Stratonovich transformation Stratonovich (1957); Hubbard (1959) and a subsequent time evolution of the system. We benchmark this approach against two paradigmatic models: the Ising model in transverse field and the Hubbard model. Here we focus on one-dimensional systems, however the method is in principle not bound to this limitation, provided an appropriate TN ansatz state is used.

TN methods are based on an efficient encoding of the many-body quantum state that allows simulations of large system sizes, unattainable via direct diagonalisation, while at the same time guaranteeing a remarkable precision. A vast class of TN algorithms have focused on an accurate determination of ground-state properties Schollwoeck (2011); Orús (2014). Prominent examples beyond the density matrix renormalisation group White (1992), the workhorse and inspiration of TN methods, are algorithms based on different one- and two-dimensional tensor network structures like Matrix Product States (MPS) Östlund and Rommer (1995); Verstraete et al. (2004a); Schollwoeck (2011), MERA Vidal (2007), PEPS Verstraete et al. (2006) or tree-TN Gerster et al. (2014). Their main common feature is the ability to minimise an energy functional over a tailored variational class of states and hence reconstruct, for instance, the many-body ground states of lattice Hamiltonians Verstraete et al. (2004a); Verstraete et al. (2008), and some elementary excitations Haegeman et al. (2013, 2012). TN techniques have been successfully extended to describe the evolution of isolated Vidal (2004); Daley et al. (2004); White and Feiguin (2004) and open quantum systems Zwolak and Vidal (2004); Verstraete et al. (2004b); Werner et al. (2016). The same techniques can be repurposed to reconstruct many-body states at finite temperature Verstraete et al. (2004b); White (2009); Molnar et al. (2015); Bañuls et al. (2015); Chen et al. (2017), thus proving that tensor networks are tools also capable of reconstructing efficiently canonical (Gibbs) ensembles Binder and Barthel (2015).

We introduce the method in Sec. I. Also, we put forward two proposals on how to adapt the standard tensor network methods for real-time evolution Vidal (2004); White and Feiguin (2004) to this purpose. In Sec. II, we test the approach on two paradigmatic models in one-dimension. Additionally, we show how the method can be extended to compute the microcanonical average of a generic observable. Finally, we draw the concluding remarks in Sec. III.

I Method

The calculation of the Density of States $\text{D}(E)=\sum_{\nu}\delta(E-E_{\nu})$ requires in principle the whole spectrum of the Hamiltonian $H$, as the sum $\sum_{\nu}$ runs over all the many-body eigenstates $|\nu\rangle$ of $H$ (with $H|\nu\rangle=|\nu\rangle E_{\nu}$). Here we show how to reconstruct, via TNs simulations, the $\text{D}(E)$ without the need to explicitly diagonalise $H$, in the same spirit of Ref. Osborne (2006). As a side note, once $\text{D}(E)$ is known, the free energy ${\cal F}$ can be easily extracted via ${\cal F}(T)=-T\ln\int e^{-E/T}\text{D}(E)dE$, where $T$ is the temperature (we set $k_{B}=\hbar=1$).

We consider a Gaussian broadening of the Dirac delta $\delta(x)\to\delta_{\eta}(x):=\sqrt{\frac{\eta}{\pi}}e^{-\eta x^{2}}$, which satisfies $\lim_{\eta\to\infty}\delta_{\eta}(x)=\delta(x)$ in the sense of distributions. Accordingly, we define a broadened density of states $D_{\eta}(E)=\sqrt{\frac{\eta}{\pi}}\sum_{\nu}e^{-\eta(E-E_{\nu})^{2}}=\sqrt{\frac{\eta}{\pi}}\mbox{Tr}[e^{-\eta(H-E)^{2}}]$. The corresponding `canonical' quantity ${\cal F}_{\eta}(T)=-T\ln\int e^{-E/T}D_{\eta}(E)dE$ is a good approximation of the free energy of the system ${\cal F}_{\eta}(T)\approx{\cal F}(T)$ as long as the broadening is sharper than the exponential decay introduced by temperature, i.e. $T\gg\eta^{-1/2}$. An important check in this respect is to guarantee that the $D_{\eta}(E)$ preserves the spectral measure, i.e. that $\int D_{\eta}(E)dE=\int\text{D}(E)dE$, which is, in turn, equal to the total number of states ($d^{L}$ where $d$ is the local space dimension and $L$ is the system size).

By means of Fourier analysis, it is possible to rewrite

[TABLE]

usually known as Hubbard-Stratonovich transformation Stratonovich (1957); Hubbard (1959). This effectively recasts the problem of characterising the DOS into a problem of simulating the real-time evolution of the quantum system. Indeed, the $D_{\eta}$ can be obtained via the numerical integral

[TABLE]

which represents the core of the Hubbard-Stratonovich Tensor Network (HS-TN) method introduced here. Notice that the Gaussian modulation $e^{-\frac{t^{2}}{4\eta}}$ drastically reduces the relevance of the long-time contributions to the $D_{\eta}$. As a consequence, it is sufficient to compute $\mbox{Tr}\!\left[e^{-itH}\right]$ for reasonably short times $t$, typically within a few orders of magnitude of $\sqrt{\eta}$ (fine-tuning details of the time-integration are discussed in Appendix A. This enables in practice the HS approach with tensor networks, since short-time evolutions are generally efficient and accurate in TN simulations.

In this paper we consider one-dimensional quantum systems with open boundary conditions and compute the time evolution using MPS Östlund and Rommer (1995); Verstraete et al. (2004a); Schollwoeck (2011) and the Time-Evolving Block Decimation (TEBD) algorithm Vidal (2004); White and Feiguin (2004). Alternative techniques to perform real-time evolution of TN, such as the Time-Dependent Variational Principle are similarly viable for this task Haegeman et al. (2011, 2016). In order to compute the trace $\mbox{Tr}\!\left[e^{-itH}\right]$ we propose two alternatives, or pathways:

(1) Random MPS sampling $-$ In this approach, the identity operator in $\mbox{Tr}\!\left[e^{-itH}{\mathbb{1}}\right]$ is approximated as the mixture of a large number $N_{R}$ of random states (see Appendix B for how to generate them). $\mathbb{1}\simeq\frac{d^{L}}{N_{R}}\sum_{R}^{N_{R}}|\psi_{R}\rangle\langle\psi_{R}|$ where the $|\psi_{R}\rangle$ are $N_{R}$ random states. This results in $\mbox{Tr}\!\left[e^{itH}\right]=\frac{d^{L}}{N_{R}}\sum_{R}^{N_{R}}\langle\psi_{R}(0)|\psi_{R}(t)\rangle$, where $|\psi_{R}(t)\rangle=e^{-itH}|\psi_{R}(0)\rangle=e^{-itH}|\psi_{R}\rangle$. The number of random states $N_{R}$ used to compute the average is increased until the convergence of $D_{\eta}(E)$ is observed. This method could be in principle partially biased as the random sampling is limited only to states which contain an amount of entanglement upper-bound by $\log m$, where $m$ is the MPS auxiliary bond dimension. However, this is not a fundamental limitation since it is conceptually possible to build many-body bases composed by states with low entanglement content only (e.g. canonical product bases).

(2) Local space doubling $-$ We briefly put forward an alternative route to the HS-TN calculation. In order to calculate the trace of a many-body operator $\mbox{Tr}[A]$, one can double the number of sites and consider the relation $\langle\Phi^{+}|^{\otimes L}(A_{\text{odd}}\otimes{\mathbb{1}}_{\text{even}})|\Phi^{+}\rangle^{\otimes L}=\mbox{Tr}[A]$, where $A$ acts only the odd sites, and $|\Phi^{+}\rangle^{\otimes L}=\bigotimes_{j=1}^{L}(\sum_{s=1}^{d}|s\rangle_{2j-1}|s\rangle_{2j})$. The (non-normalised) state $|\Phi^{+}\rangle^{\otimes L}$ can be easily expressed as an MPS with bond dimension $1$ and physical dimension $d^{2}$ (after merging together sites $2j-1$ and $2j$). This approach does not require a statistical average, however the drawback sits in the increase of the local space dimension to $d^{2}$. We realize that this route is more similar in spirit to Refs. Osborne (2006); Chen et al. (2017), but with the advantage that it uses the same pure-state time-evolution methods as pathway (1).

II Results

Hereafter we benchmark the presented method applying it to the spin-$\frac{1}{2}$ Ising model in a transverse field and the spin-$\frac{1}{2}$ Hubbard model, following pathway (1). We address Abelian global symmetries explicitly Singh et al. (2010a, b), and this allows us to reconstruct sector-specific broadened DOS, where now $D_{\eta,q}(E)=\sqrt{\frac{\eta}{\pi}}\sum_{\nu(q)}e^{-\eta(E-E_{\nu(q)})^{2}}$ for an arbitrary quantum number $q$, and the sum $\sum_{\nu(q)}$ runs over all the eigenstates $|\nu(q)\rangle$ of $H$ with a fixed quantum number $q$. Embedding such symmetry content is equivalent to substituting $\mbox{Tr}[e^{-itH}]\to\mbox{Tr}[e^{-itH}P_{q}]$ in Eq. (2), where $P_{q}=\sum_{\nu(q)}|\nu(q)\rangle\langle\nu(q)|$ is the projector over all the $|\nu(q)\rangle$ states. Notice that $[e^{-itH},P_{q}]=0$ since $q$ are good quantum numbers, connected to a symmetry of $H$. Calculating such trace is relatively straightforward in pathway (1), where one needs simply to sample random states $|\psi^{q}_{R}\rangle$ having quantum number $q$. The sector-specific densities of states can be eventually summed up to recover the full DOS, i.e. $D_{\eta}(E)=\sum_{q}D_{\eta,q}(E)$.

The spectrum of the Ising model, $H=-\sum_{j}^{L-1}\sigma^{x}_{j}\sigma^{x}_{j+1}+h\sum_{j}^{L}\sigma^{z}_{j}$, can be computed exactly since it is a free fermion theory Lieb et al. (1961) ($\sigma^{\alpha}$ are the Pauli matrices, $j$ labels the lattice sites, $h$ is the transverse field, and $L$ the number of sites). As a figure of merit for our technique, we consider the $D_{\eta,+}(E)$, restricted to the even `$+$' sector of the parity symmetry $\Pi=\sigma^{z\otimes L}$. We compare the $D_{\eta,+}$ calculated from the exact spectrum with the one obtained via the HS-TN method. Such comparison is shown in Fig. 1 (using $h=1,\eta=10L^{2}$) for different lattice sizes. Our analysis shows clearly that the HS-TN data converge to the exact solution. In Fig. 2 (main panel) we further compare the free energy per site ${\cal F}_{\eta}(T)/L$ calculated with the HS-TN method to the exact solution. We also show that we recover correctly the large $T$ limit behaviour, which predicts ${\cal F}(T)\approx-T(L-1)\ln 2+2^{1-L}\mbox{Tr}[H]$ (where $2^{L-1}$ is the dimension of the whole even symmetry sector). We observe numerical convergence with increasing bond dimension $m$, which we report in the top inset of Fig. 2 by plotting the deviation in $\mathcal{L}_{1}$ norm. Convergence with respect to the ensemble average $N_{R}$ is reported in Appendix B.

To demonstrate that the HS-TN method is not limited to the study of free-particle models, we apply it to the study of the spin-$\frac{1}{2}$ Fermi-Hubbard model, $H=U\sum_{j}^{L}n_{\uparrow,j}n_{\downarrow,j}-J\sum_{j}^{L-1}\sum_{\sigma}c^{\dagger}_{\sigma,j}c_{\sigma,j+1}+\text{h.c.}$, restricted to the symmetry sector of $N_{\uparrow}=L/2$ spin-up particles and $N_{\downarrow}=L/2$ spin-down particles (addressing explicitly the U(1) $\times$ U(1) subgroup of the full symmetry group for the model). The $D^{\prime}_{\eta}(E):=D_{\eta,N_{\uparrow}=N_{\downarrow}=L/2}(E)$ are plotted in Fig. 3 for different system sizes: $L=6$ which can be exactly diagonalized, and $L=50$ where convergence is observed. In the lower panels of Fig. 3, we report the entropy per site $S/L=-L^{-1}\partial{\cal F}/\partial T$, again confirming convergence to the exact solution where available. As a side note, we remark that when we adopt a broadening $\eta$ able to resolve the energy gap (such as in Fig. 3 for $L=6$), we can in turn access temperatures sufficiently low to display the insulating behaviour (signalled by the plateau in the bottom-left panel).

Spectral quantities $-$ The HS-TN method is not limited to the calculation of the DOS, but it can be applied to compute the spectral properties related to arbitrary observables $\Theta$. It is possible to calculate the expectation value $\bar{\Theta}_{\eta}(E)=\mbox{Tr}[\Theta\hat{\rho}_{\eta}(E)]$ over the broadened microcanonical ensemble $\hat{\rho}_{\eta}(E):=D_{\eta}^{-1}(E)\sqrt{\frac{\eta}{\pi}}e^{-\eta(E-H)^{2}}$. This allows to calculate the thermodynamical properties $\bar{\Theta}(T)=e^{{\cal F}_{\eta}/T}\int e^{-E/T}\bar{\Theta}_{\eta}(E)D_{\eta}(E)dE$ for temperatures $T\gg\eta^{-1/2}$. Analogously to Eq. (2), the spectral quantity can be recast into the integral

[TABLE]

via the HS equation. In this picture, numerically calculating $\mbox{Tr}[\Theta e^{-itH}]$ using TEBD requires little additional effort than calculating simply $\mbox{Tr}[e^{-itH}]$, provided that $\Theta$ can be efficiently expressed as a Matrix Product Operator Pirvu et al. (2010). In practice, this requires the computation of the delayed matrix elements $\langle\psi_{R}(0)|\Theta|\psi_{R}(t)\rangle$ for pathway (1), since $\mbox{Tr}\!\left[\Theta e^{itH}\right]=\frac{d^{L}}{N_{R}}\sum_{R}^{N_{R}}\langle\psi_{R}(0)|\Theta|\psi_{R}(t)\rangle$, or alternatively of $\langle\Phi^{+}|^{\otimes L}(\Theta e^{-itH})_{\text{odd}}\otimes{\mathbb{1}}_{\text{even}}|\Phi^{+}\rangle^{\otimes L}=\mbox{Tr}[\Theta e^{-itH}]$ for pathway (2).

Similarly, it is possible to reconstruct spectral probability distributions $p_{\eta}(E)=D_{\eta}(E)\left\langle\Psi|\hat{\rho}_{\eta}(E)|\Psi\right\rangle$ for a given MPS $|\Psi\rangle$. Figure 4 shows an example of such treatment, performed on the even sector of the Ising model, on which we spectrally probe the nearest-neighbour correlations of the paramagnetization $\langle\sigma^{z}_{j}\sigma^{z}_{j+1}\rangle$.

III Conclusions and Outlook

We introduced the HS-TN method to reconstruct the complete broadened microcanonical ensemble of one-dimensional quantum many-body models. It relies on performing real short-time evolution of Matrix Product States, which can be performed via standard TEBD algorithm, and then manipulates the information of the dynamics into spectral quantities by the Hubbard-Stratonovich transformation. The method we discussed has the advantage of accessing system sizes far beyond the exact diagonalization, while at the same time capable of simulating systems at finite densities or fermionic systems without incurring into any sign problem. Finally, we remark that techniques for performing time-evolution of MPS (such as TEBD or TDVP) are well established and largely available as free-released software, ultimately making our method a ready-to-use toolkit. We benchmarked the method by explicitly calculating the broadened density of states and some corresponding thermodynamical quantities for free and interacting fermion models.

We conjecture that the good accuracy of the HS-TN method even with relatively small bond dimensions $m$ is due to mixing effects. Indeed, on the one hand, high-energy eigenstates might be highly-entangled. On the other hand, it is often possible to express mixtures (in this case Gaussian ensembles) of states (at fixed resolution $\eta$, in numbers exponentially increasing with $L$) with a relatively small amount of quantum correlations, and thus to represent such mixtures efficiently with a tensor network. In fact, this effect has been already observed and discussed in other scenarios Verstraete et al. (2004b); Werner et al. (2016); las Cuevas et al. (2013); White (2009); Binder and Barthel (2015). More quantitative observations regarding such behaviour are reported in Appendix C.

We showed explicit examples of the HS-TN technique applied to one-dimensional models in open boundary conditions. However, we stress that this method extends to all those models whose out-of-equilibrium, short timescale dynamics can be numerically simulated efficiently and precisely, thus effectively extending the scope of our treatment beyond 1D.

This manuscript opens wide new perspectives in the framework of tensor networks methods, as we have shown that they are suitable tools to study complex quantum systems well beyond the low-temperatures limit. The HS-TN technique can be applied to a plethora of cases, e.g. the evaluation of thermodynamical and transport properties of quantum systems (spin chains, fermionic and bosonic models, nanotubes, small molecules, etc.) whose detailed knowledge has been so far precluded.

Acknowledgements $-$ The authors acknowledge support from the EU via RYSQ, UQUAM and QUIC projects, the DFG via SFB/TRR21 project, the Baden-Württemberg Stiftung via Eliteprogramm for Postdocs and the Carl Zeiss Foundation. SM gratefully acknowledges the support of the DFG via a Heisenberg fellowship. RF kindly acknowledges support from the National Research Foundation of Singapore (CRP - QSYNC) and the Oxford Martin School.

Appendix A Parameters and computational complexity

We now discuss in detail fine-tuning issues of the simulation parameters and the related computational complexity. Hereafter, we consider a MPS ansatz of bond dimension $m$ for the real-time dynamics, on a 1D system of size $L$ with open boundary conditions. The parameter $\eta$ introduces a spectral broadening of magnitude $\delta E\sim\eta^{-1/2}$ We numerically evaluate the HS-integral given in Eq. (2) by introducing a discrete time-step $\mathrm{dt}$ and finite evolution-time bound $\bar{t}$ reducing the integration interval $\int_{0}^{\infty}$ to $\int_{0}^{\bar{t}}$. We tune the integration parameters according to

[TABLE]

where $\Delta E=E_{\mathrm{max}}-E_{\mathrm{min}}$ is the estimated width of the complete spectrum, and $\mathrm{erfc}(z)=\frac{2}{\sqrt{\pi}}\int_{z}^{\infty}e^{-t^{2}}dt$ is the complementary error function. The additional fine-tuning parameter $\chi\ll 1$ controls the truncation-error of the Gaussian modulation (we typically set $\chi$ of the order of $10^{-10}$).

As the models we typically consider are short-range, we can safely assume a linear bound with $L$ to the spectral width $\Delta E\propto L$. Under this condition, we observed that we can employ dt as time-step for the TEBD algorithm in order to obtain a target precision independent of $L$. Consequently, the computational complexity of the algorithm equals the cost of performing ${\cal N}$ time-steps of TEBD, where ${\cal N}=\bar{t}/\mathrm{dt}\propto\Delta E/\delta E\propto\sqrt{\eta}L$. As a single TEBD time-step carries a cost of the order $\mathcal{O}\left(Lm^{3}\right)$ [21$-$23] (including measurements), the overall cost of the HS-TN algorithm with MPS is

[TABLE]

The computation of the $D_{\eta}(E)$ using $\mathcal{O}(\mathcal{N})$ points $E\in\left[E_{\mathrm{min}},E_{\mathrm{max}}\right]$ is simply a (fast) Fourier-transformation, of marginal computational cost.

Two noteworthy remarks follow: First, the bond dimension $m$ is expected to increase with $L$, however our test-cases shows that a quite small $m$ over various system-sizes is sufficient to obtain accurate results. Second, in pathway (1), the random sampling introduces an additional factor $N_{R}$ to the overall cost. However, we observed that often we do not need to increase $N_{R}$ for larger $L$ to maintain the fixed precision (on the contrary, the precision can even improve with $L$ at fixed $N_{R}$).

Appendix B Convergence with $N_{R}$ and generation of symmetric random states

Here we report the convergence of the ensemble reconstruction using pathway (1) (random initial state sampling), as a function of the number $N_{R}$ of initial random states employed. To properly test convergence with $N_{R}$, we consider results at convergence in the other refinement parameters (e.g. $m$). Figures 5 and 6 show convergence of the free energy calculated via HS-TN method, as a function of $N_{R}$, for the Ising model (even symmetry sector) and the Hubbard model (half-half filling symmetry sector), respectively. In the Ising model, we observe once more a power-law convergence to the exact result of the error $\epsilon^{\prime}$ as a function of $1/N_{R}$.

In all cases, the initial random MPS with bond dimension $m$, and quantum number $q$ of a global Abelian symmetry, were sampled as follows: We randomly draw quantum numbers $q$, for each bond-index $1,\dots,m$ and each bond, from the maximal set of all nontrivial contributions. This is equivalent to picking a random $m$-dimensional sub-representation of the largest possible symmetry representation for each bond link. Then, as an MPS in a fixed global sector $q$ has all block-diagonal tensors, we sample the elements of the nontrivial blocks from a normal distribution. Such sampled states were then normalized and gauged properly for the TEBD evolution.

Appendix C Resolution and scaling with $L$

In this section, we now focus on the interplay between the system-sizes $L$ and the broadening-parameter $\eta$.

As a measure of precision, we consider the error $\epsilon_{0}(m,\eta)$, defined as in Fig. (1), between the normalized broadened $D_{\eta,+,\text{exact}}$ (exact result) and $D_{\eta,+,m}$ computed with the HS-TN method, over various bond-dimensions $m$ and broadening-parameters $\eta$ in the Ising model. Our results are summarized in Fig. 7, and lead to the following observations: It emerges that for a given system size $L$, there is an inherent broadening threshold $\bar{\eta}(L)$, or equivalently a spectral resolution threshold $\overline{\delta E}=\bar{\eta}^{-1/2}$ accessible via HS-TN. The top panel of Fig. 7 shows that the error $\epsilon_{0}(m,\eta)$ converges to zero rapidly in $m$ as long as $\delta E>\overline{\delta E}$ (i.e. $\eta<\bar{\eta}(L)$). Conversely, when a resolution higher than this threshold is selected ($\delta E<\overline{\delta E}$) a non-zero plateau in the error emerges for large $m$, until we reach entanglement saturation of the TN ansatz. This behaviour is highlighted in the inset (red data points). In turn, in all the cases we considered, such resolution threshold $\overline{\delta E}(L)$ improves with the system size, i.e. it decreases exponentially with $L$ (bottom panel).

We interpret this behaviour with an argument based on the energy spacings of the Hamiltonian $H$, i.e. the differences between subsequent eigenvalues of $H$: It appears, in fact, that the spectral resolution threshold is well approximated by the typical energy spacing $\overline{\delta E}\approx\Delta E/d^{L}$ (or more precisely $\overline{\delta E}\approx\sqrt{\left\langle dE_{\nu}^{2}\right\rangle}$, compare red and purple data sets in the bottom plot). If the spectral broadening we introduce is larger than the typical energy-spacing, then we do not need to reconstruct single eigenstates of $H$ (expensive to represent with tensor networks when they exhibit a volume-law of entanglement) but mixtures of them, which instead may be efficiently approximated by TN. According to this argument, the HS-TN method will converge rapidly in $m$ if we set $\delta E>\overline{\delta E}$, as observed above.

Finally, since the total spectral width is usually extensive $\Delta E\propto L$, the average energy gap decays roughly exponentially with $L$ ($\overline{\delta E}\propto e^{-L}$). From this follows that, as the thermodynamical limit is approached $L\to\infty$, the resolution threshold vanishes $\overline{\delta E}\to 0$, and we expect the HS-TN method to converge regardless of the chosen broadening.