Optical properties of periodic systems within the current-current response framework: pitfalls and remedies

TL;DR

This paper compares density-density and current-current response functions for optical absorption in extended systems, highlighting pitfalls and proposing remedies to ensure physical accuracy in calculations.

Contribution

It identifies issues in current-based response calculations at zero momentum and offers corrected methods for including lifetime effects and electron-hole interactions.

Findings

Incorrect velocity operator definitions cause unphysical spectra.

Proper inclusion of lifetime effects improves spectral accuracy.

Methodology applies to semiconductors, metals, and insulators.

Abstract

We compare the optical absorption of extended systems using the density-density and current-current linear response functions calculated within many-body perturbation theory. The two approaches are formally equivalent for a finite momentum $\mathbf{q}$ of the external perturbation. At $\mathbf{q}=\mathbf{0}$, however, the equivalence is maintained only if a small $q$ expansion of the density-density response function is used. Moreover, in practical calculations this equivalence can be lost if one naively extends the strategies usually employed in the density-based approach to the current-based approach. Specifically we discuss the use of a smearing parameter or of the quasiparticle lifetimes to describe the finite width of the spectral peaks and the inclusion of electron-hole interaction. In those instances we show that the incorrect definition of the velocity operator and the violation of the conductivity sum rule introduce unphysical features in the optical absorption spectra of three paradigmatic systems: silicon (semiconductor), copper (metal) and lithium fluoride (insulator). We then demonstrate how to correctly introduce lifetime effects and electron-hole interactions within the current-based approach.

Equations

$$\hat{h}(t)=\left\{\frac{1}{2}[-i{\nabla }_{\mathbf{x}}+\frac{1}{c}\mathbf{A}(\mathbf{x},t){]}^2+v(\mathbf{x},t)\right\}$$

$$v(\mathbf{x},t)$$

$$\mathbf{A}(\mathbf{x},t)$$

$$\rho (\mathbf{x},t)$$

$$\mathbf{j}(\mathbf{x},t)$$

$$\hat{\mathbf{v}}(t)$$

$${\chi }_{\rho \rho }(\mathbf{q},\omega )$$

$${\chi }_{jj}(\mathbf{q},\omega )$$

$$\epsilon [{\chi }_{jj}^{\mathrm{IP}}](\mathbf{q},\omega )=1-\frac{4\pi }{{\omega }^2}{\chi }_{jj}^{\mathrm{IP}}(\mathbf{q},\omega ).$$

$$\epsilon [{\chi }_{\rho \rho }^{\mathrm{IP}}](\mathbf{q},\omega )=1-\frac{4\pi }{q^2}{\chi }_{\rho \rho }^{\mathrm{IP}}(\mathbf{q},\omega ),$$

$$q^2{\chi }_{jj}(\mathbf{q},\omega )={\omega }^2{\chi }_{\rho \rho }(\mathbf{q},\omega ),$$

$${\chi }_{aa}^{\mathrm{IP}}(\mathbf{q},\omega )=-\frac{1}{V}\sum _{nm}\frac{a_{\mathbf{q},(nm)}^{\mathrm{IP}}{\tilde{a}}_{\mathbf{q},(mn)}^{\mathrm{IP}}}{{\omega }^{+}-\Delta {ϵ}_{mn}^{\mathrm{IP}}(\mathbf{q})}$$

$$\Delta {ϵ}_{mn}^{\mathrm{IP}}(\mathbf{q})$$

$$\Delta f_{mn}(\mathbf{q})$$

$$a_{\mathbf{q},(nm)}^{\mathrm{IP}}$$

$${\tilde{a}}_{\mathbf{q},(mn)}^{\mathrm{IP}}$$

$${\omega }^{+}=\underset{\eta \to 0}{\lim }(\omega +i\eta ),$$

$${\chi }_{jj}^{\mathrm{IP}}(\mathbf{q},\omega )={\chi }_{j^p{j}^p}^{\mathrm{IP}}(\mathbf{q},\omega )+\frac{N}{V}.$$

$$v_0(\mathbf{r})=v_N(\mathbf{r})+v_H(\mathbf{r})+v_{xc}(\mathbf{r}),$$

$$v_{0,\sigma {\sigma }^{\prime }}(\mathbf{x},{\mathbf{x}}^{\prime },\omega )=v_N(\mathbf{r})+v_H(\mathbf{r})+{\Sigma }_{xc,\sigma {\sigma }^{\prime }}(\mathbf{x},{\mathbf{x}}^{\prime },\omega ),$$

$$H^{{\rm QP}}_{nm}(\omega)=\delta_{n,m}\Big{[}\epsilon^{{\rm KS}}_{n}+\langle n|\Sigma_{xc}(\omega)-v_{xc}|n\rangle\Big{]}.$$

$${ϵ}_n^{\mathrm{QP}}={ϵ}_n^{\mathrm{KS}}+⟨n|{\Sigma }_{xc}\left({ϵ}_n^{\mathrm{QP}}\right)-v_{xc}|n⟩,$$

$${\chi }_{aa}^{\mathrm{IP}}(\mathbf{q},\omega )=-\sum _{nm}\frac{K_{nm,\mathbf{q}}^{aa}}{{\omega }^{+}-\Delta {ϵ}_{mn}^{\mathrm{IP}}(\mathbf{q})},$$

$$F(\omega )=\frac{1}{\omega -\Delta ϵ}=-\frac{1}{\Delta ϵ}-\frac{\omega }{\Delta {ϵ}^2}+\frac{{\omega }^2}{(\omega -\Delta ϵ)\Delta {ϵ}^2},$$

$$\epsilon [{\chi }_{jj}^{\mathrm{IP}}](\mathbf{q},\omega )=1+\frac{A^{\mathrm{IP}}(\mathbf{q})}{{\omega }^2}+\frac{B^{\mathrm{IP}}(\mathbf{q})}{\omega }+C^{\mathrm{IP}}(\mathbf{q},\omega ),$$

$$A^{\mathrm{IP}}(\mathbf{q})$$

$$B^{\mathrm{IP}}(\mathbf{q})$$

$$C^{\mathrm{IP}}(\mathbf{q},\omega )$$

$${\chi }_{jj}^{\mathrm{IP}}(\mathbf{q},0)={\chi }_{j^p{j}^p}^{\mathrm{IP}}(\mathbf{q},0)+\frac{N}{V}=0,$$

$$\epsilon [{\chi }_{jj}^{\mathrm{IP}}](\mathbf{q},\omega )=1+C^{\mathrm{IP}}(\mathbf{q},\omega )=\epsilon [{\chi }_{\rho \rho }^{\mathrm{IP}}](\mathbf{q},\omega ).$$

$${\rho }_{\mathbf{q},(nm)}^{\mathrm{IP}}=-iq{d}_{\mathbf{q},(nm)}^{\mathrm{IP}}+O(q^2),$$

$${\chi }_{dd}^{\mathrm{IP}}(\mathbf{q},\omega )=-$$

$$\epsilon [{\chi }_{dd}^{\mathrm{IP}}](\mathbf{q},\omega )=1-4\pi \sum _{nm}\frac{K_{nm,\mathbf{q}}^{dd}}{{\omega }^{+}-\Delta {ϵ}_{mn}^{\mathrm{IP}}(\mathbf{q})},$$

$$\underset{\mathbf{q}\to 0}{\lim }\epsilon [{\chi }_{jj}^{\mathrm{IP}}](\mathbf{q},\omega )$$

Full text

I Introduction

The electromagnetic linear response of solids can be measured experimentally by applying a small external perturbation-or probe-which induces a small change in the sample material. This change is the response of the material to the perturbing field and it can have both longitudinal and transverse components, depending on the experimental setup and the inhomogeneity of the system. Well-known examples of such experiments are the measurements of absorption, electron energy loss, Kerr and Faraday effects, and dichroism.

From the theoretical point of view the linear response of a system to longitudinal fields can be obtained from both the density-density response function $\chi_{\rho\rho}$ and the current-current response function $\underline{\underline{\chi}}_{jj}$, which is a tensor Strinati (1988). Instead, the response to transverse fields can only be described using $\underline{\underline{\chi}}_{jj}$. This is due to the fact that the density determines the longitudinal current—through the continuity equation— but not the transverse current. An important instance where transverse electric fields come into play is the recent experimental progress on topological insulators Checkelsky et al. (2012); Chang et al. (2013). Therefore approaches based on $\underline{\underline{\chi}}_{jj}$ are more general than those based on $\chi_{\rho\rho}$. On the other hand the current–density based approach is more susceptible to numerical issues and instabilities. For example the paramagnetic and diamagnetic contributions to $\underline{\underline{\chi}}_{jj}$ must be treated on equal footing Kootstra et al. (2000); Raimbault et al. (2015) otherwise divergencies can arise.

The differences between the two response functions are intrinsically related to the different gauge used to define the coupling of the external field with the electrons. While transverse fields can be described only in terms of the transverse component of the vector potential, for longitudinal fields there are two options. One can use a scalar potential, which couples to the electron density of the system, or the longitudinal component of the vector potential, which couples to the current density. In the optical limit, i.e., in the limit where the momentum $\mathbf{q}$ carried by the external perturbation is negligible, these two gauges are called length and velocity gauge, respectively. In this limit the distinction between longitudinal, or parallel to ${\mathbf{q}}$, and transverse, or perpendicular to ${\mathbf{q}}$, fields vanishes. However the optical limit is also the case where many subtle differences between the two approaches arise, fundamentally because the density-based approach is ill-defined at $\mathbf{q}=0$.

Various works in the literature have discussed and compared the two approaches either within the independent-particle (IP) picture Sipe and Ghahramani (1993); Lamb et al. (1987); Sangalli et al. (2012), in which optical properties can be described in terms of a simple sum-over-states expression, or within a generalized Kohn-Sham approach. Springborg and Kirtman (2008) These studies focus mostly on the limitations of the density-based approaches and the avoidance of divergencies in the current-based approaches. An important point of debate is whether the two approaches give the same longitudinal macroscopic dielectric function. At the linear-response level the equivalence has been shown, but only in some specific cases, i.e., for semiconductors and insulators, for Hamiltonians with only local operators, and for absorption at resonance Schafer and Wegener (2002). The comparison for metals, for Hamiltonians with nonlocal operators, and for absorption out of resonance due to smearing,  111 The differences out of resonance can be of small importance in the linear regime but they become crucial for any non-linear phenomena because they play an important role in the construction of high order response functions. has received less or no attention. More importantly, no systematic comparison of the density-based and the current-based approaches has been carried out so far at the many-body level.

The objectives of this work are (i) to better elucidate the differences between the density and the current-based approaches at the IP level; (ii) to extend the discussion to the many-body framework. In particular, at the IP level, we compare the two approaches when lifetime effects are introduced (either from first-principles or by introducing a smearing parameter). At the many-body level we consider the Bethe-Salpeter equation (BSE). Based on Green functions theory, the BSE is the fundamental equation of Many-Body Perturbation Theory for the description of linear response properties Fetter and Walecka (2003). In particular we address some important and rather subtle aspects related to the definition of the velocity operator in the current-based BSE approach Rohlfing and Louie (2000) which did not receive proper attention so far. Strinati (1988); Rohlfing and Louie (2000); Del Sole and Fiorino (1984)

We note that alternative approaches exist in which the explicit calculation of response functions is avoided, for example by calculating the induced current-density in the frequency domain Kootstra et al. (2000); Romaniello and de Boeij (2005); Berger et al. (2005, 2006) or by calculating the time-dependent density or current-density using real-time propagation Attaccalite et al. (2011); Castro and Appel (2006); Bertsch et al. (2000). Many of the issues discussed in this work are also pertinent to these methods.

The manuscript is organized as follows. In Sec. II we compare the density and the current-based approaches at a formal level, within both the independent-particle approximation and the BSE. We compare the two approaches in the optical limit $\mathbf{q}\rightarrow{\bf 0}$ and elucidate the origin of their fundamental difference at $\mathbf{q}={\bf 0}$. In Sec III we show how an incorrect implementation of smearing and an erroneous definition of the velocity operator can lead to different spectra in the two approaches. We then show how these problems can be solved. For sake of completeness, in Sec. IV we briefly discuss alternative approaches to calculate optical spectra. Finally, in Sec. V we draw our conclusions.

II Formal equivalence of density and current-based approaches

In this section we review the basic equations which describe optical properties within both a density-based and a current-based approach and we compare these approaches in the optical limit $\mathbf{q}\rightarrow{\bf 0}$. We first consider the case of noninteracting electron-hole pairs referred to in the following as IP level, 222In condensed matter the IP approximations is often called RPA to be distinguished from RPA plus local field effects. and then the case in which the electron-hole interaction is treated within the BSE of Many-Body Perturbation Theory. In both cases we demonstrate that the two approaches should lead to the same results for the longitudinal macroscopic dielectric function (see App. B for its definition).

II.1 Noninteracting electron-hole pairs

Let us consider a collection of electrons moving in a periodic potential ${v_{0}({\mathbf{x}})=v_{0}({\mathbf{r}})}$, with ${\mathbf{x}}={\mathbf{r}}+{\mathbf{R}}$ and ${\mathbf{r}}$ inside the unit cell volume $V$, whose position is identified by the Bravais lattice vector ${\mathbf{R}}$. We are interested in the response of the system to a macroscopic time-dependent electromagnetic field characterized by a set of scalar and vector potentials, $\{\delta\phi({\mathbf{x}},t),\delta{\mathbf{A}}({\mathbf{x}},t)\}$. A macroscopic quantity is defined as an average taken over the unit cell whose location is given by ${\mathbf{R}}$ Wiser (1963) (see App. A).

The motion of the system is then governed by the following one-particle Hamiltonian,

[TABLE]

where ${\mathbf{x}}={\mathbf{r}}+{\mathbf{R}}$, and

[TABLE]

for which the macroscopic external perturbations $\delta\phi({\mathbf{R}},t)$ and $\delta{\mathbf{A}}({\mathbf{R}},t)$ vanish identically for $t\leq 0$. The $v_{0}({\mathbf{r}})$ term instead describes a periodic potential (i.e. it does not depend on ${\mathbf{R}}$) and can be seen as the mean field felt by the electrons. Here and throughout the article we use atomic units $(\hbar=m_{e}=e=4\pi\varepsilon_{0}=1)$ with the gaussian convention for electromagnetism. It is important to notice that the perturbing potentials $\delta\phi$ and $\delta{\mathbf{A}}$ are external potentials and they do not take into account the contributions due to the response of the system.

The density and the current-density of the system are defined in terms of the Bloch wave functions $\psi_{n}({\mathbf{x}},t)$, which are eigenstates of the one-particle hamiltonian $\hat{h}(t)$, as

[TABLE]

Here $n=\{\tilde{n}{\bf k}\}$ is a generalized index that comprises the band index $\tilde{n}$ and the wave vector ${\bf k}$, and $f_{n}$ are the occupations factors. In the collinear case $\tilde{n}=\{i\sigma\}$ can be further specified as the band with the corresponding spin index $\sigma$, while in the spinorial case it is the spinorial band index. The density operator is the identity in real space $\hat{\rho}=\hat{1}$. The current-density is given in terms of the velocity operator Landau and Lifshits (1965).

[TABLE]

where $\hat{h}_{\pm}=(\hat{h}\pm\hat{h}^{\dagger})/2$, $[\hat{a},\hat{b}]$ is the commutator, while $\{\hat{a},\hat{b}\}$ is the anticommutator. With the Hermitian Hamiltonian defined in Eq. (1) ${\hat{\mathbf{v}}(t)=(-i\nabla_{\mathbf{x}}+\mathbf{A}({\mathbf{x}},t)/c)}$. The momentum operator $-i\nabla_{\mathbf{x}}$ gives the paramagnetic current density, i.e. ${\hat{j}^{p}=-i\nabla_{\mathbf{x}}}$, and the potential-dependent term gives the diamagnetic current density. We note that the velocity operator, and therefore the current operator, depends on the Hamiltonian. We will discuss this important point later in the paper. Density and current-density variations are induced as a response to the perturbing potentials in Eqs. (2a) and (2b).

We restrict ourselves to macroscopic longitudinal perturbations with a small transferred momentum ${\mathbf{q}}$. In the following, therefore, we consider the longitudinali, i.e., parallel to ${\mathbf{q}}$, component only of the vector potential and the current (see App. B). Longitudinal perturbations can be described by a scalar potential $\phi$, for example in the Coulomb gauge \mbox{\boldmath\nabla}\cdot{\mathbf{A}}=0, or by a longitudinal vector potential ${\mathbf{A}}$, for example in the (incomplete) Weyl gauge $\delta\phi=0$. In the optical limit the two gauges generate the so-called length and velocity gauges (see App. C).

The macroscopic density-density,and the longitudinal current-current response functions, $\chi_{\rho\rho}$ and $\chi_{jj}$, respectively, can be written in reciprocal space as (see App. D)

[TABLE]

When the induced density and current density are calculated from the Hamiltonian in (1), then the response functions in Eqs. (5) are the independent particle ones, $\chi^{{\rm IP}}_{\rho\rho}$ and $\chi^{{\rm IP}}_{jj}$. We can now define the longitudinal component of the macroscopic dielectric tensor $\varepsilon({\bf q};\omega)$, which conveniently describes the optical properties of semiconductors in the long wavelength limit. It can be obtained, for $\omega\neq 0$, in terms of the longitudinal current-current response function as Strinati (1988)

[TABLE]

For ${\bf q}\neq{\bf 0}$, the dielectric function can also be expressed in terms of the density-density response function as

[TABLE]

with $q=|{\bf q}|$. Equation (7), in the limit ${\bf q}\rightarrow{\bf 0}$, is the expression commonly used for the calculations of optical properties in solids. For both expressions (6) and (7) the non-analytic point $({\bf q},\omega)=({\bf 0},0)$ can be described only via a limiting procedure and the final result depends on the direction of the limit.

Both $\varepsilon[\chi_{jj}]$ and $\varepsilon[\chi_{\rho\rho}]$ should lead to the same result. Indeed the two response functions are related by the expression Giuliani and Vignale (2005); Schafer and Wegener (2002)

[TABLE]

which follows from the continuity equation ${{\bf q}\cdot\mathbf{j}=\omega\rho}$, which guarantees local charge conservation. Thus at ${{\bf q}\neq{\bf 0}}$ and ${\omega\neq 0}$ relation (8) ensures $\varepsilon[\chi_{jj}]=\varepsilon[\chi_{\rho\rho}]$. We note that although we have established formal equivalence this does not guarantee numerical equivalence, as we show in Section III.

II.1.1 Choice of the reference one-particle hamiltonian

Within the IP picture defined by the Hamiltonian in Eq. (1) the response functions can be written as

[TABLE]

with

[TABLE]

(the spin index in $n$ and $m$ is identical in the collinear case) and

[TABLE]

where it is understood that the limit $\eta\rightarrow 0$ is taken at the end of the calculation. Here $\hat{a}$ is either the one-particle density operator $\hat{\rho}$ or the paramagnetic current-density operator $\hat{j}^{p}$, and $|n\rangle$ are Bloch states, which are eigenstates of the equilibrium one-particle Hamiltonian $\hat{h}(t=0$), with corresponding energies $\epsilon^{{\rm IP}}_{n}$ and occupation numbers ${0\leq f_{n}\leq 1}$.333The inclusion of $\sqrt{\Delta f_{mn}({\bf q})}$ in the definition of the matrix elements of Eq. (13) allows to define in the following sections (where the the case of interacting electron-hole pairs is considered) an excitonic matrix which remains Hermitian also in the general case of fractional occupation numbers. See also Ref. Sangalli et al., 2016. It is important to notice that the transitions with ${\Delta\epsilon^{{\rm IP}}_{mn}({\mathbf{q}})=0}$ give no contribution since at equilibrium $\Delta{f_{mn}({\mathbf{q}})=0}$. Thus all the summations in the present manuscript are intended without the zero energy poles. The full current-current response function is obtained via

[TABLE]

We shall now choose the stationary part of the one-particle Hamiltonian (1). A common choice is the Kohn-Sham Hamiltonian from density functional theory (DFT) where

[TABLE]

is the equilibrium Kohn-Sham (KS) potential. It is the sum of the potential generated by the nuclei, $v_{N}({\mathbf{r}})$, the Hartree potential, $v_{H}({\mathbf{r}})$, and the exchange-correlation potential, $v_{xc}({\mathbf{r}})$. We refer to the IP response function derived from the KS hamiltonian as $\chi^{{\rm KS}}_{aa}$. However, usually DFT-KS band structures are not a good starting point for response calculations, since, for example, the fundamental band gap is systematically underestimated in semiconductors and insulators. To improve over DFT one can introduce a more general nonlocal and frequency dependent quasiparticle (QP) potential

[TABLE]

where $\Sigma_{xc}$ is the exchange-correlation many-body self-energy, defined within the formalism of Many-Body Perturbation Theory (MBPT), evaluated at equilibrium. In particular we consider the (first-order) QP hamiltonian $H^{{\rm QP}}$ defined as the first-order correction in the perturbation ${\Sigma_{xc}-v_{xc}}$ of the KS hamiltonian:

[TABLE]

We have here emphasized the $\omega$ dependence of $\Sigma_{xc}$, and dropped the ${\mathbf{r}}$ space and $\sigma$ spin dependence of both $\Sigma_{xc}$ and $v_{xc}$. The QP eigenvalues are thus corrected KS energies

[TABLE]

while we keep the same wave functions ${\psi^{{\rm QP}}_{n}=\psi^{{\rm KS}}_{n}}$. From the quasiparticle energies $\epsilon^{{\rm QP}}$, solution of Eq. (19), we can define the QP response function, $\chi^{{\rm QP}}_{aa}$, which differs from the KS response function by the replacement $\epsilon^{{\rm KS}}_{n}\rightarrow\epsilon^{{\rm QP}}_{n}$.

II.1.2 Expansion of the current-based approach at finite momentum

We now focus on the relation between the density-based and current-based formalism by extending the approach of Ref. Sipe and Ghahramani, 1993 to ${\bf q}\neq\mathbf{0}$. To simplify the notation we rewrite the response function (9) as

[TABLE]

where $K^{aa}_{nm,{\bf q}}=a^{{\rm IP}}_{{\bf q},(nm)}\,\tilde{a}^{{\rm IP}}_{{\bf q},(mn)}/V$. By using the exact relation 444which can be obtained by Taylor expanding ${F(\omega)}$ around ${\omega=0}$ up to the first order. The last term on the right-hand side is then defined from the difference ${F(\omega)-F(0)-F^{\prime}(0)\,\omega}$.

[TABLE]

in Eq.(20) and inserting the result into Eq.(6), we can decompose $\varepsilon[\chi^{{\rm IP}}_{jj}]({\bf q},\omega)$ as

[TABLE]

with

[TABLE]

where in Eq.(23a) we used Eqs.(15) and (20). From the conductivity sum rule (CSR) for $\chi^{{\rm IP}}_{jj}$, given by

[TABLE]

we see that $A^{{\rm IP}}({\bf q})=0$. Moreover, if time-reversal symmetry holds, also $B^{{\rm IP}}({\bf q})=0$. Therefore, in case of time-reversal symmetry, only the $C({\bf q},\omega)$ term survives and we have, at finite momentum, the general result

[TABLE]

The last equality 555Eq. (25) implies $q^{2}K^{j^{p}j^{p}}_{nm,{\bf q}}=\Delta\epsilon_{nm}^{2}({\bf q})K^{\rho\rho}_{nm,{\bf q}}$, which is similar, but not the same, to Eq. (8). Notice the non trivial replacement of $\omega^{2}\rightarrow\Delta\epsilon_{nm}^{2}({\bf q})$. holds thanks to Eq. (8). In the next section we will use (23a)-(23c) to discuss the optical limit.

II.1.3 The optical limit

We are interested in computing optical properties, for which $\omega=c\,q$, with $\omega$ in the order of a few eV. Therefore we consider $q\approx 0$ compared to the size of the Brillouin zone. However, $\varepsilon[\chi^{{\rm IP}}_{\rho\rho}]$, given in Eq. (7), is not defined for $q=0$ because of the $1/q^{2}$ term. Therefore, to obtain an explicit expression for $\varepsilon[\chi^{{\rm IP}}_{\rho\rho}]$ for small $q$, we Taylor expand $\rho^{{\rm IP}}_{{\bf q},(nm)}$, defined in Eq. (13), around $q=0$,

[TABLE]

where 666The position operator is ill-defined when periodic boundary conditions (PBC) are imposed. Here we implicitly work in the crystal momentum representation in which the matrix elements of the position operator are redefined consistently with the PBC Blount (1962); Souza et al. (2004) ${d^{{\rm IP}}_{{\bf q},(nm)}=\langle n|\hat{\bf q}\cdot{\mathbf{x}}|m\rangle\sqrt{\Delta f_{nm}({\bf q})}}$ with $\hat{\bf q}={\bf q}/q$ the direction of ${\bf q}$. Substitution into Eq. (20) leads to $\chi^{{\rm IP}}_{\rho\rho}({\mathbf{q}},\omega)=q^{2}\chi^{{\rm IP}}_{dd}({\mathbf{q}},\omega)+O(q^{3})$ with

[TABLE]

the longitudinal dipole-dipole response function. Thus, for small $q$, we can rewrite Eq. (7) as

[TABLE]

which is well defined at ${\bf q}={\bf 0}$.777Here $d_{{\mathbf{q}}}$ at ${\bf q}={\bf 0}$ is the longitudinal dipole, defined by the direction of the ${\bf q}\rightarrow{\bf 0}$. See also App. B. It is Eq. (28) that is usually implemented to compute absorption spectra of cold semiconductors and insulators in the density-based approach.

However, Eq. (28) is formally exact only in the limit ${\bf q}\rightarrow{\bf 0}$, since at ${\bf q}={\bf 0}$ contributions from intraband transitions, i.e., transitions within a single band, are excluded 888For finite ${\bf q}$ Eq. (28) is exact to first order in ${\bf q}$. Instead, no such problem exists for the current-based approach. We can summarize this difference between the two approaches as

[TABLE]

We now discuss the cases ${\bf q}={\bf 0}$ and ${\bf q}\rightarrow{\bf 0}$ in more detail by separating $A^{{\rm IP}}({\bf q})$, $B^{{\rm IP}}({\bf q})$ and $C^{{\rm IP}}({\bf q})$ into inter- and intraband contributions. Moreover, $A^{{\rm IP}}({\bf q})$ also contains the constant diamagnetic term $A^{d}=-4\pi N/V$.

The ${\bf q}={\bf 0}$ case: interband transitions.

At ${{\bf q}={\bf 0}}$ only interband transitions, i.e., transitions between two different bands, contribute to the summations in Eqs. (23a)-(23c), since for intraband transitions $\Delta f_{nn}=0$. Because of the missing intraband contributions in $\chi^{{\rm IP}}_{jj}$ the CSR in Eq. (24) in general does not hold at ${\bf q}={\bf 0}$. Let us define

[TABLE]

One can verify 999Here $B^{\textrm{IP,inter}}=0$ because we are considering the longitudinal term only. Indeed the mixed longitudinal-transverse terms can instead be different from zero and describe the Anomalous Hall effect Sangalli et al. (2012). that $B^{{\rm IP}}({\bf 0})=0$ by exchanging $n$ and $m$ in Eq. (23b).

Let us first consider systems with a gap. Then $A^{\textrm{IP,inter}}+A^{d}=0$ since the dielectric function must go to a constant Sipe and Ghahramani (1993) as $\omega\rightarrow 0$. We can thus focus on $C^{\textrm{IP,inter}}(\omega)$. Using Eq. (12) for the paramagnetic current-density, we can write

[TABLE]

where we used $\hat{\mathbf{v}}=-i[\hat{h},\hat{\mathbf{x}}]$ and $d^{{\rm IP}}_{{\bf 0},(nm)}$ given by Eq. (12) with $\hat{a}=\hat{\mathbf{x}}$. From this relation we deduce that

[TABLE]

Substitution of this identity into Eq. (28) shows that at ${\bf q}={\bf 0}$ both $\varepsilon[\chi^{{\rm IP}}_{jj}]$ and $\varepsilon[\chi^{{\rm IP}}_{dd}]$ can be expressed in terms of $C^{\textrm{IP,inter}}$ for systems with a gap:

[TABLE]

Equation (33) proves the equivalence between the density-based and the current-based approaches for cold semiconductors and insulators.

For metals, however, we also need to describe the divergent Drude-like term. Since the current-based approach is exact also at ${\bf q}={\bf 0}$, this term must be described by $(A^{\textrm{IP,inter}}+A_{d})/\omega^{2}$, i.e. the Drude tail originates from the breaking of the CSR. Thus we have

[TABLE]

Instead, the density-based approach does not contain any extra term beyond $C^{\textrm{IP,inter}}(\omega)$ and thus cannot describe metals at ${\bf q}={\bf 0}$. The Drude-like tails in the density-based approach can be obtained only explicitly dealing with the ${\bf q}\rightarrow{\bf 0}$ limit as explained in the next subsection.

The ${\bf q}\rightarrow{\bf 0}$ limit: intraband transitions.

In the ${\bf q}\rightarrow{\bf 0}$ limit also intraband transitions contribute to the summations in Eqs. (23a)-(23c). One can thus define

[TABLE]

where we used Eq. (30b) and the fact that $A^{{\rm IP}}({\bf q}\neq{\bf 0})={\bf 0}$ owing to the CSR in Eq. (24). If time-reversal symmetry holds $B^{{\rm IP}}({\bf q})=0$ and hence $B^{\textrm{IP,intra}}=0$. We note that here and in the rest of the paper the ${\bf q}\rightarrow 0$ limit is taken at finite $\omega$. Using the results of Eqs. (35a)-(35c) in Eq. (22) we obtain

[TABLE]

where we used Eq. (30c). Owing to Eq. (29a) we can compare this result to Eq. (34) and conclude that

[TABLE]

This means that, within the current-based approach, the Drude-like tail, which at ${\bf q}={\bf 0}$ is described by $(A^{\textrm{IP,inter}}+A^{d})/\omega^{2}$, in the limit ${\bf q}\rightarrow{\bf 0}$ is described via $C^{\textrm{IP,intra}}(\omega)$. Indeed, the exact expression for $\varepsilon[\chi^{{\rm IP}}_{jj}]$ at finite ${\bf q}$ in Eq. (25) only depends on $C^{{\rm IP}}({\bf q},\omega)$. Thus $C^{\textrm{IP,intra}}(\omega)$ must describe all intraband transitions in metals when ${\bf q}\rightarrow{\bf 0}$.Marini (2001).

To have an explicit expression of such intraband contribution in the density-based approach, one needs to Taylor expand the energy and occupation number differences as Romaniello and de Boeij (2005)

[TABLE]

where $\varv_{n{\bf k}}=\nabla_{\bf k}\epsilon_{n{\bf k}}$. At the IP level there is no Drude tail in the absorption (see App.  E for details). Only when introducing a smearing a Drude-like peak in the absorption appears Romaniello and de Boeij (2005); Berger et al. (2005). Experimentally such smearing also exist, because of the interaction between electrons and thus in practice a peak is always measured.

In conclusion, the current-based approach is exact both in the limit ${\bf q}\rightarrow{\bf 0}$ and at ${\bf q}={\bf 0}$. In particular, if metals are considered, the Drude-like tail is present in both cases, but it is described by different terms, namely by $C^{\textrm{IP,intra}}(\omega)$ in the limit ${\bf q}\rightarrow{\bf 0}$ and by $(A^{\textrm{IP,inter}}+A^{d})/\omega^{2}$ at ${\bf q}={\bf 0}$. Also the density-based approach is exact in the limit ${\bf q}\rightarrow{\bf 0}$. This can be summarized by the following set of relations

[TABLE]

II.2 Interacting electron-hole pairs

So far we have only considered independent particles. We now also take into account the electron–hole interaction by considering the variations induced in the Hartree and the exchange–correlation self–energy by the external potential. Within MBPT such variations can be conveniently described using the BSEStrinati (1988) for the (time-ordered) two-particle propagator $L$, which reads

[TABLE]

where $i\equiv{{\mathbf{r}}_{i},t_{i}}$. In Eq. (40) $L_{0}$ is given in terms of QP Green’s functions as ${L_{0}(1,2,1^{\prime},2^{\prime})=-iG^{{\rm QP}}(1,2^{\prime})G^{{\rm QP}}(2,1^{\prime})}$ and the four-point kernel $\bar{\Xi}(3,5,4,6)$ is given by

[TABLE]

where $\bar{v}_{H}$ is the Hartree potential without the long-range component of the Coulomb potential, $v_{c}(\mathbf{G}=0)$. We are now looking at the self-consistent response of the system to a macroscopic field composed of both the external and the induced macroscopic field. Therefore $v_{c}(\mathbf{G}=0)$ is effectively removed Bussi (2004). We use the $GW$ self-energy with an instantaneous screened Coulomb potential $W$. In this case the propagator $L$ depends only on the time difference $\tau=t_{1}-t_{2}$101010As commonly done in the literature, we neglect the term $i\delta W/\delta G$ in the kernel.. A Fourier transformation to frequency space leads to

[TABLE]

where, for notational convenience, we dropped the spin and space arguments 111111for a detailed treatment of the space indexes at finite momentum see Ref. [Gatti and Sottile, 2013]; for the treatment of occupations factors see Ref. [Sangalli et al., 2016]).

The BSE can be mapped onto an effective two-particle equation, written in a basis of electron-hole transitions $(nm)$, according to  Onida et al. (2002); Gatti and Sottile (2013)

[TABLE]

where $H^{Exc}({\bf q})$ is the excitonic Hamiltonian with eigenvectors $A_{\lambda,{\bf q}}$ and eigenvalues $E_{\lambda}({\bf q})$. The (retarded)121212The exact Dyson equations holds only for the time-ordered two-particles propagator, which is formally derived assuming zero temperature (i.e. integer occupation numbers also in case of metals). To consider fractional occupations one would need to introduce a finite temperature formalism. However, using a static kernel, the Dyson equation at finite temperature reduces to a Dyson equation identical to Eq. (42) but for the retarded propagator, which is indeed what is needed to define the dielectric function. Thus from now on we can consider all quantities as retarded functions and forget about the time-ordered formalism. particle-hole propagator $\bar{L}$ can be obtained via

[TABLE]

where $B_{\lambda,{\bf q}}$ and the overlap matrix $S_{\lambda\lambda^{\prime},{\bf q}}$ are defined by

[TABLE]

Upon solving Eq. (43), one can obtain the density-density and the current-current excitonic response functions $\bar{\chi}^{Exc}_{aa}$ as

[TABLE]

where $a^{Exc}_{{\bf q},\lambda}$ and $\tilde{a}^{Exc}_{{\bf q},\lambda^{\prime}}$ are defined analogously to $a_{{\bf q},\lambda}$ and $\tilde{a}_{{\bf q},\lambda^{\prime}}$ in Eqs. (12) and (13) but the expectation value is with respect to the excitonic wave function. As for the IP case, the zero-energy transitions give no contribution to the summation 131313We assume that zero energy transitions at the BSE level originates from zero energy transitions at the IP level. Zero energy poles at the BSE level may originate from finite energies transition at the IP level as well. However this case would point to an instability of the ground state which would be degenerate to an excited state. We exclude this possibility in the present work. since $\Delta f_{nm}({\bf q})=0$ and thus ${B^{(nm)}_{\lambda,{\bf q}}=0}$. Note that neglecting the kernel $\bar{\Xi}$ in the BSE (42), Eq. (46) reduces to the IP response function given in Eq. (9).

Substitution of $\bar{\chi}^{Exc}_{jj}$ and $\bar{\chi}^{Exc}_{\rho\rho}$ in Eqs. (6) and (7), respectively, yields an expression for the macroscopic dielectric tensor in the current- and density-based approach. As in the IP case the two expressions are identical for ${{\bf q}\neq{\bf 0}}$ and ${\omega\neq 0}$ thanks to Eq. (8).

II.2.1 The optical limit

The analysis for the optical limit in Sec. II.1.2-Sec. II.1.3 also applies to the excitonic case. In particular, by rewriting the expression for the excitonic response function as

[TABLE]

where $K^{aa}_{\lambda\lambda^{\prime},{\bf q}}=a_{{\bf q},\lambda}\,S^{-1}_{\lambda\lambda^{\prime},{\bf q}}\,\tilde{a}_{{\bf q},\lambda^{\prime}}/V$ one obtains

[TABLE]

The interband dielectric function in terms of the excitonic dipole-dipole response function $\chi^{{\rm Exc}}_{dd}$ reads

[TABLE]

where now $K^{dd}_{\lambda\lambda^{\prime},{\bf q}}$ contains the longitudinal excitonic dipoles $d^{{\rm Exc}}_{{\bf q},\lambda}=\sum_{nm}A^{(nm)}_{\lambda,{\bf q}}x_{(nm)}\,\sqrt{\Delta f_{nm}({\bf q})}$.

Finally, the analogous of the relation (32) exists also for the exitonic case:

[TABLE]

from which expressions analogous to Eqs. (34) and (35) can be written for $\varepsilon[\chi^{{\rm Exc}}_{jj}]$.

III Issues in the practical implementation of current-based approaches

In the previous section the equivalence between the current-based and density-based approaches was established for the longitudinal macroscopic dielectric function. In the present section we compute optical absorption in a semiconductor (bulk silicon), an insulator (LiF), and a metal (copper) comparing the two approaches numerically, i.e., using both Eq. (6) for the current-based approach and Eq. (7) for the density-based approach. In doing so we will show potential pitfalls in the implementation of current-based approaches due to: (i) the use of a smearing parameter (Sec. III.2); (ii) the violation of the conductivity sum rule (Sec. III.3); (iii) the inclusion of lifetimes (Sec. III.4); (iv) the inclusion of QP energy corrections and excitonic effects with an incorrect velocity operator (Sec. III.5). Moreover, we demonstrate how those pitfalls can be avoided.

III.1 Computational details

The response functions entering Eqs. (6) and (7) are constructed starting from ground state DFT calculations performed with either the Quantum-Espresso codeGiannozzi et al. (2009) (Si and Cu) or the Abinit codeGonze et al. (2002) (LiF) using norm conserving pseudopotentials. The ground-state of all three materials have been calculated by using the local density approximation (LDA) for the exchange-correlation functional. For Si we used 4 valence electrons (${3s^{2}3p^{2}}$ configuration), a face-centered cubic (FCC) cell with a two atoms (diamond structure) and the experimental lattice parameter of $10.18$ Bohr. For calculating the ground state density we use an energy cutoff of 10 Ha and a $6\times 6\ \times 6$ sampling of the Brillouin zone (BZ). For the calculation of the response function we used instead a $16\times 16\times 16$ sampling of the BZ, resulting in 145 k-points in the irreducible BZ (IBZ) and 4096 in the BZ. For LiF we used 8 valence electrons (${2s^{1}}$ configuration for Li and ${2s^{2}2p^{5}}$ for F)Goedecker et al. (1996), a FCC cell with two atoms (sodium-chloride structure) and the experimental lattice parameter of $7.70$ Bohr. For calculating the ground state density and the response function we use an energy cutoff of 40 Ha with a $8\times 8\times 8$ sampling of the BZ, resulting in 29 (512) k-points in the IBZ (BZ). For Cu we used with 11 valence electrons (${3d^{10}4s^{1}}$ configuration), a FCC cell with a single atom and the experimental lattice parameter of $6.82$ Bohr. For calculating the ground state density and the response function we use 32.5 Ha and a $16\times 16\times 16$ sampling of the BZ, resulting in 145 (4096) k-points in the IBZ (BZ).

We computed the dielectric function starting from the DFT-KS wave functions and energies using the Yambo code Marini et al. (2009) where we implemented the equations for the current-based approach.141414Yambo standard implementation uses a density-based approach. For Si and LiF, we considered electron-hole pairs built from the top 3 valence and the bottom 3 conduction bands. For Cu we included 30 bands in the band summation. A scissor operator is used to mimic the effect of the $GW$ quasiparticle corrections in Si ($0.8$ eV) and LiF ($5.8$ eV) consistently to what already reported in the literature. For the BSE calculations the static screening in the random-phase approximation is computed using 50 bands in the band summation and 2.3 Ha energy cutoff in the summation over reciprocal lattice vectors for Si, and 30 bands and a 3.6 Ha energy cutoff for LiF.

We also consider QP lifetimes by introducing an imaginary part in the QP energies. To mimic the effect of the electron-phonon Fan self-energy Giustino (2016); Bernardi et al. (2014), we use a term proportional to the density of states. To mimic the effect of the $GW$ self-energy Aryasetiawan and Gunnarsson (1998) we use a term which grows quadratically in energy, from the Fermi level in Cu and from the conduction band maximum (valence band minimum) plus (minus) the band gap in Si.

III.2 Smearing parameter and conductivity sum rule

The macroscopic density-density and paramagnetic current-current response functions in Eq. (9) contain the infinitesimal $\eta$ which ensures causality and avoids having poles on the real axis. Numerically $\eta$ can be used as a smearing parameter to simulate finite lifetime effects of the excitations. As a result each peak in the absorption spectrum acquires a finite width. This is done by replacing

[TABLE]

in Eq. (9) for both the density-based and the current-based approach, where $\eta$ is now a finite number. In the current-based approach however the frequency $\omega$ enters also in the definition of the dielectric function, Eq. (6), as a factor $1/\omega^{2}$. Moreover it is common to numerically impose Berger et al. (2006); Raimbault et al. (2015) the CSR (Eq. 24) replacing the diamagnetic term, $N/V$, with minus the paramagnetic term, $-\chi^{{\rm IP}}_{j^{p}j^{p}}(\omega)$, evaluated at $\omega=0$. It is thus natural to wonder whether we should replace also $1/\omega^{2}$ by $1/z^{2}$ and if we should use $-\chi^{{\rm IP}}_{j^{p}j^{p}}(z)$ evaluated at $\omega=0$ or $z=0$ while imposing the CSR.

In Fig. 1 we plot the optical spectra (${\bf q}=0$) obtained for LiF and Si inserting $\chi^{{\rm IP}}_{\rho\rho}(z)$ in Eq. (7) and $\chi^{{\rm IP}}_{j^{p}j^{p}}(z)$ in Eq. (6). For the current-based approach we consider both the factors $1/\omega^{2}$ and $1/z^{2}$ in the definition of the dielectric function (Eq. (6)). In both cases the diamagnetic term is replaced by $-\chi^{{\rm IP}}_{j^{p}j^{p}}(i\eta)$, i.e. the CSR is imposed with $\omega=0$. The optical spectra obtained within the current-based approach are different from those obtained within the density-based approach. In particular they present clearly unphysical features: the case with the $1/\omega^{2}$ factor in Eq. (6) shows a divergent low energy contribution which resembles the Drude like behaviour of metals; the case with the $1/z^{2}$ factor shows a negative peak at $\omega=\eta$. Both unphysical features are related to the existence of a finite smearing in the low frequency region of the spectrum.

One possible solution is to adopt the recipe proposed by Cazzaniga et al. Cazzaniga et al. (2010) in another context, i.e. to choose $\tilde{z}=\sqrt{\omega^{2}+2i\eta\omega}$, instead of $z=\omega+i\eta$. The latter choice implies no smearing at $\omega=0$ and a CSR uniquely defined by $N/V=-\chi^{{\rm IP}}_{j^{p}j^{p}}(0)$. For the case of Si we show this recipe cures both unphysical features although some residual numerical noise remains.

A more rigorous solution requires to consider again the expansion defined by Eq.(21). Having a smearing parameter we can now expand either around $\omega=0$ or around $z=0$ (i.e. $\omega=-i\eta)$). The expansion around $\omega=0$ is more suited to analyze the case with the $1/\omega^{2}$ factor in Eq. (6) and yields

[TABLE]

The paramagnetic term entering $A_{\eta}^{{\rm IP}}$ is correctly balanced by replacing $N/V\rightarrow-\chi^{{\rm IP}}_{j^{p}j^{p}}(i\eta)$ in the diamagnetic term. However $B_{\eta}$ is not zero anymore (it is zero only for $\eta=0$) and thus leads to a divergence in the spectrum at $\omega=0$.

The expansion around $z=0$ is more suited to analyze the case with the $1/z^{2}$ factor in Eq. (6) and yields

[TABLE]

In this case $B^{IP}$ is numerically zero as expected theoretically, however the the paramagnetic term entering $A^{{\rm IP}}$ is not correctly balanced (since we are using $N/V\rightarrow-\chi^{{\rm IP}}_{j^{p}j^{p}}(i\eta)$) and the CSR is broken. $A^{{\rm IP}}$ is here multiplied by $1/z^{2}$ leading to a negative energy peak around $\omega=\eta$. To summarize, in order to avoid unphysical divergencies and negative peaks there are two options:

• •

the $1/\omega^{2}$ factor can be used in Eq. (6) together with the CSR imposed by $N/V\rightarrow-\chi^{{\rm IP}}_{j^{p}j^{p}}(i\eta)$ and balancing $B_{\eta}^{{\rm IP}}$ by a proper counter term;

• •

a $1/z^{2}$ term can be used in Eq. (6) together with the CSR imposed by $N/V\rightarrow-\chi^{{\rm IP}}_{j^{p}j^{p}}(0)$.

The latter option is the most straightforward to implement and we have tested that cures the unphysical negative peak at $\omega=\eta$ and gives the same spectra as within the density-based approach.

III.3 The conductivity sum rule and the Drude term

The current-current response function in Eq. (15) consists of two terms, a constant diamagnetic term $N/V$ and a paramagnetic term given by Eq. (9). In practice the sum over states in Eq. (9) is truncated. As a consequence the CSR in Eq. (24) is no longer satisfied, and it becomes impossible to converge the optical spectra for small frequencies. To solve these problems one can impose the CSR by replacing Kootstra et al. (2000); Romaniello and de Boeij (2005); Berger et al. (2005) the diamagnetic term $N/V$ with $-\chi_{j^{p}j^{p}}({\bf q},0)$ (or $-\chi_{j^{p}j^{p}}({\bf q},i\eta)$ as we did in the previous section). Thus the diamagnetic and paramagnetic contributions are treated on equal footing and no convergence problems occur.

At ${\bf q}={\bf 0}$ this strategy poses no problems for systems with a gap. For metals, instead, it suppresses the Drude tail, which is described by the term $A^{IP}/\omega^{2}$ on the right-hand side of Eq. (34). This occurs because the CSR in general does not hold at ${\bf q}={\bf 0}$ since intraband transitions are excluded in the sum over states in $\chi_{j^{p}j^{p}}$. We thus consider here the case of a metal, Cu, and compute its spectrum without imposing the CSR. Since the diamagnetic term is purely real, we need to use the $1/z^{2}$ strategy as discussed previously. However, without imposing the CSR, calculations never converge. This is shown in Fig. 2: the Drude tail in Cu has a wrong behaviour, whereas in Si, where converged calculations should give no absorption at low energy, an artificial Drude-like tail appears.

We found that the best solution is to enforce the CSR and to calculate the Drude term through the explicit inclusion of intraband transitions as described for the density-based approach. A faster convergence of the latter contribution can be obtained using the tetrahedron method for the integration in the Brillouin zone. Blöchl et al. (1994); Romaniello and de Boeij (2005)

III.4 Many-body lifetimes

In Sec. III.2 the finite width of the peaks was obtained by introducing an ad hoc smearing parameter. A more physically motivated approach is to consider finite lifetimes $\gamma$ originating from the imaginary parts of the many-body self-energy. For example, it has been shown that in a semiconductor such as Si the finite width of the peaks is well-described by means of the Fan and $GW$ self-energies, which account for electron-phononGiustino (2016); Marini (2008) and electron-electron Aryasetiawan and Gunnarsson (1998) scattering processes, respectively. To illustrate the effects of such lifetimes in the following we will consider the contributions due to the Fan self-energy. We define modified KS energies $\tilde{\epsilon}^{{\rm KS}}_{n}$ as

[TABLE]

and use them to compute the absorption both in the current-based and in the density-based approach. The KS energies $\tilde{\epsilon}^{{\rm KS}}_{n}$ correspond to the eigenenergies of the Hamiltonian (18), in which only the imaginary part of the Fan self-energy is considered. The underlying Hamiltonian, which we will indicate as $\tilde{H}_{{\rm KS}}$, is diagonal in the same basis set in which $H_{{\rm KS}}$ is diagonal, i.e. $\tilde{\psi}^{{\rm KS}}_{n}\approx\psi^{{\rm KS}}_{n}$ and thus $\tilde{x}^{{\rm KS}}_{nm}=x^{{\rm KS}}_{nm}$. The velocity matrix elements however change because the velocity operator is proportional to the hamiltonian itself. Using Eq. (4), the velocity matrix elements read

[TABLE]

Such expression generalizes the result by Del Sole et al.Del Sole and Girlanda (1993) (that we use in the next section) to the case of complex energies and is consistent with the findings of Tokman. 151515Ref. Tokman, 2009 discusses the connection between the time derivative of the dipole operator and the velocity operator in case of a dephasing of the polarization (i.e. an imaginary term in the Hamiltonian) is considered.

Moreover, since $\tilde{H}_{{\rm KS}}$ is non-Hermitian (although it is still a normal matrix), the velocity operator is also non-Hermitian. This means in practice

[TABLE]

and thus the numerator of the response function (Eq. (9)) in the current-based approach cannot be written as a square modulus $|\tilde{v}^{{\rm KS}}_{nm}|^{2}$ anymore. One may wonder if—having a non positive defined numerator—the spectrum may become negative.

In Fig. 3 we show for Cu and Si that this is not the case: the spectrum is well defined and matches the one obtained in the density-based approach as long as the non-Hermiticity of the velocity operator is correctly taken into account. If instead we ignore the non-Hermicity of the velocity operator and use $|\tilde{v}^{{\rm KS}}_{nm}|^{2}$ in the numerator of Eq. (9), a spurious divergence appears at low energy for Si and the Drude tail is not correctly described for Cu.

III.5 Quasiparticle energies and excitonic effects

We finally consider the dielectric function beyond the IP approximation by using the Bethe-Salpeter equation. We include both the $GW$ corrections to the KS band structure and the effect of the electron-hole interaction in the absorption. As for the case of the QP lifetimes we have to consider a renormalized velocity operator, Del Sole and Girlanda (1993)

[TABLE]

which corresponds to the quasiparticle Hamiltonian (Eq.(18)) with a $GW$ self-energy. However this is not the only correction to be considered. If the dielectric function is computed by diagonalizing the excitonic Hamiltonian given by Eq. (43), then it is expressed in terms of the excitonic dipole matrix elements $x^{\rm Exc}_{\lambda}$ in the density-based approach,

[TABLE]

or via the excitonic velocity matrix elements $v^{\rm Exc}_{\lambda}$, in the current-based approach. One might be tempted, in analogy with Eq. (59), to define

[TABLE]

However this expression—which has been reported in the literature Rohlfing and Louie (2000); Spataru et al. (2005); Palummo et al. (2015)—is not correct since it does not correspond to the underlying excitonic Hamiltonian $H^{{\rm Exc}}$ given by Eq. (43). The excitonic velocity hence must read

[TABLE]

which in general is different from Eq. (60).

In Fig. 4 we compare the absorption spectrum in the density-based approach with the current-based approach results with and without the renormalized velocity operator, i.e. considering the corrections due to Eqs.(58) and (62). We notice that the two corrections in general partially cancel each other since $\epsilon^{{\rm QP}}_{n}-\epsilon^{{\rm QP}}_{m}>\epsilon^{{\rm KS}}_{n}-\epsilon^{{\rm KS}}_{m}$ and $E_{\lambda}<\epsilon^{{\rm QP}}_{n}-\epsilon^{{\rm QP}}_{m}$ and the resulting error in the absorption intensity is proportional to $|E_{\lambda}-(\epsilon^{{\rm KS}}_{n}-\epsilon^{{\rm KS}}_{m})|/|\epsilon^{{\rm KS}}_{n}-\epsilon^{{\rm KS}}_{m}|$. There is thus an error cancellation in the intensity due to the fact that the onset of the KS absorption is often not too different from the onset of the BSE absorption. The opposite sign of the two corrections is evident when they are considered independently. If only Eq. (58) is considered, the absorption spectrum intensity is strongly overestimated. Conversely the spectrum intensity is strongly underestimated if only Eq. (62) is considered. Only when including both the renormalization of the velocity due the quasiparticle corrections and the nonlocal operator we recover the results obtained within the density-based approach.

Regarding the computational cost, we note that if Eq. (62) is used to evaluate the excitonic velocity operator $v^{{\rm Exc}}_{\lambda}$, then iterative schemes to calculate the dielectric function, such as the Lanczos-Haydock method Haydock (1980); Benedict and Shirley (1999); Grüning et al. (2009), are precluded since $v^{{\rm Exc}}_{\lambda}$ in Eq. (62) explicitly depends on the BSE energies $E_{\lambda}$. Instead, one should use Eq. (61) and explicitly calculate the commutator of the dipole operator and the excitonic Hamiltonian $H^{{\rm Exc}}$ which is cumbersome since $H^{{\rm Exc}}$ contains nonlocal kernels. This issue also arises in approaches in which the QP corrections are calculated in a non-perturbative way, i.e., when a non-diagonal self-energy is considered, and thus the commutator with such a self-energy should in principle be computed. Bruneval et al. (2006)

IV Direct calculation of the induced density, polarization and current

For sake of completeness we extend the discussion by considering alternative approaches in which optical response functions are not computed as sum over states (Eq. 9), but rather using the right-hand side of Eqs.(5): the change in the density or in the current-density induced by a perturbing field Bertsch et al. (2000); Grüning et al. (2016); Yabana et al. (2012) is computed first and the macroscopic response functions are then obtained by dividing the induced quantities by the perturbing fields.

A way to access the density and the current-density is through the one particle time-dependent density-matrix of the system $\varrho({\mathbf{x}},{\mathbf{x}}^{\prime},t)$.

The general equation of motion (EOM) for the density matrix can be written in the Von Neumann form Von Neumann (1927). It reads

[TABLE]

where the Liouvillian operator $\hat{L}(t)$, defined by

[TABLE]

describes the coherent evolution in terms of the Hermitian part of the Hamiltonian $\hat{h}_{+}(t)$. The relaxation operator $\hat{R}$, defined by

[TABLE]

describes relaxation processes in terms of the skew-Hermitian part of the Hamiltonian $\hat{h}_{-}(t)$.161616This EOM can also be derived as an approximation to the non-equilibrium Green’s function theory or quantum kinetic theory Kadanoff and Baym (1962); Bonitz (1998) where the density matrix is the time diagonal of the lesser Green function $\varrho(t)=-iG^{<}(t,t)$. We can solve the equation for $\varrho(t)$ in the equilibrium one-particle basis set $\varrho_{nm}(t)=\langle\psi_{m}|\varrho(t)|\psi_{n}\rangle$ with the initial condition $\varrho_{nm}(t=0)=\delta_{nm}f_{n}$. The density and current-density can then be obtained as

[TABLE]

If the Hamiltonian is expressed in terms of quasiparticle energies and only the variation of the static screened exchange self-energy is considered, this approach has been proven Attaccalite et al. (2011), to linear order, to be equivalent to the BSE. In this case the use of a smearing parameter $\eta$ in the BSE corresponds to setting the non-Hermitian part of the Hamiltonian proportional to $\eta$.

Since the Hamiltonian is gauge dependent, so is the EOM in Eq. (63) and its solution $\varrho(t)$. Using the gauge function $\Lambda({\mathbf{x}},t)$ defined in App. C, which transforms the potentials from one gauge to another, $\hat{L}(t)$, ${\hat{R}}(t)$, and $\varrho(t)$ transform according to

[TABLE]

where $\chi=c\,[\Lambda({\mathbf{x}},t)-\Lambda({\mathbf{x}}^{\prime},t)]$. One can show that a gauge transformation preserves the current, the density, the polarization and the electromagnetic energy Tokman (2009) only if all quantities are gauge transformed together. We note that in the velocity gauge the relaxation operator acquires a dependence on the perturbing potential, not present in the length gauge.Tokman (2009); Lamb et al. (1987). At finite momentum Eqs. (67)-(69) are all is needed to prove the equivalence between the density- and current-based approaches for the description of optical properties, similarly to Eq. (8) for the response functions. Propagating the density-matrix in the length gauge (with potential $\phi(t)$) and then computing the density variation the macroscopic density-density response function is obtained. Propagating the density-matrix in the velocity gauge (with potential ${\mathbf{A}}(t)$) and then computing the current-density variation the macroscopic current-current response function is obtained.

Again difficulties arise if one considers the optical limit ${\bf q}\rightarrow{\bf 0}$. Indeed at ${\bf q}={\bf 0}$ the length gauge cannot be formulated in terms of the coupling with the density only Gonze et al. (1995); Martin and Ortiz (1997). The induced macroscopic polarization needs to be defined. Up to first order in the perturbing field, it can be constructed 171717for a discussion and a general expression, valid at all orders, see Ref. Attaccalite and Grüning, 2013 from the density matrix as

[TABLE]

The longitudinal 181818While the density-density response function can be used only to compute the longitudinal response, the dipole-dipole response function, like the current-current one, can also describe transverse and mixed (longitudinal-transverse) terms. dipole-dipole response function can then be defined as ${\chi_{dd}(\omega)=\delta P(\omega)/\delta E(\omega)}$, where $E$ is the total macroscopic electric field. Thanks to Eq. (70) it is possible to show the formal equivalence Tokman (2009) between the two formalisms also at ${\bf q}=0$ for cold semiconductors. As previously, $\chi_{dd}(\omega)$ describes however only the inter-band and not the intra-band contribution. The latter can be obtained, in the length gauge, only by explicitly considering the ${\bf q}\rightarrow{\bf 0}$ limit. The direct numerical calculation of $\delta\rho({\bf q},\omega)$ at very small ${\bf q}$ is however not feasible, and one would need to analytically deal with the small ${\bf q}$ dependence as we do in App. E for the response function. Instead, the current-based formalism can be directly formulated also at ${\mathbf{q}}={\bf 0}$ and the intraband contribution is also included. However, similar to what discussed previously, one must be careful because of potential breaking of the CSR (Eq. 24).

V Conclusions

We compared the optical absorption of extended systems calculated from the density-density and current-current linear response functions obtained within many-body perturbation theory. We did this by studying the longitudinal macroscopic dielectric function both on a formal and on a numerical level.

We showed that for a finite momentum ${\bf q}$, carried by the external perturbation, the two approaches are formally equivalent thanks to the continuity equation, which relates the density to the longitudinal current.

For the case of vanishing momentum, ${\bf q}={\bf 0}$, the optical absorption is not well defined in terms of the density-density response function. A small ${\bf q}$ expansion is needed, which leads to a formulation in terms of the dipole-dipole response function plus a divergent term which describes intraband transitions for metalic systems. The current-based approach is instead exact at ${\bf q}={\bf 0}$. In practice, however, one needs to impose the conductivity sum rule in a way which suppresses the intraband transitions at ${\bf q}={\bf 0}$, thus making the small ${\bf q}$ expansion needed for metalic systems also in the current based approach.

When gapped systems are considered there are no intraband transitions and the two approaches are formally equivalent once the conductivity sum rule is imposed in the current-based approach. On the other hand we showed that the formal equivalence of the approaches may be lost in practical calculations when the strategies usually employed in the density-based approach to include smearing, quasiparticle lifetimes and the electron-hole interaction are naively applied to the current-based approach. The smearing is included straightforwardly in the density-based approach by replacing the real frequency with an imaginary frequency. However, a careless extension of this recipe within the current-based approach leads to unphysical features in the optical spectrum. We showed how correctly include the smearing by redefining the dielectric function and the conductivity sum rule. The inclusion of lifetimes, quasiparticle corrections, and excitonic effects correspond to a change in the underlying Hamiltonian. Therefore the velocity operator, which enters the definition of the current-current response function, has to be modified accordingly. We noted instead that the expression for the excitonic velocity operator reported in several published works is incorrect. In this work we thus report the correct general definition for the velocity operator when complex energies, quasiparticle corrections, and excitonic effects are taken into account.

Acknowledgments

The authors thank the EUspec COST Action for Short Term Scientific Missions that allowed to initiate and carry out this project. Discussion within the Collaboration Team on Correlation of the European Theoretical Spectroscopy Facility (ETSF) is greatly acknowledged. DS acknowledge Giovanni Onida for the computational time provided on the HPC cluster “ETSFMI” in Milano and Paolo Salvestrini for the support on the cluster; computing time has been also provided by the French national GENGI-IDRIS supercomputing centers at Orsay under contract $n^{o}$ t2012096655; financial support by the Futuro in Ricerca grant No. RBFR12SW0J of the Italian Ministry of Education, University and Research MIUR; the funding received from the European Union project MaX Materials design at the eXascale H2020-EINFRA-2015-1, Grant agreement n. 676598 and Nanoscience Foundries and Fine Analysis - Europe H2020-INFRAIA-2014-2015, Grant agreement n. 654360. JAB and PR acknowledge support from the IDEX Emergence (project no. 2016-075/CNRS)

Appendices

Appendix A Fourier Transform in Periodic Boundary Conditions

Here we explicitly introduce the Fourier transform in Periodic Boundary Conditions (PBC). To make the derivation more clear, we use, just for the present appendix, different symbols for the different functions in real and reciprocal space associated to a given observable $O$. We start from the general definition:

[TABLE]

We then divide the space in a series of microscopic unit cells with volume $V$ which are periodically repeated and we use ${\mathbf{x}}={\mathbf{r}}+{\mathbf{R}}$ with ${\mathbf{r}}$ restricted to the volume $V$ centered at ${\mathbf{R}}={\bf 0}$. Accordingly the reciprocal space results separated in two parts with ${\mathbf{G}}$ defined by ${\mathbf{G}}\cdot{\mathbf{R}}=2\pi\,n$, where $n$ is any integer number, and ${\mathbf{q}}={\mathbf{k}}-{\mathbf{G}}$ restricted to the first Brillouin Zone with volume ${\Omega=(2\pi)^{3}/V}$. Then the Fourier transform becomes

[TABLE]

A.0.1 Periodic quantities

A periodic function $O({\mathbf{r}})$ is characterized by the property ${O({\mathbf{R}}+{\mathbf{r}})=O({\mathbf{r}})}$. Substitution of this identity into Eq. (72b) leads to

[TABLE]

i.e., the Fourier transform of a periodic function has non-zero components only for ${\mathbf{q}}={\bf 0}$. We can use this result in Eq. (72a) to obtain

[TABLE]

Thus the Fourier transforms in Eq. (72) reduce to

[TABLE]

and the macroscopic part of a periodic quantity, i.e., its average over the unit cell, is $\tilde{\mathcal{O}}({\mathbf{G}}={\bf 0})/V$.

A.0.2 Macroscopic quantities

Here we are interested in processes in which the transferred momentum ${\mathbf{q}}+{\mathbf{G}}$ is small. Let us therefore consider a (non periodic) function for which, in reciprocal space, only the ${\mathbf{G}}={\bf 0}$ component is non-zero, i.e., ${\mathcal{O}({\mathbf{q}}+{\mathbf{G}})=\delta_{{\mathbf{G}},{\bf 0}}\overline{\mathcal{O}}({\mathbf{q}})}$. Its Fourier transform in real space, given by

[TABLE]

depends both on ${\mathbf{R}}$ and ${\mathbf{r}}$, but with a smooth dependence on ${\mathbf{r}}$ since it has no fast oscillating ${\mathbf{G}}$ component. Its inverse Fourier transform is

[TABLE]

If $e^{-i{\mathbf{q}}\cdot{\mathbf{r}}}\simeq 1$ we can neglect the ${\mathbf{r}}$ dependence and consider

[TABLE]

We thus obtain

[TABLE]

The relation $e^{-i{\mathbf{q}}\cdot{\mathbf{r}}}\simeq 1$ holds for small ${\mathbf{q}}$, i.e. in the “macroscopic limit”, where one can consider the volume $V$ as infinitesimal ($V\rightarrow d^{3}{\mathbf{R}}$ and $\Omega\rightarrow\mathbb{R}^{3}$) and ${\mathbf{R}}$ becomes a continuous variable. Indeed we can assume ${\mathbf{q}}$ or small $V$ if $2\pi/q\gg V^{1/3}$. In this limit materials are considered as a continuum, their atomistic structure is neglected and the macroscopic integrals used to describe electromagnetism in classical media are recovered. Instead for short wavelengths, i.e. $2\pi/q\ll V^{1/3}$, the atomistic structure can never be neglected and one is forced to use Eqs. (76)-(77) when computing the response induced at ${\mathbf{G}}={\bf 0}$ in PBC.

Appendix B Full dielectric tensor

In general the dielectric tensor $\underline{\underline{\varepsilon}}^{M}({\bf q},\omega)$ can be decomposed into a longitudinal component, a transverse component, and mixed components as

[TABLE]

For isotropic systems, in particular, there exist only two independent components, i.e., the longitudinal one ($\varepsilon^{L}$) and transverse one ($\varepsilon^{T}$), so that the dielectric function reads

[TABLE]

The longitudinal component describes the response to longitudinal fields, which are involved, for example, in electron energy-loss experiments, where the scattering cross-section of an electron traversing a medium is proportional to $-Im\left\{1/\varepsilon^{L}({\bf q},\omega)\right\}$. The transverse component describes the response to optical fields, which are characterized by small $q\approx\omega/c\approx 0$. In the long-wavelength limit $q\rightarrow 0$ the two quantities are equal Nozières and Pines (1999); Adler (1962), thus optical and energy-loss measurements contain the same physical information. In particular for cubic symmetry we have

[TABLE]

Since in the present manuscript we only deal with longitudinal perturbations and longitudinal external fields, when we write a scalar dielectric function this is understood to be its longitudinal component, i.e. ${\varepsilon({\bf q},\omega)=\varepsilon^{L}({\bf q},\omega)}$. The same applies to other quantities such as the current-current response function ${\chi_{jj}({\bf q},\omega)=\chi^{L}_{{\mathbf{j}}{\mathbf{j}}}({\bf q},\omega)=\chi_{j^{L}j^{L}}({\bf q},\omega)}$, the external potential ${A({\bf q},\omega)={\mathbf{A}}^{L}({\bf q},\omega)}$, the current ${j({\bf q},\omega)={\mathbf{j}}^{L}({\bf q},\omega)}$, the polarization ${P({\bf q},\omega)={\mathbf{P}}^{L}({\bf q},\omega)}$, and the dipoles ${d_{\bf q}={\mathbf{d}}^{L}_{\bf q}}$.

Appendix C Gauges

The scalar and vector potentials $\phi({\mathbf{x}},t)$ and ${\mathbf{A}}({\mathbf{x}},t)$, respectively, describe a general electromagnetic field by

[TABLE]

The electromagnetic field is invariant under the gauge transformations

[TABLE]

Here $\Lambda(\mathbf{x},t)$ is a differentiable, but, otherwise, arbitrary function of $\mathbf{x}$ and $t$. Notice that in quantum mechanics the gauge transformation also modifies the wave function phase as:

[TABLE]

One can use the gauge freedom to map a problem in an equivalent one, which is maybe easier to solve. For example, one can completely gauge transform the scalar potential $\phi$ into a vector potential of the form

[TABLE]

using $\Lambda({\mathbf{x}},t)=c\int_{0}^{t}\phi({\mathbf{x}},t^{\prime})dt^{\prime}$. This is the Weyl gauge. Such a vector potential, being expressed as the gradient of a scalar, is longitudinal, i.e. it describes a longitudinal vector field, since its Fourier transform is parallel to $\mathbf{q}$ for any $\mathbf{q}$. Using $\Lambda({\mathbf{r}},t)=-\int A^{L}({\mathbf{x}}^{\prime},t)d^{3}{\mathbf{x}}^{\prime}$, one can gauge transform the longitudinal component of ${\mathbf{A}}({\mathbf{x}},t)$ to a scalar potential as

[TABLE]

This is the Coulomb gauge defined by $\mathbf{\nabla}\cdot{\mathbf{A}}({\mathbf{x}},t)=0$. In the dipole approximation, i.e. $\mathbf{E}({\mathbf{x}},t)\approx\mathbf{E}({\bf 0},t)$, the two gauges reduce to the length and the velocity gauge.

Gauge transformation affects only the potentials describing longitudinal fields as it is clear if Eq. (LABEL:Eqn:gauge_transformation) is written in Fourier space. Indeed the transverse part of the vector potential is gauge independent and can never be described in terms of a scalar potential.

Appendix D Macroscopic response functions

In a system with translation invariance symmetry only for ${\mathbf{x}}={\mathbf{R}}$, the longitudinal potentials $\delta\phi^{{\mathbf{G}}^{\prime}}({\mathbf{q}},\omega)$ and $\delta A^{{\mathbf{G}}^{\prime}}({\mathbf{q}},\omega)$ in general induce variations at any ${\mathbf{q}}+{\mathbf{G}}$ component with ${\mathbf{G}}$ also different from ${\mathbf{G}}^{\prime}$ (we write here the dependence on ${\mathbf{G}}$ as a superscript for convenience). Indeed we can formally write the component of the induced density and current density linear in the perturbing potentials of form (2a)-(2b) (using either the length or the velocity gauge) as

[TABLE]

where ${\mathbf{q}}$ is now restricted to the first Brillouin zone, ${\mathbf{G}}$ is a reciprocal lattice vector, and

[TABLE]

If we consider a perturbation with only the ${\mathbf{G}}^{\prime}={\bf 0}$ component, i.e. $\delta\phi^{{\mathbf{G}}^{\prime}}({\mathbf{q}},\omega)=\delta\phi({\mathbf{q}},\omega)\,\delta_{{\mathbf{G}}^{\prime},{\bf 0}}$ and $\delta{\mathbf{A}}^{{\mathbf{G}}^{\prime}}({\mathbf{q}},\omega)=\delta{\mathbf{A}}({\mathbf{q}},\omega)\,\delta_{{\mathbf{G}}^{\prime},{\bf 0}}$ , and look for the variation of the macroscopic induced density and current density, i.e. their ${\mathbf{G}}={\bf 0}$ Fourier component, we arrive at (5a) and (5b).

Appendix E intraband contribution to C

Following similar steps as in Ref. [Romaniello and de Boeij, 2005] the intraband contribution to $\bar{\chi}^{IP}_{dd}$ becomes

[TABLE]

where $\sum_{{\bf k}}$ has been replaced by $V/(2\pi)^{3}\int d{\bf k}$ and the integration over the $k$-space reduced to an integral over the sheets $S_{i}$ of the Fermi surface originating by the partially occupied bands $i$. Here $\varv^{F}_{i{\bf k}}$ is the Fermi velocity. For the frequency-dependent factor we can use the Cauchy theorem and write

[TABLE]

In optical experiments $\omega/q$ is of the order of the velocity of light $c$. 191919The direction of the ${\bf q}\rightarrow{\bf 0}$ limit is a delicate point. Indeed the $({\bf 0},0)$ point in the $({\bf q},\omega)$ plane is non analytic, i.e. the value of the dielectric function depends on the direction of the limit $({\bf q},\omega)\rightarrow({\bf 0},0)$. This direction is determined by the experiment we would like to describe. In optical experiments the direction of interest is $\omega=c{\bf q}$. Other directions in the $({\bf q},\omega)$ plane may be of interest. For example in electron-energy loss experiments one measures the inverse dielectric function $\varepsilon^{-1}({\bf q},\omega)$ at fixed momentum ${\bf q}_{exp}$. The line ${\bf q}={\bf q}_{exp}$ always crosses the $\varv^{F}({\bf q})$ line, thus a peak must always appear if we derive $\varepsilon({\bf q},\omega)$ from the computed electron energy loss function at ${\bf q}_{exp}$. Therefore for ${\bf q}_{exp}\rightarrow{\bf 0}$ we always obtain the Drude-like tail in the absorption spectrum, also without smearing. Therefore, a Drude-like peak in the absorption can be described only if $|\varv^{F}_{n{\bf k}}|\approx c$. At the IP level, using real energies, this never happens; therefore the imaginary part in Eq. (90) is zero and the real part reduces to -1 in the limit of ${q\rightarrow 0}$ (and finite $\omega$). In this case the intraband contribution to the dielectric function is real and reads

[TABLE]

In the excitonic case the treatment of the intraband contribution is more involved and requires also the Taylor expansion of $E_{\lambda}({\bf q})$ and $A_{\lambda,{\bf q}}$.