Efficient numerical method for evaluating normal and anomalous time-domain equilibrium Green's functions in inhomogeneous systems

In this work we develop EPOCH: Equilibrium Propagator by Orthogonal polynomial CHain, a computationally efficient method to calculate the time-dependent equilibrium Green's functions, including the anomalous Green's functions of superconductors, to capture the time-evolution in large inhomogeneous systems. The EPOCH method generalizes the Chebyshev wave-packet propagation method from quantum chemistry and efficiently incorporates the Fermi-Dirac statistics that is needed for equilibrium quantum condensed matter systems. The computational cost of EPOCH scales only linearly in the system degrees of freedom, generating an extremely efficient algorithm also for very large systems. We demonstrate the power of the EPOCH method by calculating the time-evolution of an excitation near a superconductor-normal metal interface in two and three dimensions, capturing transmission as well as normal and Andreev reflections.


I. INTRODUCTION
The study of time-evolution in quantum mechanical systems provides fundamental insights into the system's underlying structure, allowing for the development of applications harnessing dynamical quantum effects, with prominent examples in quantum information processing [1][2][3][4]. Notably, even quantum systems without any external time-dependent drive usually present a highly complex dynamical behavior, exhibited in its time-dependent correlation functions without the need to go out of equilibrium. The dynamical behaviors range from the transmission and reflectance across interfaces to more general quantum transport properties, including also quantum electronics with coherent single electron excitations in solids [1,[5][6][7]. In superconductors, another remarkable time-dependent phenomenon occurs when electrons pair at unequal times, enabling odd-frequency pairing, present in for example superconductor-ferromagnetic and multiband systems [8][9][10][11][12][13][14]. Moreover, through the fluctuation-dissipation theorem, equilibrium correlations even predict the response to time-dependent external perturbations [15][16][17].
At the same time, many systems with highly nontrivial quantum dynamics also lack translation invariance, such as those dominated by the presence of junctions, interfaces, edges, or disorder. Therefore, these systems also see a dramatic growth in the number of degrees of freedom that must be treated simultaneously. This fact makes the use of analytical methods, as well as brute force eigenvalue diagonalization, completely unrealistic, even for systems with a noninteracting quasiparticle description and in equilibrium. Faced by this very challenge, we seek a method that efficiently and accurately compute time-domain properties with a low computational cost even as the number of degrees of freedom * Corresponding author: jorge.cayao@physics.uu.se proliferate.
To extract interesting dynamical properties, the correct objects are usually the time-domain Green's functions, encoding the physical observables of any system. Fortunately, for any generic non-interacting system, all higher order Green's functions reduce to the computation of the single particle Green's functions by virtue of Wick's theorem [17]. The Green's functions are also by themselves inherently interesting objects, representing the direct physical probability amplitude of finding a particle at a point in space-time (x 2 , t 2 ), after an earlier insertion at (x 1 , t 1 ). In superconductors, where particle number is not conserved, the inserted particle may later be found as a hole, and the amplitude for this conversion process is properly handled by the anomalous Green's functions. Here, as an example, the unequal time anomalous Green's functions may display exotic odd-frequency dynamical pairing [12,13], and are therefore also of large interest.
In terms of evaluating any physical observables in the time-domain, we first note that, for pure quantum states, the time-domain naturally lends itself to a step-by-step propagation in time without the need of any costly diagonalization. Simply because the Hamiltonian H is precisely the time-evolution generator in any quantum system. Therefore, the time-dependence of a pure quantum state, i.e. single wave-packet propagation, is directly achievable with only matrix-vector multiplication; for instance, by Taylor expanding the unitary time evolution operator U (t) = e −iHt , by naive finite differences of the Schrödinger differential equation, or by somewhat more involved algorithms [18].
The single wave-packet propagation methods are vastly improved, both in terms of accuracy and efficiency (exponential over power law convergence), by the Chebyshev method of time-evolution [19], a method that is already firmly established in quantum chemistry and molecular dynamics [20][21][22][23]. The basis of the Chebyshev time-evolution method is the expansion of the evo-lution operator in terms of Bessel functions J n (t) and the Chebyshev polynomials T n (x) as: U (t) = ∞ n=0 (2 − δ n0 )(−i) n J n (t)T n (H) [24], whereH = H/ H andt = H t to fit the standard domain of T n (x). The resulting series has the advantageous property that it converges rapidly, because fort n the Bessel J n (t) is exponentially suppressed. Therefore, before a given finite timet, it is only the first n terms that contribute significantly. A second underpinning of the method's efficiency is that at each expansion order, the next order is obtained via a recursive relationship connecting the T n (H) polynomials of neighboring degree. This limits the computational cost to just one additional (sparse) multiplication by the HamiltonianH per expansion order. As a result, the Chebyshev wave-packet propagation method achieves both accuracy and stability with a low computational cost.
In a broader context, the efficiency of the Chebyshev wave-packet propagation method relates directly to the advantageous properties of the underlying Chebyshev polynomials T n (x), which, like all orthogonal polynomials, form a complete basis set. It is therefore not surprising that orthogonal polynomials have not just been found to be useful for wave-function propagation, but also have a long history as a versatile technique in approximation theory [25] and also applied to electronic structure methods, where they have developed into a prominent linearscaling method known as the kernel polynomial method (KPM) [26,27]. This method has notably been used to compute spectral functions and expectation values in very large systems [28,29] and underpin several new software packages [30][31][32][33]. Despite this, it was only recently that KPM was extended to also tackle superconducting systems [34,35], but still not in the time-domain, which is necessary when studying, for instance, time evolution.
While the Chebyshev wave-packet propagation method has already been highly successful to evaluate observables in the time-domain, a fundamental limitation for condensed matter system is that it is only applicable to pure quantum states, i.e. single electrons, and not to many electron states, i.e. mixed states. Thus it cannot evaluate, for instance, thermodynamic expectation values, which are central concepts in condensed matter. Therefore, in order to achieve the same efficiency and accuracy as the Chebyshev wave-packet propagation method also in condensed matter systems, we need to incorporate Fermi-Dirac statistics alongside the time-evolution.
In this work, we develop the Equilibrium Propagator by Orthogonal polynomial CHain, or EPOCH, method, to efficiently calculate the time-domain equilibrium Green's functions for both normal and anomalous correlations in generic quantum condensed matter systems. Generalizing the Chebyshev wave-packet propagation method from quantum chemistry, the EPOCH method efficiently incorporates the Fermi-Dirac statistics that is needed for the equilibrium physics of quantum condensed matter systems.
To derive the EPOCH method, we first strategically depart from the Chebyshev method and instead use the Legendre polynomials, since they in contrast allow us to derive analytical expressions for key quantities in the time-domain. This changes our reference starting point, from the expansion of the evolution operator U (t) in terms of the Chebyshev polynomials T n (x) to an analogous expansion in terms of the Legendre polynomials P n (x) (see [22], Eq. (7.12)) where j n (t) are the spherical Bessel functions [36]. Next, in order to move beyond the wave-packet propagation and towards a Green's function formalism, we in Section II first review the standard Bogoliubov-de Gennes (BdG) formalism for a generic time-independent Hamiltonian H with an emphasis on the time-domain [17,37], since this formalism is applicable to systems both with and without superconductivity. Thereafter, we show how, at any inverse temperature β = 1/T , the thermal equilibrium Green's functions (the lesser G < (t), greater G > (t), and anomalous Green's functions F(t)) are all found as the matrix elements of a quantity which we refer to as the Equilibrium Propagator (EP) given by The key component of the EPOCH method then boils down to a general series expansion of the EP, which we derive in Section III. When specifically applied to the Legendre polynomials P n (x) in Section IV, we arrive at, whereH = H/ H andt = H t to fit the standard domain of P n (x). This expansion can be thought of as a generalization to many-body fermionic systems of the Chebyshev wave-packet propagation method already widespread in quantum chemistry [19][20][21][22][23] or, likewise, as the series expansion for the time-evolution operator in Eq. (1). Within this generalization, the coefficients as shown in Section IV, the solutions to a closed recurrence relationship with an inhomogeneous source term, and we also show how to calculate them with a numerically stable and efficient method. The computational cost in the EPOCH method is kept extremely low as the method scales linearly with the degrees of freedom (system size) and with the longest evolved time, i.e. running the time evolution for twice as long will only take twice as long to compute, all else being equal. This is because each expansion order of Eq. (2) requires only one additional multiplication by H (on very general grounds is a sparse matrix), which follows from the three-term recursion relationship of the Legendre polynomials. Another strength of the EPOCH method is that any evolved time (or arbitrary set of times) at all temperatures is directly and easily computable. Moreover, the formalism offers a detailed account of any error introduced by the inevitably finite truncation of the EP series expansion. In Section V we establish the accuracy of the expansion, and also show how Gibbs phenomenon, a possible error source due to the sharp Fermi surface in systems without an energy gap, is effectively avoided by going to small but finite temperatures. For easy use and summary, we end the formal development of the EPOCH method in Section VI by providing a step-by-step outline of the numerical implementation together with the error bounds.
As a demonstration of the capabilities of the EPOCH method, we first use it in Section VII to calculate the particle propagation both across and around a junction between a superconductor and a normal metal in both two and three dimensions. Despite the very large system sizes involved, especially for the three-dimensional problem, we capture all the transmission and reflection amplitudes, including the Andreev processes. While these examples are simple, they nonetheless clearly illustrate the capabilities of the EPOCH method, and how it enables efficient computation of the complex quantum dynamics even on a standard laptop. As a second demonstration, we show in Section VIII how the dynamical correlations between observables and the general time-dependent linear response to external probes within the Kubo formalism can be computed directly within the EPOCH method. In this application, a clear advantage of EPOCH is that it automatically gives all time and temperature dependence, without the need to recompute the quantum propagation. EPOCH therefore enables calculation of the response to pulses probes directly in the time-domain, predicting observables and material properties measured by e.g. scattering, polarizability, and transport [38].

II. TIME-DEPENDENT BOGOLIUBOV-DE GENNES FORMALISM
To provide a simple starting point and ensure that our approach is applicable to condensed matter systems with or without superconductivity, we review the standard BdG formalism. We consider a fully general quadratic time-independent Hamiltonian H and separate it into two partsĤ =Ĥ 0 +∆ where the first partĤ 0 conserves particle number, thus capturing the normal part, while the second part∆ contains all terms that break the gauge invariance, notably the terms of a superconducting condensate, if present. Such a Hamiltonian models any device, with both a spatially varying normal state or a superconducting order in any part of the system.
Adopting the block Nambu-spinor X = c c † T of dimension 2N for the N degrees of freedom in the system (lattice sites, spin, and orbital degrees of freedom), the Hamiltonian takes, up to a constant shift, the BdG bilinear form [17,37] Using the matrix block structure, the eigenvectors of the BdG Hamiltonian can be written as, H u v * T = E u v * T , with amplitudes u and v.
The particle-hole symmetry inherent to the BdG Hamiltonian implies that each eigenvector also has a symmetry companion state of the opposite energy, The accompanying canonical transformation X = U Y defines the eigenstates Y = γ γ † T having definite energies, stored on the diagonals of E in E = diag (E, −E). Notably, the eigenstates satisfy the fermionic anticommutation rules: {γ i , γ j } = 0 and γ i , γ † j = δ ij . In the Heisenberg picture they are also the eigenmodes of the time evolution operator i∂ t γ s (t) = γ s (t) ,Ĥ = E s γ s (t) and therefore we find their time-dependence simply from γ s (t) = e −iEst γ s , for each eigenstate s.
While each eigenstate has a definite time-dependence, the state of a condensed matter system is more complicated, consisting of mixed state that is not reducible to a single eigenstate. The computation of observables therefore requires not only that the eigenvalues themselves are found, but the summation over the correct overlap functions of the many involved eigenstates. An alternative is offered by the time-dependent Green's function, as they relate directly to physical observables.
Directly corresponding to either transition or pairing amplitudes, the two-point single-particle Green's functions are the inhomogeneous and homogenous solutions to the equation [17,39] Since this equation is linear, any linear combination of the solutions is also a solution, allowing for many possible definitions, satisfying different boundary conditions. We choose to work with the lesser G < ij (t 1 , t 2 ) and greater G > ij (t 1 , t 2 ) Green's functions as they relate directly to the the local electronic density, and as such naturally includes the quantum statistics via the Fermi-Dirac function. The same holds true for the anomalous Green's function capturing the pairing amplitudes F ij (t 1 , t 2 ). These Green's functions are homogenous so-lutions, following from their definitions [17,39], and through linear combinations including appropriate factors of i and the Heaviside step function θ(t) they also directly give other conventionally defined Green's functions, including the retarded, advanced, and timeordered versions [17,39], Thus, from the lesser, greater, and the anomalous Green's functions, which we will calculate, we can easily find the other Green's functions.
Since the eigenstates with their definite timedependence completely specify the system within the BdG formalism, they also determine the timedependence of the Green's functions. This is because, between any two generic times t 1 and t 2 , a transformation to the eigenstates via the transformation U gives the explicit time-dependence through the matrix elements v ij and u ij as where F β (E s ) = γ † s γ s is the Fermi-Dirac function at the inverse temperature β = 1/T and the summation runs over all eigenstates indexed by s.
Having established the explicit expressions of the timedependent Green's functions in Eq. (4), we next show that they are given by the matrix elements of the EP L β (H, t) given by Eq. (2). In fact, just like the Green's functions, the EP viewed as a matrix equation is also a homogeneous solution to [i∂ t − H]L β (H, t) = 0, while, in addition, also the unique solution that for equal times reduces to the Fermi-Dirac projection F β (H) and capturing the equilibrium statistics. This fundamentally es-tablishes a direct connection between the EP and the Green's functions given by Eqs. (4). In fact, the explicit eigenbasis forms of the Green's functions, also show that the Green's functions are the matrix elements of the EP. Because in the same eigenbasis, the Hamiltonian has a diagonal form H = U EU † and the matrix function L β (H, t) = e −iHt F β (H) is found by simply applying the same function to the real eigenvalues on the diagonal entries of E: Multiplying out the matrix blocks shows that each block corresponds to one of the Green's functions when compared to the explicit expressions of Eq. (6), Thus, the upper (lower) diagonal block in Eq. (7) is the lesser (greater) Green's functions and the upper off-diagonal block correspond to the anomalous Green's functions with its conjugate in the lower off-diagonal block, as defined in Eq. (4). We therefore conclude that, by calculating the EP, we have full access to all two-point Green's functions, or equivalently, to a singleparticle/hole excitation over time in a system with a time-independent Hamiltonian. We, therefore, focus on efficiently calculating the EP in the subsequent sections.

III. ORTHOGONAL POLYNOMIAL EXPANSION OF THE EQUILIBRIUM PROPAGATOR
It is clear from Eq. (7) of the previous section that the time-dependent Green's functions are numerically attainable from diagonalization of the Hamiltonian. Diagonalization however incurs a O(N 3 ) computational complexity cost for a system with N degrees of freedom, thus prohibiting its use for large inhomogeneous systems. Fortunately, on very general physical grounds, the full Hamiltonian matrix of any typical quantum system is highly sparse [26] due to short-ranged hopping amplitudes and interactions. This sparseness opens up for linear-scaling iterative methods. In particular, orthogonal polynomial expansion readily produces O(N ) methods because of the three-term recursion relationship connecting successive polynomials, while at the same time being numerically stable [26,27,40]. Thus, spurred by the impasse at diagonalization, we proceed with the use of orthogonal polynomials for computing the EP matrix elements.
Generally, a set of orthogonal polynomials {φ n } are defined from an inner product with a weight function w(x), on an interval I where they form a complete set of functions with the completeness relationship thus allowing for a generalized Fourier series expansion of generic functions on the interval. Because the BdG Hamiltonian matrix is Hermitian and assumed to have a bounded spectrum, we can employ a set of orthogonal polynomials to expand matrix functions of the BdG matrix, too. Because H is in principle diagonalizable with real eigenvalues, the action of the matrix function defining the EP, L β (H, t) = e −iHt F β (H), would in a diagonal representation be equivalent to applying the associated function L β (E, t) = e −iEt F β (E) to the real eigenvalues. Thus, to expand the EP, we can work directly with the real-domain function L β (E, t).
Notably, inserting the above polynomials completeness relationship and integrating over the kernel K of Eq. (8) produces the series expansion, The immediate benefit of this step is that all of the dependence on time and temperature are conveniently isolated to the mode transients defined as The transience of l n β (t) for any given n is a direct result of the Riemann-Lebesgue lemma, guaranteeing that l n β (t) → 0 in the limit t → ∞ [41]. As stated, the expansion in Eq. (9) is directly carried over to the corresponding well-defined EP matrix function as meaning that the convenient separation of variables also extends to this orthogonal polynomial expansion of the EP.
Independently of what set of polynomials are used, the expansion in Eq. (11) is, as shown, also a a series expansion of the Green's functions, as these are the matrix elements of the EP. The choice of the polynomials {φ n }, does however affect the definitions of the functions l n β (t) that depend on the set used.
To clearly show the equality between the (anomalous) Green's functions and the matrix elements of L, we introduce new basis vectors in the 2N dimensional vector space for the particle and hole part of the BdG Hamiltonian: [e i ] j = δ ij and [h i ] j = δ (i+N )j . For with these definitions, the explicit equality is: It is this series mode expansion, which explicitly separates time and temperature effects, that is the foundation for the EPOCH method for computing time-dependent Green's functions.

IV. LEGENDRE POLYNOMIAL EXPANSION OF THE EQUILIBRIUM PROPAGATOR
To use the mode expansion of Eq. (12) we need to evaluate both the elements of the matrix polynomial φ n (H) and the time-dependent mode transients l n β (t) given by Eq. (10). However, attempting the direct evaluation of either one of these is both numerically unstable and inefficient. For instance, we have not found any simple closedform analytical expressions for l n β (t). Moreover, the integrand that defines l n β (t) in Eq. (10) becomes for large n or t highly oscillatory from the energy phase factor and the high degree polynomial, thereby hindering direct numerical integration beyond the first few terms. The key to surmounting these difficulties is to use the recursive relationship between subsequent orthogonal polynomials to derive recurrence relationships for both the matrix elements φ n (H) and the mode transients l n β (t). In this way, both of these quantities are readily computed iteratively order by order, which achieves both numerical stability and efficiency in one sweep. This is a key ingredient of the EPOCH method.
For definite expressions suitable for numerical treatment, we first have to choose an appropriate set of orthogonal polynomials. We choose the Legendre polynomials P n (x) defined on the interval I = [−1, 1] having P n 2 = 2/(2n + 1). Notably, the weight function w(x) = 1 associated with the Legendre polynomials gives equal weight and therefore also equal representational accuracy to the whole interval and spectrum. This is in contrast to the Chebychev polynomials, which add an unphysical weighting issue, see Appendix A. Moreover, the straightforwardness of the Legendre polynomial weight function also streamlines all ensuing calculations thereby allowing us to derive key analytical expressions, as we show below.
Since any finite dimensional Hamiltonian has a bounded spectrum, we can without loss of generality assume that the spectrum is contained within the interval I, since the Hamiltonian can always be rescaled by a constant: H →H = H/λ. [42] For a bandwidth W , the spectrum ofH will be entirely contained in In what follows and for ease of presentation, we therefore assume that such a scaling has been performed and work directly with the dimensionless rescaled Hamilto-nianH, which has its spectrum entirely contained in I. Similarly, all quantities with energy dimensions are also subsequently rescaled, including the inverse temperature β →β = λβ and time t →t = λt, which all then become dimensionless quantities. Proceeding with the Legendre polynomials, below we calculate the matrix elements of φ n (H) = P n (H) and the mode transients l ñ β (t) in order to generate the series mode expansion in Eq. (12).

A. Matrix elements
The matrix elements e † i P n (H)e j = (P n (H)e i ) † e j and e † j P n (H)h j = (P n (H)e i ) † h j that enter the mode expansions of the normal and anomalous Green's functions in Eq. (12) are efficiently computed by using the threeterm recurrence relationship for the Legendre polynomials (n + 1) P n+1 (x) = (2n + 1)xP n (x) − nP n−1 (x). For this relationship connects sequential orders of the lefthand vector e n i = P n (H)e i via the recursion relation As a result, the matrix elements for all orders can be computed iteratively as the inner products e † i P n (H)e j = (e n i ) † e j and e † i P n (H)h j = (e n i ) † h j . Because these inner products are computationally cheap, the only significant computational cost associated with the matrix elements is, therefore, the one sparse matrix multiplicationHe n i for each use of the recursion. In addition, once the orders of a given initial state e n i have all been computed, the inner product with any end state, particle or hole, carries essentially no additional cost, making all Green's functions of that initial state available. Consequently, the complexity cost for the whole mode expansion is just O(N ), making it very efficient.

B. Mode transients
The second ingredient needed for the mode expansion of Eq. (12), in addition to the matrix elements, are the mode transients l ñ β (t). Their defining equation Eq. (10) can be seen as the Fourier transform of the product between the Legendre polynomials and Fermi-Dirac function on the interval I, since for the Legendre polynomials the weight function is identically one. Without the Fermi-Dirac function, this Fourier transform has a simple analytical close form, because the Fourier coefficients of the Legendre polynomials are directly related to the spherical Bessel functions j n (t) through [43], Eq. 7.243.5), which is also the basis for constructing Eq. (1).
The transform of the mode transient including the Fermi-Dirac function, is however not so simple. Yet, parts of the previous result can be salvaged, and the the above connection to the Bessel functions is maintained even for the non-trivial mode transients l ñ β (t). For because of the odd/even symmetry of the Legendre polynomials, P n (−Ẽ) = (−1) n P n (Ẽ) and the partition of Consequently, for all inverse temperaturesβ Thus for all n either the real or imaginary part of the complex function l ñ β (t) are already given by the spherical Bessel functions j n (t).
To find the remaining complementary parts of l ñ β (t) that is not given by the Bessel functions, we derive a new recursion relationship, by once more exploiting the recursive properties of the Legendre polynomials. Using that (2n ) and partial integration, we find, We recognize Eq. (15) as intimately connected to the well-known recurrence relationship for spherical (+) or modified spherical (−) Bessel functions z n+1 (x) ± z n−1 (x) = (2n + 1)z n (x)/x (see Ref. [43], Eq. 8.4718.1), but with an inhomogeneous source term (16) To further isolate the remaining part of l ñ β (t), we define each mode transient to be the sum of a unitary part that is simply given by the spherical Bessel function j n and one as of yet unknown projective part f ñ β (t): By this definition all of the temperature dependence is completely relegated to the projective part f ñ β (t). From Eq. (15), it follows, as we find, that the projective part f ñ β (t) satisfies the inhomogeneous spherical Bessel recurrence relation: Notably, by exploiting the dependence of the integrand in Eq. (16) with respect toẼ and n, we find that the inhomogeneous part in Eq. (17), i n−1 Sβ ,n (t)/t, is real for all n. Therefore all f ñ β (t) are also real. To summarize, the results of this section show that when the Legendre polynomials P n (x) are used as the orthogonal polynomial modes in the expansion of the EP in Eq. (11), the mode transients l ñ β each decompose into a unitary and a projective transient. First, the unitary transients are given by the spherical Bessel functions j n (t). In standard wave-packet propagation it is only this unitary part that contributes, giving the expansion presented already in Eq. (1). Second, the projective transients f ñ β (t) encode the equilibrium statistics and, therefore, account for all the temperature dependence of the Green's functions. Like the unitary transients, the projective transients are real functions satisfying the spherical Bessel recursion relationship, but for f ñ β (t) the recursion is inhomogeneous because of the source term Sβ ,n (t). With the Bessel functions j n (t) well-known, we only have left to solve for the f ñ β (t) to complete our task of computing the EP and thereby all two-point Green's functions. This task requires us to evaluate the source term Sβ ,n (t) in Eq. (16), which we do next.

C. Inhomogeneous source term
At first it might not seem that the integral defining the source term Sβ ,n (t) in Eq. (16) is any easier to evaluate than the integral originally defining the general mode transient in Eq. (10). But in analogy to the Sommerfeld expansion [44], the source term Sβ ,n (t) in Eq. (16) is readily computed as a low-temperature expansion in 1/β, by simply exploiting that for low temperatures F β (Ẽ) is a sharply peaked function that is exponentially localized toẼ = 0. As a consequence, the sole contribution to the integral defining Sβ ,n (t) is only coming from a small neighborhood around zero energy of width 1/β. The energy scale associated with this width should be compared to the overall bandwidth W of the system, and we can therefore extend the limits of integration to the whole real line, since for any realistic condensed matter system βW 1. Further, because integration is a linear operation, we can additively distribute the integration over the Legendre polynomials appearing in the integral defining Sβ ,n (t) in Eq. (16). Sβ ,n (t) is therefore always a linear combination of the dynamical Fermi moments Iβ ,k (t) = 1 −1Ẽ k −F β (Ẽ) e −iẼt dẼ. This allows us to first consider the integral of Iβ ,k (t) and afterward take the appropriate linear combination of the result. By introducing y =βẼ, we find where we have first rewritten the integral as a derivative over the integral; next, because the integrand is exponentially suppressed in both directions, we have extended the integration to the whole real line allowing us then to perform the contour integration in the last step. In restored units, the contributions from outside the original interval are exponentially suppressed in the ratio of the bandwidth W to the temperature by a factor of e −βW .
Thus, the approximation introduced by extending the integral is inconsequential for any realistic condensed matter system where βW 1. For very large moment ordinals k, however, the integrand eventually grows outside even this interval and the assumption in Eq. (18) breaks down for very large k eβW . Still, using the expression of Eq. (18) for the Fermi moments produces an asymptotic expansion that is very accurate for all moments with k βW .
The source term Sβ ,n (t) is now a linear combination of the dynamical Fermi moments Iβ ,k (t) of Eq. (18), where the coefficients are given by the explicit representation of the Legendre polynomials, P n (x) = 2 n n k=0 n k n+k−1 2 n x k . Thus, we finally arrive at a lowtemperature asymptotic expansion of the source term, where we use the gamma function Γ. In the zero temperature limit,β → ∞, it is clear that all the dynamical Fermi moments become time-independent and only the first term is non-zero, Iβ ,k (t) → δ k0 . In this case therefore S ∞,n (t) is also time-independent and nonzero only for odd n, and reduced simply to: S ∞,n = (2 n (2n + 1)Γ n 2 )/(Γ − n 2 Γ (2 + n)). To summarize, because the derivative of the Fermi function F β (Ẽ) is sharply peaked around zero for all temperatures that are small compared to the bandwidth W , we can compute the inhomogeneous source terms Sβ ,n (t) as an asymptotic series expansion in the dynamical Fermi moments Iβ ,k (t), which is rapidly converging. As a consequence, we can either use simply the zero temperature source term S ∞,n or the first terms of the low temperature expansion in Eq. (19) to always find a good numerical approximation Sβ ,n (t) for realistic temperatures.

D. Calculating mode transients
With the inhomogeneous source term Sβ ,n (t) determined in the previous subsection, we find the total mode transients l ñ β (t) by solving the recurrence relation Eq. (17) for the projective transients f ñ β (t). However, when attempting to numerically solve this recurrence relationship, care must be taken.
In fact, if one attempts to solve it by simple forward propagation, calculating the next term from the previous two, then rounding errors inevitably and disastrously accumulate with each step. This is a well-known danger of second order recurrence relationships, including that of the Bessel functions [45]. For our purposes, the instability of Eq. (17) is most readily understood starting from its underlying homogeneous Bessel recurrence relationship with its two linearly independent solutions, commonly called J n and Y n . Here, the massive domination of one solution over the other (|Y n (x)/J n (x)| ∼ 2(n!) 2 /(x/2) 2n → ∞, as n goes to infinity) implies that even if just a small rounding error falls outside the subspace of J n during numerical stepping procedure, then this spillover will quickly grow and overtake the calculated solution such that any resemblance to the desired J n is disastrously lost.
Fortunately, there already exists a stable algorithm for solving linear recurrence relationships, including inhomogeneous relationships. Central to this algorithm is the reformulation of the recurrence relationship as an auxiliary linear boundary value problem; the left boundary is set by the true initial-value of f 0 , but the right boundary, at a sufficiently large M , is instead set to an artificially imposed value of f M = 0. Solving this boundary value problem has been shown to converge to the minimal solution, as the right bound is moved further away with larger M [46,47]. The recurrence equation of Eq. (17) for the projective transients f ñ β (t) is therefore easily solved as a boundary value problem, given that we provide the initial condition f 0 To illustrate the mode transients we plot in Fig. 1 the real and imaginary parts of l ñ β=∞ (t) for n = 11, where the real part is the projective transient f ñ β (t) generated by Eq. (17), while the imaginary part is the unitary transient, the spherical Bessel function j n (t). As seen, the projective transient resembles a phase shifted version of its companion unitary transient, characterized by an oscillating tail, decaying with time ∼ 1/ √t for larget. Similarly, the envelopes of both f ñ β (t) and j n (t) are of maximal height in the neighborhood oft ∼ n. Thus for a giveñ t, it is precisely the orders with n ∼t that mainly contributes to the time dependence in the series expansion of the EP.
Generally therefore, it is the envelopes of the two transient functions that dictate which terms of the expansion make a significant contribution to the mode expansion and in turn also sets the accuracy of the Green's functions when computed from the truncated mode expansion. In this regard, there is a distinction between the unitary and projective transients. For while j n (t) is exponentially suppressed fort n, the f ñ β (0) can take on finite initial values, meaning that even for short evolution times a finite number of terms, to be specified in Sec. V, are needed to capture the projective part.

E. Static limit
While our main focus in this work is on the dynamical properties and the time-domain Green's functions, the static limitt → 0 is also of interest for two reasons. First, whent → 0, the EP reduces to the density matrix Lβ(H, 0) = Fβ(H) and its elements are thet = 0 Green's functions, measuring the densities of the system: the local electron density and equal-time pair amplitudes, respectively. To have direct access to these with in the same formalism is of course of great utility. Secondly, as seen in the previous subsection, the finite initial values of f ñ β (0) inform the rate of convergence of the mode expansion; because, unlike the case with wave-packet propagation, even att = 0 the Fermi-Dirac statistics is encoded in f ñ β (0). Conveniently, in thet → 0 limit the mode transient recursion Eq. (15) decouples and consequently the mode transients are proportional to the source terms: l ñ β (0) = i n−1 f ñ β (0) = Sβ ,n (0)/(2n + 1) for all n ≥ 1, whence j n (0) = 0. In addition, the static limit of the dynamical Fermi moments of Eq. (18) where B k are the Bernoulli numbers (follows from the Taylor expansion x/ sinh(x) = 1 + ∞ n=1 2(1 − 2 2n−1 )B 2n x 2n /(2n)!). Inserting these static Fermi moments in Eq. (19) to get the source terms Sβ ,n (0), the static limit of the mode expansion of EP in Eq. (11) gives a polynomial expansion of the density matrix:

Lβ(H, 0) = Fβ(H)
where l 0 β (0) = 1 is the origin of the n = 0 term. In the zero temperature limit, T = 0, only the first k = 0 term remains Both Eqs. (20) and (21) are particularly useful expansions of the density matrix (Fermi operator) even at finite temperatures compared to the pioneering expansions of Ref. [48,49]. Thus from Eqs. (20) or (21) the density expectation values of the system are given by Lβ(H, 0)]e j and c i c j = e † i [Lβ(H, 0)]h j . Because the assumption underlying the low-temperature expansion of the Fermi moments in Eq. (18) eventually breaks down for very large moment ordinals k, the absolute values of the static moments Iβ ,k (0) pass through a minimum at k ∼ βW before eventually growing again. In turn, the finite temperature series of Eq. (20) is not formally convergent to all orders. Still, because the minimum of Iβ ,k (0) is exponentially suppressed in the ratio of the bandwidth to the temperature as e −βW , the series of Eq. (20) is an exponentially accurate asymptotic expansion for the density matrix. In the zero temperature limit, the asymptotic limit of the coefficients in Eq. (21) decay as n −3/2 , demonstrating the norm convergence of this series.

V. TRUNCATION OF EQUILIBRIUM PROPAGATOR EXPANSION
In the previous sections we have derived an exact expansion of the time-dependent two-point Green's functions using Legendre polynomials. While this series is, in principle, exact, in practice the expansion of the EP Lβ(Ẽ,t) in Eq. (9) has to be truncated at a finite number of terms in order to produce a numerically useful algorithm. Because the unitary part j n (t) only contributes significantly aftert n, the truncation of this series to order M yields a remarkably accurate approximation of the time-dependence for allt M . We illustrate this in Fig. 2(a,b) where we plot the numerical error L ∞ (Ẽ,t) − L M ∞ (Ẽ,t) for M = 60 and M = 1000, respectively. The only exception is a narrow energy window δẼ ∼ 1/M closest to the Fermi surface atẼ = 0, where the representation of the step-like Fermi-Dirac function in the projective transients f ñ β (t) starts to deviate because of Gibbs phenomenon, see below in Sec. V A. This is similar to the residual norm of the truncated static mode expansion, which also vanishes as L ∞ (Ẽ, 0) − L M ∞ (Ẽ, 0) ∼ 1/M , because the coefficients of Eq. (21) vanish as n −3/2 . However, this error source at low energies is inconsequential for any gapped system, such as insulators or superconductors with an energy gap ∆. Simply put, since there are no contributing states in the misrepresented region around the Fermi level, there can be no numerical distortion of the Green's functions. For example, if for a superconductor we assume a temperature below the transition temperature, we are guaranteed that the time-dependent Green's functions computed using the truncated mode expansion of Eq. (12) have excellent fidelity for allt M , as long as ∆/W 1/M , where W is the bandwidth.

A. Gibbs phenomenon
The numerical errors found at the Fermi level and at low temperatures in Fig. (2) can be understood as a man- ifestation of the Gibbs phenomenon; general ringing or oscillations found when an orthogonal series expansion is used to capture a step function. In short, therefore, the Gibbs phenomenon appears near the sharp Fermi surface in the finite series representation. As such, it is however only relevant for systems with a finite density of states near the Fermi energy. Thus, the Gibbs phenomenon is not relevant for an insulating or gapped superconducting system, but it could matter for a metallic systems.
The Gibbs phenomenon comes about because truncating Lβ(Ẽ,t) in Eq. (9) at the finite order M , amounts to replacing the underlying true integration kernel of Eq. (8) with the truncated Dirichlet kernel, K D M (x, x ) = M n=0 P n (x) P n (x ) / P n 2 . This replacement is a controlled approximation that does narrow in on the Dirac delta function with a diminishing width of δẼ ∼ 1/M [50], as illustrated by the black curve in Fig. 3(a), where we plot K D M (0,Ẽ) for M = 60. Nonetheless, it is also apparent that this kernel is subject to Gibbs oscillations. A direct consequence of these oscillations are the overshooting ripples in the finite series Dirichlet kernel approximation of the discontinuous T = 0 Fermi-Dirac function seen in Fig. 3(b).
Echoing standard Fourier theory [41], the standard approach for combating the Gibbs phenomenon is to replace the naive Dirichlet kernel with a suitable summability kernel [27].
For instance, the Jackson kernel is a common choice, which is an everywhere positive optimal kernel with a minimal squared peak width, as displayed by the red curve in Fig. 3(a). It is defined by reweighing the terms of the series curve shows how adopting the Jackson kernel results in an oscillation-free approximation of the T = 0 Fermi-Dirac distribution function. While this is a satisfactory solution for the static representation of the Fermi-Dirac function, the diminished weight given to the high-order terms by the Jackson kernel (and also by other traditional summability kernels such as the Fejér Kernel or the Lorentz kernel [27]) simultaneously means that the convergence of the representation is severely damaged when extended to time-evolution. For when evolved to a timet, then it is the terms of order n ∼t that are the largest and have the most significant contribution to the representation. Thus, the diminished weight given to the higher order terms by the summability kernels, causes the approximation of Lβ(Ẽ,t) to deviate even after just a short time-evolution, as we display in Fig. 3(c) along with this effect. In contrast, the Dirichlet kernel is, also in Fig. 3(c), seen to accurately represent Lβ(Ẽ,t) for allt n, while the Jackson kernel fails even for small t. Thus, the standard approach for circumventing the Gibbs phenomenon is not useful when going beyond the static Fermi-Dirac function to the time-evolution.
Although we find that the widely accepted technique of kernel resummation is not appropriate for time-evolution, we nonetheless find a working solution to the Gibbs phenomenon by transitioning to finite temperatures. At nonzero temperature, the step discontinuity of the Fermi-Dirac distribution is replaced by a smooth transition of width δẼ ∼β. The distribution therefore becomes amenable to an accurate representation by a polynomial approximation of order M ∼β, as in Eq. (20). We illustrate this in Fig. 3(b) where theβ = 60 distribution is plotted with M = 60 moments.
Importantly, by constructing the finite temperature representation of the EP by solving the recurrence relationships of Eq. (17) for the projective mode transients with the first few terms of the low-temperature expansion Eq. (19), we avoid the problems arising from Gibbs phenomenon, and the resulting representation of the EP has all the accuracy of its zero-temperature version at all timest n away from the Fermi surface. But in addition, the representation is now also accurate around the step of the Fermi-Dirac distribution. Meaning that, Lβ(Ẽ,t) is accurate for allt n/2, including the energies near the Fermi surface, as is illustrated in Fig. 2(c). As a consequence, even quantitative calculations in systems with a finite density of states near the Fermi energy are possible at any finite temperature.

VI. SUMMARY OF EPOCH METHOD
In summary, we have shown that time-domain Green's functions of a time-independent Hamiltonian H are efficiently and accurately computed as the matrix elements of the EP by the orthogonal polynomial expansion Eq. (11) using Legendre polynomials. We call this method EPOCH, standing for Equilibrium Propagator by Orthogonal polynomial CHain. For additional clarity, we here summarize the steps of the EPOCH method:  (17) for the projective transients f n λβ (λt) as a boundary value problem with the method of Ref. [46,47]. For the source term S λβ,n use either the zero temperature source terms S ∞,n or the first terms of the lowtemperature expansion Eq. (19).
[(e n i ) † e j ]l n λβ (λt) , The mode expansion in Eq. (22) using M modes is able to faithfully represent the Green's functions for λt M , except in a narrow energy interval of width λ/M around the Fermi energy. If the system has an energy gap ∆E, then the Green's functions are faithfully represented for all λt M as long as ∆E λ/M . Even for systems with a finite density of states near the Fermi energy, it is possible to faithfully represent the finite temperature Green's functions down to temperatures of λβ M and out to λt M/2 in time, as illustrated in Fig. 1.
The computational cost of the EPOCH method is mainly set by the computation of the matrix elements of P n (H), where each additional moment requires a matrixvector multiplication byH. For a sparse matrix, this is an O(N ) operation, that furthermore need only be done for the sites of interest and is also easily parallelized over multiple CPU or GPU cores. Moreover, the EPOCH methods also scales linearly in propagated time, i.e., running the time evolution for twice as long will, all else FIG. 4. Snapshots in time of the evolution of both the particle-particle |G < ij (t)| (a-d) and the particle-hole |Fij(t)| (e-h) correlation amplitudes for a particle inserted on the normal side near the interface of a 2D NS junction, with interface indicated by dashed white line. The square lattice is of dimensions L 2 = (201a) 2 and extends beyond the regions shown and the spatial scale in indicated by bar in (h). The evolution was computed using the EPOCH method at zero temperature with M = 1000 moments applied to the Hamiltonian in Eq. (23) at the chemical potential µ = 2γ and pairing potential ∆ = 0.25γ inside the superconductor.
being equal, take twice as long to compute. Thus, even for systems with large dimensions N the time-domain Green's function for large relative time differences are efficiently and easily computed computed with the EPOCH method.

VII. EXAMPLE: SUPERCONDUCTOR-NORMAL METAL INTERFACE
To illustrate the vast possibilities of the EPOCH method, while still studying a well known, but large system, we calculate the time-dependence of the normal and anomalous Green's functions in a superconductor-normal metal (SN) junction in both 2D and 3D. At the SN interface, an incident particle can, in addition to normal transmission and reflectance, also exhibit Andreev reflection [51], where a particle is converted to a hole and a Cooper pair is transferred to the superconductor. Andreev reflection is thus the process responsible for the superconducting proximity effect, which not only changes the properties of SN junctions but also directly leads to the finite Josephson current in junctions with two superconductors leads. Andreev refection is thus experimentally detectable in tunneling experiments [52], and with dual contact measurements the Andreev reflections are even by themselves directly observable [53,54]. Scattering at SN interfaces is therefore an interesting fundamental and experimentally measurable process, for which EPOCH is ideally suited. Specifically, we consider a finite sized tight-binding model on the square and cubic lattice, respectively, described bŷ where the particles (electrons) are created by the operator c † iσ at the lattice site i with spin σ =↑, ↓ at the position r i = a(i xx + i yŷ ) in 2D and r i = a(i xx + i yŷ + i zẑ ) in 3D, where a is the lattice constant and the total length of each sides of the lattice is L. With respect to the BdG form of Eq. (3), the diagonal normal partĤ 0 includes the first two terms ofĤ with the nearest-neighbor hopping amplitude γ and the chemical potential µ. Similarly, the last term is the off-diagonal part∆ of the BdG form, representing a constant superconducting order parameter of finite amplitude ∆ in half of the system x < 0, given by the Heaviside step function. This generates an NS interface at x = 0.
We are here interested in the equilibrium timeevolution of a particle (electron) that is initially created at time t = 0 close to the SN interface on the N side, and then allowed to propagate in time throughout the whole system. The physics of this evolution is captured by the Green's functions G < ij (t) and F ij (t), where the former captures the time-evolution of the intact particle, while the latter represents the amplitude associated with a particle at t = 0 appearing as a hole at a later time, t, i.e. a pair amplitude. Because the Hamiltonian in Eq. (23) is isotropic in spin, we can without loss of generality set G < ij (t) = i c † jσ (0)c iσ (t) to be between equal spins and the anomalous correlations are those of a spin-singlet superconducting state between opposite spins F ij (t) = i c j↓ (0)c i↑ (t) . We show below that the underlying physics is the same in both 2D and 3D, but because the results are easier to illustrate in 2D we treat this case first.
In Fig. 4 we show four snapshots in time of the evolution of both |G < ij (t)| (top row) and |F ij (t)| (bottom row) for a 2D NS junction of dimensions L 2 = (201a) 2 . The dashed white line marks the interface with the N side on top. Because the superconducting gap also extends into the normal metal through the proximity effect, we can compute the time-evolution using the EPOCH method even at zero temperature to high accuracy, here using a total of M = 1000 moments.
As seen in Fig. 4(a-d), the particle propagates outwards in approximately circular waves with some anisotropy due to the the anisotropy of the Fermi surface at µ = 2γ. As seen from the last two snapshots in Fig. 4(c-d) the propagation is largely unimpeded and unaffected even into the superconductor, but on the normal side there are clear interference patterns emanating parallel to the interface from a partially reflected wave at the SN interface. Studying the anomalous Green's function |F ij (t)|, we see that even at the earliest snap-shot in Fig. 4(e) there are finite particle-hole correlations on the normal side of the interface, thus appearing even before the main wave of the particle has reached the SN interface. This is a manifestation of the superconducting proximity effect. In the next snap-shot, the particle wave in Fig. 4(b) has reached the interface where it is transmitted but also Andreev reflected back in the N region as a hole. This shows up in the anomalous part in Fig. 4(f) as growing hole correlations at and beyond the interface. As the particle wave further propagates and extends along the SN interface, a part is transmitted as a hole into the S region due to the finite order parameter ∆, resulting in |G < ij (t)| and |F ij (t)| showing matching patterns in the S region in Fig. 4(c,g) and Fig. 4(d,h). At the SN interface, the intersecting wave crest of |G < ij (t)| is also Andreev reflected and the correlations propagate back into the N region, producing the triangular hole correlation wave front of |F ij (t)| in Fig. 4(h,g) in the N region.
Because of the O(N ) scaling of the EPOCH method there are no limitations in computing the time-evolved Green's functions also for a 3D bulk system. Thus, similarly to the 2D NS junction, we show in Fig. 5 snapshots of the time-evolution of |G < ij (t)| and |F ij (t)| for a 3D NS junction with the system dimensions L 3 = (125a) 3 . The dimension of the underlying BdG Hamiltonian matrix is therefore over 3.9 × 10 6 . Still, using the EPOCH method implemented on a standard laptop we generated Fig. 5 in a matter of minutes. This clearly illustrates the power and versatility of the EPOCH numerical method.
Specifically, in Fig. 5(a,f) we show the spatial extent of |G < ij (t)| and |F ij (t)| at t = 20γ −1 , respectively throughout the 3D volume surrounding the NS interface resulting from the propagation of a particle inserted on the normal side (blue). To illustrate more details, we also show snapshots at different times of the propagation of |G < ij (t)| in Fig. 5(b-e) |F ij (t)| in Fig. 5(g-j) on a vertical plane intersecting the creation point. From Fig. 5 it is clear that the same physical processes are present in both 2D and 3D. The main difference is only that the amplitudes of |G < ij (t)| in Fig. 5(b-e) rapidly diminish on the vertical intersection as the wave's amplitude now spreads out in all three dimensions, consistent with the propagation of spherical waves in 3D. Still, it is possible to discern from Fig. 5(g-j) that the approximately spherical particle amplitude wave is both Andreev reflected at the interface and hole transformed inside the superconductor, just as in 2D.
The results in Figs. 4-5 illustrate that the EPOCH method readily computes the time-domain Green's functions even for large 2D and 3D superconducting systems. The method therefore gives direct access to the dynamic equilibrium correlations of any normal or superconducting system, enabling the fully quantum mechanical timeevolution of condensed matter systems to be investigated. Notably, analyzing an NS junction in the time-domain offers a transparency that clearly contributes to an intuitive understanding of the underlying physical processes. Therefore, the EPOCH method provides a new tool for investigating phenomenon occurring in large highly inhomogeneous systems.

VIII. APPLICATION TO LINEAR RESPONSE
In this section we highlight yet another potential application of EPOCH, using it in the calculation of correlation functions directly in the time-domain, such as current-current J µ ( r, t)J ν ( r, t) and spin-spin S µ ( r, t)S ν ( r, t) to density-density n( r, t)n( r, t) correlation functions. All these dynamical correlation functions are interesting in their own right and, through the fluctuation-dissipation theorem, they also govern a system's response to external perturbations (not driving the system out of equilibrium), such as magnetic and electric fields [38]. For the latter, it is well established that calculations within the Kubo formalism can be carried out using the time-domain equilibrium Green's functions [38], with the susceptibility of an ob-servableÂ = ij A ij c † i c j to a field coupling opera- . Thus the expression for χ AB (t − t ) only requires the Green's functions that are directly computed in EPOCH, given by Eq. (2). As a result, dynamical correlation functions and linear response to external perturbations can be computed efficiently directly in the time-domain using the EPOCH method. Moreover, because the time and temperature dependencies of the response enters only via the transients l n β (t) = (−i) n (j n (t) + if n β (t)), a distinctive advantage of EPOCH is that the response through χ AB (t − t ) can be evaluated at any time or temperature without having to recompute any matrix elements of the quantum propagation. Thus, EPOCH applied to compute the time-domain linear response is more versatile and efficient than previous approaches that rely on a separate finite-difference stepping of the underlying Schrödinger equation, which is either only applicable at zero temperature [55] or, if extended to finite temperatures through the Boltzmann-weighted method, requires that the time stepping is redone at each temperature [56]. To conclude, EPOCH is a general method of obtaining the time-domain equilibrium Green's functions in large inhomogeneous systems from which any linear response function can be computed and it automatically gives all temperature and time dependence in a single computation.
In terms of the computational efficiency of EPOCH we note that, with response functions being central to many experiments, they have been a prime target for also other linear scaling methods [27,57]. This is especially the case in the field of quantum transport [58], where several new software packages have methods that also allows the conductivity to be computed in large systems [30,32,59,60]. Since the conductivity is given by the current-current correlation within the Kubo formula [61][62][63][64][65], it can also be directly computed using EPOCH. In fact, the direct time-domain approach of EPOCH is particularly interesting for quantum transport [58], because tracking of wave-packet diffusion in time allows for the identification of ballistic, diffusive, and localized regimes. This explains why the Chebyshev single wavepacket propagation method is still often used in quantum transport works [30,[66][67][68][69][70][71][72], where the conductivity can be computed from the velocity auto-correlation function [30,71] or the mean-square displacement [66][67][68][69]. However because single wave-packet propagation does not include the full quantum statistics, unlike EPOCH, any such wave-packet propagation relies on a separate projection on the Fermi energy, as well as normalization by a separately calculated density of states [58,72]. With EPOCH we entirely avoid these approximations and separate calculations, because the full quantum statistics is always included at the outset, thus allowing the conductivity tensor to be computed at any temperature and directly in time-domain, from which localization regimes, DC, and AC conductivities can easily be investigated.
In relation to the application in quantum transport, we also note that a significant simplification used in many other methods is that the important contributions are usually only from a narrow energy window around the Fermi energy, the so-called Fermi window. This simplification makes it possible to compute the conductivity using only a few low energy states, as e.g. in the Landauer-Büttiker formalism [73]. Using sparse matrix algorithms, the (M ) lowest energy states, or generally the scattering states of semi-infinite leads, can then be constructed with a linear scaling (O(M N )) in the system size (N ) [74]. One such approach is offered by the Kwant package [59]. A further benefit of using wave functions is that external time-dependent perturbations can be incorporated by directly propagating the wave functions using the time-dependent Schrödinger equation. Numerically the propagation can be done either with the unitary Crank-Nicolson for stability, or with Runga-Kutta (Linear Multistep Methods) as in the recent TKwant code [60,75,76]. Thus, for quantum transport problems, an approach based on scattering states is highly efficient, achieving a linear scaling on par with EPOCH. This favorable scaling is however only achieved when considering the narrow low-energy window, and therefore does not extend beyond this ideal transport regime. Since EPOCH always includes all occupied states, EPOCH does not only avoid the cost that comes from working with M wave functions, including computing the necessary overlap integrals and energy integration, but EPOCH can also be directly applied with linear scaling to problems where more than a few occupied states have a finite contribution.

IX. CONCLUDING REMARKS
We develop a computationally efficient method, EPOCH: Equilibrium Propagator by Orthogonal polynomial CHain, to extract the equilibrium thermal timedependent Green's functions directly in the time-domain that excels for any large inhomogeneous system described by a time-independent Hamiltonian. The method effectively generalizes the widely adopted and successful Chebyshev wave-packet propagation method already widespread in quantum chemistry for single particle pure quantum states to also handle many-body fermionic system in thermal equilibrium.
The centerpiece of the EPOCH approach is a quantity which we refer to as the equilibrium propagator (EP), whose matrix elements may be used to construct the wellknown single particle Green's functions of any fermionic state (such as electron, hole, or Cooper pair) at equilibrium at any fixed temperature. The main advantage of the EPOCH method lies in the efficiency with which this EP may be calculated using an expansion in orthogonal Legendre polynomials. All the time dependence in EPOCH is relegated to the coefficients of these polynomials, the mode transients, which may be written in terms of two parts, a unitary part, simply given by the predefined spherical Bessel functions, and a projective part. The projective transient encodes the quantum statistics by encapsulating all temperature dependence and is readily computed from an efficient recursive relationship that we derive. We arrive at the Green's functions by summing together the mode transients with matrix elements of the polynomial powers of Hamiltonian. These powers are efficiently calculated by single matrix vector multiplications due to the recursive relations of Legendre polynomials. Consequently, the EPOCH method scales linearly in the degrees of freedom of the system, while still both handling time-dependence and equilibrium fermionic statistics. The EPOCH method thus paves the way for investigating quantum phenomenon in large, inhomogeneous, even superconducting, systems, where time-evolution or dynamical aspects are of interest, ranging from quantum transport to odd-frequency superconductivity.
We illustrate the extreme efficiency of the EPOCH method by computing both the normal and anomalous Green's functions for a large 3D superconductor-normal state junction within minutes on a standard laptop. The results show how an electron created on the normal side of the junction propagates towards the interface where it can be either transmitted, normal reflected, or Andreev reflected, and finally builds up a complex pattern of particle-particle (normal) and particle-hole (anomalous or Cooper pair) waves in the junction. We have also recently employed EPOCH to calculate the propagation of even-and odd-frequency pair amplitudes in a disordered normal metal to superconductor junction to uncover unexpected robustness of odd-frequency superconductivity [77].
In summary, EPOCH is an extremely computationally effective method to calculate the equilibrium Green's functions directly in the time-domain to capture dynamics in quantum systems, with a wide range of possible applications.
to unequal, and unphysical, treatment of different regions. For example, the weight function of the Chebyshev polynomials is 1/ √ 1 − x 2 = 1/ (1 − x)(1 + x) and thus has diverging poles at the end of the interval I = [−1, 1]. A finite generalized Fourier series in the Chebyshev polynomials therefore sacrifices representational accuracy throughout the interval for a tighter fit around the end points, where the weight function adds most of the cost to any deviation.
A more quantitative statement of this weight function effect can be made if the generalized Fourier series is combined with the Jackson kernel, introduced in Sec. V A. In the Jackson resummation, the truncated Kernel of Eq. (8) is everywhere positive, approximating a Gaussian curve that narrows in on the ideal Dirac delta function, as the number of terms are increased. At the center of the interval I = [−1, 1] i.e. at E = 0, which is arguably the most important part of the spectrum of a condensed matter system, the width of this Gaussian kernel is π/N using N Chebyshev polynomials [27], but only π 2/3/N when using N Legendre polynomials. But even more importantly in time-evolution the whole spectrum contributes. Therefore, the main advantage of the Legendre polynomials is that they avoid placing any arbitrary extra weight to any particular region of the spectrum at the expense of an other region.