F-electron spectral function of the Falicov-Kimball model and the Wiener-Hopf sum equation approach

We derive an alternative representation for the $f$-electron spectral function of the Falicov-Kimball model from the original one proposed by Brandt and Urbanek. In the new representation, all calculations are restricted to the real time axis, allowing us to go to arbitrarily low temperatures. The general formula for the retarded Green's function involves two determinants of continuous matrix operators that have the Toeplitz form. By employing the Wiener-Hopf sum equation approach and Szeg\"o's theorem, we can derive exact analytic formulas for the large-time limits of the Green's function; we illustrate this for cases when the logarithm of characteristic function (which defines the continuous Toeplitz matrix) does and does not wind around the origin. We show how accurate these asymptotic formulas are to the exact solutions found from extrapolating matrix calculations to the zero discretization size limit.


Introduction
The Falicov-Kimball model [ 1] has often been studied as one of the simplest models of strong electron correlations. Indeed, the exact solution in the limit of large dimensions displays charge density wave order, the Mott metal-insulator transition, and phase separation [ 2]. It is particularly relevant to present a paper on this model in an issue dedicated to Prof. Stasyuk, as he has worked on this model and closely related ones throughout his career; this began with his first article on the Falicov-Kimball model [ 3] while his most relevant work to this contribution involves approximate solutions for the f -electron spectral functions [ 4,5,6].
Our focus in this work is on the spectral function of the f -electrons in the exact solution of Falicov-Kimball model with dynamical mean-field theory. This problem was first investigated by Brandt and Urbanek [ 7] and by Janis [ 8]. The numerical approach of Brandt and Urbanek was extended by Zlatić et al. [ 9] and Freericks, Turkowski and Zlatić [ 10,11,12] within the original Brandt-Urbanek framework. An alternative approach was taken by employing the numerical renormalization group by Anders and Czycholl [ 13]. Some analytic work was completed by Liu [ 14] which was related to early work on this problem [ 15] and proceeded from the perspective of the X-ray edge problem [ 16,17]. Here, we develop a new approach which is numerically more tractable than other approaches and allows us to develop asymptotically exact formulas for the Green's function in the time representation.
The single-site Hamiltonian of the Falicov-Kimball model [ 1] involves two types of electronsconduction electrons, denoted by the letter d and localized electrons, denoted by the letter f where n d = d † d and n f = f † f are the number operators for the dand f electrons, respectively, U is the on-site Coulomb interaction between the d and f electrons, and E f is the site-energy of the f state. The full Hamiltonian on the lattice includes a repeat of this local Hamiltonian for each lattice site and a hopping of the conduction electrons between nearest-neighbor sites. The density matrix of the single-impurity problem is equal to

2)
where the time-ordering and integration are performed over the Kadanoff-Baym-Keldysh contour [ 18,19] and β = 1/T is the inverse temperature; the Kadanoff-Baym-Keldysh contour is shown in Fig. 1-it starts at t = 0, runs out along the real axis to a maximal time t max , then returns along the real axis back to t = 0, and finally runs along the negative imaginary axis down to −iβ.
Here, the dynamical mean-field λ c (t ′ , t ′′ ) and chemical potential µ are taken from the equilibrium solution of the conduction electron problem via dynamical mean-field theory [ 2]. In other words, the dynamical mean field satisfies (1.5) where λ(ω) is the ordinary dynamical mean field on the real axis, extracted from the dynamical mean-field theory solution of the model, f (ω) = 1/[1 + exp(βω)] is the Fermi-Dirac distribution function, and Θ c (t, t ′ ) is the Heaviside function on the contour, equal to 1 when t is ahead of t ′ on the contour, equal to 0 when t is behind t ′ on the contour and equal to 1/2 when t = t ′ on the contour (note that this definition of the λ c field is i times the definition in Ref. [ 2]). The partition function for the single-site problem contains two terms: which correspond to the different occupations of the f -particles, as follows: Here, iω m = iπT (2m + 1) is the fermionic Matsubara frequency and is the dynamical mean field evaluated at the mth Matsubara frequency [the λ c field is a function only of the difference of the two time arguments when both times lie on the imaginary part of the Kadanoff-Baym-Keldysh contour and we use the notation λ(τ ) = −iλ c (−iτ, 0), which corresponds to the conventional definition of the dynamical mean field for imaginary time].
is the f -electron spectral density. The corresponding Fourier transforms of the retarded and timeordered Green's functions are the following: (2.11) and We use the Kadanoff-Baym-Keldysh approach to calculate the real-time Green's function because there is no simple way to perform the analytical continuation of the Matsubara frequency Green's function to the real axis [ 7]. We are interested in the retarded Green's function, so we consider the greater Green's function for positive time t > 0 and the lesser Green's function for negative time t < 0 (ort = −t > 0) we use the fact that f † (0)f (t) can be replaced by f † (t)f (0) in the trace due to the time-translation invariance that we have in equilibrium. The greater and lesser Green's functions can be found for other time values by using the relations in equation (2.8). The first step in solving for the Green's function is to write equations of motion for the foperators: and substitute their solutions into equations (2.13) and (2.14), which yields the following expressions for the greater Green's function (t > 0) and for the lesser Green's function (t = −t > 0) where U c (t, t ′ ) is a new time-dependent field, which is nonzero only on the upper branch of the Kadanoff-Baym-Keldysh contour and is equal to U for t ′ ∈ [0, t] and is zero otherwise. It is because this U c field is not time-translation invariant that we need to use the Kadanoff-Baym-Keldysh formalism for the analytic continuation. Now we take into account the fact that the f -particle occupation number n f = f † f is a conserved quantity of the Hamiltonian; for the greater Green's function we have a projection onto states without f -particles (n f = 0) while for the lesser Green's function we have a projection onto states with one f -particle (n f = 1), which yields: and where the trace is now over the d-electrons only. We next expand the contour-ordered-exponent into a series that we can ultimately sum exactly. To begin, we need to recall Wick's theorem, where the expectation value of any number of pairs of fermionic operators can be expanded in terms of products of the two-operator contractions when the Hamiltonian is quadratic in the fermionic operators. Here we have ǫ 0 = −µ for sites with n f = 0 (α = 0) and ǫ 1 = U − µ for sites with n f = 1 (α = 1). The traces in equations (2.17) and (2.18) can then be written in the following diagrammatic expression 1 + e −βǫα exp where the arrows denote the zero-order Green functions g α (t ′ , t ′′ ) in equation (2.19) and the wavy lines denote a generalized time-dependent hopping (λ-field) for the greater Green's function and we have to replace for the lesser Green's function. Motivated by approaches that involve the integration over a coupling constant, we modify our time-dependent field by introducing a dependence on some new parameter that we denote by x with the following limiting behavior Next, we take the derivative of the diagrammatic series in equation (2.20) with respect to x and find: where we introduced a parameter-dependent Green's function G U (t ′ , t ′′ |x), which is the solution of the following Dyson equation: This expression holds for the greater Green's function (α = 0). For the lesser Green's function, one has to use the substitution in equation (2.22) and set α = 1. Finally, we find explicit expressions for the greater and lesser Green functions, where are the average densities of sites without (w 0 ) and with (w 1 ) f -electrons. We can rewrite equation (2.26) in the following two forms: The first form emphasizes the lambda field as an effective self-energy, while the second emphasizes the U field as the effective self-energy. In the top equation, we introduced the auxiliary Green's function as first defined by Brandt and Urbanek [ 7], while we introduced the bare time-ordered Green's function for the effective medium in the bottom equation. The bare time-ordered Green's function for the effective medium can be determined directly on the real axis via being the retarded effective medium on the real frequency axis. The first representation (2.30) was originally used by Brandt and Urbanek [ 7] and improved by Zlatić et al. [ 9] and Freericks, Turkowski and Zlatić [ 10] (see for review Freericks and Zlatić [ 2]) and requires the calculation of continuous matrix operator determinants over the entire contour.
In what follows, we shall instead use the second representation in equation (2.31) because in equations (2.27) and (2.31) all times in G α U (t 1 , t 1 |x) are only on the real axis (0 ≤ t 1 ≤ t) [or (0 ≤ t 1 ≤t)] and we can replace the contour integral in equation (2.31) by an integral over the upper real time branch of the Kadanoff-Baym-Keldysh contour where the time-dependent field Integrals over the entire Kadanoff-Baym-Keldysh contour survive only in the definition of G α (t) in equation (2.33), which can instead be calculated directly on the real axis by picking t and t ′ on the upper branch of the contour in equation (2.34) Also, for further discussions we will need its Fourier transform We are now ready to summarize the final formulas for the Green's function by evaluating all of the terms in equation (2.31) and then performing the integration over the parameter x. We begin with the formal matrix solution for G α Here δ c (t, t ′ ) is the Dirac delta function along the contour defined by After substituting this result into equation (2.27) and then integrating over the parameter x, we obtain for the greater Green's function the following results: where Tr c denotes the trace of a continuous matrix operator expressed as a line integral over the contour and det c is the determinant of the continuous matrix operator defined on the contour. However, in our case, the time-dependent field U c (t, t 2 ) (and hence all nondiagonal components) are nonzero only on the upper real time branch of the contour for 0 ≤ t 2 ≤ t, which allows us to reduce the problem to the determinant of a continuous matrix operator over the finite real time interval: It is the restriction to this finite real-time interval that makes the numerical calculations much more efficient than other methods used previously. Note that all of the temperature dependence is in the Green's function G 0 , which is the time-ordered Green's function for the effective medium of the dynamical mean-field theory.
In a similar fashion, we can derive the expression for the lesser Green's function (t = −t > 0): Here, the Green's function G 1 (t 1 − t 2 ) is defined in equation (2.37) and is the time-ordered Green's function for the effective medium when there is an f -electron on the impurity site. When we perform numerical calculations to evaluate the continuous matrix determinants, we first discretize the time contour with a discretization step ∆t, and then discretize the continuous matrix operator into a discrete matrix, which can be manipulated by standard computational libraries like LAPACK. The discretized representation for the Green's functions is then and where W i are the quadrature weights for the discrete set of times {t i }; for a uniform grid W i = W = ∆t, we have to calculate the determinant of a Toeplitz matrix. This will allow us to use Szegö's theorem and the Wiener-Hopf approach to find analytic expressions of the continuous matrix operator determinant which are accurate for long times. Note, that while it appears that these expressions are equally valid for all temperatures, the numerical solution, in the ∆t → 0 limit, becomes more difficult at lower temperatures. The retarded Green's function is found from the greater and lesser Green functions via Using the relation in equation (2.8), then yields for real times, and for real frequencies.

The Wiener-Hopf sum equation approach and Szegö's theorem
We follow closely the development of the Wiener-Hopf sum equation approach in McCoy and Wu [ 20], which is based on the original integral equation development by Wiener and Hopf [ 21], as summarized in Krein [ 22]. The Wiener-Hopf integral equation approach has been used widely in the X-ray edge problem and in the Falicov-Kimball model [ 8] primarily at T = 0; we will show that the Wiener-Hopf sum equation approach is very powerful for finite temperature calculations.
Before we can apply this technique to the problem of calculating the f -electron spectral function, we must begin with some mathematical preliminaries (further details can be found in Ref. [ 20]).
Consider a sequence of numbers c n which satisfy for 0 ≤ n ≤ N with c n and y n known. We will eventually be taking the limit where N → ∞, so we will require ∞ n=−∞ |y n | < ∞ too. Note that the matrix c n−m that appears in equation (3.1) is a Toeplitz matrix, because it depends only on the difference n − m and not on n and m separately. Our job is to find the solution x (N ) m to these equations. We define for n < 0 and n > N . In addition, we need to define two more sets of coefficients: We consider a variable ξ = exp[iθ] which lies on the unit circle, so that |ξ| = 1. We multiply both sides of equation (3.3) by ξ n and sum over n to find It may appear that we have made the problem more complicated, because we have replaced a problem with one unknown X(ξ) by a problem with three unknowns X(ξ), U (ξ) and V (ξ). But it turns out that this more complicated problem can actually be solved by carefully analyzing its analytic structure. The crucial step of Wiener and Hopf is to find a unique factorization of C(ξ) into two factors, one analytic inside the unit circle and continuous on the circle, and the other analytic outside the unit circle and continuous on the circle. Before we show how to do this, we need some more mathematical definitions. We call a function a + function if it can be expanded as a Laurent series for |ξ| = 1 ( ∞ n=0 a n ξ n ) and the coefficients satisfy ∞ n=0 |a n | < ∞. In this case, the function is analytic within the unit circle and continuous on the unit circle. Similarly, we call a function a − function if it can be expanded in a Laurent series for |ξ| = 1 ( −1 n=−∞ a n ξ n ) and the coefficients also satisfy −1 n=−∞ |a n | < ∞. This function is analytic outside the unit circle and continuous on the unit circle; it vanishes as |ξ| → ∞. Note that any function defined on the unit circle via a Laurent expansion f (ξ) = ∞ n=−∞ a n ξ n with ∞ n=−∞ |a n | < ∞ can obviously be uniquely decomposed into the sum of a + and a − functions via f (ξ) = f + (ξ) + f − (ξ).
The function C(ξ) is factorized in a so-called canonical factorization, where and both P (ξ) and Q(ξ) are nonzero for |ξ| ≤ 1; the factorization is made unique by requiring Q(0) = 1. The functions P and Q can then be found in terms of + and − functions as P (ξ) = e G+(ξ) and Q(ξ −1 ) = e G−(ξ) , (3.7) which automatically satisfies G − (∞) = 0 [and Q(0) = 1] and makes both P (ξ) and Q(ξ) nonvanishing inside the unit circle. An explicit calculation of G ± (ξ) is nontrivial, but it is easy to write down a formal solution via where the ± subscripts denote a restriction of the Laurent series for the logarithm to its nonnegative or negative power terms, respectively. Note that the condition for properly defining the + and − functions requires the logarithm to be single-valued after ξ winds around the unit circle, or, in other words, we require Ind C(ξ) = 0, where We will also need to consider situations where Ind C(ξ) = −1. In this case, a new functionC(ξ) = −ξC(ξ) will have zero index, and the canonical factorization can be applied toC(ξ) instead.
We are now ready to solve equation (3.4) by substituting in the factorization from equation (3.6). Next, multiply both sides of the equation by Q(ξ −1 ) to get (3.10) The term on the left hand side is a + function, because P (0) = 0, and [P (ξ)] −1 can be expanded in a power series with just positive powers for |ξ| < 1 and X(ξ) is a + function. Similarly, Q(ξ −1 )V (ξ −1 ) is a − function. The other two terms are neither + nor − functions, but they can be uniquely decomposed into + and − function pieces. We do this, and move all + functions to the left hand side and all − function pieces to the right hand side. We are left with The left hand side of equation (3.11) is an analytic function for |ξ| < 1 and is continuous on the unit circle, the right hand side is an analytic function for |ξ| > 1 and is continuous on the unit circle, and both functions agree on the unit circle, due to the equality, so we can define a function that is analytic in the entire complex plane, and hence is an entire function. But this function vanishes as |ξ| → ∞, and the only entire function that vanishes for large argument is the function that is identically equal to 0. Hence we learn that and Using the fact that X(ξ −1 )ξ N is a + function [recall X(ξ) is an order N polynomial in ξ], we can substitute ξ → ξ −1 in equation (3.11), multiply both sides of the equation by ξ N , and then separate into the + and − function pieces to create another vanishing entire function, and learn that both and hold. These are all of the necessary relations we will need for determining the Toeplitz determinants. Now we must take another mathematical interlude to state Szegö's theorem. This theorem determines how a Toeplitz determinant exponentially decays for large matrix size [ 23] and even finds the constant prefactor for the exponential decay [ 24,25]. The proof can be carried out using the Wiener-Hopf sum equation approach [ 20], but it would take too much space for us to repeat it here, and it is rife with many technical mathematical manipulations that would take us away from the physical phenomena we wish to discuss here. So we will instead simply state the result, and then show how to evaluate asymptotic expressions for the f -electron Green's function in real time.
If we consider the determinant D N of the N × N Toeplitz matrix defined by The conditions that Ind C(ξ) = 0, C(ξ) is continuous on the unit circle, and ∞ n=−∞ |c n | < ∞ are all required for the theorem to hold. Now, with all of the mathematical preliminaries complete, we are ready to apply these techniques to finding the Toeplitz determinants needed for the f -electron Green's function. We need to evaluate the two determinants in equations (2.43) and (2.44). We will show explicitly how to calculate the determinant for the greater Green's function. Modifications for the lesser Green's function are straightforward.
Case 1: no winding. We start by discretizing the time axis [0, t] with N + 1 points given by t n = n∆t, with ∆t = t/N . Then we define coefficients c n to satisfy c n = δ n0 − ∆tU G 0 (n∆t). (3.18) Since G 0 (t) is bounded, and decays exponentially fast in time, it is easy to see that ∞ n=−∞ |c n | < ∞. We will assume for the moment that Ind C(ξ) = 0 (which turns out will be true at half filling when U 0.866 on the hypercubic lattice, see below). Using the definition for c n , we find If we let n∆t = t ′ and ω∆t = θ, we recognize that in the limit where ∆t → 0, the sum in the right hand side of the above equation approaches the Fourier transform of G 0 (t ′ ), so we write (3.20) where G 0 (ω) is defined by equation (2.38). The coefficients g n in equation (3.17) can then be written as where we again used t ′ = n∆t. Now, defining D N to satisfy which is the asymptotically exact result in the limit of large t when the index of C(ξ) is equal to zero. Since the coefficient of t in the exponent can be a complex number, the determinant can oscillate while exponentially decaying at long times. This occurs when the system has undergone the Mott transition and the f -electron spectral function displays a gap at low frequency. But we can improve upon the result in equation (3.23) to provide finite-time corrections. We illustrate next how to do this using the Wiener-Hopf sum equation approach. We start with our original N + 1 simultaneous equations we need to solve in equation (3.1), with the coefficients c n defined in equation (3.18). Next, we choose y n = δ n0 . Cramers rule, applied to the first column of the c matrix immediately tells us that which further implies that If we assume that x (M) 0 = (1 − ∆tδ M ) exp(−g 0 ) (which we will verify below), then we find This will then provide the next order of corrections to the asymptotic limit of the determinant. We need to find x (M) 0 for large M . Our strategy is to set U (ξ) = 0 in equation (3.13) to find V (ξ −1 ), then solve for U (ξ −1 ) with equation (3.15) and the approximate V (ξ), and finally substitute U (ξ) into equation (3.12) to find X(ξ). The term x (M) 0 then follows from taking the ξ → 0 limit. Using the facts that Y (ξ) = 1 for our equation and Q(ξ −1 ) − 1 is a − function, we immediately find that equation (3.13) can be solved by Substituting into equation (3.15) (after replacing ξ → ξ −1 ) then yields (3.28) Now we must replace ξ → ξ −1 in the above equation. This requires care because we will convert the − function into a + function, but with no constant term in the power series expansion. So we write with the prime indicating there is no constant term in the + function. This result can now be substituted into equation (3.12) to yield where we used the fact that [Q(ξ −1 )] + = 1. We will need to evaluate this in the limit ξ → 0. Note that we have the following two identities for P and Q P (ξ) = exp − ∞ n=0 g n ξ n and P (ω) = exp The limit of ∆t is to remind us that there should be no constant term in the expansion. Using this strategy, we can immediately write down the final answer for x where t = N ∆t. Putting this all together finally yields the asymptotic expression for D(t)

36)
Case 2: winding of -1. Szegö's theorem fails when the index of C(ξ) is nonzero, and ln C(ξ) winds around the origin. This occurs on the hypercubic lattice at half filling when U 0.866, see below. We can still derive asymptotic formulas for the Green's function by solving an auxiliary problem, which has the winding removed. We now show how this is done. In all numerical cases we have examined, the index of C(ξ) satisfies Ind C(ξ) = −1 when U is large enough. Such a winding can be removed by consideringC(ξ) = −ξC(ξ) which has index zero (the minus sign is introduced for convenience, as will be clear below). The coefficients obviously satisfyc n = −c n−1 . (3.37) We examine the auxiliary Toeplitz determinant We need to calculatex (N ) N to find the determinant D N . First note that equations (3.12-3.15) hold for the barred functions when we have a winding around the origin. We will once again solve forV (ξ −1 ) settingŪ (ξ) = 0, then substitute into equation (3.14) to findX(ξ −1 )ξ N . Taking the limit ξ → 0 will then give usx (3.43) We take the limit ξ → 0 by integrating the function on the right hand side over the unit circle and dividing by 2π. Using θ = ∆tω, then yields where we used t = N ∆t again. Note that the limits on the integration cannot be simply extended to ±∞ as we did before. Since the function ln C(exp[iθ]) winds once around the origin, it displays a discontinuity of −2iπ to its imaginary part as θ runs from −π to π. We remove this discontinuity by adding a linear function that is equal to 0 at θ = −π and is equal to 2iπ at θ = π (the minus sign in C is needed to move the branch cut of the logarithm from the positive real axis to the negative real axis). Since this linear function cannot be extended to infinity, we need to work with a finite value of ∆t in the calculations (for a further discussion, look at Fig. 2 and the corresponding discussion below). We simply take ∆t small enough, that we see the numerical results do not change, and the system has approached its ∆t → 0 limit. The final result for the determinant is then The limits on the first integral have been extended to ±∞, because it turns out that the linear piece gives no contribution to the integral. The definitions ofḡ(t),P (ω) andQ(ω) are obvious from the above discussion. It should be noted that, in contrast to the case without winding where the third term in the exponent of equation (3.35) gives finite time corrections to the large t exponential asymptotics of equation (3.23), in the case with winding, equation (3.45) has a time dependent prefactor with a nontrivial ∆t → 0 and t → ∞ limit which can not be considered as a vanishing correction to the exponential decay of equation (3.23).

Numerical results
We begin our numerical discussion by describing the different behavior of C(ξ) when it has no winding or when it winds with an index of −1 around the origin. In Fig. 2, we show a plot of ln C(ξ) for two cases: (a) the case with no winding, and (b) the case with an index of −1. These two cases are distinguished by the sign of the real part of C(ξ) at ξ = 0, where its imaginary part changes sign. On the hypercubic lattice at half filling we have Re C(0) > 0 for U < 0.866 and there is no winding and Re C(0) < 0 for U > 0.866 and now the index is equal to −1. Notice how the imaginary part has a steep change in its value near ξ = 0 [panel (a)], which evolves into a discrete jump by 2πi when U 0.866 [panel (b)], the jump at the origin is compensated by the linear shift which runs from the minimal to maximal values of the frequency used in the calculations (|ω| = 10π here). While it is obvious that integration over ω can easily be extended to ω = ±∞ when there is no winding, one needs to carefully include contributions present at the given value of ∆t when there is winding, and one cannot extend the integration limits in this case. Next we compare our different asymptotic expansions to the results derived from a direct numerical calculation with the matrix formalism from equation (2.43) when T = 0.1. This involves a straightforward approach where we use three different discretization sizes (∆t = 0.1, 0.0666, and 0.0333) which we extrapolate to ∆t → 0, and we compare those results to calculations with different approximations. In the case with no winding, we use both the asymptotic form for large t in equation (3.23) and the more correct form for smaller t in equation (3.35). In Fig. 3, we show the results for the greater Green's function at half filling on the hypercubic lattice for U = 0.5 [panel (a)] and U = 0. 7 [panel (b)]. Note how the exponential form, that comes from Szegö's theorem, is quite accurate for large times, but becomes poor for small times, and how the corrections arising from the Wiener-Hopf approach produce remarkable agreement with the numerically exact results for essentially all t. The errors are the largest near t = 0, but even there, they lie below the few percent level. This shows that the analytic formulas are quite accurate for all t when there is no winding.
Next we examine what happens in the case with winding. We examine two values of U in   √ 2, and U = 2 which is a small gap insulator. We show the results from the scaled matrix calculations with a solid line (black) and the asymptotic Wiener-Hopf results with a (blue) dashed line. In this case, the Green's function has behavior that does not approach the asymptotic exponential behavior for short times, and the analytic formula is less accurate for small times (with maximal error is on the order of 10%).  Now we can specify the differences between the cases without winding and with winding. For U < 0.866 there is no winding and our exact results of the functional determinants rapidly approach the asymptotic result (with finite-time corrections) in equation (3.35) and we do not observe a crossing of the zero axis at any temperature and at T → 0 the long-time behavior is replaced by a power law. On the other hand, for 0.866 < U < U c we always observe a crossing of the axis at high temperatures, which originates from the time dependent prefactor in equation (3.45). With decreasing temperature, the crossing point shifts to larger times producing some kind of exponential decay for intermediate values of t and at T = 0 we observe a crossover to power law behavior when the zero crossing is pushed out to infinity. For U > U c the crossing of the axis is observed at finite time values for any temperature and does not transform into power law at zero temperature.
The Fourier transforms to the spectral formula produces results essentially equivalent to those shown already in Refs. [ 10,11,12], so we do not repeat them here.

Conclusions
In this work, we have shown a new representation for the f -electron Green's function, which allows for efficient numerical computation. By using the property that U c (t, t ′ ) is nonzero only on the interval [0, t], we restrict the determinant of the continuous matrix operators to a finite time interval instead of over the three branches of the Kadanoff-Baym-Keldysh contour. This produces a huge savings in the computational effort, because the size of the matrices do not grow with temperature and they are less than one half the size of the matrices used in previous numerical calculations. As a result, we can examine much wider regions of parameter space.
In addition, we used the Wiener-Hopf sum equation approach, along with Szegö's theorem, to derive exact analytical expressions for the Toeplitz determinants required to find the Green's functions. While these expressions are asymptotically exact in the limit where t → ∞, we find they have errors less than 10% all the way down to t = 0. This then allows us to use the matrix results for small times, and then append them by the asymptotic expressions for large times in order to find accurate results for the Green's function at all times. These are then Fourier transformed to real frequencies to determine the f -electron spectral function.
The f -electron spectral function displays the expected behavior as well. For U values smaller than the critical U for the Mott transition (U = √ 2), the spectral function develops a sharp peak at low T , which ultimately diverges as an inverse power law in ω as T → 0 [rigorously speaking, our asymptotic formulas are not valid at T = 0, because the function ln C(ξ) becomes discontinuous due to the jump in the Fermi factor, but the power law behavior can be extracted via other techniques-our focus here was on the finite-T behavior]. When U is larger than the critical U , a gap forms in the f -electron spectral function (rigorously speaking, it is a pseudogap on the hypercubic lattice) but subgap states rapidly enter as a function of T for low T .