{'fraction_non_alphanumeric': 0.028102090686939994, 'fraction_numerical': 0.0067336410534890034, 'mean_word_length': 4.666307692307693, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Spreading of correlations in the Falicov-Kimball model

We study dynamical properties of the one- and two-dimensional Falicov-Kimball model using lattice Monte Carlo simulations. In particular, we calculate the spreading of charge correlations in the equilibrium model and after an interaction quench. The results show a reduction of the light-cone velocity with interaction strength at low temperature, while the phase velocity increases. At higher temperature, the initial spreading is determined by the Fermi velocity of the noninteracting system and the maximum range of the correlations decreases with increasing interaction strength. Charge order correlations in the disorder potential enhance the range of the correlations. We also use the numerically exact lattice Monte Carlo results to benchmark the accuracy of equilibrium and nonequilibrium dynamical cluster approximation calculations. It is shown that the bias introduced by the mapping to a periodized cluster is substantial, and that from a numerical point of view, it is more efficient to simulate the lattice model directly.


I. INTRODUCTION
Physical systems are characterized by their response to external perturbations. Slow or weak perturbations probe the equilibrium state of the system through response functions, which can be classified in terms of the low-energy excitations. Outside of this regime the nonequilibrium dynamics mixes excitations at different energy scales and can be complicated. In correlated quantum systems non-equilibrium studies reveal a plethora of new phenomena 1,2 and novel states of matter. [3][4][5][6] A general classification of universal features of nonequilibrium transport in quantum systems is currently lacking. Nevertheless, theoretical predictions exist for the spreading of correlations. Lieb and Robinson 7 showed that for interactions of finite range, there is a maximum velocity associated with this spreading. The resulting light-cone dynamics manifests itself in the commutators of observables, which are related to physical response functions, while anticommutators can exhibit algebraic tails that extend beyond the light-cone. 8,9 In the case of a noninteracting Fermion model the spreading velocity is determined by the maximum Fermi velocity. The effects of interactions and disorder modify this velocity and (in a localized phase) may limit the range of the correlations. A relatively simple model which allows to explore the influence of disorder and correlations on the light-cone dynamics is the Falicov-Kimball model. 10,11 This model describes mobile c electrons which interact with immobile f electrons. In equilibrium, this interaction produces a self-consistently determined disorder potential for the c electrons. Recent lattice Monte Carlo simulations 12,13 showed that the half-filled Falicov-Kimball model on the square lattice exhibits a rich phase diagram with metallic, weakly localized, Anderson insulating, Mott-like insulating, and charge ordered insulating (CDW) phases, Phase-diagram of the spin-less particle-hole symmetric Falicov-Kimball model on a two-dimensional square lattice, in the space of local interaction U and temperature T . The unit of energy is the nearest-neighbor hopping. The noninteracting system is metallic (M) while the interacting system at low temperatures is in a charge density wave (CDW) phase. At temperatures > ∼ 0.15, the system goes through weakly localized (WL), Anderson insulator (AI), and Mott insulator (MI) phases with increasing U . The WL phase appears in finite size simulations (here 48 × 48) and is replaced by the AI in the thermodynamic limit. (Lattice Monte Carlo results by Antipov et al. 13 ) see Fig. 1. The metallic (M) and weakly localized (WL) phases are conducting. Both the weakly localized and Anderson insulating (AI) phases are characterized by a nonzero density of states at the Fermi level and a finite localization length. The crossover between the AI and WL regimes occurs when the localization length reaches the system size. In the thermodynamic limit, the WL phase is replaced by the AI phase with a vanishing static conductivity. The Mott-like insulating (MI) and charge ordered insulating (CDW) phases have no available states in the proximity of the Fermi energy, either due to c-f electronic interactions (MI) or the presence of a static checkerboard order (CDW).
For a given f configuration, the dynamics of the noninteracting c electrons can be computed explicitly. Therefore the dynamics of the Falicov-Kimball model can be studied using a sign-problem free lattice Monte Carlo simulation of real-time correlation functions.
In this paper, we use this method to study timedependent correlation functions in the different equilibrium phases of the Falicov-Kimball model, and after interaction quenches between different physical regimes. We connect the maximum spatial extent of the spreading correlations and their phase and group velocities to characteristics of the underlying equilibrium phases. Specifically, Anderson and Mott-like localized phases show a finite range of correlations, and an increasing (decreasing) phase (group) velocity with increasing interaction strength. The onset of charge order changes this behavior and the correlations spread across the whole system.
We will furthermore use the numerically exact lattice Monte Carlo data to benchmark the dynamical cluster approximation (DCA) results in and out of equilibrium. [14][15][16] The dynamical mean field theory (DMFT) and its cluster extensions provide a numerically tractable tool for the study of nonequilibrium problems, 2 such as parameter quenches and electric field excitations. [16][17][18] However, relatively little is currently known about the effect of nonlocal correlations on the nonequilibrium properties of perturbed lattice systems, and on the ability of cluster DMFT methods to capture this physics. In the case of the one-dimensional Hubbard model, a systematic comparison between cluster DMFT and time-dependent DMRG has been presented in Ref. 17 for interaction quenches from a noninteracting initial state. It was found that the dynamics of local observables converges much faster with cluster size than the dynamics of nonlocal observables, and that a proper averaging over different cluster geometries can considerably improve the convergence with cluster size. The use of a weak-coupling perturbative impurity solver however limited these benchmarks to the weakly correlated regime. In the case of the Falicov-Kimball model, the cluster impurity problem can be solved exactly, so that the nonequilibrium dynamics can be simulated for arbitrary interaction strength. 16 This, in combination with exact Monte Carlo benchmark results, allows us to reveal the errors introduced by neglecting long-range correlations, or by assuming translational symmetry on the cluster.
The paper is organized as follows. Section II details the Monte Carlo simulation method and the dynamical cluster approximation for the Falicov-Kimball model. Section IIIA investigates the convergence of the lattice Monte Carlo and DCA simulations with system size, and provides benchmark results for equilibrium and nonequilibrium DCA. Section IIIB investigates the spreading of correlations in the one-dimensional and two-dimensional Falicov-Kimball model, while Section IV contains a discussion and outlook.

A. Lattice Monte Carlo
The lattice Hamiltonian of the Falicov-Kimball model in d = 1, 2 dimensions is given by where t is the hopping between nearest neighbor sites, U is the interaction parameter, and µ is the chemical potential (µ = U/2 corresponds to half-filling). The c electrons are itinerant, while the f electrons are localized, and the number of both species is conserved.
In the Monte Carlo treatment, we restrict the system to a d-dimensional cluster of linear size L with periodic boundary conditions. The partition function of the thermal system at inverse temperature β is given by where the trace over the f electron configurations enumerated by α is taken explicitly and n f α i ∈ {0, 1} denotes the f electron occupation in configuration α at lattice site i. H α = α|H + i µf † i f i |α = c † H α c is the c electron Hamiltonian with a site-dependent potential term defined by α. Here c † is the vector of creation operators (of size S = L d ) and H α an S × S matrix.
Since H α is quadratic in the fermionic operators, we can express the partition function as and we may interpret P α as an unnormalized probability amplitude for the f electron configuration α. The thermal expectation value of an observable O = c † O c can be calculated as where p α = P α /Z is a normalized probability density, and O α is the expectation value for the f electron configuration α.
With the above weights we perform Marcov chain Monte Carlo sampling on the f electron configurations. 19 For a given f electron configuration we use exact diagonalization to treat the remaining c electron problem. The expectation value of the observable O is expressed as where |ν is an eigenstate of H α with eigenvalue ν and n F ( ) = [exp(β ) + 1] −1 is the Fermi function.

B. Nonequilibrium
To study nonequilibrium phenomena we apply an instantaneous perturbation to the system at t = 0. In particular, we consider a global interaction quench where H is the equilibrium Hamiltonian as defined in Eq. (1) and is the time-dependent interaction parameter, defined in terms of the Heaviside function θ(t).
For ease of presentation we will refer to the equilibrium Hamiltonian as H − and to the post-quench Hamiltonian as H + . The lattice Monte Carlo method described above extends to the nonequilibrium case. For each f electron configuration we diagonalize the c electron Hamiltonians H ± α . The nonequilibrium Hamiltonian H + α is timeindependent so that the time-propagation can be performed analytically.
For example the lesser and greater components of the nonequilibrium Green's functions between sites i and j for a fixed f electron configuration α are given by where the time-propagators are defined as and the equal-time equilibrium Green's functions are given by where |ν − are the eigenstates of H − α and ν − are the corresponding eigenvalues.

C. Local and Nonlocal Correlation Functions
The lattice Monte Carlo procedure gives us access to various local and nonlocal observables in nonequilibrium. In particular we study the local c electron density ρ i (t) and double occupancy D i (t), which are given by We also study two-point density-density correlations between c electrons on the sites 0 and i, and at times 0 and t, For a fixed f electron configuration α we can apply Wick's theorem to evaluate the correlation function in terms of the single-particle Green's function Furthermore, we study the charge susceptibility, which can be defined as the Fourier transform of the commutator of the density-density correlations in both space and time and performing the Fourier transform analytically yields transforms to reciprocal space, and We only define the charge susceptibility in equilibrium, which is why we do not distinguish the quenched and equilibrium Hamiltonian in this case.

D. Spectral function and optical conductivity
The c electron spectral function can be obtained by a sampling over the eigenvalues αν of the c electron Hamiltonian matrix H α where we replace the δ-function by a finite width Lorentzian. For the zero frequency value it is desirable to avoid Lorentzian broadening in order to reliably identify a gap opening. In this case we can take a histogram over all occurring energy eigenvalues in a finite window around ω = 0.
Furthermore, we study the equilibrium optical conductivity, as previously described by Antipov et al., 13 obtained by linear response of the current to an applied infinitesimal electric field. Appendix A contains a derivation of the more general time-dependent optical conductivity.

E. Dynamical Cluster Approximation
In Ref. 16 we studied local and nonlocal correlations of the Falicov-Kimball model in nonequilibrium after an interaction quench using the dynamical cluster approximation (DCA). 20 Since we are interested in comparing the exact (i. e. converged in lattice size) results obtained using the above Monte Carlo method to results from that previous study, we briefly review the implementation of the nonequilibrium DCA formalism here.
DCA is a cluster extension of the dynamical mean field theory (DMFT). 2,14,15 A cluster of sitesr on the lattice is selected and defines an impurity model with a selfconsistently determined bath. Additionally, translation invariance and periodic boundary conditions are imposed on the cluster. In reciprocal space this yields patches of reciprocal vectorsk of the super lattice around reciprocal vectors K of the cluster. The corresponding effective impurity Hamiltonian reads where µ is the chemical potential (here µ = 0 corresponds to half-filling), N cl is the number of cluster sites, K are the reciprocal cluster vectors,¯ is the coarse grained dispersion,r are the cluster vectors, V K,p is the hybridization amplitude, ε K,p are the bath energy levels, and a ( †) are the bath annihilation (creation) operators. The impurity problem is solved by explicit summation over all the f electron configurations. For each such configuration α we obtain a c electron cluster Green's function R αKK . The final impurity Green's function is then given by the weighted sum over all f electron configurations where the w α are the equilibrium weights of the f electron configuration. The self-energy is approximated as constant on each patch in reciprocal space, Σ K+k = Σ K , and the selfconsistency condition requires the cluster Green's function to be equal to the coarse grained lattice Green's functionḠ The nonequilibrium problem is then solved on the Kadanoff-Baym contour. 16 We can determine the time-dependent double occupancy straight-forwardly on the impurity model. By applying Wick's theorem in a similar form as shown above on the cluster impurity we can also measure nonlocal density-density correlations It is important to note that this is an approximation and does not equal the corresponding observable on the whole lattice. However, for larger cluster sizes these correlations will converge towards the corresponding lattice observable.

A. Convergence of lattice Monte Carlo and DCA
We start by demonstrating the convergence of the lattice Monte Carlo results with increasing lattice size, and then use the converged Monte Carlo results to benchmark the DCA calculations. We consider the twodimensional model with nearest neighbor hopping, and use this nearest-neighbor hopping as the unit of energy.
1. Equilibrium Figure 2 shows equilibrium spectral functions for U = 2, 4 and temperatures T = 0.1, 0.2. The left panels show the lattice Monte Carlo results for indicated linear size L of the cluster, while the right hand panels compare different DCA spectra to the lattice Monte Carlo result for L = 32. At the higher temperature the system is in the disordered phase 13 for all interactions (see Fig. 1), and the Monte Carlo spectra are essentially converged for L = 24. At U = 4, the system is in the crossover region to the Anderson insulator regime, and a pseudogap opens in the spectral function. The comparison with the DCA spectra (right panels) shows an overall good agreement already for small cluster size. While the 2 × 2 and 2 × 8 clusters overestimate the pseudo-gap due to strong charge-order tendencies, the diagonal 8-site cluster underestimates it.
At the lower temperature T = 0.1, the infinite lattice system is in a charge-ordered insulating phase (see Fig. 1), with the U = 2 case close to the phase boundary to the Anderson insulator. Accordingly, there is a large gap in the spectral function. As demonstrated in the left panels of Fig. 2 the convergence of the Monte Carlo spectra with cluster size is slower, due to spiky features. Nevertheless it is clear that the L = 32 result is up to small oscillations converged. The DCA results qualitatively differ from the lattice Monte Carlo spectra, and rather resemble the high-temperature results. This is because of the translation invariance which is enforced on the cluster, and the suppression of long-range order in the self-consistency. The comparison between the DCA and lattice Monte Carlo results shows how the appearance of strong charge order correlations opens a gap of size U and shift the spectral weight from the gap region to sharp peaks at the gap edge.

Nonequilibrium
We next consider the time evolution after an interaction quench. To reduce discrepancies originating from charge order, we set T = 0.2. Figure 3 shows the simulation results for a quench from U = 3 to U = 4. The left panels illustrate the convergence of the lattice Monte Carlo results with the linear lattice size L (the width of the curves corresponds to the Monte Carlo error). Looking at the initital value, one sees that the result for the local observable N c 0 N f 0 converges very rapidly with lattice size, and that even the nonlocal quantities N c 0 N f 1 and N c 0 N c 1 are converged already for L = 8. While the time evolution for L = 2 and L = 4 exhibits spurious oscillations, for L ≥ 8 the dynamics shows a rapid damping and the curves are converged up to time t = 5. We can thus use these converged Monte Carlo data to benchmark the quench dynamics predicted by nonequilibrium DCA.
The second row compares the DCA evolution for different cluster geometries to the exact Monte Carlo result (black dashed curve). We first of all note that for the moderate cluster sizes considered, there is a strong cluster size dependence in the DCA results, especially for the nonlocal observables. While the results for the largest clusters tend to be relatively close to the benchmark curve, and the damping behavior is qualitatively well reproduced for the larger clusters, there is no systematic convergence with cluster size. The situation can be improved by calculating averages over different cluster geometries (patch layouts), as discussed in Ref. 16. In the third column, we compare these averaged curves  to the lattice Monte Carlo result. One now obtains a systematic convergence towards the exact result with increasing lattice size, although the off-set to the Monte Carlo curve remains substantial for cluster size 8 and in particular larger than what one might have guessed based on the difference between the "4×2" and "8" simulations. These deviations, which are also evident in the equilibrium spectra plotted in Fig. 2, may have several origins: the suppression of charge-order correlations in DCA, the lack of vertex corrections outside the cluster, or the fact that Anderson insulator behavior cannot be captured on small-size clusters (the U value after the quench is at the border of the Anderson insulator regime).
To better judge the accuracy of the damping behavior, we remove the offsets by subtracting the value at time t = 0 from all the curves. The corresponding results are shown in the right panels. We see that the damping dynamics is qualitatively well reproduced by the larger clusters, but differences remain at longer times and, for the nonlocal observables, even the initial response to the quench is not quantitatively accurate. While the absolute changes of the correlations are small, our results indicate that clusters with substantially more than 8 sites are needed to fully converge the DCA calculations. Comparing the 8-site DCA results to the 4 × 4 sites lattice Monte Carlo curve, we conclude that the DCA construction speeds up the convergence of the damping behavior, but results in a significant offset of the local and nonlocal correlation functions, so that the convergence to the exact infinite-lattice result is faster in the lattice Monte Carlo approach.

One-dimensional model
We start by discussing the properties of the equilibrium one-dimensional Falicov-Kimball model, where lattices up to L = 128 can be treated. Here we expect metallic, WL/AI and MI behavior with increasing U , even at low temperature. The simulation results are presented in Fig. 4 for different interactions U and two values of the temperature, T = 0.1 and 1.
In the right hand panels, we plot the density of states and the real part of the optical conductivity. The U = 2 system has a finite density of states at ω = 0, but a vanishing conductivity. According to the analysis of Ref. 13, this indicates an AI or WL state. (To distinguish the two, one would have to study the scaling of the conductivity with system size.) The states at U = 5 and 9 are characterized by a vanishing density of states and correspond to the MI phase.
To study the spreading of the density-density correla-tions with time we calculate the commutator The left two columns of Fig.4 show the results on a linear and logarithmic scale. The top four panels correspond to low temperature (T = 0.1) and the bottom four panels to high temperature (T = 1). At U = 0, the correlations spread with the Fermi velocity v F = 2 up to distances comparable to the system size.
For U > 0, the spreading of correlations is not linear in time any more. The correlations only extend up to some maximum distance, which decreases with increasing interaction strength, indicative of localization behavior. The localization length is the smallest in the MI phase. This behavior is in contrast to the Mott phase of the Hubbard model, which is characterized by freely propagating spin modes. In the Falicov-Kimball model, the spreading of the charge excitations is quenched in the disordered MI, in a way analogous to the AI phase.
We also notice that the spreading velocity is reduced with increasing interactions, which is most clearly visible in the low-temperature data plotted on the linear scale. To determine the spreading velocity, we look at the charge susceptibility χ k (ω) defined in Eq. (17), which is plotted in the third row. The maximum slope in the ω-versus-k curves defines the spreading (group) velocity v and the corresponding light cone is plotted (with an arbitrary offset) in the left-hand panels. While v can be rather easily determined in the low-temperature simulation results, the χ k (ω) plots for T = 1 show both dispersing and flat features, so that the definition of v becomes ambiguous. One noteworthy point is that in the high-temperature system, the Fermi velocity of the noninteracting model controls the spreading of correlations at short times, even in the AI and MI regimes.
In contrast to the group velocity, the phase velocity v , i. e. the velocity of the wave fronts, increases with increasing interaction strength. While there is a small speed-up of the phase velocity with time, it can be unambiguously determined at a given time in the low-temperature data. We extracted v by fitting the wave front emerging from t ≈ 5 as shown by the yellow line in the second row. If one indicates the corresponding velocity in the plot of χ k (ω), it matches the maximum intensity point in the susceptibility. At T = 1, the definition of v is more difficult, since at early times, the phase velocity is essentially given by v F , while at some later time, an enhanced phase velocity reminiscent of the low-temperature data appears, at least at large U . In the χ k (ω) plot, there are correspondingly two branches -the upper branch looks similar to the lowtemperature susceptibility, while the lower branch resembles the noninteracting dispersion. The large v roughly explains the edge of the upper branch.
The top panel of Fig. 6 shows the U -dependence of the velocities v and v for the one-dimensional model at T = 0.1. For U > ∼ 2, i. e. in the MI regime, the phase velocity v scales linearly with U , while the spreading velocity v scales roughly like 1/U .

Two-dimensional model
We now turn to the two-dimensional model, where simulations up to linear size L = 24 are possible with modest resources. The results for T = 0.1 and 1 are shown in Fig. 5. At the higher temperature, the U = 0, 2, 5, and 9 panels correspond to the M, WL, AI, and MI regimes, as evidenced by the density of states and conductivity data (see also Fig. 1). At the lower temperature, U = 2 is at the border between WL and CDW, while U = 5 and 9 are in the CDW phase. In the metallic phase (U = 0), the correlations spread at the fastest Fermi velocity in the x-direction, which is v F = 2. (Because of the small system size, the left and right wave fronts cross around t ≈ 6, which is evident at later times.) For U > 0, the low-temperature data reveal a reduced spreading velocity v , but up to time t = 10, there is no evidence for a finite localization length. The spreading behavior in the CDW phase is thus similar to the expected result for the Hubbard model. At T = 1, on the other hand, one observes a localization behavior in the WL/AI and disordered MI phases. Compared to the one-dimensional case, the phase velocity v is substantially larger, and near the edge of the light cone almost independent of U . However, at shorter distances, one can also identify slower phase velocities, as indicated by the yellow lines in the second column.
We plot the U -dependence of this slower phase velocity v together with the light cone velocity v in the bottom panel of Fig. 6, for T = 0.1 and 0.2. The phase velocity shows a jump between U = 5 and 6, which can be associated with the AI to MI transition at the higher temperature. Interestingly, the same feature persist even at T = 0.1, in the CDW insulating phase, which has longrange ordered f particle cofigurations. This indicates that the phase velocity is influenced more by correlation effects than by disorder effects. (Note that correlation induced changes near the Mott transition value of U occur also inside the CDW phase, similar to the crossover from weak-coupling to strong-coupling antiferromagnet in the Hubbard model. 21 ) In the high-temperature system (T = 1), the spreading at short times is controlled by the v F of the noninteracting system, even for large U . This is the same behavior as already observed in the one-dimensional case, and it is again reflected in the susceptibility in the form of two branches -a lower branch with a maximum slope given by v F and a weakly dispersing upper branch whose maximum slope is related to the spreading velocity at later times.
We have also studied the spreading of correlations after an interaction quench between the CDW/MI and WL phases of the two-dimensional model. The simulation results show a spreading behavior similar to a "cold" system at the final interaction, even though a substantial amount of energy is injected by the quench. They also provide further support for our observation that CDW correlations in the disorder potential enhance the spread- ing range. The details are presented in Appendix B.

IV. DISCUSSION AND CONCLUSIONS
We have studied dynamical properties of the Falicov-Kimball model using Monte Carlo simulation with exact time propagation. In particular, we focused on the spreading of density-density correlations in the half-filled one-and two-dimensional model and related the observed behavior to the equilibrium phase diagram. There are three quantities which characterize the spreading: the maximum range of the correlations, the spreading (or group) velocity v , and the phase velocity v . The value of these quantities depends on the interaction U , which determines the strength of the disorder potential, and on the CDW correlations in the disorder potential. We summarize the general trends in Tab. I.
In the U = 0 metallic regime, the correlations are not bounded, and they spread with the maximum Fermi velocity. For U > 0, in the disordered phase, the range of the correlations is limited, in accordance with Anderson localization. This holds both for the WL/AI phase at weak and intermediate U and the MI phase at large U . There is a systematic trend of decreasing correlation range with increasing U .
On the other hand, in the CDW phase of the twodimensional model, we do not observe any indications of localization behavior. While the accessible system size is limited, one can clearly conclude (by comparing the behavior above and below the CDW transition temperature) that the correlations in the disorder potential allow the density-density correlations to spread farther. In fact, in the CDW phase, the correlations spread without apparent bound, independent of U , and by fixing the f particle configuration to a perfect CDW pattern, one can furher increase the correlations, compared to the thermal ensemble. We note that in the thermodynamic limit of the model 22 multiple CDW phases can be present at different values of the interaction strength and temperature, and that their dynamical signatures may vary. In the CDW phase at T = 0.1 all possible localized states are gapped out.
We next consider the spreading (or light-cone) velocity v . It is obvious from the simulation data that this velocity decreases with increasing U , and that in this case, too, there is no dramatic change at the transition from the WL/AI phase to the MI phase. Even within the CDW phase, there is a systematic trend of decreasing v with increasing U . In fact, in this case one finds that the effects of U (disorder strength) and disorder correlations cooperate in reducing v . This is for example evident by comparing v above and below the CDW transition temperature, or the fact that a simulation with a perfect CDW disorder potential leads to a slower spreading than the simulation for the thermal ensemble.
As for the phase velocity v one finds a systematic increase with increasing U . This trend is most evident in the data for the one-dimensional model, where the wave fronts can be clearly identified. In the two-dimensional case, the phase velocity at the edge of the light cone becomes very large. However, at short distances and later times, one can also identify "slow" phase velocities, which resemble the behavior in the one-dimensional case. These v exhibit a jump at the transition from the WL/AI to the MI phase. Interestingly, the same jump is found even below the CDW phase transition temperature, which indicates that CDW correlations have no important effect on the phase velocity.
While the above analysis holds for low temperatures, the spreading behavior at higher temperatures is more complicated. Here, we find that the initial spreading is determined by the Fermi velocity of the noninteracting system, while at later times, it is impossible to measure v due to strong localization effects. Also the phase velocity becomes difficult to measure for U > 0. The charge susceptibility of the high temperature systems is characterized by two dispersing features, one corresponding to the initial spreading with velocity v F , and the other to localized charge excitations.
We also used the quench set-up to study the convergence properties of lattice Monte Carlo and DCA simulations. The Monte Carlo results converge rapidly with lattice size. Simulations on an 8 × 8 system are sufficient to produce essentially exact results for times up to t = 5, which is enough to observe a complete damping of transient oscillations in local and nearest-neighbor correlation functions. The DCA results, on the other hand, exhibit a very strong dependence on the cluster size and geometry on the small clusters (≤ 8 sites) that are accessible with our implementation. While averaging over different patch layouts improves the convergence with cluster size, substantial deviations from the exact result remain for these small clusters, and a reliable extrapolation to the thermodynamic limit is not yet possible. Given the slow convergence of DCA with cluster size, it is worthwhile to discuss the computational effort of DCA compared to direct lattice Monte Carlo simulations.
The DCA implementation used in this work and in Ref. 16   not suffer from such an exponential scaling. Experience shows that even a 24 × 24 cluster only requires about 2 13 Monte Carlo measurements. Additionally, the Monte Carlo procedure is trivially parallelizable, requiring no synchronization except for the final statistical analysis. It is possible to also implement DCA using Monte Carlo sampling over the f particle configurations, and it would be interesting to study how far such a Monte Carlo based scheme can be pushed.
Another essential difference is that the nonequilibrium DCA employs a time-stepping algorithm 16 that scales cubically in the number of time steps for computation and quadratically for memory requirements. Furthermore, the time discretization cannot be chosen arbitrarily, but has to be chosen small to ensure a converged solution. The Monte Carlo method, on the other hand, scales linearly in the number of time steps, for both computation and memory, and the time grid can be chosen arbitrarily due to the analytic time-propagation.
For comparison, simulating an 8-site cluster using DCA on 16 cores requires more than 100 GB of memory and takes about eight hours, while a Monte Carlo simulation of an 8 × 8 cluster on the same machine requires about 200 MB of memory and takes about ten minutes. In view of these considerations, it seems that the cluster DMFT approach does not offer any particular advantages in the study of the Falicov-Kimball model, and that the direct lattice simulation is the better strategy, even in the nonequilibrium case, or for the calculation of dynamical response functions.
where the interaction picture operators are given by The current expectation value is where is the current-current correlation function in the realspace components a, b ∈ {x, y, z}.
We define the conductivity σ as the system's response to the electric field where V is the volume. Taking into account that one finds We define the time and frequency dependent conductivity as σ(ω, t) = ∞ 0 ds σ(t + s, t)e iωs , and separate the Drude weight (or stiffness) D(t) as the dissipationless component such that Re σ(ω, t) = 2πe 2 D(t)δ(ω) + Re σ reg (ω, t).
When evaluating the Fourier transform we use the following transform of the Heaviside function where P indicates the Cauchy principal value. Finally, we arrive at the following expressions for the time and frequency dependent conductivity and the Drude weight at a fixed f electron configuration: To investigate the spreading of correlations in a nonequilibrium situation, we consider an interaction quench in the two-dimensional model from U = 9 to U q = 2, as described in Eqs. (6) and (7), and an initial equilibrium temperature of T = 0.1 (CDW phase). The energy injected into the system can be obtained from the instantaneous change in the local energy, i. e. the interaction and chemical potential contribution. From the total energy of the system after the quench, one may then calculate an effective temperature T eff , which corresponds to the temperature of an equilibrium system with U = U q and the same total energy. For the above set-up, one finds T eff = 0.791. 25 In Fig. 7 we plot the spreading of charge correlations after the interaction quench (top panels) and compare the result to an equilibrium simulation at U q and T eff (middle panels). The spreading velocity is somewhat lower in the quenched system than in the thermal system at T eff , and the correlations are less confined. This indicates an effectively "cold" (trapped) state of the quenched system. In fact, comparing to the T = 0.1 data in Fig. 5, it seems that the correlations spread in a way analogous to the cold CDW system at U = U q , despite the much higher total energy after the quench. This is possible because . Two-dimensional model: Spreading of two-point density-density correlations in a nonequilibrium system after an instantaneous interaction quench from U = 9 to Uq = 2 starting from an equilibrium system at temperature T = 0.1 (top row). Spreading in the corresponding equilibrium system at the temperature T eff defined in the text (center row).
Spreading for a fixed f particle configuration in a perfect charge-density-wave pattern. (bottom row).
the f particle distribution of the initial (U = 9) CDW state cannot adjust after the quench, so that the disorder potential corresponds to a CDW state.
As further support of the dominant role of the disorder potential, we have calculated the spreading behavior at U = U q and T eff for a fixed f configuration in a perfect CDW pattern (bottom panels). This system shows no localization, despite the elevated temperature. We also notice that the spreading velocity in the CDW potential is reduced compared to the quenched system and the equilibrium system at U q . This shows that CDW correlations in the disorder potential reduce the spreading velocity v .