1. Introduction

Single-photon avalanche diodes (SPADs) fabricated as In_0.53Ga_0.47As/InP heterostructures are compact and reliable sensors for single-photon detection in the near-infrared (NIR) range (up to 1.7 µm), with low dark count rate (DCR), good photon detection efficiency (PDE) and low timing jitter [1,2]. Many applications, from eye-safe light detection and ranging (LIDAR) [3] to quantum computing and communications [3], would greatly benefit from having high PDE, so the design of new InGaAs/InP SPADs is often oriented towards this aim. In this paper, we present a method for performing a reliable and comprehensive estimation of the photon detection efficiency of InGaAs/InP SPADs at different wavelengths and temperatures. We selected the most suitable models for ionization coefficients and complex refractive indexes of the different materials in the layer stack among those reported in the literature, and fitted them to empirical data. Additionally, starting from the temperature dependence of complex refractive indexes taken from the literature, we extended the incomplete datasets in order to perform optical simulations at different temperatures. Electrical and optical simulations were performed employing a commercially-available technology computer-aided design (TCAD) simulator (Synopsys Sentaurus [4]), to investigate the device behavior when varying parameters such as operating temperature, cross-section geometry and doping levels. The simulation results were subsequently processed using custom MATLAB [5] codes to estimate the PDE. Thanks to our accurate selection of parameter models and the implementation of their dependences on doping and temperature, the PDE model we developed is a valid and innovative tool for supporting the design of new high-PDE devices.

2. Device structure

Planar InGaAs/InP SPADs are typically fabricated as a separate absorption, grading, charge and multiplication heterostructure, as shown in Fig. 1(a) [2]. Low-energy photons are absorbed inside the InGaAs absorption layer, which has a small energy gap (∼ 0.75 eV at room temperature), while carrier multiplication takes place in the InP multiplication region (where the energy gap is 1.35 eV at room temperature). The active area of the device is defined by the shape of the double zinc diffusion: the high electric field needed to promote carrier multiplication by impact ionization is confined below the deep diffusion, but a wider shallow diffusion is also present to prevent premature edge breakdown at the border of the deep diffusion, and thus enhance the electric field uniformity inside the multiplication region.

 

Fig. 1. (a) Typical cross-section of an InGaAs/InP SPAD heterostructure and electric field profile. Examples of: (b) simulation of avalanche triggering probability at V_EX = 5 V; (c) FDTD simulation of absorbed photon density (normalized to its peak value) with 1550 nm photons flood illuminating the full front surface, at a temperature of 200 K.

Download Full Size | PDF

The electric field profile along the vertical direction must be high only in the InP multiplication region, while it has to be low enough in the low-energy gap InGaAs absorption region in order to reduce dark count rate due to field-enhanced carrier generation. This is achieved by interposing a highly doped charge layer between such two regions. In order to facilitate the drifting of photo-generated holes, a few quaternary InGaAsP grading layers with different molar fractions have been interposed between InP and InGaAs ones. Finally, a silicon nitride layer is deposited on top of the active area for acting as anti-reflection coating.

3. Photon detection efficiency simulation

The absorption of a single photon in the InGaAs layer induces the generation of an electron-hole pair and the electric field drifts the hole towards the anode. When the device is reverse biased above its breakdown voltage (V_BD), the electric field in the multiplication region is high enough (> 500 kV/cm) that a free carrier has a finite probability of triggering a self-sustaining avalanche, which can be easily read by simple external electronics and, eventually, converted into a digital output pulse. In order to accurately estimate the photon detection efficiency of the device, two contributions have to be considered: i) the probability that an impinging photon is absorbed inside the device, i.e., the absorption probability, P_abs; ii) the probability that a photo-generated carrier triggers a self-sustaining avalanche, i.e., the avalanche triggering probability, P_trig. Eventually, the photon detection efficiency is calculated as PDE = P_abs·P_trig.

3.1 Electrical simulations

Electrical simulations are used to compute the avalanche triggering probability inside the device at a given bias voltage and operating temperature. The simulation procedure can be divided in two subsequent steps. First, a complete current-voltage (I-V) characteristic of the device is simulated, in order to assess its breakdown voltage: a proper V_BD estimation depends on the ionization coefficients model. Then, the device is biased at a specific excess voltage (V_EX) above V_BD, to compute the avalanche triggering probability. For such a simulation, the impact ionization is disabled, so that the avalanche triggering probability is computed assuming no avalanche is ignited. Therefore, the simulation reflects the condition where a SPAD is biased above breakdown and no avalanche was triggered yet, so that the electric field is not modified by the current flow. In this condition, the avalanche triggering probability of a carrier inside the device is linked to the local electric field where it is photo-generated, as well as to the electric field it experiences while drifting towards the metal contact. To compute the avalanche triggering probability at a given V_EX and operating temperature, we followed the theory of Oldham et al. [6], after choosing the most appropriate ionization coefficient model. Furthermore, in order to properly account for photon absorption outside the high-field region, we added a model to include carrier diffusion towards the depleted region from the low field quasi-neutral regions, i.e., the Zn diffusion on one side and the InP buffer layer on the other side. Indeed, there is a non-negligible contribution from short-wavelength (< 900 nm) photons absorbed in the upper quasi-neutral region, which might generate electron-hole pairs, whose electrons might diffuse towards the multiplication region, where they can trigger avalanches. Figure 1(b) shows an example of the resulting avalanche triggering probability distribution at V_EX = 5 V and T = 200 K, including the diffusion contribution, plotted over one half of the device. Any modification in the cross-section geometry affects the electric field distribution inside the device, thus affecting the avalanche triggering probability and the PDE. For instance, changing the relative width and depth of the two zinc diffusions affects the electric field uniformity in the active area, which impacts the average avalanche triggering probability. Bidimensional electrical simulations yield more detailed and accurate simulation results compared to mono-dimensional ones, thus being necessary to assess the optimal ratio between the internal dimensions of the device, and to assess the excess bias required to achieve the desired electric field distribution inside the SPAD. It also follows that the accuracy of a PDE prediction is linked to the reliability of the electrical simulations (including the models adopted therein), and that a bidimensional simulation sheds light over details that would be lost with mono-dimensional ones, as will be shown later.

3.2 Optical simulations

Bidimensional optical simulations based on geometrical optics, such as raytracing or transfer matrix methods, are not well suited for estimating the photon absorption in InGaAs/InP SPADs, since some features of the device have dimensions comparable to the wavelength of interest (0.8–1.7 µm). So, we developed a simulation approach that employs the finite difference time domain (FDTD) method [7] to directly solve Maxwell’s equations inside the device. Complete FDTD simulations would require the whole 3D SPAD structure, thus leading to resource-hungry and time-consuming simulations. However, we verified that the results of a bidimensional optical simulation (with radial symmetry) are in good agreement with the ones obtained with a 3D simulation, at a fraction of the computational cost. The small features of the device, that require the adoption of FDTD simulations instead of geometrical optics, are well represented in the 2D device model, so using a 3D model does not add critical information to the simulation. On the other hand, the bi-dimensional cross-section of the device shows several details that would be lost in a 1D analysis, such as the zinc diffusion curvature mentioned above, or other differences linked to active area diameter modifications. The combination of bidimensional optical and electrical simulations allows to accurately estimate the PDE of InGaAs/InP SPADs with different diameters. Figure 1(c) shows the normalized absorbed photon density as computed employing the FDTD method. The device model was flood illuminated with 1550 nm wavelength photons (only one half of the device is shown) travelling orthogonally towards the device surface. Our device model also includes the complete back-end material stack of the SPAD, such as metal lines, nitrides, and oxides, to emulate a real fabricated detector even better. In the FDTD simulation results, interference patterns are clearly visible, because the device is illuminated by a plane wave as wide as the SPAD itself, and its diffraction through the active area aperture (between the metal contacts) causes such patterns.

3.3 Combination of electrical and optical simulation results

After both electrical and optical bidimensional simulations are completed, a plot of the PDE as a function of wavelength, at a given temperature and excess bias voltage, is obtained by exploiting a custom MATLAB code. The absorbed photon density and avalanche triggering probability are both imported as matrixes in the MATLAB environment. Then, we perform an element by element multiplication in every cell of the mesh grid, thus obtaining the local photon detection efficiency in every point of the device. Finally, the PDE of the device is evaluated as the integral of the local PDE over the whole device area. This approach allows to easily estimate the contribution of each region of the device to the total PDE by simply performing the bidimensional integral only over the region of interest. A process map illustrating the main steps of the simulation work, from to literature research to the elaboration of simulation results using MATLAB, is shown in Fig. 2.

 

Fig. 2. Illustration of the workflow for the proposed PDE model.

Download Full Size | PDF

3.4 Model construction and parameter calibration

The simulation of InGaAs/InP SPAD photon detection efficiency relies on the choice of electrical and optical models for the involved materials. In this work, we selected the best fitting models among those available in the literature, based on the comparison between simulation and experimental results. For the InP ionization coefficients, the most suitable one was selected and adopted in the simulation workflow. The minority carrier diffusion model was extracted from the literature and fitted to experimental data. The missing data for the optical models were obtained by linear interpolation of the available ones, as this approximation proved to yield accurate results.

3.4.1 InP ionization coefficients

Concerning electrical simulations, the InP ionization coefficients play the most critical role, since carrier multiplication takes place in the InP layer. Among the various models for the ionization coefficients of InP reported in the literature [8–10], we selected the one described in [10]: it includes the temperature dependance, it was extracted from samples with a structure very similar to our devices and grown in similar conditions, and it is verified for electric fields up to 550 kV/cm. Then, it was validated by comparing the simulated I-V curve to measured ones. At 225 K, the measured breakdown voltage is typically around 66 V, while simulations results with [8,9] and [10] were 71.8 V, 72.5 V and 67.5 V, respectively. The simulated temperature dependence was roughly the same for all datasets, giving a breakdown temperature coefficient similar to the measured one, i.e. ∼ 0.1 V/K. The main effect of a temperature variation is a shift of the breakdown voltage, while the bidimensional profile of avalanche triggering probability is almost unchanged.

3.4.2 Minority carrier diffusion

In order to account for carriers generated by photons absorbed outside the depleted region, which diffuse towards the high-electric field region, minority carriers’ lifetime (τ_e) in the p-doped InP region must be adequately estimated, as well as its dependence on zinc doping concentration. This is especially relevant for wavelengths shorter than 900 nm: since the upper InP layer is more than five times thicker than its absorption length (∼ 300 nm) in such wavelength range, most of the short-wavelength photons are absorbed in the InP top layer, mainly in the quasi-neutral region, while very few are absorbed in the underlying layers. Consequently, the PDE at short wavelengths is mainly due to the absorption of photons in the Zn-doped InP region. So, τ_e can be extrapolated from PDE measurement in this short-wavelength range. From our PDE measurements, we estimated that τ_e is around 500 ps at 200 K, in good agreement with the data reported in [11], which show that minority carrier lifetime in Zn-doped InP is below 500 ps for p-InP doped higher than 4·10¹⁸ cm⁻³, a doping level similar to the one present in the double zinc diffusion in our devices. Furthermore, we used additional data from [11] to compute the electrons’ lifetime as a function of doping inside the device, and we included in the TCAD model the doping profiles, for both the shallow and the deep diffusions, as measured by secondary ion mass spectrometry on purposedly-designed test devices. We fitted the doping-dependent electrons’ lifetime model to our experimental data to include it in the electrical simulations. The simulation results were then exploited to compute the avalanche triggering probability of a photon absorbed in the quasi-neutral region with a custom code. Regarding the bottom InP buffer layer, even if diffusion of carriers from the InP buffer quasi-neutral region were considered, their contribution to PDE would be negligible, because photons at wavelengths shorter than the InP cut-off are all absorbed in the shallower layers before reaching the InP buffer.

3.4.3 Optical parameters

The complex refractive index for both InP and InGaAs have been measured and reported in the literature in the past decades. However, datasets are not always complete, since they only span over a limited wavelength range and, in many cases, values reported by different authors do not agree in the overlapping regions. This is because measured values depend on the specific measurement technique, sample preparation and surface conditions. Moreover, reported data are typically taken from samples at room temperature, whereas InGaAs/InP SPADs are usually operated at about 200–240 K, so the selected model had to be adapted for temperature variations according to the temperature dependence taken from different sources, as explained in the following.

For InP complex refractive index, we used the data reported by E. D. Palik in [12], which is a compendium of optical constants data. The author compared the results of several different sources and reported the most reliable data, considering the sample preparation and measurement conditions, as well as the accuracy, precision, and resolution of the adopted measurement technique. A plot of the complex refractive index of InP at room temperature is reported in Fig. 3 (left).

 

Fig. 3. Complex refractive index of InP at room temperature (left) and temperature dependence of absorption coefficient around cut-off wavelength (right).

Download Full Size | PDF

However, since the data reported in the handbook are referred to samples at room temperature, we had to integrate it with additional data taken from [13], which show the trend of the absorption coefficient of InP around its cut-off wavelength at different temperatures. Temperature affects the energy gap of the material: at lower temperature, the energy gap is increased, implying that only photons with a higher energy (i.e., shorter wavelength λ) can be absorbed and generate an electron-hole pair. The impact on the absorption coefficient of the material (which is obtained from the imaginary part k of the complex refracting index as: $\alpha = 4 \cdot \pi \cdot k/\lambda $) can be seen as a shift of the cutoff knee towards shorter wavelengths of about 2 nm every 10 K (see Fig. 3, right). Therefore, only the absorption of photons whose energy is directly comparable to the energy gap is affected by such bandgap narrowing. As a result, the temperature dependence of the absorption coefficient is accountable for by shifting the cutoff wavelength of the spectrum towards shorter wavelengths for lower temperatures, while we kept the remaining part of the absorption coefficient spectrum unchanged. With this observation, the temperature dependence of InP absorption coefficient can be included in the model by exploiting the measurement results reported in [13]: as a result, we were able to estimate the absorption coefficient as a function of wavelength at a given temperature. As an example, Fig. 3 (right) illustrates the position of the absorption coefficient cut-off at 225 K computed by interpolation of the ones at the available temperatures (298 K and 77.4 K), which we then used to update the complete spectrum taken from [12].

Regarding InGaAs, several different authors reported the dependence of the complex refractive index over the wavelength range between 0.8 and 1.7 µm [14–17], extracted from data obtained either by measurements or by simulations (see Fig. 4, left), at room temperature. We chose the model proposed by Adachi [14], who reported a generalized method for the calculation of the optical constants of III-V compound semiconductors, since the simulation results yielded by Adachi’s dataset closely matched our PDE measurement at 200 K, except for the long-wavelength range (namely above 1500 nm), which is the one most affected by the temperature variation. With the same approach earlier described for InP, the complex refractive index model was adapted in order to account for temperature variations by integrating it with measurements (taken from [15]) at different temperatures for the range of wavelengths around the cut-off (see Fig. 4, right). In the Zielinski dataset [15], for temperatures between 162 K and 300 K, the dependence of the cut-off wavelength on temperature is almost linear, thus allowing us to perform optical simulations at temperatures not originally included in the datasets by estimating the points at the target temperature with a simple linear interpolation. However, as shown in Fig. 4 (left), the absorption coefficient reported by [14] is quite different from the one reported by [15] at 300 K in the wavelength range where they overlap. This is due to the difficulties in estimating the exact absorption coefficient profile around the cut-off wavelength, either by measurements or simulations. We verified that none of them gave optimal match with measurement data for the overlapping wavelength range. So, we chose to minimize the amount of interpolation points required to have a full dataset by exploiting the Adachi dataset up to its intersection with the Zielinski one. As a result, the simulated PDE at long wavelength shows some unsolved discrepancies with the simulated one, as will be discussed in the following. As an alternative, dataset from [16] also yielded similar results, although the PDE was slightly overestimated for wavelengths shorter than 1.2 µm. Reference [17] reported an absorption coefficient higher than the others, so it is expected to overestimate PDE even more.

 

Fig. 4. Absorption coefficient of InGaAs at room temperature according to different sources (left) and temperature dependence around cut-off wavelength (right).

Download Full Size | PDF

Finally, we calculated the complex refractive index of the quaternary InGaAsP layers as a function of their individual molar fraction, with a linear interpolation between the ones of InP and InGaAs, as it is reported to be a good approximation in [18]. Although the overall thickness of all the grading layers is, at most, 10% of the absorption layer thickness, electric field in the grading layers is still significatively high, therefore avalanche triggering probability in the grading region is high as well. For this reason, the contribution of the grading layer to the overall PDE cannot be neglected.

4. Results

Figure 5 (left) shows the measured photon detection efficiency as a function of wavelength of a 25 µm-diameter SPAD operated in gated mode with V_EX = 5 V at T = 200 K, compared to the PDE model resulting from the simulation of a nominally identical SPAD at the same temperature and bias conditions. The measurement was performed using a custom set-up where the source is a broadband quartz tungsten halogen lamp, which is coupled to a monochromator, a filter wheel and a narrow slit to select the output wavelength (line width is about 10 nm). Then, an integrating sphere receives the monochromatic light at its input port and has two output ports (0.5-inch diameter): the first one is used to monitor the optical power with a calibrated optical power meter, while the SPAD is placed at 10 cm from the second port. In this way, the actual optical power at each wavelength is constantly monitored. Moreover, given the geometry of the setup, with good approximation the photons are impinging the active area of the SPAD almost orthogonally and the optical beam is uniform (within 1%) over a diameter of 6 mm at the location where the SPAD is placed. Therefore, the measurement conditions are quite similar to what represented in the simulation environment, where the device is flood-illuminated by a plane wave and all photons imping the detector orthogonally.

 

Fig. 5. Measured vs. simulated PDE (left); individual contributions to PDE from the different regions (right). In both cases, flood illumination was employed, with T = 200 K and V_EX = 5 V.

Download Full Size | PDF

Despite some minor local inaccuracies, the PDE model is in good agreement with the experimental measurement in the wavelength range between 950 nm and 1550 nm. The main discrepancies are present around 900 nm and 1600 nm wavelengths, which correspond to the cut-off wavelengths of InP and InGaAs at 200 K, respectively, and are a critical aspect of the model, where even a small inaccuracy has a great impact on the results. Therefore, these differences may be ascribed to a slightly imprecise estimation of the absorption cut-off wavelengths for both InP and InGaAs. However, the main trend of the PDE curve is well represented, and the cut-off wavelength is correctly located at about 1600 nm at 200 K, so the choice of the models and their extension through linear interpolation has proven to be a valid approach for the simulation of PDE. It must be noted that the matching between simulations and measurements can be improved by deeply revising the models, which goes beyond the scope of this work.

Additionally, since the main contribution to PDE in the 0.8–1 µm wavelength range is carrier diffusion, simulation results may be further improved with a finer tuning of the diffusion parameters, that may be pursued in the future by performing targeted measurements on specifically designed devices. Anyhow, at 1550 nm, which is a key wavelength for many applications, the simulated PDE (26.7%) is very close to the measured one (25.7%).

In Fig. 5 (right) we also show the individual contribution of each of the main device regions to the overall PDE curve, computed with our simulation model.

As expected, for wavelengths shorter than 900 nm, the main contribution is given by the diffusion of carriers generated upon photon absorption inside the InP quasi-neutral region, namely within the zinc diffusion.

The depleted InP multiplication region also contributes to PDE with the absorption of photons at wavelengths right before the InP cut-off (∼ 900 nm): such photons have lower energy with respect to the ones at 800 nm, therefore they can travel more deeply through the top InP layer, and are absorbed below the p-n junction, in the depleted region. In accordance with this analysis, it can be observed that the contribution from the depleted InP layer increases with wavelength, an opposite trend with respect to the one shown by the diffusion contribution curve of the Zn doped region. However, the increasing trend of the curve referred to InP-top abruptly stops at the InP cutoff wavelength.

As previously mentioned, the quaternary grading layers also contribute appreciably to the overall PDE throughout the whole spectrum, thus causing a vertical shift of the PDE curve almost at all the wavelengths. Indeed, our simulations showed that, in the wavelength range between 950 nm and 1550 nm, the photons absorbed in the grading layers contribute to about 10% of the total PDE.

The main contribution to PDE is given by the InGaAs absorption layer, as expected, with a sharp cut-off around 1600 nm.

As another example of the application of our PDE simulation model, we show the simulated PDE of two different devices, with active area diameters of 25 µm and 10 µm, designed for operation at 275 K, compared to the measured one (Fig. 6). As in the previous example, the simulation results well match the measured ones, proving to be a valid tool for the estimation of PDE before fabrication. The only exceptions are the regions around the cut-off wavelengths, where mismatches are more pronounced for the same reasons reported above. Nevertheless, the position of the cut-off wavelength is correctly estimated.

 

Fig. 6. Comparison between simulated and measured PDE of two SPAD devices designed for operation at 275 K, with an active area diameter of 25 µm (left) and 10 µm (right).

Download Full Size | PDF

Oscillations in the 0.9–1.5 µm range are present for both the simulation and the measured data, but do not always match. This might be due to the exact composition of the Si₃N₄ antireflection coating, with the fabricated one being slightly different from the nominal one, possibly due to either a fabrication inaccuracy or the presence of impurities. To simulate this, as an example, we computed the PDE curve of a device employing a complex refractive index of Si₃N₄ deviating of ± 10% with respect to nominal value (see Fig. 7, left). The simulation result suggests that a variation in the composition of the antireflection coating does affect the local peaks of the PDE curve for wavelengths between 0.9 and 1.4 µm.

 

Fig. 7. Left: simulated PDE with SiN_x variations. Right: PDE of devices with different diameter, for a comparison between bidimensional and mono-dimensional simulations. In both plots, we used T = 200 K and V_EX = 5 V.

Download Full Size | PDF

Figure 7 (right) shows the PDE from bidimensional simulations of SPADs with a diameter of 5, 10, 25 and 50 µm, compared with the one resulting from a one-dimensional simulation in correspondence of the device center. The bidimensional simulations predict an increase in PDE as the diameter is decreased. Indeed, the peak avalanche triggering probability in correspondence of the deep diffusion edge is the same for all diameters, but for a smaller diameter the edge region spans over a greater portion of the device, with respect to a larger diameter, therefore the average avalanche triggering probability is higher. On the other hand, a mono-dimensional simulation only shows the value of avalanche triggering probability in a section of the device, so all information about edge effect is lost. Moreover, shrinking the diameter affects the diffraction of the impinging electromagnetic waves, which appears as a variation of the oscillations pattern and magnitude in the PDE plot; a mono-dimensional optical simulation, instead, is performed with raytracing rather than FDTD, so no oscillations are visible.

5. Conclusions

Starting from optical and electrical models taken from the literature, we presented a comprehensive model for the estimation of the photon detection efficiency of InGaAs/InP SPADs at different temperatures, aimed at supporting the design of high-efficiency detectors. Exploiting bidimensional optical and electrical simulations, our model includes the temperature dependence of the ionization coefficients and complex refractive indexes of all the materials employed in the heterostructure, as obtained after an extensive bibliographic research, expanded with the linear interpolations needed to overcome the lack of data at different temperatures. We also included a doping-dependent model for minority carrier lifetime fitted to experimental measurements, in order to account for the contribution given by the diffusion of carriers generated by short-wavelength photons. Bidimensional simulations, performed on accurate models of the devices (which also include all the backend stack and doping profile), highlight the distribution of the avalanche triggering probability and photon absorption probability throughout the device, thus allowing us to compute the PDE with higher accuracy with respect to one-dimensional models. The simulated PDE showed good agreement with the measured one, proving that our model is a valid and versatile tool to perform reliable simulations of InGaAs/InP SPADs. The accuracy of the simulation results may be improved by developing custom models for all the parameters involved, thus requiring the fabrication of specifically designed devices. As this remains a possible future development beyond the present work, we believe that the method presented here, based on data available in the literature, proved to properly simulate the PDE of InGaAs/InP SPADs.