Monte Carlo Solution of High Electric Field Hole Transport Processes in Avalanche Amorphous Selenium.

Amorphous selenium lacks the structural long-range order present in crystalline solids. However, the stark similarity in the short-range order that exists across its allotropic forms, augmented with a shift to non-activated extended-state transport at high electric fields beyond the onset of impact ionization, allowed us to perform this theoretical study, which describes the high-field extended-state hole transport processes in amorphous selenium by modeling the band-transport lattice theory of its crystalline counterpart trigonal selenium. An in-house bulk Monte Carlo algorithm is employed to solve the semiclassical Boltzmann transport equation, providing microscopic insight to carrier trajectories and relaxation dynamics of these non-equilibrium "hot" holes in extended states. The extended-state hole-phonon interaction and the lack of long-range order in the amorphous phase is modeled as individual scattering processes, namely acoustic, polar and non-polar optical phonons, disorder and dipole scattering, and impact ionization gain, which is modeled using a power law Keldysh fit. We have used a non-parabolic approximation to the density functional theory calculated valence band density of states. To validate our transport model, we calculate and compare our time of flight mobility, impact ionization gain, ensemble energy and velocity, and high field hole energy distributions with experimental findings. We reached the conclusion that hot holes drift around in the direction perpendicular to the applied electric field and are subject to frequent acceleration/deceleration caused by the presence of high phonon, disorder, and impurity scattering. This leads to a certain determinism in the otherwise stochastic impact ionization phenomenon, as usually seen in elemental crystalline solids.

Introduction

Since the discovery of photoconductivity in amorphous selenium (a-Se) films in the 18th century,[1] devices based on a-Se have garnered many applications such as digital radiography,[2−5] solar cells,[6] threshold and memory switching,[7] and Xerography.[8] This unique disordered material is a large-area, room-temperature, direct band gap semiconductor that has shown ultralow dark current density (∼pA/cm2),[9] and has a wide detectable wavelength range (∼90% in the blue wavelength), which, in combination with scintillators, covers a significant part of the electromagnetic spectrum that includes visible,[2] ultraviolet (UV),[10] and X-ray applications.[11] Two essential features of the avalanche phenomenon are that first, in a-Se, only holes get hot and undergo impact ionization (vacuum tube a-Se devices have achieved gains as high as ∼1000),[12] and, second, the avalanche process is deterministic and non-Markovian, thus leading to a very low excess noise factor.[13,14]Figure compares the impact ionization coefficients for electrons βe and holes βh in different materials. A k-value (k = βe/βh) of ∼1 as seen in Ga0.52In0.48P[15] and GaAs[16] means that hole/electron ionization feedback will be unavoidable during electron/hole avalanche, which can have a deleterious impact on the excess noise. Although c-Si shows a higher asymmetry between βe and βh,[17] it is not enough to guarantee electron-only impact ionization, and excess noise still dominates at higher electric fields. Some wide band gap materials like 4H-SiC[18] and GaN[19] exhibit a low k value; nevertheless, their wide band gap leads to a very high breakdown voltage and their signal response is confined to the UV region of the electromagnetic spectrum. In comparison, a-Se with its hole-only impact ionization (k = βe/βh ≈ 0 for electric fields E < 1.05 MV/cm)[20] exhibits noiseless avalanche multiplication gain even at high electric fields.[2,9,13,14] Thus, understanding the dynamics of high field hot carrier transport in a-Se is of paramount importance. Recently, crystalline compound semiconductors InAs[21] and HgCdTe[22] have shown substantial promise and evidence of a true single carrier (electron only) impact ionization process with k = βh/βe ≈ 0. However, low yield, low gain (1–50), and fabrication challenges invigorated by high dark current density (∼μA/cm2) in these narrow band gap materials with limited responsivity in the infrared spectrum, and incompatibility with room-temperature operation, have hindered their potential growth to establish a true solid-state photo-multiplier with very high gain, low excess noise, high dynamic range, and linear mode of operation.

Figure 1

Change in impact ionization coefficient for electrons βe (stars) and holes βh (open circles) with inverse electric field.

Although glassy/amorphous materials are intrinsically disordered on the atomic scale, many of their electrical and semiconducting properties originate from the local structural characteristics. This observation by Ioffe and Regel states that medium- to short-range order must be maintained for amorphous solids, for a continuum of semiconducting properties to be observed across its allotropic forms.[23] Mott employed this observation and predicted the occurrence of a minimum conductivity in disordered systems considering delocalized extended-state charge transport.[24,25] The Mott criterion has been corroborated for a-Se using time-of-flight (TOF) measurements, when the activated hole mobility saturated with a value of ∼1.5 cm2 V–1 s–1,[26] beyond the onset of impact ionization avalanche once activated trap-limited transport shifts to band-like transport via delocalized extended states.

The multiscale simulation approach we have adopted in this work uses density functional theory (DFT) to calculate the electronic structure, band energies, valence band density of states (VB-DOS), Hamiltonians, Hamiltonian derivatives, dynamical matrix, and phonon dispersion in the crystalline counterpart of a-Se, that is, t-Se (as seen from molecular dynamics simulations; refer to the Supporting Information, Supplementary S1 for further details). The simulated reduced radial distribution function of a-Se compared well with both m-Se and t-Se (showed strong correlations up to 10 Å, indicating a similar short-range order). Yet, the intrinsic metastability of the monoclinic phase results in a lack of experimentally verifiable results, thus leading us to choose t-Se as the crystalline counterpart of selenium that was modeled using DFT.[27] The parameters calculated using DFT are coupled to an in-house Monte Carlo (MC) simulation framework that solves the semiclassical Boltzmann transport equation (BTE) to gain insight into the extended-state high field hole-transport process in a-Se. (refer to the Supporting Information, Supplementary S2 for further details ).

In the past, we considered an MC-BTE solution using a parabolic band approximation to the VB-DOS, to model acoustic and non-polar optical phonon-limited hole transport in t-Se.[27−29] In contrary to previous considerations,[30] we showed how holes in selenium can undergo both elastic and inelastic collisions and yet get hot, thus gaining energy at a higher rate from the electric field than they lose to the lattice in the form of phonon scattering.[27] While our model was good enough to study low electric field drift mobilities, which compared well to experimental results, in this work, we had to extend the model to examine high-field transport in a-Se by accounting for the lack of long-range order that exists in the amorphous phase.

A feature common to all amorphous chalcogenide solids is the presence of bonding defects that occurs due to atoms being over- or under-coordinated. The most studied defect in a-Se is the so-called valence alternation pair (VAP) formed via an exothermic reaction with a negative effective correlation energy.[31,32] VAP defects are represented by two selenium atoms in close proximity, found as a combination of a positively charged three-fold coordinated atom, Se3+, and a negatively charged one-fold coordinated atom, Se1–.[33] Experimental observations of thermally activated hole and electron drift mobility indicate the presence of traps in the mobility gap of a-Se, whose exact nature, though inconclusive, has been widely attributed to VAP defects (a large concentration of 1018–1020 cm–3) in the atomic structure of a-Se.[34] In the regime of extended-state transport at high electric fields in a-Se, these VAP defects will not act as trap states but instead associate themselves into pairs, thus acting as scattering centers for holes causing hole–dipole scattering. Such an approach limits the validation of our simulation results with experiments to high electric fields only, when the transport has shifted from activated localized hopping to non-activated extended-state band-like transport. The thermal agitation of the a-Se lattice causes oscillations of the VAP type dipoles, resulting in a second source of perturbation, modeled using polar optical phonon scattering. Lucovsky et al. performed infrared reflectivity measurements,[35,36] and arrived at the conclusion that the chalcogenide family, including the homopolar selenium and tellurium, has fairly strong optical phonon coupling.[37] After adding polar optical phonon scattering, disorder was introduced as an additional scattering process whose primary effect is to produce elastic and isotropic hole-lattice scattering. The strength of this scattering mechanism is inversely proportional to the magnitude of the short-range order parameter.[38] Additionally, a non-parabolic band approximation to the VB-DOS is used in the MC-BTE simulation to better stabilize the hole energy distributions and eliminate the artificial polar runaway effect. Finally, we model the hole impact ionization avalanche in a-Se as a separate inelastic but isotropic scattering process.

Computational Methods

Experimental measurements of neutron inelastic scattering, X-ray ultraviolet, and inverse photo-emission in a-Se and t-Se have shown an almost identical phonon and electronic density of states.[39−41] Moreover in a-Se, the extended-state wave-function involved in the impact ionization process has a comparable dispersion and phase-coherence to that of t-Se and thus allows this quantum mechanical MC-BTE high-field modeling using the wave-vector k (k is not considered a good quantum number in disordered materials due to the lack of periodic potentials).[27]

To improve accuracy of the MC-BTE model, an analytical approach, using the k·p method, is implemented to obtain the non-parabolic equation:[42]E(1 + αE) = ℏ2k2/2mc, where mc is the conductivity mass and α is the non-parabolicity factor from the k·p method, which depends on the material as α = 1/Eg(1 – mc/mo)2, where Eg is the band gap for a-Se (2.1 eV) and mo is the hole rest mass. The value of mc has been reported in t-Se to be 0.29mo in the direction parallel to the c-axis (m∥) and 0.75mo in the direction perpendicular to the c-axis (m⊥).[43] Given that we are ultimately interested to study the extended-state hot hole transport in the amorphous phase with topological disorder, we treat all solids as an isotropic continuum without any directional dependence. Acoustoelectric current saturation observations of hole drift mobility[44] backed up by our previous MC-BTE simulations in t-Se[27] show that the mobility ratio between the ⊥ and ∥ directions to the c-axis in t-Se is μ∥/μ⊥ ≈ 4. The higher scattering rates and lower mobility of holes (closer to that of a-Se: μ ≈ 1.5 cm2 V–1 s–1) along the ⊥ direction to the c-axis in t-Se lead us to the assumption that transport and hole-phonon coupling parameters in the ⊥ direction to the c-axis in t-Se should be used in the scattering rate calculations for a-Se (Table ).

Table 1

Hole-Phonon Coupling Parameters Used in our Calculation of Scattering Rates in a-Se

mechanism	type	parameter	value	exp conditions/method
acoustic phonons	elastic/isotropic	acoustic deformation potential Ξac (eV)	6 (t-Se⊥c)	computational DFT[27]
		sound velocity vs (m s–1)	2150 (t-Se⊥c)	comp. DFT slope of acoustic modes of vibration[27]
		DOS mass md	1.4mo (t-Se⊥c)	estimated from thermoelectric power with an isotropic single valence band maximum[46]
		density (kg/m3)	4819 (t-Se⊥c)	T = 298 K[47] calc. From X-ray data
non-polar optical phonons	inelastic/isotropic	optical phonon energy ℏωo (meV)	28.9 (t-Se⊥c)	computational DFT[27]
		optical deformation potential D0 (eV/Å)	3 (t-Se⊥c)	computational DFT[27]
polar optical phonons	inelastic/anisotropic	low ε0 frequency dielectric constant	7.35 (t-Se⊥c)	oscillator fit-IR data[48]
		high ε∞ frequency dielectric constant	6.66 (t-Se⊥c)	oscillator fit-IR data[48]
		non-parabolic factor α (eV–1)	0.15	analytical approach to the k·p method[42]
disorder	elastic/isotropic	short-range order	∼10 Å (a-Se)	molecular dynamics sim.;[49] constant matrix element[38]
VAP dipole	elastic/anisotropic	density of scattering dipole pairs NT (cm–3)	8 × 1019 (a-Se)	density of VAP defects[34]
		dipole radius ao (Å)	17.32 (a-Se)	order of nearest neighbor distance[50]
		dielectric constant ε	7	average of low ε0 and high ε∞ in t-Se⊥c
		Debye length LD (Å)	6.6	calculated analytically

The effective conductivity mass in the amorphous phase was calculated using the Herring–Vogt transformation, which resulted in an α of ∼0.15.[45] In this work, the density of states mass md (1.4mo)[46] is used for the calculations of the scattering rates, and mc (m⊥ = 0.75mo) for the equation of motion and the E–k dispersion relation. Figure a compares the DFT-calculated VB-DOS with its non-parabolic band approximation. The use of α ≈ 0.15 that arises from the simplistic assumption of a single isotropic band prevents the hole energy and drift velocity distribution from running away to higher values at electric fields substantially lower than those measured experimentally.

Figure 2

(a) DFT calculations of the VB-DOS are shown by the solid blue line. The dotted red line represents the non-parabolic band approximation (α = 0.15) to the VB-DOS. (b) Scattering rates of the mechanisms relevant for a-Se. Elastic and isotropic mechanisms (acoustic and disorder) correspond to the dashed lines. Inelastic and anisotropic polar optical phonon scattering is denoted by dashed dotted lines. The dotted lines denote isotropic but inelastic scattering caused by non-polar optical phonon vibrations. Impurity scattering from VAP dipoles is denoted by a solid line. (c) The probability of the scattering angle is shown for the three anisotropic scattering mechanisms in the simulation. Scattering from polar optical phonons and the VAP type dipoles become more anisotropic at high electric fields when the hole energy increases, thus favoring small angle forward scattering. (d) Real space trajectories of seven holes comparing lateral spread at low (100 kV/cm) and high (1000 kV/cm) electric fields. (e) The impact ionization Keldysh fit used. The distribution of the hole impact ionization band shows a narrow normal distribution (full width at half maximum ≈ 0.45 eV), an indication toward the deterministic nature of the avalanche process in a-Se.

In the regime of extended-state transport at high electric fields in a-Se, we assume that VAP defects act as a single source of perturbation causing hole–dipole and polar optical phonon scattering. For the first time, we have formulated the hole–dipole scattering interaction for non-parabolic bands (Supporting Information, Supplementary S3 and Supplementary S4). To account for the lack of long-range order, we assume that the disorder existing in the amorphous phase causes a hole–lattice elastic/isotropic interaction, which is only governed by the VB-DOS and a constant matrix element (Supporting Information, Supplementary S3).[38]

Figure b shows the scattering rates obtained from the non-parabolic band approximation discussed above. Collisional broadening effects come into play at scattering rates of ∼1015 s–1, when the average time between two successive collisions become too short.[51] However, these higher order effects are not incorporated at present and will be a topic of future investigation. Figure c shows the anisotropic nature of polar optical phonon scattering and hole–dipole scattering, thus favoring forward angle scattering, which further increases with an increase in hole energy. It is worthwhile to notice that at small energies/low electric fields, the VAP dipoles increases the lateral spread of holes.

This idea is illustrated in Figure d where a comparison between low and high electric field trajectories shows that while the actual hole path lengths are substantially increased, the VAP defects lead to more lateral spread at lower electric fields when the hole energies are comparable to the thermal energy at room temperature (∼38 meV). The longer effective path length of holes leads to an enhanced probability of energy loss via polar and non-polar phonons, which helps in stabilizing the ensemble energy distributions. Furthermore, the probability of a sharp rise in carrier energies decreases and a ″delay time″ arises before a hole–hole impact ionization event occurs, which is expected to reduce (and potentially eliminate) the excess noise in a-Se.

The only hole–hole interaction that has been accounted for in this work is impact ionization avalanche, which has been modeled using a single term power law fit, shown in Figure e (Supporting Information, Supplementary S3). We have neglected Coulomb scattering, which occurs due to charged impurities and defects. Although the VB-DOS for t-Se indicates the existence of a second band at energies greater that 10 eV, as investigated by the empirical pseudo-potential method (EPM),[52] when the VB-DOS increases abruptly, high energy processes such as intervalley scattering and interband impact ionization have been neglected. While such a process could significantly alter the high-energy tail of the hot hole distribution, in this work, we shall restrict our attention to the gross features of high electric field transport in a-Se when carrier transport shifts from hopping conduction via localized states to band-like transport via extended states.[27]

Results and Discussion

The bulk MC-BTE is used to follow the real-time scattering processes and calculate average transport characteristics on a picosecond time scale. Figure a shows the time-averaged ensemble drift-velocity calculated as a function of electric field. The saturated drift velocity is an important limiting parameter for the semiconductor industry. Experimental studies of the hot hole drift velocity in a 0.4 μm-thick a-Se HARP (high-gain avalanche rushing amorphous photoconductor) target at electric field strengths of 1000 kV/cm was measured to be 1.87 × 104 m/s.[55] The drift velocity calculated using the non-parabolic model shows a gentle peak of 1.83 × 104 m s–1 at 1000 kV/cm and saturates thereafter. This velocity saturation effect could not be observed with a parabolic bulk MC-BTE model where the drift velocity increases monotonically as a function of electric field. Figure b shows the average hole energy in a-Se as a function of electric field. The hot carrier energy in the 0.4 μm-thick a-Se HARP increases linearly as the electric field increases (as shown by solid marks in Figure b). The increase in the simulated hot hole energy in the non-parabolic MC-BTE simulation compares well to these experimentally measured values.[53] The ensemble energy of holes calculated under the parabolic approximation rapidly spreads to higher energies at electric fields much lower than the threshold for impact ionization (1000 kV/cm), a result due to polar runaway.[56] The parabolic band approximation to the DFT-calculated valence-band DOS overestimates the steady-state drift velocity and energy values at fields beyond 600 kV/cm.

Figure 3

Steady state simulation results of ensemble time-averaged (a) drift velocity and (b) average energy of holes as a function of electric fields ranging from 1–1200 kV/cm. The hollow square markers show velocity and energy spread to higher values under the parabolic band approximation. The hollow circles show average drift velocity and average energy simulated under the non-parabolic band approximation. The solid line joining solid marks shows experimentally reported values of hot hole energy in a-Se.[53] (c) TOF non-parabolic MC-BTE calculated mobility (hollow circles) compared with experimental measured saturated and electric field-independent mobility (solid markers) in a-Se.[26] The error bars indicate the statistical errors on the TOF mobility simulated for 0.5–35 μm thick a-Se bulk device lengths. (d) Impact ionization gain calculated using non-parabolic MC-BTE and compared with experimentally measured gain for 0.5–35 μm-thick a-Se films.[20,54] (e) Theoretically modeled hole impact ionization coefficient compared with measured values from Tsuji et al.(20) (f) A simulation of hole energy distribution in a 1 μm-thick a-Se film calculated using the non-parabolic MC-BTE model at different electric fields.

Figure c compares TOF-measured saturated mobility in a-Se[26] where transport has shifted from activated (trap-limited) transport at lower fields to extended-state (band-like) transport. The activated experimentally measured mobility for a-Se increases with increase in electric field and finally saturates at very high electric fields, around 900 kV/cm. The saturated extended-state experimental TOF mobility has a value of 1.5 cm2 V–1 s–1 and matches the non-activated MC-BTE calculated mobility. The MC-BTE TOF mobility is calculated as μ = L/(E × ttr), where L is the the length of the device, E is the applied electric field and ttr is the transit time of holes to cross the length of the device. The TOF mobility was calculated across device lengths of 0.5–35 μm. Although the length dependence of the TOF mobility was minute, thicker selenium simulations, in general, showed smaller mobilities as compared to thinner samples. However, there were outliers, which lead us to believe that deviation in mobility as a function of device length originates from the statistical nature of the Monte Carlo simulations, denoted by the error bars in Figure c.

Our MC-BTE calculated impact ionization gain in Figure d compares well to experimental results for 0.5–35 μm-thick a-Se samples across a wide range of electric fields. The MC-BTE simulation results, as shown in Figure d, correctly predict the decreasing onset of the avalanche electric field with increasing device thickness. Figure e compares MC-BTE-calculated hole impact ionization coefficients with experimentally measured values (refer to the Supporting Information, Supplementary S6 for the calculation of the hole impact ionization coefficient).

The Ergodic theorem states that in a statistical dynamic system, the time and space average shows similar behavior. Figure f shows the hole energy distribution for various electric fields. The histories of 1000 independent holes were simulated over a fixed distance, sufficiently far from the starting position so as to guarantee the achievement of a steady-state regime (Supporting Information, Supplementary S7). The energy of the carriers was recorded at the end of the travel. Holes starting from thermal energy (38 meV), on average, gain significant energy in the first 10 nm. At high electric fields, corresponding to high scattering rates, a short ″heat-up″ distance[51] (10–20 nm) is sufficient to achieve a steady-state distribution. At lower electric fields, this heat-up distance increases to 50–100 nm. In Figure f the carrier energies were frozen after they crossed a distance of 1 μm. At high electric fields (1200 kV/cm), the energy tail (hole energies of 2.1 eV and beyond; the impact ionization band) existing in lower field distributions undergoes the impact ionization process. At 1200 kV/cm, the total number of carriers at the end of the simulation was 2979 and 6061 (avalanche gains of 2.98 and 3.03) for 1000 and 2000 initial carriers, respectively. The MC-BTE gain is calculated by dividing the final number of carriers collected at the end of travel by the initial number of carriers.

Experimental results have previously shown that the excess noise in a-Se ≈ 1 because of the non-Markov branching of hot holes.[13] This might be a result of the high scattering rates existing in the disordered phase of selenium, which leads to a non-zero “dead space” distance, defined as the minimum distance of travel in the direction of electric field before a carrier can attain the ionization threshold energy required for avalanche.[21,57] This, in turn, can average out the noise arising due to the stochastic avalanche process and increase determinism, thus leading to a very low excess noise factor in a-Se.[13]

Conclusions

A modeling scheme was presented and discussed, which is used to study charge carrier transport in disordered structures, which lack the long-range order present in their crystalline counterparts. In other words, we have studied the transport of holes in a-Se by using an MC-BTE technique in which the carrier free-flights are interrupted by scattering from acoustic, polar and non-polar optical phonons, disorder, and dipole scattering. Material parameters for selenium’s crystalline allotropic form t-Se were calculated using DFT, and the approximation of a single isotropic non-parabolic band with non-parabolicity factor α = 0.15 was used in the MC-BTE. The saturated drift velocities and ensemble hole energies obtained with our model matched closely with experimentally measured values in a 0.4 μm-thick a-Se HARP target.[55] Calculations of our MC-BTE TOF mobility matched exactly to experimentally measured non-activated extended state mobility.[26] Moreover, the calculated impact ionization gain matched closely to experimentally measured values in 0.5 to 35 μm a-Se HARP samples.[20,54] Our MC-BTE model correctly predicted the decreasing onset of the avalanche electric field with increasing device thickness, a characteristic of the avalanche phenomenon in a-Se. We speculate that the deterministic nature of impact ionization avalanche in a-Se occurs due to the existence of a non-zero dead space and the non-ballistic nature of hole transport in a-Se, fostered by frequent scattering at high electric fields. Consequently, there is an increased spread of holes in the direction lateral to the applied electric field, seasoned by frequent acceleration and deceleration caused by the cumulative affect of the electric field and lattice/impurity/disorder scattering, which results in averaging of the distance traveled by the hot holes over a finite delay time before an impact ionization event occurs. Moving ahead, this multiscale approach of combining molecular dynamics, DFT, and MC-BTE will be used in calculating the excess noise and spatial resolution in a-Se-based devices.