# Three-dimensional dynamics of a fermionic Mott wedding-cake in clean and disordered optical lattices

###### Abstract.

Non-equilibrium quantum phenomena are ubiquitous in nature. Yet, theoretical predictions on the real-time dynamics of many-body quantum systems remain formidably challenging, especially for high dimensions, strong interactions or disordered samples. Here we consider a notable paradigm of strongly correlated Fermi systems, the Mott phase of the Hubbard model, in a setup resembling ultracold-gases experiments. We study the three-dimensional expansion of a cloud into an optical lattice after removing the confining potential. We use time-dependent density-functional theory combined with dynamical mean-field theory, considering interactions below and above the Mott threshold, as well as disorder effects. At strong coupling, we observe multiple timescales in the melting of the Mott wedding-cake structure, as the Mott plateau persist orders of magnitude longer than the band insulating core. We also show that disorder destabilises the Mott plateau and that, compared to a clean setup, localisation can decrease, creating an interesting dynamic crossover during the expansion.

A. Kartsev, D. Karlsson, A. Privitera and C. Verdozzi

## Introduction

A large part of our understanding of the physical world concerns the equilibrium state, where macroscopical observables are constant in time. However, most phenomena surrounding us are instead non-equilibrium phenomena, where systems evolve towards new equilibrium states once perturbed away from their initial conditions. For condensed matter systems, ultrafast techniques are now becoming available to probe the non-equilibrium dynamics at the typical electronic scales [1, 2]. Even so, the intrinsic complexity of real materials still makes their theoretical understanding extremely difficult, because the electronic effects are usually interwoven with phononic contributions, inhomogeneities, etc.

In contrast, idealised solids realised by loading fermionic ultracold atoms into optical lattices feature several advantages over the corresponding solid-state systems. The accurate control of the system parameters permits experimental realisations of model fermionic Hamiltonians, e.g. the single-band Hubbard model [3, 4], where higher-bands contributions, phonons, lattice defects are absent or can be introduced separately and in a controlled way. In fact, most system parameters can be tuned independently in time and in a very wide range, allowing for the exploration of regimes beyond what is feasible in condensed matter [5, 6]. Moreover, the typical timescales of ultracold-gases dynamics are much longer () than those in condensed matter systems [7]. This has made possible, for example, the observation of interesting phenomena like the crossover from ballistic to diffusive behaviour in a fermionic cloud expanding into an optical lattice [8]. Although the physics in these experiments is much simpler than in condensed-matter ones, the theoretical understanding is still in general impaired by the absence of reliable tools to address the non-equilibrium behaviour of these simple model fermionic systems.

Indeed, among the well established approaches, perturbative techniques are only suitable for weak- and strong-coupling regimes, while non-perturbative schemes scale unfavourably with system size and simulation time. On the whole, it is fair to say that using any of the methods above for large inhomogeneous samples does not appear immediate, because of rather prohibitive Êcomputational costs. Lastly, exact techniques like the time-dependent density-matrix renormalisation group are currently mainly restricted to one dimension, while Quantum Monte-Carlo approaches can be severely limited by the fermionic sign-problem. Therefore, progress in the field requires new methods to cope with higher dimensions, large inhomogeneous systems irrespective of the interaction strength.

For these challenging situations, we advocate here the use of a recently developed approach [9], i.e. the combination of Lattice Time-Dependent Density-Functional Theory and Dynamical Mean-Field Theory (Methods). Within this method, one can properly describe the Mott transition within a density-functional theory framework [9]. Here, we show its potential by addressing the real-time dynamics of a large inhomogeneous three-dimensional system (up to lattice sites) in different regimes of interaction strengths and also in presence of disorder.

Our setup resembles a recent experiment[8] on ultracold Fermi gases, consisting of a confined cloud, which expands into an optical lattice after the trap is removed. Our initial density profile is the ground-state density of the trapped interacting cloud. The corresponding initial state is the (Kohn-Sham) ground state of density-functional theory (Methods). To describe the expansion of the system, we study the time-dependent one-particle density, which can in principle be obtained exactly within our approach (Methods).

We considered two interaction regimes, corresponding to a strongly-correlated metal and to a Mott insulator in the homogeneous case, and different protocols for switching off the trapping potential. At strong coupling, the ground-state density exhibits the peculiar wedding-cake structure due to the presence of the Mott plateau where the local density is commensurate [10]. As we release the trapping potential, the high-density metallic domain immediately melts, while the Mott plateau remains remarkably stable against the expansion, over much longer timescales than below the Mott transition. Thus the intrinsic inhomogeneity of confined ultracold gases displays the multi-scale dynamics of different phases in a single experiment.

We also considered the role of disorder on the cloud expansion. Compared to the static case[11, 12, 13, 14, 15], the effect of disorder on interacting Fermi systems out-of-equilibrium is much less understood, especially for dimensions larger than one. Our results show that disorder earlier makes the Mott plateau less stable, decreasing the melting time, while slowing down the expansion at long times. This induces a noteworthy crossover in real-time.

As commonly done within the ultracold-atom community, the system is described in terms of an inhomogeneous and time-dependent Hubbard model. In standard notation,

(1) |

with the hopping parameter and the contact interaction strength. The trap strength determines the switching-off protocol in time, and is nonzero in the disordered case (Methods). We take as energy unit, thus is in units of .

In our simulations, we consider a simple-cubic lattice with and . In the homogenous case, these values correspond respectively to a strongly-correlated metal and to a Mott insulator (the critical interaction ). Starting from the trapped system in the ground state, we examine three different expansion scenarios: For , we study a sudden expansion (), as well as a slow one (), where controls the trap removal speed (Methods). For , instead, we choose to consider only the slow expansion.

Indeed, the current implementation of our approach is expected to give a reliable description for slow and moderately-fast expansions, even at strong coupling. However, for fast trap removals, its performance significantly deteriorates for increasing interaction, and large deviations occur well above the Mott (gapped) regime. Instead, in the metallic (gapless) regime at , our results should remain robust in a qualitative sense. This is further addressed in the Methods section.

## Results

Expansion into a homogeneous lattice. The qualitative behaviour in time of the density profile in the plane for the chosen setups is shown in Fig. 1. Even within a cubic lattice, the initial cloud profiles appear to be in very good approximation spherical, due to the symmetry of the trapping potential. In Fig. 1a-b, for , the initial density profile smoothly decreases as a function of the distance from the centre in every direction. At strong-coupling instead (, Fig. 1c), the repulsion is large enough to induce an insulating phase in the homogeneous system at half-filling. Accordingly, the trapped system develops a Mott plateau, i.e. a region where the density is to a very good approximation, due to the incompressible nature of the Mott phase.

It is immediately manifest that the Mott physics has striking effects on the dynamics. Indeed, for the cloud expands smoothly (Fig. 1a-b), while for (Fig. 1c) the band insulating core immediately melts to get rid of the large interaction energy, and the upper part of the wedding-cake structure at collapses over the Mott plateau. At the same time, the underlying Mott plateau is remarkably stable against the trap opening and, after some time, the system exhibits a large Mott region surrounded by a metallic domain (Fig. 1c, ). Only on a much larger timescale (), also the Mott plateau melts and the system fully relaxes into the lattice.

Different switching-off protocols are compared in Fig 1a-b. For the sudden case, the system evolves only according to the homogeneous Hamiltonian, which (together with the initial state) sets the expansion timescales. As a consequence, the metallic dome melts much faster, and the shape of the cloud develops clear signatures of the lattice symmetry (). Conversely, in the slow case, the expansion rate is controlled by the trap-removal speed, and thus the cloud anisotropy remains small.

In Fig. 2, we analyse the expanding cloud in a more quantitative way. We start by looking at the density as a function only of the radius (Fig. 2a-c). For the slow expansion (Fig. 2b-c), at the cloud core is basically spherical (i.e. the density is a single-valued function of ) and tends to maintain this symmetry. Significant anisotropy is instead observed in the low-density region around the cloud boundaries for the sudden setup (Fig. 2a). In this case one expects the tail expansion to be ballistic [8], with the cubic symmetry of the underlying lattice.

These considerations are further supported by the analysis (not shown) of an Óasymmetry parameterÓ , where and are the maximum and minimum radial size of the cloud (). The cloud radius is defined as the largest (shortest) distance from the centre where the density is above ( was also used to check that the expanding cloud did not reach the boundaries of the simulation box). The maximum asymmetry we found is for the sudden case, and for the slow expansions. Finally, a notable feature in Fig. 2c is that the Mott plateau, before collapsing, widens over a significant time interval ( ). This appears to be independent of dimensionality [16], and thus related to the intrinsic energetics of the Mott phase and its rigidity.

Another useful quantity to analyse is the cloud expansion velocity. This was done in ref. [8], where, starting from an initial band-insulator state, the dual nature of the expansion was characterised as ballistic at the cloud edge, and, due to interactions, as hydrodynamic in the cloud core. In our analysis, we compare the maximum and average radius of the cloud , where is the number of particles. Their different behaviour (Fig. 2d) can give indications on the presence of multiple domains expanding at different rates. For and in the case of a sudden expansion, and were fitted according to the expression , where generically denotes () at , and the speed is a fitting parameter. The fit of confirms that the edge of the cloud expands ballistically [8]. At the same time, taking into account the slower interacting core, increases at a lower rate than . For the slow expansions, also shown in Fig. 2d, no fit was attempted, since the the slow trap opening hinders the expected ballistic behaviour at the cloud tails.

The effect of different trap-removal protocols on the dynamics in the Mott regime is shown in Fig. 2e for . We find it informative to scale time according to . In this way, simulations for different :s appear very similar, and can be discussed together. This also suggests that the intrinsic cloud relaxation timescales are much faster than . Furthermore, due to the presence of the Mott plateau, we find it insightful to describe the cloud as consisting of three domains, naturally defined in terms of the density: a low-density () and high-density () metallic region and the Mott plateau (). For each of them, we consider a distinct and the related expansion speed. For , we find a positive speed (for , the low-density region is the same as the original in Fig. 2d). Nevertheless, the speed is negative for the other domains, which means that their outer radius is shrinking in time (note, however, that the Mott plateau is growing inwards). Even after the high-density core has disappeared, the Mott region has only just started to shrink, thus confirming that the different domains have qualitative different behaviours.

Finally, in Fig. 2f we present results for the Kohn-Sham bond-currents at time . Although their physical meaning is rather indirect, such currents can offer additional insight (Methods). From their intensity and direction, we see that the particles are flowing out of the high-density region in the centre. However, no particles accumulate in the Mott region, consistently with the rigidity of the latter.

The foregoing discussion for a clean lattice reveals several interesting dynamical features of strongly correlated fermions. Another important element that we wish to bring into our analysis is the effect of disorder, which we discuss next.

Expansion into an inhomogeneous lattice. To illustrate the role of disorder, we consider one particularly interesting situation, namely the dynamical competition between disorder and interaction in the Mott regime. To this end, we prepared the system for in the same initial state as in the previous section and let it expand into a disordered environment according to the slow trap-opening protocol used before. Disorder was introduced via a single special quasi-random structure (Methods), which mimics a 50% binary alloy ( in Eq. 1). The :s, as shown in Fig. 3b, are chosen to be non-zero only where the initial density is negligible, and thus we can specifically address the role of disorder on the cloud dynamics.

In Fig 3a, we show the density profile along the plane for the homogeneous () and the disordered case . For this setup, the effect of disorder on the density profile is hardly visible on the scale of the figure for , due to the peculiar Mott plateau dynamics. This is easily understood by noting that, initially, the expansion is mostly characterised by an internal rearrangement of the high-density region and the formation of the extended Mott plateau, as in the clean case. This has a relatively small effect on the low-density region.

The influence of disorder is instead much more pronounced for when the Mott plateau begins to collapse and the particles around its boundaries flow towards the low-density region. A first more evident effect is that the cloud expansion at large times is now hindered by the scattering against the binary disorder, which causes an irregular density accumulation close to the inner edge of the disordered region. At large times , the disorder induces a kind of dynamical localisation of the cloud as the density profile is significantly reduced in the trap centre (in fact more than in the case) but at the same time the expansion far away from the centre is considerably slowed down. The final density profile at shows a large particle accumulation in a roughly annular region at the beginning of the disordered zone, in contrast to the clean case, which is almost uniform (Fig. 3a). A second and more subtle effect, hardly visible on the plot scale and discussed further below, is that the disorder accelerates the melting of the Mott plateau.

To quantify the effect of disorder on the size and symmetry of the expanding cloud, we also analysed both the minimum and maximum radial size of the cloud with and without disorder. Not unexpectedly, in presence of disorder and rapidly diverge from each other, since disorder destroys the spherical symmetry observed in the clean case. separates from its clean-setup counterpart in the first stages of the cloud expansion, indicating that disorder earlier favours some particles to flow away from the cloud in a faster way. Interestingly, disorder induces a larger and more asymmetric cloud at every step of the expansion.

The inverse participation ratio , shown in Fig. 3d for and as a function of time, is a convenient indicator of the degree of localisation resulting from competing disorder and interactions (Methods). Since a smaller implies a smaller degree of localisation for the particle cloud, in the first stage of the expansion () we observe a rather interesting and counterintuitive result: for , stays smaller than for the homogeneous case (whereas disorder in general is expected to increase the degree of localisation), resulting in a more delocalised cloud at small expansion times. The dynamical localisation in Fig. 3a takes over at large times and the trend is inverted for .

This dynamical crossover is an intriguing consequence of the competition between disorder and interaction in the Mott regime. Indeed the disordered lattice offers in the early stages of the expansion additional energetic pathways compared to the case, which accelerate the melting of the Mott plateau. This is also why the expansion initially occurs earlier and faster than in the clean case (thus we observe a smaller and a larger ). However, after the Mott plateau has melted, the cloud expands faster (on average) in the clean setup, since the particle flow is not hindered by disorder (hence the larger at for the disordered case). We verified that the overall features are robust against distinct quasi-random structures.

## Discussion

We have studied the dynamics of an interacting fermionic cloud expanding into clean and disordered optical lattices by using a recently developed approach which merges lattice Time-Dependent Density-Functional Theory with Dynamical Mean-Field Theory. With this method, we have been able to describe the real-time dynamics of large three-dimensional (in)homogeneous Fermi systems up to an unprecedented size ( lattice sites, i.e. comparable with current experiments in cold gases), in different regimes of interaction and even in presence of disorder. Our work unveils important aspects of the fermionic dynamics, such as the central role of the Mott physics at strong-coupling and the interplay between disorder and interaction in 3D.

Above the Mott threshold, the timescales and features of the clean expansion are markedly different from the metallic regime. We observe an earlier increase in size of the Mott plateau against the trap opening at the expense of the metallic domain. Compared to lower dimensions, this finding is even more surprising, due to the larger number of runaway paths in 3D: this mainly arises from the universal features of Mott physics. The 3D nature of the system is instead crucial i) in the observed dichotomy between the weak- and strong-coupling dynamics, which reflects the presence of a finite Mott threshold in the homogeneous system. Our findings also suggest a convenient description of the cloud expansion in terms of multiple domains with positive and negative expansion speeds ii) in the rich phenomenology observed in presence of disorder.

Disorder introduces notable changes in the dynamics above the Mott threshold, altering expansion time-scales, but also resulting in interesting temporal patterns. For example, we observe a dynamic crossover in the localisation properties, as the disordered system is less localised than the clean one at the beginning of the expansion, whilst disorder increases in general localisation, as expected, at large times.

Our results shown here are just an example of how rich and interesting the real-time dynamics in the presence of disorder and interactions is. We have already considered (not shown) different disorder strengths and interactions, and observed intriguing features like, e.g., density revivals due to dynamical localisation. These findings, as well as the inclusion of proper configuration averaging and disorder in the initial state, are deferred to future publications.

In conclusion, we hope our work will stimulate further improvements of our approach, as well as new experiments on ultracold gases in optical lattices, to deepen the understanding of many-body quantum systems in many regimes of interest.

## methods

To describe the properties of the inhomogeneous Hubbard model in Eq. (1), we use static [17, 18] and Time-Dependent [19] Density-Functional Theory (DFT and TDDFT, respectively). These are in-principle-exact reformulations of the many-body problem, and we use them in their (static [20, 21, 22, 23, 24] and time-dependent [25, 26, 27, 28, 29, 30, 31, 16, 9, 32, 33]) lattice versions.

Here, we briefly recall the essentials of our treatment (for a recent review, see Ref. [32]). In this approach, the (time-dependent) number of particles per lattice site , is the basic variable and the physical observables of the system are functionals of . In operational terms, one introduces a non-interacting image system, the so-called Kohn-Sham (KS) system, and the exact many-body density is then obtained from the KS single-particle states. A key ingredient in the KS system is the exchange-correlation potential, , incorporating exchange and correlation effects. In general, is not known exactly and approximations are used. A simple but effective one (used here, and further discussed below) is the Local Density Approximation (LDA) for the static case, where at site depends locally on the site occupation , and correspondingly the adiabatic LDA (ALDA) for the time-dependent case, where depends instantaneously and locally on .

In a recent work, a lattice DFT treatment of simple-cubic Hubbard model[9] was proposed, where the pivotal ingredient is an adiabatic LDA based on Dynamical Mean-Field Theory (DMFT) [34, 35, 36]. There, for the reference homogeneous system was obtained within DMFT according to

(2) |

where is the ground-state energy, is the non-interacting kinetic energy and is the Hartree energy. DMFT properly describes the Mott metal-insulator transition [35], and gives a good variational estimate of the energy [37], although the self-energy only depends on the frequency and not on quasi-momentum.

A crucial feature of the method is the occurrence of a discontinuity in at above a critical value of the Hubbard interaction. This is how the Mott-Hubbard metal-insulator transition manifests within a DFT framework. The discontinuity is the origin of the Mott plateaus in Fig.1, for .

In our calculations, we considered simple-cubic clusters of lattice sites with open boundary conditions. We chose particles when to avoid ground state degeneracies in the density region of interest. The ground state was computed by solving self-consistently the KS equations:

(3) |

where the effective potential ( labels the lattice site) , with the kinetic energy operator on the lattice, and the Hartree contribution, with (the sum is over all occupied KS orbitals, and the factor 2 accounts for spin degeneracy). is the external trapping potential. We chose for . In the LDA, . The large-scale self-consistent eigenvalue problem of Eq. (3) was made computationally manageable by using symmetry-adapted orbitals via the point group . Also, for , was slightly smoothened for numerical convenience (see e.g. ref. [16]). As a result, the Mott domain is not perfectly flat. The threshold for self-consistency in the density was .

On removal of the trap, the system was evolved in time via the time-dependent KS equations

(4) |

where , with (in the ALDA) , and . The time-dependent external potential is . When present, the disorder is static, whilst the trap is removed according to . We chose for and when ensuring a smooth time dependence. The numerical time propagation was performed via the short iterated Lanczos algorithm [38], with a time-step .

The (A)LDA neglects memory effects and the (dynamical) broadening of the discontinuity in for inhomogeneous systems (in general, these limitations become more severe at lower dimensions). As shown by recent benchmarks ,[27, 29, 39, 32, 16, 40], where non-locality in time and space is fully taken into account, the (A)LDA performs poorly for strong interactions (e.g. in gapped Mott systems) when the perturbations are fast. However, the same benchmarks have shown that, for fast fields and not very strong interactions or slow fields irrespective of the interactions, these shortcomings appear to be much less important, and ALDA can give reliable results.

To simulate a disordered 50% binary alloy, we used a single disorder configuration chosen according to the special quasi-random structure approach, to effectively describe the random arrangements of sites at short range [41]. This still provides significant insight in localisation trends, while avoiding demanding numerical averaging over many configurations [42]. As an indicator of de/localisation, we used a modified inverse participation ratio :

(5) |

exactly accessible within lattice (TD)DFT.

Finally, we examined the KS currents. In TDDFT, the meaning of the KS current-density is rather indirect. It is in fact its divergence which, via the continuity equation, can be exactly determined (and thus the exact current out a region can be obtained). The same holds in lattice TDDFT: rigorous physical content should not be assigned to the KS bond currents

(6) |

but rather to their divergence on the lattice ( label the nearest neighbour sites defining the bond ). Even so, it may be instructive to consider as an auxiliary, albeit non-rigorous tool to illustrate the dynamics.

## References

- [1] Fausti, D. et al. Light-Induced Superconductivity in a Stripe-Ordered Cuprate Science 331, 6014 189–191 (2011)
- [2] Dal Conte, S.et al. Disentangling the Electronic and Phononic Glue in a High-Tc Superconductor Science 335, 6076 1600–1603 (2011)
- [3] Jördens R., Strohmaier N., Günter K., Moritz H., Esslinger T. A Mott insulator of fermionic atoms in an optical lattice. Nature Physics 455, 204-207 (2008)
- [4] Khatami, E., Rigol, M Accessing the Mott Regime in 2D Optical Lattices with Strongly Interacting Fermions. J. Supercond. Nov. Magn. 25, 2145 (2012)
- [5] Bloch I., Dalibard J., Zwerger W. Many-body physics with ultracold gases. Rev. Mod. Phys. 80, 885 (2008)
- [6] Bloch, I., Dalibard, J., Nascimbene, S. Quantum simulations with ultracold quantum gases. Nature Physics 8, 267-276 (2012)
- [7] Jreissaty M., Carrasquilla J., F. A., Rigol M. Expansion of Bose-Hubbard Mott insulators in optical lattices. Phys. Rev. A 84, 043610 (2011)
- [8] Schneider, U. et al. Fermionic transport and out-of-equilibrium dynamics in a homogeneous Hubbard model with ultracold atoms. Nature Physics 8, 213-218 (2012)
- [9] Karlsson, D., Privitera, A., Verdozzi, C. Time-Dependent Density-Functional Theory Meets Dynamical Mean-Field Theory: Real-Time Dynamics for the 3D Hubbard Model. Phys. Rev. Lett. 106, 116401, 213-218 (2012)
- [10] Jaksch D., Bruder C., Cirac J. I., Gardiner C. W., Zoller P. Cold Bosonic Atoms in Optical Lattices. Phys. Rev. Lett. 81, 3108 (1998)
- [11] Abrahams E. 50 Years of Anderson Localization. World Scientific, Singapore (2010)
- [12] Guarrera V., Fallani L., Lye J. E., Fort C., Inguscio M. Inhomogeneous broadening of a Mott insulator spectrum . New J. Phys. 9, 107 (2007)
- [13] Gao, Xianlong; Polini, M.; Tosi, M. P.; Tanatar, B. Effect of Disorder on the Interacting Fermi Gases in a One-Dimensional Optical Lattice. Int. J. Mod. Phys. B 22, 4500–4510 (2008).
- [14] Gao, Xianlong Ground-state phases of interacting Fermi gases in disordered one-dimensional lattices. J. Phys. B: At. Mol. Opt. Phys. 45, 225304 (2012).
- [15] Jendrzejewski, F., et al. Three-dimensional localization of ultracold atoms in an optical disordered potential. Nature Physics 8, 398-403 (2012)
- [16] Karlsson D., Verdozzi C., Odashima M. M., Capelle K. Dynamical self-stabilisation of the Mott insulator: Time evolution of the density and entanglement entropy of out-of-equilibrium cold fermion gases. EPL 93, 23003 (2011)
- [17] Hohenberg P., Kohn W. Inhomogeneous Electron Gas. Phys. Rev. 136, 864 (1964)
- [18] Kohn W., Sham L.J Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 140, 1133 (1965)
- [19] Runge E., Gross E.K.U Density-Functional Theory for Time-Dependent Systems. Phys. Rev. Lett. 52, 997 (1984)
- [20] Gunnarsson O., Schönhammer K. Density-Functional Treatment of an Exactly Solvable Semiconductor Model. Phys. Rev. Lett. 56, 1968 (1986)
- [21] Schönhammer K., Gunnarsson O., Noack R. M. Density-functional theory on a lattice: Comparison with exact numerical results for a model with strongly correlated electrons. Phys. Rev. Lett. 52, 2504–2510 (1995)
- [22] Lima N. A., Silva M. F., Oliveira L. N. , Capelle K. Density Functionals Not Based on the Electron Gas: Local-Density Approximation for a Luttinger Liquid. Phys. Rev. Lett. 90, 146402 (2003)
- [23] Xianlong G., et al. Bethe ansatz density-functional theory of ultracold repulsive fermions in one-dimensional optical lattices. Phys. Rev. B. 73, 165120 (2006)
- [24] Akande A., Sanvito S. Electric field response of strongly correlated one-dimensional metals: A Bethe ansatz density functional theory study. Phys. Rev. B. 82, 245114 (2010)
- [25] Aryasetiawan F., Gunnarsson O., Rubio A. Excitation energies from time-dependent density-functional formalism for small systems. Europhys. Lett. 57, 683 (2002)
- [26] Magyar R. J. Ground and excited-state fermions in a one-dimensional double-well: Exact and density-functional solutions. Phys. Rev. B. 79, 195127 (2009)
- [27] Verdozzi C. Time-Dependent Density-Functional Theory and Strongly Correlated Systems: Insight from Numerical Studies. Phys. Rev. Lett. 101, 166401 (2008)
- [28] Baer R. On the mapping of time-dependent densities onto potentials in quantum mechanics. J. Chem. Phys. 128, 044103 (2008)
- [29] Li W., Xianlong G., Kollath C., Polini M. Collective excitations in one-dimensional ultracold Fermi gases: Comparative study. Phys. Rev. B 78, 195109 (2008)
- [30] Stefanucci G., Perfetto E., Cini M. Time-dependent quantum transport with superconducting leads: a discrete basis Kohn-Sham formulation and propagation scheme Phys. Rev. B 81, 115446 (2010)
- [31] Stefanucci G., Kurth, S. Towards a Description of the Kondo Effect Using Time-Dependent Density-Functional Theory. Phys. Rev. Lett 107, 216401 (2011)
- [32] Verdozzi C., Karlsson D., Puig von Friesen M., Almbladh C.-O., von Barth U. Some open questions in TDDFT: Clues from lattice models and Kadanoffâ-Baym dynamics. Chem. Phys. 391, 1, 37-49 (2011)
- [33] Farzanehpour M., Tokatly I. V. Time-dependent current density functional theory on a lattice. Phys. Rev. B 86, 125130 (2012)
- [34] Metzner W., Vollhardt D. Correlated Lattice Fermions in Dimensions Phys. Rev. Lett 62, 324 (1989)
- [35] Georges A., Kotliar, G., J. Rozenberg M., Krauth, W. Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions. Rev. Mod. Phys. 68, 13 (1996)
- [36] Kotliar G., Vollhardt D. Strongly Correlated Materials: Insights from Dynamical Mean-Field Theory. Physics Today 57, 3, 53-59 (2004)
- [37] Kozik E., et al. Diagrammatic Monte Carlo for correlated fermions. EPL 90, 10004 (2010)
- [38] Park T. J., Light J. C. Unitary quantum time evolution by iterative Lanczos reduction. J. Chem. Phys. 85, 5870 (1986)
- [39] Puig von Friesen M., Verdozzi C., C.-O. Almbladh, C.-O. Kadanoff-Baym dynamics of Hubbard clusters: Performance of many-body schemes, correlation-induced damping and multiple steady and quasi-steady states. Phys. Rev. B. 82, 155108 (2010)
- [40] Uimonen A.-M., et al. Comparative study of many-body perturbation theory and time-dependent density functional theory in the out-of-equilibrium Anderson model. Phys. Rev. B 84, 115103 (2011)
- [41] Zunger A., Wei S.-H., Ferreira L. G., Bernard J. E. Special quasirandom structures. Phys. Rev. Lett. 65, 353 (1990)
- [42] Verdozzi C., Durham P. J., Cole R. J., Weightman P. Correlation and disorder effects in photoelectron and Auger spectra: The late transition metals and their alloys. Phys. Rev. B. 55, 16143 (1997)

## addendum

Acknowledgments. We thank C.-O. Almbladh and J. Hein for valuable discussions.
We gratefully acknowledge LUNARC for computational resources and technical support.
A.K., D.K. and C.V. thank the EOARD (grant FA8655-08-1-3019) and the ETSF (INFRA-2007-211956) for
financial support. A.P. thanks M. Capone for the financial support by the European Research Council
under FP7/ERC Starting Independent Research Grant ÒSUPERBADÓ (Grant Agreement n. 240524).

Author contributions. A.K. developed the computer codes to perform the numerical simulations,
with contributions from C.V. and D.K. The numerical calculations were
performed by A.K. The data analysis was performed by A.P. and D.K., with
contributions by C.V. All the authors participated to the discussion of the results and to the
writing of the paper. The project was supervised by C.V.

Competing interests statement. The authors declare that they have no competing financial interests.

## The Authors

Alexey Kartsev,
Mathematical Physics and European Theoretical Spectroscopy Facility (ETSF), Lund University, 22100 Lund, Sweden

Daniel Karlsson, Mathematical Physics and European Theoretical Spectroscopy Facility (ETSF), Lund University, 22100 Lund, Sweden

Antonio Privitera, Democritos National Simulation Center, Consiglio Nazionale delle Ricerche, Istituto Officina dei Materiali (IOM) and Scuola Internazionale Superiore di Studi Avanzati (SISSA), Via Bonomea 265, 34136 Trieste, Italy

Claudio Verdozzi, Mathematical Physics and European Theoretical Spectroscopy Facility (ETSF), Lund University, 22100 Lund, Sweden (Corresponding author: cv@teorfys.lu.se)