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. Modeling the excess electrons as interacting electrons confined to a sphere, we find that the excitation evolves from single-particle 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 and electron–hole attraction. In addition, the excitation passes through an intermediate regime where it is best characterized as an intraband exciton. We compare equation-of-motion coupled-cluster theory with those of more affordable single-excitation theories and identify the inclusion of electron–hole interactions as essential to describing the evolution of the excitation. Despite the simplicity of our model, the results are in reasonable agreement with the experimental spectra of doped ZnO nanoparticles at a doping density of 1.4 × 10^{20} cm^{−3}. Based on our quantum chemistry calculations, we develop a schematic model that captures the dependence of the excitation energy on nanoparticle radius and electron density.

## I. INTRODUCTION

Metallic nanoparticles are an important optoelectronic platform because of their strong surface plasmon resonance at visible energies, which can be tuned by the size, shape, and environment;^{1,2} however, the accessible carrier densities are limited to those of the parent metals and are typically 10^{22} cm^{−3} or higher. Recently, doping of semiconductor nanoparticles has enabled access to much lower electron densities.^{3–9} 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.^{10–13}

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.”^{14–22} Furthermore, the collective plasmon picture becomes dubious for systems containing only 1–10 excess electrons, suggesting a transition to “single-particle” excitations^{9} 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. Specifically, the focus of this work is to understand how the optical properties of doped nanocrystals are influenced by electron correlation, which is usually completely neglected or treated approximately in an uncontrolled manner.

The classical Mie theory predicts a localized surface plasmon resonance (LSPR) when *ε*(Ω_{LSPR}) = −2*ε*_{m}, where *ε*(*ω*) is the complex bulk dielectric function for the nanoparticle and *ε*_{m} is the dielectric constant of the medium.^{1,23} If we assume Im*ε* ≈ 0, then the Drude plasmon pole approximation to the dielectric function, $\epsilon (\omega )=\epsilon \u221e\u2212\Omega p2/(\omega 2+i\gamma \omega )$, leads to the expression for the plasmon frequency $\Omega LSPR=\Omega p2/(\epsilon \u221e+2\epsilon m)$, where $\Omega p=4\pi \rho /m*$ is the bulk plasma frequency, *ρ* is the free charge carrier density, *m*^{*} is the effective mass, *ε*_{∞} is the high-frequency dielectric constant, and we have taken the limit *γ* = 0. In the simple case of *ε*_{∞} = *ε*_{m} = 1, the LSPR occurs at $\Omega LSPR=\Omega p/3$. For a given material and medium, this classical LSPR 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 still neglects interparticle interactions.^{5,7,24–26}

The prevailing interacting quantum mechanical theory of plasmons in metals is the random-phase approximation (RPA).^{27,28} 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 bulk plasmon excitation, including its dispersion and strong oscillator strength.^{29,30} Despite its success for simple metals, the RPA (also known as 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 apply equation-of-motion coupled-cluster theory with single and double excitations (EOM-CCSD),^{31–36} which contains the physical ingredients needed for an accurate treatment of excitons and plasmons.^{36–38} Based on this high-level theory, we 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. We apply these quantum chemistry techniques to interacting electrons confined in an infinitely deep spherical well. 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* → ∞ 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.

By analyzing the underlying physics and comparing to 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 use of time-dependent density functional theory (TDDFT) with the local density approximation^{22,39–41} predicts results that are similar to TDHF, but due to the cancellation of errors. The intraband excitation in our model evolves from confinement-dominated and single-particle in character to excitonic and to plasmonic, when the number of electrons is increased. The transition from a confinement-dominated to plasmon-like excitation can be driven by both the number and density of electrons; we find that for a fixed number of electrons, increasing the density decreases the plasmonic character of the excitation due to increasing quantum confinement and electron–hole attractions. While this is the same qualitative trend observed in noninteracting theories,^{26} the latter conclusions must be considered carefully because noninteracting theories alone do not predict plasmons. The exciton binding in small doped nanocrystals will impact their applications in charge transfer, sensing, and catalysis. The character of the excitation across the entire range of the 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

### A. 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. This neglect of atomistic detail is an approximation, but the model parameters are chosen to ensure that the energy scales are appropriate, leading to correct qualitative conclusions about the nature and behavior of collective excitations in small assemblies of confined, interacting electrons. A direct comparison to experiment is precarious, but the results in Sec. III D suggest that our model captures the essential physics. This “jellium sphere” model and various levels of theory have been used to successfully describe the structure and excitations of nuclei (including the giant dipole resonance, a collective harmonic mode^{42}) and especially the optical properties of alkali clusters.^{43–49}

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

where *m*^{*} is the effective mass of the conduction band electrons, *ε*_{∞} is the dielectric constant of the nanoparticle, and *v*(*r* < *R*) = 0 and *v*(*r* ≥ *R*) = ∞. 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^{6,7}) do *not* preserve charge neutrality.

Throughout this manuscript, we will use a Coulomb interaction in our Hamiltonian equation (1) that is unscreened, i.e., *ε*_{∞} = 1. In a first-principles quantum chemistry approach, screening arises naturally from higher-order excitations and not from a separate calculation (such as in the GW/Bethe–Salpeter equation approach^{50}), and the comparison between single-excitation theories and EOM-CCSD provides some insight into the role of screening (see below). However, the neglected valence band electrons may also contribute to screening, which we can account for with a static dielectric constant because the semiconductor bandgap is relatively large. This is only done in Sec. III D, when we compare to experimental results on ZnO nanoparticles. We note that the simple use of a screened 1/*r* interaction neglects the dielectric contrast between the nanoparticle and its environment.

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

Further computational details are provided in the Appendix. In particular, we describe the choice of basis sets and provide expressions for the two-electron integrals ⟨*pq*|*rs*⟩, which are not analytic but can be reduced to a two-dimensional quadrature along the radial axis. All electronic structure calculations are performed by defining a custom Hamiltonian for use in the PySCF software package.^{51}

### B. Excited states

Due to the potential importance of electron correlation, we focus on quantum chemical theories of excited states. The highest level of theory we apply is equation-of-motion coupled-cluster theory with single and double excitations (EOM-CCSD).^{31–36} To the best of our knowledge, this represents the first application of EOM-CCSD to the study of plasmons in finite systems. Given a set of one- and two-electron integrals, the formalism is standard and we refer the interested reader to the literature.^{36} However, in terms of its physical content, we note that EOM-CCSD includes ground-state correlation and excited-state double excitations, which both contribute to screening and lifetime effects in an effective single-excitation theory such as the Bethe–Salpeter equation. Furthermore, EOM-CCSD completely subsumes all single-excitation theories considered below. Within EOM-CCSD, we find the plasmon state by using an iterative eigensolver that targets states having a large overlap with the CIS plasmon state.

The high cost of EOM-CCSD, which scales as *N*^{6} where *N* is a measure of the system size, will preclude its future applications to large, atomistic nanocrystals. Therefore, we also consider more affordable theories that consider only single excitations (and potentially de-excitations). Specifically, we consider excited states of the form

where here and throughout *i*, *j*, *k*, *l* and *a*, *b*, *c*, *d* index the occupied and unoccupied HF orbitals, respectively, and *X*_{ai} and *Y*_{ai} correspond to coefficients for the excitation and de-excitation of an electron from orbital *i* to *a*, respectively. The de-excitation operator implies that the ground state |Ψ_{0}⟩ is potentially correlated, though unspecified. The amplitudes *X*_{ai} and *Y*_{ai} are obtained from the eigenvalue problem,^{42,52}

where **Ω** is a diagonal matrix of excitation energies. The single-excitation theories considered here correspond to the specific choices of the **A** and **B** matrices. The most “complete” theory is TDHF for which

and the antisymmetrized integrals are defined as ⟨*pq*$\u2225$*rs*⟩ = ⟨*pq*|*rs*⟩ − ⟨*pq*|*sr*⟩. 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 **B** = 0, i.e., **AX** = **XΩ**, 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 TDA provides a minor reduction in computational cost, but the accessible system sizes are similar for methods with or without the TDA. We note that the strongly absorbing state of interest in doped nanoparticles is rarely the lowest-energy excited state, and therefore, we fully diagonalize the Hamiltonian in Eq. (3).

The RPA is the minimal quantum mechanical theory necessary for the description of plasmons. For the UEG, the RPA predicts a plasmon dispersion Ω(*q*) that has the correct long-range limit, Ω(*q* → 0) = Ω_{p}, where $\Omega p=4\pi \rho /m*$ is the classical plasma frequency.^{29,30} The TDA (including CIS) predicts a collective excitation, but one whose energy unphysically diverges in the long-range limit, Ω^{TDA}(*q* → 0) → ∞. Here, we will see the same behavior in the *R* → ∞ limit, highlighting the failure of the TDA for large plasmonic nanoparticles.

For comparison, we also present the results obtained at lower levels of theory; in particular, for the primary resonance, we consider the orbital energy differences from the noninteracting 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 molecular orbitals (HOMOs) and lowest unoccupied molecular orbitals (LUMOs).

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

where $\mu ^x=\u2211n=1Nx^n=\u2211pqxpqap\u2020aq$ is the *x* component of the dipole operator. The polarizability of metallic nanoparticles was theoretically shown to be dominated by the surface plasmon resonance with a minor contribution from the bulk plasmon.^{53} This behavior is in qualitative agreement with classical electrodynamics, a local theory that only attributes intensity to the surface plasmon. Due to its dominance, even in the quantum spectrum, the surface plasmon is the state of interest throughout this manuscript. We note that momentum-resolved techniques such as energy loss spectroscopy, which enable the tunable excitation of both surface and bulk plasmons,^{5,54} can be straightforwardly simulated using the formalism described here.

## 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.^{7} As discussed in Sec. II A, we consider systems containing 2, 8, 18, 32, 50, 72, and 98 electrons and study densities from the experimentally measured *ρ* = 1.4 × 10^{20} cm^{−3} up to 1 × 10^{22} cm^{−3},^{4} which defines the radius *R* for a given number of electrons *N*. The highest density studied is comparable to that of metal nanoparticles, where quantum confinement effects were studied early on with noninteracting theories.^{25} In this sense, our work studies the evolution of collective excitations from high densities achievable in metal nanoparticles down through many of the low densities achievable in doped semiconductor nanoparticles. In this range, the radii of the nanoparticles are on the order of 1–10 nm, where the effective mass approximation for ZnO is reasonable.^{55} We set *ε*_{∞} = 1, which ignores screening from the valence electrons, but will return to this topic in Sec. III D.

### A. Spectral properties and peak position

First, we analyze the single-particle orbital energies. The top two panels of Fig. 1 show the density of states (DOS) predicted by the noninteracting (left within a column) and mean-field (HF, right within a column) level of theory, where blue lines indicate the occupied orbitals and red lines indicate the 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 × 10^{20} cm^{−3} and 1 × 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/(2*m*^{*}*R*^{2}). In the noninteracting results, the energy of the HOMO (roughly the Fermi energy) decreases with increasing *R*, which is typical particle-in-a-box type physics. 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 Fig. 1, we show the RHF gap (black), TDHF absorption spectrum (red), and density of excited states *D*(*ω*) = ∑_{m}*δ*(*ω* − Ω_{m}) (gray, 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;^{37,42} 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. In the metallic limit, this peak corresponds to the surface plasmon. An additional small intensity can be seen in the TDHF spectrum, which corresponds to the bulk plasmon.

In the bottom two panels of Fig. 1, we compare the energy of this dominant absorption peak vs radius, for these same two densities, as predicted by various theories (Fig. 4 shows additional densities for TDHF and RPA). We present results for TDHF, 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 bandgap. 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; the exciton binding energy, defined as the difference between the gap and the excitation energy, can be as large as 1 eV even at the lowest electron densities. At large *R*, the excitation energy increases with increasing *R* and goes *above* the HF gap. While the relative cost of EOM-CCSD precludes the study of very large *R*, we see that the peak position approaches an energy that is above the classical LSPR energy. This blueshift may be due to the double excitations, as observed for the bulk plasmon energy,^{37} the strong surface confinement and level quantization, as observed in TDDFT calculations on jellium spheres,^{53} or the absence of a neutralizing background for the *charged* nanoparticles studied here.

We observe that TDHF is the most accurate single-excitation theory, based on comparison to EOM-CCSD, which provides the highest-quality treatment of electron correlation for our model Hamiltonian. 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 **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 excitation at any radius due to the lack of electron–hole interactions, and CIS contains the exchange Coulomb interaction but neglects the **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 that goes above the HF gap and diverges at a large radius, which is precisely analogous to its *q* → 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, where electron correlation is weak, 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, which reduces the accuracy of the RHF ground-state solution.^{56–59} While this physics could be approximately addressed via spin-symmetry breaking, we do not pursue this approach here.

### B. 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 Fig. 2 summarize the energy of the peak absorption for TDHF and the RPA, and the bottom two plots of Fig. 2 show the sum of the *Y* coefficients for the peak absorption, for all densities studied here. The value ∑_{ai}|*Y*_{ai}|^{2} measures the correlations in the ground state: as the HOMO–LUMO gap closes with increasing *R*, the ground state can be improved by mixing in excitations with the nearly degenerate virtual states [Eq. (3)], which corresponds to the approach to a metallic state and a plasmon excitation. The plasmonic character is sensitive to both the absolute number of electrons and the density. At a fixed number of electrons, the plasmonic character increases with decreasing density (increasing radius). At fixed density, the plasmonic character increases with increasing electron number (increasing radius). The sensitivity of the plasmonic character to density is exacerbated by electron–hole interactions (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 is close to the classical LSPR energy, 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,

where *χ* 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 × 10^{20} cm^{−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 (identified by the previous discussion). 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 noninteracting 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 ⟨*i*|$z^$|*a*⟩ 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,

The second line in Eq. (7), containing the “direct” ⟨*ib*|*aj*⟩ integrals that describe the interaction energy between single-particle excitations, has been termed the “plasmonicity”^{22} and is nonzero for TDHF, RPA, CIS, and RPA(TDA). The third line, containing the “exchange” ⟨*ib*|*ja*⟩ integrals from the **A** matrix, provides a measure of the exciton binding energy and is nonzero for TDHF and CIS. Finally, the fourth line, containing the “exchange” ⟨*ij*|*ba*⟩ integrals from the **B** matrix, has no classical interpretation and is only nonzero for TDHF elements.

Figure 3 (bottom) shows the plots of the above contributions vs radius for a doping density of 1.4 × 10^{20} cm^{−3}, for all four singles theories. CIS and the RPA(TDA) do not account for ground-state correlations (**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 excitations.^{42} The RPA and RPA(TDA) overpredict the excitation energy at small radii due to the lack of the excitonic interaction (see Fig. 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 **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 (Fig. 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*.

### C. Schematic model

In Secs. III A and III B, we showed the existence of an excitonic state in between the transition of a confinement-dominated excitation to a surface plasmon. The existence of such an excitonic state is significant as it reacts differently to the dielectric environment compared to a LSPR. Here, we derive a schematic model that gives the dependence of the energy of the excitation on the radius and separates the excitonic state from the plasmon based on simple scaling arguments. This treatment is motivated by the presentation in Ref. 42. The results of the schematic model at all levels of theory except EOM-CCSD are compared to the numerical results at four densities in Fig. 4.

We will first consider the RPA and RPA(TDA) theories, i.e., those without antisymmetrized integrals. Consider the factorization of the direct two-electron integrals ⟨*ib*|*aj*⟩ ≈ *λρ*_{ai}*ρ*_{bj}, with *λ* > 0. This factorization can be seen as a crude density fitting but is also a simple tool to enable analytic insights into the behavior of single-excitation theories. In this approximation, the RPA(TDA) equation for an excited state *n* can be written as

which leads to^{42}

The latter equation can be solved graphically for Ω_{n}, leading to a number of single-particle excitations with energies approximately given by *ε*_{a} − *ε*_{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 *ε*_{a} − *ε*_{i} = *ε*, which yields for the plasmon state

where *N*_{trans} is the number of transitions and $\rho \xaf$ 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., 1s^{2}, 1p^{6}, 1d^{10}, 1f^{14}, …, and so on. The angular momentum of the HOMO is thus defined by the number of electrons, $N=2(lmax+1)2$ or $lmax=N/2\u22121$, at fixed density. To a good approximation, we find that the zeros of the spherical Bessel functions can be written as *k*_{nl} = *nπ* + 1.32*l* (in particular, the value 1.32 is empirically better than the asymptotic value *π*/2). This yields a noninteracting bandgap from 1*l* to 1(*l* + 1) of

At the HF level, we find that the form of the one-electron contribution to the bandgap is very similar and most significantly modified by the exchange contribution, which we m2odel with the form 1/*R* for all densities, $\epsilon gHF(R)=\epsilon gNI(R)+1/R$.

The HOMO and LUMO each has a degeneracy proportional to *l*(*l* + 1), which, combined with the dipole selection rule Δ*m* = 0, ±1, leads to a number of transitions $Ntrans\u221dlmax\u221d\rho R3$. The Coulomb interaction has a scaling *λ*(*R*) ∝ *R*^{−1}. Therefore, within the RPA(TDA), the schematic model predicts an excitation energy

where *c* is a free parameter; we find *c* = 1 provides a good fit to our numerical data. At small *R*, the excitation energy of a doped nanoparticle is given by the kinetic-energy-determined bandgap; 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* → 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

The same approximations as above for a spherical nanoparticle lead to

This excitation energy has the same kinetic-energy-determined bandgap at small *R* but now has a finite *R* → ∞ limit, $\Omega RPA(R\u2192\u221e)=5.04\rho /m*$, in good agreement with the classical surface plasmon energy of $(4\pi /3)\rho /m*$.

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

where *C* = ∑_{ai}|*ρ*_{ai}|^{2} is a normalization constant, and in the RPA, the *X* and *Y* coefficients are

respectively, where Ω is the energy of the collective state and

We note that as *ε*_{g} → 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 *ε*_{g} increases, the *X* coefficients dominate and *Y*_{ai} → 0. This behavior is observed numerically, as shown in Fig. 2. The sum of de-excitation coefficients for the collective state, ∑_{ai}|*Y*_{ai}|^{2}, is therefore another measure of plasmonic character. From Fig. 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 *ε*_{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 ⟨*ia*|*jb*⟩ ≈ *λρ*_{ij}*ρ*_{ab} does not facilitate the solution of the eigenvalue problem. However, the largest-in-magnitude element ⟨*ia*|*ia*⟩ can be included exactly as it just shifts the HF bandgap *ε*_{g} → *ε*_{g} − ⟨*ia*|*ia*⟩. Within the schematic model, we take all such excitonic Coulomb integrals to be equal and obeying the scaling ⟨*ia*|*ia*⟩ = *c*/*R*, with *c* = 1.4. This approximation gives the CIS and TDHF excitation energies as

We note that the remaining excitonic Coulomb integrals ⟨*ia*|*jb*⟩ could be included via perturbation theory, though we do not pursue this here.

In Fig. 4, we show the performance of the analytic solutions of this schematic model compared to our numerical calculations for two additional densities in between those shown in Fig. 1. 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 behavior of our calculations based on simple scaling arguments on the density and radius and supports the physical interpretations given to the various single-excitation theories.

### D. Comparison to experiment and DFT-based approaches

Before concluding, in this subsection, we make a comparison with the experimental results and with a popular DFT-based approach. Here, we approximate the screening induced by the valence band electrons in ZnO nanoparticles by setting *ε*_{∞} = 3.72 in Eq. (1). 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 × 10^{20} cm^{−3}, against the experimental results from Ref. 7. We see again that EOM-CCSD identifies TDHF (left panel) as the most accurate singles theory and TDHF compares favorably with the experimental result over the full range of radii. Based on the comparison with the RHF gap, we propose that the primary transition is excitonic in nature below 3.5 nm and plasmonic in nature above 3.5 nm.

In Fig. 5 (right), we compare the experimental results to those predicted by density functional theory within the local density approximation (LDA). Specifically, we show the LDA Kohn–Sham gap and the excitation energy predicted by LDA plus the RPA. The LDA+RPA approach has been applied previously to our model^{20,22} and has found success in describing the experimental ZnO results.^{41,60,61} As usual, the LDA gap is smaller than the RHF gap due to the latter’s treatment of exchange, but the LDA+RPA result is in reasonable agreement with experiment, comparable to the result from TDHF. However, RPA only acts to increase the small LDA gap at all *R* due to the plasmonic effects; in contrast, TDHF lowers the large RHF gap at small *R* (due to excitonic effects) and increases the RHF gap at large *R* (due to the plasmonic effects). Therefore, ignoring exchange, as in the LDA+RPA approach, can produce an accurate result due to the cancellation of errors. It is known that in the condensed phase, the TDDFT excitation energy of semilocal functionals is equal to the Kohn–Sham gap and unable to describe excitonic effects,^{62} which is closely related to the description of excitations with charge-transfer character.^{63} In our doped nanoparticle calculations, we see that the accuracy of LDA+RPA with respect to TDHF is worst where there is strong quantum confinement and electron–hole attraction; at the smallest value of *R* considered, the LDA+RPA excitation energy is more than 0.2 eV higher than that of TDHF.

As discussed in Sec. II A, our simple treatment of screening ignores the dielectric effect of the solvent and ligands, which can be included by an electrostatic perturbation on the Coulomb interaction.^{64,65} A simple estimate based on the Drude model (*ε*_{m} = 2.25 for toluene) suggests that this effect could change our results by about 20%. However, we emphasize that the model still lacks atomistic detail and so a precise agreement is not to be expected.

## 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 (and therefore some hybrid functionals) 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 and surface, ligand, or solvent effects. The model is also unaware 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 toward atomistic or tight-binding 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;^{66} Delerue^{67} demonstrated the proper use and success of the RPA for an atomistic tight-binding model of doped nanoparticles. Second, the retention of attractive electron–hole “exchange” integrals is essential for an accurate wavefunction description of excitonic states, which may or may not be important depending on the regime studied. 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.^{68}

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## ACKNOWLEDGMENTS

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

### APPENDIX: COMPUTATIONAL DETAILS

#### 1. Basis set

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

where *r* is the radial coordinate, *α*_{nl} = *k*_{nl}/*R*, *k*_{nl} is the *n*th zero of the spherical Bessel function *j*_{l}(*r*), and the normalization constant is $Nnl=R3/2jl+1(knl)$. Each orbital is characterized by three quantum numbers *n*, *l*, *m*, with the limits *n* ≥ 1, *l* ≥ 0, and *m* = −*l*, …, *l*. The noninteracting orbital energies are given by

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

The naive way to grow the basis set is to add particle-in-a-sphere 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*_{max} is sufficient to converge the ground-state calculation. In order to capture the singly excited states, we add an additional shell *l*_{max} + 1 based on the dipole selection rule Δ*l* = ±1. The rapidly increasing degeneracy of the basis functions limits the number of electrons to 98, which we converge with 483 basis functions (*n*_{max} of [10, 9, 9, 8, 8, 7, 7, 7] for *l* = 0, …, *l*_{max} + 1 = 0, …, 7). For the EOM-CCSD calculations, we use a basis set that adds orbitals based on increasing energy and include over 350 basis functions for 72 electrons.

#### 2. Two-electron integrals

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

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

where ** x** = (

**,**

*r**σ*) is a combined space and spin variable.

The spherical harmonics $Ylm$ are generally complex. To maximize the symmetry of the two-electron integrals, we use the real form of the spherical harmonics, $yl\mu =\u2211mUm\mu lYlm$. With this choice, the two-electron integrals ⟨*pq*|*rs*⟩ are given by the Laplace expansion of the Coulomb potential,

where the angular integrals are

a linear combination of the integrals of three complex spherical harmonics,

which vanishes unless *l* + *l*_{1} + *l*_{2} = 2*g*, $g\u2208Z$, 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 factors.^{69} The radial integral of the normalized spherical Bessel functions *R*_{nl} is

and can be computed numerically.