Determination of the responsivity of a predictable quantum efficient detector over a wide spectral range based on a 3D model of charge carrier recombination losses

Charge carrier transport in 3D model is simulated by using Cogenda Genius Device Simulator [13]. Figure 3 shows an illustration of the simulation structure with the intention of replicating the actual photodiode structure. Only 1/8 of the real device is used in the simulation model due to limited computational resources. Based on the actual size of photodiode, the simulation structure has the size of 5500 μm × 5500 μm. The structure includes a p-doped silicon substrate, p-doped and n-doped 100 μm wide contacts for electrical connection separated 300 μm apart. A SiO₂ layer is placed on top of the silicon substrate to create trapped surface charge area. The rear side of the silicon substrate is heavily p-doped. Detailed information about the photodiode structure can be found in [1]. The whole structure is illuminated at the corner with 1/4 of the radiant power, which represents an illumination at the center of the photodiode. A flat top uniform illumination area is calculated from the measured Gaussian-like light beam used in the experiments with an incident angle of 45 degrees generating an elliptical spot on the photodiode.

Zoom In Zoom Out Reset image size

The IQD caused by charge carrier recombination is calculated by:

$$\begin{equation}\delta (V,\lambda )=\frac{{R}_{\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{f}}+{R}_{\mathrm{b}\mathrm{u}\mathrm{l}\mathrm{k}}}{{G}_{\mathrm{o}\mathrm{p}\mathrm{t}}},\end{equation} \tag{ 2 }$$

where R_surf is the total surface recombination rate, R_bulk is the bulk recombination rate of electron–hole pairs and G_opt is the total rate of optical generation of charge carriers. The charge carriers are generated within the photodiode volume according to the illumination area and spectral attenuation of the light beam in silicon [14]. Even though measured at a 45 degrees angle of incidence, the beam direction will deviate only about 10 degrees from the surface normal inside the photodiode due to silicon's high refractive index. This corresponds to a penetration difference to normal incidence of less than 1.6% and is not corrected for in the generation distribution.

The simulation software is not capable of implementing a Gaussian illumination beam profile, so uniform illumination beam conversion is necessary to generate an approximation to a Gaussian beam. A Gaussian beam may be defined in various ways. We define throughout the manuscript 2σ as the beam radius at which the light intensity falls down to 1/e² of the beam axial value. A 2D Gaussian cross-section of the intensity profile of a light beam with intensity amplitude I_o and standard deviation σ as shown in figure 4 is described as

$$\begin{equation}f(r)={I}_{\mathrm{o}}\,\mathrm{exp}\left[\frac{-{r}^{2}}{2{\sigma }^{2}}\right],\end{equation} \tag{ 3 }$$

where r is the radial distance from the center of the beam. The volume of a 3D Gaussian beam represents the total number of photons in the beam

$$\begin{equation}{V}_{\mathrm{G}}=2{\int }_{0}^{\infty }{I}_{\mathrm{o}}\,\mathrm{exp}\left[\frac{-{r}^{2}}{2{\sigma }^{2}}\right]\pi r\,\mathrm{d}r\end{equation} \tag{ 4 }$$

$$\begin{equation}=2\pi {I}_{\mathrm{o}}{\sigma }^{2}.\end{equation} \tag{ 5 }$$

In a flat top cylindrical beam approximation, we want the cylindrical beam amplitude to match the Gaussian beam amplitude and calculate the radius such that it has the same volume as the Gaussian beam. The volume of a cylinder with radius R with the same amplitude I_o as the Gaussian light beam intensity profile is

$$\begin{equation}{V}_{\mathrm{C}}=\pi {R}^{2}{I}_{\mathrm{o}}.\end{equation} \tag{ 6 }$$

This means that the volume of a cylinder, representing the approximation of the profile of a Gaussian beam with the same amplitude, has a beam radius R relative to the Gaussian beam width σ as

$$\begin{equation}R=\sqrt{2}\sigma .\end{equation} \tag{ 7 }$$

With this approximation, we have a conversion of a flat top simulation beam to a measured Gaussian beam with any σ value. For elliptical Gaussian beams, we use the standard deviations σ_x and σ_y along the principal axes.

Zoom In Zoom Out Reset image size

In a seven reflections PQED trap detector, the first photodiode has an oxide thickness of about 301.6 nm and the second photodiode has an oxide layer of about 220.6 nm [1, 15]. The first photodiode absorbs 85.7% and 85.9% of the total radiant power at 476 nm and 647 nm, respectively. This means that the first photodiode absorbs about six times the power of the second photodiode, as also confirmed by the photocurrent ratio measurements between the two photodiodes [3]. Thus, when measuring the photocurrent as a function of bias voltage, the I–V curve is defined by the response of the first photodiode as for instance in a 700 μW laser power measurement 600 μW is absorbed by the first photodiode and 100 μW is absorbed by the second photodiode. Unless operated close to the bias voltage where saturation in the response occurs, the IQD of the first photodiode is orders of magnitude higher than that of the second photodiode. This means that it is the first photodiode that defines the structure and shape of the response in a relative I–V measurement. In our simulations, we have used the approximation that about 86% of the applied laser power is absorbed by the first photodiode and this represents the I–V measurement at each power level for the PQED.

The aim of the work is to fit a 3D simulation model to experimental data to extract the value of IQD. An illustration of a photocurrent measurement with increasing reverse bias voltage is given in figure 5. In order to perform a fit of an IQD model to experimental data, the experimental data has to be converted to an IQD value according to [16]

$$\begin{equation}\delta (V,\lambda )=1-\frac{I(V,\lambda )}{{I}_{\mathrm{i}\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{l}}},\end{equation} \tag{ 8 }$$

where I(V, λ) is the measured photocurrent of a photodiode at bias voltage V and I_ideal is the photocurrent from an ideal photodiode with vanishing IQD at a given absorbed radiant power Φ. I_ideal can be related to other parameters as

$$\begin{equation}{I}_{\mathrm{i}\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{l}}=\frac{e\lambda }{hc}{\Phi}.\end{equation} \tag{ 9 }$$

The expected measured current from a photodiode with an unknown internal loss δ(V, λ) is given by

$$\begin{equation}I(V,\lambda )=[1-\delta (V,\lambda )]\frac{e\lambda }{hc}{\Phi}.\end{equation} \tag{ 10 }$$

Zoom In Zoom Out Reset image size

In the experiment, the ideal photocurrent I_ideal is unknown and a correct representation of δ(V, λ) of equation (8) is possible only after I_ideal is determined. However, it is known that the expected IQD at high bias voltages is small. We use the maximum measured photocurrent I_max to facilitate the fitting process and a parameter describing the small difference between I_ideal and I_max is introduced by = I_ideal − I_max as shown in figure 5. Equation (8) then becomes

$$\begin{equation}\delta (V,\lambda )=1-\frac{I(V,\lambda )}{{I}_{\mathrm{max}}+{\epsilon}}=1-\frac{I(V,\lambda )}{{I}_{\mathrm{max}}}\frac{1}{1+{\epsilon}/{I}_{\mathrm{max}}}.\end{equation} \tag{ 11 }$$

By defining offset parameter ξ(λ) as ξ(λ) = /I_max, equation (11) can be rearranged as

$$\begin{equation}1-\frac{I(V,\lambda )}{{I}_{\mathrm{max}}}=1-[1-\delta (V,\lambda )][1+\xi (\lambda )],\end{equation} \tag{ 12 }$$

$$\begin{equation}=\delta (V,\lambda )-\xi (\lambda )+\delta (V,\lambda )\xi (\lambda ).\end{equation} \tag{ 13 }$$

Since ξ(λ) ≪ 1, the term δ(V, λ)ξ(λ) can be neglected relative to δ(V, λ). Equation (13) becomes

$$\begin{equation}\delta (V,\lambda )=1-\frac{I(V,\lambda )}{{I}_{\mathrm{max}}}+\xi (\lambda ).\end{equation} \tag{ 14 }$$

Equations (2) and (14) form a consistent set of simulation values and experimental data to be used in presenting the figures of the next sections. The offset has a significant effect only at high bias voltages where the IQD is low. The offset parameter ξ(λ) = /I_max can be positive or negative depending on noise level in the measurement and IQD of the PQED.

Measured photocurrents I(V, λ) at different power levels from 100 μW to 1000 μW as a function of bias voltage V are used for determining IQD from the experimental data. The wide dynamic range is used because the fitting parameters have various sensitivities at different power levels. The I(V, λ) curves of two different PQEDs are measured by two different laboratories at wavelengths of 647 nm and 488 nm. PQED 1 is characterized at 647 nm. The experimental Gaussian beam profile is circular with a 1/e² diameter of 2680 μm (σ = 670 μm), while PQED 2 is characterized at 488 nm with beam profile that is slightly elliptical with σ values of 267 μm (σ_x ) and 242 μm (σ_y ) in two principal axes. The equivalent flat top beam sizes used in the simulation are calculated from the experimental beam sizes by equation (7) as explained in section 3.2.

Figures 6(a) and (b) show experimental data and fitting results of PQED 1 and PQED 2, respectively. All simulations use materials data values for 300 K. The simulation result is fitted to experimental data by rewriting equation (14) as

$$\begin{equation}\delta (V,\lambda )-\xi (\lambda )=1-\frac{I(V,\lambda )}{{I}_{\mathrm{max}}}.\end{equation} \tag{ 15 }$$

The left-hand side of equation (15) is determined by the simulation results using fitted parameters. The right-hand side of the equation is calculated from experimental data. The fitting is a manual trial and error process based on the known sensitivity of the IQD prediction to the fitting parameters. Typically, the simulation is started by adjusting the values of fixed charge density Q_f and surface recombination velocity (SRV) to match the corner points of the highest power curves as indicated with dashed circles in figure 6. Change of the fixed charge Q_f shifts the corner points horizontally, whereas SRV shifts the corner points vertically. Then the remaining parameters defining the simulated IQD, i.e., bulk lifetime τ_bulk and substrate doping concentration p, are adjusted to all experimental curves to form a consistent set between simulation and experimental data, as different parameters dominate in different parts of the curves.

Zoom In Zoom Out Reset image size

The offset ξ(λ) is determined by matching the saturated IQD at high bias voltages between simulation and experimental data of the right-hand side of equation (15). Noise may affect the offset via = I_ideal − I_max (see figure 5). However, the results of the fitting process are insensitive to the exact value of I_max. If the data point corresponding to I_max is removed and the fitting process is carried out using the second largest experimental current value, the values of ξ(λ) and change in such a way that the left and right side of equation (15) agree. The saturation level IQD as well as values of the fitted parameters fixed charge Q_f, SRV, bulk lifetime and doping concentration remain unchanged.

Taking into account the rough beam shape approximations we needed to do, remarkably good agreement between experiment and simulation is obtained in figure 6 for both PQEDs, with radiant power varying by one order of magnitude and δ(V, λ) varying by more than four orders of magnitude. The effect of approximating the Gaussian laser beam profile with the flat top profile is considered in more detail in appendix A. It should be noticed that only one set of parameters defining the simulated IQD is used for each PQED as listed in table 1. The extracted doping concentration of the two PQEDs are 1.60 × 10¹² cm⁻³ and 1.95 × 10¹² cm⁻³. These values are close to the nominal value of 2 × 10¹² cm⁻³ from wafer manufacturer specifications [1]. The values of fitted fixed charge density Q_f are 4.5 × 10¹¹ e cm⁻² and 3.0 × 10¹¹ e cm⁻² for the two PQEDs while the experimental Q_f result of the photodiode manufacturer more than ten years ago is 6.2 × 10¹¹ e cm⁻² [1]. The same bulk lifetime τ_bulk value of 2.5 ms is obtained from the fitting for both PQEDs. PQED 1 has extracted SRV S₀ value of 5000 cm s⁻¹ while this value is 3000 cm s⁻¹ for PQED 2. An offset of ξ(λ) = −15 ppm is determined for PQED 1 at 647 nm, while the offset is found to be close to zero for PQED 2 at 488 nm.

In order to investigate the influence of beam size on PQED responsivity, a test on PQED 2 with a larger beam size than beam 1 used for the measurement result of figure 6(b) is performed. The new illumination beam 2 has σ values of 490 μm in both principal axes. This means that the diameter of beam 2 is approximately twice that of beam 1. Other parameters of beam 2 are close to those of beam 1 as described in table 2.

Effect of the beam size on the IQD is shown in figure 7. By using the same fitted parameters to the experiment as with beam 1 on PQED 2, the simulation model well predicts measurement result at 510 μW radiant power with beam 2 parameters. This supports the conclusion that our beam approximation method is valid as the simulation results are in good agreement despite a significant change in the photocurrent response with bias voltage in the two experiments of figure 7. A large beam size offers lower IQD at the same bias voltage and the IQD reaches the saturation level at lower bias voltages as compared to a small beam size. However, the illumination beam size does not affect the saturation IQD level at high bias voltages. As described in appendix A, the saturated IQD at large bias voltages is not sensitive to the beam shape either.

Zoom In Zoom Out Reset image size

To validate the simulation model, simulation with change in illumination wavelength is also examined. Measurements of PQED 1 are carried out at an additional wavelength of 476 nm at different radiant power levels from 100 μW to 843 μW in a similar way as the measurements at 647 nm in figure 6(a). The experimental beam has σ values of 630 μm in one principal axis and of 610 μm in the other principal axis. The beam sizes used for the measurement at 476 nm and at 647 nm (σ = 670 μm) are therefore quite close to each other. Hence, comparison of the experimental IQD in figures 6(a) and 8 reveals mainly the effect of wavelength variation. The fitted parameters of PQED 1 at 647 nm in table 1 are used for the simulation at 476 nm without any additional adjustment except for the change of offset ξ(λ) from −15 ppm at 647 nm to 0 ppm at 476 nm.

Zoom In Zoom Out Reset image size

The simulation well predicts the experimental data at 476 nm in figure 8 using the detector parameters extracted from the fit at 647 nm. Although the simulation of 843 μW radiant power shows a small deviation with experimental data, the results support the conclusion that the extracted parameters are defining the response of the detector. The slight deviation may be due to uncontrolled differences in beam shape and beam size between the two wavelengths. The used flat top approximation of a Gaussian beam profile reduces simulation accuracy at both wavelengths.

The IQD sensitivity to simulation parameters and simulation uncertainty of PQED 1 are shown in figure 9. In the sensitivity and uncertainty analysis, wavelength 476 nm is used as the predicted response is compared to measured data. The simulation is performed at 201 μW (solid lines) and 601 μW (dashed lines) of radiant power. The fitted parameters of table 1 are varied so that the experimental data $1-\frac{I(V,\lambda )}{{I}_{\mathrm{max}}}$ of equation (15) lie between the lower and upper envelope curves provided by the left-hand side of equation (15) in the bias voltage range where curves are influenced by the fitted parameters. In general, the IQD is more sensitive to the change of parameters at high radiant power levels and each parameter affects the IQD curve differently. In the case of figure 9, the nominal offset ξ(λ) is zero and the left-hand side of equation (15) is equal to the IQD.

Zoom In Zoom Out Reset image size

As shown in figure 9(a), the SRV changes the IQD at bias voltages above the corner point where the back surface recombination ceases to influence (highlighted by dashed circles in figure 6). The change in IQD with SRV at bias voltages just above this point is most sensitive at high radiant power. By changing SRV, the corner point is shifted vertically. SRV also influences the saturation level of IQD at higher bias voltages in the tens of ppm range. The value of Q_f affects the whole IQD curve. In contrast to SRV, Q_f shifts the corner point horizontally. Therefore, the corner point at high radiant power is used to determine the initial fit values of Q_f and SRV as mentioned in section 4.2. Bulk lifetime variation changes the IQD curve mostly at low bias voltages well below saturation of IQD and is more dominant at low power levels compared to SRV. Doping concentration shifts the whole IQD curve, similar to Q_f, except that just before saturation of IQD the relative sensitivity is much higher than that of Q_f and that at high bias voltages the saturation IQD stays unchanged. The change of the offset only influences the saturation level of the curves as shown in figure 9(e). The offset values are chosen so that the change in saturation levels caused by offset matches the change in saturation levels caused by SRV in figure 9(a).

As compared to the best fit S₀ value of 5000 cm s⁻¹ shown in table 1 for PQED 1, the value of S₀ ranges between 2000 cm s⁻¹ and 8000 cm s⁻¹ for the lower and upper envelope curves of figure 9(a). For the envelope curves, the Q_f varies between 3.5 × 10¹¹ e cm⁻² and 5.0 × 10¹¹ e cm⁻² while bulk lifetime τ_bulk takes values between 1.5 ms and 3.0 ms. The corresponding nominal values from table 1 are 4.5 × 10¹¹ e cm⁻² and 2.5 ms, respectively. The substrate doping varies between 1.45 × 10¹² cm⁻³ and 1.65 × 10¹² cm⁻³ in figure 9(d), while the nominal value from table 1 is 1.60 × 10¹² cm⁻³. Offset ξ(λ) varies between −20 ppm and +20 ppm in figure 9(e).

Combining the lower and upper uncertainty estimates of separate individual simulation parameters of figures 9(a)–(e) generates expanded uncertainty bands for the predicted IQD that are shown in figure 9(f). The expanded IQD uncertainty contributions at 476 nm from each of the simulation parameters as half of the difference between the upper and lower envelope curves in figure 9 at 6 V reverse bias are listed in table 3. Possible errors due to the TIA described in section 2 have only a minor contribution, because the effect of nonlinearities is largely cancelled in the ratio I(V, λ)/I_max in equation (15). From the simulation point of view, the position of the beam on the photodiode is not critical unless the beam approaches the active area boundaries of the photodiode. When the bias voltage is above the value where the back surface ceases to influence on the IQD (marked with a dashed ring in figure 6), the charge carriers are confined to the volume under the illumination beam and hence the losses are not sensitive to beam position. In addition, the most critical parameters in the uncertainty budget of table 3 are SRV S₀ and offset ξ and the envelope curves of these parameters in figures 9(a) and (e) are not sensitive to the values of these parameters at low bias voltages, where the back surface may affect the IQD.

Using the results of figure 9 and table 3, the saturation level IQD δ(λ) of PQED 1 at the 476 nm wavelength is (30 ± 27) ppm, where the expanded uncertainty of 27 ppm is calculated as the quadratic sum of the components of table 3. According to equation (1), this result determines the responsivity of PQED with similar relative uncertainty, because the reflectance and vacuum wavelength are straightforward to determine with an uncertainty of the order of 1 ppm.

As the S₀ parameter component defines the photodiode saturated IQD, experiments with higher sensitivity to S₀ parameter should be sought. From figure 9(a), it is seen that the dependence of photocurrent on bias voltage is much more sensitive to S₀ at high power levels. Further measurements at higher power levels than 1 mW could be performed. However, their relevance was revealed only after the experiments and analyses of this study.

The fitting results of figure 6(a) and table 1 can be used to calculate the spectrally dependent IQD of PQED 1 at 4 V bias and 100 μW radiant power shown in figure 10 with the best case and worst case scenarios. The best case scenario is defined as the set of parameters that makes the nominal IQD curve to shift to lower values and the worst case scenario is defined as the set of parameters that makes the nominal IQD to shift to higher values. Only dependence of IQD on fitted parameters in table 1 is taken into account for the best case and worst case scenarios. As seen in figures 9(a) and (e), there is strong correlation of the effects of SRV and offset at the saturation level IQD. For the purpose of calculating the worst and the best case scenarios, the full variation of the saturation level IQD in figure 9(f) can be scaled to SRV leading to the worst case and best case estimates of S₀ = 10 000 cm s⁻¹ and S₀ = 1000 cm s⁻¹, respectively, used in figure 10. It should be noted that at wavelengths below 450 nm, quantum yield must be taken into account to correctly calculate the effective IQD [17]. This correction is not implemented in figure 10.

Zoom In Zoom Out Reset image size

From figure 6(a), it can be seen that the IQD at 100 μW radiant power has saturated well below the bias voltage of 4 V. Thus, simulated IQDs at 4 V bias represent the saturated IQD of PQED 1. In the spectral range from 400 nm to 800 nm, the predicted IQD decreases with wavelength in figure 10. This is due to the decreasing optical absorption coefficient with wavelength in silicon. At short wavelengths, the saturated IQD is dominated by the recombination losses close to the surface of photodiode. More electron–hole pairs are generated at the surface of the photodiode due to the high absorption coefficient and it is basically the hole concentration that determines the recombination rate in the induced n-type layer. Therefore, the simulated IQD given by equation (2) increases at short wavelengths. When the wavelength is increased above 800 nm, the illumination beam penetrates deeper into the structure and beyond the depletion region where bulk recombination starts to dominate the losses and hence the IQD increases.

Although a large deviation from the nominal SRV value in table 1 is used in the best case and worst case scenarios, the change of IQD is relatively small at long wavelengths. The IQD values are 1.9 ppm, 1.2 ppm and 1.0 ppm at 600 nm, 700 nm and 800 nm for best case scenario while those values are 16.4 ppm, 8.3 ppm and 4.8 ppm for the worst case scenario. Calibration points traceable to a cryogenic radiometer are included in figure 10. The calibrated IQDs are (71 ± 160) ppm at 476 nm and (150 ± 200) ppm at 647 nm. The predicted IQDs agree with the calibration points within their combined uncertainties.