Modeling a Thick Hydrogenated Amorphous Silicon Substrate for Ionizing Radiation Detectors

There is currently a renewed interest in hydrogenated amorphous silicon (a-Si:H) for use in particle detection applications. Whilst this material has been comprehensively investigated from a numerical perspective within the context of photovoltaic and imaging applications, the majority of work related to its application in particle detection has been limited to experimental studies. In this study, a material model to mimic the electrical and charge collection behavior of a-Si:H is developed using the SYNOPSYS©Technology Computer Aided Design (TCAD) simulation tool Sentaurus. The key focus of the model is concerned with the quasi-continuous defect distribution of acceptor and donor defects near the valence and conduction bands (tails states) and a Gaussian distribution of acceptor and donor defects within the mid-gap with the main parameters being the defect energy level, capture cross-section, and trap density. Currently, Sentaurus TCAD offers Poole-Frenkel mobility and trap models, however, these were deemed to be incompatible with thick a-Si:H substrates. With the addition of a fitting function, the model was able to provide acceptable agreement (within 10 nA cm−2) between simulated and experimental leakage current density for a-Si:H substrates with thicknesses of 12 and 30 μm. Additional transient simulations performed to mimic the response of the 12 μm thick device demonstrated excellent agreement (1%) with experimental data found in the literature in terms of the operating voltage required to deplete thick a-Si:H devices. The a-Si:H model developed in this work provides a method of optimizing a-Si:H based devices for particle detection applications.


INTRODUCTION
Currently, silicon based technology is the gold standard for the majority of solid-state microelectronics and detector technology. High quality crystalline Silicon (c-Si) is typically utilized in such devices due to its well-known characteristics and technological refinement. However, fabrication of silicon sensors requires the production and implantation of wafers which limit the use of silicon to relatively small area (cm 2 ) devices. An alternative material suitable for fabrication of large area, pixellated devices with potentially a very high radiation hardness is hydrogenated amorphous silicon (a-Si:H) which is already in use in a variety of applications from photovoltaic devices [1,2], to large area x-ray imaging detectors [3][4][5].
The desirability of a-Si:H is in part owed not only to the capability of the material to be deposited over a large area (i.e., without the need for physical tiling) but also to be deposited above a variety of different substrates including flexible materials like polyamide (kapton), opening up a myriad of potential applications in radiation detection physics [6]. Furthermore, it's becoming increasingly desirable in particle detector applications, given its low cost and superior radiation tolerance [7,8]. This radiation hardness can be accounted for in the disordered structure inherent to a-Si:H, and the passivation of delocalized states or defects through the introduction of high concentrations of hydrogen. In contrast to the ordered nature of c-Si, consisting of a well-defined electronic energy structure consisting of energy bands and energy gaps, the electronic states of the highly disordered a-Si:H material can be considered to be quasicontinuous. a-Si:H is characterized by the presence of dangling bonds (DBs) due to the condition that some Si-Si bonds are not satisfied. The introduction of hydrogen (typically 4-10% atomic minimum) is used to passivate DBs, and will have a direct impact upon the width of the band-gap (E g ) and concentration of tail states depending on both the hydrogen concentration introduced and the temperature of annealing. The distribution of defects within the band structure of a-Si:H can be described by the defect pool model developed in 1990 in order to model solar cells [9]. DBs act as recombination centers or defects within the a-Si:H material and are present as a continuous distribution of states within E g [2,8]. These states can be classified as either extended "tails" of the valence and conduction bands, or as a localized distribution of states within E g . Tail states can be considered as localized states with almost zero mobility (hopping conductivity will be observed at low temperature). Conductivity at room temperature or above goes through the extended tail states (in the valence band or conduction band) through a mechanism called multiple trapping, whereby charge carriers undergoing transport can become spatially localized/trapped and excited [10]. Doped a-Si:H can be obtained through the addition of phosphine, i.e., PH 3 (n-type) or Diborane(6), i.e., B 2 H 6 (p-type) to the process gas (e.g., a mixture of silane-SiH 4 and hydrogen-H 2 ) during the deposition process.
Given the large variance in reported material properties for a-Si:H, there is at present no established commercial model of the material available to truly mimic the electrical and charge collection behavior of the unique material within the context of ionizing radiation detection applications. The vast majority of simulation based studies dealing with a-Si:H are primarily concerned with thin film devices (i.e., nm thickness range) for photovoltaic applications. These studies have been performed with a variety of commercial and non-commercial software packages. One of the first such instances of utilizing numerical techniques to investigate the performance of a-Si:H based devices was that of Dutta et al., [11] and Despeisse et al. [12]. Dutta, utilized a custom electrical-optical model known as ASDMP to investigate possible means of improving the sensitivity of the open circuit voltage for a-Si:H based pi-n solar cells [11]. Despeisse, utilized HSPICE to simulate the generated signal within an a-Si:H based particle detector due to pulsed laser light (660 nm) showing excellent agreement with experimental results. The simulation was able to mimic the experimental behavior for both the magnitude of the induced current as well as the timing related to the drift based charge collection [12]. Best results in this study were found by setting the electron induced signal time constant τ c to be 5.6 ns and the electron mobility µ e to be 2 cm 2 V −1 s −1 . In 2010, Nawaz et al. utilized the commercial Technology Computer Aided Design (TCAD) simulation tool-kit provided by SILVACO©to investigate a-Si:H solar cells [2]. Nawaz implemented perhaps the most comprehensive description of the defect distribution within a-Si:H to date. Instead of using the standard Density Of States (DOS) models, Nawaz defines defects as a combination of single energy level defects within the "mid-gap." More recently, there has been a shift to utilizing the Sentaurus TCAD tool-kit distributed by SYNOPSYS© [13]. The Sentaurus TCAD simulation tool was developed primarily for design and fabrication optimization of semiconductor electronic devices. Given that a-Si:H is not available within the standard material libraries, SYNOPSYS©users interested in simulating a-Si:H based devices are forced to develop their own material models to mimic the material behavior in the regime of interest. Lee et al. presents one of the first documented uses of Sentaurus TCAD for simulating the performance of an a-Si:H based device and provides a comparison with AMPS-1D, already showing the functionality of Sentaurus TCAD in competing with other numerical codes such as AMPS-1D, SCAPS-1D, PC-1D, SimWindows, ADEPT-F, AFORS-HET, and ASA [14]. Sentaurus TCAD was also used in a later study by Otero et al., for the purpose of modeling the performance of a-Si:H solar cells [15]. The n-i-p structure of interest in this study, although much thinner, is very similar to that which is investigated in this work. Similar to Nawaz et al. the authors of this study considered the defects present within a-Si:H as a continuous exponential tail distribution originating from the conduction and valence bands, respectively, as well as a Gaussian distribution of mid-gap trap states. Unlike Nawaz, Otero utilizes the DOS model, which while validated for crystalline silicon has yet to be validated for a-Si:H. With the exception of the study carried out by Despeisse, the majority of simulation based studies related to a-Si:H are concerned with understanding or improving the efficiency of optically stimulated current and do not consider the operating requirements necessary for detection of ionizing particles.
The renewed interest in a-Si:H n-i-p photo-diodes for radiation sensors has prompted the investigation into numerical models to mimic the electrical and charge collection response of such devices. A key focus point in this study is concerned with appropriately modeling the various deep level defects within the a-Si:H band-gap that can act as recombination centers and/or trap states. The charge collection efficiency and charge collection timing for minimum ionizing particles (MIPs) were selected as the reference figure of merit for the model validation. The realization and application of an accurate numerical model for the a-Si:H material provides an opportunity to take advantage of the predictive capability to optimize the geometries and implantation structure of such semiconductor devices.

Material Parameters
Sentaurus TCAD is a suite of tools that have been specifically developed to simulate fabrication processes used to realize semiconductor devices as well as the operating behavior under steady-state and dynamic conditions. Sentaurus TCAD is equipped with a variety of experimentally validated numerical models to describe the relevant physical and electronic behavior of commonly used materials used in the fabrication of semiconductor devices. Unfortunately, a model to mimic the behavior of a-Si:H within the Sentaurus TCAD material library framework has yet to be realized. It is therefore, necessary to implement a custom material parameter file "aSiH.par." To that end, the standard "Silicon.par" model was modified to incorporate the unique and peculiar physical and electronic attributes of the hydrogenated form of amorphous silicon. In addition to the development of this new material parameter file, additional description of the defect distribution is required to be implemented within the "sdevice" physics module native to Sentaurus TCAD.
In order to correctly model the disordered nature of a-Si:H within the framework of Sentaurus TCAD, it is necessary to implement both acceptor (electron trap) and donor (hole trap) defects at discrete energy levels within the mid-gap. To that end, an approach similar to the two level defect model described by Petasecca et al. [16], was adopted and modified to agree with the defect model proposed by Nawaz. The energy levels of the individual defects specified by Nawaz was first modified to fit within the mid-gap of the a-Si:H material used in this study. This was achieved by multiplying each energy level by the ratio of the band-gap used in this study and the band-gap defined by Nawaz. Next, the Gaussian defect states were normalized to a DB density of 4 × 10 15 cm −2 . In order to obtain agreement with the experimental results, unique scaling factors were applied to the concentration of individual acceptor/donor defects in the tails (DTi and ATi) and Gaussian distribution of defects (DGi and AGi), where the notation (i) refers to each single defect state (see Figure 1). The exact details in terms of energy level, electron/hole cross section and concentration are provided in the Supplementary Material. The energy levels described in the following tables refer to their position within the mid-gap with respect to the valence band (VB).
Despite the significant differences within the structure of c-Si and a-Si:H, only a few parameters within the new "aSiH.par" material parameter file were required to be modified. Firstly, the electron and hole mobilities were adjusted from 1,417 and 470.5 cm 2 V −1 s −1 , to 10 and 0.01 cm 2 V −1 s −1 , respectively, in accordance with the relevant literature [2,17]. The second major change to the material parameter file was the adjustment of the band-gap and is presented in the results, see Table 1.

The Physics Model
The Device Simulation for Smart Integrated Systems (DESSIS) tool kit native to Sentaurus TCAD solves the Poisson and electron/hole continuity equations in order to simulate the electrical behavior of the a-Si:H n-i-p device. The doping dependent Shockley-Read-Hall and trap assisted Auger physics models were activated to describe all carrier generation/recombination processes. The doping dependent, Enormal and high field saturation models were activated to describe the mobility of carriers within the simulation. The Poole-Frenkel model related to traps within the mid-gap was activated in order to take into account the enhanced emission probability (Ŵ PF ) for charged trap centers, due to the reduction in the potential barrier associated with the presence of high external electric fields (E). As shown in Equation (1), the role of Ŵ PF is to adjust the initial trap cross sections from σ 0 n,p to their new value of σ enh n,p as determined by the Poole-Frenkel trap model. The enhanced emission probability (see Equation 2) can be treated as a variable through the manipulation of the tuning parameter α (see Equation 3) via the adjustment of ε PF from the default value of 11.7 [30]. It should be noted that k, T, and q in Equation (3) refer to Boltzmann's constant, temperature and electronic charge.
The Poole-Frenkel mobility model was originally designed in TCAD to simulate organic semiconductors for thin film transistors and photovoltaic devices and was considered in the early phases of this work, however, its inclusion was subsequently discontinued as it was deemed to be incompatible with the distribution of defects introduced with the defect-pool model of the DBs and thus its inability to correctly mimic the E-field ( √ E) dependence upon the current density.

Geometry of the Device
In this study, we chose to model a simplified 2D version of the n-i-p (n-doped, intrinsic, and p-doped layers) a-Si:H diode structure as described by Wyrsch et al. [8]. The geometry depicted in Figure 2, features a 90 nm thick n-type layer upon a 30 µm thick intrinsic layer upon a 90 nm thick p-type layer. The intrinsic layer is simulated in two different thicknesses to evaluate the capability of the model to scale macroscopic characteristics such as leakage current, electric field distribution and charge collection efficiency with thickness appropriately. The p-type and n-type doped layers are modeled using a Gaussian analytic function with peak concentrations of 3 × 10 +18 and 8 × 10 +18 cm −3 , respectively. The modeled device is 50 µm wide and given its 2D nature, is by default 1 µm deep. Given that the model is a 2D representation of the a-Si:H ni-p device, a geometric scaling factor A F was used to allow for a comparison between experimental and simulated leakage current density. This geometric scaling factor is used to account for the considerably different areas of the experimental (2 × 2 mm 2 ) and simulated (50 × 1 µm 2 ) device structures. A discretized meshing strategy was applied to the entire geometry in order to solve and evaluate the Poisson and electron/hole continuity equations throughout the geometry. As part of the meshing process, the doping dependent mesh strategy was employed to optimize the number of nodes in regions of interest to improve simulation accuracy and minimize computational time.

Simulation Process
After appropriate specification of the adjusted material parameters and defect distribution within the mid-gap, a preliminary steady-state simulation of the electrical behavior of the 30 µm thick a-Si:H n-i-p device was performed and compared to experimental results. The DC-simulation technique was used to investigate the leakage current of a 12 and a 30 µm thick device from 0 to 200 and 0 to 400 V, respectively. The same simulation set-up was used to investigate the electric field distribution as a function of bias for verification of the depletion region achievable in the intrinsic layer of the device. Lastly, the charge collection efficiency was evaluated by means of a transient simulation of a MIP with a linear energy transfer of 1.28×10 −5 pC µm −1 traversing the center of the device.  (4) was employed to fit the field dependent components of the data. The summation of the simulated and fitted result (green circles) compared to the experimental result is presented on the right. FIGURE 5 | Comparison of experimental (black square) and simulated (red circle) current density characteristics as a function of externally applied bias for the 30 µm thick a-Si:H n-i-p device. A fitting function (Blue triangles) described by Equation (4) was employed to fit the field dependent components of the data. The summation of the simulated and fitted result (green circles) compared to the experimental is presented on the right.

Current-Voltage Characteristics
Because the band-gap of a-Si:H is dependent upon the concentration of introduced hydrogen. A preliminary study was performed to investigate the effect of the band-gap (1.70 eV ≤ E g ≤ 1.90 eV) upon the simulated leakage current density and is depicted in Figure 3. In this study, the most appropriate value of the band-gap or E g was determined to be 1.84 eV. The appropriateness of the value of E g was determined by considering the agreement between the experimental and simulated leakage current density and is described in more detail below.
Unique scaling factors were applied to the acceptor and donor defect states ATi, AGi, DGi, and DTi depicted in Figure 1, which are 0.001, 5, 1, and 0.001, respectively. The results of the leakage current simulation by means of a static DC simulation for the 12 and 30 µm thick a-Si:H n-i-p devices are presented in Figures 4, 5. A comparison of the simulated results with experimental results provided by Wyrsch et al. is made [8]. In both cases, the model employed was not able to mimic the experimental results satisfactorily. This discrepancy is a result of the transition from drift dominated transport to Poole-Frenkel transport with its well known square root of the electric field ( √ E) dependency [31]. As the magnitude of the electric field increases, charge carriers trapped within defect centers are able to acquire sufficient energy to contribute to conduction within the device. To that end, a fitting function described by Equation (4), was employed to mimic the high field dependency of the experimental leakage current density (J). The aforementioned equation includes a modified version of the Poole-Frenkel effect which describes the √ E dependency. Note in this equation, J 0 is a dimensionality constant (A cm −2 ), E is the electric field strength (V cm −1 ), d is the substrate thickness in cm and J SRH is the current density simulated by TCAD, based upon Shockley-Read-Hall statistics, taking into account the temperature, electric field, effective carrier concentration, and generation-recombination centers. It should be noted that the Poole-Frenkel trap model previously discussed was investigated in this study. However, it was determined that its activation has no impact upon the results.
The summation (green circles) of the simulated current density (red circles) and the fitting function (blue triangles) provides excellent agreement with the experimental data (black squares) for both the 12 and 30 µm thick a-Si:H n-i-p devices (see Similar to the ideality factor in the Shockley diode equation [32,33], the fitting parameters b and J 0 are to account for variations of the a-Si:H material quality and/or fabrication processes.

Electric Field Distribution
The electric field is obtained from the gradient of the electrostatic potential. In this section the distribution of the electric field for the 12 and 30 µm thick a-Si:H n-i-p structures for different applied bias is investigated. Figure 6 depicts the cross section of the electric field distribution for the 12 and 30 µm thick a-Si:H ni-p structures, respectively. The cross section is taken through the center of the device, i.e., x = 0, in order to show how the electric field changes with respect to diode thickness. The DC results of the 12 µm thick a-Si:H n-i-p structures show full depletion for

Charge Collection Study
The "HeavyIon" model native to Sentaurus TCAD was adapted to mimic a MIP traversing the 12 µm thick a-Si:H n-i-p devices at specific time (t 0 = 10 ns) and location (x = 0). The Linear Energy Transfer (LET) for the MIP was set to 1.28×10 −5 pC µm −1 corresponding to approximately 80 pairs.µm −1 . As shown previously, the large concentration of traps is associated with a large leakage current density. It also leads to substantial instability in the transient behavior of the material. In order to calculate the charge collection associated with the ion strike, the transient instability must be taken into account.  Figure 8. In order to assess the functionality of the Poole-Frenkel trap model, the simulation was performed with and without the model activated as shown in Figure 8.
Although the simulated CCE behavior with the Poole-Frenkel trap model activated shows a predictable trend at lower biases, it displays a significant deviation from the expected behavior at higher biases, with obvious non-linear trend and lack of CCE saturation indicating that the model is not presently capable of modeling the depletion of thick a-Si:H substrates. In contrast, the simulated CCE behavior without the Poole-Frenkel trap model activated manifests the expected trend with a linear increase in CCE at lower biases and saturation in the condition of over-depletion. At low biases, the CCE is hampered by the defects within the a-Si:H material limiting charge collection through trapping processes and reduction of the carrier mean free path. Once the bias is increased sufficiently, there is an increased proportion of charge carriers able to contribute to the current, until such point that the charge collection saturates. Using the CCE behavior without the Poole-Frenkel trap model, allows an alternative method of determining the operating conditions required to fully deplete the device in order to measure the signal due to incident ionizing radiation. Plotting two linear fitting functions to the CCE behavior, one in the low bias region and the second in the high bias region where the charge collection has saturated, allows determination of the point at which these lines intersect. Determination of the intersection position on the x-axis, i.e., the square root of the bias, gives a rough approximation of the depletion voltage. Using this method, a voltage of 68 ± 20 V is given, corresponding to an electric field strength of 5.7 V µm −1 . Despeisse et al. [12] and Wyrsch et al. [34] provides an empirically based formalism to predict the depletion voltage (V D ) of a-Si:H based devices with DB density of the order of 10 14 using only the device thickness in µm (D), i.e., V D ≈ 0.48 × D 2 . According to this formalism, the substitution of our device thickness of 12 µm gives a depletion voltage of 69 V, corresponding to electric field strengths of 5.75 V µm −1 . In the case of the 12 µm thick device, the two methods provide excellent agreement (≈1%) in terms of the prediction of the depletion voltage. This is in agreement with the result predicted by the DC simulation which depicts the electric field strength as a function of depth within the a-Si:H substrate and predicts depletion for electric field strengths between 4 and 8 V µm −1 for the 12 µm thick device.
Applying the same methodology to the data set obtained for the 30 µm thick device (see Figure 9), the depletion voltage, as determined by the CCE as a function of the the square root of the bias, is determined to be 300 ± 30 V, corresponding to an electric field strength of 10 V µm −1 . This result fits with the earlier prediction of depletion occurring between 300 and 360 V based upon the electric field distribution. In comparison, the depletion voltage obtained by the formalism proposed by Despeisse shows a discrepancy of approximately 20%, corresponding to an electric field strength of 14.4 V µm −1 .

DISCUSSION AND CONCLUSION
Amorphous silicon technology presents several advantages for the development of ionizing radiation detectors: it allows deposition over flexible materials such as plastics; it can be easily patterned creating pixellated large area detectors; it has an intrinsically high radiation hardness associated with its high degree of native disorder. For these reasons, a tool capable of optimizing the design of radiation detector geometries based on a-Si:H is desirable. In this study, a new material model to mimic the electrical and charge collection behavior of a-Si:H has been developed using the Sentaurus TCAD suite. Two a-Si:H based device structures, identical with the exception of their thickness (12 and 30 µm) were modeled. In order to model the highly disordered structure of an a-Si:H substrate, multiple acceptor/donor defects characteristic of the quasi-continuous conduction/valence tail states and the Gaussian mid-gap states have been introduced.
The use of a commercially available TCAD suite, largely diffused in the radiation detection scientific community, presents some advantages to develop a model for a radiation detector: the model is portable, reproducible and expandable to include further properties or parameters in the simulation of the device. A limitation associated with this strategy is the use of models which may not be developed specifically for the material adopted. At present, TCAD doesn't have a library for a-Si:H and none of the models associated with Poole-Frenkel transport have been optimized for this specific material. Because of this limitation, the leakage current density simulated as a function of the applied potential was unable to mimic the experimental behavior without application of a fitting function featuring the square-root of the electric field dependency associated with Poole-Frenkel transport in a-Si:H. Regardless, the model has demonstrated capability of determining the operating conditions, i.e., the bias required for depletion, of a-Si:H based devices. The static DC simulations were used to show the electric field distribution as a function of position through the thickness of the 12 and 30 µm devices. The model was able to show that an operational voltage for the 12 µm thick device between 48 and 96 V corresponding to an electric field strength of 4-8 V µm −1 was required to fully deplete the device; an obvious necessity considering the context of operating as a sensor for ionizing events. Likewise, the model of the 30 µm predicts an operational voltage between 300 and 360 V, corresponding to an electric field strength of 10-12 V µm −1 to fully deplete the device. A previous study by De Greef et al. proved that a model based on the combination of SRH statistics and amphoteric behavior of the defects is effective in modeling the behavior of aSi:H [35]. Whilst they did not investigate the effect of high electric fields on the dark reverse current, they proved that at low bias, recombination is the dominant source of leakage current while at high temperatures (equivalent to a high electric field condition) the formalism of amphoteric states becomes more significant. This result confirms our finding where the current density calculated by TCAD using the SRH statistics is significant only at low bias while it is dominated at high voltages by emission associated with tail states described by the Poole-Frenkel effect. The model discussed by De Greef et al. [35] calculates with great accuracy all the effects on current density due to temperature, density and spatial distribution of the defects in a device with uni-dimensional geometry (using the software D-AMPS). The Synopsys TCAD model developed in this study does not currently provide the same accuracy of the model developed in D-AMPS due to the limitations of the numerical models already embedded in the commercial tool. However, Synopsys TCAD provides the means to tune the models to achieve an acceptable accuracy. Furthermore, it provides the additional advantage of simulating 2D or even 3D models of the device. This feature allows having a better understanding of the effects on electric field distribution and on the charge collection efficiency created by specific geometries, relative position and dimensions of the electrodes which is of paramount importance for fine tuning of a position sensitive ionizing radiation detector. Transient analysis of MIPs was performed with and without activation of the Poole-Frenkel trap model for the 12 µm thick a-Si:H device structure. The results derived from the simulation with the Poole-Frenkel trap model activated depict behavior that is at odds with what has been shown experimentally within the literature [8], and so the Poole-Frenkel trap model in its current implementation is deemed to be incompatible with thick a-Si:H substrates. Conversely, the results derived from the simulation without activation of the Poole-Frenkel trap model depict the expected charge collection behavior of a device before and after being fully depleted. Analysis of the aforementioned results gives a depletion voltage of 68 ± 20 V which is within 1% of that predicted by the formalism provide by Despeisse et al. [12] of 69 V. The depletion voltage obtained from the CCE as a function of the square root of the applied bias for the 30 µm thick substrate is 300 ± 30 V, which differs from the formalism of Despeisse by approximately 20%. The model underestimates the depletion voltage in thick devices in comparison with the results obtained by experimental methods because a thicker substrate presents larger probability to have a different degree of disorder across the intrinsic layer, requiring a larger voltage to be depleted. Such behavior can be taken into account by creating successive layers of substrate with increasing concentration of tail states. The good agreement between the simulated charge collection and the aforementioned requirements to fully deplete the device indicate the suitability of the TCAD modeling approach as a predictive tool for the design and optimization of a-Si:H radiation detectors. It should be noted that the determination of depletion voltage via the charge collection efficiency as a function of voltage is not the standard method and should only be treated as an approximation. Thus the next step in this work will be to utilize the small-signal analysis simulation method, whereby the capacitance of the a-Si:H devices as a function of bias may be determined. Future work in this area will focus upon optimization of device geometries for particle detection, including charge sharing, inter-pixel capacitance, and electric field distribution, as well as continued development of the model in terms of radiation damage upon the effects and variation of the CC and depletion as a function of the temperature. In addition, the Sentaurus Physical Model Interface (PMI) will be investigated as a means of implementing modifications to the existing Poole-Frenkel mobility model to be applicable for thick a-Si:H substrates in order to improve the robustness of the numerical model.

DATA AVAILABILITY STATEMENT
A detailed list of the energy states included in the simulations is included in the Supplementary Material.

AUTHOR CONTRIBUTIONS
JD, BT, and MP contributed equally to the conduction of all simulations. ML assisted with the data analysis of the J-V characteristics of the simulation data. All other authors contributed to discussion of the results and the write-up of the manuscript.

FUNDING
This study is part of the project 3D-SiAm, funded by the National Commission 5 of the Istituto Nazionale di Fisica Nucleare-Italy.