Effective doublon and hole temperatures in the photo-doped dynamic Hubbard model.

Hirsch's dynamic Hubbard model describes the effect of orbital expansion with occupancy by coupling the doublon operator to an auxiliary boson. In the Mott insulating phase, empty sites (holes) and doubly occupied orbitals (doublons) become charge carriers on top of the half-filled background. We use the nonequilibrium dynamical mean field method to study the properties of photo-doped doublons and holes in this model in the strongly correlated regime. In particular, we discuss how photodoping leads to doublon and hole populations with different effective temperatures, and we analyze the relaxation behavior as a function of the boson coupling and boson energy. In the polaronic regime, the nontrivial energy exchange between doublons, holes, and bosons can result in a negative temperature distribution for the holes.

INTRODUCTION

Strongly correlated electron systems are often described in terms of the Hubbard model or its multiorbital extensions. This model assumes that the electrons interact only locally on each lattice site, and that this interaction is of a particle-hole symmetric form. In the case of the single-band Hubbard model, the interaction on site i, parametrized by U, can be written as , so that the energy costs for adding and removing an electron in the half-filled state are the same. As was pointed out by Hirsch and co-workers in a series of papers, this model ignores the fact that doubly occupied atomic orbitals expand as a result of Coulomb interactions and electron correlations. To capture this physics, Hirsch introduced the dynamic Hubbard model, which contains an additional coupling of the doublon operator to an auxiliary boson. The expansion of the orbital and associated reduction of the Coulomb interaction is then described by the dynamics of the boson and the resulting screening effect. Several variants of this model have subsequently been formulated and explored. These equilibrium studies have revealed interesting properties resulting from the filling-dependent particle-hole asymmetric correlation effects, which may even be relevant for the understanding of unconventional superconductivity.

More recently, a dynamic Hubbard model has been proposed to explain the behavior of organic charge-transfer salts in which molecular orbitals are periodically modulated via the selective excitation of an optical phonon mode with intense laser fields. In this context, the particle-hole asymmetric electron-boson coupling provides an experimental “knob” to tune the interaction effects in the solid, and hence raises the interesting prospect of dynamical control of material properties.

Theoretically, the nonequilibrium properties of dynamic Hubbard models are largely unexplored, but as in equilibrium, one can expect interesting effects resulting from the particle-hole asymmetric correlations. Examples are the nontrivial interplay between heat and particle currents (Seebeck effect), or the unusual nature of the photo-doped insulator with its coexisting weakly/strongly correlated doublons/holes. Motivated by these considerations, we investigate here some nonequilibrium phenomena in the dynamic Hubbard model by means of the nonequilibrium dynamical mean field theory (DMFT) formalism. In particular, we will study how the dynamic Coulomb correlations are reflected in the relaxation and recombination of doublons and holes produced by an impulsive stimulation of the Mott insulating state. We will show that the particle-hole asymmetry leads to different transient temperatures for the doublons and holes during the relaxation, and that the respective cooling rates depend in a nontrivial way on the coupling strength of the auxiliary boson.

The rest of this paper is organized as follows. Section II describes the model considered in this study and the approximate impurity solver used in the DMFT calculations. Section III presents results for the nonequilibrium spectral function and the time-evolution of the effective doublon and hole temperatures, while Sec. IV provides a short summary.

MODEL AND METHOD

We consider the dynamic Hubbard model introduced in the original paper by Hirsch where v denotes the hopping amplitude, U is the onsite repulsion, μ is the chemical potential, and measures the occupation of site i with electrons of spin σ. The coupling of the double occupancy to a harmonic oscillator with frequency ω0 allows to describe the expansion of the atomic orbital (and associated reduction of the Coulomb interaction) which occurs when the orbital is occupied by two electrons. X and P are the boson position and momentum operators, respectively, and g is the electron-boson coupling. (Here, the expansion of the orbital corresponds to a shift of X in the negative direction.) In the DMFT approximation, this lattice model is mapped to a self-consistent solution of a quantum impurity model with Hamiltonian and a local term where measures the double occupation. The hybridization term and the noninteracting electron bath are determined by a self-consistent procedure in such a way that the bath mimics the lattice environment. In the action representation of the impurity model, the coupling to the environment appears as a hybridization term , and the hybridization function is determined by the DMFT self-consistency.

To solve this impurity model, we perform a unitary transformation which transforms operators , and thus shifts . The transformed local Hamiltonian reads where we denote transformed operators by a tilde. The electrons and bosons can be decoupled by choosing , so that the unitary transformation becomes with and the terms involving X0 in Eq. (3) evaluate to . The transformed interaction and chemical potential therefore read

We also have to evaluate the transformed electron creation and annihilation operators, . The creation operator gives a nonzero result only if it acts on a state with . In this case, the operator W on the right becomes 1, and the operator on the left becomes . Hence, we find , and similarly . Writing in terms of the boson creation and annihilation operators and separating the contributions from the occupied and empty state, we finally get where we have introduced the Hubbard operators and . Using Eqs. (6) and (7), the hybridization term of the impurity action can be written as . (The remaining terms involving products and are not written explicitly because they do not contribute to the lowest order expansion considered below.)

As introduced for the Holstein-Hubbard model in Ref. 12, an approximation that captures the essential physics in the insulating phase and in the metal-insulator crossover regime at not too low temperature is to partially resum the diagrams of a hybridization expansion (non-crossing approximation, NCA), and to decouple the electron and boson terms by factorizing products in the action and in the Green function, with the bosonic function ( is the contour ordering operator) Here, denotes the earlier/later time on the Kadanoff-Baym contour. The same electron-boson decoupling has recently been shown to provide the optimal representation of an Anderson impurity model with dynamically screened interactions in terms of a model with static interactions. In the case of the dynamic Hubbard model, we follow the same procedure and write Similarly, the contour-ordered Green function is factorized within the NCA like The correlation functions and are then evaluated using the standard NCA diagrammatic expressions, as illustrated in Fig. 1.

FIG. 1.

Illustration of the NCA self-energy diagrams for the dynamic Hubbard model. An open circle corresponds to an operator , a full circle to , an open square to , and a full square to . is the hybridization function, is a pseudo-particle Green's function, and B is the bosonic propagator defined in Eq. (8). Red squares represent the “box weight” (14) in a simulation with externally driven bosons.

Because the dynamic Hubbard model is not particle-hole symmetric, an external driving of the bosons can be expected to yield nontrivial effects. Such a driving might have to be considered in the modeling of experiments similar to those recently reported by Kaiser et al., or Singla et al., where the lattice and hence the Wannier orbitals are shaken. This effect can be described within the present formalism, starting from the model with f(t) some time-dependent function. Replacing in Eq. (3), the Lang-Firsov shift becomes , and the shifts of the interaction and chemical potential are obtained from . Hence, , and there is an additional time-dependent energy shift .

According to the discussion in Section II B of Ref. 12, the time-dependent unitary transformation W(t) implies an additional term in the Hamiltonian. Including this term, the bosonic part of the transformed Hamiltonian becomes From Eq. (6) and the solution of the Heisenberg equations of motion for the b-operators, , one finds (assuming ) and similarly for (t). In the hybridization expansion, all vertices and are thus multiplied by the “box weight” with for creation and annihilation operators (see Fig. 1). If the driving is of the form , then for , and for .

RESULTS

“Photo-doping” of the Mott insulator

Since the approximate impurity solver described in Section II is expected to give qualitatively correct results in the insulating phases of the model, we will consider a “photo-doping” type perturbation in a paramagnetic Mott insulator. For simplicity, we do not simulate the excitation of carriers across the Mott-gap by an electric field pulse, but instead apply a short U-pulse to the system: For , the effect of this modulation is to promptly create a certain number of doublon-hole pairs, so that the state of the system after the U-pulse is similar to a photo-doped state after impulsive stimulation with a laser pulse with frequency comparable to U. The injected doublons and holons (or simply “holes”) are long-lived mobile charge carriers which exist on the background of the half-filled Mott insulator. The spectral weight associated with these carriers appears in the upper and lower Hubbard bands, respectively. If U is much larger than the kinetic energy, doublon-hole recombination processes are suppressed, and it then becomes meaningful to study the energy distribution and effective temperatures of these long-lived carriers. In the following calculations, we assume a semi-circular density of states of bandwidth 4, so that the DMFT self-consistency loop simplifies to , and we use as the unit of energy, and as the unit of time.

The photo-doping induced changes in the spectral function of the (particle-hole symmetric) Holstein-Hubbard model have been investigated in Ref. 18. In contrast to the Hubbard model, where the Mott insulating spectral function is little affected by photo-doping, in particular, in the gap region, the rapid cooling of the photo-doped carriers and the possibility of polaron formation in the Holstein-Hubbard model with a strong electron-phonon coupling can lead to a series of doping-induced in-gap states. The appearance of these polaronic states also has an important effect on the life-time of the photo-carriers, which becomes much shorter than in the Hubbard model, and also doping dependent (, with the doublon population of the initial state). The energy dissipation rate of photo-doped carriers in the Holstein-, Holstein-Hubbard, and Peierls-Holstein-Hubbard model has been studied using exact diagonalization, density matrix renormalization group (DMRG), and DMFT techniques. In the case of the Mott insulating Holstein-Hubbard model, it was found that the cooling rate is controlled by the effective boson coupling .

Equilibrium results

We now turn to the dynamic Hubbard model (1) and first briefly discuss the spectral properties of the Mott insulating equilibrium system. The left panel of Fig. 2 plots the spectral function for U = 10, μ = 0, and indicated values of g. The inverse temperature is β = 5. Because the chemical potential is inside a large gap, the filling is approximately n = 1 independent of g. For g = 0 (Hubbard model), the system is particle-hole symmetric and the spectral function features two Hubbard bands of approximately semi-circular in shape. As the electron-boson coupling is increased, the spectral function becomes asymmetric: while the lower Hubbard band essentially retains its semi-circular shape, the upper Hubbard band broadens and for large enough g splits into sidebands. In addition, there is a small shift of the Hubbard bands to higher ω. If the parameter is kept fixed while ω0 is varied, the separation between the sidebands changes, while the broadening of the upper Hubbard band stays roughly the same (middle panel). The boson-induced broadening of the upper Hubbard band in the small-g regime is a manifestation of the reduced doublon correlations, which is the effect described in the original work by Hirsch. The splitting into sidebands at larger g is due to the formation of polarons. We will not discuss here under which conditions this polaronic regime may be realized, but simply explore the nonequilibrium properties of the dynamic Hubbard model in the full range of boson couplings.

FIG. 2.

Equilibrium spectral functions and compressibilities in the dynamic Hubbard model. The left panel shows spectral functions for U = 10, μ = 0, , and indicated values of g, while the middle panel shows spectral functions for U = 10, μ = 0, fixed , and indicated values of ω0. The right panel plots the compressibility as a function of μ for g = 1, , U = 5, and U = 10. The inverse temperature is β = 5 for all plots.

One can also observe the asymmetry in the electron and hole correlations by computing the density or the compressibility . As shown in the right panel of Fig. 2, the compressibility is larger on the electron-doped side than on the hole-doped side, which is consistent with weaker doublon correlations.

Spectral functions of the photo-doped system

We next investigate the changes in the spectral function after the impulsive excitation defined in Eq. (15). Figure 3 shows results for U = 10, g = 1, , and pulse parameters and . The top left panel corresponds to a half-filled system. We plot both the total spectral function , the distribution of occupied states , and the distribution of unoccupied states . These time-dependent spectra were obtained by forward time integration, using the formulas We plot the results for t = 15, i.e., during the relaxation following the short 2-cycle pulse, which lasts until . As can be seen by comparing the transient spectral function to the initial equilibrium result in grey, the photo-doping induces substantial changes, especially in the gap region. In analogy to the previous results for the Holstein-Hubbard model we find that the production of doublons results in a series of in-gap states, which can be viewed as additional boson sidebands of the upper Hubbard band. However, in contrast to the particle-hole symmetric Holstein-Hubbard case, the sidebands are limited to the region , and there are no in-gap states appearing in the vicinity of the lower Hubbard band. The lower Hubbard band merely exhibits an asymmetry reminiscent of a smeared-out quasi-particle band near the upper band edge.

FIG. 3.

Left and bottom panels: Nonequilibrium spectral functions measured at time t = 15 in a half-filled, hole-doped, and electron-doped system after photo-doping. The parameters are U = 10, g = 1, , and β = 5. The frequency of the two-cycle excitation pulse Eq. (15) is and the amplitude is . The black curve shows the total spectral function, while the green and pink curves indicate the distributions of empty and occupied states, respectively. For comparison, we also plot the initial equilibrium spectral functions. The top right panel illustrates the time evolution of the distribution functions for n = 1.

The distribution of occupied states shows that the additional doublons created by the pulse are rapidly cooled and are the origin of the polaronic in-gap states. At the same time, the distribution of unoccupied states shows an accumulation of holes near the upper edge of the lower Hubbard band. The time-evolution of the doublon and hole populations for n = 1 is illustrated in the top right panel. The curves with colors ranging from green to blue show the function in the lower Hubbard band from t = 3 to 21, and the set of curves with colors ranging from red to dark-violet show the function in the upper Hubbard band for the same times. One notices that the doublon population exhibits almost no time dependence, which shows that the cooling of the doublons and the appearance and population of the polaronic in-gap states is so fast that it essentially takes place during the application of the pulse. This observation is consistent with the Holstein-Hubbard results in Ref. 18 for the same parameter regime. The redistribution of the hole population is a much slower process, which is expected because the holes do not couple directly to the bosons in the dynamic Hubbard model. Nevertheless, the cooling of the holes which is evident in this set of data is very fast compared to the relaxation in the Hubbard model (g = 0), where there would be almost no change in the energy distribution of doublons and holes on the timescales of this plot. Apparently, the scattering with “cold” doublons provides a relatively efficient dissipation mechanism for the holes. In Section III D, we will analyze this relaxation process more quantitatively.

We can repeat similar calculations also for hole and electron doped systems. In the bottom panels of Fig. 3, we show results for n = 0.85 and n = 1.15. In the hole-doped system, the spectral function features a narrow quasi-particle peak at the Fermi energy. After the pulse, this peak broadens due to the effects of heating and photo-doping. Near the upper Hubbard band, we again observe the appearance of a series of polaronic side-bands. They are less prominent than in the half-filled case, because in the doped system, the doublon-hole recombination is faster, so that the number of photo-carriers in the upper Hubbard band is substantially reduced by t = 15. In the electron-doped system, we observe a series of doping-induced side-bands already in the equilibrium state. After the pulse, the weight of these polaronic states is further increased, while in the lower Hubbard band, we observe the appearance of a small “quasi-particle” peak. Again, because of the faster doublon-hole recombination, the changes in the lower Hubbard band are smaller than in the half-filled case.

Effective doublon and hole temperatures

Except deep in the polaronic regime, the doublon and hole distribution functions quickly reach an approximately thermal shape in the energy region corresponding to the upper and lower Hubbard band. It is therefore instructive to analyze the relaxation process in terms of the effective “temperatures” of the doublons and holes. To measure these temperatures, we compute the nonequilibrium distribution function . In equilibrium, corresponds to the Fermi function for the corresponding temperature. As shown in the top left panel of Fig. 4 for the model with g = 0.5 and , the photo-doped system does not have a Fermi like nonequilibrium distribution. Rather, there are two Fermi edges, one in the energy range of the upper Hubbard band ( in the figure), and the other in the energy range of the lower Hubbard band (). By fitting the distribution functions to a Fermi function in these energy regions, one can extract an effective doublon and hole temperature, and in principle also the corresponding effective chemical potentials. For the earliest time (t = 12), these fits are indicated in the figure by dashed black lines. From these fits, one extracts and .

FIG. 4.

Top left panel: Nonequilibrium distribution function for the model with g = 0.5, , and indicated values of t. To illustrate the energy range covered by the Hubbard bands, we also plot (black line). For the t = 12 curve, we show the fits to Fermi distribution functions by dashed lines (fitting range and ). Top right panel: Effective electron and hole temperatures extracted from the Fermi function fits for g = 0.5, . The dashed line indicates the temperature of the initial state, and the arrow the end of the 2-cycle pulse. Bottom panels: Effective doublon and hole temperatures for different boson couplings and boson frequencies.

The time evolution of these effective doublon and hole temperatures is illustrated in the top right panel. Since the Fermi fits depend on the fitting range, there is some uncertainty associated with the calculation of the effective temperatures, and the error bars in the figure give a conservative estimate of this uncertainty. It is obvious that the doublons, which couple directly to the bosons, are rapidly cooled down to an effective temperature close to the initial equilibrium temperature of (horizontal dashed line), while the holes are still in the process of dissipating their energy at the longest accessible time. We note that in our closed system, the energy injected by the pulse will lead to a thermalization at a temperature above the initial equilibrium temperature, so that the effective temperature of the doublons might increase on longer timescales. However, since the energy changes associated with the redistribution of spectral weight within the Hubbard bands and the cooling of doublons and holes during the initial stage of the relaxation can apparently be easily absorbed by the bosons, such an increase of the doublon temperature is controlled by the doublon-hole recombination time, which is very long in large-gap insulators.

The bottom panels show similar results for different boson couplings and boson frequencies. Obviously, the cooling rates and the ratio between effective doublon and hole temperatures depend on these parameters in a nontrivial way.

Dependence of the cooling rate on the electron-electron and electron-boson interaction

While the coupling to the bosons provides a dissipation mechanism for the doublons, the cooling of the holes is an indirect process which involves doublon-hole scattering. Here, one could imagine scattering processes with virtual break-up of the doublon, , or the scattering between a high energy hole and a low-energy doublon, which transfers energy from the hole to the doublon without such virtual excitations. If the former processes are relevant, one should observe a dependence of the cooling rate on U, because the local energies of the states d + h and differ by U. Figure 5 shows the U-dependence of the effective doublon and hole temperatures for different boson couplings and boson energies. There is a small dependence of the effective hole temperature on U, which suggests a more rapid cooling in small-gap insulators, and thus a certain role of the first type of scattering processes in the initial stage of the relaxation. The effective temperatures at later times are however essentially independent of U. It thus appears that the dominant scattering mechanism is of the second type, i.e., energy transfer from hot holes to cold doublons without virtual excitations.

FIG. 5.

Effective doublon and hole temperatures in systems with photo-doping concentration 0.1, indicated boson couplings and boson energies, and different values of U. The upper curves show the hole temperatures and the lower curves the doublon temperatures.

We also analyzed the dependence on the parameters g and ω0 and found that the cooling rate of the doublons and holes is essentially determined by . This is illustrated in the left panel of Fig. 6, which shows the effective doublon and hole temperatures for U = 15, photo-doping concentration 0.1, and different parameter choices corresponding to . The rate with which the effective temperature decreases is similar in the three cases, even though the boson frequencies and coupling strengths are very different, and the slowest boson oscillation period () is longer than the plotted time range. The vertical offsets of the curves indicate that the effective temperatures reached in the limit (with the cooling time and the doublon-hole recombination time) depend on g and ω0, not only on the effective coupling λ. We also find that, in agreement with the Holstein-Hubbard results of Ref. 18, the data collapse for becomes worse for larger couplings and larger boson frequencies.

FIG. 6.

Effective temperatures in systems with U = 15 and photo-doping concentration 0.1. Left panel: constant effective coupling strength . Middle panel: constant boson frequency . The upper curves show the hole temperatures and the lower curves the doublon temperatures. Right panel: inverse effective temperatures for doublons and holes as a function of g. The dashed (full) curves show the result at t = 20 (t = 40). Note that for , the definition of the doublon temperature becomes ambiguous due to deviations of from a Fermi function. The uncertainty in is , and thus larger than the difference between the t = 20 and t = 40 values.

For fixed , the cooling rate of the holes depends in a non-monotonous way on the coupling strength g (see middle panel). If g is small, the energy dissipation to the bosons becomes slow. As a result, the doublon temperature decreases slowly and the effective hole and doublon temperatures approach each other. In the intermediate coupling regime, the energy dissipation is most effective and both the doublon and hole temperatures decay quickly to a value close to the initial equilibrium temperature. In the polaronic regime, the doublons cool down very rapidly (although their effective temperature is no longer very well defined due to deviations of from a Fermi distribution function, see Fig. 7), while the transfer of kinetic energy from the holes to the doublons becomes increasingly inefficient, so that the holes remain hot. The inefficiency of the cooling via doublon-hole scattering is explained by the narrow width of the polaron sidebands, and the relatively large energy separation between them.

FIG. 7.

Spectral functions and distribution functions in the polaronic regime (U = 15, , g = 2, and 1.5). In the right hand panels, we plot the distribution functions in the energy range corresponding to a sizable hole and doublon population. For the later time (t = 40), we also show the Fermi function fits from which the effective doublon and hole temperatures have been extracted. For g = 2, the holes exhibit an inverted population. The high frequency noise is due to the cutoff of the Fourier integral. The inset in the top right panel shows the time evolution of the effective doublon and hole temperatures for g = 2 (solid lines) and g = 2.25 (dashed lines).

The g-dependence of the inverse effective temperatures is shown in the right panel of Fig. 6. The most efficient cooling occurs around (corresponding to ). For smaller g, the doublon and hole temperatures increase, while for larger couplings, only the hole temperature shows a marked increase. Interestingly, at , the holes reach an infinite temperature (flat) distribution, and for even stronger couplings, a negative temperature distribution (population inversion). From the comparison of the dashed and solid lines, which show the results for t = 20 and t = 40, it furthermore follows that the holes “cool down” in this negative temperature state ( decreases), i.e., the inverted population becomes steeper with increasing time (see also the inset of Fig. 7).

Figure 7 shows the spectral functions and distribution functions for g = 1.5 and g = 2, corresponding to the polaronic regime with positive (g = 1.5) and negative (g = 2) effective hole temperature. By comparing the spectral functions for the two coupling strengths, one notices that in the g = 2 case, the doublons populate in-gap states which fill the entire gap region. Due to these in-gap states, the doublon-hole recombination rate is substantially enhanced compared to g = 1.5. A possible explanation for the inverted hole population is that the energy in a recombination processes is partially transferred to the hole subsystem, where it is stuck because holes do no longer efficiently thermalize with the doublons in the polaronic regime. The recombination process is still entropically favorable (although the “entropy” of the holes decreases with increasing negative temperature) if a sufficient amount of energy is transferred to the bosons via the doublons. The rather long-lived transient state with inverted hole population thus results from doublon-hole recombinations, and the complicated energy exchange between holes, doublons and bosons in the polaronic regime with large ω0. On longer timescales, we expect that the system will thermalize in a state with identical, positive, temperatures for the doublons and holes.

SUMMARY AND OUTLOOK

We have discussed a DMFT based method for the simulation of nonequilibrium phenomena in the dynamic Hubbard model. This class of models has originally been introduced to capture the fact that atomic orbitals expand if occupied by two-electrons, which results in particle-hole asymmetric correlations. Dynamic Hubbard models have also recently appeared in the discussion of externally driven solids, where the orbital shapes are modulated periodically. We have considered here the model originally discussed by Hirsch, where the particle-hole asymmetric correlations are generated by a term gdX, which couples the doublon number operator to the position operator of an auxiliary boson.

One can think of numerous nonequilibrium setups in which the particle-hole asymmetry leads to interesting effects, which are absent in a conventional Hubbard model or Holstein-Hubbard model description. While the external driving of the boson in a homogeneous Holstein-Hubbard model would result in a trivial time-evolution, the driving of the boson in the dynamic Hubbard model may produce nontrivial effects, ranging from particle-hole excitations to dynamically induced boson-magnon couplings. If a temperature gradient is applied across a slab, the unequal doublon and hole mobilities should result in a nonzero particle current.

In the present work, we have focused our attention on photo-doped paramagnetic Mott insulators, where the asymmetric boson coupling leads to the rapid cooling of doublons and a hot distribution of holes, since the holes can dissipate their energy to the bosons only indirectly, via doublon-hole scattering. We have shown that a short time after a photo-doping excitation, the doublon and hole distributions in the upper and lower Hubbard band exhibit an approximately thermal shape, so that one can extract an effective doublon and hole temperature. We have analyzed the time evolution of these effective temperatures as a function of the model parameters. As in the Holstein-Hubbard model, the cooling rate seems to be determined by , and essentially independent of U. The fastest decrease of the effective hole temperature is observed for intermediate coupling strength (), while in the polaronic regime, the energy transfer from the hot holes to the trapped doublons becomes slow. If the energy spacing ω0 between the polaronic states is large enough, the photo-doping of the Mott insulator can even result in a long-lived transient state with an inverted hole population.

A long-lived negative temperature state could lead to potentially interesting effects. In the weakly correlated metallic regime, it has been shown that a system with a dynamically generated population inversion behaves as if the Coulomb interaction had been switched from repulsive to attractive. While the present situation with strongly correlated photo-doped holes is quite different, one might nevertheless expect different phases to occur for T > 0 and T < 0 when becomes small. Whether or not the inverted hole population leads to symmetry breaking to unconventional ordered states is an interesting open issue.

In the future, apart from studying driven systems and inhomogeneous set-ups, it will be interesting to extend the nonequilibrium DMFT simulations to other types of electron-boson couplings, such as the gdX2 coupling of Ref. 7. The extension of the Lang-Firsov scheme to these variants of the dynamic Hubbard model is not straight-forward and may require additional approximations. An alternative route is the development of combined strong/weak-coupling impurity solvers along the lines of Ref. 25.