Cryogenic solid-state detectors are used in a number of dark matter (DM) search experiments [1][2][3][4][5][6]. In these experiments, incoming DM particles are expected to scatter off of the detector nuclei or electrons, creating phonon signals which are measured by high resolution phonon sensors. Resolution on the order of 1 eV is achieved, which allows for reduced energy thresholds and enables the detection of nuclear recoils with energies as low as ∼10 eV [1][2][3]. Low-mass DM candidates that produce small interaction energies can be probed via electron recoils by measuring the ionization signal-the number of produced e -h þ pairs in the detector [7,8]. In phonon-based crystal detectors, when a voltage bias is applied across the crystal, the ionization signal is converted into an amplified phonon signal via the Neganov-Trofimov-Luke (NTL) effect [9,10]. A charge carrier with a charge e accelerated by the electric field scatters off of the crystal lattice and produces NTL phonons with the total energy equal to the work done by the electric field to move the charge through the electric potential difference Δφ:
Normally, when an e -h þ pair is created in the crystal, each charge drifts in the electric field all the way to the corresponding electrode on the crystal surface. Together they traverse the entire voltage bias of the detector, so the total energy of the produced NTL phonons is given by
where n eh is the number of e -h þ pairs and V bias is the voltage bias. The total phonon energy produced in an event is then given by
where E dep is the energy deposited in the detector by the incoming particle. For a detector with a good phonon energy resolution σ res and a large voltage bias, where σ res ≪ eV bias , the spectrum of the phonon energy in Eq. ( 3) is expected to have quantized peaks corresponding to the integer number of created e -h þ pairs. This e -h þ -pair quantization is observed in SuperCDMS high-voltage (HV) and HV eV-scale (HVeV) detectors when operated with a voltage bias on the order of 100 V [7,8,11].
Due to the presence of impurities in the crystal, a nonquantized amount of NTL energy can be produced in an event. We distinguish two categories of effects causing nonquantized NTL energy: charge trapping (CT) and impact ionization (II). In a CT process, a charge carrier gets trapped in an impurity state in the bulk of the crystal. In an II process, a propagating charge ejects (or "ionizes") an additional unpaired charge from a shallow impurity state. Trapped charge carriers and unpaired charge carriers created in an II process terminate or start their trajectories in the bulk of the detector, respectively. As a result, they traverse only a fraction of the voltage bias, producing a nonquantized amount of NTL energy.
A proper modeling of these detector response effects is crucial for DM search analyses. In Sec. II, we develop an analytical model (the so-called "exponential CTII" model) that describes the NTL energy spectrum for events affected by the CT and II processes. We improve upon the previously used CT and II model introduced in Ref. [12] (the so-called "flat CTII" model) by taking into account the distribution of locations at which the CT and II processes occur. We demonstrate a difference between the NTL energy spectra of events produced on the detector surface and events produced in the detector bulk that can be used for statistical discrimination between surface background and bulk DM events. In Sec. III, we incorporate the CT and II model into the full detector response model, and take into account additional surface effects that may be relevant to certain calibration data. This modeling is used in Sec. IV to extract detector response parameters for HVeV detectors.
The underlying physical assumption of the exponential CTII model is that there are three possible processes that can occur to a charge carrier (an electron or a hole) when it traverses the bulk of the crystal under the influence of an electric field. It can get trapped in an impurity state, it can create a single free electron from an impurity state by promoting it into the conduction band, or it can promote an electron from the valence band to an impurity state, creating a single hole in the valance band. The probabilities for these processes to occur may differ between holes and electrons; therefore we consider in the model a total of six different CT and II processes: electron trapping ("CTe"), hole trapping ("CTh"), creation of a hole by an electron ("IIeh"), creation of an electron by an electron ("IIee"), creation of an electron by a hole ("IIhe"), and creation of a hole by a hole ("IIhh").
The model assumes that each of the six processes has a small constant probability of occurring at any point of the charge carrier's trajectory, independent of the location in the bulk of the crystal, of the path already traveled by the charge, and the presence of other charges simultaneously traversing the crystal. Additionally, it is assumed that charges propagate along some z axis that is parallel to a uniform electric field (detectors, including the HVeV detectors used in Refs. [7,8], are typically designed to have a uniform electric field throughout the bulk). We start by considering that impurities are distributed uniformly throughout the bulk of the crystal, where we let p i denote the probability for a charge to undergo a certain process i per unit of distance traveled along the z axis. Here, i refers to the specific CT or II process a charge may undergo (CTe, CTh, IIee, IIeh, IIhe, or IIhh). p i itself may depend on various factors, including the impurity density and the amount of charge diffusion. If a charge travels a distance Δz in n steps, the total probability of the charge not undergoing some CT or II process Ci ðΔzÞ is ð1p i Δz=nÞ n . In the limit of infinitesimally small step sizes, Ci ðΔzÞ becomes
where the p i term in Eq. ( 4) is replaced with 1=τ i , with τ i defining the characteristic length of that particular CT or II process. Ci ðΔzÞ is the complementary cumulative distribution function of the probability density function (PDF) that describes the probability for a charge to travel a distance Δz before a particular CT or II process occurs. This PDF is therefore given by
While these PDFs are described in terms of a distance traveled, the model also imposes the condition that the charges terminate when reaching the crystal surface.
That means for a crystal with a thickness Z, the charges are bound between z ¼ 0 and z ¼ Z. For convenience we choose Z ¼ 1, and let z describe the proportion of the crystal thickness rather than a physical distance. The six characteristic lengths τ i measured in fractions of the crystal thickness are the only fundamental input parameters of the model. We write these characteristic lengths in terms of probabilities f i defined as
Hence f i is the probability of a particular process occurring if a charge can traverse the entire length of the detector. Equations ( 5) and ( 6) are repeated for each of the six processes, and together make up the fundamental building blocks of our exponential CTII model.
The end product of the model is a PDF of the NTL energy produced in an event. This energy is proportional to the distance traveled by the charges along the field lines. We adopt an energy scale E neh such that a unit of E neh is equivalent to the amount of NTL energy produced by a charge that travels a distance equal to the thickness of the crystal. Using this energy scale, a charge going from z ¼ 0 to z ¼ 1, as well as an e -h þ pair starting at z ¼ 0.5 whereby both charges travel a distance Δz ¼ 0.5, will result in a total energy of E neh ¼ 1. With such units, there is a one-to-one correspondence between the PDFs of the total NTL energy and the total distance traveled by the charges along the electric field.
The exponential CTII model is constructed by finding the analytical solutions for the NTL energy produced by a single e -h þ pair for events of three distinct classes. The first are surface events, where a single charge is created at one of the surfaces (i.e. along the z ¼ 0 or z ¼ 1 plane) and propagates toward the opposite surface; this class of events does not include events created along the lateral surfaces of the crystal. Surface events correspond to laser or lightemitting diode (LED) calibration data, whereby optical photons are absorbed near the z ¼ 0 or z ¼ 1 surface of the crystal, as well as to charge leakage originating at the crystal surface. The second class of events are single charges produced throughout the bulk of the crystal. These events may correspond to some charge leakage process that happens throughout the detector bulk. The third class of events are bulk-e -h þ pairs produced throughout the bulk of the crystal. These events are what is expected for DM interactions. For each class of events, we consider various unique combinations of CT and II processes occurring to the charges, and solve for the probabilities of measuring an energy of E neh given those unique combinations of processes.
Modeling multiple II processes in a single event poses a significant challenge: Each additional II process allowed adds an new charge carrier, causing the number of potential combinations of CT and II processes to grow exponentially, and the complexity of each new solution greatly increases. For this reason, we limit the number of solutions to a certain "order" of processes, where the order of a process is defined as follows: For processes of order N, charges that participated or were produced in a primary II process can take part in no more than (N -1) additional II processes. For surface events and bulk-single-charge events, the solutions for processes up to second order are found, resulting in 28 unique solutions for each event type. For bulk-e -h þ -pair events, the solutions for processes up to first order are found, resulting in 16 unique solutions. When solving for these analytical solutions, we assume that any charges existing after the order limit is reached will propagate to a crystal surface with 100% probability. Appendix A provides a detailed description of how these solutions are found, with Appendixes A 1-A 3 adding further details on solving the solutions for each of the three classes of events. The full list of process combinations and the corresponding solutions are cataloged in the Supplemental Material [13] and are displayed in Fig. 1. It is immediately apparent how the computed PDFs differ for the different classes of events. Namely, the regions above and below the first e -h þ -pair peak are relatively flat for surface events, in contrast to bulk-e -h þ -pair events where the PDF in the same regions is more curved. Furthermore, the PDF for bulk-singlecharge events does not have a delta function at E neh ¼ 1.
The analytical solutions are found for when there is initially only a single charge or e -h þ pair produced. However large energy depositions in the crystal will often produce multiple charges or e -h þ pairs for a single event. Let F ð1Þ type ðE neh Þ be the probability distribution function for one charge or e -h þ pair in the n eh -energy space. The "type" refers to the specific event type to model: either surface events, bulk-single-charge events, or bulk-e -h þ -pair events. F ð1Þ type ðE neh Þ is found by summing the analytical solutions for the given event type, and examples of this function are shown by the black, dashed curves in Fig. 1. Without any additional detector response, the PDF for j e -h þ pairs F ðjÞ type ðE neh Þ is found by convolving F ð1Þ type ðE neh Þ with itself (j -1) times:
In practice, F ðjÞ type ðE neh Þ is found using numerical convolution. We can use this to construct the PDF for events that generate multiple e -h þ pairs, defined as HðE neh Þ. The solution for HðE neh Þ up to J e -h þ pairs is given by
where a j are the weights associated with producing j e -h þ pairs, which are discussed more in Sec. III. A comparison of the PDFs for single-e -h þ -pair events F ð1Þ type ðE neh Þ and multi-e -h þ -pair events HðE neh Þ for different event types is shown in Fig. 2 for arbitrary CT and II probabilities. The PDFs are convolved with a Gaussian function to emulate the energy resolution for illustrative purposes. Furthermore, the solutions are compared to the PDFs computed using the flat CTII model described in Ref. [12].
The example PDFs from Fig. 2 allow us to make some broad observations about the exponential CTII model. First, while the higher-order processes are significant for the singlee -h þ -pair solutions (as seen in the top plot of Fig. 2 above E neh ¼ 2), they generally become less significant or even negligible for multi-e -h þ -pair solutions. Second, the type of events being modeled has a significant impact on the shape of the PDFs between the e -h þ -pair peaks. Notably, the between-peak shape for the bulk-e -h þ -pair events differs greatly from that of surface events, as well as that of the flat CTII model which does not differentiate between event types.
A. Single-and multihit solutions Equation (8) describes the PDF of producing a certain amount of NTL energy for a given event type that generates FIG. 2. Example PDFs found for single-e -h þ -pair events F ð1Þ type ðE neh Þ (top) and multi-e -h þ -pair events HðE neh Þ (bottom). The PDFs are computed for surface events (solid, blue curves), bulk-e -h þ -pair events (dashed, green curves), and bulk-singlecharge events (dot-dash, orange curves) using the exponential CTII model. For comparison, the PDFs computed using the flat CTII model from Ref. [12] are shown by the dotted, purple curves. These examples are shown for arbitrary CT and II parameters:
for the flat CTII model, f CT ¼ 20% and f II ¼ 4%. Furthermore, the multi-e -h þ -pair solutions assume that the a j terms in Eq. ( 8) that describe the e -h þ -pair probabilities follow a Poisson distribution with a mean of two e -h þ pairs. For illustrative purposes, the PDFs are convolved with a Gaussian function with a width of E neh ¼ 0.05 to emulate the detector energy resolution. FIG. 1. Analytical solutions in the E neh energy space of the exponential CTII model for single-e -h þ -pair events. The unique solutions represented by the solid, colored curves are found for surface events (top), bulk-single-charge events (middle), and bulke -h þ -pair events (bottom). This example is shown for arbitrary values of the CT and II parameters:
and f IIhh ¼ 5%, and the top and middle plots assume that the initial charge is an electron. The black, dashed curves in each plot are the sums of the analytical solutions for each event type and are examples of
multiple e -h þ pairs, HðE neh Þ, which is derived from the analytical solutions of the exponential CTII model. However HðE neh Þ can be extended to include other detector response effects that are either measured or expected. These additional effects include ionization probabilities, conversion to the phonon energy scale, and continuous spectra of energy deposition. To start, we define H ð1Þ as the PDF for events resulting from a single interaction of a particle with the crystal, which we call "single-hit" events. Examples of a single-hit event include a single photon absorbed by the crystal, or a single DM particle scattering off of an electron. First we will construct a general formula for H ð1Þ , and then subsequently see how this formula is used to model specific interactions from various sources.
The first step to extend the detector response model is to replace the generic weights a j in Eq. ( 8) with the probability mass function (PMF) describing the probability of producing a given amount of ionization. The ionization PMF is specific to the detector material, and is a function of the energy deposited in the detector E dep . Let p eh ðjjE dep Þ describe the ionization probability of producing j e -h þ pairs given E dep . For silicon, results of the ionization yield at low energies can be found in Ref. [14].
Next we need to convert the PDFs to the correct energy scale. As mentioned in Sec. I, event energies are measured by the total phonon energy E ph described by Eq. ( 3), whereas the F ð1Þ type functions of the exponential CTII model are described in the n eh -energy space E neh . Using Eq. ( 3), E neh can be written in terms of E ph as
This change in energy units also changes the overall scaling of the PDFs. To account for this, the PDFs must be scaled by a factor of jdE neh =dE ph j ¼ 1=eV bias . Finally, we need to consider the general case where there is a continuum of energy depositions that can occur for a given source of events. This continuum can be described by a normalized differential rate spectrum d R=dE dep ðE dep Þ, where for a total single-hit event rate of
Putting this all together, the extended detector response model for single-hit events in the phonon energy space modeled up to J e -h þ pairs is given as
Here we assume that J is large enough such that the ionization PMF sums to unity for all E dep . As the d R=dE dep ðE dep Þ function in Eq. ( 10) is normalized, H ð1Þ ðE ph Þ is describing a PDF for single-hit events from a given source. We also consider so-called "multihit" events, which are events generated from simultaneous particle interactions in the detector. In general, the PDF solutions for multihit events are found by recursively convolving the single-hit solution from Eq. (10). An example of constructing a multihit PDF solution is shown in Sec. III A 1.
Up to this point, the detector response model has been described without considering the detector energy resolution σ res . While σ res can be incorporated into the model in different ways, this work assumes that the energy resolution is constant over E ph . Therefore the single-hit model including the energy resolution H ð1Þ ðE ph ; σ res Þ can be expressed as
where GðE ph jμ ¼ 0; σ res Þ is a Gaussian function with a mean μ ¼ 0 and width of σ res . One could also consider, for example, an energy resolution that depends on n eh . In which case, a Gaussian function with a width equal to the energy resolution at each e -h þ -pair peak σ ðjÞ res would be convolved with the corresponding F ðjÞ type function.
A typical way to calibrate the energy of HVeV detectors is to use a photon source. Specifically in Refs. [7,8], a laser source of optical photons was pointed at one of the detector surfaces. The laser was pulsed at some known frequency f γ , and, depending on the laser intensity, produced an average number of photons per pulse λ that were detected. The probability of a given number of photons per laser event is given by a Poisson distribution with a mean of λ. These photon-calibration events are therefore examples of multihit events, where the probability distribution must also account for the probability of the simultaneous absorption of multiple photons in a single event. The differential rate dR=dE ph ðE ph Þ for photon-calibration events is then given by
where H ðlÞ ðE ph Þ corresponds to the NTL energy produced in an event with l ≤ L photons absorbed. H ð0Þ ðE ph Þ corresponds to events with no photons absorbed. Such events may be present in the calibration data if the detector trigger is synchronized with the laser pulses and λ is small. In the simplest scenario, H ð0Þ ðE ph Þ ¼ δðE ph Þ. However the zeroth-e -h þ -pair peak can also take a more complex form, like the modified Gaussian noise peak described in
Ref. [12]. H ð1Þ ðE ph Þ corresponds to events with one photon absorbed and is generally given by Eq. ( 10).
For photon-calibration events, we assume that the photons are absorbed sufficiently close to a detector surface such that these events can effectively be modeled as surface events created along the z ¼ 0 or z ¼ 1 plane. Furthermore for a spectrum of photon energies E γ , d R=dE dep ¼ d R=dE γ . For an LED source, d R=dE γ is the normalized energy spectrum of the LED photons. Yet for a laser source like in Refs. [7,8], the photons are monoenergetic, and therefore d R=dE dep ¼ δðE dep -E γ Þ. Moreover the laser used for the calibration in Refs. [7,8] produced 1.95 eV photons which, for silicon, always ionize exactly one e -h þ pair per absorbed photon [14]. Photons of this energy have an absorption length in silicon of Oð10 μmÞ which, for a detector that is 4 mm thick [8], is sufficiently small to model these events as surface events. For this particular case, Eq. ( 10) reduces to
Again if we assume a constant energy resolution σ res , Eq. ( 12) is convolved with a Gaussian function with a mean μ ¼ 0 and a width of σ res in order to compute dR=dE ph ðE ph Þ with the energy resolution. The distinction between singlehit and multihit events displayed here is subtle yet important. For the case of low-energy photon-calibration events, the multi-e -h þ -pair peaks observed are not due to multiple e -h þ pairs ionized from a single absorbed photon, but rather simultaneously absorbed photons that each ionize a single e -h þ pair.
For any dark matter search experiment, a detector response model is required to determine the expected signal distribution of a DM candidate in the detector. Therefore Eq. ( 10) can also be used to compute expected DM signals in HVeV detectors. Unlike photon-calibration events, DM interactions are considered to be single-hit events; generally DM models exclude the possibility of the simultaneous interaction of multiple DM particles with a detector. While the exact signal distribution will depend on the specific DM candidate that is modeled, we will look at examples of two DM candidates commonly searched for using HVeV detectors.
The first candidate is the dark photon that is modeled, for example, in Ref. [15]. In this model, nonrelativistic dark photons with a mass m A 0 constitute all relic dark matter. The interaction rate of dark photon absorption R A 0 ðm A 0 ; εÞ depends on its mass and is proportional to the kinetic mixing parameter ε that couples dark photons to standard model photons. In this model dark photons provide a monoenergetic source of energy deposition equal to its mass such that d R=dE dep ¼ δðE depm A 0 c 2 Þ, where c is the speed of light. Substituting this into Eq. ( 10) and noting that DM interactions are modeled as bulk-e -h þ -pair events, the differential rate of dark photon absorption dR
The second candidate we consider is light DM that elastically scatters off of electrons, as described in Ref. [16]. In this model, the dark matter particle χ with mass m χ is also assumed to constitute all relic DM, and scattering interactions with electrons are mediated via a dark-sector gauge boson. The total rate of DM-electron scattering interactions R χ ðm χ ; σe Þ is dependent on the DM mass as well as the effective DM-electron scattering cross section σe . However, this DM-electron scattering process produces a spectrum of recoil energies E r . Specifically in Ref. [16], the recoil spectra are provided as rates over discrete recoil energy bins. Therefore the integral of d R=dE dep in Eq. ( 10) is replaced by a sum over weights w k corresponding to the recoil energies E ðkÞ r , where the weights are normalized such that
The different rate functions in Eqs. ( 14) and (15) do not yet include the detector energy resolution. As before we assume that σ res is constant over E ph , and therefore the energy resolution is incorporated by convolving Eqs. ( 14) and ( 15) with a Gaussian function with a mean μ ¼ 0 and a width of σ res .
The detector response model can be extended further by modeling other phenomena that are observed in the detector. One such phenomenon is the apparent deposition of nonionizing energy measured together with photoncalibration events. We surmise that this detector response effect occurs in HVeV detectors because of the observed dependence of the e -h þ -pair peak positions on λ, the average number of photons per laser or LED pulse [17]. One hypothesis is that some proportion of photons is absorbed directly into the aluminum fins of the phonon sensors. Another hypothesis is that, due to the random initial trajectory of electrons and holes after ionizing, there is some probability that both charges will happen to recombine at the nearest detector surface. This so-called surface-trapping effect has been observed in detector simulations using G4CMP [18,19], and is illustrated in the top plot of Fig. 3. In any case, these hypotheses suppose that some proportion of photons will deposit some nonionizing energy without generating a typical e -h þ pair that undergoes the bulk CT and II processes.
In the case of the hypothesized surface-trapping effect, we can include this effect in the model by modifying the single-hit PDF for photon-calibration events described by Eq. ( 13). Let α be the probability of the created e -h þ pair to undergo surface trapping, where 0 ≤ α ≤ 1. That means there is a (1 -α) probability that the e -h þ pair will propagate through the crystal, undergoing the typical bulk CT and II processes. For photons that undergo surface trapping, the deposited energy will only be the absorption energy of the photon E γ . We are then able to include the surface-trapping effect by modifying Eq. ( 13) in the following way:
The multihit solution for photon-calibration data with the surface-trapping effect is given by Eq. (12), where H ðlÞ ðE ph Þ is found by recursively convolving Eq. ( 16) with itself. The result of including the surface-trapping effect in the detector response model is illustrated in the bottom two plots of Fig. 3. Due to the presence of nonionizing photons, each peak in the spectrum splits into multiple subpeaks separated by E γ , as seen in the middle plot of Fig. 3 before resolution smearing. Each subpeak corresponds to q ionizing photons and p nonionizing photons, with the subpeak location defined as q • eV bias þ ðq þ pÞ • E γ . Note that a FIG. 3. Top: illustration of the hypothesized surface-trapping effect as observed from simulation data using G4CMP [18,19]. Two examples are shown of the trajectory of an ionized e -h þ pair in terms of the depth below the detector surface and the perpendicular x-coordinate relative to the hit position of the absorbed photon. The right example shows a typical event, where the electron eventually travels in the direction opposing the electric field. The left example shows a surface-trapped event, where the electron recombines with the detector surface before it can turn around. Middle and bottom: examples of modeling the surface-trapping effect using Eqs. ( 12) and ( 16) with V bias ¼ 100 V, E γ ¼ 1.95 eV, f γ ¼ 1 Hz, and arbitrary CT and II probabilities. The additional spikes seen in the middle plot demonstrate the contribution of nonionizing energy deposition when α > 0 which, when smeared by the energy resolution, widen and shift the e -h þ -pair peaks. This effect also causes a peak position dependence on λ, as demonstrated by the bottom plot with α ¼ 0.3. subpeak corresponding to p and q is part of the function H ðqþpÞ ðE ph Þ, rather than of the function H ðqÞ ðE ph Þ. Therefore, in order to properly model the substructure of the qth e -h þ -pair peak, modeling of higher e -h þ -pair peaks is required. When normalized, the underlying amplitudes of the subpeaks in each e -h þ -pair peak follow a Poisson distribution of the number of nonionizing photons with a mean of α • λ (see Appendix B for more details). To include all the significant subpeaks, it is recommended to set the maximum number of modeled peaks [L in Eq. ( 12)] to a number exceeding the number of peaks in the region of interest by the mean number of nonionizing photons plus at least 3 standard deviations of its distribution, i.e. by ðα
When the energy resolution is applied, the peak substructure from nonionizing photons gets smeared and appears as a shift and a widening of the e -h þ -pair peaks, as shown in the bottom two plots of Fig. 3. When λ increases, there is a greater contribution from nonionizing photons, resulting in wider peaks that are shifted by a larger amount.
To demonstrate the performance of the exponential CTII and extended detector response model described in Secs. II and III, we fit the model to laser-calibration data acquired from Ref. [8]. Specifically, the model for photoncalibration events described by Eqs. ( 12) and ( 13) is fit to laser-calibration datasets from Ref. [8] acquired with V bias ¼ 100 V and E γ ¼ 1.95 eV. The individual datasets differ by the laser intensity used during data acquisition, and thus by the value of λ. The fits of the model to two of these datasets are shown in Fig. 4. For simplicity, we reduced the number of parameters in the fits by requiring
We obtain the best-fit results from the fit to the left (right) laser-calibration dataset in Fig. 4 of λ ¼ 0.41 AE 0.01 (0.475 AE 0.005), σ res ¼ 3.30 AE 0.04 eV (3.37 AE 0.02 eV), f CT ¼ 11.6 AE 0.6% (12.1 AE 0.3%), and f II ¼ 0.7 AE 0.4% (0.9 AE 0.2%). Within uncertainties, these results are consistent with the results obtained by fitting the flat CTII model from Ref. [12] to the same datasets. The consistency of the results is expected, as the flat CTII model has previously demonstrated that it can accurately describe photon-calibration data [20]. Figure 4 shows that the exponential CTII model is able to obtain equivalent results for this relatively simple scenario. However as will be shown below, the advantages of the model presented in this work become apparent when including additional detector response effects or when modeling different event types.
We can further evaluate the extended detector response model by performing a simultaneous fit of the model to multiple photon-calibration datasets. To do this, the simultaneous fit to multiple datasets is done separately for data acquired from two different experiments. The first are three laser-calibration datasets acquired in Ref. [8]. The second are three LED-calibration datasets acquired at the Northwestern EXperimental Underground Site (NEXUS) at Fermilab (Batavia, IL). This NEXUS facility is located in the NUMI tunnel, which provides an overburden of 225 mwe [21], and hosts a Cryoconcept dry dilution refrigerator. The LED-calibration data reported in this work were acquired by one of four 1-cm-side HVeV detectors that were operated at NEXUS between May 14 and July 27, 2022. More information about the experiment design, data acquisition, and data analysis can be found in Ref. [22].
There are several key similarities and differences between laser-calibration datasets acquired in Ref. [8] and the LED-calibration datasets acquired at the NEXUS facility reported in this work. In both cases, the data were acquired using an HVeV detector with an "NF-C" sensor design [11]. Both devices are constituted by a 10 × 10 × 4 mm 3 silicon target with two channels of quasiparticletrap-assisted electrothermal-feedback transition-edge sensors [23] (QETs) patterned on the top surface to measure the phonon signal. While the HVeV detector used to acquire the NEXUS data is not the same as the one used in Ref. [8], the substrate from both detectors belongs to the same silicon wafer. This means that the impurity levels in both detectors are likely to be similar.
Furthermore, both detectors generated an electric field throughout the bulk of the crystal by applying a high voltage to an aluminum electrode deposited on the detector surface FIG. 4. Results of fitting the exponential CTII and extended detector response model to two laser-calibration datasets from Ref. [8] acquired with V bias ¼ 100 V and E γ ¼ 1.95 eV. The residuals from the fit results are shown in the bottom plots.
opposite the surface patterned with the QETs; the QET surface of the detectors was kept grounded. In Ref. [8], lasercalibration data were acquired by emitting 1.95 eV photons from a laser onto the center of the QET face of the detector. In contrast, the LED-calibration data from the NEXUS facility were acquired using an ∼2 eV LED collimated on the center of the electrode face of the detector. Yet the data from both detectors were acquired with V bias ¼ þ100 V. This means that for the laser-calibration data from Ref. [8], the initial propagating charges are electrons, whereas for the LEDcalibration data from the NEXUS facility, the initial propagating charges are holes. Moreover, because the LEDcalibration data from the NEXUS facility illuminated the electrode face of the detector, any nonionizing energy deposition due to photon absorption into the aluminum fins of the QETs is expected to be minimal, especially compared to the laser-calibration data from Ref. [8].
We fit the extended model assuming nonionizing energy deposition caused by surface trapping [Eq. ( 16)] simultaneously to three laser-calibration datasets from Ref. [8] and three LED-calibration datasets from the NEXUS facility, all acquired with V bias ¼ 100 V. Each fit includes the parameters λ 1 , λ 2 , and λ 3 corresponding to the λ value of each dataset, but includes only one value of f CT , f II , σ res , and α for all datasets. For simplicity, we again reduced the number of parameters in the fit by imposing the requirements given by Eq. ( 17).
In the fit to the data from Ref. [8], we kept the energy of the laser photons fixed at E γ ¼ 1.95 eV, whereas in the fit to the NEXUS datasets, we allowed the energy of the LED photons to float. Furthermore, a measurement of the LED wavelength spectrum at 4 K found the spread in photon energies to be ∼0.0012 eV, and therefore we can adequately treat the LED as a monoenergetic source of photons described by Eq. (13). For the NEXUS datasets, we additionally included parameters to calibrate the data. The calibration converts the pulse amplitude A (in units of μA) of each event to the total phonon energy E ph . The fit includes three calibration parameters c 0 , c 1 , and c 2 that follow the equation
where the quadratic coefficient c 2 is included to account for any saturation effects in the QET sensors that can cause a nonlinear response at higher energies [11].
The top and bottom plots of Fig. 5 show the fit results to the datasets from Ref. [8] and the dataset acquired at the NEXUS facility, respectively. The best-fit results of the fit parameters are listed in Table I. As the CT and II probabilities of the initial propagating charge have the largest impact on the expected signal for photon-calibration events, we can interpret the values of f CT and f II from the fits to the datasets from Ref. [8] and the datasets acquired at the NEXUS facility as the CT and II probabilities for electrons and holes, respectively. Therefore, these results suggest that for these detectors (that come from the same silicon wafer), the CT probability for electrons may be higher than for holes. By using Eq. ( 6) and knowing that the thickness of these detectors is 4 mm, the fitted f CT values in Table I can be converted to the characteristic lengths of CT, giving 27.6 AE 0.4 mm and 32.7 AE 0.6 mm for electrons and holes, respectively. FIG. 5. Results of simultaneous fits of the extended detector response model to three laser-calibration datasets from the HVeV detector in Ref. [8] (top) and three LED datasets acquired from a HVeV detector at the NEXUS facility (bottom). All of the datasets were acquired with V bias ¼ 100 V. The model assumes nonionizing energy deposition caused by surface trapping as described by Eq. ( 16). The inset plots show the data and fit enlarged around the first e -h þ -pair peak in order to clearly observe the peak shifts between datasets.
In both cases, the fit determined the amount of surface trapping to be ∼40%. Indeed the inset plots in Fig. 5 that are enlarged around the one e -h þ -pair peak clearly show the peak position dependence on λ-a feature predicted when assuming surface trapping. Notably, the surfacetrapping probability found for the NEXUS data (where the LED photons are illuminated on the electrode face of the detector) is slighter higher compared to the data from Ref. [8]. This strongly disfavors the hypothesis that nonionizing energy deposition occurs due to photons being absorbed directly into the aluminum fins of the QETs. While this result supports the surface-trapping hypothesis, we stress that it is just one interpretation of these data. Additional dedicated measurements are needed to confirm these detector response effects, as well as to understand the differences that these effects have on electrons and holes.
We can additionally demonstrate the differences in the detector response for expected DM signals which, as mentioned in Sec. III A, are modeled as bulk-e -h þ -pair events. To do this, we ran simulations using G4CMP [18,19] of two DM models by generating events within the bulk of a silicon HVeV detector with V bias ¼ 100 V. The first simulated signal is dark photon absorption following Ref. [15] with a dark photon mass m A 0 ¼ 10 eV=c 2 and kinetic mixing parameter ε ¼ 10 -12 , and the second is DMelectron scattering following Ref. [16] with a DM mass m χ ¼ 5 MeV=c 2 , effective DM-electron scattering cross section σe ¼ 10 -33 cm 2 , and a DM form factor of F DM ¼ 1. The total number of events in both simulations were determined assuming an exposure of 6 gram days, and the ionization PMFs in the simulations are computed using the binomial approach taken in Refs. [7,8]. Finally, we assumed an energy resolution of σ res ¼ 3 eV and the following CT and II probabilities:
The results of these simulations are shown in Fig. 6, which also shows the expected signal for the two DM models computed using the exponential CTII and extended detector response model from this work following Eqs. ( 14) and (15). It is important to note that G4CMP-based simulations model CT and II processes using the same PDF described by Eq. ( 5) and are parametrized using characteristic lengths [18]. The consistency between the solid, orange curves and simulated data shown in Fig. 6 therefore provides a verification of the analytical solutions found for the various CT and II processes modeled. For comparison, Fig. 6 also shows the expected DM signals computed using the flat CTII model from Ref. [12] which, unlike the exponential CTII model, does not distinguish between surface and bulk event types. Differences in the energy spectra for different event sources are evident by comparing Fig. 6 with the energy spectra seen in Figs. 4 and 5. Our modeling expects the signal shape in the between-peak regions to differ for a bulk-e -h þ -pair source of events, such as DM, compared to surface events, such as photoncalibration data. The signal shape in the between-peak regions has a lot more curvature for bulk-e -h þ -pair events, in contrast to the relatively straight signal shape between the peaks for surface events. This feature is not captured by the flat CTII model, as evident in Fig. 6.
The exponential CTII model introduced in Sec. II addresses several limitations of the flat CTII model from Ref. [12]. This more robust model adopts a more physically motivated approach to describe CT and II processes in phonon-based crystal detectors. Consequently, it effectively characterizes the detector response across a range of event types and allows for the differentiation of CT and II probabilities for electrons and holes. These advantages, TABLE I. Best-fit results found by fitting the extended model assuming nonionizing energy deposition caused by surface trapping [Eq. ( 16)] simultaneously to multiple datasets. The results in the left column are determined from the fit to three laser-calibration datasets from Ref. [8], and the results in the right column are determined from the fit to three LED-calibration datasets from the NEXUS facility presented in this work. The fit to the NEXUS data also allowed the energy of the LED photons E γ to float, and included the calibration parameters c 0 , c 1 , and c 2 .
Laser data (Reference [8])
LED data (NEXUS)
1 3.18 AE 0.04 2.2 AE 0.1 λ 2 0.66 AE 0.02 4.2 AE 0.2 λ 3 2.90 AE 0.03 5.7 AE 0.2 f CT (%) 13.5 AE 0.2 1 1 .5 AE 0.2 f II (%) 0.4 AE 0.2 0.0 þ0.3 -0.0 σ res (eV) 3.41 AE 0.02 2.71 AE 0.06 α (%) 36.9 AE 0.7 4 1 AE 2 E γ (eV) 2.02 AE 0.07 c 0 (μA) ð1.7 AE 0.2Þ × 10 -3 c 1 (μA=eV) ð9.102 AE 0.006Þ × 10 -4 c 2 (μA=eV 2 ) ð-242 AE 7Þ × 10 -10 coupled with the expanded detector response model detailed in Sec. III, provide a more accurate representation of the detector's response to different sources of events.
The exponential CTII model combined with the extended response model assuming nonionizing energy deposition caused by surface trapping is shown in Fig. 5 to provide an accurate description of laser-calibration data from Ref. [8] and LED-calibration data acquired from the NEXUS facility reported in this work. While these results are encouraging, we note that this model may not encompass all relevant detector response effects that may occur. Rather, these results motivate us to acquire future measurements in order to continue investigating the full extent of these effects. Apart from understanding detector response effects, the result of the fit to the NEXUS data shown in Fig. 5 and Table I additionally demonstrates how the extended detector response model can be utilized to calibrate the detector. Incorporating these additional detector response effects can help to improve the accuracy of the energy calibration.
Furthermore, Fig. 6 illustrates how the signal shape differs for DM models, where the detector response is instead modeled by bulk-e -h þ -pair events. Indeed the analytical solutions to the exponential CTII model presented in this work match the spectra of events simulated using G4CMP, whereby CT and II processes are also parametrized using characteristic lengths [18]. The ability of the exponential CTII model to describe different event types provides a large advantage over the flat CTII model from Ref. [12], and can enable statistical discrimination between different background signals and expected DM signals.
Yet it is equally important to outline the limitations of the exponential CTII model. Physically, this model provides no description of the underlying mechanisms of CT or II processes, and simply assumes that there is some probability that these processes occur due to impurities throughout the crystal bulk. The characteristic lengths of these processes, τ i , could depend not only on the density of impurities, but also on the strength of the electric field, prebiasing history, "baking" history (impurity neutralization by detector irradiation), and temperature [24,25]. Furthermore, as mentioned in Sec. II, the exponential CTII model limits the analytical solutions to second-order processes for surface and bulk-single-charge events, and first-order processes for bulk-e -h þ -pair events. This limitation is a practical necessity, as including higher-order processes would exponentially increase the number and the complexity of solutions to solve for. However, the probability of II has been found in measurements to be of the order of 1% [8,20]. Therefore second-and third-order processes are expected to be extremely subdominant with probabilities ≪ 0.01%. Nevertheless, this limitation is quantitatively assessed for each event type in Appendix C.
While the detector response effects described in this work can also be modeled using Monte Carlo (MC) simulations, the exponential CTII model also provides a considerable computational advantage. For instance, calculating a DM exclusion limit requires the calculation of the expected DM signal, which itself may need to be computed for a large selection of CT or II values. The analytical solutions found for the exponential CTII model provide a quick and easy method for generating multiple DM signals with different CT or II probabilities compared to the computationally intense method of running many MC simulations.
The bulk CT and II processes modeled in this work are also present in other solid-state DM experiments, including SENSEI and DAMIC-M [26,27]. These experiments use charge-coupled devices (CCDs) that read out the amount of charge collected at each pixel. As such, charges that become trapped within the bulk of the crystal result in a signal loss. Furthermore, because charge packets are FIG. 6. Simulated data of events generated in the bulk of a silicon HVeV detector with V bias ¼ 100 V for DM signals using G4CMP [18,19]. The simulations were run for (top) dark photon absorption following Ref. [15] with m A 0 ¼ 10 eV=c 2 and ε ¼ 10 -12 and (bottom) DM-electron scattering following Ref. [16] with m χ ¼ 5 MeV=c 2 , σe ¼ 10 -33 cm 2 , and F DM ¼ 1. For running the simulations, we assumed an exposure of 6 gram days, an energy resolution of σ res ¼ 3 eV, and CT and II probabilities of
The solid, orange curves are the expected signals computed using our model following Eqs. ( 14) and (15). For comparison, the dashed, purple curves are the expected signals computed using the flat CTII model from Ref. [12].
required to move from pixel to pixel, the drift length of charges is typically larger for CCD detectors compared to HVeV detectors. The proposed Oscura experiment aims to account for this signal loss by considering that these trapped charges may be released at a later time and measured as single-electron events [28]. This bulk trapping is distinct from charge transfer inefficiency that occurs in CCD detectors, whereby charges are lost to surrounding pixels as the charge packet is moved from pixel to pixel. As shown in this work, the advantage of using phonon-based crystal detectors, including HVeV detectors, lies in the ability to exploit the nonquantized e -h þ -pair peaks regions on the energy spectrum to extract CT and II parameters as well as to differentiate among different sources of events.
Future plans involve dedicated detector response investigations using HVeV detectors. For example, taking CT and II measurements using crystals of different impurity levels while varying the voltage bias or the amount of prebiasing applied to the crystals can help develop our understanding of the factors that contribute to CT and II. Taking measurements with different voltage polarities will allow us to probe the differences in CT and II processes for electrons and holes. Furthermore, we aim to explore additional detector response effects and phenomena, including sources of nonionizing energy deposition. While these proposed measurements can be made using a photoncalibration source, finding a source of low-energy, bulke -h þ -pair events would additionally allow us to investigate the detector response expected for DM signals. The nuclear-recoil ionization yield measurement in Ref. [29] demonstrates a method of producing such events by measuring low-energy neutron recoils off of the nuclei of an HVeV detector. The model presented in this work not only provides a more robust understanding of the detector response effects in phonon-based crystal detectors, but can be utilized to help discriminate among different sources of events in order to improve the sensitivity of DM search experiments.
This appendix provides details on how to find the unique solutions for single-e -h þ -pair events. As mentioned in Sec. II, the events are categorized into three distinct classes: surface event, bulk-single-charge events, and bulk-e -h þpair events. Recall that both propagating electrons and holes have three processes which can occur: charge trapping (CTe and CTh), creation of a charge of the same kind (IIee and IIhh), and creation of a charge of the opposite kind (IIeh and IIhe). When writing the equation for the probability of a particular scenario occurring, we must consider the probabilities of all possible processes. For example, consider the probability of CT for a propagating electron as a function of Δz, P CTe ðΔzÞ. The probability of CTe occurring at a distance Δz traveled is the probability of CTe occurring at Δz and IIee not occurring by Δz and IIeh not occurring by Δz. Using Eqs. ( 4) and ( 5), this scenario is described by
where we define
Equivalent equations to Eq. (A1) can be found for the other five unique processes and by defining
We also need to consider the probability of a charge propagating a distance Δz without any CT or II process occurring. In this context, it is evaluating the probability of a charge reaching one of the crystal surfaces after traveling a distance Δz. We denote this probability as P S ðΔzÞ. For a propagating electron, this probability is given by
Equivalently, the probability for a propagating hole to reach a surface without a CT or II process occurring as a function of Δz is P S ðΔzÞ ¼ e -ΔzT h . Equation (A4) says that if a charge travels the full length of the detector (i.e. Δz ¼ 1), the probability that it does not undergo a CT or II process is e -T e=h .
So far, we have only described probabilities of propagating charges undergoing CT and II processes as a function of the distance traveled. However, detectors do not directly measure the distance charges are able to travel. As mentioned in Sec. II, the NTL energy measured by the detector is proportional to the distance traveled by the charges. In this parametrization, the energy measured due to ionization is therefore equal to the total distance traveled by all propagating charges involved with an event. The probability of measuring some energy E (where E is in n eh -energy space denoted as E neh is Sec. II) can be described as the sum of probabilities whereby the total distance traveled by all charges z tot is equal to E. While there is no direct constraint on z tot , and thus E, individual charges are constrained by the bounds of the crystal surfaces (i.e. 0 ≤ z ≤ 1).
Using the understanding of how the measured energy is related to the distance traveled by the propagating charges, we can start by solving for the almost trivial solutions, which will also create a set of base equations in which all other solutions can be found. We consider the solutions for charges propagating toward the z ¼ 0 and z ¼ 1 surfaces separately. For a charge traveling toward the z ¼ 1 surface that starts at a position z 0 , the total energy that can be measured by the charge is E ¼ 1z 0 . Conversely for a charge traveling toward the z ¼ 0 surface that starts at a position z 0 , the total energy that can be measured by the charge is E ¼ z 0 . Using Eq. (A4), we can write the probability distribution for a charge reaching the z ¼ 0 and z ¼ 1 surfaces as
where the subscript q ¼ ðe; hÞ indicates if the charge is an electron or hole, and the superscripts 0 and 1 indicate which surface the charge is traveling toward. Next, we consider the probability distribution of measuring an energy E before some CT or II process occurs. For now, whether the process is CT or II does not matter. These base solutions are found using Eq. (A1) as a framework and solved separately for charges propagating toward the z ¼ 0 and z ¼ 1 surfaces. For charges that undergo some process i and start at a position z 0 , these probability distributions P i ðE; z 0 Þ go as
where again the superscripts 0 and 1 indicate the direction of propagation. Some specific solutions can immediately be found from Eq. (A6). Specifically when a charge undergoes CT, the charge can no longer propagate and produce more energy or undergo additional processes. Therefore Eq. (A6) are also the solutions for the probability distribution of a CT process for a charge starting at a position z 0 . This is decisively not the case for a charge that undergoes II, as both the initial charge and the newly created charge will continue to traverse the crystal and increase the energy measured. Nevertheless, Eqs. (A5) and (A6) provide the necessary foundation for calculating the probability distribution for any specific scenario. Equipped with the base equations of the probability distributions in energy for propagating charges traveling in either direction with some starting position, we can now analytically solve for the single-e -h þ -pair probability distributions corresponding to specific events and scenarios.
A surface event starts with the creation of either an e -or a h þ at the starting position z 0 ¼ 0 or z 0 ¼ 1. Knowing the polarity of the voltage bias and the starting position will necessarily decide whether the propagating charge is an e - or a h þ . We can reduce the number of solutions to solve for by accounting for the symmetries that exist in this scenario. First, if a solution is found for when an e -is the initial propagating charge, the solution for when a h þ is the initial propagating charge can be immediately found by swapping the τ and T parameters. Therefore we need only to keep the distinction between charges that are the same as the initial charge and charges that are the opposite. This is done by replacing the "e" and "h" labels in the τ and T parameters with "s" and "o" labels to indicate the same or opposite charge. Second, the solutions should be the same whether the charge is propagating toward the z ¼ 0 or z ¼ 1 surface. Therefore solutions need only be found for one direction of propagation. However for good practice, the solutions were solved for both directions of propagation and are confirmed to give matching results. Here we will solve some of the solutions for when the initial charge is propagating toward z ¼ 1.
We start by defining the probability distribution of the initial starting position of the charges. For surface charges that propagate toward z ¼ 1, the probability distribution of the initial starting position P surf ðz 0 Þ trivially goes as
Next, we can begin to solve for the solutions corresponding to specific scenarios. These probability distributions are indexed as P k ðEÞ, where k iterates through the different solutions. The first, and easiest, solution to solve for is the case where the initial charge reaches the surface
IMPROVED MODELING OF DETECTOR RESPONSE EFFECTS IN … PHYS. REV. D 109, 112018 (2024) 112018-13
without a CT or II process occurring. This probability P 0 ðEÞ is by combining Eqs. (A5) and (A7),
The solution to P 0 ðEÞ is a delta function at E ¼ 1 with an amplitude of e -T s . This makes sense, as a charge traveling from z ¼ 0 to z ¼ 1 will produce exactly one e -h þ pair worth of energy. Next, we can solve for the solution when the initial charge undergoes CT, P 1 ðEÞ, which is found by combining Eqs. (A6) and (A7),
The next scenario to consider is the case where the initial charge creates a new, like charge, and both the original charge and the new charge happen to reach the surface. We call this probability distribution P 2 ðEÞ. Let z II be the position where the new charge is created. It is important to remember that the quantity we are interested in is the total distance traveled by all the charges in the scenario. It is helpful to think about the energy gained by each "segment" of charge propagation. In this scenario, there are three such segments: the initial charge that travels from z 0 to z II , and the initial and additional charges that each travel from z II to z ¼ 1. It is additionally helpful to think of E as fixed number that constrains the problem. We want to find the combinations of segments that result in E total energy.
Let E 2a and E 2b be the energy contributions from each of the two charges after II. The combined energy contribution from both charges after II is therefore
is the energy contribution of the initial charge before II, then the total energy measured will be E ¼ E 1 þ E 2 . We also know that E 1 can be expressed as z IIz 0 . Rearranging gives z II ¼ E -E 2 þ z 0 . What we have done is expressed z II not as a position in the crystal, but rather in terms of energy contributions. The choice of expressing the parameters this way initially seems odd. After all, we already know that E 2a and E 2b are equal in this scenario. But that is not true for all scenarios, and it turns out that this way of formulating the problem provides a generic framework for solving all of the solutions.
We first need to find the probability that the energy contribution after II is E 2 . This probability is equal to the probability that one charge contributes an energy of E 2a times the probability that the other charge contributes an energy of E 2 -E 2a , given a starting position of z II and summed over all possibilities of E 2a . Expressed in terms of the base equations from Eq. (A5) gives
Next, the probability of having a total energy of E for a given starting position z 0 is the probability that the energy contribution after II is E 2 and the energy contribution before II is E -E 2 . We find this by combining Eqs. (A6) and (A10), and summing over all possibilities of E 2 ,
The final step to find P 2 ðEÞ is to integrate over all z 0
The next scenario to consider is when the initial charge creates an opposite charge, and both the original and new charge happen to reach the surface. We call this probability distribution P 3 ðEÞ. Like the previous scenario, the additional charge is created at a position z II , and we need to find the probabilities of each segment of charge propagation. Because of how we formulated the problem, the way to solve for P 3 ðEÞ is exactly the same as how to solve for P 2 ðEÞ with just two key substitutions. The II process created an opposite charge, and that opposite charge will propagate toward the z ¼ 0 surface. Therefore, one of the P 1 S;s terms in Eq. (A10) is replaced with P 0 S;o with the same inputs. And because the process in this scenario is II to an opposite charge, the P 1
IIss term in Eq. (A11) is replaced with P 1
IIso , also with the same inputs. Making these substitutions and solving for P 3 ðEÞ gives
The solutions for other scenarios can be found by employing the same logic of considering segments of energy contribution, using the correct combination of base equations, and nested integrals. In total, there are 28 solutions found for surface charge events, all of which are cataloged in the Supplemental Material [13]. The individual solutions are shown together in the top plot of Fig. 1.
A singe-charge bulk event starts with the creation of either an e -or h þ at some starting position z 0 that ranges between z ¼ 0 and z ¼ 1. As with surface charges, we can use the same symmetry arguments to reduce the number of solutions to solve for. Again, the e and h labels in the subscripts are replaced with s and o to indicate charges that are the same and opposite as the initial charge, and solutions are only needed to be found for charges propagating in one direction. For bulk-single-charge events, the probability distribution of z 0 is defined as a uniform distribution between z ¼ 0 and z ¼ 1:
We can again consider the simplest scenarios to solve for P 0 ðEÞ (the charge reaches the surface), P 1 ðEÞ (the charge undergoes CT), P 2 ðEÞ (the charge undergoes II to the same charge and both charges reach the surface), and P 3 ðEÞ (the charge undergoes II to the opposite charge and both charges reach the surface). Fortunately, these solutions are mostly solved for in Eqs. (A8)-(A13), except now P surf ðz 0 Þ is replaced with P bulk ðz 0 Þ. Making this substitution and solving for the probability distributions gives P 0 ðEÞ ¼ e -T s •E ; 0 ≤ E < 1; 0 else;
IIss e -T s •E ð2 -EÞ; 1 ≤ E < 2; 0 else;
It is evident that the solutions to these problems become complex. One way to determine if these solutions make sense is to examine the boundary conditions. For example, the probability distribution P 3 ðEÞ in Eq. (A15) ranges from one to two e -h þ pairs of energy. If a charge has an initial position of z ¼ 0 and immediately creates an opposite charge, the event will produce one e -h þ pair of energy. The same is true if the charge has an initial position of z ¼ 1 and immediately creates an opposite charge. There is no scenario where this process can produce an energy less than one e -h þ pair. Furthermore, if the charge has an initial position of z ¼ 0 and creates an opposite charge only when it reaches z ¼ 1, the event will produce two e -h þ pairs of energy. There is likewise no scenario where this process can produce an energy greater than two e -h þ pairs. The 28 unique solutions found for bulk-single-charge events are cataloged in the Supplemental Material [13], and the individual solutions are shown together in the middle plot of Fig. 1.
A bulk-e -h þ -pair event starts with the creation of both an electron and hole at some starting position z 0 that ranges between z ¼ 0 and z ¼ 1. As with bulk-single-charge events, we assume that z 0 is a uniform distribution between the surfaces of the detector and follows Eq. (A14). However unlike the solutions for single charges, the solutions for e -h þ -pair events need to keep the distinction between the parameters for electrons and holes. The initial e -and h þ will propagate in opposite directions and are treated as independent charges. The only constraint the initial starting position that they both share. Like with the single-charge events, solutions will be the same regardless of which direction of propagation is chosen for the charges.
Here we will demonstrate how to find the solutions for the simplest scenarios. Let P 0 ðEÞ be the probability that both the electron and hole reach the surface. We assume the electrons and holes travel toward the z ¼ 1 and z ¼ 0 surfaces, respectively. As before, the solution can be found by considering the segments of charge propagation in the scenario. If the total energy of the event is E and the electron contributes an energy of E 1 , then the hole must contribute an energy of E -E 1 . The probability of measuring an energy of E giving a starting position of z 0 is therefore the probability that the electron contributed an energy of E 1 starting at z 0 times the probability that the hole contributed an energy of E -E 1 starting at z 0 summed over all possibilities of E 1 . Using the base equations from Eq. (A5), this is written as
The final step to solve for P 0 ðEÞ is to multiply Eq. (A16) by Eq. (A14) and integrate over all z 0 . However this last step must be considered separately for when T e ¼ T h and T e ≠ T h in order to avoid undefined solutions. For the case where T e ≠ T h , P 0 ðEÞ is found to be
For the case where T e ¼ T h ≡ T, P 0 ðEÞ is found to be
The next scenario to consider is when the e -is trapped while the h þ reaches the surface. Let the probability distribution for this process be P 1 ðEÞ. As before, E 1 is the energy contribution from the electron, and E -E 1 is the energy contribution from the hole. The probability for measuring an energy E given a starting position of z 0 is found in the same way as in Eq. (A16), except that for the electron, the appropriate base equation from Eq. (A6) is used:
Again we can find P 1 ðEÞ by multiplying Eq. (A19) with Eq. (A14) and integrating over z 0 . For the case where T e ≠ T h , P 1 ðEÞ is found to be
For the case where T e ¼ T h ≡ T, P 1 ðEÞ is found to be
The solutions for the other scenarios can be found by considering the probabilities for the process that happens to each charge and constraining the energy contribution from each charge to the total measured energy. In total there are 16 unique solutions found for bulk-e -h þ -pair events, which are cataloged in the Supplemental Material [13].
The individual solutions are shown together in the bottom plot of Fig. 1.
Section III B introduced the phenomenon of nonionizing energy deposition and how the surface-trapping effect can be incorporated into the extended detector response model. In this model, the number of photons that hit the detector is given by a Poisson distribution with a mean of λ. Each photon will create an e -h þ pair, where there is a probability α that the e -h þ pair undergoes surface trapping. Absorbed photons that result in an e -h þ pair that undergoes surface trapping are classified as nonionizing photons, whereby the deposited energy in the detector will only be the absorption energy of the photon E γ . This effect results in the formation of a subpeak structure at each e -h þ -pair peak in the energy spectrum, as can be seen in the middle plot of Fig. 3. This appendix provides further details on the distribution of these subpeak structures and its dependence on λ and α.
Each subpeak corresponds to q ionizing photons and p nonionizing photons. The q ionizing photons will produce q e -h þ pairs that will propagate through the detector where they may undergo bulk CT and II processes. As will be discussed below, the bulk CT and II processes do not affect the shape of the underlying subpeak distributions, and thus can be ignored. For the subpeak distribution of p at the qth e -h þ -pair peak, we want to determine the probability PðpjqÞ, simply given as
To find these probabilities, we must first consider the probabilities of the separate processes. The total number of photons absorbed in the detector is (p þ q), and the probability of (p þ q) photons hitting the detector is determined from a Poisson distribution with a mean of λ. The probability of having p nonionizing photons is determined from a binomial distribution with (p þ q) trials and a probability of α. Therefore Pðq ∩ pÞ is the probability that (p þ q) photons hit the detector and p photons are nonionizing:
If the mean number of photons hitting the detector is λ and there is a (1 -α) probability that a photon will be ionizing, then the mean number of ionizing photons hitting the detector is λ • ð1 -αÞ. Therefore PðqÞ is determined from a Poisson distribution with a mean of λ • ð1 -αÞ:
Likewise, PðpÞ is determined from a Poisson distribution with a mean of λ • α:
Using Eqs. (B2) and (B3), Eq. (B1) becomes
λ ðqþpÞ e -λ α p ð1 -αÞ q λ q ð1 -αÞ q e -λð1-αÞ ¼ λ p α p e -λ•α p! ¼ Poiss:ðp; λ • αÞ:
Equation (B5) shows that the subpeak distribution of p for a given q is just the probability of having p nonionizing photons, which is a Poisson distribution with a mean of λ • α. Importantly, the subpeak distribution is independent of q, and is therefore the same for each e -h þ -pair peak. As shown in Fig. 3, when the resolution smearing is applied, these subpeak distributions shift the location of the e -h þpair peaks in the spectrum. The amount that the peaks are shifted by ΔE ph is determined from mean energy of nonionizing photons at each e -h þ -pair peak. As the subpeak distributions are the same for each peak, the amount that each peak is shifted by is also constant:
Lastly, we consider what effect the bulk CT and II processes may have on the subpeak distributions. For the qth e -h þ -pair peak there are q ionizing photons and thus q e -h þ pairs that propagate through the detector. The peaks in the subpeak distribution only arise when all of the primary charges from the e -h þ pairs reach the surface without undergoing a CT or II process, otherwise the measured energy will be in a nonquantized region of the spectrum. Appendix A showed that the probability of a charge from a surface event to traverse the detector without undergoing a CT or II process is e -T e=h , where T e=h encodes the CT and II probabilities for either the electron or hole. The probability of q charges from a surface event to traverse the detector without undergoing CT or II processes is found from a binomial distribution with q trials and a probability of e -T e=h : Binom:ðq; q; e -T e=h Þ ¼ e -qT e=h .
Therefore while the overall scaling of the subpeak bution depends on q, the shape of the underlying distribution remains for each e -h þ -pair peak.
As discussed in Sec. V, the exponential CTII model is limited by the highest order of processes that are modeled. Specifically, solutions for surface events and bulk-singlecharge events are found for up to second-order processes, whereas for bulk-e -h þ -pair events, solutions are found for up to first-order processes. In order to assess and quantify these limitations, the single-e -h þ -pair solutions are compared to simple MC simulations of the CT and II processes. The MC simulations model the CT and II processes using the same initial assumptions as the analytical model: that the probability distributions of CT or II occurring are described by Eq. ( 5), and where CT and II processes are parametrized by the characteristic lengths τ i . However unlike the analytical model, the MC simulations are able to include higher-order CT and II processes. In these MC simulations, there are no physical or detector response processes that are modeled other than CT, II, and generic resolution smearing.
We would like to determine where these higher-order processes become significant such that the analytical model is no longer a suitable description of the MC simulations, and thus of these CT and II processes. There are two main factors that will cause the analytical solutions to deviate from the MC simulations. The first is the total probability of impact ionization, and the second is the total number of events in the simulation. Increasing either the total probability of II or the total number of events will increase the number of events in the MC simulation that undergo higher-order CT or II process that the analytical solutions do not model.
The limitations of the single-e -h þ -pair solutions can then evaluated by using a simple procedure. For each event type, we scanned over the total II probability in the model and the total number of events in the MC simulation. After computing the analytical model and running the simulation for each set of parameters, we performed a Kolmogorov-Smirnov (KS) test to determine if the simulated spectrum is described by the analytical model for a given level of confidence. For simplicity, we define the total II probability of a single charge f II;tot as f II;tot ¼ f IIee þ f IIeh ¼ f IIhe þ f IIhh , where each f i is equal to f II;tot =2. For surface events and bulk-single-charge events, the tests assume that the initial charge is an electron. Furthermore, each of the solutions and simulations assume a small CT probability of f CTe ¼ f CTh ¼ 1% in order to include all possible processes. The KS tests take the null hypothesis that the MC simulation is described by the same probability distribution as the analytical model, and the results from the tests are subsequently placed into three categories: accepted (failed to reject the null hypothesis at 90% confidence level), rejected the null hypothesis at a 90% confidence level, and rejected the null hypothesis at a 99% confidence level. The results of the KS tests for each event type are shown in Fig. 7.
The test results from Fig. 7 clearly illustrate the regions of this parameter space where the analytical solutions of the exponential CTII model deviate from the MC simulations. For reference, measurements of f II;tot in HVeV detectors have been on the order of 1% [8,20]. However, these results represent a worst-case scenario for the model, as other parameters can extend the boundary of this limitation. For instance, high CT probabilities will generally lower the probabilities of high-order II processes. Furthermore, these tests were performed using the single-e -h þ -pair solutions, FIG. 7. Results of the KS tests performed using the analytical CT and II solutions and the MC simulations, where each dot corresponds to a test performed for a particular value of the total II probability for a single charge f II;tot and of the total number of events in the simulation. The tests were performed separately for surface events (top), bulk-single-charge events (middle), and bulk-e -h þ -pair events (bottom).
whereas modeling the full detector response will often require the multi-e -h þ solutions. In many cases, the multie -h þ solutions will cause the high-order II processes to be subdominant within the total probability distribution, as can be seen by comparing the top and bottom plots of Fig. 2. In these scenarios, the analytical solutions may adequately describe the CT and II processes even for higher II probabilities or for a larger number of events.
112018-12