Singlet oxygen generation by nanoporous silicon: photoluminescence dynamics in magnetic field

Singlet oxygen generation in porous silicon (PSi) was investigated by a magneto-optical experiment. Photoluminescence (PL) quenching due to an energy transfer (ET) process mediated by an exchange interaction was monitored in the spectral range 1.4–2.5 eV and in a magnetic field of 0–6 Tesla at different levels of oxygen concentration and excitation pump power. When a magnetic field was applied, both PL recovery and, for magnetic fields below 2 Tesla and high concentrations of oxygen, an unusual additional pump power dependent quenching of the PL was observed. A rate equation model describing the behavior of PL from PSi with oxygen adsorbed at cryogenic temperatures in magnetic field was developed. The model has been expanded to cover the ET process as a function of the nanoparticle size.


Introduction
Energy transfer (ET) from photoexcited silicon nanoparticles (SiNPs) to molecular oxygen generates highly reactive singlet state species that have found application in fields as diverse as medicine, photodynamic cancer therapy [1,2] and chemical engineering, in optically activated reactors [3]. Exchange coupling between SiNP excitons and the ground triplet level of oxygen (the Dexter mechanism [4]) is believed to be the dominant process. Monitoring SiNP photoluminescence (PL) has proved particularly useful in understanding this phenomenon [5][6][7][8]. In the presence of oxygen, ET will depopulate the SiNP exciton levels and hence reduce (quench) the PL. The application of a large magnetic field at low temperature will spin polarize both the exciton and molecular triplet levels. The need to satisfy spin selection rules reduces the ET rate in this circumstance, and usually results in a partial recovery and enhancement of the PL.
The dependences of the efficiency of the process on magnetic field and oxygen concentration in terms of the dynamics of the ET and other excitation and relaxation processes were investigated in our previous work [8]. However, for simplicity, that work neglected three questions: these were the effects on the ET process of (i) the NP size distribution, (ii) the pump power and (iii) the random orientation of the oxygen molecules. Here, we extend our model to take into account all of these, together with fresh experimental data to test the new model. Importantly, the revised model reproduces some experimental observations (for instance, a surprising initial quenching of the PL with increasing magnetic field) that were previously unexplained.
In detail, two major simplifications were [8]: (i) the two singlet oxygen levels, 1 Σ and 1 Δ at E Σ =1.63 eV and E Δ =0.98 eV, respectively, were not distinguished; and (ii) oxygen molecules were assumed to be oriented along the magnetic field. The first simplification does not affect the singlet-to-triplet relaxation times as long as both radiative and non-radiative processes are grouped together; however, by introducing the possibility to excite 1 Δ to 1 Σ by optical pumping it can modify PL dynamics. This simplification also does not affect ET rates as long as PL intensity close to either E Δ or E Σ is analyzed, because in these cases ET to either 1 Σ or 1 Δ level can be neglected due to the large energy gap between these two states (see section 3.1.4). However, for the NPs emitting PL at about 1.3 eV, the ET to the 1 Σ and 1 Δ levels is comparable and cannot be neglected. Furthermore, the transition rates used in the model cannot be compared straightforwardly with the literature data for those two states. Meanwhile, the assumption that the oxygen molecules are not parallel to the magnetic field has a significant effect on PL dynamics due to the relaxation of selection rules for ET as well as for intramolecular transitions.
In the present work, (i) the singlet 1 Σ and 1 Δ levels of molecular oxygen were included in the model explicitly; (ii) the model was generalized for the oxygen molecules oriented randomly in the magnetic field (however the magnetic splitting of the orbital doublet 1 Δ, m L =±2 was not included); and, further, (iii) the ET process is now modeled as a function of NP size and excitation pump rate.
The modeling requires a large set of input parameters. Many of these are available in the literature; others have been estimated as part of the present work. Although our new analysis of the literature data (section 4.5 and the appendix) suggests a set of the parameters which are more realistic than in our previous work [8], we identify some parameters that still need clarification.

Experimental
The samples were produced in the form of porous silicon (PSi) layers (thickness ∼8 μm) on bulk crystalline substrates by conventional electrochemical etching from wafers consisting typically of p-type boron-doped CZ (100) silicon with resistivities of 1-25 Ω cm. Room temperature anodization was performed in a 1:1 solution of 48% aqueous HF and ethanol; the porosity was controlled by variation of the current (10-40 mA cm −2 ) and was typically 60%-70%. The etched layers were left attached to the substrates for better mechanical strength and were glued to a copper cold finger with heater and thermometer resistors attached. The samples were held either in a continuous-flow cryostat (base temperature ∼10 K) or a superconducting magnet in superfluid helium (base temperature ∼1.5 K, magnetic field up to 6 Tesla either parallel or perpendicular to the layers). In both systems, the cold finger could be raised to the top of the cryostat to expose the cold sample briefly to oxygen gas and it could be heated whilst in vacuum to desorb oxygen. PL was excited by a continuous wave solid state diode laser (wavelength ∼450 nm, power ∼0.4−40 mW at the sample, with a weakly focused laser spot, size ∼4 mm 2 ) and detected with an intensified CCD camera and compact single-grating spectrometer. The frequency-resolved spectroscopy (FRS) [9] has been used to examine the dynamics of PSi.

Rate equations: steady state solution
The model is based on ideas about donor-acceptor recombination that were undertaken by Dunstan and Davies [10]upon an assumption that O 2 molecules will be excited to their singlet states only by adjacent excitons. O 2 molecular vibrations are neglected.
The differences between low and high O 2 spectra, as well as between those for low and high pump power, point to competition between the physical processes (light absorption, radiative recombination, spin relaxation, and ET) that control the shape of the PL spectrum. These processes are indicated schematically in figures 1 and 2, and serve as a guide to the rate equation model we develop below. In this model, the photo-excited populations of the separate spin states of the excitons and the oxygen molecules are treated explicitly, taking into account the spin-dependence of the ET to O 2 , the radiative exciton recombination rate, the processes of thermal excitation and spin-lattice relaxation that lead to population redistribution between the spin states for a given silicon NP, and the rates of relaxation from singlet to triplet and singlet to singlet oxygen states.
3.1.1. Silicon nanoparticles without oxygen. At the low measurement temperatures necessary for magneto-optical experiments, we can assume that oxygen is not mobile and we can therefore divide the NPs into two separate populations; those which are not in contact with oxygen (figure 1(c)), and those which are (represented in figure 2). S levels was taken as 0.4429 meV, and an isotropic g-value of 2.0 was used. (b) The triplet, 3 Σ, and two singlet, 1 Δ and 1 Σ, levels of O 2 in a magnetic field  B z; (c) the ground state (bottom) and triplet exciton states, X T , of a SiNP in a magnetic field. In this model, the exchange-coupling between the two triplet states results in energy transfer between three possible pairs of states, of which one example is indicated by the vertical and curved arrows.
We write the proportion of NPs which do not have adsorbed oxygen molecules and which do not currently contain an exciton as n 0 ; excitons are created in these in one of the three triplet exciton states (index i=1, 2, 3) with equal pumping rates P 1 3 to generate fractional populations u i (all populations involved in the rate equations are fractional, from here we will omit the word 'fractional'). The photoexcited NPs can depopulate only by radiative emission (rates r 1 , r 2 , r 3 for m S ={−1, 0, 1}, respectively; assuming r 3 =r 1 ), spin-lattice relaxation to spin states lower in energy (γ ij ) or thermal excitation to spin states higher in energy by ΔE ( ). Under these assumptions, the steady-state solution of four rate equations for the populations u i (i=1, 2, 3) and n 0 yields: where F is the total fraction of NPs with bound O 2 . The low temperature PL intensity in the case without oxygen will be: The influence of non-radiative transitions is found to be negligible for T<100 K but increases substantially above 150 K [11,12]; in the latter case, the radiative contributions to decay rates r 1 and r 2 must be considered in (2).
3.1.2. Silicon nanoparticles with oxygen. We now consider the second population of NPs, those which are in contact with oxygen. We write the proportions of NPs which do not contain an exciton as n j , where j=1, 2, 3 runs over the three possible oxygen triplet states m S ={−1, 0, 1}, respectively. As above, excitons are created in these NPs in one of the three triplet exciton states (index i=1, 2, 3 corresponding to m S ={−1, 0, 1}, respectively) with equal pumping rates P 1 3 to generate coupled exciton-oxygen populations n ij . The exciton radiative recombination and spin-lattice relaxation terms are as above, and we introduce β ij , a spin-lattice relaxation and thermal excitation term between the oxygen triplet states analogous to γ ij . We must account for NPs in which the oxygen is in the singlet state and no exciton is present. This is the state of an NP after ET and before relaxation of the oxygen molecule, at which point the oxygen states have populations n S ( 1 Σ) and n D ( 1 Δ). We also must account for NPs in which an exciton has been created while the oxygen is still in the singlet state (populations w j and v j of 1 Σ and 1 Δ levels, respectively). Also, we have to consider the possibility of optically allowed excitation of O 2 from its lower singlet 1 Δ to the higher singlet 1 Σ level. Finally, we introduce the ET process which is the focus of this work through the rates t ij . For the case of an oxygen molecule which is aligned along the magnetic field, the matrix [t ij ] has a simple form in order to impose the overall conservation of spin angular momentum, D = m 0 where t is the ET rate at B=0. However, for O 2 tilted with an angle θ with respect to the magnetic field, the spin states of the O 2 triplet are mixed, the selection rules (3) are relaxed and, from Fermi's golden rule, the matrix needs to be considered as follows: , and D is the zero field splitting shown in figure 1. Subscripts Σ or Δ for the rate t and rate matrix t ij used below indicate ET to 1 Σ or 1 Δ levels, respectively (e.g. S D t t , ij ij ). Figure 1(a) shows the oxygen triplet states in a magnetic field, calculated according to (5), for two orthogonal directions: B parallel and perpendicular to the molecular symmetry axis, z.
As in the previous sub-section, we present the steady state solutions of the resulting nineteen rate equations plus the SiNPs that can recombine to emit photoluminescence (red arrows, 'PL') or can transfer energy to those absorbed oxygen molecules that are in the triplet ground state (orange arrows, 'energy transfer'). Excited oxygen molecules in the singlet state can return to their ground state (green arrows, 'relaxation') via emission of luminescence and/or non-radiative relaxation processes. The details about fractional populations n, v, w, relaxation rates r, R, γ, β, energy transfer rates t (with the corresponding indices) and pump rate P are given in the text.
condition that the total number of NPs with adsorbed oxygen remains constant.
The first sets of expressions, (6), represent the generation and loss of excitons in NPs with adsorbed triplet oxygen (block 2 in figure 2); the existence of two triplet entities gives nine possible joint spin states, so that nine equations are required where i=1, 2, 3, j=1, 2, 3 1 ) are the relaxation rates from the singlet 1 Σ (or 1 Δ) level of oxygen to un-mixed levels of triplet oxygen at  B z with m S =0 and m S =±1, respectively.
The sets (9) and (10) of three equations each represent the optical pumping and de-excitation of NPs with adsorbed oxygen in its singlet 1 Σ state (block 4 in figure 2) and 1 Δ state (block 6 in figure 2), respectively; the three equations (i=1, 2, 3) in each of the two sets arise from the three exciton states Here it is assumed that Pn D goes to block 3. The set of three ( j=1, 2, 3) equations (11) represent the generation and loss of NPs with triplet oxygen but no exciton (block 1 in figure 2).
The two equations (12) and (13) represent the generation and loss of NPs with singlet oxygen in the 1 Σ level (block 3 in figure 2) and 1 Δ level (block 5 in figure 2), respectively, but with no exciton: With these equations, we have a fully-determined system of twenty linear equations with twenty variables, however, the balancing equation (14) is important for normalizing and cannot be omitted; so any one of the above equations can be made redundant. Thus, equation (14) imposes the requirement that the total fraction of NPs with adsorbed oxygen should remain constant at F: We can sum all the exciton radiative processes in order to obtain an expression for the low temperature PL intensity I PL as follows: and this expression can be evaluated as a function of magnetic field, pump rate or particle size. As before, a zero nonradiative contribution to r 1 (=r 3 ) and r 2 is assumed.
We note that the probability of relaxation of the oxygen singlet to the triplet state (and also of 1 Σ→ 1 Δ transition) with the simultaneous recombination of a Si exciton (such as transitions from blocks 4 and 6 to block 1 and from block 4 to block 5 in figure 2) was not considered. These probabilities are proportional to the squares of the corresponding fractional populations and are smaller than the linear ones, and thus can be neglected. They also make the system of equations nonlinear which requires a more laborious algorithm to solve it. Exciton hopping (tunneling) processes in interconnected SiNPs are also ignored in the model.

Effect of randomly oriented O 2 .
In assuming that the oxygen molecules are randomly oriented within PSi, the number of oxygen molecules oriented with an angle θ with respect to the magnetic field B is proportional to sinθ. The PL intensity calculation was performed for one hemisphere of solid angle, i.e. for θ from 0°to 90°(with the step of 0.1 • ), with the corresponding normalized weighting factor of q p sin 2 .
3.1.4. NP size (PL energy) dependence of the ET rate. Zero phonon and phonon assisted contributions to the ET need to be considered. The contribution of each phonon replica must be summed for each PL energy. The ET rate t is proportional to the probability of ET from confined Si excitons in the particles, where their energy levels are in resonance with E Σ =1.63 eV, and from Si excitons in smaller and bigger particles with the energies separated by multiple transverse optical (TO) phonon energies of w  =  n n , 1, 2, 3 ...
TO from E Σ . The ET rate also proportional to the number of oxygen molecules adsorbed at the surface of a NP of a certain radius R (smaller for higher energies); the surface area is proportional to R 2 . Finally, the number of oxygen molecules is proportional to frequency of appearance of this particular radius given by the radius-distribution (RD) function ν(R). Thus: where ν(E) is the RD function with the R-scale converted to the E-scale according to the following relation between the SiNP PL peak position E in eV and the particle radius R in nm [13][14][15][16][17]: is the band gap of bulk Si at T=0 K, and the parameters a and A are constants which slightly vary in the literature-we used a=1. 39 (17): According to (18), the PL peak positions of samples 1 and 2, 1.71 and 1.60 eV, and E Σ =1.63 eV, correspond to the NP radii 2.0, 2.4, and 2.3nm, respectively. As shown by Suemoto et al [18] and further discussed by Yorikawa and Muramatsu [16,19] (see also [20]), an RD function of an ensemble of NPs of different size can be deduced from the PL spectrum of the ensemble excited at an energy E exc with an intensity I(E exc ): where C is a constant including the quantum yield and α tot is the total absorption of exciting laser light by the ensemble. At low temperature and for sufficiently large gap between E exc and the PL energy E, α tot (E) is a simple parabola as a function of ( ) -E E exc 2 [18]. The function ν(R) is usually fitted by a lognormal or a Gaussian function as commonly observed for colloid or granular particle size distributions [12,16,20,21]. A lognormal function of the following form was found to be preferable in most of the cases [12,16]: where R c is the median, σ is the standard deviation and R 0 is the minimal radius. The RD maximum (mode) is at 2 and the full width at half-maximum is given by . Figure 3(c), (d) shows that the PL quenching occurs far from E Σ to the both sides of the spectra. This is because, as mentioned above, the phonon replicas of the ET resonance are involved. The probabilities of these replicas decrease exponentially with an increasing gap between the replica Figure 3. Low temperature PL of PSi containing zero, low and high concentrations of O 2 . PL spectra of PSi samples 1 (a) and 2 (b) exposed to small and larger quantities of oxygen gas are shown for magnetic fields of zero to 6 Tesla in (a) and 0, 3, and 6 Tesla (for clarity) in (b). The spectra at B=0 T without O 2 are also shown for comparison. The sample was held in superfluid helium at 1.6 K and the PL was excited with 450nm (2.76 eV) continuous wave excitation with the power of 0.4 and 17.5 mW for samples 1 and 2, respectively. The vertical dashed lines at 1.63, 1.96 and 1.74 eV indicate the threshold energies above which photo-excited excitons in the silicon nanoparticles have sufficient energy to excite the adsorbed oxygen from its triplet 3 Σ to its singlet 1 Σ state, dimer state 2 A different formula, which describes the relation between the first energy level of confined Si exciton and the size of SiNP, is given in [12]: and the resonance energy as follows from the theory by Miyakawa and Dexter [22]: .. In our model, we approximated t 1 (E) per a single oxygen molecule by sum of multiple lorentzians: the first lorentzian positioned at E Σ (or E Δ ) and the others at the distance of ±nÿω TO from the first, with the amplitudes decreasing exponentially according to (21). For brevity, only one TO phonon, namely the energy-conserving TO phonon with energy 63 meV having the highest density of states and located at the Γ point of the silicon phonon dispersion was considered in our model for NPs emitting above as well as below the triplet to singlet transition energy thresholds. However, as reported in [6,7,23] and seen in figure 3(b), the phonon replicas associated with the momentum-conserving TO phonon with energy 57 meV located at the X point of the Brillouin zone can also contribute in the ET process; this momentum-conserving phonon was reported in [6,7,23] to be dominant for the 1st replica below the triplet to singlet transition energy. However it does not seem to be a common rule. Our results (figure 3) demonstrate the gap of 63.4 meV between the resonant dip at 1.620 7 eV and the 1st dip below it at 1.5573 eV for sample 1 and the gap of 57.2 meV between the resonant dip at 1.63 eV and the 1st dip below it at 1.5728 eV for sample 2. Thus, only sample 2 reveals ET assisted dominantly by the momentum-conserving TO phonon. This 1st replica below the resonance is interesting also from the other point of view: due to the indirect character of a Si exciton in NP, ET at one phonon below the resonance occurs as no-phonon process [24].
The resulting size (PL energy) dependent ET rate was normalized to the value ( ) Typical PL spectra at 1.6 K for two PSi samples exposed to zero, low and high O 2 concentrations are shown in figure 3. PL from samples 1 ( figure 3(a)) and 2 ( figure 3(b)) were excited with a laser power I exc of about 0.4 and 17.5 mW, respectively (see details above, section 2). The broad PL band corresponding to a wide distribution of SiNP sizes [20,[25][26][27] is observed for the cases of low and zero O 2 concentration. Zero-oxygen PL for samples 1 and 2 peaked at about 1.71 and 1.60 eV have linewidths of about 0.34 and 0.36 eV, respectively. The low O 2 spectra in figure 3(a), (b) are of a similar shape to those obtained in the absence of oxygen but are lower in intensity. The high O 2 spectra demonstrate much stronger quenching with severe deformation of the PL shape. It is not possible to measure quantitatively the concentration of oxygen adsorbed on the SiNPs but the stronger quenching of the PL (with the other conditions, e.g. pump power, kept the same) gives a clear indication that the concentration is higher than in the case of the weaker quenching. The strongest quenching of the PL occurs precisely (as seen more clearly in figures 3(c), (d)) for NPs having an exciton energy equal to E Σ =1.63 eV [28][29][30]. We note that the dip position at 1.621 eV for sample 1 is lower than that expected for gaseous oxygen, 1.63 eV 3 . The spectra show a number of other sharp downward-pointing peaks or dips part of which originate from the enhanced ET to oxygen for NPs whose exciton energies differ from 1.63 eV by energies corresponding to one or more momentum-and energy-conserving phonons (located, as mentioned above, at the X and Γ points of the silicon phonon dispersion, respectively). These phonon related dips are partly masked by the dips whose energies differ by phonon energies from 1.96 and 1.74 eV, which correspond to ET to the oxygen dimer state  3 As discussed in [6], the weak van der Waals interaction of adsorbed O 2 molecules with Si surface atoms lowers the energy of the excited states of molecular oxygen and slightly broadens the transitions. on photon energy for B=0, 3 and 6 T at T=1.6 K. An abrupt (in log scale) rise near 1.63 eV is followed by a more gentle slope prolonged until about 2.4 eV. The dimer and phonon related features manifest here as peaks. QS at zero magnetic field at 1.63 eV for sample 2 ( figure 3(d)) is about 30 times weaker than that for sample 1 ( figure 3(c)).

PL intensity at E Σ : magnetic field dependence
We extract PL intensities at E Σ , 1.621 eV (sample 1) and 1.63 eV (sample 2), from the spectra of figure 3(a), (b) and plot them in figure 4 as a function of magnetic field B, normalized to the PL intensity at B=0 T. This normalization eliminates the difficulties associated with considering absolute PL intensities and will facilitate the comparison of data from different samples.
The high O 2 PL in figure 3(a) both above and below the 1.63 eV threshold shows a much stronger recovery of intensity as the magnetic field is increased, by a factor of about three times (figure 4(a), right scale) and, unlike the case of low O 2 coverage ( figure 4(a), left scale), the recovery of the PL has not saturated up to B=6 T. The enhanced quenching at B=0.5 T is not an artefact. Rather different behaviors are exhibited by sample 2, when excited by higher pump power ( figure 4(b)): (i) the level of PL recovery in a magnetic field for the specimen with high O 2 is smaller than that for the specimen with low O 2 ; (ii) in the case of high O 2 , PL is additionally quenched in the wider range (up to B=2 T) than that for sample 1 (below 1 T) and deeper 4 ; (iii) further increasing pump power leads to decreasing PL recovery at higher fields with simultaneous shallowing the dip of additional quenching at lower fields (not shown). Figure 4 also shows results calculated using (15), in which we take a set of parameters based on the recent literature and our measurements. These are summarized in table 2. For the two sets of experimental data, we maintain all parameters at the same values, except for those associated with pump power and the ET process itself: these are P, the pump rate of Si excitons, F, which expresses the fraction of NPs with O 2 , and the ET rate t, which decreases as the number of an NPs having multiple O 2 molecules available decreases. The ratio of the pump rates used in simulation of the magnetic field dependences for samples 1 and 2 ( (the exact fitting numbers are incidental, they can be scaled together with the other rates). The parameters r 1 , r 3 and r 2 in the model can be magnetic field-dependent, however, with the lack of the literature data for these quantities, we fitted the data with field-independent r 1 =r 3 and r 2 . There is also the lack of knowledge about the actual difference between r 1 and r 2 (see more details about the justification of parameters in the appendix).
Having more than a dozen parameters which can be varied, it is possible to fit the magnetic field dependence of PL recovery with different sets of parameters. Two more experimental dependences-power dependence and size (PL energy) dependence (see sections 4.3.3 and 4.4.1 and figures 6 and 7)-allowed us to iterate those sets of parameters to get overall self-consistent agreement of simulated curves with the experiment. Thus, we get not only excellent fits of I(B) dependences for two different pump rates and oxygen coverage (figures 4 and 7(c), (d)), but also prediction of QS (figure 7(b)) and reasonable good fit of size (energy) dependence ( figure 6(b)).
We note, that (i) the I(B) for sample 2 with the additional quenching below 2 T ( figure 4(b)) can be fitted well if ¹ r r ;  (table 2), thus the system is not fully thermalized; (iii) the transfer rate t Σ at E Σ is slightly faster than that previously reported [23], however, this gives the subsequent good agreement of power dependence estimation with the experimental results; (iv) some rates, such as t Δ and -R b X , are too slow in comparison with the others, thus the system is not sensitive to their variation up to some extent.
In the simulations, the temperature also can be varied, since the field at which the PL recovery approaches saturation is sensitive to the relationship between m g B B and kT (see the previous publication [8]).  (20); the parameters are summarized in table 1. We note that our samples can be fitted reasonably well by Gaussians with sample 2 fitted better than sample 1; the Gaussian parameters are given in table 1 in brackets. In figure 5(a), the function ( )| | n E dR dE (deduced from PL as well as calculated), the radius R(E) and the normalized PL of sample 1 without O 2 are shown plotted versus energy.   [31] and was attributed to strong enhancement of surface recombination due to magnetic-field-induced confinement of photoexcited free carriers near the surface.
(b) and overall spectral shape in figure 6(b) and it is close to half the energy of the TO phonon, 63 meV. The ET rate t Δ from Si excitons in NPs, emitting PL at 1.4-2.5 eV range, to oxygen 1 Δ level was calculated as above with the same χ and corresponding energy distance  from E Δ =0.98 eV. The calculation shows that, at E Σ , the rate t Δ of the non-resonant ET process is about four orders of magnitude slower than the resonant process t Σ (table 2).
In the literature, the energy dependence of t in PSi, obtained for the range 1.5-1.85 eV at 300 K and 1 bar O 2 ambient, has been reported in [23] (points in figure 5(b)). In the range 1.63-1.85 eV it is simply proportional to ( ) R E 2 . Figure 6 shows the experimental and calculated ratio ( ) ( ) = I B I B 0 for samples 1 and 2 plotted for the spectral range of 1.4-2.4 eV for 1, 3 and 6 Tesla. Each spectral point was calculated using (15) with the energy-dependent parameters being normalized by their values at E Σ , obtained from fitting the experimental results in figure 4. It is seen from the experimental spectra in figure 6(c), (d) that the PL recovery for sample 2 turns to additional quenching beginning from some energy level (for low O 2 above 1.8 eV), whereas the ratios for sample 1 demonstrate negligible additional PL quenching at high energies ( figure 6(a), (b)). Despite some quantitative discrepancies, qualitatively the simulated spectra reveal all the main features in the spectrum such as the overall shape and the additional PL quenching for sample 2 and (almost) no additional quenching for sample 1. The discrepancy between simulated and experimental curves for low O 2 is larger than that for high O 2 in figure 6-partly because the parameters were iterated by simultaneous fitting figure 6

Comparison to experiment.
The discrepancies between the experimental and simulated spectra in figure 6 may have a variety of causes: these include neglected fast non-radiative recombination channels [23], which would impact on the normalization PL PL B 0 differently at different energies; the involvement of more than one phonon type in ET; the neglect of ET to oxygen dimer states, and the breakdown of the assumption of the number of adsorbed O 2 molecules being proportional to the surface area at low concentrations. Figure 7(a) demonstrates the PL intensity, I(P), simulated at E Σ . All parameters in (15) were kept the same as for the simulation of I(B) shown in figure 4 for high O 2 with varying P. Two curves for B=0 T and B=6 T for each sample are shown. Also, I(P) for the NPs without O 2 , i.e. simulated with the parameter F=0, is shown (thin black lines). I(P) from specimens without O 2 is typical for two-level system at steady state: linear at lower P with saturation at higher P; the kink is at P∼r, and the saturation occurs when P overwhelms r. The behavior of the system with O 2 is different: (i) the linear I (P) at lower P saturates at P higher than that for the case without O 2 ; (ii) there is one more kink between the linear dependence and the saturation. Those kink positions are conditioned by the complex relation between the pump, ET and relaxation rates involved (see section 3.1). The saturation level of I(P) (with and without O 2 ) is magnetic field dependent due to ¹ r r 1 2 (table 2). In figure 7(b) at 0 and 6 T, is shown. At higher pump rates, the QS at B=0 T ( figure 7(b), solid lines) approaches unity, which means that no oxygen-induced PL quenching  occurs, while the QS at B=6 T becomes slightly less than unity at P ∼ 5×10 4 ( figure 7(b), dashed lines) before it approaches 1 at P ∼ 10 6 . However, if QS is defined as   figure 4(a) and (b), respectively, indicate that the system is in the almost linear regime with and without Table 2. Parameters used in modeling low temperature PL at E Σ for sample 1 (2); inverse rates in seconds, zero-field splitting D in meV, and F is dimensionless. The given system of the inverse rates is scalable with no effect to the modeled PL dynamics. 10 −3 10 −3 10 −5 -10 −2 [11,12,27,[32][33][34][35][36][37][38][39][40][41][42]

Pump power dependence
Oxygen

Parameters used for the calculations
The justification of the parameters used in the calculations is given in the appendix. The final set of parameters used is summarized in table 2.

Conclusions
A rate equation model was developed to describe PL dynamics in silicon nanoparticles with adsorbed O 2 in an external magnetic field. The key feature of the dynamics is the ET to oxygen molecules by an exchange mechanism, accompanied by singlet oxygen generation. This was validated by a magneto-optical experiment.
Excellent fits of the magnetic field dependence of the PL intensity at the energy = S E 1.63 eV, resonant with the transition energy from oxygen triplet ( S 3 ) to singlet ( S 1 ) state, were demonstrated for different levels of pump power and oxygen concentration.
Phonon-assisted ET processes at the energies non-resonant with S E and D E were considered. Good agreement with the theoretical model was demonstrated for the magnetic field dependence of the broad PL band generated by the distribution of nanoparticle sizes in PSi. It was found that the PL quenching (and its partial recovery in a magnetic field) far from S E was still mediated by size-dependent non-resonant phonon-assisted ET processes competing with the relaxation and pump rates.
Importantly, this new model describes well the unusual, magnetic field-induced additional PL quenching observed (i) at the high-energy wing of the PL band at higher pump levels and (ii) in the magnetic field dependence below 2 Tesla for higher as well as at lower pump levels.
The dependence of the PL intensity at E Σ on pump power was analyzed. An excellent prediction of the PL QS for the different pump rates was demonstrated. The quenching vanishes at higher pump rates at the PL saturation level.
The same set of parameters was used to fit the magnetic field, pump and O 2 concentration dependences of the PL intensity and only those parameters directly related to the dependences, i.e. magnetic field, pump power and the ET rate, respectively, were varied.
The zero-field splitting parameter D of oxygen molecule frozen on the SiNP surface was roughly estimated from the experiment to be of 0.4476 meV.
In this work (see blue points in figure A1), a PL lifetime, τ PL , of about 1 ms at T=10 K and of about 0.1 ms at T=50 K was measured in the spectral range 1.52-1.91 eV by mean of FRS (see also [42]). The low temperature τ PL of SiNPs in PSi -which at cryogenic temperatures is equal to the radiative decay time, τ r -as measured in previous works [11,27,[32][33][34][35][36][37][38][39][40] to be from 0.25 to 10 ms, with the average value of about 3 ms (in the temperature range of 1.5-50 K and in the spectral range of 1.55-2.2 eV). In the recent work by Lüttjohann et al [41], the PL decay time of SiNPs with mean radius of 2.4 nm was measured to be about 0.2 ms at 5 K; a similar value was predicted theoretically for SiNPs of 2nm radius in [12]. The low temperature   Though a relatively strong dependence of Si exciton radiative decay time on detected PL energy (NP size) at high temperatures (reduction from 200 to 10 μs within the energy range from 1.3 to 2.3 eV at 300 K) was observed, no significant change of the decay time for PL from SiNPs versus detected energy was observed at low temperatures [6,23,26,36,38,41,65] (see figure A1). Our results from PSi for the PL detection range 1.52-1.91 eV are consistent with the literature data: a weak decrease of τ PL from 1.2 to 0.7 ms at 10 K and from 50 to 12 μs at 300 K. Therefore, following the experimental data, we kept r 1 and r 2 in our calculations (at 1.6 K) constant versus photon energy.
A.3. Si exciton spin-lattice relaxation time, γ À 1 (T 1 in conventional notation) Zwanenburg et al [66] reported a spin-lattice relaxation time γ −1 =2.8 s at B=1.85 T in a Si/SiGe dot; a longer γ −1 for smaller dots and a γ that depends on magnetic field as g µ B 5 or/and g µ B 7 ; there is also B-independent g = ms 4 1 . For electrons trapped in Si lateral quantum dots, Glavin and Kim [67] calculated γ, which arises from two contributions: (i) a one phonon process with g µ -- where E Z is Zeeman splitting; this is similar to the case of donor spin relaxation [68,69]; (ii) a two phonon process with g µ -- At B=1 T and T=2 K, the summed contribution from both one and two-phonon processes is g » -3 min 1 . This is close to the result for g -1 obtained experimentally by Honig and Stupp [70] for electron in phosphorus-doped Si of 2 min at B=1 T with g µ B 4 at higher fields. For a shallow boron-acceptor in strained Si at the temperatures below 2 K, Shimizu and Nakayama [71] calculated g -1 of 0.12-2 ms.
In our fits, we used the spin-lattice relaxation time g = -0.7 ms 1 (B-independent) for the Si exciton triplet (table 2).
A.4. Transition rates between levels of the oxygen molecule: R b À X ;1 ; R b À X ;0 ; R a À X ;1 ; R a À X ;0 and R b À a Figure A2 shows the triplet (X, S 3 ) and two singlet (a, D 1 and b, S 1 ) levels of molecular oxygen. For unperturbed gaseous ,0 For perturbed O 2 in condense phase environments, the transition rates are significantly increased [48]. At 4-6 K, the transition times between the energy levels of oxygen in solid matrices of four different nobel gases was measured to be:   [30] for condensed oxygen at 77 K. t D 1 from 7.9 to 64 μs was measured for gaseous oxygen in a zeolite and porous silica [50], 3 μs in water [51], and 3.9-15 ms in PSi [51,53]. For O 2 adsorbed on PSi at 5 K, t = D 0.5 ms 1 was reported [51]. In our fits, we used (see table 2). These values are within the wide span of the literature data. As mentioned above, these rates are included in the model not necessarily as radiative only rates.

A.5. Zero field splitting D
The zero field splitting parameter D ( l 2 0 or D 3 2 zz in alternative notations) of the spin-Hamiltonian (5) was obtained in earlier works [54][55][56][57][58][59]62] to be of 0.492 16 meV for a oxygen molecule in the gas phase. However, for condensed oxygen, slightly reduced values were obtained: 0.4429 meV for oxygen molecules formed as defects in single crystals at 26 K [61], in solid N 2 at 4 K [60] and in solid air [63], and 0.4336 meV in a N 2 matrix [64].
In our calculations, we used the value D=0.4429 meV which is the closest literature value to that roughly estimated in our experiment. In the inset of figure 4, it is clearly seen that PL varies almost linearly in the regions 1.5-4 T and 4-5.5 T. The kink (intersection of these regions) at 3.867 T is attributed to the level crossing (see figure 1(a)), at which m = D g B B . This gives D=0.4476 meV. D for the Si triplet exciton was calculated to be of 0.0055 meV [43] and has been neglected in the model.
A.6. O 2 spin-lattice relaxation time, β À 1 We were unable to find the literature value for b -1 . A magnetic field independent β −1 =0.12 ms for oxygen triplet (table 2) was used. Figure A2. Molecular oxygen. Schematics of the electronic levels (shown also in figure 1(b)) and transitions under consideration at zero magnetic field. For each transition, relaxation times for oxygen molecules in a dilute gas [48] (upper value) and in a nobel gas matrix at 4-6 K [47, 52] (lover value) are given.
A.7. ET rate At 1.63 eV, the ET rate t −1 =2.56 μs at 5 K and t −1 =17 μs at 300 K from Si exciton to a single oxygen molecule adsorbed on PSi was reported [23]. The latter results are shown in figure 5(b). It was assumed in [6,23] that about 8 molecules of O 2 were adsorbed on PSi per one NP (emitting PL at 1.63 eV) at low temperature. Given the evidence that PL above the 1.63 eV threshold is less suppressed (quenched) in our case than that of the low temperature case in [6,23], we can conclude that the number of O 2 molecules per NP in our case is less than 8 (or at least less than that in [6,23]). The values for the ET rate of t −1 =0.09 μs and 0.3 μs, which we have used in our model to fit magnetic field dependence of PL from sample 1 and 2, respectively, for high O 2 coverage in figure 4 and simultaneously power dependence of QS in figure 7, suggest about 2.56/0.09=28 and 2.56/0.3=9 O 2 molecules per NP in these cases assuming t proportional to the number of molecules per NP. Thus, our fitting for sample 2 is in better agreement with the assumption of 8 molecules per NP (or with the value of 2.56 μs per molecule) than for sample 1. The values t −1 of 960 μs and 42 μs, used to fit PL for low O 2 coverage for samples 1 and 2, respectively, suggest 2.56/960=0.003 and 2.56/42=0.06 O 2 molecules per NP. With fraction of NPs with adsorbed O 2 F=0.68 and F=0.9 for sample 1 and 2, respectively, it suggests 0.003/ 0.68=0.004 and 0.06/0.9=0.07 O 2 molecules per NP in the group of NPs with oxygen. In the ideal case it must be at least 1 molecule per NP in the group of NPs with oxygen, however, further decreasing of F leads to the worse fit, so we compromised between further reducing t and/or reducing F for the best fit of I(B) for low O 2 concentration. Also, the reported in [23] rise of t −1 from 0.32 to 0.45 μs when one TO phonon involved in the ET process implies χ=5.5 eV −1 in (21), which is about 3 times smaller than that used in our model, 15 eV −1 .
We note that the system of rate equations is scalable, thus, multiplication of all reciprocal rates in table 2 by the same factor will not effect the resulting fittings. Therefore, we can scale the ET rate down by the factor of, say, four, bringing it for sample 1 to '8 molecule' limit and still keeping the other rates within the range of the literature data (e.g. ). The ET rate from SiNP to a dimer (again, by the exchange mechanism) was proposed to be similar to that to an O 2 monomer [23], however, this process is not included in the model.
The recent theoretical calculation [72,73] of the rate of ET (due to the exchange only mechanism) from 2D excitons in a 1.2μm thick Si layer to an adsorbed O 2 molecule gives the value of = t 10 10 The other ET mechanisms, e.g. charge transfer, can also be involved [74] and they may result in the experimentally observed ET efficiency, however, we believe that the exchange mechanism is dominant and we do not consider those charge transfer related processes in our model. The good agreement of our experimental results with the theoretical model validates the exchange mechanism of ET.

A.8. NP size (energy) dependence of the absorption coefficient
In the case of SiNPs, the pump rate ( ) s µ P R , where σ is absorption cross-section per NP of radius R. Kovalev et al [75] reported σ versus the detected energy E (eV) of PL from PSi excited at 2.7 eV, which can be fit by 2 . The excitation energy level in our case (2.76 eV) was higher than that in [75], so we could expect weaker dependence of P on ( ) s R and, also, constant P gave the better fit of energy (size) dependence of PL in figure 6. Therefore, we used a constant P in our simulations.