Multi-band non-equilibrium GW+EDMFT formalism for correlated insulators

TL;DR

This paper develops a nonequilibrium multi-band GW+EDMFT formalism to study charge-transfer insulators under photo-excitation, revealing significant spectral and optical property changes driven by non-local charge fluctuations.

Contribution

It introduces a detailed derivation of a nonequilibrium multi-band GW+EDMFT framework and applies it to analyze photo-induced dynamics in cuprate-like insulators.

Findings

Charge-transfer gap is strongly renormalized after photo-excitation.

Momentum-dependent broadening observed in the upper Hubbard band.

Optical conductivity indicates gap renormalization proportional to photo-doping.

Abstract

We study the dynamics of charge-transfer insulators after a photo-excitation using the three-band Emery model which is relevant for the description of cuprate superconductors. We provide a detailed derivation of the nonequilibrium extension of the multi-band GW+EDMFT formalism and the corresponding downfolding procedure. The Peierls construction of the electron-light coupling is generalized to the multi-band case resulting in a gauge invariant combination of the Peierls intra-band acceleration and dipolar intra-band transitions. We apply the formalism to the study of momentum-dependent (inverse) photo-emission spectra and optical conductivities. The time-resolved spectral function shows a strong renormalization of the charge-transfer gap and a substantial broadening of some of the bands. While the upper Hubbard band exhibits a momentum-dependent broadening, an almost rigid band shift is observed for the ligand bands. The inverse photo-emission spectrum reveals that the inclusion of the non-local and inter-band charge fluctuations lead to a very fast relaxation of holes into the lower Hubbard band. Consistent with the changes in the spectral function, the optical conductivity shows a renormalization of the charge-transfer gap, which is proportional to the photo-doping. The details of the photo-induced changes strongly depend on the dipolar matrix elements, which calls for an ab-initio determination of these parameters.

Figures

Equations

$$\begin{array}{rl}H=& H_{\text{kin}}+H_{\text{int}},\\ H_{\text{kin}}=& \sum _{ij\sigma }\sum _{AB}{t}_{ij}^{AB}{c}_{i,A\sigma }^{†}{c}_{j,B\sigma }-\mu \sum _{iA}{n}_i^A,\\ H_{\text{int}}=& \sum _{i,A}{V}_{ii}^{AA}{n}_{i\uparrow }^A{n}_{i\downarrow }^A+\sum _{i\ne j\text{leavevmode}\text{or}\text{leavevmode}A\ne B}{V}_{ij}^{AB}{n}_i^A{n}_j^B,\end{array}$$

$$\Psi [G,W]=$$

$$-{\Psi }_{\text{weak}}[G_{ij,ab},W_{ij,ab}],$$

$${\Psi }_{GW+\text{EDMFT}}[G,W]={\Psi }_{\text{EDMFT}}[G_{ii,ab},W_{ii,ab}]$$

$$+{\Psi }_{GW}[G_{ij,AB},W_{ij,AB}]-{\Psi }_{GW}[G_{ii,ab},W_{ii,ab}].$$

$$\Sigma =\frac{\delta {\Psi }_{GW+\text{EDMFT}}}{\delta G},\Pi =-\frac{\delta {\Psi }_{GW+\text{EDMFT}}}{\delta W}.$$

$$S_{\text{imp}}^e[c_{a,\sigma }^{\ast },c_{b,\sigma }]={\int }_{\mathcal{C}}dtd{t}^{\prime }$$

$$\displaystyle\Big{\{}\sum_{\sigma,ab}c_{a,\sigma}^{*}(t)\left[\left((\mathrm{i}\partial_{t}+\mu)\delta_{ab}-\delta\mu_{ab}\right)\delta_{\mathcal{C}}(t,t^{\prime})-\Delta_{ab,\sigma}(t,t^{\prime})\right]c_{b,\sigma}(t^{\prime})$$

$$\displaystyle-\frac{1}{2}\sum_{ab}n^{a}(t)\mathcal{U}_{ab}(t,t^{\prime})n^{b}(t^{\prime})\Big{\}},$$

$$[G_k^{-1}{]}_{AB}(t,t^{\prime })=[(\mathrm{i}{\partial }_t+\mu ){\delta }_{AB}-{ϵ}_{k,AB}]\delta (t,t^{\prime })-{\Sigma }_{k,AB}(t,t^{\prime }),$$

$$[W_k^{-1}{]}_{AB}(t,t^{\prime })=[V_k^{-1}{]}_{AB}\delta (t,t^{\prime })-{\Pi }_{k,AB}(t,t^{\prime }),$$

$$\begin{array}{rl}& [G_k^{-1}{]}_{AB}(t,t^{\prime })=[(\mathrm{i}{\partial }_t+\mu ){\delta }_{AB}-({ϵ}_k+{\Sigma }_k^{HF}{)}_{AB}]{\delta }_{\mathcal{C}}(t,t^{\prime })\\ & -{\delta }_{Aa}{\delta }_{Bb}[{\Sigma }_{ab}^{\text{EDMFT}}(t,t^{\prime })-{\Sigma }_{\text{loc},ab}^{GW}(t,t^{\prime })]-{\Sigma }_{k,AB}^{GW}(t,t^{\prime }).\end{array}$$

$$[G_k^{-1}{]}_{AB}=[(\mathrm{i}{\partial }_t+\mu ){\delta }_{AB}-s_{AB}-{\Sigma }_{ab}^{\text{EDMFT}}{\delta }_{Aa}{\delta }_{Bb}]-[{\eta }_k{]}_{AB}.$$

$$\begin{array}{r}{s}_{AB}=\left(\begin{array}{cc}\frac{1}{N_k}{\sum }_k({ϵ}_{k,ab}+{\Sigma }_{k,ab}^{HF})& 0\\ 0& \frac{1}{N_k}{\sum }_k({ϵ}_{k,\bar{a}\bar{b}}+{\Sigma }_{k,\bar{a}\bar{b}}^{HF})\end{array}\right).\end{array}$$

$$[{\eta }_k{]}_{AB}=$$

$$+({\Sigma }_{k,AB}^{GW}-{\delta }_{Bb}{\delta }_{Aa}{\Sigma }_{\text{loc},ab}^{GW}).$$

$${\left(\begin{array}{cc}{M}_{11}& M_{12}\\ M_{21}& M_{22}\end{array}\right)}^{-1}=\left(\begin{array}{cc}[M_{11}-M_{12}{M}_{22}^{-1}{M}_{21}{]}^{-1}& -[M_{11}-M_{12}{M}_{22}^{-1}{M}_{21}{]}^{-1}{M}_{12}{M}_{22}^{-1}\\ -[M_{22}-M_{21}{M}_{11}^{-1}{M}_{12}{]}^{-1}{M}_{21}{M}_{11}^{-1}& [M_{22}-M_{21}{M}_{11}^{-1}{M}_{12}{]}^{-1}\end{array}\right).$$

$$\begin{array}{rl}& [g_k^{-1}{]}_{ab}(t,t^{\prime })=\\ & ([\mathrm{i}{\partial }_t+\mu ]{\delta }_{ab}-s_{ab})\delta (t,t^{\prime })-{\Sigma }_{ab}^{\text{EDMFT}}(t,t^{\prime })-{\eta }_{k,ab}(t,t^{\prime }),\end{array}$$

$$\begin{array}{r}{\bar{g}}_{k,\bar{a}\bar{b}}^{-1}(t,t^{\prime })=([\mathrm{i}{\partial }_t+\mu ]{\delta }_{\bar{a}\bar{b}}-s_{\bar{a}\bar{b}})\delta (t,t^{\prime })-{\eta }_{k,\bar{a}\bar{b}}(t,t^{\prime }).\end{array}$$

$$\begin{array}{rl}& G_{k,ab}^{-1}(t,t^{\prime })\\ & =([\mathrm{i}{\partial }_t+\mu ]{\delta }_{ab}-s_{ab})\delta (t,t^{\prime })-{\Sigma }_{ab}^{\text{EDMFT}}(t,t^{\prime })-{\tilde{\eta }}_{k,ab}(t,t^{\prime }),\end{array}$$

$$\begin{array}{rl}{\tilde{\Sigma }}_{k,ab}& ={\eta }_{k,a\bar{b}}\ast g_{k,\bar{b}\bar{c}}\ast {\eta }_{k,\bar{c}b}.\end{array}$$

$$\frac{1}{N_k}\sum _k[G_k{]}_{ab}=G_{\text{imp},ab}.$$

$$\begin{array}{rl}& G_{\text{imp},ab}^{-1}(t,t^{\prime })\\ =& \left([\mathrm{i}{\partial }_t+\mu ]{\delta }_{ab}-{\Sigma }_{\text{imp},ab}^H(t)-\delta {\mu }_{ab}-{ϵ}_{\text{loc},ab}\right){\delta }_{\mathcal{C}}(t,t^{\prime })\\ & -{\Delta }_{ab}(t,t^{\prime })-{\Sigma }_{ab}^{\text{EDMFT}}(t,t^{\prime }),\end{array}$$

$$\delta {\mu }_{ab}(t)$$

$$=s_{ab}(t)-{ϵ}_{\text{loc},ab}(t)-{\Sigma }_{\text{imp},ab}^H(t).$$

$$[G_k^{-1}{]}_{AB}$$

$$=[Z^{-1}-{\eta }_k{]}_{AB}.$$

$$\begin{array}{rl}{Z}_{AB}^{-1}=& {\delta }_{Aa}{\delta }_{Bb}[(\mathrm{i}{\partial }_t+\mu ){\delta }_{ab}-s_{ab}-{\Sigma }_{ab}^{\text{EDMFT}}]\\ & +{\delta }_{A\bar{a}}{\delta }_{B\bar{b}}[(\mathrm{i}{\partial }_t+\mu ){\delta }_{\bar{a}\bar{b}}-s_{\bar{a}\bar{b}}].\end{array}$$

$$G_k=Z+Z\ast {\eta }_k\ast G_k,$$

$$\begin{array}{rl}{G}_k^{-1}=& Z^{-1}-{\tilde{\eta }}_k,\\ G_{\text{imp}}^{-1}=& Z^{-1}-\Delta +[s-\delta \mu -{ϵ}_{\text{loc}}-{\Sigma }_{\text{imp}}^H].\end{array}$$

$$\begin{array}{rl}{G}_k=& Z+Z\ast {\tilde{\eta }}_k\ast G_k,\\ G_1=& \frac{1}{N_k}\sum _k{\tilde{\eta }}_k\ast G_k=\Delta \ast G_{\text{loc}},\\ G_2=& \frac{1}{N_k}\sum _k[{\tilde{\eta }}_k+{\tilde{\eta }}_k\ast G_k\ast {\tilde{\eta }}_k]=\Delta +\Delta \ast G_{\text{loc}}\ast \Delta ,\end{array}$$

$$[1+G_1]\ast \Delta =G_2.$$

$$[W_k^{-1}{]}_{AB}=[V_k^{-1}{]}_{AB}$$

$$-[{\Pi }^{\text{EDMFT}}{]}_{ab}{\delta }_{Aa}{\delta }_{Bb}-({\Pi }_{k,AB}^{GW}-{\Pi }_{\text{loc},ab}^{GW}{\delta }_{Aa}{\delta }_{Bb}),$$

Full text

I Introduction

Mott insulators are materials which according to band theory should be metals, but where the motion of the electrons is blocked by the large Coulomb repulsion in relatively localized $d$ or $f$ orbitals. Photo-excitation of charge carriers across the Mott gap, or the dielectric breakdown induced by a strong quasi-static field, results in a nonthermal conducting state. Both types of nonequilibrium phase transitions have been extensively studied experimentally Iwai et al. (2003); Perfetti et al. (2006); Okamoto et al. (2007); Mayer et al. (2015); Yamakawa et al. (2017) and using analytical Oka and Aoki (2010); Lenarčič and Prelovšek (2012) or computational tools Eckstein and Werner (2011, 2013). While the theoretical works have focused so far mainly on the simple case of a single-band Hubbard model, the experiments are typically performed on materials which would be classified as charge-transfer insulators Okamoto et al. (2007, 2011); Cilento et al. (2018); Peli et al. (2017). In these systems, ligand bands are located in between the upper and lower Hubbard bands, or in the energy range of the lower Hubbard band, and the applied fields may induce holes in the ligand bands, rather than in the lower Hubbard band Zaanen et al. (1985). For a better understanding of photo-induced nonthermal metallic states in Mott systems and a closer connection between theory and experiment, it is thus essential to extend the previous theoretical studies to charge transfer systems.

Here, we focus on the dynamics in the $d$-$p$ model for cuprates after field-pulse excitations, which is a set-up relevant for numerous experimental studies Okamoto et al. (2007, 2011); Cilento et al. (2018). Extending the results reported in a recent short study Golež et al. (2019), we will address several relevant issues which have not been considered in previous single-band model calculations: (i) what is the effect of the pulse energy on the photodoped metal state? Does it make a difference if holes are created in the $p$ bands or in the lower Hubbard band? (ii) How do photoinduced changes depend on the excitation protocols, which may include intraband (Peierls) acceleration and interband dipolar transitions? (iii) To what extent are the band gaps and band dispersions renormalized in the photodoped state? (iv) How are the photoinduced changes in the electronic structure reflected in the optical conductivity?

In terms of the formalism and numerical effort, the study of the $d$-$p$ model is considerably more involved than a simulation of a single-band Hubbard model with $s$ or $d$ orbitals. In the $d$-$p$ model, we need to treat different bands with different correlation strengths, as well as the dipolar terms appearing in the light-matter coupling. In this work, we will consider a combination of the extended dynamical mean-field theory (EDMFT) Sun and Kotliar (2002) and the $GW$ methodHedin (1965) ($GW$+EDMFT).Biermann et al. (2003); Boehnke et al. (2016); Ayral et al. (2013); Huang et al. (2014) EDMFT is suitable for describing Mott insulators and allows us to capture strong correlation effects in the $d$ orbitals, while the $GW$ method is appropriate for the description of the more weakly correlated $p$ orbitals, as well as the nonlocal charge fluctuations in the system. In order to assess the relevance of dynamical fluctuations for the description of equilibrium and nonequilibrium states, we furthermore compare the $GW$+EDMFT results to simpler Hartree-Fock+EDMFT (HF+EDMFT) calculations, in which the $p$-orbital interactions are treated at the Hartree-Fock level Weber et al. (2008).

Apart from its immediate relevance for the theoretical study of photoexcited cuprates, the implementation of a nonequilibrium $GW$+EDMFT scheme for the $d$-$p$ model represents an important step towards the development of ab-initio simulation approaches for nonequilibrium strongly correlated systems. The formalism captures time-dependent changes in the polarization function and allows to self-consistently compute dynamically screened, effective interaction parameters.Golež et al. (2017) This is particularly relevant for the description of excited states of solids, where for example the injection of charge carriers into a Mott insulator completely changes the screening properties.Golež et al. (2015) Since $GW$+EDMFT is rooted in the diagrammatic $GW$ formalism, it can furthermore be coupled to $GW$ ab-initio simulation codes in a way analogous to what has been demonstrated for equilibrium systems in Refs. Boehnke et al., 2016; Nilsson et al., 2017.

Selected results from our $GW$+DMFT investigation of the $d$-$p$ model have been reported in Ref. Golež et al., 2019. The purpose of the present work is to provide a detailed description of the formalism, and to present and analyze polarization- and momentum-dependent data. The paper is organized as follows. In Sec. II, we derive the general formalism which combines an EDMFT and $GW$ description of different orbital subspaces, as well as the corresponding downfolding procedure. We then define the $d$-$p$ model and apply the formalism to this specific model. The coupling of electrons to light for the present multiband problem is derived in Sec. III. The equilibrium results are presented in Sec. IV, and the dynamics after photoexcitation in Sec. V. We finish with the conclusions in Sec. VI.

II Model and method

II.1 Multitier $GW$+EDMFT

The multitier approachNilsson et al. (2017) is motivated by a general property of solid-state systems, where a small subset of orbitals close to the Fermi level are strongly correlated, while higher lying bands are weakly correlated. For example, in the cuprate-inspired models studied in this work, the $p$ orbitals are weakly correlated and can be treated by perturbation theory, while the low-lying $d$ orbitals with strong correlations require a more sophisticated treatment, such as EDMFT. Therefore we separate the full orbital space into strongly and weakly correlated subspaces. To keep the procedure general, we will denote the degrees of freedom in the full space by capital letters, e.g., $A$, those in the correlated subspace by lower case letters, e.g., $a$, and those in higher lying orbitals by lower case letters with an overline, e.g., $\bar{a}$. The presented formalism is therefore applicable to the general multiorbital case in both subspaces. Only in Sec. II.3, will we restrict the discussion to the specific example of cuprates.

The formalism applies to a general Hamiltonian of the form

[TABLE]

where $i$ and $j$ denote the spatial index of the unit cell and $\mu$ is the chemical potential. We have restricted ourselves to density-density interactions and explicitly separated the local interband interaction term. Here, we have introduced the spin dependent density operators $n_{i\sigma}^{A}=c_{i,A\sigma}^{\dagger}c_{i,A\sigma}$ and the total density operators $n_{i}^{A}=\sum_{\sigma}n_{i\sigma}^{A}$. The on-site energies are fixed by the local-diagonal part of the single-particle term $t_{ii}^{AA}$, and the remaining matrix elements $t_{ij}^{AB}$ define the hopping amplitudes.

The many-body complexity of the system is captured by the Almbladh functional $\Psi[G,W]$ which contains all possible single-particle irreducible diagrams built with the single-particle propagator $G$ and the screened interaction $W$Almbladh et al. (1999).

In general, we aim to construct a conserving approximation for the total system starting from two approximate functionals, one for the correlated subsystem $\Psi_{\text{strong}}[G,W]$ and one for the higher lying orbitals $\Psi_{\text{weak}}[G,W]$ as

[TABLE]

where the last term removes the double counting. In the $GW$+EDMFT approximation, strong correlations are treated within a local EDMFT approximation $\Psi_{\text{EDMFT}}[G_{ii,ab},W_{ii,ab}]$. The higher-lying orbitals and nonlocal fluctuations are treated within the $GW$ approximation, i.e., the lowest order diagram in the weak-coupling expansion of \Psi_{\text{GW}}[G,W]. The corresponding Almbladh functional is given byBiermann et al. (2003); Ayral et al. (2013)

[TABLE]

From the stationarity condition we obtain the expression for the self-energy $\Sigma$ and the polarization $\Pi$ as

[TABLE]

Within EDMFT, the local propagators $G_{ii,ab}$ and $W_{ii,ab}$ in the correlated subspace are obtained from the solution of an effective quantum impurity model with actionSun and Kotliar (2002); Ayral et al. (2013); Golež et al. (2015)

[TABLE]

where the retarded interaction $\mathcal{U}$ is coupled to the density operator $n^{a}$. Here we have separated the self-consistent hybridization function into a time-local part $\delta\mu$, and a contribution $\Delta$ which has no instantaneous term. The origin of the additional shift $\delta\mu$ is that the lattice Hartree-Fock term includes nonlocal contributions, which are not present at the impurity level. The hybridisation function $\Delta(t,t^{\prime})$ is determined from the self-consistency condition $G_{ii,ab}=G_{\text{imp},ab}$, and the effective impurity interaction $\mathcal{U}(t,t^{\prime})$ is implicitly defined by the condition $W_{ii,ab}=W_{\text{imp},ab}.$ From now on, we will assume a translationally invariant system, although the formalism is easily extended to the spatially inhomogeneous case, and introduce momentum $k$ as a good quantum number. $G$ and $W$ in the full orbital space are given by the fermionic and bosonic Dyson equations

[TABLE]

[TABLE]

where the dispersion $\epsilon_{k,AB}$ is the Fourier transform of the hopping amplitudes $t_{ij}^{AB}.$ These relations conclude the formal derivation of the theory.

The restriction to density-density interactions considered in this work can be overcome. Here, one should distinguish two cases: i) Additional local interaction terms in the correlated subspace, such as a Kanamori interaction, can be trivially incorporated in a multiorbital extension of the present formalism, as presented in the next section. ii) A treatment beyond density-density interactions outside the local correlated subspace would require a proper EDMFT decoupling and an introduction of new auxiliary fields and two-particle selfconsistency loops. While this is in principle feasible and an important task for the realistic description of correlated materials, it goes beyond the scope of this paper.

II.2 Downfolding

While the perturbative $GW$ self-energies and polarizations are explicit functions of $G$ and $W$, $\Sigma^{\text{EDMFT}}$ and $\Pi^{\text{EDMFT}}$ are only implicitly defined through the impurity model. In order to determine the auxiliary quantities $\Delta$ and $\mathcal{U}$, we downfold the electronic and bosonic Dyson equations, see Eqs. (6) and (7), to the correlated subspace, as described in the following paragraph.

Fermionic self-consistency loop.

We first discuss the downfolding for the electronic Green’s function. We can rewrite the fermionic Dyson equation (6) by separating the local and nonlocal self-energies

[TABLE]

Here, $\Sigma^{HF}$ is the lattice Hartree-Fock contribution evaluated in the full space. The local self-energy of the correlated subspace is described within EDMFT by $\Sigma^{\text{EDMFT}}_{ab}$, and the nonlocal correlations in the entire space are described at the $GW$ level by $\Sigma^{GW}_{k,AB}$.Hedin (1965) Since instantaneous contributions can be treated at the level of the effective Hamiltonian, we explicitly separate the Hartree-Fock contribution $\Sigma_{k}^{HF}$ from the higher order EDMFT ($GW$) self-energy contributions $\Sigma^{\text{EDMFT}}$ $(\Sigma^{GW}).$ The term $\Sigma_{\text{loc},ab}^{GW}=\frac{1}{N_{k}}\sum_{k}\Sigma^{GW}_{k,ab}$ is the local component of the $GW$ self-energy in the correlated subspace, which is subtracted to avoid a double counting of local self-energy contributions in this subspace. Here and in the following we use a short notation for the sum over repeated indices, $\delta_{Aa}X_{ab}\delta_{bB}\equiv\sum_{ab}\delta_{Aa}X_{ab}\delta_{bB}$. In the following, we will separate all diagonal and local contributions to the self-energy from the rest. This step is a convenient choice for calculations based on strong coupling impurity solvers, where we want to avoid the explicit use of the self-energy in the self-consistency loop, see below for a detailed discussion. We rewrite the lattice Dyson equation as

[TABLE]

Here, we have separated the local part of the mean-field dispersion, which is now a block-diagonal matrix with strongly and weakly correlated subspaces

[TABLE]

This separation shifts all the block off-diagonal and block diagonal nonlocal contributions into

[TABLE]

We note that the local components vanish in the correlated subspace by construction, $\frac{1}{N_{k}}\sum_{k}[\eta_{k}]_{ab}=0.$

Once we have solved the lattice Dyson equation, we have to define the auxiliary impurity model (hybridization function $\Delta$ and retarded interaction $\mathcal{U}$) for the correlated subspace. The main identity to downfold the Dyson equation from the full space to the correlated subspace is the equation for the inverse of a block matrix

[TABLE]

In order to apply it to the lattice Dyson equation (9) we define an auxiliary propagator for the correlated subspace $g_{k}$,

[TABLE]

and similarly for the high-lying orbitals

[TABLE]

Using Eqs. (11) and (12), the lattice Dyson equation for the correlated subspace may now be rewritten as

[TABLE]

where $\tilde{\eta}_{k,ab}(t,t^{\prime})=\eta_{k,ab}(t,t^{\prime})+\tilde{\Sigma}_{k,ab}(t,t^{\prime}),$ and the auxiliary self-energy is

[TABLE]

Here, $*$ denotes a matrix multiplication in the orbital sector and a convolution of times. (For HF+DMFT, $\eta_{k}$ is instantaneous and the convolutions reduce to multiplications.)

We determine $\Sigma^{\text{EDMFT}}(t,t^{\prime})$ by solving the auxiliary impurity problem and imposing the self-consistency in the correlated subspace,

[TABLE]

The auxiliary impurity problem fulfills the Dyson equation, which is understood as a matrix equation in orbital space

[TABLE]

where the impurity Hartree term is given by $\Sigma^{H}_{\text{imp},ab}(t)=\delta_{ab}\int d\bar{t}\leavevmode\nobreak\ \mathcal{U}_{ac}(t,\bar{t})\langle n^{c}(\bar{t})\rangle$ and the local energy level by $\epsilon_{\text{loc}}=\frac{1}{N_{k}}\sum_{k}\epsilon_{k}.$ The additional chemical potential shift $\delta\mu$ corresponds to the difference between the Hartree-Fock contribution on the lattice and on the impurity level,

[TABLE]

Self-consistency without self-energy.

Most impurity solvers provide a way to compute the Green’s function rather than the self-energy. In particular, this is the case for the strong-coupling expansion (NCA) which will be used in this work Grewe and Keiter (1981); Eckstein and Werner (2010) but also DMRG Wolf et al. (2014); Balzer et al. (2015) and QMC Hirsch and Fye (1986); Rubtsov et al. (2005); Werner and Millis (2006). In this case, the determination of the time-dependent self-energy from an inversion of the impurity Dyson equation can be a numerically ill-conditioned problem Aoki et al. (2014), and it is preferable to close the self-consistency cycle without explicitly evaluating the self-energy. The main idea is to rewrite the lattice Dyson equation with the auxiliary function $Z$, which is purely local and diagonal and shift the rest into $\eta_{k}$. Here, we rely on the simple relation that the inverse of a block-diagonal matrix is a block-diagonal matrix, which decouples the evaluation of the auxiliary quantity $Z$ for each subsector. We should stress that any other choice would lead to a coupling between different orbitals and would require an evaluation of the impurity self-energy. To be explicit, we rewrite the lattice Dyson equation as

[TABLE]

The inverse of the local propagator $Z$ explicitly reads

[TABLE]

The elements $Z_{ab}$ are obtained from the impurity Dyson equation $G^{\text{imp}}=Z+Z*\Delta*G^{\text{imp}}$ as a matrix equation in the correlated subspace, while the rest are free propagators. The solution of the lattice Dyson equation is then obtained from

[TABLE]

which is a matrix equation in the full orbital space. This is the stable form of the Dyson equation, which is used in our numerical implementation. $\eta_{k}$ is obtained from Eq. (11) and in the case of $GW$+EDMFT has an instantaneous and a retarded part. Then, by imposing the self-consistency condition $G_{\text{loc},ab}=G_{\text{imp},ab}$ we obtain the equation for the hybridization function $\Delta$ from Eqs. (15) and (18). In terms of $Z$, these equations read

[TABLE]

The term in the square bracket is eliminated by the choice of $\delta\mu$ in Eq. (19). Then the hybridization function $\Delta$ is determined from Eqs. (23), in a way analogous to the single-band case:Eckstein and Werner (2011); Aoki et al. (2014)

[TABLE]

so that $\Delta$ can be calculated by solving

[TABLE]

All relations from Eq. 22 to Eq. 25 should be understood as block-diagonal matrix equations with blocks corresponding to the strongly correlated and weakly correlated subspaces. The last expression shows that the hybridization function $\Delta$ includes information about direct hopping events to the neighboring sites via $\eta_{k}$, as well as the hybridization between the two subspaces via $\tilde{\eta}_{k}.$ Note that $G_{2}$ does not have an instantaneous part, because $\sum_{k}\tilde{\eta}_{k}=\sum_{k}\tilde{\Sigma}_{k}$ has no instantaneous part and in the correlated subspace $\sum_{k}\eta_{k}=0$ by construction. Our numerical implementation is based on the stable version of the Dyson equation [Eq. (25)] and $\tilde{\eta}_{k}$ is evaluated using the auxiliary self-energy $\tilde{\Sigma}_{k}$ defined in Eq. (16).

Bosonic self-consistency loop

We next discuss the downfolding for the screened interaction. The bosonic lattice Dyson equation is

[TABLE]

while the impurity Dyson equation reads

[TABLE]

We will use these equations and the possibility to calculate the charge susceptibility $\chi_{\text{imp}}$ of the impurity model. The charge susceptibility acts like a $T$-matrix for the bosonic impurity Green’s function and allows us to calculate the impurity polarization $\Pi_{\text{imp}}$ as follows: $W_{\text{imp}}=\mathcal{U}+W_{\text{imp}}*\Pi_{\text{imp}}*\mathcal{U}=\mathcal{U}-\mathcal{U}*\chi_{\text{imp}}*\mathcal{U}$, or

[TABLE]

In the lattice Dyson equation, we set the elements of the lattice polarization $\Pi^{\text{EDMFT}}$ to $\Pi_{\text{imp}}$ within the correlated subspace. The rest of the lattice polarization is filled with the $GW$ polarization, where local components in the correlated subspace have been subtracted due to the double counting. It thus becomes

[TABLE]

and from $W_{k}$ we can calculate the local screened lattice interaction $W_{\text{loc}}.$ This is once again a stable form of the bosonic Dyson equation, which is employed in our implementation. The new approximation for the impurity effective interaction $\mathcal{U}$ is obtained by identifying the local screened lattice interaction in the correlated subspace with the screened impurity interaction, $W_{\text{loc}}=W_{\text{imp}}$:

[TABLE]

which is understood as a matrix equation in the correlated subspace.

Summary of the self-consistency loop

The full self-consistency loop involves the following steps:

Start with an initial guess for $\Delta$ and $\mathcal{U}$ and $s_{ab}=\epsilon_{\text{loc},ab}+\Sigma^{HF}_{\text{loc},ab}$. 2. 2.

Solve for $G_{\text{imp}}$ and $\chi_{\text{imp}}$ in the subspace with the impurity Hartree term $\Sigma^{H}_{\text{imp}}$ and $\delta\mu=\Sigma^{HF}_{\text{loc}}-\Sigma^{H}_{\text{imp}}$. 3. 3.

Invert $Z+Z*\Delta*G_{\text{imp}}=G_{\text{imp}}$ to get $Z$ in the correlated subspace, while the $Z_{\bar{a}\bar{b}}$ are obtained from the free solution $Z_{\bar{a}\bar{b}}=(\mathrm{i}\partial_{t}+\mu-s_{\bar{a}\bar{b}})^{-1}$. The total $Z$ is constructed as in Eq. (21). 4. 4.

Get the local polarization in the correlated subspace by solving Eq. (28). 5. 5.

Compute the lattice propagator $G_{k}$ ($W_{k}$) in the full orbital space for all momenta $k$ using the fermionic (bosonic) lattice Dyson equation, see Eq. (22) (Eq. (29)). These are both stable versions of the Dyson equation. 6. 6.

Optional: Evaluate the nonlocal $GW$ self-energy $\Sigma_{k,AB}^{GW}$ and polarization $\Pi_{k,AB}^{GW}.$ 7. 7.

Compute the $\bar{g}_{k}$ in the weakly correlated subspace by solving $\bar{g}_{k}=\bar{g}_{k,0}+\bar{g}_{k,0}*\eta_{k}*\bar{g}_{k}$ for $\bar{g}_{k}$, where $\bar{g}_{k,0}^{-1}=\mathrm{i}\partial_{t}+\mu-s_{\overline{a}\overline{b}}$. As the self-energy in the weakly correlated subspace is explicitly known, this is just an ordinary Dyson equation. Compute the auxiliary self-energy $\tilde{\Sigma}_{k}=\tilde{\eta}_{k}^{\dagger}\bar{g}_{k}\tilde{\eta}_{k}$ for the downfolding. 8. 8.

Solve Eqs. (24) and (25) to obtain the new approximation for $\Delta$ and calculate the updated $\delta\mu$. 9. 9.

Obtain the new approximation for the effective impurity interaction $\mathcal{U}$ by solving Eq. (30).

II.3 $d$-$p$ model

We will now apply the formalism introduced in the previous section to the two-dimensional Emery model for cupratesEmery (1987), which includes the correlated Cu $d_{x^{2}-y^{2}}$ and higher lying $O$ $p_{x}$ and $p_{y}$ orbitals. The position of the $p_{x}$ $(p_{y})$ orbital is shifted with respect to the $d$ orbital by half a unit vector in the $x$ $(y)$ direction. Since the local interactions in the $d$ orbital are stronger than in the $p$ orbitals, we restrict the correlated subspace (EDMFT treatment) to the $d$-orbital, while the $p$ orbitals are treated at the Hartree-Fock or $GW$ level.

The system Hamiltonian is parametrized by Eq. (1). The on-site energies are given by the local-diagonal component of the single-particle term $t_{ii}^{AA}=\{\epsilon_{d},\epsilon_{d}+\Delta_{pd},\epsilon_{d}+\Delta_{pd}\}.$ Here, the band separation $\Delta_{pd}$ determines the difference between the on-site energy for the $d$ orbital, $\epsilon_{d}$, and the on-site energy for the $p$ orbitals, $\epsilon_{d}+\Delta_{pd}.$ We denote the nearest neighbor hopping between the $d$ and $p_{x}$ and $p_{y}$ orbitals by $t_{ii}^{dp_{y}}=t_{ii}^{dp_{y}}\equiv t^{dp}$ and the hopping amplitude between the $p_{x}$ and $p_{y}$ orbitals by $t_{ii}^{p_{x}p_{y}}\equiv t^{pp}.$ The density-density interaction vertex is given by $V_{ii}^{dd}=U_{dd}$ and $V_{ii}^{dp_{x}}=V_{ii}^{dp_{y}}=\frac{1}{2}U_{dp}$ for nearest-neighbor $d$ and $p$ orbitals.

It is useful to introduce a spinor for the unit cell $i$,

[TABLE]

and its Fourier transform

[TABLE]

which allows to write the single particle part of the Hamiltonian in the momentum representation as

[TABLE]

Here, “h.c.” stands for the Hermitian conjugate part of the matrix. We note that in some previous studies a Fourier transform with a shifted position of the $p$ orbitals has been used.Maekawa et al. (2013); Hansmann et al. (2014)

By performing the same Fourier transform as for the single particle part of the Hamiltonian we arrive at the following momentum representation of the interaction

[TABLE]

where special care needs to be taken for the local and orbital diagonal components to fulfill the Pauli exclusion principle, see Sec. II.4. We solve this lattice problem using the $GW$+EDMFT method Biermann et al. (2003); Ayral et al. (2013); Huang et al. (2014); Golež et al. (2017) introduced above, which is based on EDMFT.Sun and Kotliar (2002); Ayral et al. (2013); Huang et al. (2014); Golež et al. (2015) The downfolding is thus performed from the three-orbital space of the full model to the $d$ orbital correlated subspace both for the electronic Green’s function and the (bosonic) screened interaction.

A nontrivial problem in practice is the convergence. Because the Hartree and Fock terms are quite large, the solution can jump from fully occupied to fully unoccupied $d$ orbitals during the self-consistency cycle. This is why we fix the total occupation to $n$ and the relative occupation between the orbitals $r=n_{p}/n_{d},$ which determines the local-diagonal part of the density matrix and the local Hartree and Fock shifts. $r$ is determined by minimizing the difference between the fixed and calculated local density matrix at a given chemical potential. Once a sufficiently good approximate solution is obtained, we release the constraint on the local block-diagonal matrix.

The electronic propagator is defined as

[TABLE]

with $\langle\ldots\rangle=\frac{1}{Z}\text{Tr}[\mathcal{T}_{\mathcal{C}}e^{\mathrm{i}S_{\text{latt}}}\ldots]$ and $S_{\text{latt}}$ the lattice action, $Z$ the partition function of the initial equilibrium state with inverse temperature $\beta$, and $\mathcal{T}_{\mathcal{C}}$ the time ordering operator on the Kadanoff-Baym contour which runs from time [math] to time $t$ along the real time axis, back to time zero and then to $-\mathrm{i}\beta$ along the imaginary-time axis.Aoki et al. (2014) The free propagator is similarly obtained from the single-particle Hamiltonian, or the corresponding lattice action. For the matrix elements of the tight-binding approximation, see Ref. Hansmann et al., 2014. Since there is hopping between all orbitals in the unit cell, the free propagator is a full matrix.

The bosonic propagator is defined as

[TABLE]

where the spinor $\phi_{k\sigma}=\{\phi_{k,d\sigma}$, $\phi_{k,p_{x}\sigma}$, $\phi_{k,p_{y}\sigma}\}$ originates from a Hubbard-Stratonovich decoupling of all the interaction terms in $H_{\text{int}}$.

II.4 Self-energy and Polarization diagrams

In a conserving approximation, the selfenergy and polarization are given by the derivative of the Luttinger-Ward functional with respect to the Green function and screened interaction, respectively. We will provide here the expressions for the mean-field (nonretarded Hartree+Fock) and the retarded $GW$ approximations.

Hartree

The Hartree expression for the Emery model is given by

[TABLE]

where we have introduced the single-particle density matrix $\rho_{k,\alpha\beta}=\langle c^{\dagger}_{\beta,k}c_{\alpha,k}\rangle$, its local value $\rho_{\alpha\beta}=\frac{1}{N_{k}}\sum_{k}\rho_{k,\alpha\beta}$, and the combined index $\alpha=(A,\sigma),$ where $A$ is the band and $\sigma$ the spin index. As marked in the sum we need to fulfill the Pauli exclusion principle and therefore we only sum over distinct indices $\beta\neq\alpha$. For the $d$ orbital the expression is given by

[TABLE]

where $V^{dd,\text{nl}}$ marks a possible nonlocal $d$-$d$ interaction. We have assumed spin symmetry and used the notation $\rho_{AB}=\sum_{\sigma}\rho_{AB,\sigma}$.

Now, we can calculate the time-dependent shift of the chemical potential for the impurity problem (see Eq. (19)), where the mean-field contributions are given by

[TABLE]

Here we have introduced the retarded part of the effective impurity interaction as $\mathcal{D}=\mathcal{U}-V^{dd}_{\text{loc}}$. We can get rid of the term $\textstyle{\int}d\overline{t}\mathcal{D}(t,\overline{t})\rho_{dd}(\overline{t})$ by performing an expansion in the density fluctuations $\bar{n}=n-\langle n\rangle$ instead of $n$, see discussion in Ref. Golež et al., 2015.

Fock

In order to fulfill the Pauli exclusion principle we should remove the local diagonal part of the Fock diagram,

[TABLE]

where $V^{\alpha\alpha}_{\text{loc}}$ is the local intraorbital interaction and repeated indices are understood as a piecewise product. To give a concrete example, for the diagonal $d$ component this becomes

[TABLE]

where we included the possibility of a nonlocal $d$-$d$ interaction. The off-diagonal parts of the self-energy are not modified, e. g.

[TABLE]

$GW$

The $GW$ self-energy is given by

[TABLE]

from which we have to subtract the local term already present in the impurity self energy. The explicit expression for the diagonal $d$ contribution is

[TABLE]

while the off-diagonal contributions are given by

[TABLE]

The polarization is calculated in an analogous fashion,

[TABLE]

where the sum is done only through spin indices $\sigma$ for combined indices $\alpha$ and $\beta$. The $d$-$d$ and $d$-$p$ components thus are

[TABLE]

and

[TABLE]

and analogous for all the other terms. From this it follows that the polarization matrix is a full matrix.

III Coupling with light

In order to simulate the photo excitation, we need to understand how light couples to matter in the multiband case. In principle, one can obtain a light-matter Hamiltonian by projecting the continuum theory onto a given set of bands. There is however a well-known gauge ambiguity:Loudon (2000) While the full theory can be described equivalently in different gauges for the electromagnetic fields, the description resulting from a projection to a given set of bands or Wannier orbitals can depend on the gauge in which the projection is performed. In particular, different descriptions are obtained by projecting the continuum Hamiltonian in the Coulomb gauge, which corresponds to the replacement $\vec{p}\rightarrow\vec{p}-e\vec{A}$, or in the dipolar gauge, where the coupling is between a polarization density and the electric field. This ambiguity is actively discussed for atomic Hamiltonians in the context of cavity quantum electrodynamics.De Bernardis et al. (2018); Di Stefano et al. (2019)

In the present case, we are interested in the semiclassical description (the fields are treated classically). Nevertheless, the question of the optimal light-matter Hamiltonian becomes relevant for strongly driven systems, beyond linear response. A common strategy is to request the multiband light-matter Hamiltonian to have a gauge structure even after projection to a given manifold of bands.Boykin et al. (2001) This leads to a Hamiltonian which in general includes both multipolar magnetic and electric matrix elements, and Peierls phases depending on the vector potential. Here we give a compact summary of the corresponding construction, which from the outset takes into account the main simplifications which are appropriate for the description of optically driven solids: In particular, all fields are assumed to be spatially homogeneous throughout the solid ($q=0$ approximation). This is suitable in the long wavelength limit; it effectively neglects magnetic dipolar matrix elements, which are relatively small even for electric fields of the order of MV/m.

Before presenting this construction, we remark that the $q=0$ approximation can be dropped relatively easily, and one can instead obtain the same light matter Hamiltonian by only neglecting the variation of the fields on the atomic scale, derive a Hamiltonian with spatially varying fields, and restrict the fields to $q=0$ in the end. Moreover, by making the semiclassical approximation from the outset one does not capture induced long-range dipolar interactions and Lamb-shifts. However, in the present case all interactions (including the long-range Coulomb interaction) are approximated by an effective short-range interaction, whose ab-initio determination is challenging even in the case without light-matter coupling. An alternative but more lengthy derivation of the Hamiltonian, which leads to the same Hamiltonian as below, without making the $q=0$ and semiclassical approximation from the outset, can be found in Ref. Jiajun et al., 2019.

For the construction of the light-matter Hamiltonian, we start with the semiclassical coupling (which is obtained from the minimal coupling Hamiltonian using a Power-Zienau-Wolley (PZW) transformation, see Ref. Loudon, 2000)

[TABLE]

where we have introduced the dipolar moment $\mathcal{\vec{P}}=e\sum_{j}\vec{r}_{j}$, $\vec{E}$ is the (classical) and spatially homogeneous transversal component of the electric field, $\vec{r}_{j}$ is the position of the $j$-th electron and $e$ the electric charge. We then make the following approximations:

Project the continuum Hamiltonian to a set of localized Wannier basis functions and explicitly separate intra- and inter-unit-cell contributions. 2. 2.

As the inter-unit-cell terms break the translational invariance, we apply a unitary transformation which restores the translational symmetry at the expense of time-dependent phase factors. This step is a generalization of the Peierls substitution. The transformation of the inter-unit cell terms will reintroduce the vector potential $\vec{A}$ into the description. The mixed representation with the electric field $\vec{E}$ and the vector potential $\vec{A}$ is the main result of this section. It allows the modeling of translationally invariant systems with intra-band acceleration and inter-orbital dipolar excitations.

In more detail, the projection of (50) to a set of Wannier orbitals is obtained as follows: The unit cell may contain several atoms described by Wannier functions, which we denote by $\langle\vec{r}|k,n\rangle=\phi_{n}(\vec{r}-\vec{R}_{k})$ for the $k$-th unit cell and $n$-th orbital in the unit cell. We thus obtain the matrix elements of the position operator as:

[TABLE]

where the vector $\vec{R}_{k}$ represents the center of mass in the unit cell $k$. We have introduced the dipolar matrix element as $\vec{D}_{ik}^{nm}=e\langle in|\vec{r}-\vec{R}_{k}|km\rangle.$

The tight-binding Hamiltonian on the lattice can now be written as

[TABLE]

The inter-unit-cell contributions to the Hamiltonian $e\vec{E}\sum_{\begin{subarray}{c}i\end{subarray}}n_{i}\vec{R}_{i}$ break the translational invariance. In order to restore the translational symmetry, we perform the unitary transformation $U=e^{\mathrm{i}e\sum_{j}\chi_{j}n_{j}}.$ This transformation removes the $e\vec{E}\sum_{j}n_{j}\vec{R}_{j}$ term and introduces time-dependent phase factors as in the Peierls substitution:

[TABLE]

from where we get the condition $\dot{\chi}=-\vec{E}(t)\cdot\vec{R}_{i}\rightarrow\vec{A}(t)=-\int^{t}\vec{E}(s)ds$ in order to cancel the spatially inhomogeneous electric field contributions. The Peierls phases are thus given by $\chi_{ij}=\chi_{i}-\chi_{j}=\int_{\vec{R}_{j}}^{\vec{R}_{i}}\vec{A}(x)d\vec{x}=(\vec{R}_{i}-\vec{R}_{j})\vec{A}$. At this point, the integral is independent of the path, as the vector potential is homogeneous and can be written as a gradient field. The final Hamiltonian is now given by

[TABLE]

The Peierls transformation yields time-dependent hoppings, but also the nonlocal dipolar matrix elements have a time-dependent phase. Up to this point, the derivation was rather general (up to neglecting quantum effects and making the $q=0$ approximation). This Hamiltonian can now be restricted to a given subset of bands. The important point is that the Hamiltonian still has an exact gauge structure, i.e., the description is invariant under the transformation $\vec{A}\to\vec{A}+\vec{\nabla}\Lambda(\vec{r},t)$, $c_{j}\to c_{j}e^{i\Lambda(\vec{R}_{j},t)}$, $\phi(\vec{r})\to\phi(\vec{r})-\partial_{t}\Lambda(\vec{r},t)$ when a scalar potential term $H_{\phi}=e\sum_{j,n}\phi(\vec{R}_{j})c_{jn}^{\dagger}c_{jn}$ is added, which is the requirement stated in Ref. Boykin et al., 2001. The real-world application however deserves several comments.

Some ambiguity in this truncation comes from the definition of the dipolar moment, which depends on the choice of the origin $\vec{R}_{k}$ in Eq. (51). This issue will be important for an ab-initio determination of the dipolar matrix elements, see the discussion below, which goes beyond the scope of the current paper. 2. 2.

In the current work, we will truncate the range of the dipolar matrix element to the nearest neighbor. However, this approximation is not crucial and can be relaxed. 3. 3.

It is a nontrivial question whether the dipolar light-matter Hamiltonian is optimal for a quantitative description of strongly driven phenomena beyond linear response, or whether, e.g., working with a projected continuum Hamiltonian in a fixed gauge can give quantitatively better results for a given number of bands. A general answer to this question is beyond the scope of this paper. An explicit comparison for sub-cycle driven dynamics in a single band shows that the dipolar Hamiltonian is not only favorable because of its gauge structure, but it also converges more rapidly with the number of bands taken into account.Jiajun et al. (2019) A future combination with ab-initio input and its comparison with experimental data could serve as a-posteriori justification.

As this is the first work which addresses the nonlinear response in correlated systems with both intraband and interband terms, we will focus on their qualitative contributions to physical observables.

Given $\tilde{H}_{ED}$, the corresponding currents can be obtained as $\vec{j}=\partial_{t}\mathcal{\vec{\tilde{P}}}(t)=-\mathrm{i}[\mathcal{\vec{\tilde{P}}},\tilde{H}_{ED}]=-\mathrm{i}[\mathcal{\vec{\tilde{P}}},H_{t}+H_{D}],$ where we have separated the pure Peierls term $H_{t}$ and the dipolar term $H_{D}.$ The polarization is given by $\mathcal{\vec{\tilde{P}}}=\mathcal{\vec{P}_{\text{dia}}}+\mathcal{\vec{P}_{\text{off}}}=e\sum_{i}\vec{R}_{i}n_{i}+\sum_{ij\\ nm}\vec{D}_{ij}^{nm}c_{in}^{\dagger}c_{jm}e^{-e\mathrm{i}\chi_{ij}}$. The commutator of the diagonal term $\mathcal{\vec{P}_{\text{dia}}}$ with the kinetic part of the interaction gives the usual contribution from the electron hopping

[TABLE]

or in momentum space

[TABLE]

The commutator of the diagonal term $\mathcal{\vec{P}_{\text{dia}}}$ with the dipolar part of the Hamiltonian evalutes to

[TABLE]

and expresses the change in the polarization due to the dipolar transition. In momentum space, this expression becomes

[TABLE]

The commutator of the off-diagonal polarization $\mathcal{\vec{P}_{\text{off}}}$ with the kinetic term $H_{t}$ reduces to

[TABLE]

If the dipolar matrix element is local in space this expression reduces to

[TABLE]

or in momentum space

[TABLE]

The Peierls $\vec{j}_{K}$ and the dipolar $\vec{j}_{D2}$ contributions to the current are independent of the electric field which makes them dominant at the linear-response level. The dipolar contribution to the current $\vec{j}_{D1}$ contains a factor of $\vec{E}$ and therefore it will only contribute at the nonlinear level. In the practical calculation, we evaluate the Peierls contribution to the current $\vec{j}_{K}$ by Eq. (III) and the total contribution to the dipolar current by a numerical derivative of the inter-band polarization $\vec{j}_{D1}+\vec{j}_{D2}=\partial\mathcal{\vec{P}_{\text{off}}}/\partial t.$

The dipolar matrix element is a material-specific quantity. In the case of the $d$-$p$ model we will treat only the nearest-neighbor $d$ to $p$ dipolar matrix elements and neglect the rest:

[TABLE]

where $R^{*}$ is the center of mass in the unit cell. For the dipolar matrix element between the $d$ and $p_{x}$ orbital the time-dependent phase is $e^{-\mathrm{i}eA_{x}a_{x}/2},$ since the distance between the $d$ and $p_{x}$ orbital is $a/2.$ The same expression holds for the $d$-$p_{y}$ element, but due to the phase dependence of the orbital it has the opposite sign.

To the best of our knowledge, the dipolar matrix element between the $d$ and $p$ orbital has not been reported in the literature. For the purpose of this work, we have performed the integral numerically using a cluster of $d_{x^{2}-y^{2}}$ and $p_{x},p_{y}$ orbitals, whose spatial profile is given by a Gaussian function with a variable width $\sigma$. The evaluation of the dipolar matrix element for a Gaussian with width $\sigma=0.01a,$ where $a$ is the lattice constant, yields a dipolar matrix element $|D|\approx 0.3$ $e$a, with $e$ the electronic charge. The matrix element strongly depends on the spatial width of the orbital and in the following we will use two characteristic values, namely $|D|=0.1,0.3$ $e$a. For a quantitative description of the photoexcitation in multiband materials, it will be crucial to determine these parameters from ab-initio calculations, taking into account material-specific properties.

IV Equilibrium

In the following we will focus on the parameter regime relevant for La2CuO4, as obtained from LDA calculations and the constrained RPA downfolding to the low-energy space composed of the $d$ and $p$ orbitals, with parameters $U_{dd}=5.0$ eV, $U_{dp}=2.0$ eV, $t_{dp}=0.5$ eV, $t_{dd}=-0.1$ eV, $t_{pp}=0.15$ eV, and $\Delta_{pd}=-3.5$ eV, see Refs. Werner et al., 2015; Hansmann et al., 2014 for ab-initio estimates. In all calculations, we set the inverse temperature to $\beta=5.0$ eV${}^{-1}\approx 1/2000$ K*-1*, which is above the Néel temperature. We have employed this high temperature due to the limitations of the NCA approximation, which prevent us from stabilizing calculations at lower temperatures. As the studied system has a 2 eV gap, the thermal fluctuations across the gap even at this elevated temperature can be neglected and we do not expect qualitative changes for lower temperatures until the Néel temperature is reached. We will compare results for two different approximations, namely, $GW$+EDMFT and Hartree-Fock+EDMFT (HF+EDMFT). As only $GW$+EDMFT includes nonlocal dynamical charge fluctuations, the comparison between the two methods will illustrate the role of nonlocal screening.

In Fig. 1, we present the orbital-resolved local spectral function $A_{A}(\omega)=-\frac{1}{\pi}\text{Im}[G_{\text{loc},AA}(\omega)],$ for $A=d,p_{x},p_{y}$ obtained with the two schemes. The lower Hubbard band and the $p$ orbital lie in the same energy range and are strongly hybridized due to the direct hopping $t_{dp}$. The spectral function below the Fermi level is split into three distinct peaks, see discussion in Ref. Zhang and Rice, 1988: (a) a peak at $\omega=-3.5$ eV of predominantly $p$ orbital character corresponding to the nonbonding orbital combination, (b) the antibonding band corresponding to the Zhang-Rice singlet around $\omega=-1.5$ eV, and (c) the bonding band pushed to lower energy, with a center at $\omega=-5$ eV. The Zhang-Rice band has been originally identified as an antibonding singlet state in a typical cuprate setup Zhang and Rice (1988). While our numerical approach properly captures the orbital nature of this state, the two-particle spin fluctuations would have to be included to capture the singlet nature of the wave function. A detailed analysis using a cluster extension of DMFT showed that the feedback of the spin fluctuations onto the single-particle spectrum leads to a narrowing of the Zhang-Rice band due to the polaronic effect Yin et al. (2008).

The presence of a large gap suppresses the screening in equilibrium, which explains the very similar $GW$+EDMFT and HF+EDMFT results. (We have explicitly checked that this is not anymore the case in doped systems.) The weak screening effect is also evident in a small reduction of the effective impurity interaction $\mathcal{U}(\omega=0)-U^{dd}\approx-0.05$ eV, so that the static approximation is justified.

The momentum- and orbital-resolved equilibrium spectral function $A_{k}(\omega)$, shown in Fig. 2, reveals that the band of predominantly $p$ character is split into two distinct bands around the $X$ point. This splitting is present already in the noninteracting spectrum (black lines) and originates primarily from the direct $p_{x}$-$p_{y}$ hopping, although also the hopping between the $d$ and $p$ orbitals leads to a slight splitting. The hybridization between the $d$ and $p$ band (see admixture of blue and red color in Fig. 2) is strongest between the $X$ and $M$ point and between the bonding band and the $p$ bands.

V Photo-excitation

V.1 Orbital occupation

Now we turn to the relaxation dynamics after photoexcitation. We excite the system with a short pulse

[TABLE]

with frequency $\Omega$ and amplitude $E_{0}$, polarized along the (11) direction. The width of the pulse $t_{0}=2\pi n/\Omega$ is chosen such that the envelope accommodates $n=2$ cycles, unless otherwise stated.

The electric field pulse creates charge-transfer excitations and results in an increase of the double occupancy $d_{\text{occ}}=\langle n^{d}_{\uparrow}n^{d}_{\downarrow}\rangle$ of the $d$ orbitals, as evidenced in the time evolution illustrated in Figs. 3(a) and (c). In the following, we will adjust the strength of the pulse ($E_{0}$) to fix the photodoping at time $t=12$ fs to be $\Delta d_{\text{occ}}$= 1$\%$, 2.5$\%$, and 5$\%$. The magnitude of the dipolar matrix element determines the ratio between direct interband excitations and Landau-Zener-like tunneling due to the intraband acceleration. The latter can be seen from the comparison of the contributions from the Peierls term in Eq. (III) and the dipolar term in Eq. (57) to the photoinduced current. Because the dipolar and the Peierls contribution have different interorbital matrix elements, it is expected that if one or the other of the two processes dominantes, the orbital population (density of holes in the $d$ and $p$ bands) or the photo-induced current (optical conductivity) may be different after the photo excitation. We investigate the effect of the dipolar matrix element by comparing two values, namely $|\text{D}|=0.3$ $e$a and $|\text{D}|=0.1$ $e$a, at fixed photo-doping. In order to reach the same photo-doping concentration, a much larger field strength $E_{0}$ is needed in the case of the smaller dipolar matrix element. Surprisingly, the final occupation of the orbitals for fixed photo-doping, see Figs. 3(b) and 3(d), is almost independent of the dipolar matrix element. This illustrates that the state reached after the photodoping excitation is essentially determined by $d_{\text{occ}}$ and does not depend on the details of the excitation process or the value of the dipolar matrix element.

To further confirm this finding, we investigate the effect of different pulse frequencies. In particular, we will focus on excitations from the ZRS to the UHB ($\Omega=4.0$), from the band with predominant $p$ character to the UHB ($\Omega=6.0$) and the bonding band to the UHB ($\Omega=7.5$). A similar analysis has already been presented in Ref. Golež et al., 2019, with the conclusion that an orbital selective excitation of charge carriers is not possible. Here we make this test stronger by using pulses with a larger number of cycles, $n_{p}=4$, corresponding to a narrower energy distribution of the pump pulse. The results are presented in Fig. 4, where we also provide a direct comparison to the HF+EDMFT results. As the absorption is reduced in HF+EDMFT, a stronger pulse strength has to be applied to reach a fixed number of doublons after the pulse. The excitation from the ZRS to the UHB produces pronounced oscillations in both approximations, see Figs. 4(a) and (b). For higher frequencies the inclusion of the nonlocal fluctuations strongly damps the long-lived oscillations in comparison to HF+EDMFT. Within $GW$+EDMFT the long-time orbital occupation of the $d$ orbital is almost independent of the pump pulse frequency, which further confirms the statements in Ref. Golež et al., 2019. In contrast, in HF+EDMFT, the $d$-orbital occupation exhibits a more pronounced pulse-frequency dependence.

The change in the $d$ occupancy, $\Delta n_{d}$, is related to the change in the $p$ occupancy, $\Delta n_{p}$, by $\Delta n_{d}=-2\Delta n_{p}$. The change in the double occupancy $\Delta d_{\text{occ}}$ corresponds to the total amount of charge transferred across the gap, and hence the change in the number of empty $d$-sites can be extracted as $\Delta h_{d}=\Delta d_{\text{occ}}+2\Delta n_{p}$. Since the $p$ orbitals are initially fully filled, the hole density on the $p$ orbitals (corresponding to singly occupied sites) can be estimated as $\Delta h_{p}=0.5\Delta n_{d}$, see Fig. 4 for the dynamics of both types of holes. Similarly to the dynamics of the occupation, the dynamics of the density of holes in the lower Hubbard band shows a strongly oscillatory behavior in both approximations only for the lowest frequency $\Omega=4.0$. For high-frequency excitations, see Figs. 4(c) and 4(g), the density of $d$ holes in the lower Hubbard band is small and the long-time dynamics is essentially constant within HF+EDMFT. In contrast, the $GW$+EDMFT result demonstrates a clear relaxation toward an enhanced orbital occupancy of the $d$ orbital. Due to the strong oscillations in the number of $d$ and $p$ holes in HF+EDMFT, it is hard to discuss the orbital selectivity of holes. However, it is clear that the deviation between different excitations is much larger than in $GW$+EDMFT. Within the latter approximation, the system relaxes to a state with practically identical occupations of the orbitals and similar hole concentrations for all excitation frequencies. The explicit relaxation dynamics of the holes can be monitored by the greater component, which will be discussed in the following section.

V.2 Spectral function and occupation

Now we turn to the dynamics of the spectrum and the population dynamics. The time-dependent spectral function is defined as a partial Fourier transform over the difference in time $A(\omega,t)=-\frac{1}{\pi}\text{Im}\int_{t}^{t+t_{\text{cut}}}dt^{\prime}e^{\mathrm{i}\omega(t^{\prime}-t)}G^{R}(t^{\prime},t),$ with a cutoff $t_{\text{cut}}=8$ fs. This allows us to compute the spectrum with a time-independent resolution up to relatively long times. We have checked that the precise choice of the cutoff does not qualitatively affect the dynamics of the spectral features discussed below.

In the following, we will compare the photoinduced changes of the spectrum for $GW$+EDMFT [Fig. 5(a)] and HF+EDMFT [Fig. 6(a)]. The band-gap renormalization is present in both cases, however it is much larger in $GW$+EDMFT. Moreover, the $GW$+EDMFT result exhibits strong broadening of all spectral features. We can distinguish two processes leading to the band gap renormalization: (a) the static Hartree-shift due to the Coulomb interaction $U_{dp}$ between holes in the $p$ orbitals and doublons in the $d$ orbital, $\Delta\Sigma^{H}_{dd}=(U_{dd}-2U_{dp})\Delta n_{d}$ Cilento et al. (2018); Peli et al. (2017); Lantz et al. (2017), where the factor of 2 originates from the number of nearest-neighbor $p$ orbitals and we have used the conservation of the total charge $\Delta n_{d}=-(\Delta n_{p_{x}}+\Delta n_{p_{y}})$, (b) the photoinduced enhancement of screening, which reduces the effective interaction. To demonstrate that the main contribution to the band-gap renormalization originates from the dynamical screening, we mark the static Hartree shift $\Delta\Sigma_{H}$ by a vertical red line in Figs. 5(a) and 6(a). Within HF+EDMFT, the Hartree shift perfectly describes the photoinduced changes, while for $GW$+EDMFT it only accounts for a small part of the observed band renormalization, see also the discussion in Ref. Golež et al., 2019.

Within $GW$+EDMFT, the local spectrum is also strongly broadened, see Fig. 5(a). The broadening can originate from two distinct mechanisms: (a) a change in the lifetime due to increased scattering and (b) photoinduced changes in the band dispersion (velocities). To distinguish these two effects, we analyze the change in the momentum- and orbital-resolved spectral function with respect to its equilibrium value $\Delta A_{k}(\omega,t)=A_{k}(\omega,t)-A_{k}^{\text{eq}}(\omega)$, see Fig. 7. The momentum-dependent information shows that the change in the region of the UHB depends on momentum: the spectral weight close to the $\Gamma$ point is shifting towards the chemical potential, while spectral weight near the $Z$ point is shifting away from the chemical potential [green arrows in the $d$-orbital panel of Fig. 7(a)]. A detailed comparison of the momentum cut at the $\Gamma$ point in Fig. 7(c) reveals asymmetric changes and a net shift of spectral weight towards the chemical potential at the lower edge of the UHB, see the green lines for the difference of the spectra in Fig. 7(c). The opposite trend is seen at the Z point. The latter difference could originate either from the photoinduced changes in the Fock contribution, which affects the bandwidth of the Hubbard bandAyral et al. (2017), or the $GW$ contribution originating from the momentum-dependent scattering with the plasmonic excitations. We have checked that the change in the bandwidth due to the Fock term after the excitation is at least an order of magnitude smaller than the shift in the momentum-dependent spectrum. For the UHB, the latter is presented in Fig. 7(b) as a photoinduced change in the maximum of the spectrum. Besides the change in the dispersion, the total change in the spectrum is accompanied by a broadening, which is roughly momentum independent. This leads to the conclusion that the photoinduced changes in the UHB and gap originate from the combination of both, namely a broadening of the spectrum, and changes in the dispersion.

In contrast, the Zhang-Rice singlet exhibits only a slight broadening, without any clear changes in the dispersion. Moreover, the main part of the $p$ band shows an almost perfectly rigid band shift with minor lifetime effects, in agreement with the change in the local spectrum, see black arrows in Fig. 7(a). Because the shift is rigid, the effect is also seen in the local ($k$-integrated) spectra, see Fig. 5(a). This is confirmed by the behavior of the photoinduced difference at $\omega=-5.0$ in Fig. 7(c), which shows a symmetric oscillation-like shape. An interesting observation is that the bonding band shows a dichotomy between the $d$ and $p$ contribution, where the latter exhibits a momentum dependent broadening, while the broadening of the former is almost momentum independent, compare the $d$ and $p$ spectra in Fig. 7(a). The comparison of the differences in the corresponding spectral peaks at the $\Gamma$ and Z points, see Fig. 7(c), confirms that these quasi-particles are only broadened without a clear shift.

To investigate the population dynamics we compute the time- and orbital-resolved lesser component $G^{<}(\omega,t)$, which measures the occupation dynamics as a function of the probe time $t$ and the energy $\omega$

[TABLE]

This expression is similar to the formula for the photoemission spectrum Freericks et al. (2009), but does not contain the pulse envelope functions. This simple forward integration allows us to analyze longer times.

We choose to excite the system with a pulse which is resonant to the transition between the bonding band and the UHB ($\Omega$ = 7.5 eV) to clearly resolve the hole dynamics. We have checked that the dynamics of holes is qualitatively similar for the resonant excitation between the predominantly $p$ band and the UHB. In Figs. 5(b) and 6(b), we compare the evolution of holes in the $GW$+EDMFT and HF+EMDFT approximation. The lesser (greater) component is plotted for $\omega>0$ $(\omega<0).$

After the pump pulse the occupation dynamics shows a relaxation of the doublons to the lower edge of the UHB within $GW$+EDMFT, see $G^{<}$ curves in Fig. 5(b). In contrast, the doublon distribution does not relax within HF+EDMFT, see Fig. 6(b). Interestingly, a small number of $p$ electrons is also present in the UHB due to the direct hybridization between the doublons on the $d$-orbital and the electrons on the $p$-orbitals.

In order to follow the dynamics of the holes, we will now analyze the greater component $G^{>}(\omega,t)=\frac{1}{\pi}\text{Im}\int_{t}^{t+t_{\text{cut}}}dt^{\prime}e^{\mathrm{i}\omega(t^{\prime}-t)}G^{>}(t^{\prime},t).$ In both approximations, the photo-excitation creates a substantial density of holes on the $p$ orbitals. Within $GW$+EDMFT, see Fig. 5(b), these $p$ holes quickly relax to the upper band edge. The weight of the $d$ holes increases, which is the result of relaxation processes including inter-band scattering and impact ionizationWerner et al. (2014). These two processes are responsible for a substantial increase in the density of $d$ holes in the LHB, compare Figs. 5(a) and 3. In contrast, the holes are almost trapped within HF+DMFT, see Fig. 6(b), as the interband relaxation is highly suppressed without charge fluctuations. Assuming that the pulse creates holes at similar energies in both approximations, the comparison furthermore shows that due to the nonlocal and interband charge fluctuations a significant fraction of holes already relaxes during the pulse.

The experimental renormalization of the nonbonding oxygen $p$ band has been recently reported in a time-resolved ARPES study of optimally doped Y-Bi2212 single crystals, see Ref. Cilento et al., 2018. As we are comparing a doped system in the experiment with a half-filled system in theory, a comparison close to the chemical potential is not possible. Therefore, we will focus on the dynamics of the nonbonding oxygen $p$ band close to the anti-nodal point. In the experiment, the dynamics of the oxygen band is characterized by a long-lived spectral weight variation and the comparison with Fig. 5(a) for the $p$ band close to the X point reveals a qualitatively similar behavior. While experimentally, the renormalization of the nonbonding $p$ band is only reported close to the X point, our modeling predicts an almost momentum-independent renormalization. A more detailed comparison between theory and experiment would be desirable to understand the momentum-dependent renormalization in the whole BZ.

V.3 Optical conductivity

Optical spectroscopy has been extensively used to study the photoinduced dynamics in charge-transfer insulators.Giannetti et al. (2016) The pioneering experiments on cuprates have shown a clear band gap renormalization, which at the time was attributed to photoinduced in-gap states.Matsuda et al. (1994) Subsequent experiments confirmed the effect and pointed toward a multiband origin of the gap renormalization.Okamoto et al. (2011, 2010); Novelli et al. (2014) Here we analyze the optical absorption by explicit simulations of the pump and probe pulse and extract the photoinduced current as the difference of the current with and without a probe pulse, $j_{\text{probe}}=j_{\text{pump+probe}}-j_{\text{pump}}.$ For a weak probe pulse, we evaluate the optical conductivity as the ratio

[TABLE]

where $X(\omega,t_{p})=\int_{0}^{t_{\text{cut}}}dsX(t_{p}+s)e^{-\mathrm{i}\omega s-\eta s}$ is the Fourier transform of $X=j_{\text{probe}}$ or $E,$ and $t_{p}$ is the start of the probe pulse. Due to the long-lived oscillations in the photoinduced current, a rather large broadening $\eta=0.2$ is used to avoid artifacts from the finite time-window Fourier transform. This procedure avoids the calculation of the current-current correlation function including vertex corrections.

In order to address the role of the dipolar matrix element D on the dynamics of the optical conductivity we compare results for $|\text{D}|=0.0$ $e$a and $0.3$ $e$a. In the following, we apply both the pump and probe pulses in the (11) direction. Different dipolar matrix elements D modify the equilibrium optical conductivity, compare Fig. 8(a) and 8(b). In both cases, the optical signal is composed of two characteristic features corresponding to either the transitions between the Zhang-Rice singlet and the UHB (smaller peak at $\omega\approx 4$ eV) or from the band with predominant $p$ character to the UHB. The position of the latter shows a substantial shift to lower energies as the dipolar matrix element increases. This suggest that the main contribution of the dipolar current comes from slightly lower energies than the Peierls contribution, which we have confirmed by analyzing both the intraband and interband contributions to the current (not shown). Moreover, the height of the peaks is strongly increased for $|\text{D}|=0.3$ $e$a compared to $|\text{D}|=0$ $e$a. This is easily understood as the bigger dipolar matrix element leads to a higher transition probability between bands with different orbital character. We have checked that the heights of the peaks do not simply scale with the size of the dipolar term. The orbital-resolved optical conductivity shows that the photo-induced current through the $p$ and $d$ orbitals is almost the same. This is a consequence of a small direct hopping between the $p_{x}$ and $p_{y}$ orbitals, namely, t${}_{pp}=0.15$ eV. Hence, the current is mainly originating from the hopping between $p$ and $d$ orbitals.

After the photoexcitation, the edges of these peaks are shifted to lower energies and furthermore, the peaks get broadened in analogy to the single-particle spectra. Both effects are clearly seen for the main peak for both values of the dipolar matrix element. While the height of the main peak is reduced in both cases, the effect is strongly enhanced for $|\text{D}|=0.3$ $e$a. This can be intuitively understood as a population effect since after the photoexcitation the presence of doublons in the UHB partially prevents a direct dipolar excitation. The peak at $\omega=3.5$ eV, which corresponds to transitions between the Zhang-Rice singlet and the UHB, is strongly smeared out and can only be recognized as a shoulder-like structure. The latter could be an effect of a rather large broadening $\eta=0.2$ and therefore we cannot resolve small photoinduced features associated with transitions between the Zhang-Rice singlet and the UHB after the photoexcitation. The orbital-resolved optical conductivity shows that the changes in the spectrum are comparable for all orbitals.

We extracted the photoinduced shifts for different photodopings from the change in the half-width at half-maximum for fixed dipolar matrix elements, see Fig. 8(c). The latter is compared to the static shift determined from the expression for the Hartree shift $\Delta\Sigma^{H}_{dd}=(U_{dd}-2U_{dp})\Delta n_{d}$, which only depends on the number of the photo-induced doublons. We conclude that the static shifts are generally much smaller and their increase slightly depends on the dipolar matrix element, see Fig. 8. In contrast, the total shift has a stronger dependence on the dipolar moment. Without dipolar moment ($D=0.0$), the bandgap shift with respect to the photo-doping density exhibits almost a linear dependence, while for finite dipolar moment the dependence becomes strongly nonlinear. Since the photo-induced changes for the different dipolar matrix elements $|\text{D}|=0.0$ and $0.3$ $e$a show a qualitatively and quantitatively different dynamics [compare Figs. 8(a) and 8(b)], it will be crucial to determine the dipolar matrix elements from ab-initio calculations for comparison with experiments.

VI Conclusions

We have discussed and demonstrated the nonequilibrium generalization of the multiband $GW$+EDMFT formalism by considering the Emery model relevant for La2CuO4. This approach is based on a multiscale modeling, where the low-energy correlated $d$ orbitals are treated within $GW$+EDMFT, while the higher-lying $p$ orbitals are described with $GW$ (in $GW$+EDMFT) or Hartree-Fock (in HF+EDMFT). We provided a detailed description of the downfolding formalism and the mapping to the effective impurity problem, which we solved within the noncrossing approximation (NCA).

We discussed the coupling of the electrons with classical light in the case of multiband systems. Starting from the continuum description and the minimal coupling, and after performing the Power-Zienau-Wolley transformationLoudon (2000), we projected the Hamiltonian to the Wannier orbitals. A Peierls-like transformation then led to the gauge-invariant formulation of the intraband acceleration described by a complex hopping parameter and additional interband dipolar transitions.

In the detailed analysis of the photodoped system, we addressed the three questions raised in the introduction. (i) We showed that the frequency of the pulse can lead to long-lived oscillations if the photoexcitation connects the Zhang-Rice singlet and the upper Hubbard band. For higher-frequency excitations, the systems’ properties depend only on the strength of the excitations and even if predominantly $p$-like charge carriers are excited, they quickly relax toward the band edge (Zhang-Rice singlet state). However, the latter process is only described within $GW$+EDMFT, while within HF+EDMFT the charge carriers are almost frozen within the band. This analysis implies that the dynamical charge fluctuations speed up not only the intraband relaxation but also the interband relaxation. (ii) The size of the dipolar matrix element determines the amount of energy absorbed from the high-frequency excitations on the fs timescale. However, if we fix the amount of photodoping by adjusting the pulse strength, the long-time relaxation dynamics and photoinduced changes are agnostic to the excitation protocol. (iii) The photoexcitation leads to modifications of the bands and the charge-transfer gap. In the local spectral function, the band-gap renormalization and lifetime effects can be clearly distinguished. Moreover, the momentum- and orbital-resolved spectral function reveals that the main origin of these features is a nontrivial interplay between a photoinduced change of the effective velocity, $k$-dependent broadening of the UHB and a rigid shift of the $p$ band. (iv) The charge-transfer gap renormalization is also evident in the optical conductivity. First, the size of the band shift depends on the dipolar matrix element. Moreover, the width of the main peak and the visibility of the Zhang-Rice like structure after the pulse qualitatively depends on the size of the dipolar moment.

An alternative approach to the formalism presented in this work is a recent hybrid implementation of the time-dependent density-functional theory (TDDFT+U), which has been applied to photoexcited NiO Tancogne-Dejean et al. (2018). Similar band shifts and effectively reduced interactions have been observed, but a direct comparison is difficult due to the lack of spectral information in TDDFT+U. As the authors in Ref. Tancogne-Dejean et al., 2018 showed an interesting spectrum of high-harmonics, such an analysis would provide a useful test for the proposed formalism. The photoinduced changes in the screening spectrum Golež et al. (2017, 2015) can lead to an enhanced recombination of charge carriers and consequently a larger high-harmonic response Murakami et al. (2018).

Apart from these insights into the nonequilibrium properties of charge-transfer insulators, our work represents an important step in the development of ab-initio simulation approaches for strongly correlated systems in nonequilibrium states. The $GW$+EDMFT method implemented here features a fully consistent treatment of correlation and screening effects, and can be combined with material-specific input obtained within a multitier approach analogous to the scheme recently demonstrated in equilibrium systems in Refs. Werner et al., 2015; Boehnke et al., 2016; Nilsson et al., 2017.

Acknowledgements.

D.G. acknowledges insightful discussions with Lewin Boehnke, Yuta Murakami and Michael Schüler. The calculations have been performed on the Beo04 and Beo05 clusters at the University of Fribourg and the Rusty cluster at the Flatiron Institute. DG and PW acknowledge support from SNSF Grant No. 200021_165539 and ERC Consolidator Grant No. 724103. ME acknowledges financial support from ERC Starting Grant No. 716648. The Flatiron Institute is a division of the Simons Foundation.