Co-doped silicon nanoparticles (NPs) are promising for the realization of novel biological and optoelectronic applications. Despite the scientific and technological interest, the structure of heavily co-doped Si NPs is still not very well understood. By means of first principles simulations, various spectroscopic quantities can be computed and compared to the corresponding experimental data. In this paper, we demonstrate that the calculated infrared spectra, photoluminescence spectra, and Raman spectra can provide valuable insights into the atomistic structure of co-doped Si NPs.

Few nanometers-sized Si nanoparticles (NPs) have been in the focus of nanoscience as their size-dependent bandgap enables the tuning of their electrical and optical properties.1 Due to the abundance and environmental friendliness of elemental Si, and the accumulation of knowledge from the semiconductor industry, Si NP-based solutions are a very attractive and promising platform over similar but less developed technologies based on compound semiconductors. An advantageous feature of Si NPs is that they can be doped by various dopants. In particular, Si NPs can be co-doped by boron and phosphorus beyond the bulk solubility of dopant atoms.2–10 These heavily co-doped Si NPs can be fabricated by the annealing borophosposilicate glass (BPSG) matrix.7 The size of the nanoparticles can be controlled by the annealing temperature in the range of 1–14 nm.7 The resulting nanoparticles can be extracted from the BPSG matrix by HF etching. According to scanning transmission electron microscopy (STEM)2,11 and X-ray diffraction measurements,12 the extracted co-doped Si NPs have highly crystalline core encapsulated by an amorphous layer.11 The concentration of B and P in the nanocrystals determines their dispersibility in solution.13 The resulting co-doped Si NPs exhibit numerous advantageous properties compared to undoped Si NPs: they are stable, resistant against oxidation, and dispersible in aqueous environments.13,14 Similarly to undoped Si NPs, they exhibit size-dependent photoluminescence (PL) in the visible and near-infrared (NIR), but the absorption edge and PL peak are generally red-shifted compared to their non-doped counterparts at identical NP diameters.10 For larger diameters, the luminescence peak falls beyond the bandgap of bulk Si.2 Overall, heavily co-doped Si NPs are very promising for the realization of novel optoelectronics applications.

Electron energy loss spectroscopy (EELS) can be utilized13 to confirm the presence of dopants in Si NPs. The dopant concentration of heavily doped Si NPs can be estimated by inductively coupled plasma atomic emission spectrometry (ICP-AES) measurements,13,15 or atom probe tomography (APT),16 which also offers some insight into the distribution of dopants within the Si NPs. According to ICP-AES measurements,10,13 the concentration of dopants can be controlled by the boron and phosphorous concentration of the BPSG matrix. The B and P concentration of extracted Si NPs generally exceeds the dopant concentration of the BPSG matrix,13 and the co-doped Si NPs typically contain significantly more B than P.10 It was also established that the annealing temperature does not have a significant impact on the dopant concentration.10 APT analysis of co-doped Si NPs can offer additional insight about the dopant distribution within the Si NPs. According to a recent study16 of Si NPs embedded in BPSG, P atoms are more likely to be found in the inner region of the Si NPs whereas B accumulates near the NP-BPSG boundary. By contrast with ICP-AES results, the P concentration was found to exceed the concentration of B. While this contradiction might be resolved by undetected B counts,17 it is not clear how these findings translate to the case of freestanding Si NPs. In addition to the radial distribution of dopants, the same study also established the abundance of B–P and P–P clusters within the Si NPs. X-ray photoelectron spectroscopy (XPS) measurements provide further information about the chemical configurations of the dopant atoms.8,18 In particular, the presence of non-oxidized B and P atoms in the proximity of the surface was established.8 

Infrared (IR) spectroscopy and Raman spectroscopy measurements can shed further light on the atomistic structure of Si NPs. Raman spectroscopy offers insights about the lower energy vibrations (300-800 cm−1),8 where IR spectroscopy is ineffective. IR spectroscopy measurements can complement these results by the identification of oxygen and hydrogen related vibrational bands in the 700-3500 cm−1 region.11,13,19 Limited amount of structural information can be also deducted by analyzing the photoluminescence (PL) spectra. Traditionally spectroscopy measurements are performed on large ensembles of nanoparticles; thus, the obtained PL spectra are featureless2–4,6,10,13 due to the slightly varying zero-phonon line (ZPL) of the different NPs. However, it is also possible to measure the PL signal of individual Si NPs,20 where the shape of the PL spectra shows some characteristic features at low temperatures. Comparing the PL spectrum of individual nanoparticles at different temperatures or single NP measurements with ensemble measurements, one can partition the broadening of PL spectra into contributions from zero-point vibrations, thermal broadening, and inhomogeneous broadening.

Despite the significant experimental and computational effort, the structure of heavily co-doped Si NPs is still not understood at the atomic level to a satisfactory degree.

Raman, infrared, and photoluminescence spectra can also be calculated utilizing first-principles computational methods. Calculating the vibrational density of states (VDOS) can be considered as routine and offers valuable information about the vibrational energies of molecules or bulk material. The IR spectrum of a system can be obtained with little additional cost by calculating the effective Born charges, and the result offers insight not only about the characteristic vibrational energies but the relative intensities of the IR bands. Our goal was to predict the experimental signals of certain structural features of co-doped Si NPs in order to improve the current understanding of their atomic structure. After the vibrational signature of clusters of few atoms is identified, the results can be compared to the experiments and the gained insight can be utilized to build accurate atomic models for co-doped Si NPs.

The IR spectra are calculated from the effective Born charge tensors and the normal modes as

IIRν=μlZl,μ*ulν2,
(1)

where l is the atomic index, uν is the νth normal mode and the Born effective charges are defined as

Zl,μλ*=2Etotul,μEλul,μ=0,Eλ=0,
(2)

where Etot is the total energy of the system and Eλ is the electric field applied in the λ = x, y, z direction.

We calculate the PL spectra within the Franck-Condon approximation using the generating function method developed by Kubo and Toyozawa21 and Lax.22 The emission spectra are calculated using the so-called partial Huang-Rhys factors defined by

Sk=ωkqk22,
(3)

where Sk is the partial Huang-Rhys factor of the kth mode, ωk is the energy of the kth mode, and qk is given by

qk=imi(RiexRigr)uk;i,
(4)

where Rgr and Rex are the atomic coordinates in the ground and excited state, m is the atomic mass of the corresponding atom, and uk is the kth normal mode. We neglect the effects of the Duschinsky rotation, and thus ωk and uk correspond to equilibrium coordinates of the electronic ground state.

From the partial Huang-Rhys factors, the emission spectra can be obtained (see Ref. 23) by constructing the so-called electron-phonon spectral function

S(ω)=kSkδ(ωωk)
(5)

and two auxiliary functions defined as

f±(t)=S(ω)nωnω+1e±iωtdω,
(6)

where nω is the Bose-Einstein statistics nω(T)=[exp(ω/kbT)1]1. At T = 0 K, Eq. (6) simplifies to

f(t)=S(ω)eiωtdω.
(7)

With these auxiliary functions, the generating function G(t) can be constructed as

G(t)=eiE0t/f+(t)+f(t)f+(0)f(0),
(8)

where E0 is the energy of the no-phonon transition or the ZPL energy. Finally, the emission line shape can be obtained as

I(ω)=nω33ϵ0πc3|μ|2G(t)eiωtγ(t)dt,
(9)

where n is the refractive index of the media (that is unity for simulation of Si NPs in the gas phase), c is the speed of light, μ is the transition dipole moment between the ground and excited electronic states, and γ(t) can be chosen appropriately to describe the inhomogeneous broadening of the spectral lines. The choices of γ(t) = γ1|t| and γ(t) = γ2t2 results in Lorentzian and Gaussian-broadening, respectively.24,25 The description of inhomogeneous broadening is outside of the scope of this work; thus, we present our results with small, artificially chosen Gaussian-broadening (2-5 meV), in order to improve the clarity of our PL plots.

In this approximation, one has to obtain the following quantities to calculate the PL emission spectra: the ground and excited state geometries, the vibrational energies and normal modes, and the ZPL energy. This method has been proven to reliable for the description of the PL spectrum emitted by color centers in solids.26,27

Resonant Raman spectra was calculated using the Savin28,29 formula, where the relative intensities of fundamental Raman intensities of totally symmetric vibrational modes are given by

Ij,10Ik,10=(ωIωj)3(ωjΔj)2(ωIωk)3(ωkΔk)2ωj2Sjωk2Sk,
(10)

where ωI is frequency of the incident light, ωj, ωk are the frequencies of the corresponding vibrational modes, and Δk is the dimensionless displacement between the ground and excited electronic state minima that is related to the partial Huang-Rhys factors defined in Eqs. (3) and (4) by Si=Δi2/2. The approximation in Eq. (10) holds when ωIωj, ωk. Generally, the Savin formula is expected to give a reasonable estimation for the resonance Raman spectra, when (i) short time dynamics govern the Raman process,24,29 (ii) the Duschinsky rotation can be neglected, (iii) the ground state vibrational frequencies are good approximates for the excited state ones, (iv) and the overtones in the Raman spectra can be neglected (true if Si ≪ 1). Similar computational approach to obtain the resonant Raman spectra was already utilized in the context of CdSe nanocrystals30 and silicon nanocrystals covered with an amorphous surface layer.31 

In addition, we comment on the excitation wavelengths in our simulation and its implications when compared to experimental data. To calculate the resonance Raman spectra, the ground and excited state geometries and the vibrational properties are determined from first principles calculations. In density functional theory (DFT) based calculations, the relaxation of the geometry in excited electronic states is very challenging. Time-dependent density functional theory (TDDFT) can be utilized to obtain excited state properties of molecules,32 but the structural relaxation of nanocrystals with ∼1000 electrons is extremely time consuming in the computations. Instead, we utilized the ΔSCF method where excited state electronic structure was mimicked by forcing non-Aufbau electron occupations. While this method is justified for geometry relaxation in the lowest lying excited state, it fails for higher excitation. This is not an issue for the calculation of PL spectra, where the excited system quickly relaxes to the lowest lying excited electronic state (Kasha’s rule) via vibrational relaxation; thus, the lowest lying electronic excited state participates in the emission. By contrast, the higher excited electronic states should be considered to calculate the resonance Raman spectra, unless the laser frequency coincides with the ZPL. This is not the case in a realistic experimental scenario, where the laser frequency is usually significantly higher than the ZPL. This may result in a larger error in the calculated Raman spectra than that in the calculated PL spectra.

In our work, we employed various projection techniques, to gain more insight about correlations between structural and vibrational properties of co-doped Si NPs. In this approach, we examine which normal modes are localized on a given set of atoms (e.g., a given atomic species, or atoms at the NP surface) by calculating

Pi,S=α,kS|ui,k,α|2,
(11)

where α = {x, y, z}, k is the atomic index, S is the chosen set of atoms, and Pi,S is the weight of the ith normal mode projected to the atoms of set S. After the calculation of the projected weights for each normal mode, the VDOS, IR spectra, and Raman spectra can be decomposed as these spectral quantities are calculated as a sum over the normal modes.

We utilized density functional theory (DFT) and time-dependent density functional theory (TDDFT) to study selected nanostructure models. We calculated the relaxation of structures and the vibrational energies by the use of the semi-local Perdew-Burke-Ernzerhof (PBE)33 exchange-correlation functional. During the structural optimization, the atoms were relaxed until all of the forces fell below 0.02 eV/Å. Before calculating the vibrational energies and normal modes, the structure was further relaxed to <0.001 eV/Å in order to achieve accurate vibrational energies. The dynamic matrix was constructed utilizing the finite-differences method with 0.015 Å displacements. Born effective charges were calculated utilizing density functional perturbation theory (DFPT).34 The vertical excitation energies and transition dipole moments were obtained from linear-response TDDFT calculations with the PBE035 hybrid-functional in the adiabatic TDDFT kernel. Excited state geometries were calculated by forcing non-Aufbau occupation numbers according to the first optical excitation. We utilized the plane-wave Vienna Ab initio Simulation Package (VASP)36,37 code for the ground state and excited state structure relaxation, the determination of vibrational energies and normal modes, and the computation of effective Born charges. VASP utilizes a plane-wave basis set where the core electrons were taken into account within the projector augmented wave (PAW)38 framework. Because of the periodic boundary conditions, at least a 10 Å vacuum was ensured between the nanoparticles and their mirror images, in order to minimize the spurious nanoparticle-nanoparticle interaction. The TDDFT calculations were carried out utilizing the Turbomole code,39,40 with a gaussian-type double-zeta basis set with additional polarization (DZP) orbitals and the valence electrons were taken into account utilizing electric core potentials (ECPs).

We found that our first principles methodology did not always give accurate results for the vibrational energies. Most notably, our calculations underestimate the energy of the transverse optical (TO) phonon peak of silicon which is a prevalent signal in Raman studies of Si NPs. On the other hand, we found that the calculated frequency of the substitutional B atom (618 cm−1) is in very good agreement with the experimental value41 for B11, 620 cm−1. As there are two common isotopes of B, one would have to consider both 10B and 11B in the simulation process. For clusters containing multiple B atoms, the number of possible combinations scales exponentially with the number of B atoms drastically increasing the complexity of the calculations. To simplify our task, we performed the calculations using the abundance-weighted average standard atomic weight for B (10.811), so the local mode energies of 10B and 11B atoms are expected to slightly differ from our results. The ratio of the vibrational frequencies of 10B and 11B in bulk Si is ≈1.03941 and the 1.03-1.05 ratios were found to hold for isotopically pure complexes of 10B/11B too.42 As a result, we expect that our results slightly overestimate the vibrational frequencies of 11B (<1% error), while underestimate the vibrational frequencies of 10B (3%-4% error).

First, we discuss our results regarding the the characteristic local mode energies of certain boron and phosphorous related clusters that are expected to occur in co-doped Si nanoparticles. During this study, our goal was not to create realistic atomistic models of co-doped nanocrystals, but to study the vibrational properties of small, few atom-sized atomic structures embedded in the Si NP, e. g., substitutional dopants, and small dopant clusters. This approach can be justified by the observation that the global structure of the Si NPs have only small effect on the vibrational energies belonging to the covalent bonds of the investigated clusters.

While it is impossible to investigate every possible structure that might possess a unique vibrational footprint in the experimental characterization of co-doped Si NPs, we believe that our extensive study is able to assist and guide the experiments in the interpretation of the observed spectra.

First, we investigate the (projected) vibrational density of states of a model Si NP that is terminated by H atoms but does not contain B or P dopants. Figure 1 shows the projected vibrational density of states (PVDOS) of a 1.6 nm Si NP model. Utilizing projection to the Si, (Si)–H, (Si)–H2, we can analyze the nature of IR bands. Apart from the characteristic stretching modes around 2100-2200 cm−1, the rolling, rocking, and twisting modes of the (Si)–H, (Si)–H2 moieties can interfere with the local modes of dopant-containing structures. The measurements indicate13,19 that H is either absent or present in diminishing quantities in co-doped Si NPs; however, we terminate the surface of our Si NP models with hydrogen in order to reduce the number of electrons entering our simulations. Fortunately, the calculation of PVDOS allows us to draw conclusions about the contributions of the local modes of specific atomic species. This is especially important for the Si–H/Si–H2 related IR band around 620 cm−1 which happens to coincide with the important Raman signal of substitutional B. Even after subtracting the H contribution from the VDOS, a small peak remains signaling that the Si atoms of the outermost layer vibrate differently compared to the Si atoms of the bulklike region. In the experiments, this IR peak is absent for undoped Si NPs, which can be explained by the smaller surface/volume ratio of experimental NP ensembles (usually 3-7 nm), the oxidation of the NP surface, or the Raman inactivity of such modes.

FIG. 1.

Vibrational density of states of a hydrogen terminated, 1.6 nm Si NP built up from 147 Si atoms. The VDOS was decomposed to show the contribution of surface (Si–H and Si–H2) and bulk (Si–Si) vibrations.

FIG. 1.

Vibrational density of states of a hydrogen terminated, 1.6 nm Si NP built up from 147 Si atoms. The VDOS was decomposed to show the contribution of surface (Si–H and Si–H2) and bulk (Si–Si) vibrations.

Close modal

We also calculated the PL spectra of H-terminated undoped Si NPs with diameteres in the range 1.1–2.8 nm and plotted the results in Fig. 2. The quantum confinement effect (QCE)43 is significant in this NP size range, as the PL peak shifts from 3.2 eV to 2.3 eV with increasing NP size. In addition, the vibration induced broadening is also size-dependent. For the smaller NPs, the calculated emission spectra show no distinguishable ZPL due to the strong vibration induced broadening. As the size of the NPs increases, the Debye-Waller factor increases, and the PL is dominated by the ZPL at low temperatures, in good agreement with the experimental data.20 We also calculated the Raman spectra, where the size-dependence of the transverse optical (TO) phonon peak of the Si crystal can be identified in accordance with earlier experimental44–47 and theoretical48–50 results that also confirmed that the shift of the TO-peak towards lower energies is caused by QCE. It is also notable that the shape of the calculated Raman spectra also changes with the diameter. For the smallest Si NP models, there is a significant feature in the 600-700 cm−1 region that gradually disappears with increasing NP size. This indicates that the 600-700 cm−1 feature is related to the surface of the NPs, which we confirmed by calculating a projected Raman spectra where the projection excluded the outermost Si layer and the H atoms.

FIG. 2.

The calculated photoluminescence emission spectra and resonant Raman spectra of hydrogen-terminated, pristine Si nanocrystals with diameters in the region of 1.1-2.8 nm. For the PL spectra, black and red curves correspond to 0 K and 300 K, respectively. For the Raman spectra, black curves correspond to the full spectra, while green curves correspond to the projected Raman spectra where the outermost layer of Si atoms and H atoms is excluded from the projection.

FIG. 2.

The calculated photoluminescence emission spectra and resonant Raman spectra of hydrogen-terminated, pristine Si nanocrystals with diameters in the region of 1.1-2.8 nm. For the PL spectra, black and red curves correspond to 0 K and 300 K, respectively. For the Raman spectra, black curves correspond to the full spectra, while green curves correspond to the projected Raman spectra where the outermost layer of Si atoms and H atoms is excluded from the projection.

Close modal

The most essential vibrational signal related to the point defects is the characteristic local modes of substitutional phosphorous and boron atoms. From an experimental viewpoint, their significance is different as the phosphorous related signal appears as a shoulder on the lower energy side of the silicon TO phonon peak;7 thus, it is rather difficult to determine its intensity. While the TO related peak could be easily subtracted from the total spectra for bulk Si, the picture is more complicated for nanocrystals as the shape of the TO peak becomes size and structural dependent. On the other hand, the local modes of substitutional boron produces a well distinguishable peak in the Raman spectra7,8 for which the intensity and shape seem to vary with doping concentration.11 Because of this fact, we focused on energy and intensity of the substitutional boron as a promising indicator of atomic-scale properties.

First, we investigated the case of a single, isolated boron atom inside a hydrogen terminated Si NP. Our goal was to determine whether the position of B atoms within the NP affects the energy of its local modes. As experiments indicate that co-doped Si NPs are compensated, we gave an overall negative charge to the B containing Si NP. Figure 3 shows the boron-projected vibrational density of states for a single B containing, 2.4 nm sized Si NP. It is clear that the position of the B atom has a significant impact on the vibrational energy of the B–Si related normal mode.

FIG. 3.

Boron-projected vibrational density of states of a 2.4 nm sized Si NPs containing a single substitutional B atom with different boron positions indicated in the corresponding pictograms. The B atoms are negatively charged and four-fold coordinated except from the uppermost model, where the B atom sits at the surface of the Si NP, is three fold coordinated, and has a neutral charge state. The numbers above the spectra correspond to the B atoms’ distance from the outermost layer of Si atoms.

FIG. 3.

Boron-projected vibrational density of states of a 2.4 nm sized Si NPs containing a single substitutional B atom with different boron positions indicated in the corresponding pictograms. The B atoms are negatively charged and four-fold coordinated except from the uppermost model, where the B atom sits at the surface of the Si NP, is three fold coordinated, and has a neutral charge state. The numbers above the spectra correspond to the B atoms’ distance from the outermost layer of Si atoms.

Close modal

For B atoms far from the surface, the B related local modes are degenerate and fall around 615 cm−1. By contrast, the B atom in the proximity of the surface exhibits significantly higher vibrational energies with broken degeneracy in the region of 630-660 cm−1. The most extreme examples are the B atoms introduced to the outermost and the second outermost atomic layer of the Si NP. For the second outermost layer, the B related peak is split to a 765 cm−1 mode and several lower energy modes in the 600-700 cm−1 region. For the outermost Si layer, the B atom is only three-fold coordinated and is in a neutral charge state. This atomic configuration results in a vibrational energy of 735 cm−1 which is much higher than local modes of the four-fold coordinated B residing in the bulk region of Si NPs. We find that the transition point between bulk-like and surface-like local modes of substitutional B falls between 3.8 and 5.4 Å distance measured from the outermost Si layer.

The case of B proximate to the NP surface is particularly important, as experiments indicate a strongly B enriched, amorphous layer encapsulating the crystalline NP core.11,51 In addition, the concentration of B in B/P co-doped Si NPs is usually higher, while the electronic structure of the co-doped Si NPs indicates that they are compensated; thus, the excess B must be inactive. This is seemingly consistent over a relatively broad range of B and P concentration, raising the questions about the nature of inactive B.

Neutral, triply coordinated B atoms are the simplest B defects that do not introduce acceptor states. In contrast to substitutional B in bulk Si (where neutral B always introduces an acceptor state), it is possible to introduce B atoms to the outermost Si layer of Si NPs that are inactive. Because of the significance of triply coordinated B atoms of the outermost atomic layer, we also investigated clusters of triply coordinated B atoms on the NP surface (see Fig. 4). These surface clusters are formed by second neighbor three fold coordinated boron atoms. Introducing a cluster of 3 three fold coordinated B atoms to the surface, the B related vibrational modes fall between 710 and 760 cm−1. This IR band becomes even wider as one of the facets of the Si NPs is fully terminated by three fold coordinated B atoms, while the shape of the TO peak remains preserved. Finally, the model where surface Si atoms are randomly substituted by three-fold coordinated B atoms at 20% total B concentration was considered. This model exhibits an even broader, featureless band falling in the 640-860 cm−1 region.

FIG. 4.

Silicon- and boron-projected vibrational density of states of Si NP models containing three fold coordinated substitutional B atoms at their surface. The hexagonal shape of the Si NP silhouettes emphasizes the faceted nature of the model Si NPs. From the top to the bottom: fully B terminated Si NP, Si NP with one of its facets terminated by B atoms, Si NP with 3 B atom on one of its facets, and finally a Si NP with a single B atom on its surface. The B-projected VDOS is scaled up to 25% concentration equivalent for sake of clarity.

FIG. 4.

Silicon- and boron-projected vibrational density of states of Si NP models containing three fold coordinated substitutional B atoms at their surface. The hexagonal shape of the Si NP silhouettes emphasizes the faceted nature of the model Si NPs. From the top to the bottom: fully B terminated Si NP, Si NP with one of its facets terminated by B atoms, Si NP with 3 B atom on one of its facets, and finally a Si NP with a single B atom on its surface. The B-projected VDOS is scaled up to 25% concentration equivalent for sake of clarity.

Close modal

Next we turn our attention towards B–P complexes within Si NPs. As the local modes of P are expected to be masked by the TO peak of the Si NP, our main focus is to determine how the proximity of P atoms affects the local modes of B. Our previous work indicated that the formation of neighboring B–P pairs is energetically preferred; thus, the vibrational energies of small B–P cluster are of interest. The concentration of B always exceeds the concentration of P in experimental samples of co-doped Si NPs, and thus we only considered B–P cluster (which are favored energetically) and B–B clusters (which are energetically not favored, but still feasible to occur due to the high concentration of B dopants).

Figure 5 shows the B-, P-, and Si-projected vibrational density of states of various Si NP models containing B and P atoms. Figures 5(a) and 5(b) show the case of individual B and P, respectively, whereas (c) and (d) are the case of separated and first neighbor B–P pair, respectively. For the single B and P dopants, a total charge of −1 and +1 was applied, respectively, in order to compensate the extra hole/electron. Figure 5(g) shows the projected vibrational density of states of a Si NP model containing a ring-like B–P cluster, while (h) shows the PVDOS of a Si NP model for which the innermost 26 atoms were substituted by a boron phospide core. Finally, Fig. 5(i) is the PVDOS of a model containing 3 fourfold coordinated P and B atoms and 3 threefold coordinated B atoms at the NP surface.

FIG. 5.

Silicon-, boron-, and phosphorus-projected vibrational density of states with a pictogram representing the dopant configurations of the corresponding co-doped Si NPs models. Green and orange dots with black rectangles represent four fold-coordinated substitutional B and P atoms, while green dots with black triangles represent threefold coordinated B atoms at the surface of the Si NPs. Adjacent dots represent covalent bonds between the relevant dopant atoms whose remaining bonds are connected to adjacent silicon atoms. For the sake of clarity, we scaled up the B- and P-projected VDOS (to an equivalent intensity of 25 B/P per Si NP) meaning that the PVDOS intensities between different subplots are not comparable.

FIG. 5.

Silicon-, boron-, and phosphorus-projected vibrational density of states with a pictogram representing the dopant configurations of the corresponding co-doped Si NPs models. Green and orange dots with black rectangles represent four fold-coordinated substitutional B and P atoms, while green dots with black triangles represent threefold coordinated B atoms at the surface of the Si NPs. Adjacent dots represent covalent bonds between the relevant dopant atoms whose remaining bonds are connected to adjacent silicon atoms. For the sake of clarity, we scaled up the B- and P-projected VDOS (to an equivalent intensity of 25 B/P per Si NP) meaning that the PVDOS intensities between different subplots are not comparable.

Close modal

The B- and P-projected VDOS of the BP pair approximates the B- and P-projected VDOS of the single-dopant model, with minor differences. For the nearest neighbor pair, a B-P mode appears around 560 cm−1, well distinguishable from the 620 cm−1 B–Si modes. For the distant B–P pair, this softer mode does not occur, but the B related modes are shifted towards higher energies due to the proximity of the NP surface. The P-projected VDOS remains scattered in a broad energy range, rendering the P-related bump on the low wavenumber side of the TO-peak less valuable for quantitative analysis. Even the introduction of a simple B–B–P cluster induces a relatively complex B-related signal in the vibrational spectrum [see Fig. 5(e)].

By calculating the B-projected VDOS and investigating the normal modes that contribute the most to the B-projected VDOS, we can identify the following vibrations: the stretching mode of the B–B bond around 850 cm−1, the rocking mode of the B–B bond around 700 cm−1, and the mode that is localized on the B atoms, qualitatively similar to the normal mode of a single fold coordinated B atom. We also find a mode that is localized on the B, P, and neighbor Si atoms around 633 cm−1 and two modes that are localized in the B atoms and the proximate Si atoms at 606 and 607 cm−1. At 588 cm−1, there is a mode that is localized on the B–P bond, some neighbor Si atoms, and also hybridized with the bending modes of the Si–H bonds at the surface. The 526 cm−1 mode is localized on the B, P dopants and the proximate Si atoms. The broad IR band below 440 cm−1 can be attributed to B and P related vibrations hybridized by rather delocalized, almost bulk-like vibrations of the host Si NP.

In the case of the B–P–B cluster [see Fig. 5(f)], the B–B related local mode is missing, and the most distinct IR band (720-750 cm−1) is related to the vibration of the threefold coordinated B atom bonding to a P atom and two Si atoms. The inner B atom contributed to the VDOS in the 580-670 cm−1 region where the vibrations are localized in the first neighbor atoms.

The 6-atomic ring-like BP cluster [Fig. 5(g)] introduces a rather broad IR band between 520 and 690 cm−1. The introduction of a small cubic boron cluster in the center of the Si NP results in a similar broad band from 530 to 730 cm−1. Interestingly, the characteristic TO and longitudinal optical (LO) lines (799 and 829 cm−1, respectively52) of cubic boron phosphide are absent from the PVDOS, indicating that the stress induced by the mismatch between the lattice constant of boron phosphide (4.54 Å) and Si (5.43 Å) changes the vibrational properties of the BP core. The lack of >800 cm−1 frequencies indicate negative pressure on the BP core53 which is consistent with differences between the lattice constants.

Finally, we investigate the case of the surface-bound BP cluster [Fig. 5(i)] that is expected to occur due to the concentration of dopants on the NP surface. In this cluster, the presence of B in both of its three- and fourfold coordinated configurations results in complex vibrational signal in the 520-800 cm−1 region. The calculated PVDOS is somewhat similar to the B–B–P cluster, with the notable exception of the missing 850 cm−1 as this model does not have B–B bonds. These results indicate that even relatively simple B–P complexes can induce rather complex IR signals in Si NPs.

As the vibrational density of states cannot be directly measured, it is beneficial to calculate the IR spectra, in order to provide a basis for comparisons with the experiments. Since bulk Si is a nonpolar crystal, we expect the IR signal of pristine Si NPs to be rather weak. Our calculations indicate that the IR spectra of the Si core do not vanish, but it is masked by the signal of the Si–H bonds. As a consequence, we subtract the contribution of the H atoms and the outermost layer of Si atoms for a better approximation of realistic Si NPs. Figure 6 shows the case of a single B–P pair introduced into two different sized host NPs (1.6 nm, 2.1 nm). It is apparent that a small peak in the 600-700 cm−1 region remains even after the subtraction of the contributions from the NP surface. The relative intensity of this small peak is expected to decrease with increasing NP size, similarly to the VDOS (see Fig. 2).

FIG. 6.

Calculated IR spectra for two different sized Si NP models, with and without a B–P pair. The contribution of the surface (H atoms and the outermost Si layer) is subtracted in order to relate the calculated spectra to the experimental data. Black, green, and orange curves correspond to the contributions of the Si atoms (without the outermost layer), B atoms, and P atoms, respectively. Note that the B- and P-related contributions to the IR spectra were scaled up for the sake of clarity.

FIG. 6.

Calculated IR spectra for two different sized Si NP models, with and without a B–P pair. The contribution of the surface (H atoms and the outermost Si layer) is subtracted in order to relate the calculated spectra to the experimental data. Black, green, and orange curves correspond to the contributions of the Si atoms (without the outermost layer), B atoms, and P atoms, respectively. Note that the B- and P-related contributions to the IR spectra were scaled up for the sake of clarity.

Close modal

A surprising result is that the introduction of the BP pair significantly changes the overall intensity and shape of the IR spectra in addition to the appearance of the B-related peak. In particular, the IR signal of the bulk-like region in the <500 cm−1 of the Si NP becomes an order of magnitude smaller. Furthermore, the contribution from the B-related local modes are negligible for the case of adjacent B–P pair. The case of the distant B–P pair is different, as the donor-acceptor pair (DAP) acts as a huge dipole, resulting in a significant contribution from the B-related local modes to the IR spectra.

The calculated IR spectra for the other model structures are plotted in Fig. 7. The IR signal of the threefold coordinated B-atom is rather weak [notice the 10× factor in Fig. 7(a)]. There is also significant difference between the IR intensities of P–B–B and B–P–B complexes. For the B–P–P cluster [Fig. 7(b)], the IR peaks related to the stretching and rocking modes of the B–B bond (at ≈850 cm−1 and ≈700 cm−1, respectively) are both relatively stronger compared to the IR intensity related to B–Si vibrations in the 600-635 cm−1 region.

FIG. 7.

The calculated IR absorption spectra for the depicted Si NP models containing a few B/P atoms decomposed for the different elements. In the case of Si, the outermost layer of Si atoms is excluded from the projection.

FIG. 7.

The calculated IR absorption spectra for the depicted Si NP models containing a few B/P atoms decomposed for the different elements. In the case of Si, the outermost layer of Si atoms is excluded from the projection.

Close modal

For the B–P–B cluster [Fig. 7(c)], the intensity of the threefold coordinated B atom’s local modes (at ≈730 cm−1) and the vibrations related to the bonds between the inner B atom and the neighboring Si atoms and P atom (≈590-670 cm−1) contribute the most to the B-projected IR intensity.

The case of the ring-like BP cluster [Fig. 7(d))] introduces a broad band in the 600-700 cm−1 region that is localized on B–Si vibrations. The model with a large BP cluster introduced into the center of a Si NP [Fig. 7(e)] has a rather broad IR band in the 520-720 cm−1 region. The absolute intensity of this band is stronger compared to the other model systems; however, the IR intensity normalized for the dopant concentration is rather small.

Finally, we discuss the case of the BP cluster introduced to the NP surface [Fig. 7(e)]. The threefold coordinated B atoms once again induce a distinct peak in the calculated IR spectra (≈780 cm−1), while the bonds between the inner B atoms and their neighboring P and Si atoms are manifested as a broad band in the 590-680 cm−1 region.

We also calculated the Raman and PL spectra of these model structures and plotted them in Fig. 8. The PL spectra of the various models do not contain substantial information about the structural properties. The Si NP models with a single B or P dopant resulted in relatively sharp emission spectra, but the spectra of the rest of the models could be only distinguished by the energy of the PL peak. For a single DAP, the absorption and emission energies depend on the donor-acceptor separation distance, resulting in lower emission energies for the NP model where the B and P are relatively far from each other. The threefold coordinated B atoms on the surface significantly reduce the emission energy for each Si NP model where they occurred [(e), (f), (i)] which is also in contradiction with experimental results. The PL spectra of the Si NP model containing a small (g) and larger (h) BP cluster is rather similar to the case of a single BP pair.

FIG. 8.

Calculated resonance Raman and PL spectra for co-doped Si NP models depicted in Fig. 5. In case of the Raman spectra, black and green lines correspond to the full Raman spectra and the Raman spectra projected to the atoms in the bulk-like region of the nanocrystal, respectively. For the PL spectra, black and red lines correspond to spectra calculated for 0 K and 300 K temperatures, respectively.

FIG. 8.

Calculated resonance Raman and PL spectra for co-doped Si NP models depicted in Fig. 5. In case of the Raman spectra, black and green lines correspond to the full Raman spectra and the Raman spectra projected to the atoms in the bulk-like region of the nanocrystal, respectively. For the PL spectra, black and red lines correspond to spectra calculated for 0 K and 300 K temperatures, respectively.

Close modal

By contrast, the Raman spectra calculated for the different Si NP models differ significantly, indicating that Raman spectroscopy is a more useful experimental tool for structural characterization. In order to differentiate between the dopant related and surface related signals, we utilize the projection technique introduced earlier. This is important as we wish to model larger Si NPs with our 1.6-2.4 nm Si NPs (see Fig. 2).

The Si NP with a single P dopant shows qualitatively similar Raman spectra to the undoped Si NP, while the models containing a single B–P pair or a single B atom show the characteristic B-related Raman peak around 600-640 cm−1 after the surface atoms are projected out.

The Raman spectra of the Si NP model with a small, ring-like B-P cluster is also quite reminiscent of the aforementioned case, but the B-related band is more diffuse with a prominent B–P peak at ≈550 cm−1.

On the other hand, the Raman spectra of the Si NP model with a large BP cluster in its center (h) differs significantly from the experiments. The B-related band extends from 500 to 720 cm−1, and the shape and relative intensity of the TO-peak is very different from that is observed in the experiments.

Finally, we discuss the case of Si NP models that contain threefold coordinated surface B atoms. The model with adjacent B-B atoms (e) is exceptional as the 850 cm−1 stretching mode of the B–B bond dominates the Raman spectra. The lack of such signals in Raman measurements could be a good indication that this configuration is indeed very rare or does not occur at all.

On the other hand, the calculated Raman spectra of model (f) seems compatible with the experimental results with a strong TO-peak and a rather weak and broad B-related IR band between 580 and 730 cm−1.

Finally we turn our attention towards the model with a surface-bound B–P cluster (i). The calculated Raman spectra for this model shows the usual TO-peak and B-related band, but exhibits two addition bands around 520-560 cm−1 and 740-780 cm−1. The latter can be identified as the characteristic local mode of the threefold coordinated B-atoms at the surface, while the former is a collective mode involving the 3 P atoms, 3 fourfold coordinated B atoms, and some neighboring Si atoms to a lesser degree.

We note that the B and P concentration of these models differs from the experimental co-doped Si NPs. Models (a)–(f) have lower while model (h) has higher B/P concentration compared to experimental samples. This means that the relative intensity of the TO-peak and the B-related band should be scaled accordingly when compared to experiments. One other question is that how these different characteristic dopant configurations interplay in the same Si NP. While the VDOS and IR spectra are expected to be superposition of the characteristic atom clusters’ signals due to the local nature of the characteristic vibrations, the situation is more complicated for Raman spectroscopy. Resonance Raman spectroscopy probes the vibrational properties by projecting the normal modes to the geometry distortion between the ground state structure and the structure of excited electronic states. The relaxation of the structure upon electronic excitation depends on the nature of the electron-hole pair induced by the absorbed photon; thus, the prediction of resonant Raman spectra for systems that contain qualitatively different dopant clusters is not straightforward. It is possible that the induced electron-hole pair (at a given excitation wavelength) becomes localized on the bulk region of the Si NP; thus, the measurement becomes insensitive to the surface region (or vice versa).

Next, we analyze the PVDOS of a couple of randomly generated B and P co-doped Si NPs. The size of the NPs is once again 1.6 nm, and we substituted 5-5 of the 147 Si atoms by B and P atoms, respectively. Our previous study demonstrated54 that the formation energy and optical properties of co-doped Si NPs are greatly influenced by the distribution of dopants even at a given doping concentration. Because of this, we have chosen the 5 most stable dopant configurations from our randomly generated ensemble of 395 co-doped Si NP models.54 In addition, we generated another 5 energetically stable configurations utilizing the empirical formula we established between the formation energy and mean dopant-dopant distances within the NP. Figure 9(a) show the calculated PVDOS averaged over the 10 different dopant configurations.

FIG. 9.

(a) The PVDOS and IR spectra averaged over 10 energetically stable, B– and P–codoped Si NPs. For the sake of clarity, the B– and P–projected VDOS were multiplied by a factor of 5 in the plot. The 1.6 nm sized host Si NP is H-terminated, and the concentration of B and P is 3% for both dopants. (b) A heatmap representation of the B-projected VDOS, where a second projection was applied which represents the involved B atoms’ distance from the NP surface. By integrating the 2D function in the y-direction, the original B-projected VDOS is obtained. This heatmap offers some insight about the correlation between the B-related normal modes’ energy and the B-atoms’ position within the Si NP.

FIG. 9.

(a) The PVDOS and IR spectra averaged over 10 energetically stable, B– and P–codoped Si NPs. For the sake of clarity, the B– and P–projected VDOS were multiplied by a factor of 5 in the plot. The 1.6 nm sized host Si NP is H-terminated, and the concentration of B and P is 3% for both dopants. (b) A heatmap representation of the B-projected VDOS, where a second projection was applied which represents the involved B atoms’ distance from the NP surface. By integrating the 2D function in the y-direction, the original B-projected VDOS is obtained. This heatmap offers some insight about the correlation between the B-related normal modes’ energy and the B-atoms’ position within the Si NP.

Close modal

The B-related peak between 570 and 720 cm−1 is familiar with the Raman measurements of co-doped Si NPs, with the exception that our results indicate a second, smaller intensity peak around 780 cm−1. We suspected that this IR band is related to B atoms in the proximity of the NP surface; thus, we calculated the correlation between the normal modes’ energy and the position of the involved B atoms. From Fig. 9(b), it is apparent that the smaller, high energy B-related IR band is indeed related to B atoms close to the surface. While the outermost layer of Si atoms were not substituted by dopants during the generation of these random structures, the layer below that was allowed to be substituted. According to our previous result (see Fig. 3), B atoms substituting the second-outermost Si layer introduce this higher energy band in the VDOS.

1. IR spectra

We also calculated the averaged IR intensities projected to B, P, and the Si atoms (excluding the outermost atomic layer). The B related IR peak is much more intense than the TO-peak, similarly to the case of single B–P pairs (see Fig. 6). In addition, the relative IR intensity of the higher energy B-Si vibrations is diminished compared to the VDOS, which might explain that this high energy peak does not appear in experimental IR spectra,11,13 even though there is significant evidence suggesting that the surface of co-doped Si NPs contains B in very high concentration. Alternatively, the amorphous B-rich surface layer encapsulating the crystalline core may alter the vibrational energies of surface-bound substitutional B atoms.

We note that some IR measurements11 of co-doped Si NPs do exhibit small absorption peaks around ≈850 cm−1 and ≈900 cm−1 that were assigned to Si–H vibrations. Given that the presence of Si–H bonds was established, this assignment is plausible, but it cannot be ruled out that B atoms also contribute to (one of) the IR peak(s) (see Fig. 1). IR measurements performed on aged samples also show a broad absorption band in the 800-900 cm−1 region,13,19 but the source of this signal could be possibly oxygen related.

2. Raman spectra

The calculated resonance Raman spectra averaged for the 10 co-doped Si NP models is shown in Fig. 10(b) compared to the experimental Raman spectra taken from Ref. 11. The B and P concentration for sample A was 0.80% and 0.29%, respectively, while 0.86% and 0.62% for sample B. After eliminating the Si–H related vibration from the Raman spectra, the result shows reasonable agreement with the observations.7,8,11 Notable differences between the simulated and measured Raman spectra are the underestimation of the TO peak energy, the raggedness of the B-related band, the small 790 cm−1 peak that is absent in the experiments, and the overall shift of the B-related band towards higher energies. The simulated B-related band falls around 680 cm−1, whereas the measured B-related band is at 635 cm−1 and 650 cm−1 for sample A and B, respectively. Finally, the relative intensity of the B-related band is slightly overestimated for our models.

FIG. 10.

(a) The experimental Raman spectra of two different co-doped Si NP ensembles (data taken from Ref. 11) compared to the (b) averaged full– and projected–resonance Raman spectra calculated for the 10 randomly generated Si NP models. (c) Measured PL spectra for a co-doped Si NP ensemble (black line) and individual NPs (colored lines) (data taken from Ref. 57). (d) The calculated PL spectra for the individual co-doped Si NP models (blue/orange thin lines), and the averaged PL spectra (black/red thick lines) at 0 K (blue/black lines) and 300 K (orange/red lines).

FIG. 10.

(a) The experimental Raman spectra of two different co-doped Si NP ensembles (data taken from Ref. 11) compared to the (b) averaged full– and projected–resonance Raman spectra calculated for the 10 randomly generated Si NP models. (c) Measured PL spectra for a co-doped Si NP ensemble (black line) and individual NPs (colored lines) (data taken from Ref. 57). (d) The calculated PL spectra for the individual co-doped Si NP models (blue/orange thin lines), and the averaged PL spectra (black/red thick lines) at 0 K (blue/black lines) and 300 K (orange/red lines).

Close modal

Before further discussion of the differences between the simulated and experimental spectra, it is important to highlight the known differences between our Si NP models and the experimental NPs. Our models are significantly smaller than the experimental NPs (1.6 nm vs. 7.0 ± 1.4 nm8), with higher dopant concentration. The sample size of our simulation is rather small, and only contains Si NPs with very low formation energies, while experimental samples are expected to contain less stable configurations due to entropy arguments (see Ref. 54). The underestimation of the TO peak energy may be attributed to the smaller NP size and the error of the DFT method. The broader TO peak can be explained by the rather small NP size, while the higher relative intensity of the B-related peak is the consequence of the larger dopant concentration. Due to the relatively small size of our NP models and the low formation energies of our sample [the dopants are typically accumulated near the surface because of the effective repulsion between B–B and P–P pairs,54 see Fig. 9(b)], B-atoms that are in the proximity of the surface are probably overrepresented compared to those in the larger NPs of the Raman measurements. This can explain the shift of the B-related band towards higher energies in our simulations. On the other hand, the experimental B-related band cuts off around 700-725 cm−1, that is, the local modes belonging to B atoms at the two outermost atomic layers of the Si NPs are missing from the Raman spectra. This might be explained by the amorphous dopant-rich layer encapsulating the crystalline core11 that could affect vibrational frequencies of these modes. It is possible, that the induced electron-hole pair is confined in the inner region of the Si NP; thus, the vibrational signature of B atoms that are proximate to the surface is not manifested in the measured Raman spectra.

We also note that boron clusters containing interstitial B or Si atoms could be also responsible for vibrational signals in the 600-720 cm−1 region.42,55 While we cannot disregard this scenario for Si NPs too, we note that most of these clusters introduce deep defect states inside the Si bandgap56 which fact is inconsistent with the electronic density of states (DOS) measurements performed for co-doped Si NPs56 where no such deep levels were found. On the other hand, a recent work12 studied the Auger interactions in the nonlinear excitation regime of co-doped Si NPs and concluded that the nonradiative recombination can be attributed to dopant-related deep defect states.

The lower energy tail of the broad B-related band is most likely the signature of existing B–P bonds in the NPs, as our calculations indicate that the vibrational energy of the B-P inside a Si NP ∼560 cm−1, while larger BP clusters introduce broader bands in the 550-700 cm−1 interval, covering the whole experimentally observed B-related band (see Fig. 5).

3. PL spectra

We also calculated the PL spectra for the aforementioned 10 Si NP models and plotted the results in Fig. 10(d). The shape of the individual PL peaks is rather similar for all 10 NPs, but the ZPL energies show a significant variance. The full width at half maximum (FWHM) of the calculated 10 spectra is 0.22 ± 0.03 eV at room temperature, while the FWHM of the 10 particle ensemble is 0.36 eV. In comparison, the FWHM of experimental samples grown at ∼900 °C (corresponding to 1.2 nm particle size) is 0.35-0.40 eV in the BPSG matrix and 0.5-0.6 eV after etching.10 Individual co-doped Si NPs have also been characterized recently.57 The investigated co-doped Si NPs were grown at 1000-1100 °C corresponding to ∼2-4 nm particle sizes, and the single-dot PL measurements were performed at 77 K. Kanno et al. found that the mean FWHM of single co-doped Si NPs is between 0.22 and 0.29 eV, in contrast to the 0.39-0.43 eV FWHM of the ensemble. They also found that the distribution of the PL FWHM changes with the growth temperature (or the NP diameter) with the smallest Si NPs showing the widest distribution of PL FWHM (75-360 meV), while the PL FWHM of the larger NPs show somewhat less variance (110-350 meV).

The most significant deviation of our results from the experiment is the energy of the PL peak. Our 10 NP ensemble has a PL peak around 2.4 eV compared to the 1.6 eV peak of the experimental sample from Ref. 57. As these Si NPs are larger than our NP models, the quantum confinement effect can be partially responsible for this discrepancy. Si NPs similar to the size of our models were also characterized by PL spectroscopy exhibiting emission peak around 1.8 eV.10 While it is plausible that the applied PBE0 TDDFT kernel is not sufficiently accurate, it is difficult to assess the error originating from the TDDFT approximation as PBE0 is known to underestimate the excitation energy of very small Si NPs,58 but overestimates the bandgap of bulk Si.59,60

On the other hand, it is possible that our Si NP models are too oversimplified for direct comparison to experimental data. Admittedly, our treatment of the NP surface is rather crude and encapsulating the crystalline core with the experimentally observed amorphous dopant-rich layer might lead to improved absolute emission energies. We also note that we observed a very strong correlation between the co-doped Si NPs formation energy and the optical gap in our earlier study.54 Due to the expensive nature of the calculation of the vibrational properties, we restricted our simulations to co-doped Si NPs with very low formation energies. We showed in our previous work54 that not only these low energy configurations occur in realistic ensembles due to the very low entropy of the stable configurations. As our sample selection favored dopant configurations with very low formation energies, we obtained a narrower, higher energy emission spectra compared to what we expect for an ensemble with broader distribution of formation energies. Extending our calculations to include less stable dopant configurations could potentially lead to larger variation of the FWHMs of individual Si NPs and a larger FWHM of the simulated ensemble.

Overall, our results compare well with the experiments apart from the absolute energy of the PL peak, despite being smaller than the experimentally studied ones and our simplistic treatment of the NP surface.

We studied the vibrational properties of B and P co-doped Si NPs by means of first principles simulations, in order to characterize co-doped Si NPs. In this work, we focused in the bulk region of the Si NPs, as the introduction of B–P clusters already creates a rather challenging problem even without considering the various surface related issues. To this end, we terminated the NP surface by H atoms, in order to reduce the computational cost. However, this compromise complicated our job as the rolling vibrational modes of the Si–H bonds overlap with the characteristic B-related local modes. In order to quantify the effect of dopants on the vibrational properties of Si NPs, we applied a projection technique to filter out the contributions of the Si–H bonds. This was especially important as most of our model Si NPs were rather small (1.6 nm) with larger surface/volume ratios compared to experimental samples. Despite of this compromise, we did identify quantifiable effects that are related to the proximity of the NP surface. We found that the vibrational signal of substitutional B is dependent on the dopants position within the NP. B atoms at least 5 Å from the NP surface show the same vibrational signature as in bulk Si, while the energies of the B related normal modes are shifted towards higher energies if the B atom is closer to the surface. In general, vibrational frequencies of a given cluster are expected to change with the distance between the cluster and the NP surface, which can lead to the incorrect interpretation of experimental spectra when the bulk data are applied.

We calculated the VDOS, IR, and Raman spectra of several BP clusters containing a few dopant atoms. We found that even small clusters containing three dopant atoms can introduce rather complex vibrational signals, especially in the proximity of the surface. B–B bonds and B atoms close to the NP surface induce the most distinct vibrational frequencies as they fall between the characteristic frequency of substitutional B (≈620 cm−1) and the frequency range associated with Si–O–Si bridges (>1000 cm−1). We characterized the size dependent PL and Raman spectra of dopant-free Si NPs and found good agreement with the experiments.

Based on these encouraging results, we utilized the same methodology to study the Raman and PL spectra of 10 randomly generated heavily co-doped Si NP models (where we have chosen stable dopant configurations), and found that the results are in good agreement with the experimental spectra, apart from the absolute energy of the PL peak which might be explained by the deficiencies of the PBE0 TDDFT kernel or the simplistic surface model applied in our work.

A.G. acknowledges the support from the National Research Development and Innovation Office of Hungary (NKFIH) Grant No. NN118161 that finances Visegrad Group (V4) and the Japan Joint Research Project on Advanced Materials (NamSeN project). This work was also financed by NKFIH within the Quantum Technology National Excellence Program (Project No. 2017-1.2.1-NKP-2017-00001). A.G. acknowledges the support from KIFÜ Supercomputer Center Grant No. 1090.

2.
M.
Fujii
,
K.
Toshikiyo
,
Y.
Takase
,
Y.
Yamaguchi
, and
S.
Hayashi
,
J. Appl. Phys.
94
,
1990
(
2003
).
3.
M.
Fujii
,
Y.
Yamaguchi
,
Y.
Takase
,
K.
Ninomiya
, and
S.
Hayashi
,
Appl. Phys. Lett.
85
,
1158
(
2004
).
4.
M.
Fujii
,
Y.
Yamaguchi
,
Y.
Takase
,
K.
Ninomiya
, and
S.
Hayashi
,
Appl. Phys. Lett.
87
,
211919
(
2005
).
5.
K.
Fujio
,
M.
Fujii
,
K.
Sumida
,
S.
Hayashi
,
M.
Fujisawa
, and
H.
Ohta
,
Appl. Phys. Lett.
93
,
021920
(
2008
).
6.
M.
Fukuda
,
M.
Fujii
, and
S.
Hayashi
,
J. Lumin.
131
,
1066
(
2011
).
7.
M.
Fujii
,
H.
Sugimoto
, and
K.
Imakita
,
Nanotechnology
27
,
262001
(
2016
).
8.
M.
Fujii
,
H.
Sugimoto
,
M.
Hasegawa
, and
K.
Imakita
,
J. Appl. Phys.
115
,
084301
(
2014
).
9.
Y.
Hori
,
S.
Kano
,
H.
Sugimoto
,
K.
Imakita
, and
M.
Fujii
,
Nano Lett.
16
,
2615
(
2016
).
10.
H.
Sugimoto
,
M.
Fujii
,
K.
Imakita
,
S.
Hayashi
, and
K.
Akamatsu
,
J. Phys. Chem. C
117
,
11850
(
2013
).
11.
H.
Sugimoto
,
M.
Yamamura
,
M.
Sakiyama
, and
M.
Fujii
,
Nanoscale
10
,
7357
(
2018
).
12.
R.
Limpens
,
M.
Fujii
,
N. R.
Neale
, and
T.
Gregorkiewicz
,
J. Phys. Chem. C
122
,
6397
(
2018
).
13.
H.
Sugimoto
,
M.
Fujii
,
K.
Imakita
,
S.
Hayashi
, and
K.
Akamatsu
,
J. Phys. Chem. C
116
,
17969
(
2012
).
14.
H.
Sugimoto
,
M.
Fujii
,
Y.
Fukuda
,
K.
Imakita
, and
K.
Akamatsu
,
Nanoscale
6
,
122
(
2014
).
15.
X. D.
Pi
,
R.
Gresback
,
R. W.
Liptak
,
S. A.
Campbell
, and
U.
Kortshagen
,
Appl. Phys. Lett.
92
,
123102
(
2008
).
16.
K.
Nomoto
,
H.
Sugimoto
,
A.
Breen
,
A. V.
Ceguerra
,
T.
Kanno
,
S. P.
Ringer
,
I. P.
Wurfl
,
G.
Conibeer
, and
M.
Fujii
,
J. Phys. Chem. C
120
,
17845
(
2016
).
17.
P. A.
Ronsheim
,
M.
Hatzistergos
, and
S.
Jin
,
J. Vac. Sci. Technol., B: Nanotechnol. Microelectron.: Mater., Process., Meas., Phenom.
28
,
C1E1
(
2010
).
18.
L. M.
Wheeler
,
N. J.
Kramer
, and
U. R.
Kortshagen
,
Nano Lett.
18
,
1888
(
2018
).
19.
M.
Sasaki
,
S.
Kano
,
H.
Sugimoto
,
K.
Imakita
, and
M.
Fujii
,
J. Phys. Chem. C
120
,
195
(
2016
).
20.
I.
Sychugov
,
R.
Juhasz
,
J.
Valenta
, and
J.
Linnros
,
Phys. Rev. Lett.
94
,
087405
(
2005
).
21.
R.
Kubo
and
Y.
Toyozawa
,
Prog. Theor. Phys.
13
,
160
(
1955
).
22.
M.
Lax
,
J. Chem. Phys.
20
,
1752
(
1952
).
23.
T.
Miyakawa
and
D. L.
Dexter
,
Phys. Rev. B
1
,
2961
(
1970
).
24.
J.
Sue
,
Y. J.
Yan
, and
S.
Mukamel
,
J. Chem. Phys.
85
,
462
(
1986
).
25.
S.
Mukamel
,
J. Chem. Phys.
82
,
5398
(
1985
).
26.
A.
Alkauskas
,
B. B.
Buckley
,
D. D.
Awschalom
, and
C. G. V.
de Walle
,
New J. Phys.
16
,
073026
(
2014
).
27.
G.
Thiering
and
A.
Gali
,
Phys. Rev. X
8
,
021063
(
2018
).
28.
F. A.
Savin
,
Opt. Spectrosc.
19
,
308
(
1965
).
29.
S.
Lee
and
S. C.
Lee
,
J. Chem. Phys.
96
,
5734
(
1992
).
30.
31.
D. C.
Hannah
,
J.
Yang
,
N. J.
Kramer
,
G. C.
Schatz
,
U. R.
Kortshagen
, and
R. D.
Schaller
,
ACS Photonics
1
,
960
(
2014
).
32.
F.
Furche
and
R.
Ahlrichs
,
J. Chem. Phys.
117
,
7433
(
2002
).
33.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
34.
S.
Baroni
,
S.
de Gironcoli
,
A.
Dal Corso
, and
P.
Giannozzi
,
Rev. Mod. Phys.
73
,
515
(
2001
).
35.
J. P.
Perdew
,
M.
Ernzerhof
, and
K.
Burke
,
J. Chem. Phys.
105
,
9982
(
1996
).
36.
G.
Kresse
and
J.
Hafner
,
Phys. Rev. B
47
,
558
(
1993
).
37.
G.
Kresse
and
J.
Furthmüller
,
Phys. Rev. B
54
,
11169
(
1996
).
38.
P. E.
Blöchl
,
Phys. Rev. B
50
,
17953
(
1994
).
39.
R.
Ahlrichs
,
M.
Bär
,
M.
Häser
,
H.
Horn
, and
C.
Kölmel
,
Chem. Phys. Lett.
162
,
165
(
1989
).
40.
R.
Bauernschmitt
and
R.
Ahlrichs
,
Chem. Phys. Lett.
256
,
454
(
1996
).
41.
A. S.
Barker
and
A. J.
Sievers
,
Rev. Mod. Phys.
47
,
S1
(
1975
).
42.
P.
Deák
,
A.
Gali
,
A.
Sólyom
,
P.
Ordejón
,
K.
Kamarás
, and
G.
Battistig
,
J. Phys.: Condens. Matter
15
,
4967
(
2003
).
43.

When the size of semiconductor nanocrystals become comparable with the De Broglie wavelength of its electrons, their electronic and optical properties start to exhibit strong size-dependence due to the surface that acts as a confining potential. The optical gap becomes strongly size-dependent, and the absorption cross section becomes enhanced for indirect bandgap materials. This phenomena is commonly referred as the quantum confinement effect (QCE).

44.
P.
Mishra
and
K. P.
Jain
,
Phys. Rev. B
62
,
14790
(
2000
).
45.
G.
Viera
,
S.
Huet
, and
L.
Boufendi
,
J. Appl. Phys.
90
,
4175
(
2001
).
46.
G.
Faraci
,
S.
Gibilisco
, and
A. R.
Pennisi
,
Phys. Rev. B
80
,
193410
(
2009
).
47.
Y.
Duan
,
J. F.
Kong
, and
W. Z.
Shen
,
J. Raman Spectrosc.
43
,
756
(
2011
), https://onlinelibrary.wiley.com/doi/pdf/10.1002/jrs.3094.
48.
J.
Zi
,
K.
Zhang
, and
X.
Xie
,
Phys. Rev. B
55
,
9263
(
1997
).
49.
W.
Cheng
and
S.-F.
Ren
,
Phys. Rev. B
65
,
205305
(
2002
).
50.
R.
Meyer
and
D.
Comtesse
,
Phys. Rev. B
83
,
014301
(
2011
).
51.
M.
Fujii
,
H.
Sugimoto
, and
S.
Kano
,
Chem. Commun.
54
,
4375
(
2018
).
52.
J. S.
Kang
,
H.
Wu
, and
Y.
Hu
,
Nano Lett.
17
,
7507
(
2017
).
53.
J. A.
Sanjurjo
,
E.
López-Cruz
,
P.
Vogl
, and
M.
Cardona
,
Phys. Rev. B
28
,
4579
(
1983
).
54.
B.
Somogyi
,
R.
Derian
,
I.
Štich
, and
A.
Gali
,
J. Phys. Chem. C
121
,
27741
(
2017
).
55.
N. X.
Truong
,
B. K. A.
Jaeger
,
S.
Gewinner
,
W.
Schöllkopf
,
A.
Fielicke
, and
O.
Dopfer
,
J. Phys. Chem. C
121
,
9560
(
2017
).
56.
P.
Deák
,
A.
Gali
,
A.
Sólyom
,
A.
Buruzs
, and
T.
Frauenheim
,
J. Phys.: Condens. Matter
17
,
S2141
(
2005
).
57.
T.
Kanno
,
H.
Sugimoto
,
A.
Fucikova
,
J.
Valenta
, and
M.
Fujii
,
J. Appl. Phys.
120
,
164307
(
2016
).
58.
R.
Derian
,
K.
Tokár
,
B.
Somogyi
,
A.
Gali
, and
I.
Štich
,
J. Chem. Theory Comput.
13
,
6061
(
2017
).
59.
J.
Paier
,
M.
Marsman
,
K.
Hummer
,
G.
Kresse
,
I. C.
Gerber
, and
J. G.
Ángyán
,
J. Chem. Phys.
124
,
154709
(
2006
).
60.
M.
Marsman
,
J.
Paier
,
A.
Stroppa
, and
G.
Kresse
,
J. Phys.: Condens. Matter
20
,
064201
(
2008
).