It has been reported in photodoping experiments that localized surface plasmonic resonances can be sustained with electrons as few as 3. We performed first principles calculations of density functional theory, with the Hubbard U correction, to see if localized surface plasmonic resonances can also be sustained by doping a wide bandgap ZnO with few shallow donors of Ga. We distributed 3–6 dopants approximately uniformly, due to quasi-spherical geometry of the quantum dot, in the dilute doping limit. The uniform distribution of dopants in quantum dots has been reported experimentally. Although the dopant configurations are limited due to computational cost, our findings shed light on absorption trends. Results for quantum dots of 1.4 nm, passivated with pseudo-hydrogens, show that localized surface plasmonic resonances can be generated in the near infrared range. The absorption linewidths for such small-sized quantum dots are broad. We find that the resonance linewidth depends on the orientation of surfaces and the number of secondary peaks on the concentration of dopants. The absorption coefficients, as functions of the principal values of the dielectric tensor, indicate that an electric field with orientation parallel to that of the most symmetric surface will produce localized surface plasmonic resonances with high quality factors.

## I. INTRODUCTION

Plasmonic resonances, collective oscillations of conduction electrons induced by an external electromagnetic radiation, have been extensively studied in metallic as well as heavily doped semiconductor nanoparticles over the past decade.^{1,2} These resonances, known as localized surface plasmonic resonances (LSPRs), lead to the confinement of light to nanoscale dimensions and open novel potential applications in nanophotonics and the concomitant experimental and theoretical investigations leading to new physics. For large nanoparticles, with dimensions greater than the excitonic Bohr radius, the classical Maxwell's equations are solved for simple geometries such as spheres or numerically for more complex ones to determine the optical properties of a dielectric material surrounded by a dielectric medium.^{2} But when dimensions less than or equal to the excitonic Bohr radius are involved, quantum mechanics must be invoked, and the emerging field is known as quantum plasmonics.^{3}

According to conventional wisdom, LSPRs are a phenomenon of large electron density of metals or heavily doped semiconductors. Recently, however, photodoping experiments have been reported in which LSPRs were sustained with electrons as few as four per nanocrystal by Faucheaux *et al*.,^{4} and as low as three electrons by Gamelin and co-workers.^{5} This raises a fundamental question about the nature of LSPRs. Jain applied a simple model of free electrons confined to a potential well to explain his results.^{6} Gamelin *et al.* employed a semiclassical approach that includes quantum corrections to explain their findings.^{5} Concurrently, Govorov and co-workers developed a quantum theory of LSPR using time dependent density functional theory (DFT) based on the jellium model that ignores the discreteness of the impurities.^{7} Their model is applicable to quantum dots with many electrons that screen the impurity potential seen by a delocalized electron. In this work, in contrast, we focus only on a few electrons excited from impurity levels of shallow donors, taking the discreteness of the QD into account.

It has been underscored that investigations of the electronic structure of doped metal oxide nanocrystals constitute a frontier in a fundamental understanding of LSPRs.^{2,8} Here, we investigated the possibility of plasmonic oscillations by a few electrons in a nanocrystal of ZnO doped with Ga and dressed with pseudo-hydrogens to simulate ligands in a colloid. We employed the DFT + U method to find out if the plasmonic resonances are sustained with few electrons. The DFT + U is often used to obtain an experimental bandgap which is severely underestimated by the conventional DFT alone. To our knowledge, this is the first DFT + U calculation of LSPRs of its kind that takes into account the discreteness of the nanostructure as well as simulates ligands through passivation by pseudo-hydrogen atoms to mimic a quantum dot in a colloidal environment. Although computationally intensive, such a DFT calculation can serve at least two purposes. First, our work shows that it is possible to compute optical properties of dressed QDs of ∼1.5 nm, therefore, such a simulation can complement experiments in understanding their results as well as designing them. Furthermore, it can also serve to validate the more approximate calculations that combine DFT with less computationally expensive approaches such as tight binding (TB), for instance, DFT + TB, to access larger sized QDs. Even more exciting would be to use it for validation of codes that are actively being developed, known as linear-scaling DFT, provided the dopants are sparsely distributed. Sparsity can be justified in the diluted doping range by cutting off the tail end of the exponentially decreasing wavefunctions of localized basis functions. Linearized DFT computation scales as the number of atoms, N, instead of N^{3} in the conventional DFT, and will enable investigations of LSPRs across the whole size range, involving thousands of atoms, for quantum plasmonics.

## II. COMPUTATIONAL DETAILS

DFT is the most widely used and successful computational method to study the geometric, electronic, and magnetic structure of materials.^{9} The exchange-correlation functional of DFT cannot be calculated exactly, therefore, two approximations are known as the local density approximation (LDA) and the generalized gradient approximation (GGA) are often used in its implementation. The LDA and GGA approaches are well known to underestimate the bandgap. In the case of metallic oxides such as ZnO, this is due to the underestimation of the binding energy of the semi-core d states, which leads to larger hybridization with p valence states of oxygen resulting in the valence band edge maximum being pushed up.^{10} One way to correct this is to account for the strong correlation of the localized electrons by introducing the Hubbard U parameter. When this correction is applied only to the d orbitals of ZnO, however, the bandgap is still underestimated even at high U values.^{11,12} More recent theoretical studies have included the U correction to the p state of oxygen in addition to that of the d state.^{13,14} In this work, we followed the latter approach, and the U parameters for d and p orbitals from Ma *et al* are applied.^{15}

We performed GGA + U calculations as implemented in the Vienna *ab initio* simulation package (VASP) designed for a periodic structure.^{16,17} The calculation of the optical properties in VASP proceeds in two or three stages depending on the level of approximation: the independent particle approximation (IPA) or the random phase approximation (RPA). In either approximation, first, the geometry is optimized to determine the geometric and electronic structure of the ground state. The coordinates for the optimization were created as follows. First, we obtained the coordinates for bulk unit cells of ZnO from the Materials Project site^{18} and re-optimized them with U parameters included. The re-optimized unit cell was then enlarged to cut out a spherical quantum dot (QD). In a colloidal synthesis of QDs, the surface is passivated with ligands to arrest growth and prevent agglomeration of the QDs. To simulate the ligands, the surface of the QD was passivated with pseudo-hydrogen atoms.^{19} The passivated QD is then optimized to get the ground state energy and structure.

We chose the size of the QD based on computational feasibility as such computations involving plane-wave basis sets, such as in VASP, require a large amount of memory. The combination of ZnO and pseudo-hydrogen atoms increases rather fast, for example, from ∼217 atoms for a QD of 1.4 nm to ∼ 600 atoms for a QD of ∼2 nm. This is due to the shorter bond length of ZnO compared to, for instance, CdSe where a 2 nm packs approximately 300 atoms (including the pseudo-hydrogens) in our recent calculations.^{20} DFT, as noted above, scales as N^{3}, where N is the number of atoms, which makes the computation of large systems computationally too demanding. Experimentally, the approximately 2 nm QD is often cited as the lowest experimentally accessible size, although in some cases as low as 1.5 nm has been reported in the literature.^{21} Here, we chose QD of 1.4 nm due to computational cost, however, it is large enough to simulate the QDs in the lowest accessible range of ∼ 1.5 nm experimentally. Recently, LSPRs were reported for metallic nanoparticles of sizes less than 2 nm.^{22}

To employ the VASP code for periodic structure for the spherical QD, the QD was surrounded by a vacuum region of 12 Å that is large enough to avoid interaction with images in adjacent supercells. The 3d electrons in Zn and Ga were treated as valence electrons on a similar footing as 4 s electrons. The electronic wavefunctions were expanded in a plane-wave basis with an energy cutoff of 550 eV, which is 4/3 times the default value of VASP. The core electrons were replaced by a projector augmented wave (PAW) pseudopotentials supplied by VASP. The Perdew–Burke–Ernzerhof (PBE) GGA was used for the exchange-correlation functional.^{23} The U parameters were set to 10 eV and 7 eV for d and p orbitals, respectively.^{15} The Brillouin zone integration was performed at the gamma-point only for such a large supercell. The total energy is obtained by optimizing the geometry with the conjugant gradient scheme until a convergence criterion of 0.01 eV/Å for forces on the ions and energy convergence of 10^{−8} eV for self-consistent calculations are reached. Such very stringent constraints were required to reach convergence due to oscillations of energy values near the potential minimum.

The results from geometry optimization, namely, the coordinates and wave functions, were used as inputs for optical calculations. The eigenfunctions from the geometry optimization were used as wavefunctions for the optical calculation as it involves a new Hamiltonian that includes the interaction of the QD with radiation. The Hamiltonian matrix was diagonalized self-consistently until the energy as well as force tolerances, stated above, were met. The size of the matrix is determined by the number of unoccupied orbitals above the gap. A parameter of significance in optical calculations that controls the number of empty bands or unoccupied orbitals in the case of QDs is NBANDs in VASP. It must be kept large enough, subject to memory constraints of the computing hardware, to obtain valid results. In this work, NBANDS is taken to be five times the default value suggested by VASP.^{24} Once the orbitals of the Hamiltonian are determined, the frequency dependent dielectric matrix is calculated.

The linear optical properties in VASP are obtained from the complex dielectric function *ɛ*(*ω*) given by^{25}

$\epsilon \alpha \beta =\epsilon \alpha \beta 1+i\epsilon \alpha \beta 2$, where the imaginary part is

The indices *c* and *v* refer to conduction and valence band states, respectively. In this work, instead of the valence band, the indices *v* will correspond to impurity levels and the indices c to low-lying orbitals mentioned above. The angular frequency $\omega $ has the units of energy, and *u _{ck}* is the cell periodic part of the orbitals at k-point

**k**. The indices α and β denote Cartesian components, Ω is the volume of the unit cell, $\omega k$ are the k-point weights, and $e\alpha ,e\beta $ are Cartesian unit vectors. The long wavelength limit $q\u21920$determines the optical properties in the wavelength regime accessible to optical or electronic probes. The real part $\epsilon \alpha \beta 1$ of the dielectric tensor is obtained from the imaginary component by the usual Kramers–Kronig transformation. The absorption coefficients are then computed from the real and imaginary components of the dielectric function.

^{26}

## III. RESULTS AND DISCUSSION

Determining the electronic structure of finite systems such as QDs which involves supercells is challenging at the level of accuracy of the conventional DFT. The bulk unit cell, in contrast to the supercell of a QD with a much larger number of atoms, is much less costly to compute and can provide an approximate result for the QDs of sizes 1 nm or larger which have a bulk-like structure in the core.^{27} We optimized the bulk unit cell with DFT + U to obtain a bandgap of 3.37 eV in agreement with the experimental gap of ∼3.40 eV. In contrast, DFT calculation, without the U parameters, in GGA yield a bandgap of 0.6–0.8 eV.^{15}

In Fig. 1, the density of states for bulk wurtzite ZnO is presented for Brillouin zone integration involving 25 × 25 × 25 points. The opening up of the gap is clearly seen. The uppermost two valence bands, the first consisting mostly of the s and p orbitals and below it one consisting mostly of the localized d orbitals are seen. The conduction band states gradually rise from the gap and then more or less fluctuate around a constant value except around 10 eV.

Next, we computed the ground state of the undoped ZnO QD to get a HOMO–LUMO gap of 4.3 eV. This is consistent with what one should expect for a small-sized QD of 1.4 nm in diameter relative to the 3.37 eV obtained from the bulk unit cell above. The gap of 4.3 eV restricts the range of energies of photons that excite the electrons from the impurity levels if excitations from the occupied states below the gap are to be avoided. In subsequent optical calculations involving impurity levels, an upper bound of 4 eV was placed.

To proceed with doping, we first determined the formation energy of Ga to be significantly lower at the Zn site than at the O site, as expected. This is due to Ga and Zn being approximately of the same size and being much larger than O. Once the doping site is determined, we looked into the effect of the dopant environment on impurity levels. There are three distinct sites of Zn in the QD: a core site, a singly passivated surface site, and a doubly passivated surface site. The impurity levels of Ga at three sites were found to be 0.5618, 0.5677, and 0.4495 eV, respectively, below the LUMO, showing that Ga is a shallow donor. The impurity levels, as well as the energy spectrum around the gap, are shown in Fig. 2. The spread in the levels is mainly due to differences in the coordination numbers at the different sites. The results for the singly passivated donor at the surface being less than at the center are consistent with similar results for CdSe QD.^{28} The higher value for a doubly passivated donor is due to electron–electron repulsion from neighboring pseudo-hydrogen atoms. As seen in Fig. 5(b), the variations in the levels contribute to the broadening of the plasmonic resonance for the absorption coefficients.

Now that the doping site preference is determined, the question remains as to where the 3 or 4 electrons could originate from. For optical excitations of the order of impurity level, that is ∼0.5 eV, the number of electrons is equal to the number of Ga dopants. For that matter, it is the case for photon energies less than the gap. There are many possible configurations of the few Ga dopants given the 57 Zn sites in the 1.4 nm QD. The possibilities, however, can be cut down on an experimental basis. First and foremost, to get the few electrons on the order of 3 or 4, the doping concentration must be within the dilute limit. How would then these few dopants distribute themselves in the QD in such a limit? The distribution of dopants has been an issue as to whether they distribute uniformly or self-purify to the surface.^{29} Recently, however, it has been reported that uniform distribution can be realized,^{30} and on the basis of the experimental result, we assumed uniform distribution subject to the geometry of the QD which, in reality, is quasi-spherical. The dopants will diffuse during synthesis such that they are as distant from each other as possible to minimize the total energy before they settle into the Zn sites. The configurations for 3–6 dopants are given in Fig. 3. The dopants are placed at least two lattice constants away from each other to avoid nearest neighbor dopant ion—ion interactions involving the minor distortions in the positions of O atoms as will be shown below.

We optimized the structures to determine the ground state energy and configuration for each of the QDs. The results for the total energies with respect to that of the undoped QD and the energies per dopant are given in Table I. The dopants in the 5% and 8.8% QDs can be grouped together in the sense that they both have one dopant at the center of the QDs. The energy of the 8.8% is higher by 0.14 eV than the 5%, most likely due to electron–electron interaction between the hydrogenic electrons of Ga. Similarly, the 7% and the 11% can be grouped together; their atoms are distributed around the surface. The energy for 11% went up by ∼8 meV, much smaller than that of the first group. This is due to surface relaxation that gives more degrees of freedom for the dopants to lower their energy. The lowest energy per dopant is 5%. While the intent of the work is to see if LSPRs can be formed in such a dilute doping limit with few electrons, we note in passing that the limited result here suggests the more likely optimal structure would be one with fewer electrons based on energy comparisons.

. | 5% . | 7% . | 8.8% . | 11% . |
---|---|---|---|---|

Total energy | −8.625 51 | −11.073 55 | −13.658 55 | −16.559 81 |

Energy/dopant | −2.875 17 | −2.768 39 | −2.731 71 | −2.759 97 |

. | 5% . | 7% . | 8.8% . | 11% . |
---|---|---|---|---|

Total energy | −8.625 51 | −11.073 55 | −13.658 55 | −16.559 81 |

Energy/dopant | −2.875 17 | −2.768 39 | −2.731 71 | −2.759 97 |

We also looked into the geometrical structure around the Ga dopant. To see the effects of the substitution into the Zn site, we plotted the input and optimized output geometries together for 5% concentration as shown in Fig. 4. It is clearly seen the green Ga atoms cover up the Ga atoms except the one at the bottom that covers some portion of Ga. This suggests that the Ga dopant fits in at the Zn site and being slightly larger than Zn covers up Zn completely, but for the one at the bottom surface. The surrounding O atoms, on the other hand, clearly show that they have slightly shifted from each other. We present the shifts in one of the Ga dopants in Table II. The shifts range from 0.0 to 0.209 45 along the coordinate axes.

. | Δx
. | Δy
. | Δz
. |
---|---|---|---|

O1 | 0.000 44 | 0.209 45 | −0.051 71 |

O2 | 0.182 03 | 0.104 46 | −0.052 34 |

O3 | −0.181 05 | 0.104 70 | −0.052 27 |

O4 | −0.000 34 | 0.000 04 | 0.187 98 |

. | Δx
. | Δy
. | Δz
. |
---|---|---|---|

O1 | 0.000 44 | 0.209 45 | −0.051 71 |

O2 | 0.182 03 | 0.104 46 | −0.052 34 |

O3 | −0.181 05 | 0.104 70 | −0.052 27 |

O4 | −0.000 34 | 0.000 04 | 0.187 98 |

Compared to the Zn–O bond length of 1.835 Å, the shifts are relatively small. Moreover, the shifts are negative and positive in different directions such that the net shift is negligible as can be seen in Fig. 4.

After the ground state structure is calculated first, the resulting optimized structure and the electronic wavefunctions are used as inputs for the optical calculation in what is known as the independent particle approximation (IPA) according to VASP. The low-lying 1088 orbitals from the ground state calculations, determined by the tag NBANDS, were used in the self-consistent electronic structure optimization. As noted above, the NBANDS value applied is five times the suggested default value in VASP and was assumed more than enough to generate the low-lying unoccupied states of the Hamiltonian, with optical interaction included, after diagonalization of its matrix. After convergence of the self-consistent loop, ionic optimizations were repeated until the force tolerance stated above was reached. Following the geometric and energy optimization, the dielectric matrix is computed. Our computation was stopped at the level of IPA and did not proceed further to the next stage known as the random phase approximation (RPA). In the VASP context, the IPA does not take into account what is referred to as local effects, that is, variations in the periodic orbital $uk$ in Eq. (1). If local effects are considered, the wavefunction and its derivatives from the IPA calculation are used as input for the RPA. RPA requires a large amount of memory and is computationally very expensive for such a large supercell used in this work.

So, the question arises as to whether the IPA is justified? IPA, in the context of many body physics, is the neglecting of the electron–electron correlation. In DFT, the electron—electron correlation is considered through the exchange-correlation functional. In this calculation, the correlation is considered through the GGA of PBE. The IPA in the optical calculation of VASP is different. In the context of Eq. (1), IPA means ignoring the change in *u _{ck}*, the cell periodic part of the orbital, due to local effects. To minimize the local effects due to dopant, we chose Ga instead of Al, the commonly used dopant in ZnO, as the Ga ion is almost identical to Zn except for one protonic charge in the nucleus. This extra charge is screened by core electrons and would not be seen by optically excited, oscillating electrons in the low-lying orbitals above the gap. In addition, as seen in Fig. 4, the distortion at the dopant sites is minimal; therefore, we assume the potential seen by the oscillating electrons is identical within the core of the QD. As to the surface, there is variation if one takes the whole surface of the QD. However, oscillation of the electron due to the oscillating electric field is a local phenomenon and locally the electron sees the same average potential near or around the surface. In our construction of the surface with passivating pseudo-hydrogens, the two electrons along the bonds are maintained. Moreover, the distortions around Ga on the top surface in Fig. 4 are similar to that of Ga in the center. Ga at the bottom surface shows a very slight difference due to the relaxation of the surface, but the surrounding O atoms show similar distortions that are very minimal. Therefore, we assume that

*u*are almost the same as those near the core atoms at least to a first order approximation.

_{ck}The dielectric tensor is calculated according to Eq. (1). The reciprocal space calculations were performed at the gamma-point only for such a large supercell of ZnGaO. Then, the sum in Eq. (1) reduces to a single k-point and the computation is highly simplified. It also means that only inter-band optical transitions are considered in this work. For intra-band transitions, one could compute with a few more k points for low-lying states, but we found out that only gamma-point calculations, for which VASP has a separate code, were feasible. The imaginary and real parts of the dielectric tensor were processed using the vaspkit code^{26} to get the absorption coefficients plotted in Figs. 5 and 6.

Figure 5 is plotted as a scale for Fig. 6. For the three distinct types of Zn atoms stated above, tetrahedrally coordinated Zn atoms in the core (4C), triply coordinated atoms at the surface (3C), and doubly coordinated atoms at the surface (2C), we performed optical calculations. The amplitudes of the plots will serve as a scale, that is, to determine whether absorption coefficients in Fig. 6 are due to the single electron excitation or multiple ones. The plots show that the absorption coefficient peaks occur at ∼0.5 eV consistent with the impurity level distance from LUMO for single electron excitation. It is worth noting that the amplitudes for yy—principal values remain the same for all the sites while xx and zz are affected. This point will further be discussed below for multiple excitations.

Next, in Fig. 6, absorption coefficients for multiple electron excitations are plotted. The percentages refer to the number of Zn atoms (out of 57 in total) that were replaced by Ga atoms. Clearly, the amplitudes are larger than those in Fig. 5. For α_{yy}, there is a pronounced absorption peak at ∼1 eV in all three doping concentrations. There is a secondary peak at ∼1.75 eV for 5% and 7%, which is washed away for 8.8% and 11%. For α_{zz}, the first peak which is less than that for α_{yy} appears at ∼0.8 eV and the second one at ∼1.75 eV for 5% only and loses the peak features to broaden into a single band for 7%. The peak gets washed away at 8.8% and 11%. For α_{xx}, the spectrum is identical to that of α_{zz} for 5% and almost identical to α_{yy} for 11%. The first peak appears at ∼0.8 eV for all the concentrations, while the second peak gets washed away for all but the 5%. The spectra clearly show that they depend both on the direction and concentration. α_{yy} has the spectra with the highest quality factor, relative to α_{xx} and α_{zz}, and α_{zz} has the lowest. This is due to the symmetric nature of the distribution of the atoms on different surfaces as shown in Fig. 7, which affects the potential seen by the excited electrons as they oscillate.

As stated above, we assume that the oscillating electrons see the same potential in the core except at the surface due to the termination of the structure as well as the ligands. The surfaces of the wurtzite ZnO are different along the axes as seen in Fig. 7, and their potentials are different accordingly. The difference in potential results in variation in absorption coefficients along the axes. However, on the same plane, we assume that the oscillating electrons see the same average potential locally around the surface. Therefore, if the polarization of the electric field can be related to the structure of the facet of QD, it would result in LSPRs with high quality factors. It has been reported that the orientations of facets of individual QDs can be controlled through small passivating ligands.^{31} This implies that the orientation of the oscillating electric field can be correlated with the orientation of the facets in the colloid. We also assume that the QD is free of defects such as oxygen vacancy. Even if there is a vacancy, in general, such a vacancy forms an acceptor level near the valence band and is energetically separate from the shallow donor level of interest in this work.^{20}

The surfaces along the z direction are unsymmetric, terminating with Zn or Ga atoms at (0001) and oxygen atoms at $(0001\xaf)$. The xy surfaces have a symmetric distribution of atoms, terminating with the almost flat surface in the y direction and surfaces in the x direction that edge out at the center and taper off. It is due to the symmetry and uniform distribution of the atoms along the y direction that the resonant peaks are most pronounced with relatively higher quality factors than those in the other two directions of x and z. The xz and zy surfaces terminate with polar surfaces along the z axis. The correlation between the surfaces and the sharpness of the resonance underscores the importance of choosing surfaces with a symmetric distribution of atoms in crystals with different packings of atoms along crystallographic directions as in the wurtzite structure of this work.

In addition to symmetry, the concentration and how the dopants are distributed affect the quality factor of the resonance. The concentration, as well as the distribution, did not affect the primary peak for α_{yy} indicating that the resonant oscillations are robust. On the other hand, α_{zz} with unsymmetric surfaces is most affected and loses the consistent peak feature as a function of concentration except for 5%. 7% and 11% have all the dopants on the surface (except one near the surface) and the peaks of their spectra are affected most for α_{xx} and α_{zz} relative to 5% and 8.8%. 5% and 8.8% have a dopant in the center. Although the samples are limited, we conjecture that a uniform distribution throughout the QD would produce spectral features with distinct peaks. Moreover, the size of our dopant is at the lowermost end of accessible QD sizes experimentally. Such a small-sized dopant has small facets, sharp edges, and as many surface atoms as interior atoms which results in variations of potentials seen by the oscillating electrons. Obviously, such variations in potential result in significant scattering and degradation of the quality factor.

The primary peak at ∼1 eV corresponds to near infrared radiation. It is within the range of values for excitation of 2 electrons, that is, from 0.8990 to 1.1354 eV, suggesting the oscillation of coupled two electrons may be the dominant normal mode of the oscillation. This absorption energy can be manipulated by varying the dopant, the host, or both. Thus, doping can open up opportunities to extend the absorption spectrum to the far infrared regime through the interplay of hosts and dopants. Wide bandgap semiconductors such as ZnO, therefore, are great candidates for extending the absorption spectrum into the far infrared regime.

## IV. CONCLUSIONS

Our results show how the absorption spectra are related to the geometric and electronic structure of the nanocrystal. The absorption quality factor strongly depends on the symmetric distribution of atoms on surfaces and the concentration of dopants. Although our samples are limited with respect to the distribution of dopants, we conjecture that a uniform distribution of dopants throughout the nanocrystal will produce LSPRs with high quality factors. The primary peak at ∼ 1 eV corresponds to absorption in the near infrared regime. Thus, doping can open opportunities to extend the absorption spectrum to the far infrared regime through the interplay of hosts and dopants. Finally, the LSPRs for strongly confined colloidal QDs at the lowest end of experimentally accessible QDs of ∼ 1.5 nm have LSPRs with low quality factors due to inhomogeneities in surface potentials.

## ACKNOWLEDGMENTS

D.D. and M.D.M were supported by the National Science Foundation Grant No. DMR-2013854. M.D.M. acknowledges computer time allocation (No. TG-DMR100055) from the Extreme Science and Engineering Discovery Environment (XSEDE) for Stampede at TACC, which was supported by the National Science Foundation Grant No. ACI-1548562. M.D.M. also acknowledges computer time at the carbon cluster of the Center for Nanoscale Materials at the Argonne National Laboratory and support from Dr. M. Chan.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflict of interest to disclose.

## DATA AVAILABILITY

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