Quantum plasmons and intraband excitons in doped nanoparticles: Failure of the Tamm-Dancoff approximation and importance of electron-hole attraction

TL;DR

This study investigates the electronic excitations in doped semiconductor nanoparticles, highlighting the failure of the Tamm-Dancoff approximation and emphasizing the importance of electron-hole attraction for accurate modeling of plasmons and excitons.

Contribution

It demonstrates the limitations of common single-excitation theories and introduces a schematic model that accurately reproduces numerical and experimental results for doped nanoparticles.

Findings

Time-dependent Hartree-Fock best describes excitation character.

Excitations transition from confinement to excitonic to plasmonic with increasing electrons.

Electron-hole attraction and de-excitations are crucial for accurate descriptions.

Abstract

We use excited-state quantum chemistry techniques to investigate the intraband absorption of doped semiconductor nanoparticles as a function of doping density, nanoparticle radius, and material properties. The excess electrons are modeled as interacting particles confined in a sphere. We compare the predictions of various single-excitation theories, including time-dependent Hartree-Fock, the random-phase approximation, and configuration interaction with single excitations. We find that time-dependent Hartree-Fock most accurately describes the character of the excitation, as compared to equation-of-motion coupled-cluster theory with single and double excitations. The excitation evolves from confinement-dominated, to excitonic, to plasmonic with increasing number of electrons at fixed density, and the threshold number of electrons to produce a plasmon increases with density due to quantum confinement. Exchange integrals (attractive electron-hole interactions) are essential to properly describe excitons, and de-excitations (i.e.~avoidance of the Tamm-Dancoff approximation) are essential to properly describe plasmons. We propose a schematic model whose analytic solutions closely reproduce our numerical calculations. Our results are in good agreement with experimental spectra of doped ZnO nanoparticles at a doping density of $1.4\times 10^{20}$ cm$^{-3}$.

Equations

$$\begin{array}{rl}H& =\sum _{n=1}^N\left[-\frac{1}{2m^{\ast }}{\nabla }_n^2+v({\mathbf{r}}_n)\right]+\sum _{n=1}^N\sum _{m<n}\frac{1}{|{\mathbf{r}}_m-{\mathbf{r}}_n|},\end{array}$$

$${\varphi }_{nlm}(r,\theta ,\varphi )=N_{nl}^{-1}{j}_l({\alpha }_{nl}r)Y_{lm}(\theta ,\varphi )=R_{nl}(r)Y_{lm}(\theta ,\varphi ),$$

$${\epsilon }_{nlm}=\frac{k_{nl}^2}{2m^{\ast }{R}^2},$$

$$H=\sum _p{\epsilon }_p{a}_p^{†}{a}_p+\frac{1}{2}\sum _{pqrs}⟨pq|rs⟩a_p^{†}{a}_q^{†}{a}_s{a}_r,$$

$$⟨pq|rs⟩=\int d{\mathbf{x}}_1\int d{\mathbf{x}}_2{\varphi }_p^{\ast }({\mathbf{x}}_1){\varphi }_q^{\ast }({\mathbf{x}}_2)r_{12}^{-1}{\varphi }_r({\mathbf{x}}_1){\varphi }_s({\mathbf{x}}_2),$$

$$|{\Psi }_n⟩=\sum _{ai}\left[X_{ai}{a}_a^{†}{a}_i+Y_{ai}{a}_i^{†}{a}_a\right]|{\Psi }_0⟩$$

$$\left(\begin{array}[]{cc}\mathbf{A}&\mathbf{B}\\ -\mathbf{B}^{*}&-\mathbf{A}^{*}\end{array}\right)\left(\begin{array}[]{c}\mathbf{X}\\ \mathbf{Y}\end{array}\right)=\left(\begin{array}[]{c}\mathbf{X}\\ \mathbf{Y}\end{array}\right)\mathbf{\Omega},$$

$$A_{ai,bj}$$

$$B_{ai,bj}$$

$$\alpha (\omega )=\sum _m|⟨{\Psi }_0|\hat{\mu }|{\Psi }_m⟩{|}^2\delta (\omega -{\Omega }_m)$$

$$\begin{array}{rl}\delta \rho (\mathbf{r},\omega )& =\int d{\mathbf{r}}^{\prime }\chi (\mathbf{r},{\mathbf{r}}^{\prime },\omega )v_{\mathrm{ext}}({\mathbf{r}}^{\prime },\omega )\\ & =\sum _n\frac{2{\Omega }_n}{{\omega }^2-{\Omega }_n^2}\sum _{ai}|X_{ai}^n+Y_{ai}^n{|}^2⟨i|\hat{z}|a⟩{\varphi }_a(\mathbf{r}){\varphi }_i(\mathbf{r}),\end{array}$$

$$\begin{array}{rl}{\Omega }_n& =2\sum _{ai}({\epsilon }_a-{\epsilon }_i)|X_{ai}+Y_{ai}{|}^2\\ & +4\sum _{abij}⟨ib|aj⟩\left(X_{bj}+Y_{bj}\right)\left(X_{ai}+Y_{ai}\right)\\ & -2\sum _{abij}⟨ib|ja⟩\left(X_{ai}{X}_{bj}+Y_{ai}{Y}_{bj}\right)\\ & -2\sum _{abij}⟨ij|ba⟩\left(X_{ai}{Y}_{bj}+Y_{ai}{X}_{bj}\right).\end{array}$$

$$\left[{\Omega }_n-({\epsilon }_a-{\epsilon }_i)\right]X_{ai}^n=\lambda {\rho }_{ai}\sum _{bj}{\rho }_{bj}{X}_{bj}^n,$$

$$\frac{1}{\lambda }=\sum _{ai}\frac{{\rho }_{ai}^2}{{\Omega }_n-({\epsilon }_a-{\epsilon }_i)}.$$

$$\Omega =\epsilon +\lambda \sum _{ai}{\rho }_{ai}^2\approx \epsilon +\lambda N_{\mathrm{trans}}{\bar{\rho }}^2,$$

$${\epsilon }_{\mathrm{g}}^{\mathrm{NI}}(R)=\frac{6.55}{2m^{\ast }{R}^2}+\frac{5.04\sqrt{\rho }}{2m^{\ast }\sqrt{R}}.$$

$${\Omega }^{\mathrm{RPA}(\mathrm{TDA})}(R)={\epsilon }_{\mathrm{g}}^{\mathrm{HF}}(R)+c\sqrt{\rho R},$$

$$\frac{1}{\lambda }=\sum _{ai}{\rho }_{ai}^2\frac{2({\epsilon }_a-{\epsilon }_i)}{{\Omega }_n^2-({\epsilon }_a-{\epsilon }_i{)}^2}.$$

$${\Omega }^{\mathrm{RPA}}(R)=\sqrt{{\left({\epsilon }_{\mathrm{g}}^{\mathrm{HF}}(R)\right)}^2+2\sqrt{\rho R}{\epsilon }_{\mathrm{g}}^{\mathrm{HF}}(R)}.$$

$$X_{ai}^{\mathrm{RPA}(\mathrm{TDA})}=\frac{1}{\sqrt{C}}\times {\rho }_{ai}$$

$$X_{ai}^{\mathrm{RPA}}$$

$$Y_{ai}^{\mathrm{RPA}}$$

$$C=4\sum _{ai}|{\rho }_{ai}{|}^2{\epsilon }_{\mathrm{g}}{\Omega }_{\mathrm{P}}$$

$${\Omega }^{\mathrm{CIS}}(R)={\epsilon }_{\mathrm{g}}^{\mathrm{HF}}(R)+\sqrt{\rho R}-\frac{1.4}{R}$$

$${\Omega }^{\mathrm{TDHF}}(R)=\sqrt{{\left({\epsilon }_{\mathrm{g}}^{\mathrm{HF}}(R)-\frac{1.4}{R}\right)}^2+2\sqrt{\rho R}\left({\epsilon }_{\mathrm{g}}^{\mathrm{HF}}(R)-\frac{1.4}{R}\right)}.$$

$$\begin{array}{rl}& ⟨n_1{l}_1{\mu }_1;n_2{l}_2{\mu }_2|n_3{l}_3{\mu }_3;n_4{l}_4{\mu }_4⟩\\ & =\sum _{l=0}^{\infty }{R}^l(n_1{l}_1,n_2{l}_2;n_3{l}_3,n_4{l}_4)\sum _{\mu =-l}^l⟨l\mu |l_1{\mu }_1|l_2{\mu }_2⟩⟨l\mu |l_3{\mu }_3|l_4{\mu }_4⟩,\end{array}$$

$$\begin{array}{rl}⟨l\mu |l_1{\mu }_1|l_2{\mu }_2⟩& =\int d\Omega y_{l\mu }{y}_{l_1{\mu }_1}{y}_{l_2{\mu }_2}\\ & =\sum _{m_1{m}_2{m}_3}[U_{m\mu }^l{]}^{\ast }{U}_{m_1{\mu }_1}^{l_1}{U}_{m_2{\mu }_2}^{l_2}⟨lm|l_1{m}_1|l_2{m}_2⟩,\end{array}$$

$$\begin{array}{rl}⟨& lm|l_1{m}_1|l_2{m}_2⟩=\int d\Omega {\left[Y_l^m\right]}^{\ast }{Y}_{l_1}^{m_1}{Y}_{l_2}^{m_2}\\ & =\sqrt{(2l_1+1)(2l_2+1)}(-1{)}^m\left(\begin{array}{ccc}l& l_1& l_2\\ 0& 0& 0\end{array}\right)\left(\begin{array}{ccc}l& l_1& l_2\\ -m& m_1& m_2\end{array}\right),\end{array}$$

$$\begin{split}&R^{l}(n_{1}l_{1},n_{2}l_{2};n_{3}l_{3},n_{4}l_{4})=\\ &\int_{0}^{R}dr_{1}\int_{0}^{R}dr_{2}r_{1}^{2}r_{2}^{2}\frac{r_{<}^{l}}{r_{>}^{l+1}}R_{n_{1}l_{1}}^{*}(r_{1})R_{n_{2}l_{2}}^{*}(r_{2})R_{n_{3}l_{3}}(r_{1})R_{n_{4}l_{4}}(r_{2})\\ &=\frac{1}{R}\Bigg{[}\int_{0}^{1}dx_{1}\int_{0}^{x_{1}}dx_{2}x_{1}^{2}x_{2}^{2}\frac{x_{2}^{l}}{x_{1}^{l+1}}u_{n_{1}l_{1}}^{*}(x_{1})u_{n_{2}l_{2}}^{*}(x_{2})u_{n_{3}l_{3}}(x_{1})u_{n_{4}l_{4}}(x_{2})\\ &\hskip 10.00002pt+\int_{0}^{1}dx_{1}\int_{x_{1}}^{1}dx_{2}x_{1}^{2}x_{2}^{2}\frac{x_{1}^{l}}{x_{2}^{l+1}}u_{n_{1}l_{1}}^{*}(x_{1})u_{n_{2}l_{2}}^{*}(x_{2})u_{n_{3}l_{3}}(x_{1})u_{n_{4}l_{4}}(x_{2})\Bigg{]},\end{split}$$

Full text

I Introduction

Metallic nanoparticles are an important optoelectronic platform because of their strong plasmon resonance at visible energies, which can be tuned by size, shape, and environmentWillets and Van Duyne (2007); Odom and Schatz (2011); however, the accessible carrier densities are limited to those of the parent metals and are typically $10^{22}$ cm*-3* or higher. Recently, the doping of semiconductor nanoparticles has enabled access to much lower electron densities Mocatta et al. (2011); Luther et al. (2011); Scholl, Koh, and Dionne (2012); Schimpf et al. (2014, 2015); Shen and Guyot-Sionnest (2016). These doped nanoparticles, with their tunable charge carrier density, exhibit strong intraband absorption in a wide range of the THz regime, enabling promising new infrared and plasmonic applications Scotognella et al. (2013); Faucheaux, Stanton, and Jain (2014); Kriegel, Scotognella, and Manna (2017); Lhuillier and Guyot-Sionnest (2017).

Chemical techniques have enabled semiconductor nanoparticles to be doped with as few as 1-100 electrons. In this regime, the classical electrostatic picture of plasmons breaks down, demanding a theory of so-called “quantum plasmons” Pitarke et al. (2007); Morton, Silverstein, and Jensen (2011); Bernadotte, Evers, and Jacob (2013); Tame et al. (2013); Krauter et al. (2014); Casanova, Matxain, and Ugalde (2016); Varas et al. (2016); Bursi et al. (2016); Zhang et al. (2017). Furthermore, the collective plasmon picture becomes dubious for systems containing only 1-10 excess electrons, suggesting a transition to “single-particle” excitations Shen and Guyot-Sionnest (2016) or – as we will argue – excitonic transitions. Here, we aim to present a detailed quantum mechanical understanding of the microscopic nature of these intraband excitations over a range of experimentally-relevant sizes and densities.

The classical Mie theory of plasmons in nanoparticles predicts a localized surface plasmon resonance when $\varepsilon_{1}(\Omega_{\mathrm{LSPR}})=-2\varepsilon_{\mathrm{m}}$, where $\varepsilon_{1}(\omega)$ is the real part of the complex bulk dielectric function for the nanoparticle, $\varepsilon(\omega)=\varepsilon_{1}(\omega)+i\varepsilon_{2}(\omega)$ and $\varepsilon_{\mathrm{m}}$ is the dielectric constant of the mediumWillets and Van Duyne (2007). The Drude plasmon pole approximation to the dielectric function, $\varepsilon(\omega)=\varepsilon_{\infty}-\Omega_{\mathrm{p}}^{2}/(\omega^{2}+i\gamma\omega)$, leads to the expression for the plasmon frequency $\Omega_{\mathrm{LSPR}}=\sqrt{\Omega_{\mathrm{p}}^{2}/(\varepsilon_{\infty}+2\varepsilon_{\mathrm{m}})-\gamma^{2}}$, where $\Omega_{\mathrm{p}}=\sqrt{4\pi\rho}$ is the bulk plasma frequency, $\rho$ is the free charge carrier density, $\varepsilon_{\infty}$ is the high-frequency dielectric constant, and $\gamma$ is the scattering rate of the electrons. For a given material and medium, this Mie+Drude plasmon frequency only depends on the density and not on the size, therefore failing to account for quantum confinement effects in small nanoparticles. The Drude dielectric function can be replaced with a microscopic dielectric function that accounts for quantum confinement effects, but this approach typically neglects interparticle interactionsKraus and Schatz (1983); Genzel, Martin, and Kreibig (1975); Scholl, Koh, and Dionne (2012); Jain (2014); Schimpf et al. (2014).

The prevailing quantum mechanical theory of plasmons in metals is the random-phase approximation (RPA) Bohm and Pines (1953); Gell-Mann and Brueckner (1957). As a theory of the ground-state energy density of bulk metals, the RPA famously removes the divergences encountered in finite-order perturbation theory. As a theory of the dynamical response, the RPA predicts the collective plasmon excitation, including its dispersion and strong oscillator strength Pines and Nozieres (1994); Giuliani and Vignale (2008). Despite its success for simple bulk metals, the RPA (by which we mean the time-dependent Hartree approximation) is not an accurate theory of excitation energies in molecules, casting doubt on its applicability to quantum plasmonics. In particular, the RPA fails to describe bound states such as excitons. To go beyond the RPA requires the tools of higher-level many-body theory or quantum chemistry.

In this paper, we investigate a model of interacting electrons confined to the interior of a sphere of radius $R$. This model generalizes the uniform electron gas (UEG), sometimes referred to as “jellium”, which is the canonical model of bulk metals and their plasmonic excitations. In the $R\rightarrow\infty$ limit, our model approaches the UEG (up to a background charge density). At finite $R$, quantum confinement produces a one-particle spectrum that is gapped, which alters the nature of the dominant excitations. We develop the requisite machinery, especially the two-electron integrals, and evaluate the performance of the RPA, time-dependent Hartree-Fock (TDHF), and configuration interaction with single excitations (CIS), which is the Tamm-Dancoff approximation (TDA) to TDHF. For accessible system sizes, we compare our results to higher-level equation-of-motion coupled-cluster theory with single and double excitations (EOM-CCSD) Emrich (1981); Koch and Jorgensen (1990); Bishop (1991); Stanton and Bartlett (1993); Bartlett and Musiał (2007); Krylov (2008), which has recently been used to characterize plasmons in the bulk UEG Lewis and Berkelbach (2019).

The intraband excitation in our model evolves from confinement-dominated and single-particle in character, to excitonic, to plasmonic, when the number of electrons is increased. By analyzing the underlying physics and comparing to higher-level EOM-CCSD, we argue that the evolution of the excitation is most accurately described by TDHF; the RPA and CIS are distinct approximations to TDHF and are only able to correctly describe plasmons and excitons, respectively. The transition from a confinement-dominated to plasmon-like excitation can be driven by both the number and density of electrons, and we find that for a fixed number of electrons, increasing the density actually decreases the plasmonic character of the excitation, opposite to the prediction of noninteracting models Jain (2014), due to increasing quantum confinement. The combination of local density functional theory (DFT) and the RPA, a popular technique in literatureProdan and Nordlander (2003); Zhang et al. (2014, 2017); Ipatov, Gerchikov, and Guet (2018), predicts results that are similar to TDHF, but due to a cancellation of errors. The character of the excitation across the entire range of number of electrons and density, in particular the intermediate excitonic state, can only be described by properly accounting for the attractive electron-hole (exchange) interaction, similar to the situation in molecules or semiconductors.

II Theory

II.1 Model

As a model of a doped nanoparticle, we treat the conduction band electrons as particles in an infinitely deep spherical well, where the atomic details of the nanoparticle are represented by the effective mass, radius, and dielectric constant. When used, the dielectric constant approximates the effect of the ignored valence electrons and higher excitations, which screen the Coulomb interaction. A more sophisticated model of the surface would include a finite or stepped barrier, however we do not expect our qualitative conclusions to be sensitive to the details of the surface. Furthermore, we neglect the dielectric contrast with the environment, which alters both the single particle band gap and the optical gap, although these two effects partially cancel in low-dimensional semiconductors Cho and Berkelbach (2018). We will treat the two-body Coulomb interactions between the conduction band electrons, which represents the focus of this work and goes beyond simple models of noninteracting electrons under confinement. We note that this “jellium sphere” model and various levels of theory have also been used to describe the structure and excitations of nuclei (including the giant dipole resonance Ring and Schuck (2004)) and especially the optical properties of alkali clusters Puska, Nieminen, and Manninen (1985); Guet and Johnson (1992); de Heer (1993); Brack (1993); Koskinen, Manninen, and Lipas (1994); Madjet, Guet, and Johnson (1995); Patterson (2019).

In first quantization and atomic units, the total $N$-electron Hamiltonian is

[TABLE]

where $m^{*}$ is the effective mass of the conduction band electrons and $v(r<R)=0$ and $v(r\geq R)=\infty$. Although charge neutrality would imply an additional $r$-dependent harmonic potential, here we neglect this potential because many experimental procedures for nanoparticle doping (e.g. photodoping with hole scavengers Schimpf et al. (2014)) do not preserve charge neutrality.

We use an orthogonal one-particle basis of eigenfunctions of the one-electron part of the above Hamilonian, corresponding to the well-known particle-in-a-sphere (PIS),

[TABLE]

where $r$ is the radial coordinate, $\alpha_{nl}=k_{nl}/R$, $k_{nl}$ is the $n$th zero of the spherical Bessel function $j_{l}(r)$, and the normalization constant is $N_{nl}=\sqrt{R^{3}/2}\left|j_{l+1}(k_{nl})\right|$. Each orbital is characterized by three quantum numbers $n,l,m$, with the limits $n\geq 1$, $l\geq 0$, and $m=-l,...,l$. The noninteracting orbital energies are given by

[TABLE]

and are $(2l+1)$-fold degenerate. Unlike the hydrogen atom, the orbital energies of the PIS are not degenerate with respect to the principle quantum number $n$, but are $m$-fold degenerate for a given $n$ and $l$.

In this orthogonal single-particle basis, the interacting, second-quantized Hamiltonian is

[TABLE]

where the indices $pqrs$ run over spin-orbitals, i.e. $p=(n,l,m,\sigma)$ and the two-electron integrals are given by

[TABLE]

where ${\bm{x}}=({\bm{r}},\sigma)$ is a combined space and spin variable.

We study nanoparticles containing $N=2$, 8, 18, 32, 50, 72, and 98 electrons, which correspond to closed-shell solutions of restrictedf Hartree-Fock (RHF). For these system sizes, we find that the RHF solution occupies the orbitals 1s,1p,1d,… up to 1$l_{\mathrm{max}}$, and therefore these closed-shell fillings correspond to $N=2\sum_{l=0}^{l_{\mathrm{max}}}(2l+1)=2(l_{\mathrm{max}}+1)^{2}$, where the factor of 2 accounts for spin.

The naive way to grow the basis set is to add PIS orbitals based on increasing energy; however, the RHF orbitals are pure eigenfunctions of $l$, so increasing the number of basis functions $n$ for each $l=0...l_{\mathrm{max}}$ is sufficient to converge the ground state calculation. In order to capture singly excited states we add an additional shell $l_{\mathrm{max}}+1$ based on the dipole selection rule $\Delta l=\pm 1$. The rapidly increasing degeneracy of the basis functions limits the number of electrons to 98, which we converge with 483 basis functions ($n_{\mathrm{max}}$ of $[10,9,9,8,8,7,7,7]$ for $l=0...l_{\mathrm{max}}+1=0...7$) requiring about 50 GB to store the two-electron integrals.

Expressions for the two-electron integrals $\langle pq|rs\rangle$, which are not analytic but can be reduced to two-dimensional quadrature along the radial axis, are given in the Appendix. All electronic structure calculations are performed by defining a custom Hamiltonian for use in the PySCF software package Sun et al. (2018).

II.2 Excited states

We focus on quantum chemical single-excitation theories due to their favorable $O(N^{4})$ scaling with system size, which makes them practical for future atomistic studies. Specifically, we consider excited states of the form

[TABLE]

where here and throughout $i,j,k,l$ and $a,b,c,d$ index occupied and unoccupied HF orbitals, and $X_{ai}$ and $Y_{ai}$ correspond to coefficients for the excitation and deexcitation of an electron from orbital $i$ to $a$, respectively. The deexcitation operator implies that the ground state $|\Psi_{0}\rangle$ is potentially correlated, though unspecified. The amplitudes $X_{ai}$ and $Y_{ai}$ are obtained from the eigenvalue problem Ring and Schuck (2004); Scuseria, Henderson, and Sorensen (2008)

[TABLE]

where $\mathbf{\Omega}$ is a diagonal matrix of excitation energies. The single excitation theories considered here correspond to specific choices of the $\mathbf{A}$ and $\mathbf{B}$ matrices. The most “complete” theory is TDHF, for which

[TABLE]

and the antisymmetrized integrals are defined as $\langle pq||rs\rangle=\langle pq|rs\rangle-\langle pq|sr\rangle$. The RPA is obtained by neglecting this antisymmetrization (consistent with time-dependent Hartree theory), and retains only the “direct” Coulomb matrix elements. The neglected “exchange” Coulomb matrix elements (also sometimes called “direct electron-hole interactions”) are responsible for the electron-hole attraction and the formation of bound excitons. The Tamm-Dancoff approximation (TDA) corresponds to $\mathbf{B}=0$, which neglects potential correlation in the ground state. When applied to TDHF and the RPA, the TDA leads to theories we will refer to as CIS and the RPA(TDA).

The RPA is the minimal theory necessary for the description of plasmons. For the UEG, the RPA predicts a plasmon dispersion $\Omega(q)$ that has the correct long-range limit, $\Omega(q\rightarrow 0)=\Omega_{\mathrm{p}}$, where $\Omega_{\mathrm{p}}=\sqrt{4\pi\rho}$ is the classical plasma frequency Pines and Nozieres (1994); Giuliani and Vignale (2008). The TDA (including CIS) predicts a collective excitation, but one whose energy unphysically diverges in the long-range limit, $\Omega^{\mathrm{TDA}}(q\rightarrow 0)\rightarrow\infty$. Here, we will see the same behavior in the $R\rightarrow\infty$ limit, highlighting the failure of the TDA for large plasmonic nanoparticles.

For comparison, we also present results obtained at lower and higher levels of theory, where available. At a low level, we consider simply the orbital energy differences from the noninteracting (PIS) and mean-field (HF) theories, i.e. the difference in energies of the occupied and unoccupied orbitals with the largest transition dipole matrix element. For HF, these orbitals are always the highest occupied and lowest unoccupied molecular orbitals (HOMO and LUMO). At a high level, we present results from equation-of-motion coupled-cluster theory with single and double excitations (EOM-CCSD) up to 50 electrons. Importantly, we note that EOM-CCSD includes ground-state correlation and excited-state double excitations, which both contribute to screening in an effective single-excitation theory like the Bethe-Salpeter equation. We use a basis set that adds PIS orbitals based on increasing energy to converge the EOM-CCSD calculations.

Since we are interested in the absorption properties of doped nanoparticles, we calculate the dynamical polarizability

[TABLE]

where $\hat{\mu}=\sum_{n=1}^{N}\hat{r}_{n}=\sum_{pq}r_{pq}a_{p}^{\dagger}a_{q}$ is the dipole operator.

III Results and Discussion

We study a model nanoparticle with $m^{*}=0.28$; which is characteristic of the conduction band of ZnO, whose plasmonic properties under doping have been experimentally studied Schimpf et al. (2014). As discussed in Sec. II.1, we consider systems containing 2, 8, 18, 32, 50, 72, and 98 electrons and study experimentally relevant densities from $\rho=1.4\times 10^{20}$ to $1\times 10^{22}$ cm*-3*, which defines the radius $R$ for a given number of electrons $N$. In this range, the radii of the nanoparticles are on the order of 1-10 nm.

III.1 Spectral properties and peak position

First we analyze the single-particle orbital energies. The top two panels of Figure 1 show the the density of states (DOS) predicted by the noninteracting (PIS, left) and mean-field (HF, right) level of theory, where blue lines indicate occupied orbitals and red lines indicate unoccupied orbitals. Results are shown for $N=2,8,32,$ and 98 electrons (bottom to top) at the smallest and largest densities considered here, $1.4\times 10^{20}$ cm*-3* and $1\times 10^{22}$ cm*-3* (left and right), the former of which corresponds to the ZnO system. In both cases, the energy spacings are reduced with increasing $R$ due to the scaling of the kinetic energy term, $1/(2m^{*}R^{2})$. In the nonineracting results, the energy of the HOMO (roughly the Fermi energy) decreases with increasing $R$, which is typical particle-in-a-box type physics. However in the HF results, the spectrum is further compressed and the energy of the HOMO is shifted to higher energies due to the mean-field effects of the Coulomb interaction.

In the middle two panels of Figure 1, we show the RHF (black) and TDHF absorption spectrum (red), and the density of excited states $D(\omega)=\sum_{m}\delta(\omega-\Omega_{m})$ (grey, filled) for the same systems as above. From the spectrum of excited states, which is becoming gapless in the limit of large $N$, we observe a dominant peak in the absorption spectrum that arises from a state that is typically not the lowest in energy. This redistribution of oscillator strength from a low-energy continuum into a single high-energy state is reminiscent of the plasmon peak in the dynamical structure factor of the UEG Ring and Schuck (2004); Lewis and Berkelbach (2019); the bright state is only the lowest in energy for $N=2$. Henceforth, we focus on the excitation energy of the dominant bright peak, which is the energy of the excited state with the largest transition dipole matrix element. The spectra in Figure 1 also highlight a computational challenge, as the strongly absorbing state of interest is buried in the interior of the eigenvalue spectrum of the Hamiltonian.

In the bottom two panels of Figure 1, we compare the energy of this dominant absorption peak versus radius, for these same two densities, as predicted by various theories. We present results for TDHF, the RPA, CIS, RPA(TDA), and EOM-CCSD, along with the noninteracting and HF transitions for comparison. The EOM-CCSD result, the most accurate solution here, demonstrates the qualitative evolution of the excitation, which we separate into three regions. At small $R$, the excitation energy is dominated by the kinetic (confinement) energy, and scales with the band gap. At intermediate $R$, the EOM-CCSD result is below the HF gap and exhibits a minimum, which we attribute to the formation of intraband excitons. At large $R$, the excitation energy increases with increasing $R$ and goes above the HF gap before reaching a plateau near the classical plasma energy. This latter behavior is consistent with the $R\rightarrow\infty$ limit of our model, i.e. the UEG, which has a well-known $q\rightarrow 0$ plasmon at the classical plasma frequency.

We observe from Figure 1 that TDHF is the most accurate single excitation theory. For a given density, as the number of electrons (or radius) increases, the energy of the TDHF excitation crosses from below to above the HF gap, which follows the EOM-CCSD result. This behavior is consistent with the physics embodied in TDHF, which contains the ingredients necessary to capture the three regimes described above, i.e. confinement-dominated excitations, excitons (due to exchange integrals, i.e. direct electron-hole interactions), and plasmons (due to the nonzero $\mathbf{B}$ matrix). These claims are validated by comparison with the “approximations” to TDHF, which contain only a subset of these ingredients: the RPA correctly predicts the evolution of the excitation to the classical plasmon (always above the HF gap), but cannot lower the energy of the excitaton at any radius due to the lack of electron-hole interactions, and CIS contains the exchange Coulomb interaction but neglects the $\mathbf{B}$ matrix, so it is accurate in describing bound excitons at small radii, but fails to correctly describe the plasmonic state at large radii. In particular, CIS predicts an excitation energy which goes above the HF gap and diverges at large radius, which is precisely analogous to its $q\rightarrow 0$ behavior in the UEG. Finally, the RPA(TDA) result is not accurate at any radius.

Despite the good overall agreement of TDHF with EOM-CCSD at high electron densities, the TDHF solution suffers at low density and large $R$ (always underpredicting the EOM-CCSD result), which is precisely where electron correlation is expected to be strongest. In this regime, one may anticipate the onset of Wigner crystallization and multireference character Thompson and Alavi (2002); Jung and Alvarellos (2003); Thompson and Alavi (2004, 2005). While this physics could be approximately addressed via spin-symmetry breaking, we do not pursue this approach here.

III.2 Characterization of the excited state

We now discuss a microscopic characterization of the wavefunction of the strongly absorbing excited state. The top two plots of Figure 2 summarize the energy of the peak absorption for TDHF and the RPA, and the bottom two plots of Figure 2 plot the sum of the $Y$ coefficients for the peak absorption, which is a measure of the plasmonic character of the excitation. The plasmonic character is sensitive to both the absolute number of electrons and the density, and the emergence of a plasmon is not solely a density-driven transition. For a fixed number of electrons, the plasmonic character decreases with increasing density, opposite to the noninteracting result Jain (2014); at higher densities, a greater number of electrons is needed to reach the same plasmonic character as a fewer number of electrons in a lower density system. This decrease in plasmonic character at smaller $R$ is due to the quantum confinement of the one-electron kinetic energy that scales as $1/(2m^{*}R^{2})$, and further exacerbated by the strong electron-hole interaction (compare TDHF to the RPA). When the energy of the peak absorption is compared to the plasmonic character, we find the surprising result that even though the energy lies within the range of classical plasma frequencies, the excitation can be far from plasmonic, and is instead excitonic or single-particle-like.

Figure 3 is a series of plots of the induced charge density,

[TABLE]

where $\chi$ is the density-density linear response function, for an external electric field oriented along the $z$-axis. We evaluate the induced charge density at the peak absorption energy for the noninteracting, HF, TDHF, and CIS levels of theory. The doping density of $1.4\text{\times}{10}^{20}\text{\,}\mathrm{c}\mathrm{m}^{-3}$ and 2 (top), 8 (middle), and 98 (bottom) electrons correspond to nanoparticle radii of 1.5 nm, 2.4 nm, and 5.5 nm, and excitation character of confined, excitonic, and plasmonic. As the number of electrons increases, the induced density concentrates at the surface, in agreement with classical theory. The addition of mean-field interactions (HF) repels the non-interacting induced density to the surface. The induced density at the mean-field level is not changed under configuration mixing, because the HOMO-LUMO transition has the largest dipole matrix element $\langle i|\hat{z}|a\rangle$ and this configuration contributes with a large weight in the bright state. In other words, the excited-state wavefunctions are all qualitatively similar, despite predicting very different energies.

In light of the above observation, we next quantify the contributions of the various singles theories to the total excitation energy,

[TABLE]

The second line in Eq. (11), containing the “direct” $\langle ib|aj\rangle$ integrals that describe the interaction energy between single-particle excitations, has been termed the “plasmonicity” Zhang et al. (2017) and is nonzero for TDHF, RPA, CIS, and RPA(TDA). The third line, containing the “exchange” $\langle ib|ja\rangle$ integrals from the $\mathbf{A}$ matrix, provides a measure of the exciton binding energy and is nonzero for TDHF and CIS. Finally, the fourth line, containing the “exchange” $\langle ij|ba\rangle$ integrals from the $\mathbf{B}$ matrix has no classical interpretation and is only nonzero for TDHF elements.

The bottom of Figure 3 plots the above contributions versus radius for a doping density of $1.4\text{\times}{10}^{20}\text{\,}\mathrm{c}\mathrm{m}^{-3}$, for all four singles theories. CIS and the RPA(TDA) do not account for ground state correlations ($\mathbf{B}=0$) and predict a diverging plasmonicity. The RPA adds a nearly constant plasmonicity at all $R$, such that the excitation energy within the RPA is only modulated by the mean-field gap, which is consistent with the prediction of a higher energy collective excitation that is a combination of degenerate single-particle excitationsRing and Schuck (2004). The RPA and RPA(TDA) overpredict the excitation energy at small radii due to the lack of the excitonic interaction (see Figure 2, bottom). CIS fails at large $R$, but when quantum confinement dominates at small $R$, both CIS and TDHF correctly describe the bound excitonic state: notably, for the most confined case of two electrons, the exciton binding energy lowers the excitation energy by roughly 1.5 eV. The exciton binding energy decreases with increasing $R$ as the electron-hole spatial overlap becomes smaller. The nonclassical exchange contribution from the $\mathbf{B}$ matrix makes a small positive contribution at all values of $R$, never exceeding 0.2 eV.

TDHF mixes the correct large $R$ limit of the RPA and the small $R$ limit of CIS. Importantly, TDHF inherits the RPA prediction of a relatively constant plasmonicity at all radii, Figure 3 bottom left. Thus, at least at this density, there is no clearly distinguishable single-particle excited state, and the exciton binding energy and plasmonicity make non-negligible contributions at all values of $R$.

III.3 Schematic model

We can understand the dependence of the energy of the excitation on radius by considering a schematic model of the RPA and RPA(TDA) theories, i.e. those without antisymmetrized integrals. This treatment is motivated by the presentation in Ref. Ring and Schuck, 2004. The results of the schematic model at all levels of theory will be compared to the numerical results at four densities in Figure 4.

We consider the factorization of the direct two-electron integrals $\langle ib|aj\rangle\approx\lambda\rho_{ai}\rho_{bj}$, with $\lambda>0$. In this approximation, the RPA(TDA) equation for an excited state $n$ can be written as

[TABLE]

which leads to Ring and Schuck (2004)

[TABLE]

The latter equation can be solved graphically for $\Omega_{n}$, leading to a number of single-particle excitations with energies approximately given by $\varepsilon_{a}-\varepsilon_{i}$ and one higher-energy collective excitation (the plasmon). For illustrative purposes, we consider the subspace containing only the HOMO and LUMO (each potentially degenerate), such that $\varepsilon_{a}-\varepsilon_{i}=\varepsilon$, which yields for the plasmon state

[TABLE]

where $N_{\mathrm{trans}}$ is the number of transitions and $\bar{\rho}$ is an average quantity.

We now seek to understand the $R$ dependence of this excitation energy at fixed density. For all of the system sizes studied here, the closed-shell RHF solution is only stable for those configurations for which all occupied orbitals have $n=1$, i.e. 1s2, 1p6, 1d10, 1f14, …, and so on. The angular momentum of the HOMO is thus defined by the number of electrons, $N=2(l_{\mathrm{max}}+1)^{2}$ or $l_{\mathrm{max}}=\sqrt{N/2}-1$ at fixed density. To a good approximation, we find that the zeros of the spherical Bessel functions can be written as $k_{nl}=n\pi+1.32l$ (in particular, the value 1.32 is empirically better than the asymptotic value $\pi/2$). This yields a noninteracting band gap from $1l$ to $1(l+1)$ of

[TABLE]

At the HF level, we find that the form of the one-electron contribution to the band gap is very similar and most significantly modified by the exchange contribution, which we model with the form $1/R$ for all densities, $\varepsilon_{\mathrm{g}}^{\mathrm{HF}}(R)=\varepsilon_{\mathrm{g}}^{\mathrm{NI}}(R)+1/R$.

The HOMO and LUMO each have a degeneracy proportional to $l(l+1)$ which, combined with the dipole selection rule $\Delta m=0,\pm 1$, leads to a number of transitions $N_{\mathrm{trans}}\propto l_{\mathrm{max}}\propto\sqrt{\rho R^{3}}$. The Coulomb interaction has a scaling $\lambda(R)\propto R^{-1}$. Therefore, within the RPA(TDA), the schematic model predicts an excitation energy

[TABLE]

where we take $c=1$. At small $R$, the excitation energy of a doped nanoparticle is given by the kinetic-energy-determined band gap; at large $R$, the excitation energy diverges due to the Coulomb interaction. This divergence is unphysical and analogous to the behavior of the TDA in the $q\rightarrow 0$ limit of the three-dimensional uniform electron gas.

The divergence is fixed in the full RPA, which is given in the schematic model by

[TABLE]

The same approximations as above for a spherical nanoparticle leads to

[TABLE]

This excitation energy has the same kinetic-energy-determined band gap at small $R$, but now has a finite $R\rightarrow\infty$ limit, $\Omega^{\mathrm{RPA}}(R\rightarrow\infty)=\sqrt{5.04\rho/m^{*}}$. Importantly, the schematic model recovers the exact limiting form of the classical plasmon energy, up to constants of order 1, i.e. $\sqrt{5.04}\approx 2.24$ compared to $\sqrt{4\pi}\approx 3.54$.

Before continuing on to theories with exchange (CIS and TDHF), we first turn to an analysis of the excitation coefficients. Again within the degenerate schematic model, the $X$ coefficients in RPA(TDA) are

[TABLE]

where $C=\sum_{ai}|\rho_{ai}|^{2}$ is a normalization constant, and in the RPA are

[TABLE]

where $\Omega$ is the energy of the collective state and

[TABLE]

We note that as $\varepsilon_{\mathrm{g}}\rightarrow 0$ and the model becomes more metallic, the excited state is a plasmon and the $X$ and $Y$ coefficients become equal in magnitude. In contrast, as the gap $\varepsilon_{\mathrm{g}}$ increases, the $X$ coefficients dominate and $Y_{ai}\rightarrow 0$. This behavior is observed numerically, as shown in Figure 2. The sum of deexcitation coefficients for the collective state, $\sum_{ai}|Y_{ai}|^{2}$, is therefore another measure of plasmonic character. From Figure 2, it is clear to see that this character increases with increasing $R$ at fixed $N$ or increasing $N$ at fixed $R$, because both situations correspond to reducing the gap $\varepsilon_{\mathrm{g}}$. By this measure, a low plasmonic character does not imply that the excitation is single-particle-like, because the excitation can still be delocalized over the $X_{ai}$ coefficients.

The inclusion of antisymmetrized integrals in CIS and TDHF prevents an analytic treatment of the schematic model because the integral factorization $\langle ia|jb\rangle\approx\lambda\rho_{ij}\rho_{ab}$ does not facilitate the solution of the eigenvalue problem. However, the largest-in-magnitude element $\langle ia|ia\rangle$ can be included exactly as it just shifts the HF band gap $\varepsilon_{\mathrm{g}}\rightarrow\varepsilon_{\mathrm{g}}-\langle ia|ia\rangle$. Within the schematic model, we take all such excitonic Coulomb integrals to be equal and obeying the scaling $\langle ia|ia\rangle=c^{\prime}/R$, with $c^{\prime}=1.4$. This approximation gives the CIS and TDHF excitation energies as

[TABLE]

[TABLE]

We note that the remaining excitonic Coulomb integrals $\langle ia|jb\rangle$ could be included via perturbation theory, though we do not pursue this here.

In Figure 4, we show the performance of the analytic solutions of this schematic model compared to our numerical calculations. The noninteracting, HF, RPA, and RPA(TDA) schematic models fit the data remarkably well, with the only parameterizations being the fit of the spherical Bessel function zeros and the proposed $1/R$ form of the HF exchange contribution to the gap. The CIS and TDHF schematic models fit the data reasonably well, with the main difficulty being the estimation of the contribution of the antisymmetrized two-electron integrals. Importantly, the schematic model captures the overall behaviour of our calculations based upon simple scaling arguments on the density and radius, and supports the physical interpretations given to the various single-excitation theories.

III.4 Effect of screening and comparison to DFT

In our model, the effect of screening from the ignored valence electrons can be approximately accounted for by a scaling of the Coulomb interaction, $r_{12}^{-1}\rightarrow(\varepsilon r_{12})^{-1}$. We can approximate ZnO nanoparticles by setting $\varepsilon=3.72$. Figure 5 plots the energy of the intraband absorption as a function of radius for this “ZnO” nanoparticle, at the experimental doping density of $1.4\text{\times}{10}^{20}\text{\,}\mathrm{c}\mathrm{m}^{-3}$, against the experimental results from Ref. Schimpf et al., 2014. We compare TDHF (left) with DFT plus the RPA (right), where the RPA excited state is calculated from the solution of the restricted Kohn-Sham equations with the local density approximation (LDA). The LDA+RPA approach has been applied previously to our model Varas et al. (2016); Zhang et al. (2017) and has found success in describing the experimental ZnO results Goings et al. (2014); Ipatov, Gerchikov, and Guet (2018); Gerchikov, Ipatov, and Guet (2018).

From Figure 5, both TDHF and DFT+RPA compare favorably with the experimental result, although they embody different physical effects as explored in this paper. As usual, the LDA gap is smaller than the RHF gap, due to the latter’s treatment of exchange. The inset on the right of Figure 5 plots the sum of all contributions given in Eq. (11) except for the orbital-energy contribution, from LDA+RPA and RHF+TDHF. The RPA only acts to increase the small LDA gap at all $R$; by contrast, the TDHF lowers the large RHF gap at small $R$ (due to excitonic effects) and increases the RHF gap at large $R$ (due to plasmonic effects). Therefore, ignoring exchange in the ground state and excited state, as in the LDA+RPA approach, can produce an accurate result due to cancellation of errorsJung et al. (2004).

IV Conclusions

In summary, we provide a fully quantum mechanical study of a confined, interacting electron gas as a model for doped semiconductor nanoparticles. We observe strongly absorbing excited states whose wavefunction character can be classified as single-particle-like (confinement dominated), excitonic, or plasmonic. Within the framework of the most computationally affordable single-excitation theories, only TDHF is capable of capturing the qualitative behavior at all studied densities and particle sizes. We also present a schematic model of the strongly absorbing excited state that reproduces the $R$-dependence observed in our simulations.

Our model is simple in order to focus on the essential features of electronic interactions in the excited states of confined systems. The model neglects atomistic details as well as surface, ligand, or solvent effects. The model is also ignorant of the doping mechanism and neglects the atomic defect potential that is introduced by impurity doping (but not by electron transfer or photodoping). Nonetheless, the results of our calculations argue strongly against the interpretation or simulation of doped nanoparticle spectra based on single-particle transitions between orbitals, and we propose an interpretation of intraband excitons as the primary excitations at low doping or small nanoparticles.

Looking forward towards atomistic or tight-binding Pi and Delerue (2013) simulations, our work has two important ramifications. First, the TDA fails spectacularly and should be avoided in all simulations seeking to address the possibility of plasmonic excitations. Second, the retention of attractive electron-hole “exchange” integrals is essential for an accurate wavefunction description of excitonic states. With these criteria in mind, we suggest that the most promising and affordable ab initio methods are TDHF (as explored here), TDDFT with hybrid functionals, or the $GW$+Bethe-Salpeter equation approach without the TDA Sander, Maggio, and Kresse (2015).

Acknowledgments

T.C.B. thanks Philippe Guyot-Sionnest for early conversations related to this work. All calculations were performed with the PySCF software package Sun et al. (2018), using resources provided by the University of Chicago Research Computing Center. This work was supported by the Air Force Office of Scientific Research under AFOSR Award No. FA9550-18-1-0058 and by the National Science Foundation CAREER program under Award No. CHE-1848369. The Flatiron Institute is a division of the Simons Foundation.

Appendix A Two-electron integrals for particle-in-a-sphere orbitals

The spherical harmonics $Y_{l}^{m}$ are generally complex. To maximize the symmetry of the two-electron integrals, we use the real form of the spherical harmonics, $y_{l\mu}=\sum_{m}U_{m\mu}^{l}Y_{l}^{m}$. With this choice, the two-electron integrals $\langle pq|rs\rangle$ are given by

[TABLE]

where the angular integrals are

[TABLE]

a linear combination of the integrals of three complex spherical harmonics

[TABLE]

which vanishes unless $l+l_{1}+l_{2}=2g$, $g\in\mathbb{Z}$, and $m_{1}+m_{2}=m$, thus truncating the infinite sum over $l$. The integral of three real spherical harmonics is invariant under all permutations of the order of the functions and can be simplified into a single complex integral times the appropriate factorsHomeier and Steinborn (1996). The radial integral of the normalized spherical Bessel functions $R_{nl}$ is

[TABLE]

and can be computed numerically.