{'fraction_non_alphanumeric': 0.049807428683334425, 'fraction_numerical': 0.011717475031007245, 'mean_word_length': 3.940180586907449, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fast Algorithm for Finding the Eigenvalue Distribution of Very Large Matrices

A theoretical analysis is given of the equation of motion method, due to Alben et al., to compute the eigenvalue distribution (density of states) of very large matrices. The salient feature of this method is that for matrices of the kind encountered in quantum physics the memory and CPU requirements of this method scale linearly with the dimension of the matrix. We derive a rigorous estimate of the statistical error, supporting earlier observations that the computational efficiency of this approach increases with matrix size. We use this method and an imaginary-time version of it to compute the energy and the specific heat of three different, exactly solvable, spin-1/2 models and compare with the exact results to study the dependence of the statistical errors on sample and matrix size.


I. INTRODUCTION
The calculation of the distribution of eigenvalues of very large matrices is a central problem in quantum physics. This distribution determines the thermodynamic properties of the system (see below). It is directly related to the single-particle density of states (DOS) or Green's function. In a one-particle (e.g., one-electron) description knowledge of the DOS suffices to compute the transport properties [1].
The most direct method to compute the DOS, i.e. all the eigenvalues, is to diagonalize the matrix H representing the Hamiltonian of the system. This approach has two obvious limitations: The number of operations increases as the third power of the dimension D of H and, perhaps most importantly, the amount of memory required by state-of-the-art algorithms grows as D 2 [2,3]. This scaling behavior limits the application of this approach to matrices of dimension D = O(10000), which is too small for many problems of interest. What is needed are methods that scale linearly with D.
There has been considerable interest in developing "fast" (i.e. O(D)) algorithms to compute the DOS and other similar quantities. One such algorithm and an application of it to electron motion in disordered alloy models was given by Alben et al. [4]. In this approach the DOS is obtained by solving the time-dependent Schrödinger equation (TDSE) of a particle moving on a lattice, followed by a Fourier transform of the retarded Green's function [4]. Using the unconditionally stable split-step Fast Fourier Transform (FFT) method to solve the TDSE, it was shown that the eigenvalue spectrum of a particle moving in continuum space can be computed in the same manner [5]. Fast algorithms of this kind proved useful to study various aspects of localization of waves [6][7][8] and other oneparticle problems [9,10]. Application of these ideas to quantum many-body systems triggered further development of flexible and efficient methods to solve the TDSE. Based on Suzuki's product formula approach, an unconditionally stable algorithm was developed and used to compute the time-evolution of two-dimensional S=1/2 Heisenberg-like models [11]. Results for the DOS of matrices of dimension D ≈ 1000000 where reported [11]. A potentially interesting feature of these fast algorithms is that they may run very efficiently on a quantum computer [12,13].
A common feature of these fast algorithms is that they solve the TDSE for a sample of randomly chosen initial states. The efficiency of this approach as a whole relies on the hypothesis (suggested by the central limit theorem) that satisfactory accuracy can be achieved by using a small sample of initial states. Experience not only shows that this hypothesis is correct, it strongly suggests that for a fixed sample size the statistical error on physical quantities such as the energy and specific heat decreases with the dimension D of the Hilbert space [12].
In view of the general applicability of these fast algorithms to a wide variety of quantum problems it seems warranted to analyze in detail their properties and the peculiar D dependence in particular. In Sections II and III we recapitulate the essence of the approach. We present a rigorous estimate for the mean square error (variance) on the trace of a matrix. In Section IV we describe the imaginary-time version of the method. The statistical analysis of the numerical data is discussed in Section V. Section VI describes the model systems that are used in our numerical experiments. The algorithm used to solve the TDSE is reviewed in Section VII. In Section VIII we derive rigorous bounds on the accuracy with which all eigenvalues can be determined and demonstrate that this accuracy decreases linearly with the time over which the TDSE is solved. The results of our numerical calculations are presented in Section IX and our conclusions are given in Section X.

II. THEORY
The trace of a matrix A acting on a D-dimensional Hilbert-space spanned by an orthonormal set of states {|φ n } is given by Note that according to (1) we have Tr 1 = D. If D is very large one might think of approximating Eq. (1) by sampling over a subset of K (K ≪ D)"important" basis vectors. The problem with this approach is that the notion "important" may be very model-dependent. Therefore it is better to sample in a different manner. We construct a random vector |ψ by choosing D complex random numbers, c n ≡ f n + ig n , with mean 0, for n = 1 . . . D, so and calculate If we now sample over S realizations of the random vectors {ψ} and calculate the average, we obtain Assuming that there is no correlation between the random numbers in different realizations and that the random numbers f n,p and g n,p are drawn from an even and symmetric (both with respect to each variable) probability distribution (see Appendix A for more details), we have where E (.) denotes the expectation value with respect to the probability distribution used to generate the c n,p 's. In the r.h.s of (5) the subscripts of c n,p have been dropped to indicate that the expectation value does not depend on n or p. It follows immediately that showing that we can compute the trace of A by sampling over random states {ψ p }, provided there is an efficient algorithm to calculate ψ p |Aψ p (see Section VII). According to the central limit theorem, for a large but finite S we have meaning that the statistical error on the trace vanishes like 1/ √ S, which is not surprising. What is surprising is that one can prove a much stronger result as follows. Let us first normalize the c n,p 's so that, for all p, D n=1 |c n,p | 2 = 1.
This innocent looking step has far reaching consequences. First we note that the normalization renders the method exact in (the rather trivial) case that the matrix A is proportional to the unit matrix. The price we pay for this is that for fixed p, the c n,p are now correlated but that does not cause problems (see Appendix A). Second it follows that E (|c| 2 ) = 1/D. Obviously the error can be written as where is a traceless (due to Eq. (8)) Hermitian matrix of random numbers. We put X = Tr RA and compute E(|X| 2 ). The result for the general case can be found in Appendix A. For a uniform distribution of the c n,p 's on the hyper-sphere defined by D n=1 |c n,p | 2 = 1 the expression simplifies considerably and we find an exact expression for the variance in terms of the sample size S, the dimension D of the matrix A and the (unknown) constants Tr A † A and | Tr A|. Invoking a generalization of Markov's inequality [14] where P(Q) denotes the probability for the statement Q to be true. We find that the probability that | Tr RA| 2 exceeds a fraction a of | Tr A| 2 is bounded by or, in other words, the relative statistical error e A on the estimator of the trace of A is given by if | Tr A| > 0. We see that e A = 0 if A is proportional to a unit matrix. From (14) it follows that, in general, we may expect e A to vanish with the square root of SD. The prefactor is a measure for the relative spread of the eigenvalues of A and is obviously model dependent.
The dependence of e A on S, D and the spectrum of A is corroborated by the numerical results presented below. It is also of interest to examine the effect of not normalizing the c n,p 's. A calculation similar to the one that lead to the above results yields Clearly this bound is less sharp and does not vanish if A is proportional to a unit matrix.

III. REAL-TIME METHOD
The distribution of eigenvalues or density of states (DOS) of a quantum system is defined as where H is the Hamiltonian of the system and n runs over all the eigenvalues of H. The DOS contains all the physical information about the equilibrium properties of the system. For instance the partition function, the energy, and the heat capacity are given by respectively. Here β = 1/k B T and k B is Boltzmann's constant (we put k B = 1 and = 1 from now on).
As explained above the trace in the integral (16) can be estimated by sampling over random vectors. For the statistical error analysis discussed below it is convenient to define a DOS-per-sample by where the subscript p labels the particular realization of the random state |ψ p . The DOS is then given by Schematically the algorithm to compute d p (ǫ) consists of the following steps: 1. Generate a random state |ψ p (0) , set t = 0.
4. Solve the TDSE for a small time step τ , replacing |ψ p (t) by |ψ p (t + τ ) (see Section VII for model specific details).

Repeat N times from
Step 3.
6. Perform a Fourier transform on the tabulated result and store d p (ǫ).
In practice the Fourier transform in Eq. (16) is performed by the Fast Fourier Transform (FFT). We use a Gaussian window to account for the finite time τ N used in the numerical time-integration of the TDSE. The number of time step N determines the accuracy with which the eigenvalues can be computed. In Section VIII we prove that this systematic error in the eigenvalues vanishes as 1/τ N.
Since for any reasonable physical system (or finite matrix) the smallest eigenvalue E 0 is finite, for all practical purposes d p (ǫ) = 0 for ǫ < ǫ 0 < E 0 . The value of ǫ 0 is easily determined by examination of the bottom of spectrum. To compute Z, E, or C we simply replace the interval [−∞, +∞] by [ǫ 0 , +∞].

IV. IMAGINARY-TIME METHOD
The real-time approach has the advantage that it yields information on all eigenvalues and can be used to compute both dynamic and static properties without suffering from numerical instabilities. However for the computation of the thermodynamic properties, the imaginary-time version is more efficient. We will use the imaginary-time method as an independent check on the results obtained by the real-time algorithm.
Repeating the steps that lead to Eq. (17) we find with similar expressions for E and C. Furthermore we have assuming H is Hermitian. Therefore we only need to propagate the random state for an imaginary time β/2 instead of β. Furthermore we do not need to perform an FFT. Disregarding these minor differences, the algorithm is the same as in the real-time case with τ replaced by −iτ .

V. ERROR ANALYSIS
Estimating the statistical error on the partition function Z is easy because it depends linearly on the trace of the (imaginary) time evolution operator. However the error on E and C depends on this trace in a more complicated manner and this fact has to be taken into account.
First we define for the real-time method and for the imaginary-time method.
For each value of β we generate the data {z p }, {h p }, and {w p }, for p = 1, . . . , S. For both cases we have The standard deviations on z, h, and w are given by where var(x) ≡ x 2 − x 2 denotes the variance on the data {x p }. However the sets of data {z p }, {h p } and {w p } are correlated since they are calculated from the same set {|ψ p }. These correlations in the data are accounted for by calculating the covariance matrix M k,l (k, l = 1, . . . , 3) the elements of which are given by are a shorthand for {z p }, {h p }, and {w p } respectively. The estimates for the errors in Z, E and C are given by where E = h/z and and C = β 2 (w/z − h 2 /z 2 ).

VI. EXACTLY SOLVABLE SPIN 1/2 MODELS
The most direct way to assess the validity of the approach described above is to carry out numerical experiments on exactly solvable models. In this paper we consider three different exactly solvable models, two spin-1/2 chains and a mean-field spin-1/2 model. The former have a complicated spectrum, the latter has a highly degenerate eigenvalue distribution. These spin models differ from those studied elsewhere [11,12] in that they belong to the class of integrable systems.

A. Spin chains
Open spin chains of L sites described by the Hamiltonian where σ x i , σ y i , and σ z i denote the Pauli matrices and J, ∆ and h are model parameters, can be solved exactly. They can be reduced to diagonal form by means of the Jordan-Wigner transformation [15]. We have where c + i and c i are spin-less fermion operators and are L × L matrices. By further canonical transformation this Hamiltonian can be written as where n k is the number operator of state k and the Λ k 's are given by the solution of the eigenvalue equation In the general case this eigenvalue problem of the L × L Hermitian matrix ( is most easily solved numerically. In the present paper we confine ourselves to two limiting cases: The XY model (∆ = 1) and the Ising model in a transverse field (∆ = 0).

B. Mean field model
The Hamiltonian of the mean-field model reads and can be rewritten as The single spin-L/2 Hamiltonian has eigenvalues with degeneracy This rather trivial model serves as a test for the case of highly degenerate eigenvalues.

VII. TIME EVOLUTION
For the approach outlined in Sections III and IV to be of practical use it is necessary that the matrix elements of the exponential of H can be calculated efficiently. The purpose of this section is to describe how this can be done.
The general form of the Hamiltonians of the models we study is where the first sum runs over all pairs P of spins, σ α i (α = x, y, z) denotes the α-th component of the spin-1/2 operator representing the i-th spin. For both methods, we have to calculate the evolution of a random state, i.e. U(τ )|ψ ≡ exp(−iτ H)|ψ or U(τ )|ψ ≡ exp(−τ H)|ψ for the real and imaginary time method respectively. We will discuss the real-time case only, the imaginary-time problem can be solved in the same manner.
Using the semi-group property Then the main step is to replace U(τ ) by a symmetrized product-formula approximation [17]. For the case at hand it is expedient to take where Other decompositions [11,18] work equally well but are somewhat less efficient for the cases at hand. In the real-time approach U (τ ) is unitary and hence the method is unconditionally stable [17] (also the imaginary-time method can be made unconditionally stable). It can be shown that U(τ ) − U (τ ) ≤ sτ 3 (s > 0 a constant) [19], implying that the algorithm is correct to second order in the time step τ [17]. Usually it is not difficult to choose τ so small that for all practical purposes the results obtained can be considered as being "exact". Moreover, if necessary, U (τ ) can be used as a building block to construct higher-order algorithms [20][21][22][23]. In Appendix B we will derive bounds on the error in the eigenvalues when they are calculated using a symmetric product formula. As basis states {|φ n } we take the direct product of the eigenvectors of the S z i (i.e. spin-up |↑ i and spin-down |↓ i ). In this basis, e −iτ Hz/2 changes the input state by altering the phase of each of the basis vectors. As H z is a sum of pair interactions it is trivial to rewrite this operation as a direct product of 4x4 diagonal matrices (containing the interaction-controlled phase shifts) and 4x4 unit matrices. Still working in the same representation, the action of e −iτ Hy/2 can be written in a similar manner but the matrices that contain the interactioncontrolled phase-shift have to be replaced by non-diagonal matrices. Although this does not present a real problem it is more efficient and systematic to proceed as follows. Let us denote by X(Y ) the rotation by π/2 of each spin about the x(y)-axis. As it is clear that the action of e −iτ Hy/2 can be computed by applying to each spin, the inverse of X followed by an interaction-controlled phase-shift and X itself. The prime in (53) indicates that J z i,j and h z i in H z have to be replaced by J y i,j and h y i respectively. A similar procedure is used to compute the action of e −iτ Hx . We only have to replace X by Y .

VIII. ACCURACY OF THE COMPUTED EIGENVALUES
First we consider the problem of how to choose the number of time steps N to obtain the DOS with acceptable accuracy. According to the Nyquist sampling theorem employing a sampling interval ∆t = π/ max i |E i | is sufficient to cover the full range of eigenvalues. On the other hand the time step also determines the accuracy of the approximation U (τ ). Let us call the maximum value of τ which gives satisfactory accuracy τ 0 (for the imaginary-time method, this is the only parameter). For the examples treated here τ 0 < ∆t), implying that we have to use more steps to solve the TDSE than we actually use to compute the FFT. Eigenvalues that differ less than ∆ǫ = π/N∆t cannot be identified properly. However since ∆ǫ ∝ N −1 we only have to extend the length of the calculation by a factor of two to increase the resolution by the same factor.
At first glance the above reasoning may seem to be a little optimistic. It apparently overlooks the fact that if we integrate the TDSE over longer and longer times the error on the wave function also increases (although it remains bounded because of the unconditional stability of the product formula algorithm). In fact it has been shown that in general [17] where t = mτ , suggesting that the loss in accuracy on the wave function may well compensate for the gain in resolution that we get by using more data in the Fourier transform. Fortunately this argument does not apply when we want to determine the eigenvalues as we now show. As before we will discuss the real-time algorithm only because the same reasoning (but different mathematical proofs) holds for the imaginary-time case. Consider the time-step operator (52). Using the fact that any unitary matrix can be written as the matrix exponential of a Hermitian matrix we can write U (τ ) = e −iτ Hz/2 e −iτ Hy/2 e −iτ Hx e −iτ Hy/2 e −iτ Hz/2 ≡ e −iτ H(τ ) .
It is clear that in practice the real-time method yields the spectrum of H(τ ), not the one of H. Therefore the relevant question is: How much do the spectra of H(τ ) and H differ? In Appendix B we give a rigorous proof that the difference between the eigenvalues of H(τ ) and H vanishes as τ 2 . In other words the value of m (or t = mτ ) has no effect whatsoever on the accuracy with which the spectrum can be determined. Therefore the final conclusion is that the error in the eigenvalues vanishes as τ 2 /N where N is the number of data points used in the Fourier transform of Tr e −it H(τ ) .

IX. RESULTS
In all our calculations we take J = 1 and h = 0, except for the Ising model in a transverse field, where we take h = 0.75. The random numbers c n,p are generated such that the conditions Eqs. (A3) and (A4) are satisfied. We use two different techniques to generate these random numbers: 1. A uniform random number generator produces {f n,p } and {g n,p } with −1 ≤ f n,p , g n,p ≤ 1. We then normalize the vector (see Eq. (8)).
2. The c n,p 's are obtained from a two-variable (real and imaginary part) Gaussian random number generator and the resulting vector is normalized.
Both methods satisfy the basic requirements Eqs. (A3) and (A4) but because the first samples points out of a 2D-dimensional hypercube and subsequently projects the vector onto a sphere, the points are not distributed uniformly over the surface of the unit hypersphere. The second method is known to generate numbers which are distributed uniformly over the hyper-surface. Although the first method does not satisfy all the mathematical conditions that lead to the error (14), our numerical experiments with both generators give identical results, within statistical errors of course. Also, within the statistical errors, the results from the imaginary and real-time algorithm are the same. Therefore we only show some representative results as obtained from the real-time algorithm.
In Fig. 1 we show a typical result for the DOS D(ǫ) of the XY model, the Ising model in a transverse field and the mean-field model, all with L = 15 spins and using S = 20 samples. Because of the very high degeneracy we plotted the DOS for the mean-field model on a logarithmic scale.
In Fig. 2 we show the relative error δZ/Z based on Eq. (36) for the three models of various size, as obtained from the simulation (symbols). For these figures we used the imaginarytime algorithm, because then the statistical error can be related to e A directly (see Eq. (14) with A = exp(−βH)). The theoretical results (lines) for the error estimate, obtained by a direct exact numerical evaluation of (14) are shown too. We conclude that for all systems, lattice sizes and temperatures there is very good agreement between numerical experiment and theory.
Results for the energy E and specific heat C are presented in Fig. 3 (XY model), 4 (Ising model in a transverse field), and 5 (mean-field model). The solid lines represent the exact result for the case shown. Simulation data as obtained from S = 5 and S = 20 samples are represented by symbols, the estimates of the statistical error by error bars. We see that the data are in excellent agreement with the exact results and equally important, the estimate for the error captures the deviation from the exact result very well. We also see that in general the error decreases with the system size. Both the imaginary and real-time method seem to work very well, yielding accurate results for the energy and specific heat of quantum spin systems with modest amounts of computational effort.

X. CONCLUSIONS
The theoretical analysis presented in this paper gives a solid justification of the remarkable efficiency of the real-time equation-of-motion method for computing the distribution of all eigenvalues of very large matrices. The real-time method can be used whenever the more conventional, Lanczos-like, sparse-matrix techniques can be applied: Memory and CPU requirements for each iteration (time-step) are roughly the same (depending on the actual implementation) for both approaches.
We do not recommend using the real-time method if one is interested in the smallest (or largest) eigenvalue only. Then the Lanczos method is computationally more efficient because it needs less iterations (time-steps) than the real-time approach. However if one needs information about all eigenvalues and direct diagonalization is not possible (because of memory/CPU-time) there is as yet no alternative to the real-time method. The matrices used in this example (up to 32768 × 32768) are not representative in this respect: The real-time method has been used to compute the distribution of eigenvalues for matrices of dimension 16777216 × 16777216 [11].
Once the eigenvalue distribution is known the thermodynamic quantities directly follow. However if one is interested in the accurate determination of the temperature dependence of thermodynamic (and static correlation functions) properties but not in the eigenvalue distribution itself, the imaginary-time method is by far the most efficient method to compute these quantities. For instance the calculation of the thermodynamic properties for βJ = 0, . . . , 10 of a 15-site spin-1/2 system (i.e. implicitly solving the full 32768×32768 eigenvalue problem) takes 1410 seconds per sample on a Mobile Pentium III 500 Mhz system.
Finally we remark that although we used quantum-spin models to illustrate various aspects, there is nothing in the real or imaginary-time method that is specific to the models used. The only requirement for these methods to be useful in practice is that the matrix is sparse and (very) large. ACKNOWLEDGMENTS Support from the Dutch "Stichting Nationale Computer Faciliteiten (NCF)" and the Dutch "Stichting voor Fundamenteel Onderzoek der Materie (FOM)" is gratefully acknowledged.

APPENDIX A: EXPECTATION VALUE CALCULATION
In this appendix we calculate the expectation value of the error squared, as defined in Section II. By definition we have where p and p ′ label the realization of the random numbers c n,p ≡ f n,p + ig n,p . First we assume that different realizations p = p ′ are independent implying that Second we assume that the random numbers are drawn from a probability distribution that is an even function of each variable P (f 1,p , g 1,p , f 2,p , g 2,p , . . . , f k,p , g k,p , . . . , f D,p , g D,p ) = P (f 1,p , g 1,p , f 2,p , g 2,p , . . . , −f k,p , g k,p , . . . , f D,p , g D,p ) = P (f 1,p , g 1,p , f 2,p , g 2,p , . . . , f k,p , −g k,p , . . . , f D,p , g D,p ), and symmetric under interchange of any two variables P (f 1,p , g 1,p , . . . , f i,p , g i,p , . . . , f j,p , g j,p , . . . , f D,p , g D,p ) = P (f 1,p , g 1,p , . . . , f j,p , g i,p , . . . , f i,p , g j,p , . . . , f D,p , g D,p ) = P (f 1,p , g 1,p , . . . , g i,p , f i,p , . . . , f j,p , g j,p , . . . , f D,p , g D,p ), for all i, j, k = 1, . . . , D. This is most easily done by drawing individual numbers from the same even probability distribution i.e. . . , f D , g D ) we find that where in the last equality we omitted the subscripts of c m,p to indicate that the expectation value does not depend on m or p. An expectation value of a product of two c * 's and two c's can be written as Furthermore for m = n we have E(c * m,p c n,p c n,p c * m,p ) =E((f 2 m,p − 2if m,p g m,p − g 2 m,p )(f 2 n,p + 2if n,p g n,p − g 2 n,p )) =E(f 2 m,p f 2 n,p ) + 2iE(f 2 m,p f n,p g n,p ) − E(f 2 m,p g 2 n,p ) − 2iE(f m,p g m,p f 2 n,p ) + 4E(f m,p f n,p g m,p g n,p ) + 2iE(f m,p g m,p g 2 n,p ) − E(g 2 m,p f 2 n,p ) − 2iE(g 2 m,p f n,p g n,p ) + E(g 2 m,p g 2 n,p ) =E(f 2 m,p f 2 n,p ) − E(g 2 m,p f 2 n,p ) − E(f 2 m,p g 2 n,p ) + E(g 2 m,p g 2 n,p ) =0. (A8) By symmetry E(|c m,p | 2 |c n,p | 2 ) does not depend on m, n or p and the same holds for E(|c m,p | 4 ). The fact that the vector (c 1,p , . . . , c D,p ) is normalized yields the identities and where c and c ′ refer to two different complex random variables. Therefore we have and (δ m,n δ k,l (1 − δ m,k ) + δ m,k δ n,l (1 − δ m,n )) + δ p,p ′ δ m,l δ n,k δ m,n E(|c| 4 ). (A13) and the final result for the variance reads An expression for the fourth moment E(|c| 4 ) cannot be derived from general properties of the probability distribution or normalization of random vector. We can only make progress by specifying the former explicitly. As an example we take a probability distribution such that for each realization p the random numbers f n,p and g n,p are distributed uniformly over the surface of a 2D-dimensional sphere of radius 1. This probability distribution can be written as where we omitted the subscript p because it is irrelevant for what follows. The even moments of |c n | = (f 2 n + g 2 n ) 1/2 are defined by It is expedient to introduce an auxiliary integration variable X by We can perform the integration over X last and regard (A17) as the M-th moment of the variable X with respect to the probability distribution The calculation of P (X) amounts to computing integrals of the form Changing to spherical coordinates we have yielding The moments E(|c| 2M ) are given by and the values of interest to us are where the first two results provide some check on the procedure used. Substituting (A23) into (A14) yields

APPENDIX B: ERROR BOUNDS
Here we prove that the difference between the eigenvalues of the Hermitian matrix A + B and those obtained from the approximate time-evolution exp(zA/2) exp(zB) exp(zA/2) (z = −iτ, −τ ) is bounded by τ 2 . In the following we assume A and B are Hermitian matrices and take τ a real, non-negative number. We start with the imaginary-time case.
We define the difference R(τ ) by a well-known result [21]. We have [23] ||R(τ )|| ≤ Hence the bound in R(τ ) does not depend on the sign of τ so that we can write where C(τ ) is Hermitian. Clearly we have We already have an upperbound on R(τ ) and now want to use this knowledge to put an upperbound on the difference in eigenvalues of C(τ ) and A+B. In general, for two Hermitian matrices U and V with eigenvalues {u n } and {v n } respectively, both sets sorted in nondecreasing order, we have [2] |u n − v n | ≤ ||U − V ||, ∀ n.
Clearly (B14) proves that the differences in the eigenvalues of A + B and C(τ ) vanish as τ 2 .
We now consider the case of the real-time algorithm (z = −iτ ). For Hermitian matrices A and B the matrix exponentials are unitary matrices and hence their norm equals one. This simplifies the derivation of the upperbound on R(−iτ ). One finds [17] where ||A|| 2 E ≡ Tr A † A denotes the Euclidean norm of the matrix A [2]. In general the eigenvalues of a unitary matrix are complex valued and therefore the strategy adopted above to use the bound on R(τ ) to set a bound on the difference of the eigenvalues no longer works. Instead we invoke the Wielandt-Hoffman theorem [24]: If U and V are normal matrices with eigenvalues u i and v i respectively, then there exists a suitable rearrangement (a permutation ̺ of the numbers 1, . . . , n) of the eigenvalues so that Let U and V denote the exact and approximate real-time evolution operators respectively. The eigenvalues of A+B and C(τ ) are x n (0) and x n (τ ) respectively. All the x n 's and x n (τ )'s are real numbers. According to the Wielandt-Hoffman theorem where y j (τ ) = x ̺(j) (τ ), ̺ being the permutation such that inequality (B17) is satisfied. We see that Eq. (B17) only depends on (τ x j (0) mod 2π) and (τ y j (τ ) mod 2π), but so does the DOS (see Eq. (16)). Since the inequality (B17) and the DOS only depend on these "angles" modulo 2π, there is no loss of generality if we make the choice Rewriting the sum in (B17) we have As we have the restriction Eq. (B18) allows us to write implying In summary, we have shown that in the real-time case there exists a permutation of the approximate eigenvalues such that the difference with the exact ones vanishes as τ 2 .
Finally we note that both upperbounds (B22) and (B14) hold for arbitrary Hermitian matrices A and B and are therefore rather weak. Except for the fact that they provide a sound theoretical justification for the real and imaginary-time method, they are of little practical value.