{'fraction_non_alphanumeric': 0.026359576145021456, 'fraction_numerical': 0.006422045129462591, 'mean_word_length': 4.789758323474734, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Inhomogeneous spectral moment sum rules for the retarded Green function and self-energy of strongly correlated electrons or ultracold fermionic atoms in optical lattices

Spectral moment sum rules are presented for the inhomogeneous many-body problem described by the fermionic Falicov-Kimball or Hubbard models. These local sum rules allow for arbitrary hoppings, site energies, and interactions. They can be employed to quantify the accuracy of numerical solutions to the inhomogeneous many-body problem like strongly correlated multilayered devices, ultracold atoms in an optical lattice with a trap potential, strongly correlated systems that are disordered, or systems with nontrivial spatial ordering like a charge density wave or a spin density wave. We also show how the spectral moment sum rules determine the asymptotic behavior of the Green function, self-energy, and dynamical mean field, when applied to the dynamical mean-field theory solution of the many body problem. In particular, we illustrate in detail how one can dramatically reduce the number of Matsubara frequencies needed to solve the Falicov-Kimball model, while still retaining high precision, and we sketch how one can incorporate these results into Hirsch-Fye quantum Monte Carlo solvers for the Hubbard (or more complicated) models. Since the solution of inhomogeneous problems is significantly more time consuming than periodic systems, efficient use of these sum rules can provide a dramatic speed up in the computational time required to solve the many-body problem. We also discuss how these sum rules behave in nonequilibrium situations as well, where the Hamiltonian has explicit time dependence due to a driving field or due to the time-dependent change of a parameter like the interaction strength or the origin of the trap potential.


I. INTRODUCTION
Spectral moments are integrals of powers of frequency multiplied by the corresponding spectral function and integrated over all frequency. They can reveal important information about the structure and spread of the spectral function and, in some cases, can also reveal interesting information about different many-body correlation functions. In theory, knowledge of all spectral moments allows one to reconstruct the function, but that procedure is well-known to be unstable and is not commonly used in numerical calculations. The spectral moments also correspond to derivatives of the Green functions with respect to relative time (either real time or imaginary time), evaluated at the point where the relative time is zero.
As such, the spectral moments provide information about the relative-time dependence, when expanded as a Taylor series in time.
Spectral moment sum rules for the many-body problem were investigated in 1967 by Harris and Lange 1 shortly after the Hubbard model 2 was introduced. In that work, one can find the moment sum rules for the first three moments of the retarded Green function of the Hubbard model and various approximations like the alloy analogy problem (which is equivalent to an inhomogeneous Falicov-Kimball model 3 ). They also developed a strong-coupling projection method to find spectral moments within each of the different Hubbard bands.
This approach is an approximate one, as the projection operator is developed in a power series of the hopping divided by the interaction strength. Shortly thereafter, Nolting 4 applied the spectral moment sum rules to develop approximations for the momentum-dependent Green function that have two poles, with the weights and locations of the poles fixed by the corresponding sum rules. This approach has been developed quite extensively, and applied to a variety of different problems including dynamical mean-field theory 5,6,7,8,9 . That work examined ferromagnetic and antiferromagnetic long-range order in the Hubbard model, different approximation schemes for perturbation theory that produce the correct strongcoupling limit, and extended the sum rules to the retarded self-energy. The approach has also been applied to photoemission and inverse photoemission 10 , where it was recognized that the moments of the so-called lesser and greater Green functions play an important role.
Finally, Harris and Lange's projection technique was applied to the lesser and greater Green functions to examine spectral properties of the Hubbard model 11 .
Most of that work has had as its focus using the sum rules to develop different approx-imations, or to learn semiquantitative features of the many-body problem. But there is another application of spectral moment sum rules that is quite important for computational work. The sum rules allow, in certain circumstances, for numerically exact computations to be benchmarked against the sum rules. Steven White applied this to the Hubbard model in two dimensions 12 , and Deisz, Hess, and Serene applied it to the fluctuation-exchange approximation 13 at the Matsubara frequencies. The sum rules have been generalized to nonequilibrium situations for the Green functions 14 and self-energies 15 . In the case when one applies a spatially uniform, but time-dependent, electric field to the system, it turns out that many of the low-order moments are time-independent, even though the Hamiltonian has explicit time-dependence due to the field being turned on at a specific time (or because the field has nontrivial time dependence). Nonequilibrium dynamical mean-field theory 16,17,18 can solve the many-body problem exactly, and the sum rules are employed to benchmark the quality of the solutions 19 .
Finally, spectral moment sum rules, when viewed as a Taylor series expansion in time, are now being employed to improve both the speed and the quality of Hirsch-Fye quantum Monte Carlo approaches 20 for solving the impurity problem in dynamical mean-field theory 21 . In particular, recent work 22,23,24 has shown that employing the exact Taylor series expansions for short imaginary times allows one to use much larger Trotter error, yet achieve high accuracy with the Hirsch-Fye algorithm, so that it is competitive with continuous-timebased approaches 25 .
In this work, we want to extend the spectral-moment sum rules for the retarded Green function and retarded self-energy to inhomogeneous cases, which are becoming more important, and where one can find similar improvement in the efficiency of impurity solvers for use in inhomogeneous dynamical mean-field theory problems. There are four main thrusts of work in the inhomogeneous many-body problem currently: (i) examining the properties of strongly correlated multilayers because of their potential for quantum-mechanical engineering of device properties 26 ; (ii) examining the properties of ultracold atoms on optical lattices but spatially confined by a magnetic or optical trap 27 ; (iii) examining the properties of a strongly correlated material that is disordered; and (iv) examining properties of strongly correlated materials with inhomogeneous spatial ordering (like a charge or spin density wave). In the multilayer problem, most solutions have relied on approximate techniques for the Hubbard model, or have examined simplified models like the Falicov-Kimball model (recently, however, the numerical renormalization group has been applied to multilayer Hubbard systems 28,29 ). In the cold atom problem, little work has been applied to Green function-based techniques (although this is increasing); here the spectral-moment sum rules could aid both in benchmarking and in improving the accuracy and efficiency of numerical algorithms. The strongly correlated material with disorder problem has had many different techniques applied to it, but most of them are approximate in one way or another, so understanding the quality of the approximations is important. Less work has taken place in ordered phase calculations, especially for sum rules when the system is spatially ordered; these results can be immediately examined with the results presented here. In all cases, use of sum rules can help provide quantitative data to analyze how accurate the numerical solutions are. Finally, inhomogeneous systems are likely to be studied within the nonequilibrium context. Here, we also generalize the nonequilibrium sum rules to the inhomogeneous environment and we consider a wide range of different nonequilibrium contexts.
The remainder of this contribution is organized as follows: in Section II, we discuss the formalism for deriving the spectral moments in equilibrium; in Section III, we generalize to nonequilibrium cases; in Section IV, we discuss different applications of sum rules within dynamical mean-field theory; and in Section V, we present the conclusions and summary.

II. FORMALISM FOR DERIVING SPECTRAL MOMENTS IN EQUILIBRIUM
The equilibrium Hamiltonian for the inhomogeneous Falicov-Kimball 3 and Hubbard 2 models can be written in the following unified form: it is in the "slow limit" if the potentials correspond to an electric field, where the charge rearranges itself into a static redistribution in response to the potential rather than allowing current to flow). We will also examine a wide class of nonequilibrium cases below.
We do not make any assumptions about the translational invariance of any of the parameters that enter the Hamiltonian, so the Green function will generically depend on two spatial coordinates rather than on their difference. But in equilibrium, the system does have time-translation invariance, so we can describe the Green functions with a single frequency by making a temporal Fourier transform. The retarded Green function is defined to be The spectral function is then defined to be and the spectral moments become Similar to the homogeneous case 14,15 , one can easily show that note that this representation uses the fact that θ(0) = 1/2, so the factor of 2 is needed in order to make 2θ(0) = 1 and have the correct normalization. Using the Heisenberg equation of motion, one can relate the time derivatives to commutators with the Hamiltonian, to where In Eq. (7), we used a shorthand notation, where the angular brackets denote the trace over all states weighted by the density matrix exp(−βH)/Z.
Evaluating the commutators is straightforward, but tedious. The results are where n f i = f † i f i ; this result corrects a factor of 2 error in the third line of Eq. (29) of Ref. 15. Note, that in cases where we have inhomogeneity arising from the spatial long-range order (say charge or spin density wave order, for example), then the Hamiltonian is actually translationally invariant, but the sum rules have inhomogeneous results due to the explicit evaluation of the different expectation values (the charge density will vary from site to site in a charge density wave, for example). This then gives rise to spatially inhomogeneous moments.
Next, we examine the retarded self-energy moments. To begin, we start with the Dyson equation in the real space for the frequency-dependent Green function and self-energy: where G R0 ij (ω) is the noninteracting retarded Green function on the lattice. When the frequency ω is large enough, all Green functions and self-energies become purely real (on the infinite-dimensional hypercubic lattice, there can be an exponentially small imaginary part since the bandwidth is infinite, but this plays little role for large frequencies), and by using the spectral formulas, we can expand the functional dependence of the Green function and self-energies in terms of the corresponding moments (since the denominators in the integrals never vanish when the numerators are nonzero) yielding where are the spectral moments for the self-energy [Σ R ij (ω = ∞) is a real constant equal to the large-frequency limit of the self-energy]. The expansions in Eqs. (15) and (16) (and a similar expansion for the noninteracting retarded Green function) are substituted into the Dyson equation in Eq. (12) and then one equates powers of 1/ω to find where the matrixμ R n (R i , R j ) is the nth spectral moment of the noninteracting retarded Green function on the lattice. Those noninteracting moments are found from Eqs. (8)- (11) with U i = 0. Substituting the explicit values of the moments into Eqs. (18)-(21) finally yields the self-energy moments: Note that the algebra required to arrive at these results is nontrivial. In cases where i and j are farther apart than the range of the hopping matrix (which is often chosen to be nonzero only for nearest neighbors) many moments are identically zero, but nontrivial cancellations are required to ensure that all of the off-diagonal moments vanish (no local approximation has been made for the self-energy here-these results hold in all dimensions). The first self-energy moment is: which contains an off-diagonal term when one is in finite dimensions (it is local in infinite dimensions due to the scaling of the hopping matrix element with dimension).
The expression for the nonhomogeneous Green function moments Eqs. and then perform the averaging. If the system tends to self-average, not too many specific configurations would be needed, but if there is interesting physics arising from rare regions of the distributions, many calculations would be needed, and these calculations can get to be rather lengthy.

III. GENERALIZATION TO NONEQUILIBRIUM SITUATIONS
One rather general form for the nonequilibrium Hubbard-Falicov-Kimball Hamiltonian is In this case, we have added time dependence to all of the parameters in the Hamiltonian, but have not introduced any additional forms of interaction within the Hamiltonian. Nevertheless, this generalization allows for a rather rich class of nonequilibrium problems to be studied. For example, if we are examining a multilayered device with electronic charge reconstruction 30 (where the potentials V i and V f i are determined by an additional semiclassical Poisson equation), and we use a vector potential to describe an external electric field that drives current through the system 26,31 , then we would have a time dependent hopping determined by the Peierls substitution 32 . If we want to examine an interaction quench, as is often studied in cold atom systems, we would have a (typically harmonic) trapping potential V i and V f i and the interaction U i would become time dependent switching from one value for early times to another value for later times 33 such as would occur near a Feshbach resonance if the bias magnetic field is changed from one value to another (for some experiments, the potentials or hopping could also change when the coupling changes); the switching could be sudden as in a rapid quench, or adiabatic, with a slowly varying change, or anything in between. In addition, within the cold-atom picture, we could imagine creating time-dependent trap potentials V i (t) and V f i (t). This would allow us to examine what would happen if we applied an impulse to the atomic cloud, or if we shifted the origin of the harmonic potential from one spatial location to another, and then examined how the center-of-mass oscillates and damps back to the thermal state, or to some nonthermal steady state. Finally, we could examine the so-called Bragg spectroscopy experiment, where the optical lattice potential amplitude is oscillated with some set frequency and one observes things like the change in the momentum distribution after the system is driven for a certain period of time, or the change in the double occupancy. In this case, the hopping, the local potentials, and the interaction could all become time dependent, but can still be described by the general form of our Hamiltonian.
There are a few subtle points to keep in mind with these nonequilibrium problems. For example, in the case of a multilayered system with electronic charge reconstruction, the charge reconstruction is created in the distant past, so it corresponds to an equilibrium static potential which does not contribute to any flow of current. The field that drives the current is described by a time-dependent vector potential, which is turned on at a particular time.
One can examine the transient current flow or the steady-state current flow by examining how the system responds to the external field.
Given this general form for the Hamiltonian, we next need to derive the sum rules.
We work in a Heisenberg picture, because the operator algebra for the time-dependent creation and annihilation operators, at equal times, is unchanged from the standard fermionic commutation relations. The only difference from the nonequilibrium derivations worked out previously 14,15 is that here we need to work in real space for all calculations, because there is no translational invariance. In order to evaluate the expression for the nonequilibrium spectral moments, it is convenient to introduce the relative t = t 1 − t 2 and the average T = (t 1 + t 2 )/2 time coordinates for the retarded Green function Eq. (2). In this case, the physical time at which one would likes to calculate the moments will correspond to the average time T , and the Fourier transform to frequency space must be performed with respect to the relative time t. Then, one can define the nonequilibrium Green function moments by generalizing the expression in Eq. (5): In a similar way, one can define the nonequilibrium moments for the self-energy: (for more details, see Ref. 15). By using these equations, one can easily show that the expressions for the nonequilibrium retarded Green function moments in Eqs. Similarly, one can show that the expressions for the large frequency self-energy in Eq. (22) and for the zeroth self-energy moment in Eq. (23) remain the same in the nonequilibrium case. However, the expression for the first moment in Eq. (24) will acquire an additional term proportional to the second derivative of the large frequency limit of the self-energy Σ ij (T, ω = ∞). Namely, in the nonequilibrium case one finds: For most nonequilibrium problems we would consider, except for an adiabatically changing interaction, or an amplitude oscillation of the optical lattice, the derivative term would vanish almost everywhere.

MEAN-FIELD THEORY
There are two immediate applications of spectral moment sum rules within dynamical mean-field theory. The first is to use them to evaluate the high-frequency asymptotic behavior exactly and then supplement the high frequency results by numerical calculations at low frequencies, and the second, closely related, is to use them to evaluate the short imaginary time behavior exactly and supplement with long-time numerical calculations.
The latter has been already discussed within the context of the Hirsch-Fye quantum Monte Carlo algorithm for solving the impurity problem in dynamical mean-field theory and the results there show great promise as a means to improve the accuracy and the efficiency of calculations 23,24 . The basic idea (at half-filling) is that the curvature, which grows sharply with increasing U, immediately determines the short-time behavior for the Green function.
Because of the sharp dependence on τ , one would need to use a very small discretization step for the QMC to accurately describe such behavior, which would then be very costly in terms of computational time. Instead, one uses a coarser grid, but before performing the Fourier transform to Matsubara frequencies, one simply creates a finer grid and uses the short time relations to find the behavior close to τ = 0, and uses simple interpolation for other time values (a shape preserving spline would work well in this context). Then the Fourier transformation will much more accurately reflect the true behavior of the system, and all of the high-frequency structure will be properly recovered, so that the QMC can be used to determine the low-frequency data where it is most accurate. (It has already been demonstrated that this approach is competitive with other QMC techniques such as the continuous-time algorithm.) Issues of accuracy and efficiency will become increasingly important for inhomogeneous dynamical mean-field theory problems (like multilayers or ultracold atoms in a trap) because one needs to solve an impurity problem at each inequivalent lattice site of the inhomogeneous system. This can range from tens to hundreds of impurity solvers for multilayered systems to many thousands or more for ultracold atomic systems in a trap. We won't discuss the application within quantum Monte Carlo approaches further here, and instead will concentrate on examining a different application, which is to solve the Here ω n = π(2n + 1)/β is the fermionic Matsubara frequency. When the system is homogeneous, a Fourier transformation allows this problem to be solved immediately. When one has inhomogeneity in one dimension, as in multilayered structures, the local Green function (of an infinite device) can be found from the so-called quantum zipper algorithm 26,37 which expresses the Green function in terms of continued fractions. For finite cold atom systems, one can use LAPACK (or sparse matrix) routines to perform numerical matrix inversions to find the local Green function 38,39 . Whatever the technique, we assume that one can solve the Dyson equation to determine the local Green function. Next, the effective medium G 0 i and dynamical mean field λ i are extracted via the scalar equations and [All quantities in Eqs. (30) and (31) and Now we find the new local Green function and finally extract the local self-energy These equations are then iterated until they converge. The conduction electron filling sat- this result is not part of the dynamical mean-field theory iteration, but it is needed if one wants to update the chemical potential during the iterations to achieve a particular electron filling. Note that this summation is ill-defined and needs to be properly regularized (see below).
Typically, one chooses a set number of Matsubara frequencies, usually with an energy cutoff many multiples of the noninteracting electron bandwidth, and solves the (now finite) set of equations for the Green functions and self-energies by iteration starting from Σ i = 0.
One can try to include the effects of the neglected tails of the summations and infinite products to improve the accuracy and minimize the effect of the energy cutoff. By employing the exact sum rules, one can make this procedure work well. We describe this process next.
First recall that we have already proven that the Matsubara frequency Green function and self-energy satisfy and when the Matsubara frequency is large |ω n | ≫ |U max |+W int , where W int is the half bandwidth (in real frequency) of the interacting density of states (valid only for finite-dimensional systems). These results follow from the definition of the spectral moments, and the spectral formula for the Green function and self-energy. Using these two relations, substituting into Eq. (30), and recalling the definition in Eq. (31), produces This asymptotic expansion along with the expansion in Eq. (37) will allow us to treat the tails in the summation for the conduction-electron filling in Eq. (36) and in the infinite products needed for the localized-electron filling in Eqs. (32) and (33).
We imagine taking an energy cutoff E c which is larger than the interacting density of states half bandwidth. Using that cutoff to determine the explicit Matsubara frequencies solved in the numerical implementation of dynamical mean-field theory, we examine only Matsubara frequencies with |ω n | < E c . This defines a cutoff integer n c corresponding to the Matsubara frequency closest to E c but lying below it. The tail for the conduction electron /β can now be evaluated analytically for the first four moments that we have calculated. Define the approximate Green function via Then, the identity allows us to evaluate infinite sums of inverse powers of iω n noting that the sum of 1/iω n requires special regularization which is given by the result in Eq. (41). By taking derivatives, and evaluating at µ = 0, we immediately learn that Substituting into the expression for the conduction-electron filling then yields The summation, which needs to be evaluated numerically, vanishes rapidly for large ω n , so the filling can be computed quite accurately.
Now we show how to evaluate the infinite products by properly taking into account the asymptotic limits using the spectral sum rules. Substituting the results from Eq. (39) into Eqs. (32) and (33) allow for the local electron filling to be computed. First note that and The infinite products that have an infinite number of terms are approximated by rewriting the infinite product as the exponential of the sum of the logarithm of the individual terms.
Replacing the sum by an integral (valid when the temperature is much smaller than the interacting half bandwidth) and converting the integral over frequency to an integral over where we use a = µ i for Z 0i , a = µ i − U i for Z 1i , b = −1/2 for both, and c = −(U i n f i − µ i )/2 for both. This then allows for an accurate evaluation of the filling using just the small set of numerical data generated for |n| < n c . Note that our use of the asymptotic expressions for the different many-body functions aided us in reducing the effort of computation only for problems that can be solved along the imaginary axis. This includes the Falicov-Kimball model (for static properties) and Hirsch-Fye quantum Monte Carlo techniques for other models (like the Hubbard model).
There does not appear to be any simple use of these sum rules within real-frequency-based approaches like the numerical renormalization group. Of course, all of these results can also be applied to the homogeneous case.
In Fig. 1, we plot false grayscale images of two moment sum rules for each lattice site of a 51 × 51 square lattice. This is a system close to phase separation with parameters U = 5, the computational time by about a factor of 10 for the same accuracy. In this real axis calculation, we use a frequency grid with a step size of 0.004t that runs from −9.6t to 9.6t.
This step size does not allow us to pick up fine structure smaller than the step size. It also cannot pick up spectral weight lying outside of the bounds. Both of these issues can cause inaccuracies in calculations, especially when a system begins to order or phase separate; they don't enter too significantly for this example though; when we repeat the calculations with a smaller step size of 0.0003t, we find the error in the zeroth self-energy moment is reduced from 1.5% to 0.025% for the central site of the lattice. The temperature we use here is high enough that there is no order, and the calculations are under good control. In any case, we show the spectral function and −ImΣ i (ω)/π for the central site of the lattice in Fig. 2, where the relative error for the zeroth moment of the self-energy is about 1.5%, and the relative error for the Green function is about 0.01%. Notice the sharp peak in the self-energy which leads to the higher error, while the Green function is quite smooth, and hence has very small errors. As the step size is reduced, the errors in the moments are also reduced, indicating that these errors are arising predominantly from the discretization size of the real frequency axis and the structure of the sharp features in the functions. We also examine the sum rules for inhomogeneous multilayered systems. In cases where there is no electronic charge reconstruction, we have found that our data satisfies all of the sum rules to high accuracy (typically better than 0.01%) except for cases with an insulating barrier where the self-energy develops a sharp peak at low frequency. Our frequency grid in previous calculations was sometimes too large to properly extract the zeroth moment sum rule for the self-energy, and errors could become very large because the numerical quadrature is greatly overestimating the weight within the sharp peak near ω = 0. In cases where the peak is not so sharp, we once again find excellent agreement. A more challenging case, though, is a case where there is an electronic charge reconstruction, because the calculations become much more difficult numerically in this case, and we usually need to introduce a finite broadening into the calculation to be able to estimate the local DOS on each plane.
Hence, it is much more interesting to examine these cases for calculations of the sum rules. Note how the DOS is once again smooth and rather featureless, so the sum rules work to high accuracy, while the self-energy has a sharp peak (maximum amplitude is 8000), which leads to less accurate moments.
In Fig. 3, we plot the absolute errors of a self-energy moment and a spectral moment ranged ω from −11t to 11t. One can see that, similar to the cold atom example above, here we also see errors which are much larger for the self-energy than for the spectral function (we only show the barrier planes for the self-energy moment, since the self-energy vanishes in the metal). This arises for the same reason as before, but is more acute here, since the step size in frequency is larger. This is shown in Fig. 4 for the plane with the largest error in the self-energy. Note how once again the DOS is smooth, which is why the moments are so accurate, but the self-energy has a narrow peak, whose weight is overestimated with the coarse grid used in the calculation.
where the sum over momentum is over the CDW Brillouin zone, which satisfies ǫ k ≤ 0, and the plus sign is for i ∈ A sublattice and the minus sign for i ∈ B sublattice. The time evolution operator U is a time-ordered product We used the Peierls' substituted band structure with ǫ k = lim d→∞ t * d i=1 cos(k i a)/ √ d the such long-ranged tails in time, one can see an effect even before the average time where the field has been turned on. When the field is turned on, the DOS naturally changes shape (and the square-root singularity appears to be smoothed into a finite peak), but the changes are much smaller than in the normal phase. Surprisingly, we predominantly see structure at either odd Bloch frequencies or even Bloch frequencies, but not both, and the structure does not look like a broadened, and split delta function as in the normal phase, but instead looks more like some kind of ordered phase gap structure, but the features can be quite small in some cases. It is clear we have not yet fully reached the steady state, but it also appears clear what the steady state DOS will eventually look like. Finally, we comment that if we take the DOS as functions of frequency (our frequency range is chosen to run from −10 < ω < 10), multiply by the appropriate power of frequency and integrate to check the sum rules, then all sum rules except the second moment are satisfied to better than 1%.
The second moment is worse, because we get significant contributions from the noisy tails of the DOS which don't cancel as they do for the odd moments. If we put the frequency cutoff closer to |ω| < 5, then the second moment accuracy increases dramatically. But we know that this is the least accurate way to check the moment sum rules.

V. CONCLUSIONS
In this work, we have shown how to determine exact spectral moment sum rules for the retarded Green functions and self-energies of inhomogeneous strongly correlated systems.
While our analysis is quite general, for concreteness, we provided explicit results for a combined Hubbard-Falicov-Kimball model. We envision these results can be applied to strongly correlated multilayers, ultracold atomic systems in traps, or strongly correlated materials with disorder (although we did not discuss the last case in much detail). The sum rules can be used to gain qualitative information about the many-body solutions, benchmark numerical calculations, be used to evaluate the tails of infinite sums (or infinite products) allowing for smaller energy cutoffs, or be used to improve the accuracy of Hirsch-Fye quantum Monte Carlo approaches by determining the short-imaginary-time behavior exactly. We provided a full derivation of the sum rules in equilibrium, and then discussed a number of different nonequilibrium situations appropriate for multilayers, for cold atom systems, and for ordered phase systems (such as a CDW). In the case of multilayers, the change in the local electrical potential energy, induced by an electronic charge reconstruction, modifies the sum rules, but the additional scalar potential, that creates an electric field to drive current through the device, does not. This motivates one to decouple the description, using a scalar potential and Poisson's equation in a semiclassical analysis to determine the modified Hamiltonian due to the electronic charge reconstruction, but use a time-dependent vector potential to describe the electric field that drives current through the system. In the case of cold atoms, we discussed nonequilibrium situations corresponding to pulling the lattice in the presence of a static trap, rapidly changing the location of the trap and examining the transient response, and changing the interaction strength (say due to a Feshbach resonance, by changing the magnetic field) and examining the response to an interaction quench. For the ordered phase, we examined the CDW state of the Falicov-Kimball model at zero temperature, with an uniform electric field turned on at a specific time. We feel these sum rules will have a broad application across a number of different systems and calculations.
It is obvious that a similar exercise can be carried out for bosonic systems. These will have most relevance for cold atom problems, where bosonic atoms are readily available within the alkali family. The sum rules for bosons become more complicated, especially at higher orders, because we used identities such as c † i c i c † i c i = c † i c i which hold only for fermionic systems. We are currently working on the generalization of these sum rules to bosonic problems which will be presented elsewhere. Of particular interest is a Falicov-Kimball model with light fermions and heavy bosons, which is realized in K-Rb mixtures used to make dipolar molecules if the system is placed on a deep enough optical lattice.