Invited Review Article: Physics and Monte Carlo Techniques as Relevant to Cryogenic, Phonon and Ionization Readout of CDMS RadiationDetectors
Abstract
This review discusses detector physics and Monte Carlo techniques for cryogenic, radiation detectors that utilize combined phonon and ionization readout. A general review of cryogenic phonon and charge transport is provided along with specific details of the Cryogenic Dark Matter Search detector instrumentation. In particular this review covers quasidiffusive phonon transport, which includes phonon focusing, anharmonic decay and isotope scattering. The interaction of phonons in the detector surface is discussed along with the downconversion of phonons in superconducting films. The charge transport physics include a mass tensor which results from the crystal band structure and is modeled with a Herring Vogt transformation. Charge scattering processes involve the creation of NeganovLuke phonons. Transitionedgesensor (TES) simulations include a full electric circuit description and all thermal processes including Joule heating, cooling to the substrate and thermal diffusion within the TES, the latter of which is necessary to model normalsuperconducting phase separation. Relevant numerical constants are provided for these physical processes in germanium, silicon, aluminum and tungsten. Random number sampling methods including inverse cumulative distribution function (CDF) and rejection techniques are reviewed. To improve the efficiency of charge transport modeling, an additional second order inverse CDF method is developed here along with an efficient barycentric coordinate sampling method of electric fields. Results are provided in a manner that is convenient for use in Monte Carlo and references are provided for validation of these models.
Contents
 I Introduction
 II Radiation Modeling
 III Phonon Simulation
 IV Quasiparticle Down Conversion

V Charge Monte Carlo
 V.1 Introduction
 V.2 Holes: propagation and scattering with isotropic bands and isotropic phonon velocity
 V.3 Electrons: propagation and scattering with anisotropic bands and isotropic phonon velocity
 V.4 Charge Time Steps, First Order
 V.5 Charge Time Steps, Second Order
 V.6 Select constants for charge Monte Carlo
 VI ElectricField Lookup
 VII Transition Edge Sensor Simulations
 VIII Final Remarks
I Introduction
Cryogenic radiationdetectors that utilize ionization, phonon and / or scintillation measurements are being used in a number of experiments. Both the Cryogenic Dark Matter Search (CDMS) Ahmed2010 (); Ahmed2011 () and EDELWEISS Armengaud2010 () dark matter search utilize silicon and / or germanium targets to detect recoils of radiation in the target masses. A combination of ionization and phonon readout is used to provide discrimination of gamma and neutronrecoil types. The CRESST dark matter search utilizes CaWO targets and readout scintillation and phonon signal to discriminate between recoil types. The advantage of reading out both phonon and ionization (or scintillation) signals comes about from the differing ratios of ionization and phonon energy or scintillation and phonon energy created in electron and nuclearrecoils in the detectors. The ratio of these two energies leads to a powerful discriminator for the experiment’s desired recoil type.
Both the ionization and phonon readout can be used to generate position estimators for the initial radiation interaction, leading to fiducial volume selection. In the ionization signal this is generally accomplished by instrumenting different parts of the detector with independent readout channels and vetoing events with large signal contribution outside of the desired fiducial volume. In the phonon signal it is generally required to measure the early, athermal component of the phonon signal which still retains a position dependent component.
The physics required to accurately model these detectors is presented in this paper along with appropriate numerical tricks that are useful for an efficient detector Monte Carlo. This paper proceeds with a review of radiation interactions, charge transport physics, phonon transport physics, instrumentation. Monte Carlo techniques and relevant physical constants are included where appropriate.
This paper will focus on the use of silicon and germanium detector masses, both of which are group IV semiconductors. However there are other relevant materials in use such as calcium tungstate (CaWO) which leads to a small loss of generality.
i.1 The CDMS Experiment and Detectors
The Cryogenic Dark Matter Search Ahmed2010 (); Ahmed2011 () utilizes silicon and germanium detectors to search for Weakly Interacting Massive Particle (WIMP) dark matter Spergel2007 (); Tegmark2004 () candidates. The silicon or germanium nuclei provide a target mass for WIMPnucleon interactions. Simultaneous measurement of both phonon energy and ionization energy provide a powerful discriminator between electronrecoil interactions and nuclearrecoil interactions. Background radiation primarily interacts through electronrecoils whereas a WIMP signal would interact through nuclearrecoils. The experiment is located in the Soudan Mine, MN, U.S.A.
The most recent phase of the CDMS experiment has involved fabrication, testing and commissioning of large, 3 inch diameter, 1 inch thick [100] germanium crystals. The CDMSiZIP (interleaved Z–dependent Ionization and Phonon) detectors are 3 inches in diameter and 1 inch thick with a total mass of about 607 grams Brink2006 (). The iZIP detector utilizes both anode and cathode lines on the same side of the detector similar to a MicroStrip Gas Chamber (MSGC) Knoll (); Oed1988 (); Luke1994 (); Brink2006 () as shown in Figure 2 and 2. Unlike an MSGC however, there is a set of anode and cathode lines on both sides of the detector. This ionization channel design is used to veto events interacting near the detector surfaces. An amorphous silicon layer, deposited under the metal layers, increases the breakdown voltage of the detectors. The total iZIP aluminum coverage is 4.8% active and 1.5% passive per side.
Ii Radiation Modeling
When using a Monte Carlo of a detector, it is often helpful or necessary to have a numerical model of radiation interactions in the detector. Many readers will find it valuable to use separate modeling software such as GEANT4 Agostinelli2003 (). A brief description of these interactions follows.
Low energy gammarays (xrays) predominantly interact via photoelectric absorption in which all of the gammaray energy is deposited in a single interaction location. High energy gammarays interact via Compton scattering in which some of the gammaray’s initial energy is transferred to an electron and the gammaray continues along a different trajectory with reduced energy. The gammaray will generally continue to scatter until it leaves the detector volume or terminates with a photoelectric absorption. In silicon (germanium), for photon energies greater than 60 (160) keV, Compton scattering dominates Knoll (); RPP ().
Both of these electron interactions result in a high energy electron being produced which then undergoes a rapid cascade process resulting in a large number of electronhole pairs RPP (); Cabrera1993 (). This initial cascade process ceases around the scale of the mean electronhole pair creation energy () resulting in an expected number of electronhole pair . Due to correlations in the cascade process, the variance in the number of electronhole pairs is reduced, relative to Poisson statistics, and given by , where is the Fano factor Fano1947 (). These high energy electronhole pairs will then shed phonons until they reach the semiconductor gap energy which results in the fraction of energy in prompt phonons and the remainder in the electronhole pair system.
Neutrons that interact in the detector bulk will knock an ion out of its lattice site and displace it to some other location in the crystal. This high energy ion will interact with both the lattice ions and / or valence electrons, with competing cross sections, before reaching some other location in the crystal. Interactions with the lattice ions can be described to first approximation as Rutherford scattering Rutherford1911 () with differential energy loss per unit length , where is the ion’s velocity Bonderup1978 (). For interactions with valence electrons, the number of electron states which can be excited to an accessible state outside of the Fermi sphere scales like velocity hence the differential energy loss per unit length scales like Bonderup1978 (). A description of this process is shown in Figure 3. There are important screening and velocity dependent cross sections in both energy loss lengths Bonderup1978 (); Lindhard1954 (); Lindhard1963 (); Lindhard1963_2 (); Jones1975 () but to first approximation these equations show the velocity dependent competition in the two scattering rates. Whichever interaction occurs first, there is generally a cascade of many scattering events resulting in a larger amount of energy being deposited in the ion lattice compared to gammaray interactions. For historical reasons, the reduced amount of energy in the electronhole (ionization) system, compared to an equal energy deposition by a gammaray, has brought about the term nuclear quenching to describe the reduced ionization signal. The details of the cascade and reduction in the number of electronhole pairs is described by Lindhard theory Lindhard1963 () and given as where , and Lewin1996 (). The recoil energy is and given in units of keV, is the atomic number and is the atomic mass.
Beta radiation represents another class of electronrecoil interactions. The attenuation lengths are much shorter, however, resulting in a class of events sometimes referred to as surface interactions. These surface interactions can result in signals that differ from bulk events of the same energy and electron / nuclearrecoil type. For example, the events may be located in regions of large electric fringing fields, directing charges away from readout electrodes or near mounting hardware that absorbs scintillation light.
Iii Phonon Simulation
iii.1 Introduction
Phonons, in the context of this review, are quantized vibration modes that exist in periodic structures such as silicon and germanium crystals. They are excitations which, in the lattice, and in materials sufficiently cooled that charge carriers are frozen out, mediate thermal transport. They are described by the Hamiltonian
(1) 
where is the mass of the lattice atoms and is the frequency of oscillations about the center of the harmonic potential between an atom and its nearest neighbor atom Kittel (); Ashcroft (); Wolfe ().
In one dimension, the solution to the timeindependent Schrodinger equation yields the different oscillation modes at atom
(2) 
where here, the wave number , is the lattice spacing and are the number of lattice sites.
In Monte Carlo, phonons are simply treated as non interacting particles with decay and mass defect scattering properties described in the remainder of this chapter. In general they have a nonlinear dispersion relationship; however, due to rapid down conversion to lower frequencies they spend most of their time in an energy region with linear dispersion relationship. Hence, a linear dispersion relation is sufficient for phonon transport modeling in these detectors.
iii.2 Prompt Phonon Distributions
Immediately after the recoil, a population of prompt phonons exist in the immediate vicinity of the interaction point. Particular details of the frequency distribution and mode population are not well known but we can make a few deductions, explained in more detail in the following sections, that will lead us an initial distribution for Monte Carlo.
Anharmonic decay, due to nonlinear terms in the elastic coupling between adjacent lattice ions, causes the phonons to rapidly down convert into a lower frequency distribution. This process allows us to start a Monte Carlo with any high frequency distribution and details of the distribution will rapidly be lost; we use the Debye frequency as a naive starting point. Isotope scattering, which also occurs at a high rate for a high frequency phonons, causes the phonons to obtain their equilibrium mode density; we use the equilibrium mode density as a naive starting point. The approximations are good in the sense that the detector’s phonon response is insensitive to variations in the distributions.
They are valid since the detectors are large compared to the initial characteristic interaction lengths of phonons. Furthermore, later generations of phonons are ballistic and timing information in the measured phonon pulses is determined by the detector geometry and the loss rate of phonons at surfaces.
iii.3 Phase Velocities and Polarization Vectors
The so called phase velocity surfaces represent the direction dependent phonon phase velocity . In general they are given by the eigenvalue Equation 3.
(3) 
where is the crystal’s mass density,
is the phonon frequency,
is a component of the polarization vector
,
are components of the elastic constant tensor and,
is a component of the phase velocity vector
Kittel ().
Not all of the elastic constants are independent, reducing via a Voigt contraction Nye () the number that we need to keep track off. Additionally, symmetries in a cubic crystal allow for further reduction in components and we can define three independent constants , and . This contraction simplifies Equation 3 significantly Ashcroft () and we are left solving matrix III.3 for its eigenvectors and eigenvalues.
(4) 
The eigenvectors represent the three polarization vector directions and the three eigenvalues equal for the longitudinal, slowtransverse and fasttransverse modes. The anisotropy in the cubic silicon crystal leads to the phase velocity surfaces being nonspherical (see Figure 4). The phase velocities are used to determine both the group velocities (Section III.4) and also the isotope scattering rates.
iii.4 Group Velocities
Phonon group velocities are found by solving
(5) 
The slight lack of sphericity in the phase velocity surfaces (see Figure 4) has a very dramatic effect on the transverse phonon group velocities Northrop1980 (); Tamura1991 (); Maris1993 (); Tamura1993LT () (see Figure 5). The longitudinal phonon’s group velocity is only mildly affected. Energy is focused in the direction of heavy banding and leads to the term phonon focusing. The point density in the plots is misleading as the three modes are shown to be equally populated. Isotope scattering including anisotropic scattering rates (see Section III.5) leads to the phonon modes in silicon being populated as follows: slowtransverse (55%), fasttransverse (35%) and longitudinal (10%).
iii.5 Anisotropic Isotope Scattering
Phonons scatter off mass defects in the crystal (see Figure 6). Additionally, they can change modes. The bulk scattering rate is given by
(6) 
where is the phonon frequency and is a scattering rate constant Tamura1991 (); Tamura1993LT (); Maris1990 (); Msall1993 (); Tamura1993 () (see Table 1). The scattering rate for individual phonons is given by
(7) 
where is the final state phonon frequency in Hz, is the polarization vector, represents the initial phonon and represents the outgoing phonon Tamura1993LT (); Tamura1993 (). It is the dot product in Equation 7 which allows mode mixing and the denominator which ensures the correct populations in the ratios. In silicon the populations when including anisotropic scattering rates are slowtransverse (55%), fasttransverse (35%) and longitudinal (10%). The standard treatment is to determine if a phonon isotope scatters via and then determine its polarization and direction via . After the initial anharmonic decay has settled down, isotope scattering dominates.
This process is unfortunately computationally expensive due to samplerejection techniques NumericalRecipes () and after several iterations an isotropic scattering process can be used for individual phonons with little loss of modeling accuracy.
iii.6 Anharmonic Decay
iii.6.1 General Case
Nonlinear terms in the elastic coupling constants cause a longitudinal phonon to down convert to two lower energy phonons (see Figure 7). The bulk decay rate is given by
(8) 
where is the phonon frequency and is a decay rate constant Tamura1993LT (); Maris1990 (); Msall1993 (); Tamura1993 () (see Table 1). The decay rate for transverse phonons is negligible Tamura1985TA (). The three body problem requires that both energy and momentum are conserved. These conditions make an exact solution computationally prohibitive for large amounts of phonons.
iii.6.2 Isotropic Approximation
To allow computations to proceed in a finite time, an exact solution to anharmonic decay is abandoned and an isotropic approximation is used (the full anisotropic phase velocities and group velocities are still easily used for isotope scattering and phonon transport). The energy distribution calculations are still difficult but fortunately have already been carried out Tamura1985 (); CabreraTamura (). Once the energies have been determined, calculating the resultant scattering angles based on energy and momentum conservation is fairly straightforward. Due to the different energymomentum dispersion relations for the longitudinal and transverse phonons, there are two different decay branches, and . The energy distribution is given by
(9) 
where ,
and,
.
This approximation results in the outgoing phonons having an angular displacement from the initial phonon given by
(10) 
(11) 
The energy distribution is given
(12) 
where ,
,
,
,
,
and,
.
The constants and are the Lamé constants and and are third order elastic constants in an isotropic model (there is additionally a third independent elastic constant but it drops out of the equations).
This approximation results in the outgoing phonons having an angular displacement from the initial phonon given by Equations 13 and 14.
(13) 
(14) 
where the trivial substitution is made for the second phonon. With these closedform energy densities and scattering angles, plots can be generated to aid understanding of these events (see Figures 9 and 9).
These angles are relative the initial momentum vector and need to be converted into the Monte Carlo coordinate system. Polar coordinates are useful and angles are provided in this system. In addition to rotation angles and (or and ) an additional azimuth angle, relative to the initial momentum vector , and in the isotropic approximation is randomly distributed from , is specified as .
In terms of initial elevation and azimuth angles and , scattering angles and and azimuth scattering angle , the final angles and that describe the phonon momentum vectors are
(15) 
and
(16) 
The final angles for the other phonon momentum vector are found by replacing with and with .
iii.7 Phonon Losses at Surfaces
Eventually the phonons will interact at surfaces where they are instrumented, reflected back into the crystal, down converted to a lower energy or are lost to the environment.
Phonons can reflect either specularly or diffusively from the surfaces. Specular reflection can occur on smooth, untreated semiconductor wafer surfaces. It is the simplest to describe with incident and reflected angles relative to the normal equal, .
Diffusive scattering is also common on surfaces that have been damaged or roughened during fabrication. In the ideal case the scattering angle is described by Lambert’s cosine law where the angular distribution scales like where the angle is measured relative to the normal. Scattering surfaces that satisfy Lambert’s cosine law scatter phonons isotropically regardless of their incident angle. Generally diffusive scattering has been found to be a good model for phononsurface reflections, likely due to some small roughness in the surface.
Phonons are strictlyspeaking eigenstates of a Hamiltonian that describes an infinitely large, periodic lattice. This description necessarily breaks down at the detector surfaces and there is some probability that phonons down convert to lower energy daughters at the surface. The details of this process would be highly material / surface treatment dependent but the probability of this occurring for a particular phononsurface interaction will be small for a high purity crystal. It is generally easiest to tune this probability by running numerous Monte Carlo to find the best value. In the CDMS detectors we have measured a loss of 0.1% for each phononsurface interaction Leman2011 ().
It is the goal of the experimental setup to absorb the phonons into some sensor and provide instrumentation into a data acquisition system. The details for the phononsensor interaction probability are again complicated and highly depend on the type of absorber and attachment / fabrication details. Acoustic mismatch theory provides a good starting point for analytic calculations and can be performed over both normal and nonnormal incidence angles. The relevance of these calculations can be lost however when detector to detector variations are considered. Additionally the angular dependence of such calculations can be washed out when integrating over a distribution of phonon incident angles and phonon energies. In practice it is usually again easiest to tune this probability by running numerous Monte Carlo and identifying the best fit value. In the CDMS detector there is additionally an amorphous silicon dielectric in between the crystal and thin metal films; we have empirically matched a phononaSialuminum interaction (details of the interaction are discussed in Section IV) probability of 33% with the remaining 67% diffusively scattering back into the crystal Leman2011 (); Leman2011_4 (). For any well designed detector, phonon absorption into the instrumentation sensors will dominate over other loss processes allowing the probability to be tuned by matching pulse decay times.
iii.8 Time Steps
It would be grossly inefficient to run all phonons with the same time step. Therefore we generate scattering times according to the distributions in Sections III.5 and III.6. The scattering and decay probabilities go like where is a combination of isotope scattering and anharmonic decay rates as given by Matthiessen’s rule .
This scattering time has to be compared with the time that the phonon will take to interact with a surface and the event that occurs soonest will be chosen for each individual phonon. If it is determined that the bulk interaction time is less than the surface interaction time, then it must be determined which event occurs based on their relative, frequency dependent rates. This can be done by drawing a uniform random number and comparing the rate of the process in question to the total rate. For example, if then an anharmonic decay is selected.
iii.9 Random Number Sampling
Only in rare circumstances will a uniform random number be needed in Monte Carlo without some transformation. Often we are trying to draw the number out of the probability distribution function (PDF) . An efficient method for transforming to the desired probability distribution function is to integrate to find the cumulative distribution function . The cumulative distribution function (CDF) has the desirable property that it is bounded by as is . The CDF is then inverted and solved at to determine , NumericalRecipes ().
As an example, we first considering the bulk interaction rate where the probability of having an interaction is . After integration and inversion, the randomly generated scattering time is given by
(17) 
As a second example we can consider diffuse scattering off the detector walls. In a spherical coordinate system where represents the angle from the surface normal, the PDF in this coordinate is . This is also easily integrated and inverted to yield . Care must be taken in this example however since is defined over the domain [1,1] which modifies to .
There are times when a PDF cannot be analytically integrated to yield a CDF or the CDF cannot be inverted. This is an unfortunate situation since an expensive rejection technique is required. This technique involves drawing a pair of uniform random numbers where and . If then is retained, otherwise is rejected and the process is repeated. The inefficiency of this method is related to the area coverage and will be successfully drawn in one of attempts.
The rejection method can be improved however for static distributions that do not change during the Monte Carlo run. An example includes diffusive scattering off of the side walls of the detector when using a spherical coordinate system. In this case, the Jacobian results in the PDF for which cannot be found analytically. The PDF can be integrated numerically to generate a CDF which is subsequently inverted. The process lacks a certain degree of elegance but is significantly more efficient than using a rejection method.
iii.10 Numerical Constants for Phonon Simulations
Table 1 lists numerous constants that are used in the phonon simulations. They define the propagation dynamics, scattering and decay rates and energy carrier statistics.
Parameter Description  Symbol  Units  Silicon  Germanium  Reference 

Isotope Scatter Rate  ()  [s]  Tamura1985 ()  
Anharmonic Decay Rate  ()  [s]  Tamura1985 ()  
Longitudinal  ()  [m/s]  9000  5310  Tamura1985 () 
Velocity  
Transverse Velocity  ()  [m/s]  5400  3250  Tamura1985 () 
Decay Rate Constant  ()  0.429  0.732  Tamura1985 ()  
Decay Rate Constant  ()  0.945  0.708  Tamura1985 ()  
Decay Rate Constant  ()  0.524  0.376  Tamura1985 ()  
Decay Rate Constant  ()  0.680  0.561  Tamura1985 ()  
Density  ()  [g/cm]  2.33  5.32  Tamura1985 () 
Decaying Ratio  0.204  0.260  Tamura1985 ()  
Lattice Constant  ()  1.66  1.29  Ashcroft ()  
Lattice Constant  ()  0.64  0.48  Ashcroft ()  
Lattice Constant  ()  0.80  0.67  Ashcroft ()  
Debye Frequency  ()  [1/s]  Ashcroft ()  
ElectronHole  [eV]  1.17  0.75  Ashcroft ()  
Energy  
Mean ElectronHole Creation Energy  [eV]  3.81  2.96  Pehl1968 () 
Iv Quasiparticle Down Conversion
Phonons have some probability of entering and interacting in the thin aluminum films that are patterned on the surface. These aluminum films make up both the phonon collecting films and ionization ground lines; it is interactions in the former that are measured in the phonon sensors. If the phonons have energy greater than or equal than twice the superconducting gap then Cooper pairs can be broken creating two quasiparticles. These quasiparticles will be of high kinetic energy and contain some probability of scattering off phonons. These daughter phonons thereby introducing a population of down converted phonons back into the metal film. These phonons could break additional Cooper pairs in a cascade process that ceases when all phonons have energy below or have a probability of being reintroduced back into the crystal. The key points in the cascade process are summarized in the following list Kaplan1976 (); Kurakado1982 (); BrinkThesis ().

Quasiparticle recombination lifetimes are long compared to quasiparticle decay and quasiparticle absorption into the aluminum films and therefore the recombination processes can be ignored.

Quasiparticle decay via absorption of a phonon is suppressed at low temperatures due to a phonon density of states term , where is the phonon energy, in the Green’s function and therefore can be ignored.

Quasiparticle decay via emission of a phonon results in phonons with an energy distribution given by , where is the quasiparticle energy and the quasiparticle density of states at temperature =0 goes like and .

Phonons break Cooper pairs producing quasiparticles with an energy distribution given by , where . The phonon is completely absorbed so that the second quasiparticle has energy .

Phonons are lost to the crystal if they reach the aluminum / crystal interface.
iv.1 Monte Carlo process ordering
In the physical processes can be ordered as follows

If the phonon energy is sufficient , a quasiparticle pair is created with the distribution previously described.

If there is a quasiparticle with energy then the quasiparticle emits a phonon with energy distribution . Quasiparticles with energy would shed a phonon with and therefore provide an endpoint for quasiparticle generation. They may be removed from the Monte Carlo.

The probability of a phonon escaping the crystal is a function of the distance to reach the aluminum / crystal interface and the phonon / Cooper pair interaction length. Given the large number of phonons it is generally not necessary to track this process in detail and instead a simple model is sufficient. On average, phonons are assumed to populate the center of the aluminum film (), where is the aluminum thickness, and the phonons have 1/2 probability of traveling upwards and 1/2 probability of traveling downwards. For the downward going phonons there is an probability of reaching the aluminum / crystal interface before scattering, where nm is a characteristic phonon interaction length BrinkThesis (). The factor of is provided to integrate over different phonon incidence angles. The factor is replaced by for upward going phonons.

If the phonon has not been removed from Monte Carlo in step 2 or reintroduced into the crystal in step 3 then the process repeats at step 1.
V Charge Monte Carlo
v.1 Introduction
Accurate modeling of charge propagation is included in Monte Carlo for numerous reasons. First, the ionization signal, compared to the phonon signal, provides a discriminator between electronrecoil and nuclearrecoil events in the silicon and germanium detectors. Second, electron transport is described by a mass tensor, leading to electron transport which contains components oblique to the applied field. This description is necessary to explain and interpret signals in the primary and guardring ionization channels, which function as a fiducial volume cut Sasaki1958 (); Jacoboni1983 (). Third, for electron recoils in the germanium bulk, charges drifting through the detector produce a population of Luke phonons which contribute 56% of the total phonon signal at 3 V bias. This fraction is understood by considering that for every of gamma energy, an electronhole pair is created which contributes of phonon energy; is the contribution of Luke phonons to the total phonon signal. These phonons’ spatial, time, energy and emitteddirection distributions should therefore be properly modeled in Monte Carlo of the detector response. Fourth, phonons created during electronhole recombination at the surfaces contribute 13% of the total phonon signal but in a low frequency, ballistic regime that is used to provide a surfaceevent discriminator. This fraction is understood by considering that for every of gamma energy, of phonon energy is released at the surface; is the contribution of electronhole gap energy phonons to the total phonon signal. These phonons also need to be properly modeled in a Monte Carlo.
Germanium has an anisotropic band structure described schematically in Figure 10 and shows energy band structure in which the hole ground state is situated in the band’s [000] direction and the electron ground state is in the Lband [111] direction. Hole propagation dynamics are relatively simple due to propagation in the band and the isotropic energymomentum dispersion relationship . Electron propagation dynamics are significantly more complicated due to the band structure and anisotropic energymomentum relationship. At low fields and low temperatures, electrons are unable to reach sufficient energy to propagate in the or Xbands, and are not considered necessary to consider in Monte Carlo. The electron energymomentum dispersion relationship is given by , where the longitudinal and transverse mass ratio .
This chapter will proceed by describing hole propagation and scattering, electron propagation and scattering utilizing a HerringVogt transformation and finally electronhole recombination. Higher order mass terms and scattering processes which occur at high electric fields are discussed elsewhere in the literature AubryFortuna2010 ().
v.2 Holes: propagation and scattering with isotropic bands and isotropic phonon velocity
Hole propagation dynamics are described by momentum evolution in an electric field and propagation in position space , where is the effective carrier mass.
Charge carriers cannot accelerate indefinitely however, and the shedding of NeganovLuke phonons Neganov1985 (); Luke1988 () limits their speed to around the longitudinal phonon phase velocity . As described in Figure 11, chargephonon scattering is an elastic processes, conserving both momentum (where and are the initial and final hole momentum vectors and is the phonon momentum vector) and energy (where and are the initial and final hole energies and is the phonon energy). The phonon energymomentum dispersion relationship is given by . Due to the low carrier energy, Umklapp processes Kittel () in which , where is a reciprocal lattice vector, are suppressed.
Energymomentum conservation coupled with the previous dispersion relationships leads to the final states and and
(18) 
where is the angular displacement between and , is the angular displacement between and , and is defined as . If we can determine a scattering rate and phonon angular displacement , then we can use these formulae to find the final states.
Fermi’s Golden Rule provides the transition probability per unit time per unit energy as
(19) 
For phonon emission processes,
(20) 
where is the deformation potential constant and is the phonon occupation number given by . A characteristic length can be defined as . The transition probability can be integrated over and to obtain a scattering rate
(21) 
The angular distribution then follows to be
(22) 
where .
The phonon scatter azimuthal rotation angle is uniformly distributed about . The charge carrier azimuthal rotation angle is required to be .
v.3 Electrons: propagation and scattering with anisotropic bands and isotropic phonon velocity
Electron propagation and scattering is complicated by the anisotropy of the electron bands but can be simplified by performing first a transformation into a space defined by the vectors (where is aligned with [111] and the other two are perpendicular) and then a HerringVogt transformation into a space where the electron bands are isotropic Herring1956 (); Jacoboni1983 (). The HerringVogt transformation is nonunitary and in the space is given by
(23) 
but the speed of sound remains unchanged and isotropic.
In this space, , , , and the effective mass is given by . The change in velocity, in position space, is found by a back transform, incorporating both the mass and momentum transforms, and is given by .
After the HerringVogt transform is used to find the electricfield induced acceleration, the electrons shed phonons via the same prescription given to the holes. The phonon and electron momentum is found first in the HerringVogt space and the back transform is applied to return to position space. The only additional concern is correctly back transforming the phonon momentum due to the nonunitarity of the HerringVogt transform. To handle this, we maintain the phonon momentum magnitude (ie, conserve energy), but use the back transform to find the correct angular distribution.
v.4 Charge Time Steps, First Order
Determining an efficient time step size, for charge transport is complicated since the scattering time varies as the charge carrier accelerates. In this charge transport model it was decided to use sufficiently small such that is relatively constant at each iteration. The requirement, sufficiently small, in practice means that energy is conserved in the transport process and the following is a detailed description of how this requirement is implemented.
By observing the charge Monte Carlo over numerous field strengths and in germanium, it was observed that scattering limits the maximum possible charge momentum magnitude to about
(24) 
(25) 
for electrons and holes respectively. These momenta are then used to determine stepping times which conserve energy; the shed Luke phonon energy must equal the change in potential energy at the carrier drifts. Again running numerous Monte Carlo it is determined that a stepping time of is sufficiently small to conserve energy. Larger stepping times will result in a deficit of NeganovLuke phonons being created.
v.5 Charge Time Steps, Second Order
The earlier described first order method can be improved upon by developing a second order method. The challenge is to efficiently and accurately determine the time until a charge sheds a NeganovLuke phonon, which is challenging due to the changing interaction time, as the charge is accelerated by the field. An inverse CDF technique would be advantageous and one is developed here that adapts to the changing interaction time.
This is done by sampling the scattering rate at two different times and . We start with the differential equation and expand to next order in time
(26) 
Integrating, we can obtain
(27) 
This continues with the standard technique of solving for the CDF and inverting to obtain which can be solved for the scattering time . The positive root is retained as physical which provides a scattering time of
(28) 
This is completed by recognizing that
(29) 
and
(30) 
There is a maximum sampling time step () that can be used before the linear interpolation of scattering rates is inaccurate. As in the first order case, this results in lack of energy conservation. It is found that the electron sampling time step can be a factor of 15 greater than the time step shown in Equation 24 and the hole sampling time step by a factor of 20 greater than the time step in Equation 25. Given these sampling time steps, most NeganovLuke phonons are produced in a time , hence this method is much more efficient than the first order method. The efficiency of this second order method also implies that pursuit of higher order methods will not yield much additional improvement in computational efficiency.
This method also couples well to a second order spatial transport method. The velocity form of the Verlet algorithm Verlet1967 (); Swope1982 () is convenient and given by a description that should look familiar.
This can be easily modified to incorporate the second order inverse CDF sampling method via the following procedure:

Make a guess for the step size , which we will call





From second order inverse CDF method, determine the randomly distributed scattering time

If




Save and for use in next iteration


Else, save and for use in next iteration
Since holes are described by a scalar mass, this procedure is straightforward. There is however a slight modification to this procedure that is useful for electrons. Due to the use of the Herring Vogt transform, it is generally easier to keep track of momentum in the Herring Vogt space, , rather than velocity. We also identify that the acceleration is given by , where the transform from the cartesian space to the space defined by basis vectors is shown explicitly.
v.6 Select constants for charge Monte Carlo
Silicon  Germanium  
Electrons  Holes  Electrons  Holes  
  0.5    0.35  
0.91    1.58    
0.19    0.081    
0.26    0.12    
(km/s)  9.0  5.4  
(g/cm)  2.335  5.323  
(eV)  9.0  5.0  11.0  4.6 
(eV)      11.0  3.4 
(m)  16.9  7.5  257  108 
Vi ElectricField Lookup
vi.1 ElectricField Lookup from Triangulated Mesh
A numerical electric field model is necessary for the charge transport described in Section V. The simplest model is a constant, longitudinally directed field. However it may be desirable to include fringing fields and details from the electrode structure. A more accurate model will utilize a triangulated mesh. A 3d mesh contains nodes, with each mesh node containing an associated electric potential . At points other than a mesh node, the potential must be interpolated. The MATLAB programming language offers a few options for this interpolation, the fastest of which, utilizing a barycentriccoordinates linear interpolation via the TriScatteredInterp class. The Computational Geometry Algorithms Library (CGAL) is available for C++
cgal () though this paper will be presented for a MATLAB implementation. Furthermore, the barycentric transformation involves a linear transformation which implies that the electric field is constant within a triangulation, a property which can be exploited to speed up computation.
MATLAB’s method of looking up the potential, while very convenient, is not efficient considering the number of repeated field queries that occur before the carrier has moved to a location with a differing field. Efficiency can be improved by exploiting the fact that a charge remains within its triangle for numerous iterations and that the field is constant with the triangulation. On the contrary, MATLAB solves for the potential at every iteration in its lookup procedure, which is a bit slow. These repeated searches can be avoided but at the expense of significant code complexity. The effort is justified however as charge transport imposes a dominant computational expense in Monte Carlo.
vi.2 Barycentric Coordinates
Given a triangulation (the mesh is made of tetrahedra in 3space but the term triangulation often persists) with four node points , and , the arbitrary point can be described by the barycentric coordinates , and where
(31) 
The additional constraint is imposed that
(32) 
The barycentric coordinates become more intuitive when thought of as area (volume in 3d) coordinates (see Figure 12). In this paradigm, consider the 2d node points , and along with the probe point . To start with, let’s normalize the area enclosed by , and to (this normalization is equivalent to the constraint ). Then we can consider the three different areas enclosed by 1) , and (), 2) , and () and 3) , and (). It turns out that these areas (, and ) are identically equal to the barycentric coordinates , and , providing a quick and intuitive interpretation of the barycentric coordinates. The process and interpretation is the same in 3d when volume is substituted for area. It is not actually recommended to calculate through this procedure but to instead follow the procedure in Section VI.3.
vi.3 Barycentric Coordinate Formulae
In this section we derive formulae useful for solving the barycentric coordinates and electricpotential. We start again with the definitions given by equations 31 and 32. After separating equation 31 into the , and components we can solve for the through the following linear procedures and the formula is written out explicitly below.
(33) 
(34) 
where . Explicitly we can write out
(35) 
where .
The potential is simple to solve for in the barycentric coordinate system and equal to . This procedure is more obvious when one considers the connection to area coordinates.
vi.4 Barycentric Coordinate Procedures and Shortcuts
The procedure for finding which tetrahedra a charge resides is not unique and a canned MATLAB procedure (or CGAL library in C++
) can be utilized for this step. After the tetrahedra in which the carrier resides is determined, the electric field, is computed. This can be performed by probing the potential at four locations ( and ) and computing gradients. The drawback is that this procedure requires conversion to barycentric coordinates and electric potential lookup for three additional points. These steps can be eliminated with some simple derivations that are outlined below. Most of these steps provide a conceptual framework and only the last step is actually computed.
First we consider two points and and find their associate barycentric coordinates and .
(36) 
and
(37) 
where, as always, is constrained by . Next we solve for the potentials and
(38) 
and