Advances in Optical Technologies

Advances in Optical Technologies / 2008 / Article
Special Issue

Silicon Photonics

View this Special Issue

Review Article | Open Access

Volume 2008 |Article ID 279502 |

V. A. Belyakov, V. A. Burdov, R. Lockwood, A. Meldrum, "Silicon Nanocrystals: Fundamental Theory and Implications for Stimulated Emission", Advances in Optical Technologies, vol. 2008, Article ID 279502, 32 pages, 2008.

Silicon Nanocrystals: Fundamental Theory and Implications for Stimulated Emission

Academic Editor: D. Lockwood
Received18 Feb 2008
Accepted09 May 2008
Published29 Jun 2008


Silicon nanocrystals (NCs) represent one of the most promising material systems for light emission applications in microphotonics. In recent years, several groups have reported on the observation of optical gain or stimulated emission in silicon NCs or in porous silicon (PSi). These results suggest that silicon-NC-based waveguide amplifiers or silicon lasers are achievable. However, in order to obtain clear and reproducible evidence of stimulated emission, it is necessary to understand the physical mechanisms at work in the light emission process. In this paper, we report on the detailed theoretical aspects of the energy levels and recombination rates in doped and undoped Si NCs, and we discuss the effects of energy transfer mechanisms. The theoretical calculations are extended toward computational simulations of ensembles of interacting nanocrystals. We will show that inhomogeneous broadening and energy transfer remain significant problems that must be overcome in order to improve the gain profile and to minimize nonradiative effects. Finally, we suggest means by which these objectives may be achieved.

1. Introduction

Silicon is the most widespread semiconductor in modern microelectronics technologies. Its natural abundance, low cost, and high purity, as well as the high electronic quality of the Si/SiO2 interface, have led to its overwhelming dominance in microelectronic devices. Nevertheless, the use of silicon in optoelectronics remains highly limited. This state of affairs has remained, in fact, almost unchanged because of a fundamental property of the silicon band structure—the indirect band gap.

The indirect radiative interband transitions in bulk Si are strongly suppressed because an emitted photon cannot satisfy the momentum conservation law for transitions from the conduction-band minimum (Δ-point) to the top of the valence band (Γ-point). The photon wave vector is about three orders of magnitude less than that required for the transition between the Δ- and Γ-points. This difference in 𝐤-space is 𝑘Δ=0.86×2𝜋/𝑎0, with 𝑎0 being the lattice constant of silicon, equal to 5.43 Å. The electron-hole radiative recombination in the bulk material is exactly forbidden unless additional mechanisms allowing the momentum to be conserved are involved in the recombination process. The most probable means to have a radiative indirect transition without breakdown of the momentum conservation law is via phonon absorption or emission. However the electron-phonon interaction is weak; consequently, phonon assistance is a low-probability (and hence slow) process. This leads to a substantial increase of the total recombination time, and a decrease of the recombination probability compared to the direct no-phonon ΓΓ radiative transitions in direct-gap semiconductors. In this sense, we call such transitions in Si “strongly suppressed.”

The discovery of visible-range emission from nanocrystalline [1] and porous [24] silicon in the early 1990s suggested that some of the problems associated with the silicon band structure might be overcome in nanoscale crystallites hosted in a widegap dielectric matrix like SiO2, in order to create a strong confining potential for carriers inside the nanocrystal. Electronic states become localized within the NC and the momentum distribution spreads due to the Heisenberg uncertainty relations. In other words, the wave functions consist of plane waves with all possible wave vectors including both 𝑘𝑘Δ for holes and 𝑘0 for electrons, respectively. Thus, the momentum conservation law is not violated in this case, which yields a nonzero probability of the ΔΓ radiative transition even in the absence of phonons.

Indeed, some time later, efficient visible photoluminescence (PL) from silicon nanocrystals was experimentally demonstrated, (e.g., [5]) and attributed to the transitions between confined electron and hole states inside the nanocrystal [6, 7] (the so-called quantum confinement effect) or between interface states [8, 9]. However, the emission quantum efficiency in Si nanocrystals remains low compared to the direct gap III-V or II-VI materials. This is naturally explained by the small “weights” of the plane waves with 𝑘𝑘Δ for the valence states and 𝑘0 for the conduction states. Thus, improvement of the light emission efficiency of Si quantum dots remains a challenge for optoelectronic technologies.

As a means to modify the optical properties of silicon crystallites, doping with shallow impurities has been suggested [1018]. In some cases (depending on the conditions and methods of preparation), the emitting properties of the dots were improved significantly. In particular, the PL intensity was several times greater when the nanocrystals were doped with phosphorus [1012, 14] or codoped with phosphorus and boron [15, 17]. The origin of this phenomenon is not fully understood at the present time, and we will touch on the problem of impurity states in silicon nanocrystals in this review.

PL experiments are not usually carried out with a single quantum dot but rather with large ensembles. As a result, interpretation of the experimental data is attendant with difficulties because of the associated inhomogeneous broadening and various collective effects that can occur as a result of the mutual influence of the nanocrystals in the ensemble. Some recent studies have reported on the PL spectra of individual silicon quantum dots [19, 20] and this is shedding new light on the physics involved in the light emission process. A principal distinction in the emission of a single NC and a nanocrystal ensemble lies in the various mechanisms of interaction amongst the nanocrystals in an ensemble. In particular, in solid nanocrystal arrays, some collective effects caused by electron, photon, and phonon transfer between the dots can strongly influence the luminescence dynamics of the nanocrystals in comparison with the case of isolated NCs. The size distribution of the nanocrystals can play an essential role in the excitation exchange between the clusters.

One possible mechanism of the NC-NC interaction is the direct tunneling of excited carriers from one quantum dot to another [2123]. As an example, we can imagine two (or more) closely separated quantum dots with different sizes. Presumably, excited carriers in the smaller nanocrystal, having higher quantized energy levels, may relax either to the valence band of this nanocrystal or to the conduction band of the adjacent nanocrystal with a larger size. The transition to the valence band is indirect, and therefore suppressed. Meanwhile, the transfer to the neighboring quantum dot occurs as if within the same (conduction) band, and may be more probable in the case where the nanocrystals are sufficiently close. If so, then the smaller nanocrystals will inject their own excited carriers into the larger nanocrystals, in which radiative interband transition will subsequently take place. This idea implies the possibility of an optical nanofountain [24]—a device emitting photons from the area where quantum dots with larger sizes are concentrated.

An additional nonradiative mechanism of energy transfer between adjacent nanocrystals becomes possible in dense arrays. This is the so-called Forster mechanism originating from the dipole-dipole (or higher-order multipole) interaction between excitons in different quantum dots [2527]. Because of the dipole-dipole interaction, the electron-hole excitations can “travel” throughout the nanocrystals without charge transfer. Such transitions have been observed in solid arrays of CdSe nanocrystals [28, 29]. Some theoretical aspects of the Forster transfer in silicon quantum dots were recently discussed by Allan and Delerue [30].

Various collective effects in quantum-dot arrays can result in a complicated time dependence of the experimental PL decay. As a rule, the decay is described by the stretched exponential function exp{(𝑡/𝜏)𝛽} with 𝛽 less than unity [3135]. Evidently, the value of β can be due to a size distribution of the nanocrystals (and hence a distribution of the recombination rates). Theoretical explanations of the stretched exponential from the point of view of the multiexponential decay and associated Laplace transforms have been given in [36].

This is, in short, the framework in which we will address the present review. At first we will briefly describe the experimental situation in this field in Section 2. Then the single-particle electronic structure and many-body corrections will be considered in Section 3 for a single silicon quantum dot as a basis for treatment of more complicated problem of solid quantum-dot ensembles. The latter will be reviewed in Section 4. Finally, Section 5 presents a general conclusion and a commentary on the potential for stimulated emission in silicon nanocrystal ensembles.

2. Experimental Results: Overview

Different methods can be used for preparation of silicon nanocrystals, for instance, Si ion implantation [3741], chemical vapor deposition [42], magnetron sputtering [43, 44], colloidal synthesis [45], electron beam evaporation [46, 47], and some others. A high-temperature thermal treatment is generally required in order to precipitate the crystallites. All these techniques allow one to form silicon NCs with sizes predominantly ranging from 2–6 nm. Their electronic structure and optical properties depend, of course, on the preparation conditions and method of fabrication. However, there are some common properties typical for silicon NCs, independent of the fabrication technique employed. In particular, the nanocrystals' surroundings, either vacuum or some host material like silicon dioxide, represent a high potential barrier for carriers of both kinds. Such a barrier is often referred to as a confining potential that mainly defines the energy spectrum of the nanocrystal. In what follows we will discuss some manifestations of the quantum-confinement effect in experiments and theoretical models. First, let us briefly discuss experimental data on the PL spectra in silicon quantum dots.

2.1. Size Dependence

As has been mentioned above, silicon NCs are capable of emitting electromagnetic energy in the visible spectrum. This is in contrast with bulk silicon, in which energy of the interband transition corresponds to the silicon bandgap energy of 1.12 eV. The increase (decrease) of the photon frequency (wavelength) in nanocrystals compared to the bulk material is a universal phenomenon taking place in various semiconductor materials and quantum dots. In the general case one may say that the energy of the emitted photon increases as the nanocrystal size decreases. Such an increase is usually called a “blueshift” because the photon energy shifts toward the shorter-wavelength side of the visible spectrum. This blueshift is illustrated for Si NCs in Figure 1. Here, the mean NC size is controlled via the excess Si concentration, with the smallest NCs occurring in the most silicon deficient samples. The PL intensity has not been normalized here; the drop in intensity on the silicon-poor side of the compositional map is due to the lower number density of NCs, and on the silicon-rich side it is due to the opening of nonradiative pathways in large and highly interconnected nanoclusters.

In the simplest model of an infinitely strong confining potential (i.e., infinitely high potential barriers at the dot boundary) it is possible to estimate the energies of electrons and holes localized inside the nanocrystal as proportional to 𝑅2, where 𝑅 is the nanocrystal radius. Thus, the optical gap may be calculated as 𝜀𝑔(𝑅)=𝜀𝑔+𝐴/𝑅2, where 𝜀𝑔 is the bandgap of bulk silicon, 𝐴 is some positive constant, and 𝐴/𝑅2 represents the total energy of the non-interacting electron-hole pairs inside the dot. It is obvious that: (i) the nanocrystal gap must be always greater than the bulk one due to the additional positive term 𝐴/𝑅2, and (ii) the photon energy, equal to the energy of the lowest electron-hole transition, increases as the dot size decreases.

In experimental work carried out over the past fifteen years with silicon nanocrystals, the optical-gap dependence on the dot size was measured and discussed extensively. Although it is impractical to cite all the papers dealing with this topic, a sampling is given in [8, 23, 4856]. The data in these papers is summarized in Figure 2. There is a fairly large spread in the calculated values of the optical gaps as a function of NC diameter. Presumably, several factors influencing the accuracy of the optical-gap measurements are as follows. First, it is difficult to determine exactly the dot sizes and the size distribution in a luminescent ensemble of NCs. Second, using the mean size in an ensemble of NCs in a diagram like Figure 2 can be misleading, since it is possible that the observed PL peak does not correspond exactly to the mean size but instead to the largest PL rate. Third and more practically, the nanocrystals studied by different research groups have been prepared by different methods. As a result, the NCs have different surroundings, surface bonds, and shapes, all of which could lead to scatter in the experimental data. Finally, as shown in later sections, NC-NC interactions can play a dominant role in the emission spectrum.

Furthermore, the blue-shift energy Δ𝜀𝑔(𝑅), determined from the experiments, does not obey the law 𝐴/𝑅2 following from the simplest quantum-mechanical model. This dependence is rather Δ𝜀𝑔(𝑅)𝑅𝑏 with 1<𝑏<1.5 or an even weaker dependence on 𝑅 in some cases. Such behavior has been a reason for supposing a key role for interface states in the radiative recombination process [51]. The question about the origin of the radiative electron-hole transitions in the nanocrystals is remains under extensive debate even today.

Nevertheless, employing the principle of Occam's razor, the deviations from the simplest model predictions for high-energy luminescence (when the peak-position energy exceeds the bulk gap 𝜀𝑔) may be explained within the quantum confinement framework without necessarily the need to invoke sub-gap radiative centers. Let us, first, take into account the finiteness of the potential barriers, leading to some weakening of the quantum confinement and, as a consequence, to a more complicated and less steep size-dependence of Δ𝜀𝑔(𝑅) other than 𝑅2. In addition, electron-hole Coulomb interaction contributes a further 𝑅1 dependence in the blue-shift energy. Both these mechanisms result in a more gradual dependence of the optical gap on the dot radius, which is in fact observed in experiments.

At the same time, low-energy PL with photon energies less than the gap of bulk silicon are also observed in experiments, see, for example, [57]. In this work an extra luminescence peak arose at about 0.9 eV at low temperatures, and its position remained almost unchanged for nanocrystals of various sizes. This may be indeed treated as an attribute of the surface states which can appear inside the band gap of bulk silicon. This possibility and the associated pair-wise trapping rates have been discussed and derived theoretically by Lannoo et al. [58].

2.2. Decay Rate

Time-resolved studies of the PL of silicon NCs demonstrate near-exponential decays of the photoluminescence intensity. Typical decay times 𝜏PL vary within a wide range of about 1 to1000 𝜇s depending on the particle size, temperature, detected wavelength, method of preparation, and so forth [31, 32, 51, 5962]. Such lifetimes are indeed large, and this is an explicit indication of an indirect-band-gap material. Contrary to direct-gap III-V or II-VI compounds, the drop of luminescence intensity for silicon crystallites is sufficiently slow to result in a low PL irradiance compared to direct-gap NCs, even for comparable quantum efficiencies.

The characteristic time of the PL decay is determined by two different recombination mechanisms: radiative, with typical time 𝜏𝑅; and nonradiative, with typical time 𝜏NR. The rate of the photoluminescence decay is the sum of the radiative and nonradiative recombination rates: 1/𝜏PL=1/𝜏𝑅+1/𝜏NR. In the case where one of the two times is much less than the other, the photoluminescence lifetime coincides with the smallest time. Usually in silicon, the nonradiative channel turns out to be faster compared to the radiative one (see, e.g., [31]), so that the PL lifetime equals 𝜏NR. The quantum yield 𝜂, which is proportional to the PL intensity, may be defined as the “weight” of the radiative channel in the recombination process: 𝜂=𝜏𝑅1/𝜏1PL. If 𝜏𝑅𝜏NR, then 𝜂=𝜏NR/𝜏𝑅1. Obviously, one can use two physically different ways to increase 𝜂. The first one is an increase of the radiative transitions, while the second is a decrease of the nonradiative processes. The second way appears easier in practice because the nonradiative processes can be influenced by preparation conditions, as has been recently demonstrated by Miura et al. [62]. On the contrary, according to [62], radiative recombination cannot be controlled by the preparation conditions. Therefore, control of the radiative channel efficiency seems to be fairly difficult to achieve in practice.

Sometimes, both radiative and nonradiative channels contribute comparably in the interband recombination in silicon nanocrystals. In these cases the quantum efficiency becomes very high—possibly on the order of tens of percent [62]. It should be noted, however, that these situations are possible, presumably, at pumping levels corresponding to no more than one excited carrier in every quantum dot. Then, traps such as surface defects may be the main sources of nonradiative recombination. Alternatively, when the excitation power is very high and several electrons are excited in the dots, the Auger process becomes possible. Since Auger interactions are fast, the radiative channel will be “shunted” in this case. Accordingly, 𝜂 tends to some small value, on the order of a fraction of a percent.

In cases when the nonradiative channel is mostly “closed,” that is, 𝜏𝑅𝜏NR, silicon nanocrystals demonstrate strong temperature dependence of the PL lifetime [23, 33, 57, 63, 64]. In particular, the lifetime becomes smaller as the temperature decreases from ~300 to 4 K. This behavior is explained within the framework of the singlet-triplet two-level model suggested by Calcott et al. [3, 4] for porous silicon. According to this model the exciton state splits due to the electron-hole exchange interaction into the upper singlet with relatively short lifetime 𝜏, and lower triplet with about 2-3 orders-of-magnitude longer lifetime 𝜏. The decay rate is then defined by 1𝜏PL1𝜏𝑅=3/𝜏+1/𝜏ΔexpS-T/𝑘𝑇Δ3+expS-T/𝑘𝑇,(1) where ΔS-T is the energy of the singlet-triplet splitting, 𝑘 is the Boltzmann constant, and 𝑇 is the temperature. The energy of the singlet-triplet splitting depends on the nanocrystal size and increases from a few meV to 10–20 meV for photon energies increasing from 1.4 eV to 2.2 eV [57, 64]. At low temperature, when ΔS-T𝑘𝑇 the decay of the singlet state is strongly suppressed, and the total lifetime 𝜏PL coincides with the longer time 𝜏. In the opposite case (when ΔS-T𝑘𝑇), 𝜏PL is close to 4𝜏. As a result, the PL decay occurs substantially faster at low temperature. More rigorous theoretical analysis [6568] reveals a rich fine structure in the excitonic spectrum of silicon NCs. Exciton states with different symmetry may be ascribed to so-called bright and dark excitons which have essentially different recombination lifetimes.

2.3. Enhancement and Quenching of Photoluminescence due to Impurities

During the past decade, doping Si NCs with shallow impurities has been explored as a potential route for modifying the luminescent properties of silicon nanocrystals. The main questions are the following: are we able to improve the emission of silicon quantum dots by incorporating donors or acceptors (or both)? If so, then the next question is what is the mechanism for the improvement? The results obtained by different research groups indicate that the PL spectra of silicon nanocrystals are, indeed, sensitive to doping. However, both an increase [1013] and a decrease [18] of the PL intensity have been observed. This, again, depends on the sample fabrication method and the type of doping.

Silicon NCs doped with phosphorus [1013, 16, 69], boron [18], or co-doped with P and B [14, 15, 17] have recently been investigated. At the same time, the authors of [7072] reported on experiments with silicon NCs doped not only with phosphorus or boron but also with hydrogen and nitrogen. Summarizing the results, one can conclude that doping with phosphorus or codoping with phosphorus and boron can enhance the PL intensity of the nanocrystals by several times, as initially reported in [10]. Meanwhile, no enhancement is observed for NCs doped with boron, hydrogen, or nitrogen. Moreover, the degree of quenching of the luminescence depends on the annealing procedure [18, 7072].

In accordance with the results of Miura et al. [12], Fujii et al. [13], and Mikhaylov et al. [72] reported that an increase of the phosphorus concentration leads to an initial rise of the photoluminescence peak, whereupon the PL intensity subsequently decreases with a further increase in doping. Moreover, the P concentration corresponding to the maximum intensity becomes greater with decreasing NC size [1113]. Therefore, in some cases, where the NCs had small sizes of about 3-4 nm [11, 70], only a monotonic rise of the PL peak was observed as P concentration increased. One possible reason for this phenomenon lies in the passivation of dangling bonds (i.e., neutral or charged Pb centers) at the NC surfaces. The absence of the enhancement effect for heavily P-doped crystallites at a given size has been explained by more extensive Auger recombination [14, 15, 17] or a coalescence mechanism [72]. The authors of [72] also pointed out the essential distinction in the photolumine scence of the P-doped nanocrystals formed by ion implantation at 1000°C and 1100°C. The former demonstrate up to 5 to 6 times intensity enhancement, while for the latter ony PL quenching was found (Figure 3).

As a method to reduce the Auger recombination, simultaneous codoping with phosphorus and boron, thereby compensating the numbers of donors and acceptors inside the crystallite, has been proposed [14, 15, 17]. PL enhancement was successfully obtained for larger nanocrystals of ~6–8 nm in diameter. In this case, a large shift of the PL peak, even below the optical gap of bulk silicon, takes place when increasing P concentration with fixed concentration of boron [14]. This may be evidence of the formation of bulk-like impurity states. At the same time, codoping with high concentration of both donor- and acceptor-type impurities result in a strong chemical shift of the ground-state energy levels in the conduction and valence bands [73], which may substantially reduce the NC optical gap. In the latter case no assumptions on the bulk-like form of electron states in silicon quantum dots are required, because the effect of the strong chemical shift is exclusively due to the confining potential of the dot.

3. Single-Dot Problem

We next turn to the quantum-mechanical analysis of various aspects of light emission in silicon nanocrystals. It is quite clear, of course, that a consideration of electronic structure, electron-hole recombination, and the role of impurities in a single quantum dot provides a basis for understanding optical properties of the quantum-dot ensembles. Therefore the present chapter first discusses some theoretical models concerning the interband transitions in an individual quantum dot. First, we will discuss the optical gap and energy spectra of undoped nanocrystals. Then the “band” structure and charge distribution in doped dots will be analyzed. Finally, the recombination lifetimes will be calculated for undoped and doped crystallites.

3.1. Electronic Structure of Perfect Silicon Nanocrystals

The first step is a calculation of the ground-state energies both for the conduction and valence bands of a nanocrystal. Note that here and throughout the paper we will use terminology that is typical for the bulk material. Of course, no genuine energy bands exist in the nanocrystal because of its finite size. Nevertheless, the states above and below the NC optical gap originate from the size-quantized states of conduction and valence bands of the bulk crystal, respectively. Moreover, in the following we mainly intend to employ the envelope-function approximation (𝐤𝐩 method) to describe electronic states of the nanocrystals. Consequently, such the terminology is convenient and should not lead to any misunderstandings. The states above and below the optical gap will be further referred to as conduction and valence bands, as in the bulk. Different approaches are used when calculating the gap of silicon nanocrystals. Among them, the simplest is the single-particle approximation. Usually, the single-particle gap exceeds real optical gap because the electron-hole Coulomb interaction reduces the total energy of the system. For this reason, the values of the single-particle gap may be treated as an upper bound for the optical gap. In this section we discuss electron and hole states of a “pure” crystallite, which has no defects and a perfect diamond lattice all over its whole volume. The results are compared with the simpler approximations discussed in the introduction.

3.1.1. Optical Gap. Numerical Results

It is possible to separate all calculation methods applied for computation of the electronic spectra in pure NCs into two groups. In the first group we have the numerical first-principles, empirical, and semi-empirical methods, such as density functional theory (DFT) [68, 7478], pseudopotential (PP) method (or various combinations of PP and DFT) [7982], tight-binding (TB) models [8386], combined Green's function (GW-approximation), and the Bethe-Salpeter equation (BSE) technique [8789]. All these methods require considerable computational time, a problem that worsens as the nanocrystal size is increased. Analytical methods such as the effective mass approximation or 𝐤𝐩 method [67, 9095] may be ascribed to the second group. In following we combine these two methods into a single envelope function approximation (EFA). The latter allows one to perform calculations without any restrictions on cluster sizes. Additionally, this approach can be used to find the higher excited states of the nanocrystal.

In the present subsection, we summarize the results of some of the methods discussed above, which illustrate dependence of the optical gap on the nanocrystal size (Figure 4). In all calculations, any dangling bonds at the nanocrystal surface were considered to be H-passivated. As is seen in the figure, there is some scatter of data due to the different computational methods and different approximations utilized in the computational procedure.

In particular, taking into account the self-energy correction or electron-hole interaction can considerably change the magnitude of the nanocrystal optical gap, in opposite directions, as has been pointed out by many authors (see, e.g., [81, 88, 89]). In particular, it is known that the local density approximation of the DFT tends to underestimate the optical-gap (similar results happen for silicon nanowires) [96, 97] compared to those obtained in experiments. For a more accurate account of the self-energy corrections and excitonic effects, the GW+BSE approach is often applied. The combination of the GW and BSE methods yields almost complete compensation of the exciton and the self-energy corrections to the single-particle optical gap [88]. The same effect was also found by Franceschetti and Zunger [81] within the framework of the PP method. Unfortunately, the GW+BSE approach is employed for nanocrystals whose diameters do not exceed 1 nm [8789]. Therefore, we did not depict the results of this approach in Figure 4.

3.1.2. Electron and Hole Spectra in a Dot

Let us now discuss electronic spectra in both the valence and conduction bands of the quantum dot. As has been pointed out, for this purpose it is convenient to use the EFA. Here we describe in more detail not only the lowest but also the higher excited electronic states as a function of size.

We assume the total potential energy of the electron inside the quantum dot to be of the following form𝑈𝐫=𝑈0𝑟+𝑉sp𝑟.(2) Here 𝑈0(𝑟) is the confining potential equal to zero inside and infinity outside the dot. The second part 𝑉sp(𝑟) describes an interaction between the electron and its image, arising due to the charge polarization on the boundary between the silicon nanocrystal and its dielectric surrounding. 𝑉sp(𝑟) is often referred to as a self-polarization term. It can be represented as𝑉sp𝑒(𝑟)=2𝜀𝑠𝜀𝑑2𝜀𝑠𝑅𝑙=0𝑙+1𝑙𝜀𝑠+(𝑙+1)𝜀𝑑𝑟2𝑙𝑅2𝑙,(3) where 𝜀𝑠 and 𝜀𝑑 are the static dielectric constants of silicon and the dielectric matrix, respectively. In order to find the electronic states in the dot, we have to solve the single-particle Schrödinger-like equation for the envelope function vector |Φ:𝐫||||𝐻+𝑈Φ=𝐸Φ.(4) Here, 𝐻 is the bulk 𝐤𝐩 Hamiltonian operator acting on the six-dimensional (6D) envelope-function vectors |Φ, and 𝐸 is the electron energy. The components of the 6D-vector |Φ are slowly varied expansion coefficients Φ𝑗(𝐫) of the total wave function in the Bloch-state basis.

The electronic states in the valence and conduction bands can be determined separately. Starting with the valence band, the Bloch-state for the Γ-point is |𝑌𝑍|, |𝑋𝑍|, |𝑋𝑌|, |𝑌𝑍|, |𝑋𝑍|, |𝑋𝑌|, where | and | are “up” and “down” spinors, respectively, and the Bloch functions |𝑌𝑍, |𝑋𝑍, |𝑋𝑌 belong to the irreducible representation Γ25. The 𝐤𝐩 Hamiltonian matrix 𝐻 in 4 is the sum of three parts: 𝐻=𝐻(0)+𝐻(1)+𝐻(so). Here 𝐻(0)=(𝐩2/2𝑚)×𝐈 is the isotropic part obtained as the average of the total 𝐤𝐩 matrix over all angles in 𝐩-space. 𝐻(1) is the anisotropic part that can be represented by two equivalent 3×3 blocks situated on the main diagonal of the total 6×6  𝐤𝐩 matrix:𝐻1=𝐻300𝐻3,𝐻31=2𝑚0𝑄1𝑁𝑝𝑥𝑝𝑦𝑁𝑝𝑥𝑝𝑧𝑁𝑝𝑥𝑝𝑦𝑄2𝑁𝑝𝑦𝑝𝑧𝑁𝑝𝑥𝑝𝑧𝑁𝑝𝑦𝑝𝑧𝑄3,(5) where 𝑄1 denotes ((𝐿𝑀)/3)(3𝑝2𝑥𝐩2), 𝑄2 denotes ((𝐿𝑀)/3)(3𝑝2𝑦𝐩2), and 𝑄3 denotes ((𝐿𝑀)/3)(3𝑝2𝑧𝐩2).

Finally, the term𝐻so=130𝑖Δ000Δ𝑖Δ0000𝑖Δ000Δ𝑖Δ000Δ0𝑖Δ000𝑖Δ𝑖Δ00Δ𝑖Δ0000(6) describes the spin-orbit interaction. We have introduced above the spin-orbit energy Δ equal to 44 meV for silicon, the 6×6 identity matrix 𝐈, and the hole effective mass 𝑚=3𝑚0/(𝐿+2𝑀), where the numbers 𝐿, 𝑀, 𝑁 are dimensionless empirical parameters in the 𝐤𝐩 Hamiltonian operator in the valence band. For silicon they equal 5.8, 3.43, and 8.61, respectively [98].

Because of the isotropic and diagonal form of the operator 𝐻(0)+𝑈0(𝑟), it is possible to classify its eigenstates similarly to atomic systems using common terminology such as 𝑠-, 𝑝-, and 𝑑-type states, and so forth [92]. Accordingly, one may expand the components of the envelope-function vectors over these eigenstates asΦ𝑗(𝐫)=𝛼𝐶𝑗𝛼||𝛼,(7) where |𝛼 stands for the 𝑠-, 𝑝-, 𝑑-, states and 𝐶𝑗𝛼 are the expansion coefficients. For instance, the 1𝑠-, 1𝑝-, 1𝑑-, (in the following simply 𝑠-, 𝑝-, and 𝑑-) and 2𝑠-type states are described by the following functions:||𝑠=𝜋2𝑅3𝑗0,|||𝑝𝜋𝑟/𝑅𝑎=32𝜋𝑅3𝑗1𝜇1𝑟/𝑅𝑗0𝜇1𝑥𝑎𝑟,|||𝑑𝑎𝑏=152𝜋𝑅3𝑗2𝜇2𝑟/𝑅𝑗1𝜇2𝑥𝑎𝑥𝑏𝑟2|||𝑑,𝑎𝑏,𝑥2𝑦2=158𝜋𝑅3𝑗2𝜇2𝑟/𝑅𝑗1𝜇2𝑥2𝑦2𝑟2,|||𝑑3𝑧2𝑟2=158𝜋𝑅3𝑗2𝜇2𝑟/𝑅𝑗1𝜇23𝑧2𝑟2𝑟2,||2𝑠=𝜋2𝑅3𝑗0,2𝜋𝑟/𝑅(8) where 𝑗𝑛(𝑥) are the spherical Bessel functions of argument 𝑥, and 𝜇𝑛 are the first roots of 𝑗𝑛(𝑥). Below, we restrict the basis of envelope functions |𝛼 with these states only, and determine the electron and hole spectra. Substitution of Φ𝑗(𝐫) into (4) yields algebraic equations for 𝐶𝑗𝛼:𝐸𝐸𝛼𝐶𝑖𝛼=𝛽||𝐻𝛼1𝑖𝑗+𝐻so𝑖𝑗+𝑉sp𝐫𝛿𝑖𝑗||𝛽𝐶𝑗𝛽,(9) where 𝐸𝛼 denotes the energy of the state |𝛼, 𝛿𝑖𝑗 stands for the Kronecker delta, and the Einstein convention has been used when summing over 𝑗. Considering all operators in the right side of (9)) as perturbations [92], one can solve the equation and find the energies in the valence band.

Restricting the basis of the envelope states |𝛼 by 1𝑠, 2𝑠, 1𝑝, and 1𝑑 functions, and neglecting the spin-orbit interaction, we obtain the following “hierarchy” of energies (for convenience, we write down the hole energies differing in the sign from the electron energies, which are determined with (9))). The ground-state energy𝐸𝑠=2𝜋22𝑚𝑅2(10) is triply degenerate (or sixfold, if spin is taken into account). It is the energy of the 1𝑠 state. Hybridization of the 𝑝-type states yields four different levels with the energies𝐸𝑝1=2𝜇212𝑚𝑅213𝑁+2𝐿2𝑀,𝐸5𝐿+10𝑀𝑝2=2𝜇212𝑚𝑅213𝑁4𝐿+4𝑀,𝐸5𝐿+10𝑀𝑝3=2𝜇212𝑚𝑅21+3𝑁2𝐿+2𝑀,𝐸5𝐿+10𝑀𝑝4=2𝜇212𝑚𝑅21+6𝑁+4𝐿4𝑀,5𝐿+10𝑀(11) where 𝜇1=4.4934. The energies of the 𝑑-𝑑 hybridized states are as follows:𝐸𝑑1=2𝜇222𝑚𝑅216𝑁,𝐸7𝐿+14𝑀𝑑2=2𝜇222𝑚𝑅213𝑁+4𝐿4𝑀,𝐸7𝐿+14𝑀𝑑3=2𝜇222𝑚𝑅21+3𝑁2𝐿+2𝑀,𝐸7𝐿+14𝑀𝑑4=2𝜇222𝑚𝑅21+6𝑁4𝐿+4𝑀,7𝐿+14𝑀(12) where 𝜇2=5.7634. Finally, the energies of the 1𝑑-2𝑠 hybrids are written as𝐸2𝑠𝑑1=22𝑚𝑅24𝜋2+𝜇222+3𝜇2214𝑁+2𝐿2𝑀𝐿+2𝑀1265𝜋𝜇24𝜋2𝜇22𝑁,𝐸𝐿+2𝑀2𝑠𝑑2=2𝜇222𝑚𝑅2,𝐸2𝑠𝑑3=22𝑚𝑅24𝜋2+𝜇222+3𝜇2214𝑁+2𝐿2𝑀+𝐿+2𝑀1265𝜋𝜇24𝜋2𝜇22𝑁.𝐿+2𝑀(13) The contribution of the self-polarization term from (10)–(13) has been omitted since it represents, in fact, a common shift of all the levels by 𝑒2(1/𝜀𝑑1/𝜀𝑠)/2𝑅.

In Figure 5 the energy levels from (10)–(13) have been plotted as functions of the NC radius. The two lowest hole levels have very similar energies. The lowest level, being sixfold degenerate, corresponds to six electron states with an 𝑠-type envelope function. These states have an isotropic macroscopic (averaged over the unit-cell volume) charge distribution inside the dot. The second level corresponds to six electron states with 𝑝-type envelope functions. Correspondingly, the macroscopic charge has an anisotropic distribution for these states (see [92] for details).

In the conduction band the Bloch-state basis consists of six Bloch functions of three different 𝑋-points in the Brillouin zone [99]. We denote three pairs of these functions as |𝑋, |𝑋; |𝑌, |𝑌; |𝑍, |𝑍. Each pair belongs to the twofold degenerate irreducible representation 𝑋1 of the corresponding 𝑋-point. The unprimed Bloch functions have a nonzero value at the lattice site, while the primed functions are zero at those points. The choice of the 𝑋1-functions as the Bloch basis allows us to describe correctly a dispersion-law nonparabolicity that originates mainly from the energy-branch crossing in the 𝑋-points.

The 𝐤𝐩 Hamiltonian operator in the conduction band represents a 6×6 matrix having the form𝐻𝐻=𝑥000𝐻𝑦000𝐻𝑧.(14) Each element of the matrix 𝐻 is a block 2×2 matrix defined by the following expressions [99]:𝐻𝑎=𝑝2𝑎2𝑚𝑙+𝐩2𝑝2𝑎2𝑚𝑡1𝑚𝑡1𝑚0𝑝𝑏𝑝𝑐𝑝+𝑖0𝑚𝑙𝑝𝑎1𝑚𝑡1𝑚0𝑝𝑏𝑝𝑐𝑝𝑖0𝑚𝑙𝑝𝑎𝑝2𝑎2𝑚𝑙+𝐩2𝑝2𝑎2𝑚𝑡.(15) Here 𝑝00.14(2𝜋/𝑎0) is the distance in 𝐩-space from any of the energy minima to the nearest 𝑋-point, 𝑎0 is the lattice constant of silicon. 𝑚𝑙=0.92𝑚0 and 𝑚𝑡=0.19𝑚0 are the “longitudinal” and “transverse” effective masses. The origin of the 𝐸-axis coincides with the 𝑋-point energy. Each of the small indices (𝑎, 𝑏, 𝑐) runs over the values 𝑥, 𝑦, or, 𝑧 and always differs from the others.

Solving (4) for the conduction band is done in the same way as for the valence band, as in (9)–(13) above. One can then obtain the following groups of energies. The ground states, as well as the second excited ones, are the 1𝑠-1𝑝 hybrids. Their energy levels are given by𝐸𝑠𝑝𝑒1=𝐸𝑠+𝐸𝑝2𝐻𝑝𝑝2𝐸𝑝𝐸𝑠2𝐻𝑝𝑝22+𝐻2𝑠𝑝,𝐸𝑠𝑝𝑒2=𝐸𝑠+𝐸𝑝2𝐻𝑝𝑝2+𝐸𝑝𝐸𝑠2𝐻𝑝𝑝22+𝐻2𝑠𝑝,(16) where𝐻𝑠𝑝=𝑖𝑝0𝑚𝑙𝑠|||𝑝𝑧|||𝑝𝑧=2𝜋𝜇1𝑝03𝑚𝑙𝑅𝜇21𝜋2,𝐻𝑝𝑝=1112𝑚𝑡1𝑚𝑙𝑝𝑧|||𝐩23𝑝2𝑧|||𝑝𝑧=2𝜇21𝑚𝑙𝑚𝑡15𝑚𝑡𝑚𝑙𝑅2(17) are 𝑠-𝑝 and 𝑝-𝑝 type matrix elements of anisotropic part 𝐻𝑎, 𝐸𝑠=2𝜋2/2𝑚𝑒𝑅2, and 𝐸𝑝=2𝜇21/2𝑚𝑒𝑅2 are the energies of the 𝑠- and 𝑝-type states of the isotropic unperturbed Hamiltonian operator. The isotropic average electron effective mass is defined by 𝑚𝑒1=(2𝑚𝑡1+𝑚𝑙1)/3(0.26𝑚0)1. The 𝑝-𝑝 hybridization forms states with energies:𝐸𝑝𝑒1,2=𝐸𝑝+𝐻𝑝𝑝𝐻𝑥𝑦.(18) Here 𝐻𝑥𝑦=2𝜇21(𝑚0𝑚𝑡)/5𝑚𝑡𝑚0𝑅2. The next three levels are the result of the 𝑑-𝑑 hybridization:𝐸𝑑𝑒1=2𝜇222𝑚𝑒𝑅2217𝑚𝑙𝑚𝑡2𝑚𝑙+𝑚𝑡67𝑚𝑙𝑚0𝑚0𝑚𝑡2𝑚𝑙+𝑚𝑡,𝐸𝑑𝑒2=2𝜇222𝑚𝑒𝑅241+7𝑚𝑙𝑚𝑡2𝑚𝑙+𝑚𝑡,𝐸𝑑𝑒3=2𝜇222𝑚𝑒𝑅2217𝑚𝑙𝑚𝑡2𝑚𝑙+𝑚𝑡+67𝑚𝑙𝑚0𝑚0𝑚𝑡2𝑚𝑙+𝑚𝑡.(19) Finally, the 2𝑠-𝑑 hybridized energies are obtained numerically:𝐸2𝑠𝑑𝑒1=22𝑚𝑒𝑅24𝜋2+2𝜇2230.477𝜇22,𝐸2𝑠𝑑𝑒2=22𝑚𝑒𝑅24𝜋2+2𝜇2230.455𝜇22,𝐸2𝑠𝑑𝑒3=22𝑚𝑒𝑅24𝜋2+2𝜇223+0.932𝜇22.(20) The results for the single-particle electron spectrum in the conduction band are shown in Figure 6.

The wave function of the ground state turns out to be a superposition of the terms with 𝑠- and 𝑝-type envelope functions (see [92] for details). The “weight” of the 𝑠-type envelope function is greater than that of the 𝑝-type function. Therefore the macroscopic charge for the ground state has a distribution close to isotropic.

From the above calculations, the optical gap of the NC can be obtained as the sum of the lowest energies in both bands and the energy difference between 𝑋- and Γ-points in bulk Δ𝑋Γ=1.24 eV:𝜀𝑔𝑅=Δ𝑋Γ+𝐸𝑠𝑝𝑒𝑙+𝐸𝑠.(21) This energy somewhat overestimates value compared to that obtained with the numerical methods, as indicated by the solid line in Figure 7. This is a consequence of two main factors.

The first factor is that we have assumed infinitely high potential barriers at the interface. The authors of [91, 93, 95] calculated the carrier energies supposing the potential barriers to be finite and equal to approximately 3 eV (or a little more) for electrons and 4-5 eV for holes. Moreover, the discontinuity of the effective mass has been taken into account. In Figure 7 we have plotted the optical gap of silicon nanocrystal versus the confining parameter 1/𝑅 for three models differeing according the nanocrystal’s surroundings. The solid line corresponds to infinitely high potential barriers. Circles represent the finite barriers and constant effective mass throughout the sample. Finally, in case of the mass discontinuity with effective mass outside the dot coinciding with free electron mass, the optical gap has been represented by disks. For the last two cases, the dependence Δ𝜀𝑔(𝑅) becomes weaker than 𝑅2, especially for the case of the discontinuous effective mass where the dependence is closer to 1/𝑅. In the last case the gap dependence is well fitted by the function:𝜀1𝑔𝑅=𝜀𝑔2+𝐷1𝑅2(22) with 𝐷1=4.8 eV2·nm2. In general, one may conclude that the calculated values of the optical gap are reduced due to the barrier finiteness and mass discontinuity, and agree better with those shown in Figure 4.

The second factor, strongly influencing the dot gap, is electron-hole Coulomb interaction. Its contribution to the optical gap will be discussed in the next subsection.

3.1.3. EFA Calculations: Exciton Corrections

The electron-hole interaction reduces the energy of the electron-hole pair. As a result, the exciton optical gap is also reduced. In order to estimate the exciton correction to the total electron-hole single-particle energy one must solve the two-particle problem. This has been done for silicon NCs embedded in an infinitely wide-band-gap material [67, 90, 94]. The exciton corrections were found to be about 0.4 eV [67], 0.25 eV [90], and 0.3 eV [94] for 2-nm-diameter NCs, with the correction decreasing approximately as 𝑅1 with increasing size.

Meanwhile, as has been pointed out in the preceding section, the finiteness of the potential barriers and effective mass discontinuity can appreciably influence the electron and hole energies. Therefore, it is possible to expect the same for the excitonic energies. Ferreyra and Proetto [91] have studied excitonic states for a quantum dot in vacuum. Accordingly, they accepted the barrier heights equal to the electron affinity of the corresponding bulk material for electrons and infinity for holes, and the effective mass being 𝑚0. They found an exciton correction of the order of 0.2 eV for a 2-nm-diameter dot. As the nanocrystal size increases, the magnitude of this correction decreased as 𝑅0.7.

The authors of some of the previous work (e.g., [67, 90, 91]) treated the Coulomb potential energy as 𝑉𝐶=𝑒2/𝜀𝑠𝑟eh similar to bulk silicon. However, in a quantum dot the electron and hole interact not only with each other but also with their “images” due to the difference in the permittivity of the materials inside and outside the dot. These are the so-called polarization fields arising because of the charge polarization at the dot boundary. The interaction with the polarization fields has been taken into account in [94]. Nevertheless, the polarization fields did not lead to any significant corrections to the optical gap. As mentioned above, the correction does not exceed 0.3 eV, which is close to the values obtained by other authors.

Evidently, the calculations carried out with the EFA [67, 90, 91, 94] provide sufficiently good coincidence with both experimental and theoretical data. For comparison, we have depicted in Figure 4 with the solid line the EFA single-particle gap according to (22). The gap values 𝜀𝑔(1)(𝑅) agree well with those computed by PP, TB, and DFT methods. Obviously, a small exciton correction of about 0.3–0.1 eV for 2 to 5 nm diameter crystallites, respectively, does not significantly change the single-particle gap values.

It should be noted, in conclusion, that all the authors employed the perturbation theory considering the electron-hole interaction to be weak compared to the typical quantum confinement energies in the dot. Such an approach is justified if the exciton Bohr radius is significantly larger than the dot radius. This requirement remains valid for small quantum dots with the sizes of about 5-6 nm or less.

3.1.4. Interface States

We have already emphasized in the introduction that radiative transitions between interface states are also often considered as the origin of NIR, or even visible, emission from silicon nanocrystals. In a certain sense, this point of view is alternative to the idea of quantum confinement. Description of the interface states with EFA is difficult because the real potential existing in the vicinity of the nanocrystal surface has in principle a microscopic nature that manifests itself at a distance on the order of a bond length. Meanwhile, the EFA is, in fact, a “macroscopic” method which is not able to distinguish the spatial structure on scales smaller than the size of the primitive cell. Therefore TB and DFT methods are essentially more suitable for this goal.

As a rule, interface states originate from various-type defects existing at the surface of the crystallite. These include, mainly, dangling bonds (Pb centers) which can be neutral, positively, or negatively charged, and different SiO bonds (backbonded and double bonded oxygen) arising at Si/SiO2 interface. Charged dangling bonds corresponding to zero or two electrons in the bond level, respectively, produce deep levels inside the band gap of bulk silicon. Their energies are approximately 0.3 eV below the conduction band minimum for the negatively charged bond, and 0.3 eV above the valence band maximum if the bond is positively charged [100]. Consequently, the effective gap between these levels is about 0.5-0.6 eV [101]. For nanocrystals, the effective gap gradually rises as the nanocrystal size decreases. In particular, for a 5-nm-diameter NC this shift equals approximately 0.6 eV, while for a 2 nm crystallite it is 1.3-1.4 eV [100]. Thus, we see that the dangling bond levels are sensitive to the NC size.

Nanocrystal oxidation has been studied theoretically by TB [8, 102], DFT [103, 104], and GW+BSE [89] methods. The role of the so-called backbonded (Si–O–Si) and double bonded (Si=O) oxygen atoms has been examined for the case where the spherical NC surface is passivated with hydrogen. In the first case (Si–O–Si), oxygen forms two bonds with two different silicon atoms which already have four bonds with other silicon or hydrogen atoms. In the second case (Si=O), oxygen replaces two hydrogen atoms bonded with one silicon atom. It is natural to suppose that the system perturbation in the second case should be more significant because of the stronger distortions of the structure. In fact, it was found that the presence of one Si=O bond at the NC surface provides an optical gap almost independent of the dot size [8]. At the same time, the nanocrystal with a Si–O–Si bond has a gap that gradually decreases with increasing size. Nevertheless, this gap remains smaller, and varies more slowly than the gap of perfect nanocrystal with hydrogen passivated surfaces. The calculations carried out by different methods [89, 102, 103] for oxidized crystallites show an essential difference in the energy of the ground electron-hole transition in crystallites with Si–O–Si and Si=O bonds. All the methods used exhibit the smallest gap for the case of double-bond oxidation, while the crystallites with backbonded oxygen have the greater gap. In turn, the latter is less than the gap of a perfect NC. The authors of [89] have found that only the bridge-type bond (Si–O–Si) can provide a satisfactory agreement with experiments on the photoluminescence in nc-Si/SiO2 system. This fact may be explained by the considerably different oscillator strengths of the ground transitions for NCs with Si=O and Si–O–Si defects. As has been shown by Nishida [102], the oscillator strength for Si–O–Si case is several orders of magnitude greater than that for Si=O case, especially for NCs smaller than 1 nm in diameter. This is precisely the case studied by Luppi et al. [89], where NCs of 0.5 nm and 0.9 nm in diameter were considered. However, according to [102] the difference between the size-quantized energy levels and those of the Si=O and Si–O–Si defects, as well as the difference in the oscillator strengths of the ground electron-hole transitions, are essential only for NC diameters less than 2 to 2.5 nm [8, 102]. For larger sizes, the energies and the transition intensities for perfect and oxidized crystallites almost coincide.

Finally, nonradiative trapping rates are also critical in the light emission properties of silicon NCs. The group of neutral or charged Pb centers at silicon dangling bonds is the most important nonradiative trap in bulk silicon [105], and acts also in the case of Si NCs [58]. In a detailed theoretical investigation, Lannoo et al. [106] found that the carrier trapping at neutral Pb centers is nonradiative, whereas capture at charged centers can lead to photon emission at energies smaller than the bandgap of bulk silicon. The rate for carrier trapping at a neutral dangling bond defect at the Si–SiO2 interface in a NC with a single Pb center was established theoretically as a function of the NC size [106]:𝑤𝑛𝑟=𝜎0𝑣𝑉1(𝐸2𝜋0𝜔2+𝑧2)0.25×exp[𝑆coth𝜔+𝐸2𝑘𝑇0𝐸2𝑘𝑇+(0𝜔2+𝑧2)0.5𝐸0𝐸𝜔𝑎sinh0𝜔𝑧].(23) In (23), 𝜎0 is the capture cross section given by 𝑐0/𝑣 where 𝑐0 is the capture coefficient and 𝑣 is the thermal velocity equal to 8𝑘𝑡/𝜋𝑚, 𝐸0 is the ionization energy of the defect (approximately equal to the Franck-Condon energy plus the carrier confinement energy), 𝑆 is the Huang-Rhys factor (𝑆=15) [107], 𝜔 is the average phonon energy, 𝑉 is the nanocrystal volume, and 𝑧=𝑆/[sinh(𝜔/2𝑘𝑇)]. The capture rate increases strongly as a function of NC size due to the decreasing defect ionization energy, until eventually the volume term begins to dominate the trapping rate, which then begins to decrease again. Since, the calculated rate is for a single isolated defect, in practice the rate should be modified for the likelihood of multiple defects on larger clusters, as discussed in Section 4.

3.2. Shallow Impurities in Silicon Nanocrystals

There is one more type of defect which can play an important role in the optical and transport properties of bulk semiconductors and their low-dimensional counterparts (in particular, silicon). These are impurity centers, of which some experimental results were discussed previously. Let us now consider NCs doped with shallow impurities and discuss impurity states in silicon quantum dots.

Various aspects of the problem have been explored. One of them—the central-cell effect [108111]—causes splitting of the sixfold degenerate ground energy level in bulk silicon into a singlet, doublet, and triplet with a typical energy splitting of about 10–20 meV for various donors [112]. According to a number of theoretical studies [73, 113119], in quantum dots doped with donors or acceptors the central-cell potential strengthens the level splitting compared to (i) bulk silicon, and (ii) undoped dots [66, 85, 120, 121]. Investigations of the spatial charge distribution [122, 123], the formation of impurity centers inside silicon nanocrystals from the energy point of view [114117, 124, 125], intervalley scattering [126], hyperfine splitting and optical gap effects [127], and the screening of the point-charge field [123, 128136] have been performed. Below, we present EFA calculations of the electronic structure of both donor- and acceptor-doped silicon quantum dots and examine the dependence on the impurity position inside the dot within the framework of the hydrogenic-impurity model. Also, we discuss the effect of the central-cell potential in screening the point charge in NCs.

3.2.1. Screening in Quantum Dots

In order to solve accurately the Schrodinger-like equation for the envelope functions in the dot in the presence of an impurity ion, we need the correct expression for the Coulomb potential energy 𝑉𝐶. Its determination in the quantum dot is not a simple task because of the complicated character of screening in NCs. Some authors described the screening properties in an NC with a modified dielectric constant 𝜀(𝑅) that depends on the dot radius [137139]. Such a simple model, indeed, gives a moderate increase of the dielectric properties due to the finite size of the crystallite but it does not reflect correctly the local structure of the electric field.

The correct description of such a structure requires first-principles calculations [123, 129136]. Nevertheless, the microscopic picture permits a clear qualitative macroscopic interpretation using the local dielectric function 𝜀(𝑟) depending on electron position vector, 𝑟. Taking into account the short-range field and charge polarization at the dot boundary due to different static dielectric constants 𝜀𝑠 and 𝜀𝑑 of an NC and its surroundings, one can in principle obtain 𝜀(𝑟) (Figure 8) [73].

An overall decrease of 𝜀(𝑟) in a nanocrystal takes place compared to the bulk value of 𝜀𝑠=12. As was pointed out earlier [123, 128, 129], this decrease is mainly due to the polarization charges at the dot boundary. The monotonic decrease of 𝜀(𝑟) extends to the nanocrystal boundary, where the value of the dielectric function becomes equal to 𝜀𝑑. At the same time, the sharp reduction of 𝜀(𝑟) toward unity occurring at small 𝑟 is exclusively due to the short-range central-cell potential.

3.2.2. Donor States in Silicon Crystallite

First, we consider electronic states in a silicon NC with a hydrogenlike donor placed in some arbitrary position inside the dot. The total electron potential energy described by (2) transforms into 𝑈(𝐫,𝐡)=𝑈0(𝑟)+𝑉𝐶 where 𝐡 stands for the impurity position vector, and 𝑉𝐶 is a Coulomb potential energy, including not only the direct electron-ion interaction but also the interaction with polarized charges arising at the dot boundary. Solving the eigenstate and eigenvalue problem for the conduction band [140] yields the ground-state energy splitting shown in Figure 9 at 𝑛𝑥=0.8, 𝑛𝑦=0.5, and 𝑛𝑧=0.33, where 𝐧=𝐡/, and stands for the absolute value of 𝐡.

In Figure 10, the envelope-function correction ΔΦ is shown for the ground state due to the hydrogen-like donor, as a function of the angles 𝜃 and 𝜑 on the spherical surface 𝑟= for the former values of 𝑛𝑎. Angles 𝜃 and 𝜑 are the spherical coordinates. The maximal values of ΔΦ (light areas in the figure) are located around the donor site marked with the cross. This is a natural result in that the electron density shifts towards the donor site.

3.2.3. Acceptor States in Silicon Crystallite

It is interesting to compare the fine structure of the spectrum for donor- and acceptor-doped dots. Obviously, the fine structure in the valence band occurs because of the spin-orbit interaction and the asymmetry of the Coulomb field inside the dot, when the acceptor occupies some noncentral position. The energy splitting has been presented in Figure 11. Similarly to the donor case, the splitting turns out to be large compared to the magnitudes of chemical shifts in bulk silicon [141144].

The Coulomb and spin-orbit interactions remove the sixfold degeneracy of both lowest levels shown in Figure 5, keeping only the double spin degeneracy of each the levels. Furthermore the splitting increases with decreasing the nanocrystal size. This is evidence in favor of a quantum confinement effect, as was the case in the donor-doped dots.

The spatial distribution of the hole density in B-doped crystallite for all the six doublets is shown in Figure 12 for /𝑅=0.1 (upper panels) and /𝑅=0.46 (lower panels). We have plotted the average of the squared absolute value of the total wave function over the unit cell. In this case the Bloch-function oscillations do not appear in the electron density, which reduces to the “density of envelope function” 𝜌env(𝐫)=6𝑗=1|Φ𝑗(𝐫)|2, where Φ𝑗(𝐫) is the 𝑗th element of the 6D envelope function vector |Φ as before.

As the calculations show, for all the six doublets the envelope-function density has an axial symmetry with respect to the line drawn through the NC center and the acceptor. Therefore 𝜌env(𝐫) has the same distribution in any dot cross-section to which the acceptor position-vector belongs. Figure 12 shows such a central cross-section of the electron density averaged over the unit cell for all the doublets. Brighter areas in the density plots correspond to higher values of 𝜌env(𝐫). The circle represents the nanocrystal boundary, and the bold point situated at the vertical axis indicates the acceptor location.

It should be noted that the general trends in the location of the electron density under the action of the acceptor electric field are, in fact, similar for all the doublets in both cases /𝑅=0.1 and /𝑅=0.46. In particular, the ground-state electron density shifts to the acceptor site, while for the excited states, as the energy of the state increases and the electron density gradually moves into the areas free of the acceptor.

3.3. Interband Recombination in Silicon Nanocrystals

The PL intensity and quantum efficiency is defined by the relative contribution 𝜏𝑅1/𝜏1PL of a radiative recombination channel in an interband transition. In turn, the total decay rate is defined by both radiative and nonradiative recombination rates: 1/𝜏PL=1/𝜏𝑅+1/𝜏NR. Consequently, in order to discuss possible means of increasing the quantum efficiency of silicon crystallites, we should understand what factors influence the radiative and nonradiative rates. For this purpose in this section we analyze theoretically both radiative and nonradiative recombination processes.

3.3.1. Radiative Recombination Times in Undoped Si NCs

In bulk silicon, no-phonon radiative transitions between the conduction-band Δ-point and the valence-band Γ-point are forbidden because of indirect band-gap of silicon. In this case, only phonon-assisted transitions may take place, as shown in Figure 13. In a quantum dot no-phonon emission becomes possible but remains a low probability process for reasons discussed below. Meanwhile, the phonon-assisted radiative transitions have much greater rate. Here, we calculate both no-phonon and phonon-assisted radiative lifetimes in silicon nanocrystal.

We first calculate the rate 𝜏01 of the no-phonon radiative recombination. Using Fermi's golden rule in the first-order perturbation theory, one can write the decay rate for the transition between initial (𝐼) and final (𝐹) states in the conduction and valence bands, respectively, in the following form:𝜏10,𝐼𝐹=2𝜋𝐐,𝜎|||𝑊𝐼𝐹|||2𝛿𝜀𝑔𝑅𝜔𝜎(𝐐).(24) Here, 𝜔𝜎(𝐐) and 𝐐 stand for the photon frequency and wave vector. The Dirac delta-function reflects the energy conservation law for the electron transition. The operator 𝑊 describes the electron-photon interaction, and has the form:𝑊=𝐤,𝜎2𝜋𝑒2𝜅𝜀𝑠;𝜀𝑑𝑚20𝜔𝜎(𝐐)𝑉0𝜀𝑠𝑐𝐐𝜎+𝑐+𝐐𝜎𝐞𝐐𝜎𝐩,(25) where 𝐩=𝑖 is the electron momentum operator, 𝐞𝐐𝜎 is the polarization vector, the operator 𝑐𝐐𝜎 annihilates and 𝑐+𝐐𝜎 creates a photon with the wave vector 𝐐 and polarization 𝜎, and 𝑉0 stands for the volume of the electromagnetic resonator. The function 𝜅(𝜀𝑠;𝜀𝑑) is written as [145]𝜅𝜀𝑠;𝜀𝑑=𝜀𝑑𝜀𝑠3𝜀𝑑2𝜀𝑑+𝜀𝑠2.(26) This function represents the correction factor in the field magnitude due to the replacement of homogeneous media with bulk permittivity 𝜀𝑠 by a spherical silicon nanocrystal surrounded by silicon dioxide with a dielectric constant 𝜀𝑑.

The initial state in (24) corresponds to an electron-hole pair being in its ground state, and an ensemble of photons whose distribution over 𝐐 is described by the Bose-Einstein statistics. In the final state, the valence band is completely occupied and the conduction band is empty. The number of photons in the final state always increases by one.

As has already been mentioned (see Section 3.1.2), the ground spinless electron state in the conduction band is sixfold degenerate. However, TB [85, 120], PP [66], and DFT [121] calculations revealed a weak splitting of the ground state into the singlet, doublet, and triplet states due to the tetrahedral symmetry of the spherical silicon nanocrystal. The authors pointed out that as the NC size varies, 𝐴1-, 𝐸-, or 𝑇2-type states alternately become the ground level. Correspondingly, in order to take into account tetrahedral symmetry of the quantum dot, we build the six ground states in accordance with the symmetric transformations of the irreducible representations 𝐴1 (singlet), 𝐸 (doublet), and 𝑇2 (triplet) of the tetrahedral group. These functions are defined as follows [73]:Ψ𝑆𝜆|||𝐴=cos1||𝜆|||𝑋𝑠+sin|||𝑝𝑥+|||𝑌|||𝑝𝑦+|||𝑍|||𝑝𝑧3,Ψ𝐷(1)𝜆|||𝐸=cos(1)||𝜆|||𝑋𝑠+sin|||𝑝𝑥|||𝑌|||𝑝𝑦2,Ψ𝐷(2)𝜆|||𝐸=cos(2)||𝜆|||𝑋𝑠+sin|||𝑝𝑥+|||𝑌|||𝑝𝑦|||𝑍2|||𝑝𝑧6,Ψ𝑇(1)𝜆|||𝑋=cos||𝜆𝑠sin2|||𝐴1+3|||𝐸(1)+|||𝐸(2)6|||𝑝𝑥,Ψ𝑇(2)𝜆|||𝑌=cos||𝜆𝑠sin2|||𝐴13|||𝐸(1)+|||𝐸(2)6|||𝑝𝑦,Ψ𝑇(3)𝜆|||𝑍=cos||𝜆|||𝐴𝑠sin12|||𝐸(2)3|||𝑝𝑧.(27) Here, the parameter 𝜆 is defined by the relationships [92]=𝐸cos2𝜆𝑝𝐸𝑠2𝐻𝑝𝑝𝐸𝑝𝐸𝑠2𝐻𝑝𝑝2+4𝐻2𝑠𝑝,=sin2𝜆2𝐻𝑠𝑝𝐸𝑝𝐸𝑠2𝐻𝑝𝑝2+4𝐻2𝑠𝑝,(28) where vectors |𝑠 and |𝑝𝑎 stand for the 𝑠- and 𝑝-type envelope functions as before. The Bloch states |𝐴1=(|𝑋+|𝑌+|𝑍)/3, |𝐸(1)=(|𝑋|𝑌)/2, |𝐸(2)=(|𝑋+|𝑌2|𝑍)/6, and |𝑇2=|𝑋, |𝑌, or |𝑍 belong to the representations 𝐴1, 𝐸, and 𝑇2 of the tetrahedral group. In what follows, when computing the interband matrix elements of 𝑊 we will choose one of the states given by (27) as the initial one. For convenience, one can rewrite the wave functions of the initial states in a general form:Ψ𝐼=𝛼,𝑗𝐵𝐼𝛼𝑗𝐹𝛼𝐫𝜓𝑗𝐫,(29) where the index 𝛼 indicates the “values” of the 𝑠, 𝑝𝑥, 𝑝𝑦, and 𝑝𝑧, functions, and 𝜓𝑗(𝐫) stands for the six Bloch functions |𝑋, |𝑌, |𝑍, |𝑋, |𝑌, |𝑍 of the irreducible representation 𝑋1. The 𝑠- and 𝑝-type envelope functions are denoted as 𝐹(𝑟), and 𝐹𝑥(𝐫), 𝐹𝑦(𝐫), 𝐹𝑧(𝐫), respectively. Selecting expansion coefficients 𝐵𝐼𝛼𝑗, we can construct any of the states defined by (27).

The final state is one of the three degenerate states (neglecting the spin-orbit coupling) in the valence band [92] with the 𝑠-type envelope function:Ψ𝐹𝑟𝜓=𝐹𝐹𝐫,(30) where 𝜓𝐹(𝐫) stands for the Bloch state basis functions |𝑌𝑍, |𝑋𝑍, or |𝑋𝑌 of the irreducible representation Γ25.

After some algebra, one can obtain recombination rate 𝜏10,𝐼𝐹 in the form:𝜏10,𝐼𝐹=4𝑒2𝜀𝑠𝜅𝜀𝑠;𝜀𝑑𝜀𝑔𝑅3𝑚202𝑐3|||𝐩𝐼𝐹|||2,(31) where 𝐩𝐼𝐹=Ψ𝐼|𝐩|Ψ𝐹. In order to compute the momentum matrix element, we expand the Bloch functions 𝜓𝑗(𝐫) and 𝜓𝐹(𝐫) over the reciprocal lattice vectors 𝐠𝐧 as follows:𝜓𝑗𝐫=𝐦𝐶𝑗𝐦𝐠exp2𝜋𝑖𝐦+𝐞𝑗𝐫/𝑎0,𝜓𝐹𝐫=𝐧𝑉𝐹𝐧exp2𝜋𝑖𝐠𝐧𝐫/𝑎0,(32) where 𝑎0 is the lattice constant of silicon, and unit vectors 𝐞𝑗 define the 𝑋-point. We set 𝐞1,4=𝐞𝑥, 𝐞2,5=𝐞𝑦, and 𝐞3,6=𝐞𝑧. Expansion coefficients 𝐶𝑗𝐦 and 𝑉𝐹𝐧 are determined with the local pseudopotential method. They satisfy the normalizing condition 𝐧|𝐶𝑗𝐧|2=𝐧|𝑉𝐹𝐧|2=1. As a result, the momentum matrix element is given by𝐩𝐼𝐹=4𝜋𝑎0𝐛𝐼𝐹𝑅𝑎0𝑅4,(33) where 𝐛𝐼𝐹(𝑅)=𝛼,𝑗𝐵𝐼𝛼𝑗𝐉𝛼(𝑅) with 𝐉𝛼(𝑅) being oscillating vector functions [146] of the dot radius. For 𝑠- and 𝑝-type subscripts, respectively, these functions are𝐉𝑅=𝐧,𝐦𝐶𝑗𝐦𝑉𝐹𝐧𝐝𝐧+𝐦𝑑4𝐧𝐦cos2𝜋𝑑𝐧𝐦𝑅𝑎0,𝐉𝑎𝑅=𝑖3𝜇1𝜋𝐧,𝐦𝐶𝑗𝐦𝑉𝐹𝐧𝐝𝐧+𝐦𝑑𝑎𝐧𝐦𝑑5𝐧𝐦sin2𝜋𝑑𝐧𝐦𝑅𝑎0.(34) Here, 𝑎=𝑥,𝑦,𝑧, 𝐝𝐧±𝐦=𝐠𝐧±𝐠𝐦±𝐞𝑗, and 𝑑(𝑎)𝐧𝐦 represents an 𝑎th component of 𝐝𝐧𝐦. The nonzero value of 𝐩𝐼𝐹 is exclusively due to the Heisenberg uncertainty relations which cause the wide distribution of the wave function in 𝐩-space.

Averaging (31) over all possible initial and final states yields𝜏01=𝑒2𝜀𝑠𝜅𝜀𝑠;𝜀𝑑𝜀𝑔𝑅12𝜋2𝑚20𝑐3𝑎20𝑎0𝑅8|||𝐛𝐼𝐹𝑅|||2.(35) In Figure 14 we show the no-phonon recombination rate 𝜏01 as function of the dot size. Although the no-phonon decay rate has a nonzero value, its magnitude remains small (less than 103 s−1) within the range of the NC sizes shown in the figure.

Next, we calculate the phonon-assisted recombination rate [95, 147149]. The total phonon-assisted rate of the radiative electron-hole recombination is determined in the second-order perturbation theory as𝜏1𝑅,𝐼𝐹=2𝜋𝐐,𝜎𝐪,|𝑎𝑊𝐹𝐴𝑈𝐴𝐼+𝑈𝐹𝐴𝑊𝐴𝐼𝜀𝐼𝜀𝐴|2×𝛿𝜀𝑔𝑅𝜔𝜎(𝐐)𝜈𝜀(𝐪)+𝛿𝑔𝑅𝜔𝜎(𝐐)+𝜈.(𝐪)(36) Here the matrix elements of the electron-photon (𝑊) and electron-phonon (𝑈) interaction operators are calculated between the initial 𝐼, the final 𝐹, and an intermediate state 𝐴; 𝜀𝐴 and 𝜀𝐼 are the total energies of the intermediate and initial states, respectively, including not only the energies of the electrons (or holes) but also the energies of the photons and phonons. The phonon frequency of th polarization is denoted as 𝜈(𝐪) with 𝐪 being the dimensionless phonon wave vector taken in units of 2𝜋/𝑎0.

The electron-phonon interaction operator is treated within the framework of the rigid-ion model and is given by𝑈=𝐪,𝐧,𝑠2𝑀𝑁𝜈(𝐪)𝑉𝐧𝑠×𝐞𝐪𝑠exp𝑖2𝜋𝐪𝐑𝐧𝑎0𝑏𝐪+𝐞𝐪𝑠exp𝑖2𝜋𝐪𝐑𝐧𝑎0𝑏+𝐪.(37) Here, 𝑁 is the number of primitive cells in the crystal, 𝑀 is the mass of a silicon atom, 𝑉𝐧𝑠=𝑉at(𝐫𝐑𝐧𝝉𝑠) is the atomic potential, where 𝐑𝐧 stands for the position vector of the 𝐧th unit cell, and 𝝉𝑠 represents the position vectors of two atoms within the unit cell: 𝝉1=0, and 𝝉2=(1,1,1)×𝑎0/4, 𝑏+𝐪 and 𝑏𝐪 are the phonon creation and annihilation operators, and the phonon polarization vectors are denoted as 𝐞𝐪𝑠. Making the Fourier transformation of the atomic potential: 𝑉at(𝐫)=𝑁1𝐩𝑉𝐩exp{𝑖𝐩𝐫}, one can rewrite the operator 𝑈 in the form (see also [150, 151]):𝑈=2𝜋𝑖𝑎0𝐪,,𝑠,𝐦2𝑀𝑁𝜈(𝐪)𝐪+𝐠𝐦×𝐞𝐪𝑠𝑉𝐪+𝐠𝐦exp𝑖2𝜋𝐪+𝐠𝐦𝐫𝝉𝑠𝑎0𝑏+𝐪𝐞𝐪𝑠𝑉𝐪+𝐠𝐦exp𝑖2𝜋𝐪+𝐠𝐦𝐫𝝉𝑠𝑎0𝑏𝐪.(38) Calculations of the electron-photon and electron-phonon matrix elements yield𝜏1𝑅,𝐼𝐹=4𝜋22𝑒2𝜀𝑠𝜅𝜀𝑠;𝜀𝑑𝜀𝑔𝑅3𝑀𝑚20𝑐3𝑎40𝑎0𝑅3×𝑃𝐼𝐹𝜈coth𝜈,2𝑘𝑇(39) where𝑃𝐼𝐹=𝐵𝐼𝛼𝑚(𝐵𝐼𝛽𝑛)𝑥𝐹𝑚𝑥𝐹𝑛2𝜋2𝑅3×𝑑𝐫𝐹𝛼𝐫𝐹𝛽𝐫𝐹2𝑟,𝑥𝐹𝑚=𝐴𝜓𝐹||𝑤||𝜓𝐴𝜓𝐴||𝑢||𝜓𝑚+𝜓𝐹||𝑢||𝜓𝐴𝜓𝐴||𝑤||𝜓𝑚𝜀𝐼𝜀𝐴.(40)𝜓𝐴 represents the Bloch functions of the intermediate states and the operators 𝑤 and 𝑢 are directly proportional to 𝑊 and 𝑈: 𝑊=𝑤2𝜋𝑎02𝜋𝑒2𝜅𝜀𝑠;𝜀𝑑𝑚20𝜔(𝐐)𝑉0𝜀𝑠,𝑈=𝑢2𝜋𝑎02𝑀𝑁𝜈.(𝐪)(41) As has been shown earlier, the electron subsystem in silicon nanocrystals interacts more efficiently with TO and LO phonons [147, 148] similarly to bulk silicon [150, 151]. The contribution of TA phonons is negligibly small. Therefore, in the sum over the phonon modes in (39), only TO and LO modes may be taken into account.

For convenience, we have again averaged the calculated rates over all the possible initial and final degenerate states and obtained the recombination rate in the following form:𝜏𝑅1=𝜋22𝑒2𝜀𝑔𝑅𝜀𝑠𝜅𝜀𝑠,𝜀𝑑9𝑀𝑚20𝑐3𝑎4𝑎0𝑅3×7.527coth𝜈LO/2𝑘𝑇𝜈LO+32.768coth𝜈TO/2𝑘𝑇𝜈TO×4cos2.𝜆+9.26(42) Here: 𝜈LO=0.051eV and 𝜈TO=0.0576eV represent the energies of LO and TO phonons and the numbers 7.527, 32.768 come from numerical integration of electron-phonon matrix elements. The size dependence of the recombination rate has been presented in Figure 14. The dependence of 𝜏𝑅1 on the dot radius is close to 𝑅3. This is much slower than the rate obtained for no-phonon transitions (𝑅8).

The results presented here agree well both qualitatively and quantitatively with those obtained by other authors [95, 147, 148] for similar systems. In particular, according to their results, the recombination rates for NCs from 2 to 6 nm in diameter do not exceed 105 s−1 and are never less than 5×102 s−1. It should be noted, however, that the TB calculations [148] produce a considerably slower decrease of the radiative decay rate than the EFA (as developed here and in [95, 147, 149]) predicts.

Finally, concluding this subsection, we would like to discuss briefly the relation between our single-particle calculations and the singlet-triplet model which gives a reduction of the radiative lifetime with increasing temperature, as in (1). It is, of course, not applicable in the frame of a single-particle treatment to obtain the state splitting caused by the exchange electron-hole interaction. Therefore, the temperature dependence of 𝜏𝑅1 (42)) is defined exclusively by the phonons. As a result, no sharp rise of the decay rate appears in our model as the temperature increases up to the room values. At temperatures higher than about 150–200 K, both models nevertheless should give similar results because the upper excitonic states become highly populated. Thus, the single-particle model remains valid for 𝑇>200 K. As well, it should be noted that the separation of dark and bright excitons is not always possible. The spin-orbit coupling transforms forbidden exciton transitions into suppressed ones due to the entanglement of single-particle states with different spin projections. The entanglement can be crucially strengthened by imperfections of the crystallite structure such as shape nonsphericity or, especially, by various point defects such as Pb centers, SiO bonds, and impurities. In the presence of the latter effects, all the transitions can become allowed and consequently the single-particle model is acceptable.

3.3.2. Radiative Recombination in Doped Dots

Experimental evidence of shallow impurity effects in the PL of silicon NCs (see Section 2.3) has no corresponding rigorous theoretical grounds. Usually, when discussing various changes in optical properties of the NCs due to doping, the nonradiative recombination channel is considered more closely. For instance, doping may passivate dangling bonds, and, as a consequence, essentially reduce the effectiveness of the nonradiative channel. On the other hand, a large number of donors or acceptors in the nanocrystal leads to the appearance of extra carriers that can recombine through the Auger process. Thus, doping is an important factor influencing the nonradiative recombination (a more detailed description of the nonradiative recombination mechanisms can be found in the next subsection). The nonradiative lifetime 𝜏NR can be varied by doping, which leads to the change of the emission efficiency.

The impact of doping on the radiative channel has not yet reported. However, it is interesting to examine also the possibility of increasing the quantum yield via enhancement of the radiative rate 𝜏𝑅1. This is exclusively a quantum effect that can be explained by the reconstruction of the wave functions in different symmetry (or asymmetry), and changing selection rules for the interband transitions. We would like to touch upon this question in the present subsection and compute the radiative recombination rates in a silicon NC doped with a hydrogen-like donor and acceptor.

In order to compute the decay rates for the doped dot we follow (36). The main reason for the rate changes compared to the case of the undoped dot lies in some corrections to the electronic wave functions of both initial and final states (27), (30) due to the existence of a shallow impurity in the dot. In the following, we discuss two different cases.

The first case corresponds to an off-center position for the impurity inside the dot. In the valence band such an impurity location creates a field asymmetry in the system, which entangles the three spinless 𝑠-type states (30)) with the three spinless 𝑝-type states [92]Ψ𝑝1=𝐹𝑥𝐫||𝑋𝑍𝐹𝑦𝐫||𝑌𝑍2,Ψ𝑝2=𝐹𝑧𝐫||𝑌𝑍𝐹𝑥𝐫||𝑋𝑌2,Ψ𝑝3=𝐹𝑦𝐫||𝑋𝑌𝐹𝑧𝐫||𝑋𝑍2.(43) These states (30) and (43) correspond to the closely spaced lowest triply degenerate energy levels shown in Figure 5. As a result, the ground hole state becomes a superposition of the states with 𝑠- and 𝑝-type envelope functions. In the conduction band, the ground electron states are the 𝑠-𝑝 combinations even without impurities, as seen from (27). Embedding the impurity into some arbitrary position within the NC raises the number of the 𝑝-type envelope functions in the ground states as well as in the valence band. However, 𝑃𝐼𝐹() becomes smaller if the integrated envelope functions are 𝑝-type. Consequently, the squared transition matrix element, proportional to 𝑃𝐼𝐹(), and the rate itself decrease relative to the case of an undoped NC, as shown in Figure 15 for two different displacements of the impurity from the dot center.

The ratio /𝑅=0.46 corresponds to the greatest possible mix of 𝑠- and 𝑝-type states with maximal weight of the 𝑝-type envelope functions in the ground states both in the conduction and valence bands. The strong mixing results in maximal decrease of the decay rate independently of the impurity type. When the impurity is situated near the nanocrystal boundary, /𝑅=0.9, the weight of the 𝑝-type envelope functions reduces in the ground state of both bands. As a consequence, the rate rises with respect to the case /𝑅=0.46, but remains smaller than that for the undoped NC.

The second case we consider is the case of the central-located impurity. If =0, no 𝑠-𝑝 entanglement takes place in the valence band. The functions with the 𝑠-type envelope (30)) describe the ground states of both donor- and acceptor-doped NCs.

The situation is different for the conduction band. Formally, the wave functions of the six spinless ground states are described by (27) as before. However, for doped dots, parameter 𝜆 has a form which differs from that given by (28). In the presence of an impurity center, 𝜆 is defined as:=𝐸cos2𝜆𝑝𝐸𝑠+𝛿𝑉2𝐻𝑝𝑝𝐸𝑝𝐸𝑠+𝛿𝑉2𝐻𝑝𝑝2+4𝐻2𝑠𝑝,=sin2𝜆2𝐻𝑠𝑝𝐸𝑝𝐸𝑠+𝛿𝑉2𝐻𝑝𝑝2+4𝐻2𝑠𝑝,(44) where 𝛿𝑉=𝑉𝑝𝑝𝑉𝑠𝑠 is the difference between the 𝑝-𝑝 and 𝑠-𝑠 type matrix elements of the Coulomb potential energy 𝑉𝐶. The appearance of this term in the definition of 𝜆 is explained by the shift of the energy levels of the 𝑠- and 𝑝-states 𝐸𝑠 and 𝐸𝑝, respectively, due to the Coulomb interaction. The relative shift 𝛿𝑉 turns out to be positive for donors, and negative for acceptors. Consequently, cos(𝜆) increases and sin(𝜆) decreases for donors, and vice versa for acceptors.

This implies that the weight of the 𝑠-type envelope function in the electron ground state increases for the donor-doped dot and shrinks for the dot doped with acceptor, with respect to its value in the undoped dot. Hence, the coefficient 𝑃𝐼𝐹() for the donor-doped dot is greater than that for the undoped dot. In turn, the latter exceeds 𝑃𝐼𝐹() for the acceptor-doped dot. Precisely the same relationships take place for the recombination rates, as shown in Figure 15.

Note in conclusion, that the central-symmetric case =0 qualitatively corresponds to some real situation, where the NC has been highly doped with donors or acceptors. It is natural to assume a homogeneous impurity distribution within the nanocrystal in this case. Obviously, such a homogeneous density of charge induces a spherically symmetric electric field that cannot mix the 𝑠- and 𝑝-type states. Moreover, the relative Coulomb shift 𝛿𝑉 remains positive for donors and negative for acceptors, as in the above case of a single impurity center in the dot [73]. Therefore, one can expect similar behavior of 𝜏𝑅1 for highly doped NCs.

3.3.3. Nonradiative Recombination

In parallel with radiative transitions, nonradiative processes occur in silicon nanocrystals. As a rule, nonradiative processes play a problematic role in silicon crystallites, due to their high recombination rates which substantially exceed the rates of the radiative transitions. Among the possible nonradiative processes, capture on dangling bonds and Auger recombination are the most efficient. Both processes are phonon-assisted, possibly even multi-phonon. The first one takes place both at low and high pumping levels, while Auger recombination becomes possible exclusively at high excitation power when more than one excited electron (hole) exist in the conduction (valence) band of the nanocrystal.

The capture on a neutral dangling bond has been studied earlier by TB method [100, 106]. It was shown that excited electron-hole pairs can recombine on a neutral dangling bond in two stages. First, the electron (hole) is trapped on a neutral Pb center making it negatively (positively) charged. Then, the hole (electron) is trapped on the charged dangling bond. According to estimations [100], the capture rate on the neutral dangling bond for holes is always much greater than that for electrons. In turn, the electron capture rate at room temperature increases sharply as the nanocrystal gap decreases from 2.7 eV to 1.5 eV (see (23)). Within the same range of energy gaps, the hole capture rate also increases but less sharply. We should emphasize the presence of two main features of the capture process. First, there is a very stron dependence of the capture rates (even for holes) on the energy gap, that is, the NC size. The second feature is the overall high capture rates, which are much greater than the radiative decay rates. Evidently, the presence of a dangling bond in silicon crystallite drastically reduces the quantum efficiency of photon generation. Therefore the passivation of dangling bonds is required for achievement of relatively efficient luminescence.

Let us now touch briefly on the role of SiO defects in the photoluminescence. The calculations of Nishida [102] showed that the oscillator strength is almost independent of the existence of Si–O–Si bonds. In contrast, the Si=O bonds reduce the oscillator strength by several orders of magnitude. Recently, however, Sa'ar et al. [64] demonstrated experimentally that both back-bonded and double-bonded oxygen atoms can enhance the PL intensity due to formation of coupled states of electrons (holes) and SiO vibrations. Such vibrations induce electric polarization fields generating long-lived coupled states named “vibrons”. Usually, the surface vibrations scatter the carriers; however, vibrons, being coupled states, do not scatter the carriers and thereby exclude one of possible nonradiative channels from the recombination process.

For highly excited nanocrystals, Auger recombination may take place. It is sufficient to have only three particles (two electrons and a hole, or vice versa) for this process to occur. It is well known that Auger recombination has high efficiency in bulk silicon. Therefore, it is natural to expect a similar behavior in silicon crystallites. Calculations of the Auger-recombination rate in crystallites carried out with TB [58, 152] and EFA [153] techniques confirm, indeed, the fast character of this mechanism. In particular, the characteristic times have been found to be of the order 0.1 to 10 ns−1 and gradually rise with decreasing nanocrystal size. Thus, Auger recombination, as well as dangling bonds, can strongly quench the luminescence in silicon crystallites, as confirmed experimentally [152].

4. Ensembles of Quantum Dots

4.1. Size Distributions

As discussed in detail in the previous sections, the recombination energy for Si NCs is dependent on the NC size. In order to achieve optical gain and stimulated emission, therefore, a narrow size distribution would be beneficial in order to narrow the gain profile. In ensembles of silicon nanocrystals grown by various thin film methods or by ion implantation, it is so far not possible to obtain the narrow size distributions that can be achieved with chemical methods, where post-synthesis size selection methods are routinely employed [154]. However, there is hope for better size selectivity, as will be discussed toward the end of this section.

The nucleation and growth behavior of nanocrystals in ion implanted systems, in particular, has been extensively characterized [155, 156]. The lognormal size distribution results in general when multiple microscopic processes govern the nucleation and growth kinetics, as is typical in these systems. The lognormal distribution is found even in cases where the “memory” of the specific initial conditions is lost, for example, after annealing of ion implanted samples [156]. In such cases, the nanocrystal size distribution is governed by the standard lognormal curve:1𝑃(𝑥)=𝑆𝑥𝑥2𝜋exp[ln𝑀22𝑆2],(45) where the mean and variance can be related to the 𝑀 and 𝑆 parameters by 𝜇=exp(𝑀+𝑆2/2) and var=exp(𝑆2+2𝑀)[exp(𝑆2)1].

For thin films also, a lognormal distribution is predicted on the basis of random nucleation and growth in a homogeneous medium, regardless of the growth temperature, specific method of crystallization, and the mean grain size [157, 158], which is consistent with previous theories on cluster growth in metals and ceramics [159]. Lognormal distributions in the case of silicon grains crystallizing in amorphous silicon have been reported [157], and the lognormal size distribution has also been observed in silicon NCs over the range of annealing temperatures from 400 to 1100°C [47, 160]. Figure 16 shows the size distributions obtained from TEM data for two samples available in the literature. In Figure 16(a), the data from Vinciguerra et al. [161] are reproduced along with a Gaussian least squares fit (𝜇radius=1.68 nm, 𝜎=0.65 nm) which is similar to the fit in the original publication (quoted as 𝜇=1.7±0.6 nm), and a lognormal fit with 𝑆=0.27 and 𝑀=0.54). The data look decidedly lognormal with a skew toward large radii; the lognormal fitting matches the data better on both the high and low tails of the distribution, has a better correlation (𝑅2=0.94 versus 0.86), and is naturally equal to zero for negative radii. Figure 16(b) shows the original data from Glover and Meldrum [162], along with a lognormal fit with 𝑆=0.29, 𝑀=0.22. This latter data and fit will be used in the model in the following subsections. Obviously, the exact numbers obtained for 𝑀 and 𝑆 depend on the histogram binning and on the quantity of the data.

Size distributions can potentially be different than lognormal whenever the nucleation and growth processes are not random but are controlled in some way, such as via irradiation-induced nucleation [163], or through nonuniformities in the concentration of silicon, or when size selection is built into in the process as, for example, in the formation of silicon nanocrystals by laser pyrolysis [164], or when size selection is done after crystal formation in solution, or when there is a mechanism for producing nonuniform (i.e., bimodal) distributions. Nevertheless, in numerous studies of silicon nanocrystals a Gaussian distribution is assumed. It can be mentioned that in most cases where a Gaussian fit is used to model the silicon NC size distribution, there is in fact no theoretical basis for it. In fact, it can predict unphysical results (e.g., the Gaussian function does not go to zero when its argument—the radius or diameter—is negative).

The lognormal function should, therefore, be used in cases except those in which there are reasons to expect a different distribution. Since there is a strong dependence of the bandgap energy on the cluster size, accurate modeling of the behavior of ensembles of Si NCs does imply that the size distributions should be well known and correctly modeled. For modeling luminescent properties of ensembles of silicon nanocrystals, as discussed below, the lognormal distribution will be used in every case.

4.2. Pl Spectra and Dynamics

Silicon NCs are characterized by a broad emission spectrum typically peaked between 800 and 900 nm, with a full-width-at-half-maximum (FWHM) of as much as 200 nm, even in “monodispersed” multilayered samples [165]. Although the size distributions for the monodispersed samples in [165] were not actually reported, it seems likely that they would nevertheless be much wider than typically obtained for CdSe NCs due to the lack of post-synthesis size selection methods. Therefore, the broad emission band could simply be due to the range of confined carrier energies in clusters of different sizes, or to the random nature of the interface states that may affect the PL spectrum. Single quantum dots may show homogeneous broadening up to ~150 meV at room temperature [19], although this is still narrower than the overall inhomogeneous broadening due to the size distribution. Below, we will investigate the spectral effects of broadening and nonradiative defects on the emission spectrum of an ensemble of Si NCs.

In order to investigate the theoretical luminescence spectrum and dynamics, the NC size distribution for SiO films annealed at 1000°C from [47] (lognormal radius parameters 𝑆=0.21, 𝑀=0.74) was used along with the theory developed in Sections 13. First, the effect of inhomogeneous broadening was investigated by plotting one spectrum in which the energies of 106 NCs produced with a lognormal probability function were binned in histogram form. In this case, for simplicity the NCs were treated as isolated, and energy transfer between them was not permitted (purple curve, Figure 17). The effect of homogeneous broadening was determined by assigning each cluster a Gaussian range of emission energies with standard deviations ranging from 0 to 200 meV and adding all the resulting curves (Figure 17). As the homogeneous broadening increases, the emission spectrum became broader, effectively similar to a “Gaussian smoothening” of the spectrum. Additionally, the spectrum was observably wider even for small levels of homogeneous broadening, and the peak shifts to slightly shorter wavelengths due to the overall shape of the lognormal function. The spectrum maximizes at shorter wavelengths than actually observed for a real specimen with this size distribution because NC-NC interactions were not yet turned on in the simulation. At this stage, the intent is simply to demonstrate the effect of homogeneous broadening on the observed PL spectrum. One may furthermore propose that the term “monodisperse” may only be used rigorously in situations in which the width of the PL spectra is governed primarily by homogeneous rather than inhomogeneous spectral broadening.

Next, the effect of nonradiative interface defects on the emission spectra and dynamics was estimated. Two size distributions corresponding to those in [47] (𝑆=0.21; 𝑀=0.74) and [162] (𝑆=0.29; 𝑀=0.22) were used, and for the present evaluation NC-NC energy transfer was not permitted. Several different concentrations of nonradiative defects ranging from 0 to 1020 cm−3 were assigned to the clusters with probabilities weighted proportionally to the surface areas. Clusters were checked in random order and a defect was assigned to a cluster if the random number was smaller than 𝐴NC/𝐴allNCs. The process was continued until all defects had been assigned, allowing for multiple defects on a single particle (Figure 18). Every cluster was then given an electron-hole pair with a lowest excited-state energy calculated according to (22). The radiative and nonradiative rates were calculated from (42) and (23), respectively. Using a simple Monte Carlo procedure, the probabilities of radiative and nonradiative decay were calculated for each NC and a random number generator used to determine whether either type of decay occurred. This process continued until no excited carriers remained.

The resulting PL spectra were characterized by an asymmetric curve shape, tailing toward higher energies due to the form of the equations leading up to (22) (Figure 19). The effect of various concentrations of nonradiative surface defects is clearly observable in the simulated PL spectra. In addition to decreasing the overall PL intensity, the effect of the nonradiative centers is to blueshift the peak wavelength. The spectral blueshift occurs because the larger clusters are statistically more likely to contain defects, and the effect would be even stronger if defects were assigned with a volume rather than surface-area weighting (e.g., for volume defects). If extended annealing periods can remove nonradiative centers without increasing the NC size [166], then the effect would be to cause an increased intensity and a redshift of the PL peak, exactly as reported experimentally [53]. This redshift is not due to increasing cluster size, but is due to the removal of nonradiative traps. Finally, the high-energy skewness of the simulated spectra is due to a combination of effects, including the lognormal size distribution, and the lack of carrier transfer from small to large clusters.

As discussed in Section 1, PL decay dynamics in ensembles of Si NCs are characterized by the “stretched exponential” function 𝐼𝑡/𝐼0=exp[(𝑡/𝜏)𝛽], with values of the exponent β often between 0.7-0.8. The stretched exponential has been suggested to be due to a hopping mechanism in which excitons are temporarily trapped and delayed at sites in the oxide matrix or at the cluster-oxide interfaces [32, 167, 168]. However, the stretched exponential decay could also be due to the distribution of lifetimes in NCs of different sizes, as has been indicated by a few groups [35, 169]. In the absence of radiative centers, the recombination energy is given by (22), and the size-dependent radiative rates were developed in (42)). Here we found that when the defect concentration was set to zero (uppermost line in Figure 20), values of 𝛽=0.58 and 𝜏=33μs emerged naturally. The small lifetime is due to the small size of the clusters used in this simulation. Also, since interactions were not enabled, there is no energy transfer from small to large clusters which would also affect the time constant (and the value of β). Including defects changed these values considerably: with increasing defect concentration, the decays are more precipitous (Figure 20). The stretched exponential function with low β values can, therefore, be obtained naturally as a result of the lifetime distribution. Approximately 100-microsecond-timescale carrier hopping, as suggested in previous work [168] as an explanation for the stretched exponential decay, while not ruled out by this analysis, is not required in order to generate the stretched exponential dynamics, a conclusion that is in agreement with recent experimental results [169].

This is not meant to imply that energy transfer does not occur—in fact, it certainly does occur in most samples as discussed below—but that the shape of the decay curve is in fact governed by the size distribution even in the presence of energy transfer. However, energy transfer will also affect the shape of the decay curve—in the absence of nonradiative processes shifting the decay to longer lifetimes and larger β. In fact, in the absence of nonradiative decay, it is in theory possible to work back from the observed decay to the lifetime distribution using the inverse Laplace transform for the stretched exponential function [170]. From the lifetime distribution one could in principle then work back to the size distribution if 𝜏𝑝𝑙𝜏rad and the radiative rate model is accurate. If nonradiative processes are included, the values of β and 𝜏 are sensitive to the size distribution, filling fraction, and defect concentration. In the discussion below, the first interaction simulation including Forster transfer and tunneling will be performed and its effect on the PL spectrum investigated.

4.3. Interactions

In order to achieve better size control, several studies have formed silicon nanocrystals using a thin film multilayering approach. Instead of depositing a single thick layer of SiO𝑥, alternating thin layers of SiO𝑥 and SiO2 were used. In this way, the thin SiO𝑥 layers are isolated from one another by the SiO2 buffer layers. Upon annealing, silicon clusters crystallize in the SiO𝑥 layers, but the cluster growth occurs only by diffusion within the layer as opposed to three-dimensionally. Although the resulting cluster sizes should clearly be smaller than for three-dimensional diffusion in a thick layer of the same composition, it is not immediately clear without more detailed modeling whether the distribution should, in theory, be more narrowly distributed with respect to the mean value. However, narrower distributions have been reported in many cases [171]. Nevertheless, even in apparently monodispersed multilayered samples [165] the PL peak is nearly 200 nm wide at the half maximum.

Through investigating the effect of multilayering of the Si NCs on the luminescence spectrum, it has been observed that the PL spectrum and dynamics from silicon nanocrystals depends on the thickness of both the SiO2 buffer layers and the active nanocrystal layers, even if the size distributions remain approximately constant [23, 32, 162]. This is evidence of an interaction mechanism operating between the clusters: as the NCs become more isolated—in one direction at least—by the buffer layers, the PL peak shifts to shorter wavelengths and there is a change in the decay dynamics as well, with a possible trend toward higher values of the exponent β. As already discussed, this has been widely interpreted in terms of a kind of excitonic hopping or migration (e.g., see [23]) between clusters, although only more recently have the mechanisms for carrier transfer among Si NCs been more clearly elucidated [30, 160].

There are at least two fundamental means by which closely spaced nanocrystals can exchange charge carriers. The first mechanism is via tunneling of individual electrons or holes from one cluster to another. If the two clusters are approximately the same size (or if the higher levels of a larger cluster overlap with lower levels of a smaller one), then the tunneling is resonant, it is equally probable in either direction, and its rate depends on the negative exponent of the separation distance. Tunneling is well known to occur in silicon clusters and forms the basis of silicon nanocrystal memories [172, 173]. Accurate calculation of the tunneling rate in ensembles of nanocrystals with different energies and spacings is a difficult problem, however. The tunneling rate can be approximated by the Wentzel-Kramers-Brouillin (WKB) approximation, as has been established for the case of double quantum wells [174176]:𝑤𝑡=1𝑡tunnel=12𝑑𝑝2𝐸𝑚𝑝1/216𝐸(𝑉𝐸)(𝑚𝑝/𝑚𝑏)𝑚𝑉+𝑝/𝑚𝑏𝐸12×exp[2𝑑𝑏2𝑚𝑏(𝑉𝐸)21/2],(46) where 𝑚𝑏 and 𝑚𝑝 are the carrier effective masses in the barrier and the particle, 𝑑𝑝 and 𝑑𝑏 are the particle diameter and barrier thickness, 𝑉 is the barrier height, and 𝐸 is the carrier energy. This equation (using the corrected form shown in [160]) has previously been employed to estimate tunneling rates between asymmetric quantum wells [176], despite the fact that (46) is strictly applicable only to the resonant case.

If the energy levels of adjacent clusters are mismatched, then tunneling must be “assisted” by emission or absorption of phonons. The tunneling rate dependence on the energy gap, Δ𝐸, and the phonon energies (within the clusters and in the SiO2 matrix) for these mechanisms is not well known and may need to be estimated experimentally if sufficiently narrow size distributions can be obtained. Recently, computer simulations showed [160] that a large Stokes shift between absorption and emission in the specimen as a whole can occur as a result of tunneling between nanocrystals although, as for the quantum well case, (46) was used without accounting for phonon-assisted tunneling.

The second fundamental interaction mechanism is resonance energy transfer, or “Forster transfer” [177]. In this mechanism, electron-hole pairs may “migrate” from particle to particle by a dipole-dipole or higher-order multipole coupling, that is, the Forster process—which has already been established for the case of CdSe nanocrystals [178]. Others have shown that it applies to silicon nanocrystals also [30], although due to the indirect gap the rates can be much different (typically smaller) than for CdSe. Forster transfer is the basis for fluorescence resonance energy transfer (FRET) microscopy, and has been widely investigated in the case of interactions amongst the rare earth ions [179]. The distance dependence of the transition dipole transfer rate, 𝑤𝑑, is given by 𝑤𝑑=𝑤PL(𝑅0/𝑟)6, where 𝑤PL is the PL rate, 𝑅0 is the distance at which the transfer rate and the PL rate are equal, and 𝑟 is the distance between the edges of the particles. 𝑅0 itself depends on the spectral overlap integral between the donor (D) and acceptor (A) clusters, and is, therefore, technically different for any pair of nanocrystals undergoing this process. Previous work in CdSe nanocrystals has assumed that the energy levels must be resonant for the Forster transfer to occur [178], however this is not strictly true: phonon-assisted Forster transfer has been calculated for the rare earth ions and should apply to nanocrystals as well. In the Miyakawa-Dexter theory [180], the resonant Forster rate is modified by considering the electron-phonon coupling parameters in the density-of-states function for the donor and activator ions. In this way, one obtains [180]:𝑤dip=𝑤𝑝𝑙𝑅0𝑟6exp𝛽Δ𝐸,(47) where β depends on the electron-phonon coupling strength, and its value can be extracted from a plot of the transfer rate versus Δ𝐸 in silicate glass. From the data provided in [181], a rough estimate gives 𝛼=0.0018 for rare earth ions embedded in SiO2, where α indicates internal transitions and is about twice the value of β [180]. Unfortunately, the value of β is much more difficult to obtain experimentally for silicon NCs, as a result of the size distributions (i.e., every NC is slightly different; whereas all rare earth ions of the same type have identical properties).

The Miyakawa-Dexter model has been widely discussed and applied extensively in the case of transfer between different rare earth ions [182184]. In the present case, one may assume that the energy difference is lost to phonons in the surrounding silicate glass, as would be the case for isolated “point” dipoles such as the rare earth ions. As the energy gap grows, the transfer rate decreases exponentially, and experimentally, for rare earth ions at least, this trend shows little variation despite the many simplifications in the Miyakawa-Dexter model [185]. The effective Forster distance, 𝑅0, is more difficult to obtain, since it is not practical to calculate the overlap integral (including phonon sidebands) for every possible pair of interacting particles in a sample. Therefore, an experimental approach is taken to approximate a single 𝑅0 value for all pairs of particles. By varying the distance between layers of interacting particles and using the geometrical constraints imposed by interacting layers (as opposed to two point dipoles) [186], an estimated value of 𝑅0=5 nm was used in the simulations. Obviously, this is a significant simplification and it must be stressed that the results remain approximate until phonon-assisted multipole transfer between all possible pairs of NCs can be calculated in a reasonable way.

The simulation setup was done exactly as described previously: every particle was populated with a single electron-hole pair at time 𝑡=0. This time, the effect of nonradiative traps was not included in order to show only the energy migration effects. Next, for every nanocluster the probabilities associated with each of the four processes described above (radiative decay, nonradiative decay, tunneling, and Forster transfer) were calculated in random order (the NCs themselves are also checked randomly) and compared to a flat-probability random number between 0 and 1 to determine whether one of these four events occurred. When the calculation started, the simulation timestep was 1×1012 s. After iterating through all the particles, the time increments by one picosecond and the process repeats. If ten complete iterations through all the clusters (in randomized order each time) produce no events of any kind, the timestep is multiplied by a factor of ten. In this way, it is possible to overcome the otherwise impossible computational limits associated with the vastly different rates for the different processes. The simulation ends when there are no electron-hole pairs left. When a particle radiates, the energy of the outgoing photon is stored for the spectral output function, and the time is recorded in order to plot the decay of the luminescence. Finally, a set of tests were performed to find the minimum sample size for which edge effects (due to the finite box size) became unobservable. A simulation of 5,000 particles was found to provide sufficiently short computation time (~12 hours on a 2.6 GHz CPU) and no observable edge effects. Auger effects were not simulated, which effectively approximates the behavior of a lightly pumped ensemble.

Figure 21 shows the PL spectrum with and without the Forster and tunneling “energy migration” interactions. In the first few nanoseconds, numerous tunneling events and dipole-dipole interactions occurred. The rates for these processes for adjacent nanocrystals are many orders of magnitude faster than the radiative decay processes. Some nonradiative trapping can also occur early on in the simulation. Even after only a few picoseconds, there was already a small fraction of charged nanocrystals due to electron or hole tunneling events. Many nanocrystals had no remaining charge carriers, and many others—predominantly the larger ones—had multiple electron-hole pairs. After 500 ns, many thousands of tunneling and dipole events occurred and most of the smaller clusters were devoid of charge carriers. In the microsecond timescale, radiative decay joins the longer-range energy transfer events, and by 100 microseconds only a few excited nanocrystals remained. The long wavelength of the emitting states (compared to absorption)—as well as a decrease in the overall skewness of the PL peak—in dense ensembles of Si NCs can originate predominantly from energy migration, indicating that the large Stokes shift in ensembles of Si NCs can be attributed, at least in part, to these effects. We have also found that the degree of “stretching” of the PL decay depends on the resulting ensemble of particles with electron-hole pairs; in the absence of nonradiative effects the PL lifetime is longer but if defects are introduced the PL lifetime can be dramatically shorter in the presence of energy migration, due to the greater statistical probability for defects to occur on larger NCs. Based on the initial results presented here, it would seem that the combination of theory with computational simulations may be one of the best tools to obtain a more thorough understanding of the complicated and controversial luminescence behavior of ensembles of silicon clusters.