Position Representation of Effective Electron-Electron Interactions in Solids

TL;DR

This paper introduces a basis-independent method to compute the effective electron-electron interaction in solids, providing insights into screening effects and interactions in various materials using ab initio calculations.

Contribution

The authors develop a systematic scheme to calculate the dynamical electron-electron interaction in position space, applicable to multiple models and materials, enhancing the construction of accurate model Hamiltonians.

Findings

Screening channels vary significantly between different compounds.

The minimal models show regions with strong attractive static interactions.

Temporal interactions exhibit damped oscillations across studied materials.

Abstract

An essential ingredient in many model Hamiltonians, such as the Hubbard model, is the effective electron-electron interaction $U$, which enters as matrix elements in some localized basis. These matrix elements provide the necessary information in the model, but the localized basis is incomplete for describing $U$. We present a systematic scheme for computing the manifestly basis-independent dynamical interaction in position representation, $U({\bf r},{\bf r}';\omega)$, and its Fourier transform to time domain, $U({\bf r},{\bf r}';\tau)$. These functions can serve as an unbiased tool for the construction of model Hamiltonians. For illustration we apply the scheme within the constrained random-phase approximation to the cuprate parent compounds La$_2$CuO$_4$ and HgBa$_2$CuO$_4$ within the commonly used 1- and 3-band models, and to non-superconducting SrVO$_{3}$ within the $t_{2g}$ model. Our method is used to investigate the shape and strength of screening channels in the compounds. We show that the O 2$p_{x,y}-$Cu 3$d_{x^2-y^2}$ screening gives rise to regions with strong attractive static interaction in the minimal (1-band) model in both cuprates. On the other hand, in the minimal ($t_{2g}$) model of SrVO$_3$ only regions with a minute attractive interaction are found. The temporal interaction exhibits generic damped oscillations in all compounds, and its time-integral is shown to be the potential caused by inserting a frozen point charge at $\tau=0$. When studying the latter within the three-band model for the cuprates, short time intervals are found to produce a negative potential.

Figures

Equations

$$\delta \rho (\mathbf{r},t)=\int d{\mathbf{r}}^{\prime }d{t}^{\prime }\chi (\mathbf{r},t;{\mathbf{r}}^{\prime },t^{\prime })V_{\text{ext}}({\mathbf{r}}^{\prime },t^{\prime }),$$

$$\delta V_H(\mathbf{r},t)=\int d{\mathbf{r}}^{\prime }v({\mathbf{r}-\mathbf{r}}^{\prime })\delta \rho ({\mathbf{r}}^{\prime },t),$$

$$\delta V(\mathbf{r},t)$$

$$+\int d{\mathbf{r}}_{1}d{\mathbf{r}}_{2}d{t}_{2}v({\mathbf{r}-\mathbf{r}}_1)\chi ({\mathbf{r}}_1,t;{\mathbf{r}}_2,t_2)V_{\text{ext}}({\mathbf{r}}_2,t_2).$$

$$\delta V=(1+v\chi )V_{\text{ext}}={ϵ}^{-1}{V}_{\text{ext}}$$

$$W({\mathbf{r},\mathbf{r}}^{\prime };\tau )$$

$$+\int d{\mathbf{r}}_{1}d{\mathbf{r}}_{2}v({\mathbf{r}-\mathbf{r}}_1)\chi ({\mathbf{r}}_1,{\mathbf{r}}_2;\tau )v({\mathbf{r}}_2-{\mathbf{r}}^{\prime }).$$

$$W({\mathbf{r},\mathbf{r}}^{\prime };\omega )$$

$$+\int d{\mathbf{r}}_{1}d{\mathbf{r}}_{2}v({\mathbf{r}-\mathbf{r}}_1)\chi ({\mathbf{r}}_1,{\mathbf{r}}_2;\omega )v({\mathbf{r}}_2-{\mathbf{r}}^{\prime }).$$

$$\chi =P+Pv\chi$$

$$W=v+v\chi v=v+vPW=v+W^{\mathrm{c}}.$$

$${\chi }^{\text{RPA}}={\chi }^0+{\chi }^{0}v{\chi }^{\text{RPA}},$$

$${\chi }^0(\mathbf{r},{\mathbf{r}}^{\prime };\omega )=2\sum _{\mathbf{k}n}^{\text{occ}}\sum _{{\mathbf{k}}^{\prime }{n}^{\prime }}^{\text{unocc}}{\chi }_{n\mathbf{k},n^{\prime }{\mathbf{k}}^{\prime }}^0(\mathbf{r},{\mathbf{r}}^{\prime };\omega )$$

$${\chi }_{n\mathbf{k},n^{\prime }{\mathbf{k}}^{\prime }}^0(\mathbf{r},{\mathbf{r}}^{\prime };\omega )={\varphi }_{n\mathbf{k}}^{\ast }(\mathbf{r}){\varphi }_{n^{\prime }{\mathbf{k}}^{\prime }}(\mathbf{r}){\varphi }_{n\mathbf{k}}({\mathbf{r}}^{\prime }){\varphi }_{n^{\prime }{\mathbf{k}}^{\prime }}^{\ast }({\mathbf{r}}^{\prime })$$

$$\times \left(\frac{1}{\omega +{\epsilon }_{n\mathbf{k}}-{\epsilon }_{n^{\prime }{\mathbf{k}}^{\prime }}+i{0}^{+}}-\frac{1}{\omega -{\epsilon }_{n\mathbf{k}}+{\epsilon }_{n^{\prime }{\mathbf{k}}^{\prime }}-i{0}^{+}}\right)$$

$$W^{\text{r}}=v+v{P}^{\text{r}}{W}^{\text{r}}$$

$$W=W^{\text{r}}+W^{\text{r}}{P}^{\text{d}}W,$$

$$U(\mathbf{r},{\mathbf{r}}^{\prime };\omega )\equiv W^{\text{r}}(\mathbf{r},{\mathbf{r}}^{\prime };\omega ).$$

$${\chi }^0(\mathbf{r},{\mathbf{r}}^{\prime };\omega )=\sum _{\mathbf{k},\alpha \beta }{B}_{\alpha }^{\mathbf{k}}(\mathbf{r}){\chi }_{\alpha \beta }^{0\mathbf{k}}(\omega )B_{\beta }^{\mathbf{k}\ast }({\mathbf{r}}^{\prime }).$$

$$W^{\mathrm{c}}(\mathbf{r},{\mathbf{r}}^{\prime };\omega )$$

$${\mathcal{I}}_{\alpha }^{\mathbf{k}}(\mathbf{r})=\int d{\mathbf{r}}_{1}v(\mathbf{r}-{\mathbf{r}}_1)B_{\alpha }^{\mathbf{k}}({\mathbf{r}}_1),$$

$${\chi }_{\alpha \beta }^{\text{RPA}}(\mathbf{k};\omega )=⟨{\tilde{B}}_{\alpha }^{\mathbf{k}}|{\chi }^{\text{RPA}}(\omega )|{\tilde{B}}_{\beta }^{\mathbf{k}}⟩.$$

$$\{B_{\alpha }^{\mathbf{k}}\}$$

$$⟨B_{\alpha }^{\mathbf{k}}|{\tilde{B}}_{\beta }^{\mathbf{k}}⟩$$

$$\sum _{\alpha }|B_{\alpha }^{\mathbf{k}}⟩⟨{\tilde{B}}_{\alpha }^{\mathbf{k}}|$$

$${\mathbf{r}}_1={\mathbf{r}}_a+\mathbf{a},$$

$${\mathcal{I}}_{\alpha }^{\mathbf{k}}(\mathbf{r})={\int }_{R_a}d{\mathbf{r}}_a\sum _{\mathbf{T}}\frac{{\text{e}}^{i\mathbf{k}\cdot (\mathbf{a}+\mathbf{T})}}{|{\mathbf{r}}_a+\mathbf{a}+\mathbf{T}-\mathbf{r}|}{B}_{\alpha }^{\mathbf{k}}({\mathbf{r}}_a),$$

$$B_{\alpha }^{\mathbf{k}}({\mathbf{r}}_a)=b_{\alpha }(r_a)Y_{LM}({\hat{\mathbf{r}}}_a).$$

$$\sum _{\mathbf{T}}\frac{{\text{e}}^{i\mathbf{k}\cdot \mathbf{T}}}{|\mathbf{r}-{\mathbf{r}}_1-\mathbf{T}|}$$

$$+\gamma \sum _{\mathbf{T}}\frac{\text{erfc}(\gamma |\mathbf{r}-{\mathbf{r}}_1-\mathbf{T}|)}{\gamma |\mathbf{r}-{\mathbf{r}}_1-\mathbf{T}|}{\text{e}}^{i\mathbf{k}\cdot \mathbf{T}},$$

$${\text{e}}^{-i\mathbf{q}\cdot {\mathbf{r}}_a}=4\pi \sum _{L=0}^{\infty }(-i{)}^L{j}_L(q{r}_a)\sum _{M=-L}^L{Y}_{LM}^{\ast }({\hat{\mathbf{r}}}_a)Y_{LM}(\hat{\mathbf{q}}),$$

$${\mathcal{I}}_{\alpha }^{\mathbf{k}(1)}(\mathbf{r})$$

$$\times \sum _{\mathbf{G}}{A}_{\gamma }(q)j_L(q{r}_a)Y_{LM}(\hat{\mathbf{q}}){\text{e}}^{i\mathbf{q}\cdot \mathbf{r}}{\text{e}}^{-i\mathbf{G}\cdot \mathbf{a}}.$$

$$\frac{\text{erfc}(\gamma |{\mathbf{r}}_a-{\mathbf{r}}_{a\mathbf{T}}|)}{\gamma |{\mathbf{r}}_a-{\mathbf{r}}_{a\mathbf{T}}|}$$

$$\times \sum _{M=-L}^L{Y}_{LM}^{\ast }({\hat{\mathbf{r}}}_a)Y_{LM}({\hat{\mathbf{r}}}_{a\mathbf{T}}),$$

$${\mathcal{I}}_{\alpha }^{\mathbf{k}(2)}(\mathbf{r})=\frac{4\pi \gamma }{2L+1}{\int }_0^{R_a}d{r}_a{r}_a^2{b}_{\alpha }(r_a)$$

Full text

I Introduction

One of the most important quantities in many-electron physics is the screened Coulomb interaction between two electrons, $W$, which is a central quantity entering the Hedin equations.hedin65 Its asymptotic value ($\omega\to\infty$) equals the bare Coulomb interaction $v$, whereas its static value ($\omega\to 0$) is very much reduced compared to $v$ due to the dynamic screening of the system, embodied by the retarded response. For finite $\omega$, it becomes a complex quantity whose imaginary part can be directly related to the experimentally measured energy-loss spectra.hedin67 Many quantities and equations are intimately tied to $W$ since the electron self-energy $\Sigma$ is a functional of it. One example is Eliashberg theory of superconductivity,Eli ; Eli2 which for years has been investigated in terms of effective interactions,ashcroft and which recently was made parameter free by making use of $W$,jap just as in superconducting density functional theory.scdft0 ; scdft A quantity closely related to $W$ is the effective low-energy interaction or partially screened interaction $U$, which excludes screening from a low-energy subspace corresponding to a model Hamiltonian and may be regarded as a dynamical and non-local generalization of the Hubbard on-site repulsion.downfold ; downfold2 ; Upaper

In the position representation, $W$ and $U$ are functions of two position variables and time (or frequency): $W(\bf{r,r^{\prime}};\tau)$, $U(\bf{r,r^{\prime}};\tau^{\prime})$, but little is known about the actual shape of these functions. The focus is typically on their matrix elements in some set of orbitals, either because these are needed when calculating other quantities or because they are central objects in Hubbard-like models. However, matrix elements are basis-dependent and, since being projected quantities, do not contain complete information about the screened interaction. We therefore present a systematic scheme which allows for the computation of the position representations of the frequency-dependent $W$ and $U$, manifestly independent of any basis. This provides an unbiased tool to pin down how a suitable model can be constructed in a given periodic solid. A subsequent Fourier transform reveals the full spatiotemporal interactions $W({\bf r},{\bf r}^{\prime};\tau)$, $U({\bf r},{\bf r}^{\prime};\tau)$. A space-time point of view may furnish useful complementary insights into the physics problem at hand, like that of high-$T_{C}$ superconductivity. To illustrate the use of the developed scheme, we compute the screened interactions in the well-known high-temperature superconductor parent compounds La2CuO4 (LCO) and HgBa2CuO4 (HBCO), and for comparison in non-superconducting SrVO3, a prototype of correlated metals.

Shortly after the ground-breaking discovery of high-temperature superconductivity in the doped cuprates1986 it was realized that standard Bardeen-Cooper Schrieffer (BCS) theorybcstheory , based on electron-phonon interaction, could neither account for their elevated critical temperatures nor their anomalous and doping-dependent isotope effect.notbcs In the well-underdoped non-superconducting regime, the cuprates share an antiferromagnetic Mott insulating order caused by strong repulsion in the partially filled Cu 3$d$ band,mott and the superconducting phase emerges, as a consequence of doping, in the vicinity of a Mott transition. It was, for this reason, early pointed out that the pairing mechanism ought to be mainly of electronic or magnetic origin,andmott a viewpoint which is reinforced by the $d_{x^{2}-y^{2}}$ symmetry of the superconducting gap.scalapino ; dsym Unfortunately, despite the progress in the field of strong correlations, there is to this day still no consensus on what mechanism or, rather, interplay of mechanisms best describes this pairing.

The strong correlations of these materials explain the qualitative failure of the local density approximation (LDA), which predicts a metal for the undoped parent compounds. The deceptively simple low-energy electronic structure can be traced back to the CuO2 sheet, in which the Cu 3$d_{x^{2}-y^{2}}$ and O 2$p_{x/y}$ orbitals hybridize to form a bonding and an antibonding state.halffilled The antibonding state, which has a strong Cu 3$d_{x^{2}-y^{2}}$ weight, forms the half-filled and well-isolated narrow band across the Fermi level in LDA. Indeed, this antibonding band is commonly used to model the low-energy electrons participating in superconductivity and frequently constitutes one of the orbitals in model Hamiltonians.emery The additional low-lying oxygen $p$ bands provide a strong screening channel that causes a substantial reduction in the effective interaction.laurentium

Many pairing mechanisms have been put forward over the last three decades. Andersonanderson suggested that strong short-range repulsive interactions lead to spin-charge separation and that the immense antiferromagnetic superexchange opens up a $d-$wave spin gap, which by kinetic frustration converts to a superconducting gap. The charge fluctuation mechanism dates back to Kohn and Luttinger,KL who realized that Friedel oscillations lead to anisotropic pairing in an isotropic electron gas with short-range interactions at low temperatures. Numerical studies within the random-phase approximation (RPA) by Rietschel and Sham later confirmed this for a certain range of electron densities by solving the Eliashberg equation.RSham Since spin fluctuations are believed to completely overshadow charge fluctuations at short distances, the latter has not been extensively investigated for the cuprates. It is conceivable that the electron gas results persist in realistic materials, but that the relevant length scale is significantly reduced. Indeed, Kohn and Luttinger argued that a non-spherical Fermi surface can drastically increase $T_{C}$.KL The screened interaction in position representation may furnish a physical insight into this mechanism, not easily accessible from matrix elements alone.

For the undoped cuprates we consider the famous one- and three-band models and calculate the effective interactions $U_{1}$ and $U_{3}$ in the respective low-energy subspace. The metallic band with dominating Cu 3$d_{x^{2}-y^{2}}$ weight constitutes the one-band subspace, whereas the three-band subspace also includes two bonding and non-bonding bands of mainly O 2$p_{x,y}$ character.emery $U$ does not include the screening of the electrons of the subspace, hence also the screening from the pathological metallic band is excluded, which partly justifies the use of LDA as a starting point. It is worth noting that the charge gap in LCO, which is absent at the LDA level, is opened up within LDA+DMFT when a dynamic $U$ computed using constrained RPA (cRPA) is used, whereas when the static value is used the material remains metallic.laurentium The measured gap of 2 eV is almost perfectly reproduced in the three-band model and partly so in the one-band model,laurentium which shows that $U$, when calculated within cRPA, indeed embodies dynamical correlation effects required when modeling the undoped cuprates. We also calculate the fully screened interaction $W$ although its interpretation demands some caution. With some justification, it may be thought of as a crude estimation of the screened interaction of the metallic doped system, which could be systematically improved, for instance, by imposing rigid shifts in the LDA band filling.fermi

This paper is organized as follows: In Section II we summarize the theory of the partially and fully screened Coulomb interaction, $U$ and $W$, as well as the RPA and constrained RPA approximations. In Section III the space-time computation of $W({\bf r},{\bf r}^{\prime};\tau)$ and $U({\bf r},{\bf r}^{\prime};\tau)$ is described, and their interpretations are emphasized. In Section IV the results for SrVO3, LCO and HBCO are presented and discussed and in Section V the main findings are summarized.

II Screened Interaction

II.1 $W$ and RPA

Before describing the position-space computation of $W(\mathbf{r,r}^{\prime};t-t^{\prime})$ or $W(\mathbf{r,r}^{\prime};\omega)$ we recapitulate the definition of $W$ from linear response theory. When applying an arbitrary external perturbation $V_{\text{ext}}(\mathbf{r},t)$ the induced density is to linear order given by

[TABLE]

where $\chi$ is the linear density response function. This causes a change in the Hartree potential

[TABLE]

which screens the applied perturbation $V_{\text{ext}}$. The resulting change in the total potential $\delta V=V_{\text{ext}}+\delta V_{H}$ is given by

[TABLE]

Schematically we may write

[TABLE]

where we recognize that $1+v\chi$ is the inverse dielectric matrix $\epsilon^{-1}$. If we replace our external perturbation with the Coulomb interaction $v(\mathbf{r-r}^{\prime})\delta(t-t^{\prime})=\delta(t-t^{\prime})/|\mathbf{r}-\mathbf{r}^{\prime}|$, with $(\mathbf{r}^{\prime},t^{\prime})$ treated as a parameter, we arrive at

[TABLE]

This is the definition of the screened interaction in the Hedin equations.hedin65 The second term, $v\chi v$, which is the screening contribution to $W$, is usually denoted by $W^{\mathrm{c}}$, a notation we will adopt in the following. We have made use of the fact that $\chi$ depends only on relative time $\tau=t-t^{\prime}$ for a system with time-independent Hamiltonian. $W({\bf r},{\bf r}^{\prime};\tau)$ is the effective interaction between two electrons at ${\bf r},t$ and ${\bf r}^{\prime},t^{\prime}$ and contains a retarded contribution, $W^{\mathrm{c}}$, due to the dynamic response of all electrons in the system. Within RPA, this retarded response originates from successive particle-hole excitations caused by the instantaneous interaction between the electrons in the system. The Fourier component of the screened interaction is then calculated from the following equation:

[TABLE]

The screened Coulomb interaction, $W$, is uniquely determined by the linear density response function $\chi=\delta\rho/\delta\varphi$. We can introduce the irreducible polarization propagator $P$, which may be thought of as the linear density response function with respect to the total field, $P=\delta\rho/\delta V$. It then follows from the chain rule that

[TABLE]

and

[TABLE]

In the random-phase approximation (RPA) the polarization propagator is approximated by the response function of a noninteracting system $\chi^{0}$,hedin65 so that the response function takes the form

[TABLE]

where

[TABLE]

[TABLE]

is equivalent to the well-known Lindhard formula.lindhart Here, $\phi_{n{\bf k}}$ and $\varepsilon_{n{\bf k}}$ are paramagnetic eigenfunctions and eigenenergies, typically obtained using density functional theory (DFT). ${\bf k}$ is restricted to the first Brillouin zone. The factor of two is due to summing over the two identical spin contributions, and the two sums are restricted to occupied (occ) and unoccupied (unocc) states respectively. Note that Eq. (10) describes the time-ordered polarization function, which means that the resulting screened interaction is not the retarded, but the time-ordered $W$. One can recover the retarded $W$ by multiplying the imaginary part of the time-ordered $W$ by a factor of $\text{sign}(\omega)$.

A qualitative and simplified depiction of $W$ is presented in Fig. 1. By increasing the depth of the screening hole, the effective interaction is reduced and can even turn negative at certain distances. The terms attraction and repulsion, however, have to be used with caution since they originate from situations where the interaction is radially monotonous and thus either attractive or repulsive throughout. Still, we adopt the term attraction if we, for a given $\mathbf{r}^{\prime}$, identify a negative minimum of the interaction at ${\bf r}$ (local attraction) towards which the classical force field is pointing.

As can be seen from Fig. 1, at very short distances to ${\bf r}^{\prime}$, the force field is always pointing outwards, which gives a local repulsion. This can be understood intuitively since for $\mathbf{r}\rightarrow\mathbf{r}^{\prime}$ there is not sufficient charge in the region between $\mathbf{r}$ and $\mathbf{r}^{\prime}$ to create screening holes that could compensate or overcompensate the Coulomb repulsion. The screening inherently depends on the electron density in the solid. Different materials will have different screening properties and therefore also different shapes of $W(\mathbf{r},\mathbf{r}^{\prime})$. The placement of the point charge will therefore also matter. If it is put at the position of a nucleus, especially of an atomic species which is an effective ”screener”, a much more reduced $W$ emerges at short distances than from a point charge in between two nuclei.

II.2 $U$ and cRPA

To determine the effective interaction of a low-energy model, we use the cRPA method,downfold ; crpa1 in which the Hilbert space is divided into a low- and a high-energy subspace, $\mathcal{D}$ and $\mathcal{R}$. The polarization function is now decomposed into two terms, $P\!=\!P^{\mathrm{d}}\!+\!P^{\mathrm{r}}$. $P^{\mathrm{d}}$ describes polarization processes within the low-energy subspace $\mathcal{D}$ whereas $P^{\mathrm{r}}$ accounts for the rest of the polarizations, i.e., those within the $\mathcal{R}$ subspace as well as those between the subspaces. By defining

[TABLE]

it can be shown thatdownfold

[TABLE]

which allows us to interpret $W^{\text{r}}$ as the effective ”bare” interaction in $\mathcal{D}$, a non-local and dynamical generalization of the Hubbard on-site repulsion.downfold2 So,

[TABLE]

As in the case of $W$, we can write $U\!=\!v\!+\!U^{\text{c}}$, where $U^{\text{c}}\!=\!v\chi^{\text{r}}v$ and $\chi^{\text{r}}\!=\!P^{\text{r}}\!+\!P^{\text{r}}v\chi^{\text{r}}$. The low-energy subspace in the Hubbard model usually corresponds to a narrow band with strong correlations, so RPA is not expected to work well. However, when computing $U$ for the model, the polarization channels within the low-energy subspace are removed from Eq. (10), so that it is justifiable to constrain the RPA to compute $P^{\text{r}}=\chi^{\text{r}0}$.

The physics lies in the choice of the low-energy model subspace. For the low-energy bands of the cuprates, which are entangled, we use the ”disentanglement” schemed developed in Ref. (29) and define the $\mathcal{D}$ subspace in terms of maximally localized Wannier functionsmaxloc and the $\mathcal{R}$ subspace as the orthogonal space. Computational details for the calculation of $U$ in the cuprates and in SrVO3 are provided in App. A.

III Position Representation

This section deals with the computation of $W$ in position representation (Eq. (6)) and its interpretation in time domain. Any expression for $W$ has an analogue for $U$ obtained by replacing $\chi^{0}$ with $\chi^{\text{r}0}$.

III.1 Product Basis

To expand the polarization $\chi^{0}$ and response function $\chi^{\text{RPA}}$ we need a set of two-particle basis functions in the form of a product basis $\{B_{\alpha}^{\mathbf{k}}\}$. This basis can be tailored to give a complete representation of $\chi^{0}$ and be optimized such that a minimal number of basis functions is neededgunnarsson ; gunnarsson2

[TABLE]

From $\chi^{\text{RPA}}=\chi^{0}+\chi^{0}v\chi^{0}+...$, it is clear that the product basis is also complete for representing $\chi^{\text{RPA}}(\mathbf{r},\mathbf{r}^{\prime};\omega)$, since $v$ is always sandwiched between two $\chi^{0}$ so that it is immaterial whether the product basis is complete or not for $v$. In other words, only the projection of $v$ in the subspace of $\chi^{0}$ is needed. In fact, the product basis constructed for $\chi^{0}$ is in general far from complete for representing $v(\mathbf{r}-\mathbf{r}^{\prime})$. Since $W=v+v\chi^{\text{RPA}}v$ within RPA, this implies that the product basis in general cannot be used for a complete representation of $W(\mathbf{r},\mathbf{r}^{\prime};\omega)$. The way around this problem is explained in the following.

III.2 $W$ and $U$ in Position Space

Figure 2 shows the steps involved to obtain the matrix elements $\chi_{\alpha\beta}^{\text{RPA}}(\mathbf{k};\omega)$ within RPA and $\chi_{\alpha\beta}^{\text{r},\text{RPA}}(\mathbf{k};\omega)$ within cRPA. These matrix elements together with the product basis completely determine $W({\bf r},{\bf r}^{\prime};\omega)$. Since $W$ partly consists of the bare Coulomb interaction $v$, which is known analytically, it is sufficient to find an expression for $W^{\mathrm{c}}$. Schematically, if we let matrix elements be underlined, Eq. (14) reads $\chi^{0}=B{\underline{\chi}}^{0}B^{*}$. Similarly, within RPA it also holds that $\chi^{\text{RPA}}=B\underline{\chi}^{\text{RPA}}B^{*}$, which implies, together with Eq. (6), that $W^{\mathrm{c}}=v\chi^{\text{RPA}}v=(vB)\underline{\chi}^{\text{RPA}}(vB)^{*}$. We have now obtained a basis which is complete for $W^{\mathrm{c}}$. Explicitly,

[TABLE]

where

[TABLE]

(15)-(17) are the main equations for obtaining $W$ in position representation. In general, the set of functions $\{\tilde{B}_{\alpha}^{\bf k}\}$ is biorthogonal to the set $\{B_{\alpha}^{\bf k}\}$ and fulfills Eqs. (18)-(20).

After having obtained all matrix elements $\chi_{\alpha\beta}^{\text{RPA}}({\bf k};\omega)$, what remains is to calculate the basis-dependent integrals $\mathcal{I}_{\alpha}^{\bf k}({\bf r})$ as well as including the $\bm{\Gamma}$-point contribution in a suitable way. We will explain both steps in the following, but first we present the product basis, constructed in the SPEX code, which has been used in this work.

III.2.1 Mixed Product Basis

The mixed product basis is an extension of the optimized product basis within the full-potential linearized augmented plane-wave (FLAPW) method,takao ; FLAPW where space is separated into spherical ”muffin tin” (MT) spheres around each atom as well as the ”interstitial region” (IR), which constitutes the remaining region of space. In the MT spheres the product basis functions $B_{aLMP}^{\mathbf{k}}(\mathbf{r})=b_{aLP}(r)Y_{LM}(\hat{{\bf r}})$ are constructed from products of the MT functions of the LAPW basis. Here, $a$ is an orbital index, $L$ and $M$ denote the orbital and magnetic quantum numbers, respectively, and $P$ is an index for different radial functions. In the IR, products of plane waves, which are themselves plane waves $B_{\mathbf{G}}^{\mathbf{k}}(\mathbf{r})=\text{e}^{i(\mathbf{k}+\mathbf{G})\cdot\mathbf{r}}/\sqrt{\Omega}$ are constructed, where $\Omega$ is the unit cell volume. The resulting ”mixed product basis” functionsFLAPW

[TABLE]

are either non-zero only in the MT spheres or in the IR. Eq. (20) holds in the subspace of $\chi^{\text{RPA}}$.

III.2.2 Muffin-tin Contribution

We start our position space reconstruction by considering the MT spheres, where $\alpha=a,L,M,P$. By defining

[TABLE]

where ${\bf r}_{a}$ is confined to a MT of radius $R_{a}$ and ${\bf a}$ is the vector pointing to the atomic centre of $a$, Eq. (16) can be re-expressed as

[TABLE]

Here we made use of Bloch’s theorem and the sum runs over all lattice vectors $\{\mathbf{T}\}$. However, $\mathcal{I}_{\alpha}^{{\bf k}}({\bf r})$ does not converge for a finite sum over ${\bf T}$ due to the long-range integrand, so we perform Ewald summation to resolve this issue (red box in Fig. 2). For ${\bf k}\neq{\bf\Gamma}$ and with $\mathbf{q}=\mathbf{k}+\mathbf{G}$, where $\mathbf{k}$ is restricted to the first Brillouin zone and $\mathbf{G}$ is a reciprocal lattice vector, Ewald’s formula readsewald

[TABLE]

For $\gamma\to 0$, the real-space sum is recovered, and, for $\gamma\to\infty$, the second term vanishes and the real-space sum is replaced by a summation in reciprocal space. For a properly chosen $\gamma$, however, the expression is short-ranged in both $|\mathbf{r}-\mathbf{T}|$ and $q$.

We separate $\mathcal{I}_{\alpha}^{\mathbf{k}}(\mathbf{r})$ into $\mathcal{I}_{\alpha}^{\mathbf{k}(1)}(\mathbf{r})$ and $\mathcal{I}_{\alpha}^{\mathbf{k}(2)}(\mathbf{r})$, resulting from the sums over $\mathbf{G}$ and $\mathbf{T}$ respectively. We define $A_{\gamma}(q)\equiv(4\pi/\Omega)\text{exp}(-q^{2}/4\gamma^{2})/q^{2}$ and make a plane-wave expansion in spherical harmonics

[TABLE]

where $j_{L}$ are the spherical Bessel functions. This yields for the first term

[TABLE]

Introducing $\mathbf{r}_{a{\bf T}}\!=\!\mathbf{r}\!-\!\mathbf{a}\!-\!\mathbf{T}$, the second term, $\mathcal{I}_{\alpha}^{\mathbf{k}(2)}(\mathbf{r})$, diverges if $|\mathbf{r}_{a}-\mathbf{r}_{a{\bf T}}|\rightarrow 0$. To resolve this issue we make use of the expansion

[TABLE]

where $r_{<}=\min(r_{a},\!r_{a{\bf T}})$ and $r_{>}=\max(r_{a},\!r_{a{\bf T}})$. Note that the majority of the terms, corresponding to translations ${\bf T}$ that cause no divergence, can be integrated without the use of this expansion. For brevity, we here keep the expansion in all terms, and arrive at

[TABLE]

The coefficients $g_{L}$ are computed as

[TABLE]

using Gaussian integration, meaning that any angular integral $\int d\Omega f(\Omega)$ is replaced by $\sum_{i}w_{i}f(\Omega_{i})$ where the weights $w_{i}$ are tabulated and independent of $f$. In particular, we used 114 cubic directions $\Omega_{i}$, which yields exact results for angular momentum components $L\leq 15$.ferdibook

III.2.3 Interstitial Contribution

We now consider the IR, where $\alpha={\bf G}$. By extending $B_{\mathbf{G}}^{\mathbf{k}}(\mathbf{r})=\text{e}^{i\mathbf{q}\cdot\mathbf{r}}/\sqrt{\Omega}$ to all of space and subtracting the muffin-tin contribution, we can write

[TABLE]

where we have made use of the fact that $B_{\mathbf{G}}^{\mathbf{k}}(\mathbf{r}_{a}+\mathbf{T})=\text{e}^{i\mathbf{k}\cdot\mathbf{T}}B_{\mathbf{G}}^{\mathbf{k}}(\mathbf{r}_{a})$. The first term reads

[TABLE]

We divide the rest into $\mathcal{I}_{\mathbf{G}}^{\mathbf{k}(1)}(\mathbf{r})+\mathcal{I}_{\mathbf{G}}^{\mathbf{k}(2)}(\mathbf{r})$ from both terms in the Ewald summation in the same manner as before, and analogously we obtain

[TABLE]

where ${\bf q}^{\prime}={\bf k}+{\bf G}^{\prime}$. Terms in $\mathcal{I}_{\mathbf{G}}^{\mathbf{k}(1)}$ with $L>4$ are very small and excluded in this work.

III.2.4 ${\bf\Gamma}$-Point Contribution

What is left at this point is to calculate the ${\bm{\Gamma}}$-point contribution to Eq. (15), which requires special treatment since the bare interaction $v$ diverges as $1/k^{2}$ for $k\!\to\!0$. In SPEX, the divergence is treated analytically by rotating to the Coulomb eigenbasisCoulBas

[TABLE]

When $k\to 0$, $E_{\mu=1}^{{\bf k}}(\mathbf{r})\to 1/\sqrt{\Omega}$ corresponds to the divergent eigenvalue of $v$ and the matrix element $W_{\mu=1,\nu=1}^{\mathrm{c}}({\bf k};\omega)$, which diverges like $1/k^{2}$, just shifts $W^{\mathrm{c}}(\mathbf{r},\mathbf{r}^{\prime};\omega)$ uniformly to leading order.FLAPW $W_{\mu=1,\nu>1}^{\mathrm{c}}({\bf k};\omega)$ and $W_{\mu>1,\nu=1}^{\mathrm{c}}({\bf k};\omega)$ diverge only like $1/k$ and are much smaller and, for this reason, neglected in this work. This simplification corresponds to making $W^{\mathrm{c}}$ block diagonal in the Coulomb basis. The large block $W_{\mu>1,\nu>1}^{\mathrm{c}}({\bf k};\omega)$ does not contain any divergence, and we therefore rotate it back to the mixed product basis. We then get the ${\bf\Gamma}$-point contribution to $W^{\mathrm{c}}$:

[TABLE]

where

[TABLE]

$\tilde{\mathcal{I}}_{\alpha}^{\bf 0}$ is calculated in the same way as $\mathcal{I}_{\alpha}^{\bf k}$. Because of the divergent behavior of $W^{\mathrm{c}}\sim 1/k^{2}$, the Brillouin-zone integration cannot be approximated by a finite summation as in Eq. (15). Therefore, we have replaced the ${\bf k}$ sum by an integral $\int_{\bf\Gamma}$, which could be understood as an integration over a finite region around ${\bf k}={\bf 0}$. In practice, we use instead an integration over the whole reciprocal space, not of $1/k^{2}$ (which would yield infinity), but of $\text{e}^{-\epsilon k^{2}}/k^{2}$ with a small positive coefficient $\epsilon$, and subtract a double-counting correction given by the sum over the ${\bf k}$-point set excluding the $\Gamma$ point. For details, see Ref. (31) and in particular Eq. (34) therein.

III.3 $W$ and $U$ in Time Domain:

Impulse and Step Response

It is interesting to study the retarded interaction both related to the impulse response and the step response of a solid. The former is to linear order given by $W({\bf r},{\bf r}^{\prime};\tau)$, and we show below that the latter is accessible from the same quantity.

The interpretation of $W$ is provided in Sec. II.1. Since it was obtained from linear response theory by replacing the external potential with the instantaneous Coulomb interaction, $v({\bf r}-{\bf r}^{\prime})\delta(\tau)$, we here denote it by $W[\delta]$. $W[\delta]$ is connected to the impulse response of the system, and is obtained by a simple inverse Fourier transform of $W(\omega)$:

[TABLE]

$W(\omega)$ is here assumed to be retarded, but the $W(\omega)$ described in Sec. II.1 is time-ordered. For positive frequencies the time-ordered and retarded $W(\omega)$ are identical, but the former is an even function

of $\omega$ whereas the latter only has an even real part, but an odd imaginary part. By only calculating $W(\omega)$ for positive frequencies, the correct symmetries can easily be imposed.

As is also clear from Sec. II.1, if we instead introduce a point charge at ${\bf r}^{\prime},t^{\prime}$ kept frozen at later times, which means inserting $v({\bf r}-{\bf r}^{\prime})\Theta(\tau)$ into Eq. (3), the resulting screened potential $W[\Theta]$ is given by

[TABLE]

Here, $\chi$ is the retarded response function, which is related to its time-ordered counterpart in the same way as described above for $W$. Since the retarded $\chi$ fulfills causality, the upper limit of integration can be changed to $\tau_{2}=\tau$, and from the variable substitution $\tau^{\prime}=\tau-\tau_{2}$ we arrive at

[TABLE]

This equation establishes a connection between the dynamically screened interaction between two electrons of the intrinsic system (impulse response) and the dynamically screened potential

from an impurity added to the system (step response). It has the following limits

[TABLE]

$W[\Theta]$ has dimension energy while $W[\delta]$ has dimension energy/time.

IV Results

We will now apply our method to compute the position representation of $W$ and $U$ in LCO, HBCO and non-superconducting SrVO3. Computational details are provided in App. A. We focus on the cases with $\mathbf{r}^{\prime}$ at the transition metal nucleus (Cu or V) as well as at the O nucleus, and with $\mathbf{r}$ and $\mathbf{r}^{\prime}$ restricted to the same CuO2 or VO2 sheet. Furthermore, in all calculations, ${\bf r}$ and ${\bf r}^{\prime}$ belong to the same unit cell.

IV.1 Static $U$ in Position Space

We start by considering the static effective interaction $U(\mathbf{r,r}^{\prime};\omega\!=\!0)$ (Fig. 3 - 5). We study the 1-band and 3-band models for the cuprates and compare the results with the non-superconducting perovskite SrVO3 in the $t_{2g}$-model (see App. A).

An interesting finding, with ${\bf r}^{\prime}$ at the transition metal nucleus, is that the $t_{2g}$ interaction in SrVO3 is essentially positive in the entire unit cell while in both cuprates there is a region close to the Cu site where $U_{1}$ ($U$ of the one-band model) is significantly negative. This region, as illustrated in Fig. 4, has a shape which originates mainly from the 3$d_{x^{2}-y^{2}}$ orbital ($x^{2}-y^{2}$-derived) of the one-band subspace even though the intra-band screening from this orbital is excluded in the one-band model. Such a region does not appear in $U_{3}$ ($U$ of the three-band model) and thus originates from the hybridization between the Cu 3$d_{x^{2}-y^{2}}$ orbital and the O 2$p_{x}$ and 2$p_{y}$ orbitals. Since the $d$ orbitals are localized, this hybridization is expected to be strong only in their vicinity, which is consistent with the shape of the attractive region in $U_{1}$. However, while the $d_{x^{2}-y^{2}}$ orbital is antisymmetric with respect to a reflection of ${\bf r}=(x,y,0)$ across the line $x=y$, the attractive region is symmetric. This is physically clear, and can be understood from Eq. (15). If we let $\mathcal{R}$ be the reflection across $x=y$, we get

[TABLE]

Since

[TABLE]

and

[TABLE]

it follows that

[TABLE]

A striking difference can be seen between the cuprates in the one-band model (Fig. 4) and SrVO3 in the $t_{2g}$ model (Fig. 5). As already pointed out, in the cuprates, the region with strong one-band attraction coincides with the region with a large one-band density, which means that the electrons could feel the attraction. In SrVO3, on the other hand, the region with the modest attraction in the minimal (three-band) $t_{2g}$ model, does not coincide with the region of the important in-plane $xy$ orbital of the model. This means that the electrons most likely experience repulsion. This finding is backed by earlier workfirstHubb on the screening channels that determine $U_{3}$ in SrVO3, where if was found that O ${2p}-$V ${e_{g}}$ transitions constitute a stronger channel than O ${2p}-$V $t_{2g}$ transitions.

It is also worth stressing, with ${\bf r}^{\prime}$ at the Cu site, the negative $U_{1}$ at the next-nearest Cu site in both cuprates. The attraction is the strongest in HBCO, for which it survives in the three-band model. HBCO is also the only compound which displays attraction, though weak, at the neighboring Cu site (in the one-band model). The corresponding $t_{2g}$ interaction in SrVO3 at the nearest or next-nearest neighbor V site is significantly positive. When $\mathbf{r}^{\prime}$ is moved to the O site the only identified attraction is very weak and found in the one-band model of HBCO at the next-nearest Cu site as can be seen in Fig. 4.

The matrix elements of the static $U_{1}$ in the maximally localized Wannier orbitals are positive for both cuprateslaurentium ; yang16Direct but the observed region between the Cu and O sites with large negative $U_{1}$ opens up a possibility of having negative matrix elements in some other orbitals, with a large weight in the attractive region. It is conceivable that such a basis could be used to describe possible Cooper pairs derived entirely from charge fluctuations. Such a basis cannot be found in non-superconducting SrVO3 since the $U$ of the $t_{2g}$ model is almost entirely positive, at least in the first unit cell.

In Sec. IV-C we analyze the screening channels associated with Cu 3$d_{x^{2}-y^{2}}-$3$d_{x^{2}-y^{2}}$ as well as O 2$p_{x,y}-$Cu 3$d_{x^{2}-y^{2}}$ transitions, but first we discuss the fully screened interaction $W$.

IV.2 Static $W$ in Position Space

$W$ contains all screening channels of the system, including, in the case of the cuprates, the spurious metallic screening due to the pathological LDA band structures. The physical meaning of $W$ in this case is therefore not entirely clear. With this caveat in mind, it is nevertheless instructive to compute $W$ to understand the role of the screening within the antibonding band crossing the Fermi level, which may be thought of as modeling the screening of the doped system.

In Fig. 6 and 7 we compare $W({\bf r},{\bf r}^{\prime};\omega=0)$ in the CuO2 sheets of the cuprates with that of the VO2 sheet of SrVO3. When choosing ${\bf r}^{\prime}$ at the Cu or V site, large regions appear with negative $W$ in all of the compounds, but with a larger magnitude in the cuprates than in SrVO3 (-6 versus -3 eV). This can be understood by observing that in the case of the cuprates, $W$ is obtained by screening $U_{1}$ with Cu 3$d_{x^{2}-y^{2}}-$3$d_{x^{2}-y^{2}}$ transitions, which have the same shape as $U_{1}$ itself. The screening in the $x^{2}-y^{2}$ channel is thereby enhanced. In SrVO3, on the other hand, the screening in the $xy$ channel essentially only originates from within the $t_{2g}$ subspace, since there are no close-by orbitals outside the subspace to hybridize with.

To investigate whether the $x^{2}-y^{2}$-derived shape of $W$, that can be seen in Fig. 7, is consistent with a superconducting gap with $x^{2}-y^{2}$ symmetry, we consider the superconducting DFT (SCDFT) gap equation. When excluding the effect of phonons, the SCDFT gap equation contains only the Kohn-Sham eigenenergies $\varepsilon_{n{\bf k}}$ and the static $W$ and readsscdft0 ; scdft

[TABLE]

where $\mathcal{E}_{n{\bf k}}=((\varepsilon_{n{\bf k}}-\mu)^{2}+|\Delta_{n{\bf k}}|^{2})^{1/2}$. Furthermore, if $W$ or $\Delta$ have certain symmetries under a unitary transformation $S$ in position representation, this holds analogously in reciprocal space. Since we are only interested in the symmetry of the gap $\Delta$, we simplify the equation by linearizing it around $T=T_{C}$, where $\Delta_{n{\bf k}}$ is small. Since the ratio $\tanh(\frac{\beta}{2}\mathcal{E}_{n^{\prime}{\bf k}^{\prime}})/\mathcal{E}_{n^{\prime}{\bf k}^{\prime}}$ is a quickly decaying function we only keep the diagonal matrix element of $W$ from the band that crosses the Fermi level. We can then drop the band index completely and obtain

[TABLE]

The symmetry can now be deduced by considering a 3$\times$3 ${\bf k}$ mesh, corresponding to the first Brillouin zone, for which we make the posteriori ansatz

[TABLE]

where the mid element corresponds to the $\Gamma$ point. By inserting this ansatz in (46) and recalling that $\beta\approx 1/k_{B}T_{C}$ the relation

[TABLE]

is obtained. Note that $T_{C}$ is the critical temperature obtained from $W$, which in general is smaller or equal to the true critical temperature, depending on what correlations are included (plasmons in this work). Since the $\Gamma$-point contribution, $a$, is in general large and positive for $W$ this relation means that a nonzero $\Delta$ is possible only for sufficiently negative $b$. This simplified condition should be applicable also to the spin-fluctuation mechanism. Equation 48 confirms that the calculated shape of $W$ is consistent with a superconducting gap of $x^{2}-y^{2}$ symmetry. A similar ansatz could be made in the $xy$ channel for SrVO3, and it is plausible that the equivalent condition is not fulfilled since the strength of attraction in the $xy$ channel in SrVO3 (Fig. 7) is only half that found in the $x^{2}-y^{2}$ channel in the cuprates. An unfulfilled condition implies that $\Delta$ is zero throughout, which obviously is true for non-superconducting SrVO3.

IV.3 Screening Channels in Position Space

Different polarization channels enter $\chi^{\text{RPA}}$ in a non-linear fashion. With the definition that the ”$pd$ screening” comes from all terms in $\chi^{\text{RPA}}$ which contain O 2$p_{x,y}-$Cu 3$d_{x^{2}-y^{2}}$ transitions to linear order or higher, the resulting contribution to the effective interaction is exactly $U_{1}-U_{3}$ (Fig. 8 and 9). In the same manner, $W-U_{1}$ (Fig. 8 and 10) is the contribution from the ”$dd$ screening”. However, the Cu $d_{x^{2}-y^{2}}$ band in the 1- and 3-band models are not exactly identical. For this reason, in the computation of $U_{1}-U_{3}$, we calculate not only $U_{3}$ but also $U_{1}$ from the 3-band interpolation.

In agreement with earlier studies of LCO, laurentium the $pd$ screening has most of its weight at the Cu site. It is clear from Fig. 8 that the metallic $dd$ screening is stronger and has longer range than the $pd$ screening. The striking similarity between the results for LCO and HBCO indicate that the screening of the cuprates is generic, although the actual strength is material specific.

IV.4 $W$ and $U$ in Time Domain

The screened interaction $W({\bf r},{\bf r}^{\prime};\tau)$ in time domain ($W[\delta]$ in Sec. III.3) is presented in Fig. 11 together with $W^{\mathrm{c}}({\bf r},{\bf r}^{\prime};\omega)$ for LCO, HBCO and SrVO3, with $\mathbf{r}$ and $\mathbf{r}^{\prime}$ at the same transition metal nucleus (Cu or V).

$W$ shares a common characteristic feature in time domain in all compounds. Shortly after the instantaneous bare interaction, there is a sudden surge of screening holes, which causes the large dip seen in $W$. $W$ then starts to oscillate, with a dominating characteristic frequency corresponding to the main collective charge excitation (plasmon) of the system. This is superimposed by oscillations with different frequencies, corresponding to subplasmons of the system. Gradually, the oscillations decay and almost vanish after 2000 attoseconds. This can be understood by considering the simple model ($\omega_{n}>0$)

[TABLE]

where the imaginary part is assumed to be a series of sharp $\delta$-functions, each representing a subplasmon excitation with an appropriate weight $W_{n}>0$. Inverse Fourier transformation leads to

[TABLE]

The behavior of $W^{\mathrm{c}}(\tau)$ for small $\tau$ is governed by the high-frequency features of $W^{\mathrm{c}}(\omega)$ and the dominating oscillation is determined by the bulk plasmon of the system. This explains the similar behavior for small $\tau$ in all the compounds in Fig. 11 since the high-frequency electron gas-like bulk plasmon is usually present in real materials. Subplasmons of lower frequencies, on the other hand, are rather material specific and determine the behavior of $W^{\mathrm{c}}(\tau)$ at large $\tau$. Indeed, in the time window between 1000 and 2000 attoseconds, $W^{\mathrm{c}}(\tau)$ still displays dramatic oscillations with strong attraction in both cuprates (mainly HBCO), but not in SrVO3.

In Fig. 12 we display the behavior of $W$ and $U_{3}$ in time-domain when an impurity is added to the system at $t=0$ and then left frozen at its position (see Sec. III.3). As should be the case, the long-time limits equal the static ($\omega=0$) values of $W$ and $U_{3}$. $U_{3}$ is presented, but not $U_{1}$, because the static limit of the former is positive, whereas the static limit of $U_{1}$ is negative, just like that of $W$. The result for $U_{3}$ brings to light the presence of time intervals with a negative interaction, despite the static limit being positive. This shows the relevance of taking into account frequency dependence when utilizing $W$ or $U$ to model superconductivity.

V Summary and conclusions

We have presented a method for computing the position representation of the effective electron-electron interaction $U$ in real materials and generalized the picture in time domain to include the study of static impurities. This basis-independent space-time approach is complementary to matrix element studies and allows for an unbiased perspective on the screening in real materials. This can be used to construct more suitable models of strongly correlated materials.

As an illustration, we have applied the method within LDA cRPA to calculate the effective interactions in two well-known cuprate parent compounds, LCO and HBCO, as well as in the prototype of correlated metals, SrVO3. We first studied the ${\bf r}$-dependence of $U({\bf r},{\bf r}^{\prime};\omega=0)$, both with ${\bf r}^{\prime}$ put at a transition metal nucleus (Cu or V) and at an in-plane O nucleus. In the $t_{2g}$ model of SrVO3, with ${\bf r}^{\prime}$ at the V nucleus, only a small region with weak attraction was found, which did not match the shape of the $xy$ low-energy orbital of the model. In the one-band model of the cuprates, on the other hand, a strong attractive interaction was found at the exact region of the low-energy 3$d_{x^{2}-y^{2}}$ orbital. Although this does not imply that charge fluctuations mediate Cooper pairing in the cuprates, they may assist other agents such as phonons and spin fluctuations in inducing pairing.

The temporal interaction exhibited generic damped oscillations in all compounds. Its time integral was shown to be the potential caused by inserting an impurity at $\tau=0$, and the results for the three-band model illustrated the possibility of finite-time overscreening, with an attractive effective interaction, despite the static limit being repulsive.

Acknowledgements.

This work was supported by the Swedish Research Council.

Appendix A Computational Details

We use the DFT code FLEURspexHP which utilizes the full-potential linearized augmented plane-wave (FLAPW) method to obtain all eigenfunctions $\phi_{n{\bf k}}$ and eigenvalues $\varepsilon_{n{\bf k}}$. All calculations are performed using the LDA. The band structures of HBCS, LCO and SrVO3 are provided in Fig. 13 together with their crystal structures.latwocuofour ; hg ; srvothree In LCO and HBCO we study $U$ in the well-established 1- and 3-band models, the former with a single Wannier function at Cu with $d_{x^{2}-y^{2}}$ symmetry and the latter also with two additional Wannier functions at the in-plane O atoms with $p_{x}$ and $p_{y}$ symmetry respectively. For comparison we also study $U$ in SrVO3 in the $t_{2g}$ model, with three Wannier functions at V with $d_{xy}$, $d_{xz}$ and $d_{yz}$ symmetry. The Wannier interpolated band structures are provided together with the LDA band structures in Fig. 13.

For the calculation of the RPA response matrix elements in the mixed product basis, $\chi_{\alpha\beta}^{\text{RPA}}({\bf k};\omega)$, we employ the SPEX code,spexHP which uses the ab initio LDA eigensolution as the unperturbed mean-field reference system. The response matrix is then utilized to compute $W$ and $U$ in position representation in the way we have described in the present paper. Since the full frequency dependence is required for the calculation of the real-time dynamics, we have taken care to include all relevant screening processes, also virtual transitions from low-lying semicore states: Cu 3$p$ and V 3$p$. These states play an important role for large values of $\omega$, and it is indeed the 3$p$ local orbitals which are responsible for the large peak structures at around 100 eV in Fig. 11. In time domain, this only affects the first main interaction minimum. The interesting time interval around 1-2 fs is essentially unaffected.

Surprisingly, the calculation turned out to be well converged with a sparse $4\times 4\times 4$ ${\bf k}$-mesh. The effect of increasing the mesh-size to $8\times 8\times 8$ was minimal. All calculations are therefore performed using a ${\bf k}$ mesh of size 4$\times$4$\times$4. The CuO2 sheets are, for simplicity, assumed to be perfectly two-dimensional without any buckling.