Nonequilibrium GW+EDMFT: Antiscreening and inverted populations from nonlocal correlations

TL;DR

This paper investigates how nonlocal correlations affect screening dynamics in photo-excited Mott insulators, revealing phenomena like population inversion and antiscreening through a nonequilibrium GW+EDMFT approach.

Contribution

It introduces a nonequilibrium GW+EDMFT method to study long-range interaction effects in photo-doped Mott insulators, highlighting the emergence of antiscreening and inverted populations.

Findings

Long-lived transient states with nonthermal properties.

Population inversion of doublons and holes in systems with strong nonlocal interactions.

Detection of antiscreening via time-resolved electron-energy loss spectra.

Abstract

We study the dynamics of screening in photo-doped Mott insulators with long-ranged interactions using a nonequilibrium implementation of the $GW$ plus extended dynamical mean field theory ($GW$+EDMFT) formalism. Our study demonstrates that the complex interplay of the injected carriers with bosonic degrees of freedom (charge fluctuations) can result in long-lived transient states with properties that are distinctly different from those of thermal equilibrium states. Systems with strong nonlocal interactions are found to exhibit a self-sustained population inversion of the doublons and holes. This population inversion leads to low-energy antiscreening which can be detected in time-resolved electron-energy loss spectra.

Figures

Equations

$$\begin{gathered} H(t)=-J\sum _{⟨ij⟩\sigma }(e^{i{\varphi }_{ij}(t)}{c}_{i\sigma }^{†}{c}_{j\sigma }+h.c.)-\mu \sum _i{n}_i \\ +\sum _{i}U(n_{i\uparrow }-\frac{1}{2})(n_{i\downarrow }-\frac{1}{2})+\sum _{⟨ij⟩}V(n_i-1)(n_j-1), \end{gathered}$$

$${\Pi }_q^{\mathrm{R}}(\omega )=\sum _{k,{\omega }_1,{\omega }_2}{A}_k({\omega }_1)A_{k-q}({\omega }_2)\frac{f({\omega }_1)-f({\omega }_2)}{\omega -({\omega }_1-{\omega }_2)},$$

Full text

The development of time-resolved spectroscopic techniques provided important insights into the properties of complex materials Giannetti et al. (2016); Dal Conte et al. (2012, 2015); Kemper et al. (2016); Matveev et al. (2016), where charge, spin, orbital and lattice degrees of freedom are intertwined. A particularly exciting prospect is the nonequilibrium manipulation of material properties on electronic time scales, and the exploration of transient states that cannot be realized under equilibrium conditions. Prominent examples of this development are the laser-induced switching to a hidden state Stojchevska et al. (2014) in 1T-TaS2, and an apparent increase of the superconducting $T_{c}$ in phonon-driven cuprates and fulleride superconductors Kaiser et al. (2014); Mitrano et al. (2016).

Essential for the understanding of such experiments and phenomena is the ability to simulate relevant model systems using techniques that capture correlation effects in highly nonthermal states. Of particular importance is a proper description of the time-dependent screening processes, which determine the interaction parameters in such model Hamiltonians. The photo-induced change of screening was considered, e.g., as the cause of the collapse of the band gap in VO2 Wegkamp et al. (2014), or for an enhancement of excitonic order in Ta2NiSe5 Mor et al. (2016). Moreover, screening originates from charge fluctuations, which, similar to other bosonic modes like phonons Allen (1987); Werner and Eckstein (2013); Golež et al. (2012); Dorfner et al. (2015) or spin fluctuations Golež et al. (2014); Eckstein and Werner (2016), profoundly affect the relaxation pathway of the electronic distribution. As we will show in this paper, the fermionic dynamics and the bosonic screening modes are strongly coupled, so that their mutual interplay can lead to long-lived transient states which are entirely different from those characterizing equilibrium phases. These non-thermal states, with partially inverted populations, thus provide an intriguing pathway to novel light-induced properties.

A promising formalism to address these questions in strongly correlated solids is the combination of the $GW$ method and extended dynamical mean field theory ($GW$+EDMFT) Biermann et al. (2003); Werner and Eckstein (2016). Hedin’s $GW$ method Hedin (1965); Aryasetiawan and Gunnarsson (1998) is a weak coupling approach in which the self-energy is approximated by the product of the Green’s function $G$ and the screened interaction $W$. It captures nonlocal physics resulting from charge fluctuations, like screening, plasmonic collective modes and charge density waves. It however fails to describe strong correlation effects, like the Mott metal-insulator transition, which in turn are well described by the non-perturbative dynamical mean field theory (DMFT) Georges et al. (1996) and extended DMFT (EDMFT) Sun and Kotliar (2002). $GW$+EDMFT is a fully diagrammatic approach, which allows a self-consistent calculation of the screened interaction and its effect on the electronic properties in systems with long-ranged Coulomb interactions, and, in combination with a $GW$-based ab-initio simulation, a parameter-free simulation of weakly and strongly correlated materials. The recent equilibrium application of $GW$+EDMFT to model systems Ayral et al. (2013); Huang et al. (2014); Ayral et al. (2017) and real materials Boehnke et al. (2016) demonstrated the importance of dynamical screening effects originating from nonlocal interactions, e.g., for the proper interpretation of spectral features such as Hubbard bands and plasmon satellites. Here, we develop the nonequilibrium extension of the $GW$+EDMFT formalism and use it to study the effect of nonlocal interactions on the transient states and the relaxation dynamics of photoexcited carriers in Mott insulators.

As a simple but generic system with inter-site interactions we consider the single-band $U$-$V$ Hubbard model on the two-dimensional square lattice,

[TABLE]

where $c_{i\sigma}$ is the annihilation operator of a fermion with spin $\sigma$ on lattice site $i$, $n_{i}=n_{i\uparrow}+n_{i\downarrow}$, $\mu$ is the chemical potential, $U$ the on-site interaction, and $V$ the interaction between electrons on neighboring sites 111The extension to longer ranged interactions is straightforward Huang et al. (2014) and does not lead to qualitatively new physics.. The hopping integral $Je^{i\phi_{ij}(t)}$ (restricted to nearest neighbors) has a time-dependent Peierls phase $\phi_{ij}(t)=\int^{t}_{0}d\bar{t}\,\vec{E}(\bar{t})(\vec{r}_{i}-\vec{r}_{i})$ originating from an in-plane electric field $\vec{E}(t)$. In the following we will use the hopping amplitude $J\equiv 1$ as the unit of energy, and rewrite the interaction as $\frac{1}{2}\sum_{ij}v_{ij}\tilde{n}_{i}\tilde{n}_{j}$, where $\tilde{n}=n-1$ is the density fluctuation operator, and $v_{ij}=U\delta_{ij}+V\delta_{\langle ij\rangle}$.

The dynamics of the system is described in terms of the momentum-dependent electron Green’s function $G_{k}(t,t^{\prime})$ $=-i\langle T_{\mathcal{C}}c_{k}(t)c_{k}^{\dagger}(t^{\prime})\rangle$, and the charge correlation function $\chi_{q}(t,t^{\prime})=-i\langle T_{\mathcal{C}}\tilde{n}_{q}(t)\tilde{n}_{-q}(t^{\prime})\rangle$, which determines the (inverse) dielectric function $\varepsilon^{-1}_{q}=1+v_{q}\ast\chi_{q}$, and the screened interaction $W_{q}=\varepsilon^{-1}_{q}\ast v_{q}$, where $v_{q}$ is the Fourier transform of $v_{ij}$. In nonequilibrium, all quantities depend on two time arguments, or equivalently on time and frequency, and the $\ast$-product denotes convolution in time 222For a precise definition of time ordered correlation functions on the Keldysh contour, see the Supplemental MaterialSupplementary material which includes Refs. , .

To solve the extended Hubbard model in Eq. (1) we resort to the $GW$+EDMFT approximation Biermann et al. (2003), which can be derived using the Almbladh functional Almbladh et al. (1999). Nonlocal self-energy contributions for electrons and bosonic charge fluctuations are treated within the lowest order expansion of the functional (the $GW$ formalism), while the local contributions are included to all orders, by solving an auxiliary Anderson-Holstein impurity model with a self-consistently determined bosonic and fermionic bath. As a Green’s function based formalism, $GW$+EDMFT is not restricted to equilibrium or quasi-static problems, but can handle highly excited states. The derivation of the nonequilibrium formalism within the Keldysh framework is analogous to the equilibrium version Ayral et al. (2013); Hügel et al. (2016), and is presented in the Supplemental Material Supplementary material which includes Refs. , .

While powerful and numerically exact methods Werner and Millis (2010) exist for the solution of the $GW$+EDMFT equations in equilibrium, the application to nonequilibrium problems requires additional approximations at the level of the impurity solver. Since our goal is to study photo-doped Mott insulators, we use a perturbative solver that combines a self-consistent hybridization expansion (at first (second) order known as the non-crossing (one-crossing) approximation NCA (OCA) Eckstein and Werner (2010); Grewe and Keiter (1981); Coleman (1984)) with a weak-coupling expansion in the retarded density-density interactions. For technical aspects of the implementation, see Ref. Golež et al., 2015. As a benchmark, we show in Fig. 1(a) a comparison of the Matsubara component of the Green’s functions $G^{\text{Mat}}(\tau)$ for $U=10.5$, $V=1.5$ and inverse temperature $\beta=20$. (In the following local (nonlocal) correlators are distingushed by the absence (presence) of a subscript momentum label.) The NCA is found to overestimate the insulating nature of the solution Pruschke and Grewe (1989); Haule et al. (2001), as seen from $G^{\textrm{Mat}}(\beta/2)$, which can be taken as a measure for the spectral weight at the Fermi level. While this is a known artefact of the NCA Pruschke and Grewe (1989); Haule et al. (2001), the OCA substantially improves the accuracy of the solution compared to numerically exact Monte Carlo results Werner and Millis (2010). Furthermore, the finite-temperature metal-insulator transition is a crossover in the NCA description and becomes first order in the OCA solution, see Fig. 1(d). In the Mott phase, which we study here, NCA and OCA however yield qualitatively similar results, and we will resort to the numerically more tractable NCA in the following.

In the spectral function, shown in Fig. 1(b), the additional nonlocal $GW$ self-energy contributions in $GW$+EDMFT strongly enhance the plasmonic sideband at $\omega\approx\tfrac{3}{2}U$ and result in a slight reduction of the gap size compared to EDMFT. The inclusion of the nonlocal $GW$ diagrams in the $GW$+EDMFT approximation leads to a more metallic solution, since nonlocal correlations (in particular the nonlocal Fock term Ayral et al. (2017)) enhance the effective bandwidth. Also the local (momentum averaged) screened interaction $W$ is modified by the inclusion of the nonlocal polarization [Fig. 1(c)]. A noticeable feature is the strong enhancement of the plasmonic peak at $\omega\approx 12$ in comparison to EDMFT. A drawback of our approximate solver is evident at energies above the plasmon peak, where $\textrm{Im}[W]$ exhibits positive spectral weight, which is unphysical in thermal equilibrium. This problem arises because the NCA and OCA self-energies and polarizations are approximate strong-coupling solutions, which miss some of the local $GW$ diagrams. Numerically we found that these artefacts are most pronounced deep in the Mott phase, while close to the MIT transition and in the correlated metal $\textrm{Im}[W(\omega)]$ exhibits the expected analytical properties. Since the unphysical spectral weight appears only at very high energies, we believe that it is not crucial for the following discussion, which focuses on the low-energy screening properties of photodoped systems.

We now turn to the effect of nonlocal fluctuations on the relaxation dynamics after an electric field excitation. By applying a short pulse $E(t)=E_{0}e^{-4.6(t-t_{0})^{2}/t_{0}^{2}}\sin(\omega(t-t_{0}))$ with frequency $\omega=U$ and appropriately tuned amplitude $E_{0}$ a certain density of holon-doublon pairs is created. The width of the pulse $t_{0}=2\pi n/\omega$ is chosen such that the envelope accommodates $n=2$ electric field cycles. Deep in the Mott phase, the recombination of the holons and doublons after photo-excitation is strongly suppressed Sensarma et al. (2010); Eckstein and Werner (2010); Lenarčič and Prelovšek (2013). The photoexcited doublons can however relax within the upper Hubbard band, which manifests itself in the evolution of the kinetic energy. If the gap is small compared to the width of the Hubbard bands, the thermalization process, which involves impact ionization Werner et al. (2014), leads to an increase in the number of doublons $n_{d}$, see the EDMFT results (dashed lines) in Fig. 2(a). As already discussed in Ref. Golež et al. (2015) the inclusion of the nonlocal interactions on the EDMFT level decreases the relaxation times, due to the coupling to bosonic excitations (collective charge fluctuations). This picture remains valid if we include nonlocal self-energy and polarization effects in $GW$+EDMFT, but only if the nearest-neighbor interaction $V$ is small ($V\lesssim 0.5$). For larger values of $V$ (but still smaller than the critical value for the charge order transition), the double occupancy starts to decrease, which indicates that doublon-holon recombination occurs in the system, see solid lines in Fig. 2(a). Furthermore, the kinetic energy increases during the relaxation process, illustrated in Fig. 2(b), which is also intriguingly different from the behavior reported in previous photodoping studies Werner et al. (2014); Werner and Eckstein (2013); Golež et al. (2015).

In order to gain further insight into this intermediate $V$ regime we calculate the time and frequency-resolved spectral function of the system. After the pulse excitation of the system the spectral function $A(t,\omega)=-\frac{1}{\pi}\text{Im}[G^{R}(t,\omega)]$ remains almost unchanged, while the occupied density of states $N(t,\omega)=\textrm{Im}[G^{<}(t,\omega)]/2\pi i$ shows an increase of roughly 1% in the occupancy of the upper Hubbard band 333For all two-time quantities $O(t,t^{\prime})$ we define the partial Fourier-transform as $O(t,\omega)=\int_{t}^{t+t_{\text{cut}}}dt^{\prime}e^{\mathrm{i}\omega(t^{\prime}-t)}O(t^{\prime},t)$, with $t_{\text{cut}}=7$, which allows us to compute the spectrum with a time-independent resolution up to relatively long times. We checked that this choice of cut-off does not qualitatively affect the dynamics of the spectral properties under consideration., see Fig. 3(a). In agreement with the evolution of the kinetic energy, we observe a shift of the excited doublons toward higher energies, in contrast to previous DMFT and EDMFT studies Werner and Eckstein (2013); Werner et al. (2014); Golež et al. (2015) that consistently showed a relaxation of doublons to the lower edge of the upper Hubbard band. This $GW$+EDMFT evolution eventually results in a population inversion, as illustrated by the distribution function $f(t,\omega)=-2\text{Im}[G^{<}(t,\omega)]/\text{Im}[G^{R}(t,\omega)]$ shown in Fig. 3(b). We note again that this behavior is observed only for sufficiently large nonlocal interaction $V$.

The efficient recombination of doublon-hole pairs and the population inversion within $GW$+EDMFT can be understood by considering the two-particle properties, namely the screened interaction $W$ and the charge susceptibility $\chi_{q}.$ The time evolution of the local component of the screened interaction $W^{R,<}(t,\omega)$ for $U=10.5,V=1.5$ is shown in Fig. 3(d). In agreement with the previous EDMFT results, low energy screening channels appear as a consequence of photo-doping Golež et al. (2015). The main difference in $GW$+EDMFT is that the imaginary part of $W^{R}(t,\omega)$ changes sign as the system evolves into the population-inverted state. Since EDMFT and $GW$+EDMFT differ in the inclusion of nonlocal fluctuations we can qualitatively understand these results by evaluating the nonlocal charge susceptibility through the particle-hole bubble contribution to the polarization. In the stationary case the latter can be written as $\chi_{q}^{\mathrm{R}}=\Pi_{q}^{\mathrm{R}}[1-v_{q}\Pi_{q}^{\mathrm{R}}]^{-1}$ where the polarization $\Pi_{q}^{R}$ is given by

[TABLE]

which, in the case of well defined quasi-particles and thermal distributions $f$, reduces to the Lindhard formula. By exciting doublon-hole pairs in a Mott insulator, we temporarily create an inverted population in some energy range. Changing the Fermi-Dirac distribution function $f(\omega)$ in Eq. (2) to a partially inverted distribution function $\tilde{f}(\omega)$, we can change the sign of the numerator in $\chi_{q}^{\mathrm{R}}$ within a certain energy range. To illustrate this idea we evaluate $\chi_{q}^{\mathrm{R}}(\omega)$ using the Hubbard I approximation, where the lattice self-energy is approximated by the atomic limit self-energy. The resulting inverse dielectric function is shown in Fig. 4(a) where $\text{Im}[\varepsilon^{-1}_{q}(t,\omega)]=v_{q}\text{Im}[\chi_{q}^{\mathrm{R}}(t,\omega)]$. This leads to maximum spectral weight at the $\Gamma$-point and $\omega\approx U,$ which corresponds to charge excitations across the Mott gap. The lowest (highest) energies $U\pm W$ for which the imaginary part of the susceptibility $\text{Im}[\chi_{q}^{\mathrm{R}}(\omega)]$ has non-zero weight are at the $X$-point [$q=(\pi,\pi)$] 444An further increased nonlocal interaction strength $V$ would eventually lead to the condensation of the bosonic modes at $(\pi,\pi)$ and the formation of charge order.. In the case of the inverted population, see Fig. 4(c), the numerator in Eq. (2) becomes negative at frequencies corresponding to the energy width of the inverted regions, which leads to a negative spectral weight $-v_{q}\text{Im}[\chi_{q}^{R}(\omega)]<0$. These considerations show that the inclusion of nonlocal dynamical screening via the polarization bubble in $GW$+EDMFT is crucial for the appearance of the anti-screening phenomenon.

In contrast to the fermionic case, negative spectral weight in a steady state bosonic spectral function is not unphysical. The simplest example is a free oscillator, whose frequency suddenly turns unstable ($\omega_{0}<0$). Although there is no stable thermal equilibrium for $\omega_{0}<0$, the transient state remains well-defined, and its negative spectral weight reflects the possibility to increase fluctuations by emitting energy to the environment. The change of the sign of $\text{Im}[\chi_{q}(\omega)]$ in the photo-doped Mott insulator thus indicates a negative attenuation of charge fluctuations, which enable the system to emit low energy bosons to gain energy in the single particle sector. This also explains the unusual increase of the kinetic energy within the upper Hubbard band and the population inversion. A similar change in the sign of the susceptibilities was previously observed in models which are driven by (time-periodic) external fields Tsuji et al. (2009); Tsuji (2010). The intriguing observation in the present case is that the inverted population of the electronic states and the negative charge susceptibility mutually support each other (because the softening of charge fluctuations is caused by the change of the fermionic distribution), so that the peculiar state is self-sustained and stable as long as doublon-hole recombination processes inject energy into the bosonic subsystem.

A related population inversion was recently discussed in a study of Hirsch’s dynamic Hubbard model Hirsch (2001), although at unusually strong electron-phonon couplings. In the present case, the relevant strength $\lambda$ of the electron-boson coupling can be estimated from the density of states $D(t,\omega)$ of the bosonic modes in the auxiliary Anderson-Holstein impurity model (i.e., the boson-mediated density-density interaction interaction) as $\lambda=\int d\omega\sqrt{|\text{Im}D(\omega)|\omega}$ Ayral et al. (2013); Golež et al. (2015). As shown in Fig. 3(e), in $GW$+EDMFT, $\text{Im}|D(t,\omega)|$ features a pronounced peak at the energy of the gap size $\omega\approx 6,$ which corresponds to a very strong electron-boson coupling ($\lambda\approx 1.9$ for the largest value of $E$ plotted in the figure, if the integration range is chosen as $0\leq\omega\leq 8$).

Experimental probes which could be used to detect the peculiar charge fluctuation region are electron energy-loss spectroscopy (EELS) Kogar et al. (2014); Ibach and Mills (2013) and optical conductivity measurements Giannetti et al. (2016). The optical conductivity measures the frequency dependent optical constant near the $\Gamma$-point Eckstein and Kollar (2008); Lenarčič et al. (2014), while the EELS signal $-\text{Im}[v_{q}\varepsilon_{q}^{-1}(t,\omega)]=-v_{q}^{2}\text{Im}[\chi_{q}^{\mathrm{R}}(t,\omega)]$ measures the difference between the dielectric loss and gain (in equilibrium and at low temperatures, there is only loss). The generalization of EELS to the non-equilibrium situation, along the lines of the derivation of the time-dependent photo-emission formula Eckstein and Kollar (2008); Freericks et al. (2009); Randi et al. (2017), is presented in the Supplementary Material Supplementary material which includes Refs. , . The closely related inverse dielectric constant $\varepsilon_{q}^{-1}(t,\omega)$ shows a similar structure in $GW$+EDMFT as in the Hubbard I approximation, see Fig. 4(a) & (d). In particular, there is a pronounced maximum at the $\Gamma$-point at $\omega\approx U$ and dispersive bands with a minimal energy around the $X$-point. Immediately after the excitation the weight in the sub-gap region is increased in agreement with previous EDMFT results Golež et al. (2015), see Fig. 4(e). The initial increase in the screening in the sub-gap region however gives way to a negative spectrum as the inverted doublon population is formed [Fig. 4(f)-(h)], and the bosonic degrees of freedom also evolve into an inverted state. In this situation the energy gain for the probe electron at a certain energy in the EELS experiment is larger than the loss.

In conclusion, the nonequilibrium $GW$+EDMFT simulation revealed a self-sustained and long-lived transient population inversion as a result of the nontrivial energy exchange between doublons, holons and charge fluctuations. The existence of such a state provides an intriguing path to stabilize different types of light-induced order, which will be the subject of future investigations. Apart from these insights into the nonequilibrium properties of systems with nonlocal Coulomb interactions, our work represents an important step in the development of ab-initio simulation approaches for correlated systems in nonequilibrium states. The $GW$+EDMFT method implemented here features a fully consistent treatment of correlation and screening effects, and can in principle be combined with material-specific input from ab-initio $GW$ calculations within a multi-tier approach analogous to the scheme recently demonstrated for equilibrium systems in Ref. Boehnke et al., 2016.

Acknowledgements The calculations have been performed on the Beo04 cluster at the University of Fribourg. DG, LB, HS and PW have been supported by ERC starting grant No. 278023 and the SNSF through NCCR MARVEL and Grant No. 200021-165539. ME acknowledges support by the Deutsche Forschungsgemeinschaft within the Sonderforschungsbereich 925 (project B4).