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

###### Abstract

We present a technique to compute the microcanonical thermodynamical properties of a many-body 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.

###### pacs:

05.30.-d, 02.70.-c, 05.10.-a.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 requires in principle the whole spectrum of the Hamiltonian , as the sum runs over all the many-body eigenstates of (with ). Here we show how to reconstruct, via TNs simulations, the without the need to explicitly diagonalise , in the same spirit of Ref. Osborne (2006). As a side note, once is known, the free energy can be easily extracted via , where is the temperature (we set ).

We consider a Gaussian broadening of the Dirac delta , which satisfies in the sense of distributions. Accordingly, we define a broadened density of states . The corresponding ‘canonical’ quantity is a good approximation of the free energy of the system as long as the broadening is sharper than the exponential decay introduced by temperature, i.e. . An important check in this respect is to guarantee that the preserves the spectral measure, i.e. that , which is, in turn, equal to the total number of states ( where is the local space dimension and is the system size).

By means of Fourier analysis, it is possible to rewrite

(1) |

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 can be obtained via the numerical integral

(2) |

which represents the core of the Hubbard-Stratonovich Tensor Network (HS-TN) method introduced here. Notice that the Gaussian modulation drastically reduces the relevance of the long-time contributions to the . As a consequence, it is sufficient to compute for reasonably short times , typically within a few orders of magnitude of (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 we propose two alternatives, or pathways:

(1) Random MPS sampling In this approach, the identity operator in is approximated as the mixture of a large number of random states (see Appendix B for how to generate them). where the are random states. This results in , where . The number of random states used to compute the average is increased until the convergence of 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 , where 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 , one can double the number of sites and consider the relation , where acts only the odd sites, and . The (non-normalised) state can be easily expressed as an MPS with bond dimension and physical dimension (after merging together sites and ). This approach does not require a statistical average, however the drawback sits in the increase of the local space dimension to . 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- Ising model in a transverse field and the spin- 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 for an arbitrary quantum number , and the sum runs over all the eigenstates of with a fixed quantum number . Embedding such symmetry content is equivalent to substituting in Eq. (2), where is the projector over all the states. Notice that since are good quantum numbers, connected to a symmetry of . Calculating such trace is relatively straightforward in pathway (1), where one needs simply to sample random states having quantum number . The sector-specific densities of states can be eventually summed up to recover the full DOS, i.e. .

The spectrum of the Ising model, , can be computed exactly since it is a free fermion theory Lieb et al. (1961) ( are the Pauli matrices, labels the lattice sites, is the transverse field, and the number of sites). As a figure of merit for our technique, we consider the , restricted to the even ‘’ sector of the parity symmetry . We compare the calculated from the exact spectrum with the one obtained via the HS-TN method. Such comparison is shown in Fig. 1 (using ) 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 calculated with the HS-TN method to the exact solution. We also show that we recover correctly the large limit behaviour, which predicts (where is the dimension of the whole even symmetry sector). We observe numerical convergence with increasing bond dimension , which we report in the top inset of Fig. 2 by plotting the deviation in norm. Convergence with respect to the ensemble average 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- Fermi-Hubbard model, , restricted to the symmetry sector of spin-up particles and spin-down particles (addressing explicitly the U(1) U(1) subgroup of the full symmetry group for the model). The are plotted in Fig. 3 for different system sizes: which can be exactly diagonalized, and where convergence is observed. In the lower panels of Fig. 3, we report the entropy per site , again confirming convergence to the exact solution where available. As a side note, we remark that when we adopt a broadening able to resolve the energy gap (such as in Fig. 3 for ), 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 . It is possible to calculate the expectation value over the broadened microcanonical ensemble . This allows to calculate the thermodynamical properties for temperatures . Analogously to Eq. (2), the spectral quantity can be recast into the integral

(3) |

via the HS equation. In this picture, numerically calculating using TEBD requires little additional effort than calculating simply , provided that can be efficiently expressed as a Matrix Product Operator Pirvu et al. (2010). In practice, this requires the computation of the delayed matrix elements for pathway (1), since , or alternatively of for pathway (2).

Similarly, it is possible to reconstruct spectral probability distributions for a given MPS . 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 .

## 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 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 , in numbers exponentially increasing with ) 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.

## References

- Anderson (2007) J. B. Anderson, Quantum Monte Carlo: origins, development, applications (Oxford University Press, 2007).
- Bulla et al. (2008) R. Bulla, T. A. Costi, and T. Pruschke, Rev. Mod. Phys. 80, 395 (2008).
- Schollwöck (2005) U. Schollwöck, Rev. Mod. Phys. 77, 259 (2005).
- Wang (1994) L.-W. Wang, Phys. Rev. B 49, 10154 (1994).
- Lyubartsev et al. (1992) A. P. Lyubartsev, A. A. Martsinovski, S. V. Shevkunov, and P. N. Vorontsov-Velyaminov, The Journal of Chemical Physics 96, 1776 (1992).
- Bennett (1976) C. H. Bennett, Journal of Computational Physics 22, 245 (1976), ISSN 0021-9991.
- Langfeld et al. (2012) K. Langfeld, B. Lucini, and A. Rago, Phys. Rev. Lett. 109, 1 (2012), ISSN 00319007, eprint 1204.3243.
- Schollwoeck (2011) U. Schollwoeck, Ann. Phys. 326, 96 (2011).
- Orús (2014) R. Orús, The European Physical Journal B 87, 1 (2014).
- Stratonovich (1957) R. Stratonovich, in Soviet Physics Doklady (1957), vol. 2, p. 416.
- Hubbard (1959) J. Hubbard, Phys. Rev. Lett. 3, 77 (1959).
- White (1992) S. R. White, Physical Review Letters 69, 2863 (1992).
- Östlund and Rommer (1995) S. Östlund and S. Rommer, Phys. Rev. Lett. 75, 3537 (1995), ISSN 0031-9007.
- Verstraete et al. (2004a) F. Verstraete, D. Porras, and J. I. Cirac, Phys. Rev. Lett. 93, 227205 (2004a).
- Vidal (2007) G. Vidal, Physical review letters 99, 220405 (2007).
- Verstraete et al. (2006) F. Verstraete, M. M. Wolf, D. Perez-Garcia, and J. I. Cirac, Physical review letters 96, 220601 (2006).
- Gerster et al. (2014) M. Gerster, P. Silvi, M. Rizzi, R. Fazio, T. Calarco, and S. Montangero, Physical Review B 90, 125154 (2014).
- Verstraete et al. (2008) F. Verstraete, V. Murg, and J. I. Cirac, Advances in Physics 57, 143 (2008).
- Haegeman et al. (2013) J. Haegeman, S. Michalakis, B. Nachtergaele, T. J. Osborne, N. Schuch, and F. Verstraete, Phys. Rev. Lett. 111, 080401 (2013).
- Haegeman et al. (2012) J. Haegeman, B. Pirvu, D. J. Weir, J. I. Cirac, T. J. Osborne, H. Verschelde, and F. Verstraete, Phys. Rev. B 85, 100408 (2012).
- Vidal (2004) G. Vidal, Phys. Rev. Lett. 93, 040502 (2004).
- Daley et al. (2004) A. J. Daley, C. Kollath, U. Schollwöck, and G. Vidal, Journal of Statistical Mechanics: Theory and Experiment 2004, P04005 (2004).
- White and Feiguin (2004) S. R. White and A. E. Feiguin, Phys. Rev. Lett. 93, 076401 (2004).
- Zwolak and Vidal (2004) M. Zwolak and G. Vidal, Phys. Rev. Lett. 93, 207205 (2004).
- Verstraete et al. (2004b) F. Verstraete, J. J. Garcia-Ripoll, and J. I. Cirac, Physical review letters 93, 207204 (2004b).
- Werner et al. (2016) A. Werner, D. Jaschke, P. Silvi, M. Kliesch, T. Calarco, J. Eisert, and S. Montangero, Physical Review Letters 116, 237201 (2016).
- White (2009) S. R. White, Phys. Rev. Lett. 102, 190601 (2009).
- Molnar et al. (2015) A. Molnar, N. Schuch, F. Verstraete, and J. I. Cirac, Phys. Rev. B 91, 045138 (2015).
- Bañuls et al. (2015) M. C. Bañuls, K. Cichy, J. I. Cirac, K. Jansen, and H. Saito, Phys. Rev. D 92, 034519 (2015).
- Chen et al. (2017) B.-B. Chen, Y.-J. Liu, Z. Chen, and W. Li, Phys. Rev. B 95, 161104 (2017).
- Binder and Barthel (2015) M. Binder and T. Barthel, Phys. Rev. B 92, 125119 (2015).
- Osborne (2006) T. J. Osborne, arXiv preprint cond-mat/0605194 (2006).
- Haegeman et al. (2011) J. Haegeman, J. I. Cirac, T. J. Osborne, I. Pižorn, H. Verschelde, and F. Verstraete, Phys. Rev. Lett. 107, 070601 (2011).
- Haegeman et al. (2016) J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken, and F. Verstraete, Phys. Rev. B 94, 165116 (2016).
- Singh et al. (2010a) S. Singh, R. N. C. Pfeifer, and G. Vidal, Phys. Rev. A 82, 050301 (2010a).
- Singh et al. (2010b) S. Singh, H.-Q. Zhou, and G. Vidal, New Journal of Physics 12, 033029 (2010b).
- Lieb et al. (1961) E. Lieb, T. Schultz, and D. Mattis, Annals of Physics 16, 407 (1961), ISSN 0003-4916.
- Pirvu et al. (2010) B. Pirvu, V. Murg, J. I. Cirac, and F. Verstraete, New Journal of Physics 12, 025012 (2010).
- las Cuevas et al. (2013) G. D. las Cuevas, N. Schuch, D. Pérez-García, and J. I. Cirac, New Journal of Physics 15, 123021 (2013).

## 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 for the real-time dynamics, on a 1D system of size with open boundary conditions. The parameter introduces a spectral broadening of magnitude We numerically evaluate the HS-integral given in Eq. (2) by introducing a discrete time-step and finite evolution-time bound reducing the integration interval to . We tune the integration parameters according to

(4) | ||||

where is the estimated width of the complete spectrum, and is the complementary error function. The additional fine-tuning parameter controls the truncation-error of the Gaussian modulation (we typically set of the order of ).

As the models we typically consider are short-range, we can safely assume a linear bound with to the spectral width . 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 . Consequently, the computational complexity of the algorithm equals the cost of performing time-steps of TEBD, where . As a single TEBD time-step carries a cost of the order [2123] (including measurements), the overall cost of the HS-TN algorithm with MPS is

(5) |

The computation of the using points is simply a (fast) Fourier-transformation, of marginal computational cost.

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

## Appendix B Convergence with 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 of initial random states employed. To properly test convergence with , we consider results at convergence in the other refinement parameters (e.g. ). Figures 5 and 6 show convergence of the free energy calculated via HS-TN method, as a function of , 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 as a function of .

In all cases, the initial random MPS with bond dimension , and quantum number of a global Abelian symmetry, were sampled as follows: We randomly draw quantum numbers , for each bond-index and each bond, from the maximal set of all nontrivial contributions. This is equivalent to picking a random -dimensional sub-representation of the largest possible symmetry representation for each bond link. Then, as an MPS in a fixed global sector 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

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

As a measure of precision, we consider the error , defined as in Fig. (1), between the normalized broadened (exact result) and computed with the HS-TN method, over various bond-dimensions and broadening-parameters 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 , there is an inherent broadening threshold , or equivalently a spectral resolution threshold accessible via HS-TN. The top panel of Fig. 7 shows that the error converges to zero rapidly in as long as (i.e. ). Conversely, when a resolution higher than this threshold is selected () a non-zero plateau in the error emerges for large , 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 improves with the system size, i.e. it decreases exponentially with (bottom panel).

We interpret this behaviour with an argument based on the energy spacings of the Hamiltonian , i.e. the differences between subsequent eigenvalues of : It appears, in fact, that the spectral resolution threshold is well approximated by the typical energy spacing (or more precisely , 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 (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 if we set , as observed above.

Finally, since the total spectral width is usually extensive , the average energy gap decays roughly exponentially with (). From this follows that, as the thermodynamical limit is approached , the resolution threshold vanishes , and we expect the HS-TN method to converge regardless of the chosen broadening.