Sebastian Paeckel
email: sebastian.paeckel@physik.uni-muenchen.de

Thomas Köhler
email: thomas.kohler@physics.uu.se

Salvatore R Manmana

Benjamin Lenz
email: benjamin.lenz@sorbonne-universite.fr

Matrix-product-state-based band-Lanczos solver for quantum cluster approaches

We present a matrix-product state (MPS) based band-Lanczos method as solver for quantum cluster methods such as the variational cluster approximation (VCA). While a naïve implementation of MPS as cluster solver would barely improve its range of applicability, we show that our approach makes it possible to treat cluster geometries well beyond the reach of exact diagonalization methods.

The key modifications we introduce are a continuous energy truncation combined with a convergence criterion that is more robust against approximation errors introduced by the MPS representation and provides a bound to deviations in the resulting Green's function. The potential of the resulting cluster solver is demonstrated by computing the self-energy functional for the single-band Hubbard model at half filling in the strongly correlated regime, on different cluster geometries. Here, we find that only when treating large cluster sizes, observables can be extrapolated to the thermodynamic limit, which we demonstrate at the example of the staggered magnetization. Treating clusters sizes with up to 6 × 6 sites we obtain excellent agreement with quantum Monte-Carlo results.

I. INTRODUCTION

Exploring the quantum states of matter of two-dimensional strongly correlated electron systems is inherently difficult. Prime examples are the cuprate hightemperature superconductors [1], whose salient features are frequently described within a single-band Hubbardtype model on a square lattice geometry [2][START_REF] Leblanc | on the Many-Electron Problem)[END_REF][4]. More recently, two-dimensional materials have attracted substantial interest in a variety of systems, like bilayer graphene and other multi-layered materials [5][6][7], transfer metal dichalcogenides (TMDCs) [8], or kagomé metals [9,10]. Experiments on these materials have revealed unconventional phases of matter, e.g., a Wigner-crystal quantum Hall phase [11], or interesting superconducting phases [5][6][7][12][13][14], for which the role of flat bands and correlation effects is investigated. Hence, there is an urgent need to further develop numerical methods to treat 2D systems.

Important progress in describing strongly correlated two-dimensional systems has been achieved by the development of methods, which help to investigate quantum phases and quantum critical behavior directly in the thermodynamic limit. The fundamental insight is that the challenge of solving systems in the thermodynamic limit can be simplified by replacing it with an equivalent problem that is defined on a finite cluster, only. This is the foundation of so-called cluster-embedding techniques (CETs) such as dynamical mean-field theory (DMFT) [15][16][17], cluster-perturbation theory (CPT) [18,19] or VCA [20][21][22][23]. These approaches are often formulated in terms of Green's functions of the finite cluster and the embedding is defined such that the cluster Green's function describes best that of the actual physical system. Despite others, the success of these methods is based on the intimate relation between Green's functions and experimental observables, for instance spectral functions, which allows for the direct comparison between theory and experiment [24].

With the development of CETs, which contain the mapping of lattice problems to finite clusters, the focus has shifted to determine the cluster Green's function, which still constitutes a computational highly non-trivial problem. Important classes of methods to approach this task are wave function-based methods, like exact diagonalization (ED) [25,[START_REF] Sandvik | AIP Conference Proceedings 1297[END_REF] or tensor-network state (TNS) methods [START_REF] Cirac | [END_REF]28], and sampling-based methods, like quantum Monte Carlo (QMC) [START_REF]Quantum Simulations of Complex Many-Body Systems: From Theory to Algorithms[END_REF][START_REF]Quantum Monte Carlo Methods in Physics and Chemistry[END_REF][START_REF] Albuquerque | International Conference on Magnetism, The International Conference on Magnetism[END_REF]. ED cluster solvers are conceptionally the most simple approach and despite the fact that an extreme degree of optimization has been achieved, incorporating various symmetries and exploiting massive parallelization [START_REF] Wietek | [END_REF], they suffer drastically from the exponential growth of the Hilbert space, limiting the practically doable number of cluster sites for single band Hubbard models to ∼ O (20). On the other hand, TNS realized as MPS [33] are extremely successful in one dimension, but the area law of entanglement or the entanglement barrier in non-equilibrium situations highly limits their applicability to two-dimensional systems [34]. From a quantum information perspective, finite projected entangled-pair states (PEPS) [28,35] should be optimal. However, these algorithms have been shown to exhibit only extremly slow convergence for Hubbard-type systems [36]. Finally, QMC methods have their limitations due to the fermionic sign problem [37], which restricts the class of solvable cluster Hamiltonians.

In this paper, we introduce an efficient MPS-based scheme for calculating the cluster Green's function and apply it within the context of VCA. It turns out that a naïve treatment of established Krylov-based methods in a MPS framework would barely improve the range of applicability of the VCA. However, when combining the MPS with a band Lanczos scheme as well as an energy truncation, it is possible to treat cluster geometries well outside the reach of ED solvers and without principal limitations on the Hamiltonian.

The remainder of the paper is structured as follows: In Sec. II we introduce the single-band Hubbard model and the basic idea of cluster methods, such as the VCA. In Sec. III we discuss in detail our band-Lanczos ansatz for MPS, including the energy truncation scheme, an error estimate using the Hochbruck-Lubich bound [38], a useful energy rescaling, and a discussion of the loss of orthogonality of the Lanczos vectors within the MPS framework. Furthermore, Sec. III B discusses further aspects of the VCA as relevant for this paper, and Sec. III C discusses our choice of clusters, including Betts clusters. Section IV presents the main results of this paper obtained with our band-Lanczos MPS+VCA-ansatz, including a detailed analysis of the scaling of the results with the MPS bond dimension and a comparison of the staggered magnetization after scaling to infinite cluster size with numerically exact QMC results, as well as results for the spectral function. In Sec. V we conclude and provide an outlook to further applications of our method. The appendix presents technical details of the calculations.

II. MODEL AND NUMERICAL TECHNIQUE

The model for which we benchmark our solver is one of the most common ones to study strongly correlated electron systems in two dimensions, the one-band Hubbard model with nearest-neighbor hopping [39]. The corresponding Hamiltonian reads

Ĥ = -t ⟨i,j⟩,σ ĉ † i,σ ĉj,σ + h.c. + U i ni,↑ ni,↓ , (1)
where ĉ( †) j,σ are the annihilation (creation) operators of electrons of spin σ and nj,σ = ĉj,σ ĉ † j,σ is the particle number operator for a given spin σ on site j.

For our band Lanczos solver we have quantum cluster techniques in mind, for which a finite-size scaling in cluster size is known to be notoriously difficult in 2D [START_REF] Leblanc | on the Many-Electron Problem)[END_REF]4]. In particular embedding the cluster into an effective environment, various factors need to be considered, such as the ratio between cluster-bulk and cluster-boundary sites, the connectivity of boundary sites or the cluster geometry. Given these complications, which we discuss in more detail in Sec. III C, the main goal of our band-Lanczos solver is to provide a tool that can handle var-ious cluster geometries while increasing the cluster sizes as much as possible. A prime example for a quantum cluster technique is CPT [18,19], which can be derived from strong-coupling perturbation theory [40,41], and which allows to calculate the spectral function based on the cluster Green's function. In CPT, the lattice is tiled into a super-lattice of -most of the time identical -clusters. The inter-cluster hopping terms are then collected into a matrix V and treated in lowest-order perturbation theory when constructing the CPT Green's function:

G CPT ( k, ω) = (G cluster (ω)) -1 -V( k) -1 . (2)
Here k are reduced wave-vectors corresponding to the partial Fourier transform with respect to the superlattice vectors. In the following, we will omit the superscript and refer to the cluster Green's function as G(ω).

The error of the perturbative treatment within CPT is expected to reduce when enlarging the ratio of cluster 'bulk' to cluster boundary, for instance when scaling up clusters of fixed geometry. In addition to the cluster size, different geometries allow to precisely account for the spectral function at different wave vectors k.

CPT has been used not only for model calculations on the Hubbard [19], t-J [42] and pure spin models [43], but also for more realistic modeling of materials, e.g., in the context of Mott insulators such as NiO [44] or (doped) cuprates [41,[45][START_REF] Dahnken | High Performance Computing in Science and Engineering[END_REF][START_REF] Sénéchal | [END_REF]. It can also be extended to non-local interactions [48], electron-phonon coupling [49], non-equilibrium problems [50] or even be used to calculate two-particle responses [51]. CPT is thereby not only conceptually simple, but also one of the most versatile quantum cluster techniques. However, CPT is limited in that it does not allow to study symmetry-broken phases. An extension of the technique mends this deficiency by including symmetry-breaking Weiss fields on the cluster, which are determined according to a variational principle. The variational cluster perturbation theory, more commonly known as VCA [20,21], will be introduced in more detail in section III B.

For benchmarking purposes, we add here a Weiss field term on the cluster corresponding to a staggered magnetic field,

ĤAF = h z l e iQ•R l (n R l,↑ -nR l,↓ ), (3)
where l denotes cluster sites and the Néel ordering wave vector is Q = (π, π). A detailed VCA study of antiferromagnetism (AF) in the 2D Hubbard model using this Weiss field can be found in Ref. 52. For the following, it is important to note that the (cluster) Green's function of Ĥ + ĤAF needs to be determined.

III. METHOD

In this section we introduce the band Lanczos algorithm for calculating the cluster Green's function using MPS. Several routes have been pursued so far to employ MPS-based cluster solvers for CET. This includes time-evolution based solvers obtaining the cluster Green's function from real or in imaginary time evolution [53][54][55][56][57], expanding the resolvent in terms of orthogonal polynomials [58][59][60], reformulating the determination of the resolvent as optimization problem using correction vectors [61,62], or evaluating the Green's matrix elements explicitly by means of a global Lanzcos approach [63].

Nevertheless, despite each of these methods being considerably successful in the past, their applicability to VCA is rather limited. Standard time-evolution based solvers require considerable system sizes to avoid too fast entanglement growth, due to reflections of excitations at the cluster boundary. Evaluating the resolvent in terms of orthogonal polynomials on the other hand, typically requires a large number ∼ O(100) of basis states, each of which is obtained by applying the Hamiltonian to its predecessor, and thereby rapidly increases the state's entanglement. This limitation is partially overcome by a direct Lanczos expansion of each Green's matrix element, which typically requires only ∼ O(10) applications of the Hamiltonian, per matrix element. However, a lot of information about the low-energy states is generated repeatedly, since each matrix element is constructed independently of every other. Furthermore, the reduced numerical precision using MPS arithmetics typically causes severe orthogonality problems after more than 10 applications of the Hamiltonian. The necessary application of further reorthogonalization techniques common to global Krylov methods [64] deteriorates the computational efficiency such that cluster sizes beyond those doable using exact methods can only hardly be reached.

In sight of these complications, the band Lanczos [START_REF] Freund | Templates for the Solution of Algebraic Eigenvalue Problems: A Practical Guide[END_REF][START_REF] Baker | Direct solution of multiple excitations in a matrix product state with block lanczos[END_REF] offers a promising approach to evaluate the cluster Green's function in a global Krylov-subspace representation, employing various initial states. In the following, we first summarize the general idea of the band Lanczos and introduce the necessary steps to obtain a meaningful convergence criterium. Afterwards, we recapitulate the necessary theory of the VCA, in order to use the band Lanczos algorithm as its cluster solver.

A. Band Lanczos with MPS

The quantity of interest of most CETs is the cluster Green's function in frequency space, which can be written componentwise for electrons in real space as

G (e) µµ ′ (ω) = ⟨ψ|ĉ µ Ĥ -E 0 -ω -1 ĉ † µ ′ |ψ⟩ , (4)
where we combined site and spin labels into greek letters µ = (i, σ), µ ′ = (j, σ ′ ) and the superscript (e) indicates that here we are describing the electron part. The cluster Green's function for holes, indicated by a superscript (h), can be written in a similar manner by exchanging the creation and annihilation operators and replacing ω → -ω.

In the band Lanczos method, the Green's function is calculated by constructing a Krylov subspace using several electron and hole excitations, which serve as initial states.

Representing the initial states as MPS and the Hamiltonian as matrix-product operator (MPO), the Krylov subspace construction can be performed in principle using standard MPS arithmetics. The particular choice of the set of initial states can depend on the cluster symmetries, but for the sake of simplicity we restrict ourselves in the following to the simplest case and ignore symmetries.

Consider the ground state |ψ⟩ of a cluster, given by the Hamiltonian Ĥ acting on L ∈ N cluster sites, where each site j can host electronic degrees of freedom described by annihilation (creation) operators ĉ( †) j,σ . We define a set of initial states |φ (e) j,σ ⟩ = ĉ † j,σ |ψ⟩ and introduce the spin-σ electron Krylov subspace K D σ with D = N • L, generated from N -1 applications of Ĥ to the L initial states with spin σ:

K 0 σ = span |φ (e) 1,σ ⟩ , . . . , |φ (e)
L,σ ⟩ ,

K D σ = K D-1 σ ∪ span ĤN |φ (e) 1,σ ⟩ , . . . , ĤN |φ (5)
L,σ ⟩ . (6) Note that replacing the creation operators with the annihilation operators ĉ † j,σ → ĉj,σ , we obtain another set of initial states |φ (h) j,σ ⟩ = ĉj,σ |ψ⟩, from which we construct the spin-σ hole Krylov subspace KD σ . For each Krylov subspace, an orthonormal basis is obtained from an iterative orthogonalization scheme. Thereby, given a set of orthogonal Krylov vectors |k α,ν ⟩ ∈ K D σ with α ∈ 1, . . . , L ( L ≤ L) and ν ∈ {0, . . . , N -1}, a new set of candidate states {|k α,N ⟩} is generated by applying the Hamiltonian to all states {|k α,N -1 ⟩}. Just as in the conventional Lanczos scheme [START_REF] Cullum | Lanczos Algorithms for Large Symmetric Eigenvalue Computations[END_REF], every candidate state |k α,N ⟩ is then reorthogonalized against its 2L predecessors |k α,N -1 ⟩ and |k α,N -2 ⟩. Furthermore, candidate states have to be orthogonalized against each other, i.e., upon adding a new candidate state |k α,N ⟩, it needs to be orthogonalized against all |k α ′ <α,N ⟩. It may occur that the states |k α,N ⟩ generated by this recursion are not linearly independent. This situation is typically solved by applying a so-called deflation scheme, reducing the range L of the α's until all Krylov states |k α≤ L,N ⟩ are linearly independent [START_REF] Freund | Templates for the Solution of Algebraic Eigenvalue Problems: A Practical Guide[END_REF].

Following this strategy, a global Krylov basis is constructed that allows for a representation of the Green's function Eq. ( 4) in a Krylov subspace

G (e) µµ ′ (ω) ≈ ν,αν ν ′ ,α ′ ν ′ ⟨ψ|ĉ µ |k α ν ,ν ⟩ ⟨k α ′ ν ′ ,ν ′ |ĉ † µ ′ |ψ⟩ × ⟨k α ν ,ν | Ĥ -E 0 -ω -1 |k α ′ ν ′ ,ν ′ ⟩ . (7)
In absence of deflated states, the effective Hamiltonian ⟨k α,ν | Ĥ|k α ′ ,ν ′ ⟩ is block tri-diagonal with a maximum block size 2L + 1 and can be diagonalized easily, to evaluate the operator inverse. In contrast, the usual Lanczos recursion corresponds to the case L = 1: The block size is 2 • 1 + 1 = 3, which reflects the fact that only a single matrix element of the cluster Green's function can be evaluated per Krylov expansion.

The overall dimension of the generated Krylov subspace is given by N • L, which amounts to N -1 applications of the Hamiltonian, and we refer to N as the Krylov order. This, however, comes at the cost of additional, global MPS arithmetics when reorthogonalizing the candidate states. The resulting approximations generally introduce instabilities in the band Lanczos recursion. Instabilities most prominently manifest in form of a loss of the block tri-diagonal structure of the effective Hamiltonian. In other words, weight of the basis states accumulates in the orthogonal complement of the targeted Krylov subspace.

The loss of accuracy of the Krylov subspace expansion leads in particular to two problematic issues. First, the loss of orthogonality amongst the basis states in the Krylov subspace cannot be completely compensated in practice: Due to the large number of basis states, a (full) reorthogonalization would not be efficient numerically. Second, a convergence criterion that relates error bounds of the cluster Green's function to the approximation quality of the Krylov subspace is required: Standard convergence criteria such as the relative change in the energy spectrum are numerically unstable. Our methodical developments target both convergence issues with the goal of a numerically stable band Lanczos recursion using MPS-arithmetics that work at a finite, yet controlled approximation quality.

Energy truncation

Constructing Krylov spaces of higher order, an increasing number of highly excited states has to be captured by the MPS representations of the Krylov vectors. Unfortunately, states from the bulk of the many-body Hamiltonian's spectrum typically satisfy a volume law of entanglement [68] and are therefore the main reason for the exponentially increasing computational complexity. While this is a problem in general, here, we are interested in a faithful approximation of the poles of the Green's matrix, i.e., their positions and weights. Hence, we can exploit a strategy introduced previously in the context of expanding the resolvent in terms of Chebyshev polynomials [58]. The general idea is to construct local Krylov spaces at each lattice site and iteratively project the state's site tensors into a sequence of local Krylov subspaces. However, in contrast to projecting the site tensors into a certain energy window as described in [58], we use the local Krylov space expansion to construct a projector to the subspace containing the most relevant eigenstates, up to a defined threshold δ:

|ψ⟩ = n ⟨ Ẽn |ψ⟩ εn | Ẽn ⟩ ≈ εn>δ ε n | Ẽn ⟩ (8)
where

| Ẽn ⟩ = K l=0 c l |ν l ⟩ , lim K→dim H Ĥ | Ẽn ⟩ = E n |E n ⟩ (9)
and |ν l ⟩ ∈ H K Lanczos vectors in Krylov space H K = span |ψ⟩ , Ĥ |ψ⟩ , . . . , ĤK |ψ⟩ . A schematic illustration of the energy truncation is shown in Fig. 1a. The red shaded bars correspond to those eigenstates that have a weight ⟨ Ẽn |ψ⟩ < δ and are hence discarded. Refering to the expansion of the Green's function in the Krylov subspace Eq. ( 7), we can identify the weights with the expansion of the matrix elements ⟨k α ′ ν ′ ,ν ′ |ĉ † µ |ψ⟩ in the energy eigenbasis. The energy truncation thus simply discards those eigenstates, whose contribution to the total Green's function are O(δ 2 ).

In Fig. 1b we illustrate the effect of the energy truncation on the calculated lattice density of states using the CPT Green's function for a 4 × 4 cluster in the antiferromagnetically ordered state. Setting the truncation threshold too high results in a breakdown of the subspace projection. Here, this is the case for δ = 10 -2 , where the erroneous subspace projection leads to a density of states, which neither captures correctly the gap nor the high-frequency part of the spectrum. However, already for moderate values of δ ≲ 10 -3 the spectrum is well reproduced, in particular the gap and the low-energy part of the spectrum, and only small deviations from the exact result are visible for higher excitation energies. A truncation threshold of δ ≤ 10 -4 finally leads to results which are for practical purposes essentially exact. In this case, only states with negligible weight for the density of states are truncated. The systematic improvements of the high-energy part of the spectrum can be related to the fact that the Green's function is computed from locally exciting the ground state, such that the weight of higher-energy eigenstates decreases rapidly. For the following calculations, we avoid any ambiguity by using an energy truncation threshold of δ = 10 -5 .

The computational costs of this energy truncation are not negligible. The complexity per local update scales as O(Km 3 wd) + O(Km 2 w 2 d 2 ) where K is the number of applications of the local representation of Ĥ. Refering to standard notation, here we denote the MPS bond dimension by m, the MPO bond dimension by w and the local Hilbert space dimension by d. Thus, we use the energy truncation only for the large clusters with L ≥ 16 sites and apply it in the process of creating a new candidate state only once, namely after the global application of Ĥ. However, the goal is not to achieve a speed up at constant m. Instead, we aim for a reduction of orthogonality losses that occur when orthogonalizing candidate states with respect to each other. Here, the effect of the energy truncation is to remove highly entangled contributions allowing for signficantly higher precision in the orthogonalization procedure, while keeping the bond dimension fixed. We found that building up a reasonably large Krylov space of dimension D ∼ 100 for the larger clusters using computationally feasible bond dimensions m = 1024 -2048 was possible only when employing an energy truncation. Otherwise, orthogonality could not be maintained and the construction procedure would break down due to orthogonality losses.

Hochbruck-Lubich criterion

Controlling the convergence of the band Lanczos in the MPS representation is a suprisingly delicate problem. The unavoidable loss of numerical precision, caused by performing global operations using MPS arithmetics, as well as finite truncation errors, render standard convergence criteria such as the relative change in the groundstate energy in the constructed Krylov subspace not suitable. For instance, initially one observes a steep decrease of the relative error δE = (E D+α+1 -E D+α )/E D+α , where E D+α is the ground state energy after adding the αth candidate state to the Krylov subspace K D σ . However, δE exhibits rather drastic jumps over more than an order of magnitude, paired with iterations where δE = 0 w.r.t. the energy resolution, which is bounded by the MPS truncation error. This phenomenon is shown at the example of a 4 × 6-Hubbard cluster in Fig. 2 where the red curves represent δE.

There are two major sources for an erroneous convergence indication in δE. First, whenever Ĥ is applied to a complete sequence of candidate states (here, N • 2L = N • 48 marked by the dashed lines), the Krylov order N increases. Given a fixed Krylov order, the band Lanczos generates a sequence of states that optimize the approximation of the ground state, which, however, is constraint by the Krylov order. Therefore, δE decreases within a sequence of candidate states, which leads to a false indication of overall convergence. Second, candidate states |k α,N ⟩ have a different relevance for the ground state approximation. Most prominently, states that are created from initial excitations located at the boundary of the system converge more quickly.

Besides these technical problems, in practise the energy gain is not the most relevant quantity: In the end, the goal is to approximate the poles and weights of the cluster Green's function up to a certain precision. For that purpose, here we pursue a different approach using a convergence criterion that directly measures the desired approximation quality. The idea is to relate an error measure δR D (τ ), which can be derived for the approximation of the action of the operator exponential e -i Ĥτ |ψ⟩ to the Green's function. The detailed derivation is given in App. B, here we only quote the resulting bound for the error in approximating the Green's function in the Krylov subspace with dimension D:

|δG (e,h) µν (ω)| ≤ 2π 2 δR D (τ ) . ( 10
)
δR D (τ ) can be approximated by the generalized residual introduced by Hochbruck and Lubich [38,69,70], which bounds the actual approximation error of the operator exponential ϵ D ≤ δR D (τ ). An estimation for ϵ D is given by the Hochbruck-Lubich bound [38]

ϵ D τ ≤ 40 τ e -τ 4 eτ 4D D , (11)
relating the approximation error of the cluster Green's function to the Krylov space dimension D and the chosen time step τ . In our computations, we choose ϵ ≤ δR D (τ ) as error threshold, which can be evaluated at negligible numerical costs from the matrix representation of the exponential of the effective Hamiltonian ĤD [34,70] (see also App. B). We furthermore fix the time step by relating it to the spectral width W of Ĥ via τ = 2π κW where typically we choose the energy window fraction κ ≈ 1.

We can now argue why we expect this error measure to be more reliable than only considering the relative change in the ground-state energy: δR D (τ ) is not a relative measure and while we also found it to exhibit fluctuations in our numerics, these are less severe than those in δE. The stability of δR D (τ ) against fluctuations comes from the fact that increasing the Krylov space dimension D, the residual in approximating the time evolution is monotonically decreasing, which follows directly from Eq. (11). Thus, fluctuations can only be generated from approximation errors caused by the MPS arithmetics and truncation errors. The main effect of these errors is to destroy the block-tridiagonal form of the effective Hamiltonian, which, however, can be monitored by the precision of the MPS arithmetics and the truncation errors. We choose the corresponding thresholds to be sufficiently small and thereby obtain a significantly more stable convergence criterion than using δE. It should be noted that strictly speaking, Eq. ( 11) is valid only after a complete sequence of Lanczos states for the set of initial states |φ (e) j,σ ⟩ has been constructed, i.e., the Krylov space dimension has to fulfill D = N • L. We discuss this issue in App. B but from our numerical experience we found that the violation of monotonicity while constructing a new set of L candidate states is not too severe.

Orthogonality loss

A major source for numerical instabilities in Krylov subspace methods is the loss of orthogonality of Krylov vectors. While in the commonly used exact state-representations the loss of orthogonality is in general caused by round-off errors due to finite-precision arithmetics [71], the use of MPS introduces additional error sources. On the one hand, generating a new set of candidate states in the N th iteration of the band Lanczos recursion requires to act with the Hamiltonian on MPS representations of Lanczos vectors |k α,N ⟩: |k α,N ⟩ -→ Ĥ |k α,N ⟩. Performing this operation by a naïve MPO-MPS application is numerically very costly so that we resort to a variational application scheme with a zipup-preconditioner [34]. While this allows us to achieve convergence of the MPO-application typically after two sweeps, we are still limited by the growth of the MPS bond-dimension, requiring a finite truncated weight. On the other hand, the orthogonalization of new candidate vectors Ĥ |k α,N ⟩ against previous Lanczos states has to be done using a variational update scheme at finite truncated weight, too. This introduces a further loss of numerical precision if the variational optimization can get stuck in local minima. We preconditioned the optimization using an optimized initial guess state generation by mixing small contribution of higher and lower order Krylov states into the candidates. Nevertheless, the band Lanczos recursion inevitably is going to loose orthogonality, and the rate at which this loss occurs crucially depends on the maximally allowed truncated weight.

In order to keep the recursions stable as long as possible, we exploit the energy truncation, which removes poles, i.e., eigenstates, with small weights in the Green's function from the Lanczos states. This truncation then allows us to use small truncated weights δ trunc ∼ O(10 -10 ) at moderate bond dimensions. We typically allow for up to χ max = 1024 states and only for the very large clusters increased that value to χ max = 2048. Throughout the band Lanczos recursion, we then monitor the violation of the anti-commutation relations

I j,σ;j ′ ,σ ′ = α,ν ⟨0|ĉ † j,σ |k α,ν ⟩ ⟨k α,ν |ĉ j ′ ,σ ′ |0⟩ + ⟨0|ĉ j ′ ,σ ′ |k α,ν ⟩ ⟨k α,ν |ĉ † j,σ |0⟩ . ( 12
)
We consider a deviation

δ ol = j,j ′ σ,σ ′ ∥I j,σ;j ′ ,σ ′ -δ j,j ′ δ σ,σ ′ ∥ 2 < 10 -3 , (13)
which is of the order of the total truncation error ∼ L • √ δ trunc , as acceptable and terminate the Lanczos recursion if that threshold is exceeded. For all practical applications discussed in this paper, the energy truncation allowed us to perform ∼ O(100) Lanczos iterations without facing numerical instabilities caused by the orthogonality loss. Importantly, this number was sufficient to converge also the largest clusters and obtain faithful approximations of the relevant poles and weights of the cluster Green's function w.r.t. the convergence criterion introduced in Sec. III A 2.

Energy rescaling

In order to conveniently use the error bound Eq. ( 11) we always rescale the Hamiltonian such that the spectral width is given by W ≡ 1. Furthermore, the rescaling allows us to introduce another tool to control the quality of the constructed Krylov subspace. In fact, a severe loss of orthogonality in the band Lanczos is typically signalled by an artificial drop of the lowest eigenvalue ẼD 0 of the effective Hamiltonian in the Krylov space K D below the ground-state energy. For these reasons we rescale and shift the Hamiltonian Ĥ throughout our computations, which of course needs to be compensated for in the evaluation of the cluster Green's function. The first step of the rescaling is to obtain the spectral width

W = |E 0 -E max | (14)
of Ĥ by calculating the ground-state energy E 0 and the energy of the highest excited state E max . Afterwards, the Hamiltonian is rescaled and shifted,

Ĥ → Ĥ = Ĥ -E max |E 0 -E max | , (15)
to force the spectrum to be located within the interval [-1, 0]. Thus, the lowest eigenvalue of the effective Hamiltonian serves as additional proxy to the convergence of the Krylov-space by monitoring whether the condition ẼD 0 > -1 holds true. We also note that using the rescaled Hamiltonian, the Hochbruck-Lubich time step is related trivially to the energy window fraction τ = 2π/κ.

B. VCA

The VCA [20] is an established quantum-cluster technique, which is well-suited to probe correlated systems for symmetry-breaking fields [21,22]. It is based on the self-energy-functional (SEF) theory [72][73][START_REF] Potthoff | Self-energy-functional theory[END_REF] and consists of determining the stationary points of the SEF Ω t [Σ], which is related to the Baym-Kadanoff-Luttinger-Ward functional [START_REF] Baym | [END_REF]76], with respect to trial self-energies of a reference system Σ(t ′ ):

Ω t [Σ (t ′ )] =Ω ′ + Tr ln -G -1 0 -Σ (t ′ ) -1 -Tr ln (-G ′ ) . (16)
Here, Tr A = T ω,α A α,α (iω), Ω ′ denotes the grand potential, and G ′ the Green's function of the reference system. Since the theory requires the reference system to have the same interaction terms as the original system, the cluster self-energies can be varied via their one-body terms t ′ . At a stationary point, the SEF represents an approximation of the grand potential of the original system in the variational space of available self-energies of the reference system.

To investigate symmetry-broken phases, a (fictitious) symmetry-breaking Weiss field term can be added to the cluster Hamiltonian. The Weiss field strength is then one of the one-body terms, which are determined via the variational principle that leads to stationarity of the SEF with respect to the cluster self-energy.

Here, we focus on zero temperature, T = 0, and perform the summation over frequency ω m in Eq. ( 16) analytically. In this case one obtains [23] Tr ln

G -1 0 -Σ -1 = m ω m Θ(-ω m ) -R (17)
Tr ln

G ′ -1 = m ω ′ m Θ(-ω ′ m ) -R , (18)
with the Heaviside function Θ(ω), a contribution R from the poles of the self-energy, the poles ω ′ m of the cluster Green's function G ′ , and the poles ω m of the VCA Green's function (G -1 0 -Σ) -1 . It is therefore useful to work with a Lehmann representation of the Green's function since it explicitly includes the information about its poles ω m ; for more details see appendix A.

In order to limit the variational space of cluster selfenergies, a limited set of one-body parameters of the cluster Hamiltonian are varied in practice: The hopping strengths t cluster ij can be determined from the variational principle to account for cluster boundary effects [20], whereas the cluster chemical potential µ ′ should be optimized to obtain a thermodynamically stable electron filling n [77]. Additional one-body terms can be added to the cluster Hamiltonian to account for symmetrybroken phases such as magnetic order [52], superconductivity [78] or charge order [48]. In the half-filled two-dimensional Hubbard model with isotropic nearest-neighbor hopping, the chemical potential is at µ = U/2 and µ ′ does not need to be determined in the variational search. Furthermore it was shown that optimizing the cluster hopping terms did not lead to significant improvement of the approximation of the grand potential [20]. To capture the essential physics of antiferromagnetic ordering, it is therefore sufficient to apply a staggered magnetic field of strength h z on the cluster, see Eq. ( 3), and use it as the sole one-body parameter throughout the variational search.

C. Choice of clusters

The choice of a cluster is a centerpiece of most quantum cluster techniques be it for accessing specific cluster momenta or for performing a finite-size scaling. Besides the size of the cluster, its geometry is of crucial importance in two-dimensional systems and the seminal papers of Betts et al. [79,80] introduced a set of evaluation criteria for the suitability of clusters with periodic boundary conditions, which were based on the completeness of near-neighbor shells.

Within CPT and VCA one uses open boundary conditions on the clusters and a systematic study of Betts clusters is missing. Nevertheless, we resort here to the labelling and characterization scheme introduced by Betts et al. [81], in particular to parameters such as the geometrical imperfection J and the bipartite imperfection I B to identify clusters of high geometrical and topological quality. Apart from more traditional rectangular clusters, we also test our solver on skewed Betts clusters. We limit ourselves to clusters that can tile bipartite lattices and focus on those which have zero bipartite imperfection. By combining calculations of several cluster geometries, we can access different points of the lattice Brillouin zone via reciprocal vectors of the superlattice. In Fig. 3 we illustrate this for lattice tilings using different 24-site clusters that fulfil the condition I B = 0.

IV. RESULTS AND DISCUSSION

In the following, we study the antiferromagnetic (AF) ordering within the half-filled single-band Hubbard model to benchmark our solver. We focus here at the challenging point in parameter space U/t = 8, i.e., an on-site interaction equal to the band-width of the noninteracting dispersion, where the system is already in the Mott insulating regime.

Translated to the self-energy functional, this leads to one stationary point, a maximum, at h z = 0 corresponding to a non-magnetic solution and to two minima at h z = ±h c z ̸ = 0 corresponding to a phase with antiferromagnetic order, see Fig. 4. Since Ω is symmetric in h z , we plot it in the following only for positive values of h z . Compared to numerically exact solvers, the MPSbased solver introduces additional approximations, in particular since the MPS bond dimension is bound to a maximal value χ max and since perfect orthogonality between the Krylov vectors is lost faster at large iteration numbers within the band-Lanczos procedure. Comparing the self-energy functional computed with our MPSbased solver to the benchmark using an exact diagonalization (ED) based solver, we see excellent agreement in Fig. 4 if the maximal MPS bond dimension is chosen large enough. For the small cluster used here, χ can be increased until the MPS representation is essentially exact. However, for too small bond dimension, deviations from the exact curve are visible. In that case, the selfenergy functional is systematically evaluated to be too large. Except for extreme cases, e.g. χ = 128 in Fig. 4, Ω shows even for moderate bond dimension correct functional behavior as a function of the Weiss field, which allows for an interpolation around the stationary points. As we show in Sec. IV 1, fitting the self-energy functional to determine its stationary points allows to calculate observables for different bond dimension in order to perform an extrapolation to χ max → ∞.

To illustrate the effects of cluster geometry on Ω, we plot in Fig. 5 the self-energy functional of the 24-site clusters introduced in Fig. 3. Even though the rough functional forms of Ω do not differ qualitatively, their stationary points vary notably. We also note that clusters with similar geometric properties, as for instance quantified in terms of their geometrical and bipartite imperfection, see Fig. 3(B), generate quite similar self-energy functionals, whereas different imperfection of two clusters leads to significantly different Ω.

Another consequence of using a MPS representation of the quantum states is that even for a given cluster geometry, the obligatory mapping of the 2D cluster onto a 1D chain introduces a degree of freedom which has a direct impact on the precision with which the SEF is calculated. In the inset of Fig. 5 we illustrate this effect by calculating the cluster Green's function of the rectangular 24s0 cluster for two different 1D mappings. The corresponding self-energy functionals agree in the χ → ∞ limit, but differ for finite χ. Whereas the mapping minimizing longrange hopping ('4 × 6') leads to well converged results at χ max = 1024, the longer-ranged mapping ('6 × 4') leads to a larger build-up of entanglement and would require a χ max > 2048.

In presence of nearly degenerate quantum states, different symmetry-breaking orders on the cluster can be easily picked up within MPS based approaches when not providing a sufficiently large bond dimension. This is for instance the case for the '4 × 6' mapped 24s0 cluster at h z ≈ 0.11, see inset of Fig. 5. Hence, scaling the position of the stationary points of Ω in the maximal MPS bond dimension instead of individual points of the self-energy functional proved to be more robust. Around the stationary point Ω(h c z ) the functional can be fitted as a function of the Weiss field, shown as a line (A); the corresponding minima are indicated by black circles. For different χ max , the SEF shows similar functional dependence on hz. The staggered magnetization ⟨m⟩ on the lattice (B) and the Weiss field strength h c z (C) are then scaled as a function of inverse bond dimension 1/χ max .

SEF scaling with bond dimension

In Fig. 6, the value of the Weiss field h c z as well as the staggered magnetization on the lattice, ⟨m⟩, are scaled in the inverse bond dimension. The staggered magnetization scales linear in 1/χ max and gives a lower estimate of the expectation value for the χ max → ∞ limit. Since at some critical χ c the MPS representation is exact, we expect a constant value of ⟨m⟩ for χ max > χ c . The value of ⟨m⟩ at the largest calculated χ max therefore sets an upper bound for the staggered magnetization.

However, the dependence of h c z and ⟨m⟩ on 1/χ is mild (O(10 -3 )) and of the same order as the estimated error such that it can be neglected in practice for most applications as long as a sufficiently large bond dimension χ is used. In the case of the shown 22-site Betts cluster, this would be the case χ max ≳ 1000. Given the large computational cost of evaluating the cluster Green's function via the MPS-based band Lanczos solver, systematically scaling expectation values in the MPS bond dimension should remain the exception. Instead, it offers the possibility to calculate spectral functions or observables for large clusters (e.g. using cluster perturbation theory), for which symmetry breaking terms can be included via VCA. The rather mild dependence of the functional form of the SEF and the related expectation values on the chosen MPS bond dimension is important in this context.

Spectral functions

Next, we discuss the mild influence of the maximal MPS bond dimension χ max on the spectral function, see Fig. 7.

For comparatively small clusters, the spectral function A(k, ω) and the density of states (DOS) N (ω) can be converged in χ max . We illustrate this in the case of the 4 × 4 cluster for N (ω) in Fig. 7(a). The results obtained from restricting the bond dimension to χ max = 512 show already good agreement for excitation energies close to the gap, but the density of states still shows differences at larger excitation energies around ω ∼ ±6. For χ max ≳ 1024, the density of states and spectral function of the 4 × 4 cluster are converged and do not change any more when increasing the bond dimension.

For larger clusters like the 22h5 cluster discussed previously in Fig. 6, scaling h c z and ⟨m⟩ still leads to minor changes when increasing the bond dimension up to χ max = 2048. However, tracking the changes of the density of states as a function of χ max reveals that only small changes close to the outmost band edges far from the gap are visible. The DOS for excitation energies |ω| ≲ 5 is already converged for χ max = 1536 and even for χ max = 1024 the deviations from the converged solution are for most purposes negligible.

In Fig. 7(b) we finally show the k-resolved spectral function A(k, ω) along the high symmetry path (0, 0)-(π, 0)-(π, π)-(0, 0) for one of the largest clusters we systematically studied, the 6 × 6 cluster. The spectral function was calculated from the CPT Green's function, Eq.( 2), using the periodization scheme proposed in Ref. 19. In contrast to small and intermediate clusters, the SEF of the 6 × 6 cluster still shows notable differences as a function of the MPS bond dimension up to χ max = 2048. They translate to the spectral function in form of small differences in the low-energy excitations, most visible around the X point, k = (π, 0). Nevertheless, the most salient features of the spectral function are already converged in the bond dimension for a remarkably small χ max = 1536. This can be seen as an advantageous feature of the MPS solver, which enforces due to the energy truncation scheme first the convergence of poles of the Green's function that contribute most spectral weight. The present MPS solver therefore even allows to interpret spectral functions of large clusters for which a full convergence in terms of the MPS bond dimension cannot be achieved.

Finite size scaling

In order to perform a finite-size scaling of the staggered magnetization, calculations need to be performed using different cluster sizes. However, even for fixed cluster size, the cluster geometry has an impact on the selfenergy functional since the corresponding cluster Green's functions include long-range correlations to a different extent. To include this information for the finite-size scaling, we use the quality factor introduced by D. Sénéchal in Ref. [START_REF]An introduction to quantum cluster methods[END_REF]. It amounts to comparing the intra-cluster hopping terms N cl t to the number of hopping terms per site in the infinite lattice and therefore takes into account both the cluster size and geometry. In the square lattice, this amounts to q = N cl t /2L. The scaling factor q thereby contrasts with other commonly used scaling parameters like 1-1/L, which include only information of the cluster The staggered magnetization ⟨m⟩ and the magnetic Weiss field on the cluster hz (inset) are scaled as a function of the parameter q described in text.

The reference data from auxiliary-field quantum Monte Carlo (AFQMC) is taken from Ref. [82]. The staggered magnetization extrapolated from the band-Lanczos results for rectangular clusters (highlighted in red) using a non-linear ansatz is shown as dashed blue line and the shaded area indicates the 1σ confidence interval, containing the AFQMC result.

size.

When using this quality factor, the extrapolation of the Weiss field strength h c z to the thermodynamic limit leads to a vanishing field strength whereas the expectation value of the staggered magnetization on the lattice ⟨m⟩ extrapolates to a finite value, see Fig. 8. In the infinite-cluster-size limit, q = 1 and we will not need the symmetry-breaking field since the symmetry will be broken spontaneously. There, the value of the staggered magnetization ⟨m⟩| q=1 simply corresponds to the order parameter on the lattice. At half-filling, ⟨m⟩ can be calculated numerically exactly using auxiliary-field quantum Monte Carlo simulations (AFQMC) [START_REF] Hirsch | [END_REF][85][86][87][88][START_REF] Becca | Quantum Monte Carlo approaches for correlated systems[END_REF].

Comparing to the AFQMC result (dashed black line) of Seki and Sorella [82] we see that the linear scaling in q (dashed, violet graph) would lead to an overestimation by ∼ 20% of the staggered magnetization, i.e., the linear dependence ⟨m⟩ on q seems to hold only for small and intermediate cluster sizes. Indeed, for the largest clusters calculated with our MPS band-Lanczos solver, we already observe a significant deviation from linear scaling (c.f. the red highlighted symbols in Fig. 8). We believe that besides others, the non-linear q-dependence is generated from the fact that in the presented calculations we did not vary the hopping amplitudes at the cluster boundaries. In that case, an extrapolation would still be possible even though the exact functional dependency is unknown.

We attempt to account for the non-linear behavior by: (1) Linearly extrapolating ⟨m⟩ (q → 1, q 0 ) for the largest rectangular clusters (c.f. the red highlighted symbols in Fig. 8) on a domain [q 0 , 1] where we varied q 0 ∈ [0.775, 0.8125]. The obtained linear extrapolations are indicated by the colored dashed lines in Fig. 8. (2) We then make the simplest non-linear ansatz for the functional dependency on q 0 ⟨m⟩ (q → 1, q 0 ) = α(q 0 -1) 2 + m ∞ .

(

) 19
This ansatz has a stationary point at q 0 = 1, resembling the expected asymptotic behavior of the magnetization when scaling the number of cluster sites to the thermodynamic limit. The extrapolated value of the magnetization is m ∞ = 0.607 ± 0.009, which is shown in Fig. 8 as blue dashed line and the shaded area indicates the 1σ confidence interval of the extrapolation. Note that this value agrees with the AFQMC-result m QMC ∞ = 0.5982 within the confidence interval and the deviation is only about 1.5%. We therefore find that with the cluster sizes accessible via the VCA+MPS ansatz, it is possible to estimate values for local observables in the thermodynamic limit to a very high precision, even without exploring all possible variational parameters within VCA.

V. CONCLUSIONS AND OUTLOOK

We introduce a band-Lanczos solver based on MPS for quantum cluster methods such as the VCA. While a naïve implementation of MPS as cluster solver does not lead to a substantial improvement over using ED techniques, we showed that significantly larger cluster sizes can be treated with high accuracy by combining the MPS+band-Lanczos solver with a suitable energy truncation, controlling the loss of orthogonality of the Lanczos states and introducing a robust convergence criterion.

We demonstrate the potential of the approach by computing the staggered magnetization ⟨m⟩ of a single band Hubbard model on a 2D square lattice. Treating different cluster sizes and shapes (regular rectangular as well as Betts clusters), we extrapolate our MPS data in the bond dimension χ, which we show to yield a systematic improvement of the approximation quality of the Green's function and its derived quantities. In this way we can exploit the increased number of accessible cluster sizes and geometries to extrapolate observables such as the staggered magnetization into the thermodynamic limit. We show that at half filling and for intermediate interaction strengths U/t = 8 (i.e., in the Mott insulating regime, where the interaction strength is of the order of the bandwidth of the noninteracting system) an anomalous finite-size scaling with a crossover length scale is observed. Only for rectangular cluster sizes that are large enough, i.e., 4 × 4 and larger, a deviation from an otherwise linear scaling is obtained, which can be extrapolated towards the correct result in the thermodynamic limit. The resulting value of ⟨m⟩ in the thermodynamic limit agrees with numerically exact QMC data within an error margin of 1.5%. In contrast, with the smaller cluster sizes amenable to ED, the extrapolation in cluster size leads to a value, which overestimates the QMC result by ∼ 20%.

Given the significant increase of the number of cluster sites, we envisage that our approach has the potential to substantially improve the investigations in multi-band Hubbard-like models. Here, Hund's coupling leads to an interplay between orbital, spin, and electronic degrees of freedom [START_REF] Georges | Annual Review of Condensed Matter Physics[END_REF][START_REF] De' Medici | The Physics of Correlated Insulators, Metals, and Superconductors[END_REF], which affects material properties like superconductivity in nickelates and iron based superconductors, orbital selectivity in iron chalcogenides or optoelectronic and photovoltaic properties [START_REF] Lanatà | [END_REF][93][94][95][96][97][98][99]. Although the interactions are local in nature, the higher complexity of these 2D systems poses an even stronger challenge to their theoretical description. While it is essentially impossible for ED-based solvers to treat cluster sizes large enough to perform meaningful studies on such multi-or-bital systems, the MPS+band-Lanczos approach makes it possible to treat also systems with larger local Hilbert spaces in a controlled way, so that it appears natural to extend our approach to such multi-band situations.

  FIG. 1. (a) Schematic representation of the energy truncation. Overlaps ϵn = ⟨ Ẽn|ψ⟩ < δ are neglected when constructing the sublattice projection (marked red). (b) Effect of the truncation threshold on the (lattice) density of states N (ω) for a 4 × 4 cluster in the AF phase, see Sec.III B for details. The results for δ = 10 -2 , 10 -3 are shifted for clarity.

FIG. 2 .

 2 FIG. 2. Comparison of convergence critera. The blue curve (right y-axis) shows the lowest eigenvalue E D of the effective Hamiltonian in the Krylov subspace as a function of the dimension D for a 4 × 6-site Hubbard cluster at U = 8t and a staggered magnetic field hz = 0.09t. The conventional convergence criterion δE = E D+1 -E D is given by the red line (left y-axis), exhibiting strong fluctuations over several orders of magnitude. In contrast, oscillations in the suggested Hochbruck-Lubich critertion ϵ D shown by the green line are significantly smaller, overlaying a systematic decrease w.r.t. the number of iterations.

FIG. 3 .

 3 FIG. 3. (A) Illustration of the geometries of different clusters with N = 24 sites labeled according to Betts et al. [81]. Red arrows indicate the superlattice vectors that tile the lattice with the respective cluster, dashed lines indicate inter-cluster hopping terms and the Néel-type sublattice structure is indicated by different tones of grey. The clusters differ in their squareness σ, geometrical imperfection J and bipartite imperfection IB, see (B). (C) Points in the Brillouin zone that are multiples of reciprocal superlattice vectors, which are indicated by arrows. The corresponding reduced (superlattice) Brillouin zones are marked by parallelograms.

FIG. 4 .

 4 FIG. 4. The self-energy functional Ω as a function of staggered Weiss field strength for a 2×4 cluster. Stationary points of Ω are indicated by arrows. Dots denote results obtained with the matrix-product state (MPS)-based solver, in the inset shown for different maximal MPS bond dimension χ max around the minimum of Ω; the results obtained with an exact diagonalization (ED) solver are shown as a line.

FIG. 5 .

 5 FIG.5. Impact of the cluster geometry on the self-energy functional. Ω is calculated with χ max = 1024 for different 24-site clusters shown in Fig.3. Minima of Ω are indicated by arrows, lines are only a guide to the eye. The inset shows Ω for two different 1D mappings of the 24s0 cluster, '4 × 6' and '6 × 4', calculated with χ max = 1024, 2048.

-

FIG. 7 .

 7 FIG. 7. (a) Density of states N (ω) of the 4×4 (top) and 22h5 (bottom) clusters for different maximal MPS bond dimension χ max . (b) shows the k-resolved spectral function A(k, ω) of a 6 × 6 cluster for two different maximal bond dimensions. All calculations were done in the antiferromagnetically ordered state, i.e. at the stationary point of the SEF for the respective clusters.

  Scaling of the self-energy functional Ω in the maximal MPS bond dimension χ max for the Betts cluster 22h5.

	(A)
	z ) self-energy functional Ω(h	4.5045 -4.5044 -4.5043 -4.5042	χ max 700 800 1024 1200 1536 2048
		0.09	0.1	0.11	0.12	0.13
			staggered Weiss field strength h z
	(B)	0.77
	m	0.768
		0.766
		0	0.0005	0.001		0.0015
	(C)	0.107
	h c z	0.106
		0.105
		0	0.0005	0.001		0.0015
			1/χ max
	FIG. 6.

ACKNOWLEDGMENTS

We thank the Grand Equipement National de Calcul Intensif for providing supercomputing time at IDRIS and TGCC (Project No. A0130912043). We also acknowledge a mobility grant of the Franco-Bavarian University cooperation centre (BayFrance). Furthermore funding through the ERC Starting Grant from the European Union's Horizon 2020 research and innovation programme under grant agreement No 758935 is greatfully acknowledged. S. P. acknowledges support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy-426 EXC-2111-390814868. Support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) -207383564/FOR1807 project P7 at an early stage of this work is acknowledged.

Supplemental Materials: Matrix-product-state-based band-Lanczos solver for quantum cluster approaches

Appendix A: Q-matrix formulation of the cluster Green's function

To rewrite the cluster Green's function, we use the Q-matrix formulation of the Lehmann representation introduced in Ref. [100]:

with α, β being compound site and spin indices, and excited states m. The Q-matrix is defined as the concatenation of the electron and hole matrices

where |0⟩ is the ground state with energy E 0 . The corresponding excitation energies needed in Eq. ( 18) read

= (E r -E 0 , E 0 -E s ). To obtain the poles ω for Eq. ( 17) we rewrite the VCA Green's function as

where g -1 = ω -Λ with Λ m,l = δ m,l ω ′ m so that the poles ω m are nothing but the eigenvalues of Λ+Q † VQ. The Qmatrix and the one-particle excitations ω ′ are calculated via the band Lanczos algorithm.

Appendix B: Hochbruck-Lubich convergence criterion

We consider the Greens function in its time integral representation

for some Hamiltonian Ĥ and local perturbations of an initial state |ψ⟩ generated from operators Â, B:

where η > 0 is a finite broadening. The error estimation by Hochbruck and Lubich is based on approximating the action of an operator exponential e -t  in a Krylov subspace K D where D ∈ N denotes the dimension of the Krylov subspace

Decomposing the expectation value F (t) = ⟨φ|e -i Ĥτ |φ⟩ into its components the Krylov subspace and the orthogonal complement using Eq. (B.4)

and abbreviating Û D (t, ω) = e -i( ĤD -ω)t , we can rewrite the time integration as

with the residual

We can now apply Eq. (B.3) to estimate the residual

where we used the monotonicity of the generalized residual δR D (t) ≤ δR D (t+pτ ) ≤ 1 for t ≤ τ . In our numerical calculations we usually set τ = 2π such that we can determine the correction term J(η) ≡ 1-e -2πη (1+2πη)

and thus we arrive at the error bound

The generalized residual δR D (τ ) can be computed at no additional cost using the approximation of the exponential in the Krylov subspace

While these considerations provide a reliable bound for the usual Lanczos procedure, they are exact in case of the band Lanczos only, if the Krylov dimension D and the Krylov order N , i.e., the number of applications of Ĥ to the set of L initial states, are related via D = N • L. The reason behind this limitation is that in Eq. (B.6) we assumed that the residual Krylov vector |f D ⟩ has a Krylov order N +1 while ĤD is approximated with Krylov order N , only. In princple this means that Eq. (B.11) can only be applied after a whole sequence of L candidate states has been constructed. However, in practise we observed that also within such a sequence, Eq. (B.11) yields satisfying error estimations with fluctuations in the error estimation being of an order of magnitude, at the most.