NEST: A Comprehensive Model for Scintillation Yield in Liquid Xenon
A comprehensive model for explaining scintillation yield in liquid xenon is introduced. We unify various definitions of work function which abound in the literature and incorporate all available data on electron recoil scintillation yield. This results in a better understanding of electron recoil, and facilitates an improved description of nuclear recoil. An incident gamma energy range of (1 keV) to (1 MeV) and electric fields between 0 and (10 kV/cm) are incorporated into this heuristic model. We show results from a Geant4 implementation, but because the model has a few free parameters, implementation in any simulation package should be simple. We use a quasi-empirical approach, with an objective of improving detector calibrations and performance verification. The model will aid in the design and optimization of future detectors. This model is also easy to extend to other noble elements. In this paper we lay the foundation for an exhaustive simulation code which we call NEST (Noble Element Simulation Technique).
1 Introduction and Motivations
Liquid noble elements (LNE) have been established as an attractive detection medium in experiments searching for dark matter or neutrino-less double-beta decay: a large number of existing or planned experiments employ liquid xenon, argon, or neon [1, 2]. The direct dark matter search experiments looking for WIMPs (Weakly Interactive Massive Particles) aim to maximize their ability to discriminate between electron and nuclear recoils (ER and NR respectively) within their target volume. While the advantages of LNEs are abundant, one drawback is the non-linear energy dependence of the scintillation yield per unit energy deposited in the medium, for both ER and NR [3, 4]. Further, this phenomenon is influenced by the magnitude of the applied electric field . In this article, we have addressed these issues for ER in liquid xenon by developing a general model for scintillation and ionization yields. We have implemented this model in a Geant4  simulation and have demonstrated its ability to reproduce a wide variety of measurements, both for the field-free case and for high applied field. We have also extended our model to include NR scintillation yield in the zero-field case.
An understanding of light yield is especially important in the low-energy region of interest (5-50 keV) which contains the majority of the spectrum for NR from WIMPs in xenon ( and ref. therein). The energy dependence of light yield is non-monotonic at low energies, stemming from effects related to , incident particle type, and electric field. This poses a challenge for simulations that employ constant scintillation yields. This model provides a comprehensive Monte Carlo similar to that of Tawara et al.  for NaI(Tl). Although this paper focuses solely on liquid xenon, in Section 7 we point to simple steps needed to include other noble elements.
2 Physical and Mathematical Framework
We begin by reviewing the basic principles of scintillation physics. A particle depositing energy creates excitation and ionization, as well as heat. Initial excitation contributes to the production of light called the primary scintillation, or S1, through an intricate process involving excited molecular states (unnecessary to fully simulate to reproduce the light yield) [7, 8, 9]. The partitioning of deposited energy into excitation and ionization depends on incident particle type and possibly on electric field , but previous work indicates that the partition is independent of incident particle energy [10, 11]. This entire process is summarized in Figure 1. Some energy is lost to heat in the form of secondary nuclear recoils [12, 13] or sub-excitation electrons, which fail to lead to the generation of S1 light . For ER, the former is negligible and the effect of the latter can be absorbed into a higher value for the work function (below).
Ionization electrons either recombine and lead to additional (S1) scintillation or escape. Two-phase LNE direct dark matter detectors  drift such electrons with an electric field and extract them into a gas volume where they drift at higher velocity and generate more scintillation through more excitation and recombination in the gas. This is termed secondary, or S2, scintillation. The ratio of S2 to S1 light differs for ER and NR ([2, 7] and ref. therein), and hence, it is the most critical discriminating parameter in dark matter searches. The first step in fully understanding this phenomenon is to have a comprehensive description of ER behavior in xenon.
The distribution of the ionization electron population between recombination and escape constitutes the crux of the energy dependence of the S1 yield. The direct excitation contribution to S1, which depends on the small ratio of excitons to ions, is consequently small ( and ref. therein). Since an ionization electron lost to recombination cannot escape and one which escapes cannot recombine, S1 and S2 light yields are anti-correlated. Furthermore, scintillation yield, and hence recombination probability, is known to be a function of the energy loss per unit of path length  and of the applied electric field .
where is energy deposited by a particle in a single interaction (not its total energy deposition), is the number of excitons created per deposition, the number of ions, the work function (required energy) for exciting atoms, and the ionization work function. One sums the numbers of excitons and ions stemming from single interactions (recoils) to calculate the total yield of one incident particle. The ratio of excitons to ions, often labeled as , may differ for NR versus ER, and may change with field, but it is likely not a function of energy . The theoretical value of this dimensionless ratio is a constant 0.06 for liquid xenon [14, 15], used successfully to fit experimental data . Some measurements of in liquid xenon have been as high as 0.20, though arguably consistent with 0.06 within experimental uncertainties [4, 17].
Ignoring heat, including the kinetic energy of ionization electrons (which can be absorbed into ), and assuming near-100% efficiency for excited or recombined electrons to lead to S1 (well-established experimentally ), the numbers of quanta produced are
which in turn lead to photons and ionization electrons:
derived from Birks’ Law . Doke et al.  derived a scintillation yield for LNEs from Birks’ Law, while we adopt the same approach but instead derive the recombination probability underlying the scintillation yield by separating excitation from the recombination probability. The first term of the equation represents the recombination which occurs when a wandering ionization electron is captured by an ion other than its parent (termed “volume” recombination). The second term () represents “geminate” recombination, also known as Onsager recombination [18, 20]. It is defined as ionization electrons recombining with their parent ions, which is assigned a fixed probability for all . A low implies that most S1 light is a result of the recombination processes.
We plot Equation 6 in Figure 2 with the parameter values that we will show in Section 3 to reproduce the light yield of gamma rays, the first particles we consider in this work. The constraint is imposed so that as approaches infinity, the recombination probability approaches unity.
where is a constant dependent on ionization electron and hole mobilities and the dielectric constant (beyond the scope of this work), is the mean ionization electron velocity, and is a length scale defining ionization density volume. We can equate with , where is the electron charge, is the dielectric constant, is Boltzmann’s constant, and is the electron temperature. Using , K , and m , we estimate a value of 0.14 for zero electric field. We assume here that most electrons will thermalize prior to recombination (and for use the 1 atm boiling point, near which most detectors operate ). This is false at high electric field, which increases electron energy .
We define “short tracks” to be those that are shorter than the mean ionization electron-ion thermalization distance, 4.6 m in liquid xenon . If an incident gamma or electron produces secondary gammas or scatters multiple times, then there are multiple interaction sites. At each site, tracks are separately summed to determine whether to apply Thomas-Imel or Doke/Birks formalism. Each is based on the same recombination physics, albeit with different descriptions governed by short- or long-track limits . The Thomas-Imel understanding relies on recombination in a “box” (cubic) geometry [10, 20] while Doke constructs his formulation utilizing the cross-section of a long column of electron-ion pairs . We seamlessly integrate these two models into one unified picture of scintillation yield, which is continuous versus energy and consistent with data.
Our model determines recombination for individual recoils in a Geant4 simulation  by explicitly tracking energies and recoil ranges for secondary electrons produced by an incident gamma. Doing so preserves macroscopic characteristics of the energy-dependent light yield (Figure 5) while remaining stochastic. As demonstrated in Figure 3, Monte Carlo fluctuations in tracks result in a variation in the number of scintillation photons ultimately produced, even before fluctuation in the ratio of excitons to ions, or their sum (Fano factor), is taken into consideration.
For simplicity, our model defines a mean work function for production of either excitons or electron-ion pairs, such that . This is equally as effective as defining two distinct work functions for successfully explaining experimental results [4, 10, 27]. Different values for the two processes are difficult to extract experimentally, and using a mean value causes no loss of generality [10, 25]. Because is also equivalent to based on Equations 1 to 3, a higher is canceled out by a lower or vice versa, leaving the total numbers of both quanta unchanged. Similarly, a higher is counteracted in the above equation by a lower , making an exact experimental knowledge of less relevant.
In general, increases with decreasing energy for electrons less than 1 MeV . Hence, recombination probability (and scintillation yield) should increase with decreasing energy, as long as the Doke model applies [4, 18]. However, an interesting phenomenon occurs in a low-energy region where the scintillation yield instead decreases (Figure 4). For gamma rays below (10 keV), the slope of the light yield curve changes sign , as it does for NaI(Tl) in a different energy regime [30, 31], albeit due to a depletion of activation sites in the crystal .
One would expect recombination to be enhanced at high due to greater ionization density, but it has been shown that at low energies (high ) recombination becomes independent of ionization density [10, 13, 20]. This is accommodated in our model because low-energy particles are strictly in the Thomas-Imel regime, where the recombination probability depends on the energy, via the number of ions, and not on . In other words, for particle tracks smaller than the thermalization distance of ionization electrons, the total track length is no longer relevant. Dahl was the first to apply a modified Thomas-Imel box model to successfully explain the odd turnover in light yield at low energy, which changes position in energy for different electric fields . At higher energies, the Thomas-Imel box model had already been successfully utilized to fit data on the global scale, such as total charge yield (total number of ionization electrons extracted) versus electric field, glossing over microphysics . In Dahl’s approach, the model is applied to the scale of the lowest energy (sub-keV and eV) individual electron recoils using PENELOPE .
In our model, we apply the Thomas-Imel model to the electron-recoil energy scale (1 keV) accessible with Geant4, making the assumption that it can be adapted for zero electric field. Manalaysay et al.  demonstrate that the Thomas-Imel parameterization, even though undefined at zero field, can have a finite limit and produce realistic results. Moreover, Figure 5 in this work demonstrates the utility of this assumption by successfully reproducing light yields.
3 Zero-Field Scintillation Yield vs. Energy
The parameters of Equation 6 and of 7 (henceforth ) were treated as free in order to match our compilation of existing data on zero-field gamma ray scintillation. A concerted effort was made to include all known experimental results in our compilation and any possible omitted works are unknown to the authors rather than intentionally excluded. Figure 5 compares the results from simulations for this work to the available empirical data.
A traditional analysis with errors could not be performed because many papers do not quote error bars in either their figures or their text (e.g., ), or state that their errors are negligible . Furthermore, there are several instances of more than one data point for a given gamma energy, with one or more points lacking errors, and yields are most often quoted as relative to measurements made at 122, 511, or 662 keV. The latter is due to experimental issues, such as difficulties in knowing PMT quantum efficiency at low temperatures and UV wavelengths and/or inadequate knowledge of the reflectivity of the walls of the detector. Thus, assumptions had to be made in order to translate all results into absolute yields. We used 511 and 122 keV as benchmarks, translating uses of 662 keV or other energies by experiments as benchmarks as needed. For the low-energy data of Obodovskii and Ospanov , the 59.5 keV Am line is used.
The process of translating the data into absolute yield increases the error by an amount which can not accurately be estimated due to the varying age and quality of the data, and in some rare cases, due to variation in choice of benchmark energy, which has a profound effect. For instance, if the results of Barabanov et al.  are re-scaled assuming a benchmark of 122 keV in Figure 5, then 59.5 and 122 keV are not outliers, but all other points become too low. Choices like this needed to be made to bring the majority of a large quantity of old and new data into rough agreement. Overall, the ultimate consistency among different detectors is remarkable. Precedence for such comparison of different data sets is exemplified by Doke et al. , who were forced to assume a match in yield between one Na gamma line and a high-energy electron in order to incorporate the results of Barabanov et al. into their Figure 4. We treat data all as equally as possible in Figure 5 without averaging results except when necessary, such as when we calculated a mean 122 keV yield in order to incorporate data which referenced it but not 511 keV (which was our primary benchmark because it was as common as 122 keV in experiments, but provided the additional advantage of forcing the majority of data from Barabanov et al. to agree with other works).
|data window and other conditions||parameters for best reproduction (this work)|
|low energy points alone (15 keV)|
|high energy points treated independently|
|averaging points (raw) at same energy|
We minimize the mean-squared of the residuals in order to optimize the parameters of our model, which is founded on Equations 6 and 7. We perform a series of full Monte Carlo simulations with different sets of parameters until finding the optimal set. Then our model can successfully reproduce the available data, thus predicting yield in regions where the data is sparse. The two clear outliers from Barabanov et al.  are excluded. They may be explained by non-uniform light collection in the detector, as claimed by . Also discarded, for the purpose of our plot and the optimization, was one result considerably distant from all others, quoting less than 34 photons/keV at 22 and 88 keV . Below 15 keV, the Thomas-Imel box model dominates, and it is in this regime where we find the optimal value for . We then use that value for the data as a whole, keeping it fixed. The top of Figure 5 is a plot of the percent residuals between various data sets and our model. Given our assumptions for relative-absolute translation, most data do not differ by more than 4% from a simulation with the parameters which best reproduce the data. These are listed in Table 1, above. Throughout this work, we do not fit data on scintillation yield directly. Instead we systematically explore ansatz solutions for the recombination probability which are then entered into the simulation to reproduce some measured yields and predict others.
In addition to the list in Table 1, other parameters which we employed, but that we held fixed, included eV, , and 4.6 m as the electron-ion thermalization length scale, which we pragmatically selected to determine the cross-over energy where the Doke model would begin to dominate, because of how well it predicts the approximate location of the major change in slope in scintillation yield as a function of energy. The -value is an error-weighted average from reports by Doke et al.  and Dahl : 13.8 0.9 and 13.7 0.2 eV. At least three other values exist, but were not included in the average: 14.0 eV because it was reported without a clearly stated uncertainty , 13.46 0.29 eV  due to a calibration issue expressed by Dahl  (but it was still used to estimate a 122 keV yield), and 14.7 eV [9, 15, 41], superseded in the past thirty-five years by more measurements benefiting from technological improvements.
The y-axis labeled “effective” work function
The clear dip in light yield at 30-35 keV is present because of a resonance at the xenon K-edge for production of low-energy Auger electrons [36, 43, 44, 45]. They are low enough in energy ( keV) to be governed by the Thomas-Imel model. Furthermore, at or above about 29.8 keV, the K-shell x-ray is emitted, traveling a small distance to create a second interaction site for which is separately calculated, serving to enhance the effect of the scintillation decrease: tracks are summed separately at each site to determine which recombination treatment to use, enhancing the probability of being in the Thomas-Imel regime. We note that the simulated dip is too shallow, possibly due to an inability of Geant4 to produce electrons of low enough energy to accurately depict the necessary interactions, even with the Auger physics module activated. In Geant4, a gamma does not generate electrons whose energy is so low that their range would be below the minimum distance cut-off. In order to maintain energy conservation, an energy deposition is recorded for the gamma in the tracking output file. To partially mitigate this issue, we assume that an ER with that recorded energy occurs .
Thus far we have only discussed gamma rays. We now treat electrons as the incident particles and subject them to a similar simulation as was done for incident gammas. We obtain eV, for incident electron energies in the range 1 - 100 keV. For electrons below (100 keV) in energy, Séguinot et al. report eV  in agreement (5%) with the lower edge of our predicted range. Given the energy dependence of , Séguinot et al.’s measurement must be an effective value, most likely an average or a minimum. In their earlier work , they claim a value of 12.5 or 12.7 eV, in disagreement with their later work and our present work, but Chepel et al.  cite Séguinot et al. as 12.7 1.3 eV, where the upper bound on the uncertainty (14.0 eV) disagrees by only 5% with our prediction (14.7 eV). We cannot, however, reconcile their claim of a low ionization work function (defined later in this section) for low-energy electrons of 9.76 eV [26, 46] using any of the possible definitions. Further, it has been disputed [42, 47, 48].
Unlike for low-energy electrons, data on high-energy ones are plentiful , derived from Bi decay, which is a mixture of 0.976 and 1.048 MeV electrons, and a 1.064 MeV gamma, plus other lower-energy decay products . Authors typically quote yield as that of an approximate “1 MeV electron,” implying that the electrons dominate; cited yields are all in agreement with each other within uncertainties [3, 4]. We take two typical empirical results to compare to our simulation: 42.0 8.4 photons/keV  and 46.4 photons/keV.
Our initial attempt to match the absolute yield of 1 MeV electrons using the Doke model was not successful. Forcing a match results in too low of an absolute yield for other well-established results. Setting the (relative) 122 keV point from Barabanov et al.  equal to one of the absolute measurements of the 122 keV yield lowers the yield of the 1 MeV gamma, hence that of the 1 MeV electron as well, creating better agreement with data. Unfortunately, this makes other Barabanov et al. points too low, in direct conflict with much more recent measurements [34, 36] and barely within the quoted error of Ni et al. . (Ni et al. report absolute, not relative, yields, requiring no assumptions in interpreting their results.) The problem at hand is subtly evident in the work of Doke et al. (the basis for the curve in Figure 5 labeled Doke 2001). They use data from Yamashita et al.  at the expense of ignoring Barabanov et al.  and determine the absolute yields of 1 MeV electrons and gammas to be nearly identical ; using 122 or 511 keV as the key to translate relative into absolute yield makes the absolute yield too high again for a 1 MeV electron as in our own work. In Ni et al. , one possible solution can be found: using a universal (average) for interactions (explained in Section 4), the Doke model can be made to fit the 1 MeV electron yield along with equal or lower-energy gammas, but we wish to avoid losing the stochastic nature of the simulation. Finally, even if one were to try applying an existing model other than Doke’s, the similarity in between a gamma and an electron at 1 MeV and the direct contradiction between one older set and two or three newer data sets would remain points of contention.
|conditions (parameters)||yield (ph/keV)|
|(unchanged model parameters)||57.4|
|(geminate/Onsager recombination alone)||43.1|
|(additional initial excitation)||46.7|
|most recent experimental data, with the smallest uncertainty ||46.4|
An extreme solution to the challenge is to assume that non-geminate recombination ceases for a minimally ionizing particle. Parameter of Equation 6 represents such recombination, and is proportional to ionization density, which is in turn proportional to . Perhaps low ionization density prevents ionization electrons from recombining with non-parent ions because of distance, though Mozumder claims that this is impossible . The thermalization distance (4.6 m) is much greater than the reach of the Coulomb field of a xenon ion as defined by the Onsager radius (49 nm) , so Onsager recombination should not be dominant. For the sake of argument, we surmise that a thermalized ionization electron can eventually be re-attracted to its parent ion if no other ion is available. Alternatively, while thermalizing, an ionization electron could become trapped in a closed path within the Coulomb field of the parent ion if repelled by neighboring electron clouds . This may help explain the slow recombination time observed by Doke and Masuda . However, we do not present our solution as a definitive one, given a lack of fundamental plausibility, but instead as a pragmatic one, allowing the creation of a simulation which produces reasonable results for gammas and electrons alike in a wide range of energies. A summary of our results for the 1 MeV electron under parameter variation is presented in Table 2. Changing to correspond with the results of  in addition to setting matches data best. A higher may imply more Onsager recombination than encoded in . More such recombination is hard to distinguish from a shift in the initial ratio of excitons to ions .
After reproducing zero-field scintillation yield it is possible to make predictions of measurements other than yield. Aprile et al. report that the zero-field escape fraction for ionization electrons from a 662 keV gamma is  (also called ). From our simulation we derive . With an of 0.20 and the adjustment of explained above, we can reproduce the reported for Bi decay . In addition to we can derive a value for the effective ionization work function, defined as . The definition of is equivalent [4, 9, 41].
4 Discussion of Other Models
In this section, we discuss the advantages of our model as compared to others that are available in the literature. Figure 5 includes the calculation of Doke et al. . This approach is based on extracting an electron response from a portion of the gamma data [29, 36], applying that extracted response to early-generation secondary electrons and then reconstructing the gamma response. In contrast, our model, which uses Equations 6 and 7, is applied to each new electron created in the cascade, as low in energy as permitted by Geant4. The range of each electron is obtained from Geant4 output, and not assumed from a Bethe-Bloch calculation. Doing so allows for variations in . Such variations contribute to the spread in scintillation yield (Figures 2 and 3). This approach can be incorporated in any simulation which outputs lists of energy depositions, and utilized for predicting charge yields (Section 5) and the NR quenching factor (Section 6). Unlike some other previous works [3, 4], we do not utilize a direct fit to the energy dependence of total scintillation yield. We thus predict, not extrapolate, the response to an energy regime that is lower than what has been studied experimentally.
In earlier work, similar models have been applied to materials other than LNE, for instance, NaI(Tl), which is a well known scintillator. The recombination probability () for NaI(Tl) has the same functional form  as Equation 6, albeit with , leaving . Despite the compositional differences between mono-atomic LNE and solid crystal scintillators (often comprised of multiple elements, including dopants), or organic scintillators for which Birks’ Law was derived originally, this formula has been applied to LNE with modest success [3, 4, 18]. However, its application involved defining a universal for a gamma interaction, which had to be obtained by averaging over all interactions (ERs). Electron numbers and energies of course vary , diminishing the utility of a mean by leading to a static recombination probability. Our treatment of xenon avoids usage of such an average.
There was, however, no certainty that Equation 6, which works relatively well with global , would also work when accounting for individual . We explored numerous other possibilities with the same qualitative features (low probability at low and high probability at high ), but these were unmotivated and made for poor reproductions of the data (, exponential, Equation 6 with raised to different powers, and non-unity at infinite .) No benefit was seen in inventing a new approach, and using (no geminate recombination) did not minimize the mean-squared error, while a non-zero did. Numerous other physically-motivated models exist, Mozumder’s for example , but have more free parameters, are inapplicable at zero field, with no simple adaptation, or apply only to MeV energies . We achieve agreement with the most data and the fewest parameters. There also appeared to be no benefit to the introduction of heat loss for ER, to explain the trend of low-energy gamma data. The of a low-energy ER is comparable to that of an NR , but the Lindhard mechanism, by which NR lose energy to additional nuclear recoils, does not apply to low-energy gammas, which should exclusively produce ER [10, 12].
Our model is ultimately an approximation because Geant4 does not accurately simulate electron production (tertiary, quaternary generations, etc.) below 250 eV , unlike other software packages with more precise low-energy electron physics, such as PENELOPE [10, 32]. However, we demonstrate that working within the default Geant4 energy regime is sufficient for all practical purposes. As shown in Figure 5 (top), our work can match data in most cases to better than 5%, assuming an accurate conversion from relative to absolute gamma and electron yield at zero field. In addition, the model continues to work well for non-zero field, as shown in the next section.
There are various extensions to the existing simulation packages which address scintillation directly. For example, RAT  is a Geant4 add-on that simulates the scintillation yield of many elements, with the exception of xenon. Our modeling of ionization and our focus on xenon, both lacking in RAT , make our work complementary.
5 Scintillation Quenching with Electric Field
As electric field strength increases, ionization electrons liberated by an interaction are increasingly less likely to recombine with an ion. This is known as electric field scintillation quenching ( and ref. therein), not to be confused with the quenching factor for NR, which involves different scintillation loss mechanisms such as the Lindhard effect [8, 12, 13]. For a two-phase detector, this results in an increase in S2 light output at the expense of S1. The problem is quite complex because recombination probability does not change uniformly across all energies . Dahl demonstrates that a fundamental theory of recombination as a function of electric field is possible, which would explain the trend in recombination probability using a modified Thomas-Imel box model . For reasons delineated below, we approach the problem semi-empirically, adapting as free parameters , , and from Equation 6, and from Equation 7. We tune these on a subset of available data sets, and then accurately predict many others.
An alternative would be to attempt a first-principles approach in Geant4 by creating look-up tables for the recombination probability based on a Thomas-Imel application in PENELOPE. For three reasons, we decided against this approach. First, the introduction of tables would eliminate the benefits of the stochastic nature of simulations, an attribute that was a primary motivation for implementing our model. Even a look-up table with stochastic variations fails to capitalize on the details of track history captured by Geant4, thus making it less realistic. Second, the first-principles approach breaks down at zero field, where single-phase detectors operate. Third, it is not verified experimentally above (100 keV) as opposed to the approach employed by Doke et al. [3, 4, 18]. Energies in the 0.5 - 3.0 MeV range are important for double-beta decay experiments and use of xenon in PET scans . Thus, our hybrid model of Thomas-Imel plus Doke covers a large range in both electric field and energy.
|Electric Field (V/cm)||This Work||Dahl 2009|
We began our data-matching with one data set that is the most comprehensive in terms of simultaneous span in energy (2-200 keV) and electric field (60-4060 V/cm) . By starting with one very complete data set, we initially sidestepped the issue related to collecting data sets from different detectors with distinct systematic effects. First, we let vary, using data below 15 keV to determine the best reproductions of data. Our results, in Table 3, differ from the values used by Dahl  for the same data; despite the disagreement the data is reproduced well, as demonstrated later in Figure 8. Agreement with data in spite of disagreement in the value of a parameter from the same model is likely a consequence of not applying the Thomas-Imel formula to the lowest energy ( eV) electrons. Geant4 does not generate them, and we are forced to utilize the Thomas-Imel box model in an approximate fashion. (We adapted Thomas-Imel to zero field at low energy without an extrapolation, but with a recalculation, in Section 3). Another contributing factor is the fact that Dahl applies the Thomas-Imel model completely to all his data, while we use it only at low energy, switching to Doke/Birks for the longer tracks of particles at higher energies. The lack of tracking of the shortest tracks causes the Thomas-Imel model to break down in Geant4 at high energy, but we take advantage of the regime where it does work. After having found values for at five different fields, we fit a power law to it: (field in V/cm). This is consistent with Dahl’s expectation of .
The next step was to find the Doke parameters as a function of electric field magnitude. We began this study by looking at the well-established electric field dependence of the 122 keV gamma. The result for 122 keV has been repeated several times [10, 22, 54], with duplicate and triplicate points often taken at the same electric field to verify repeatability . We proceeded to check the variation of the Doke parameters achieved studying only this one energy against other energies and to adjust our fit for the electric field dependence of the parameters accordingly, though continuing to give 122 keV results the greatest weight.
There must be a transition between the Thomas-Imel regime and the Doke regime. At zero field we use the ionization electron thermalization length scale as the cross-over distance, and this immediately makes the transition smooth. As electric field increases a need arises to change this cross-over distance in order to avoid creating a discontinuity. We offer a power law fit as a function of electric field as a first effort without optimization: the distance in m is , where the electric field magnitude is in units of V/cm.
Monotonic variation emerged naturally in our model (Figure 6) for and , and a negligible value for , so we set . These parameters have physical meaning as discussed earlier. We interpret a near-zero value for as demonstrating that volume recombination dominates over Onsager recombination. This corresponds exactly with the findings of Dahl, who reports that Onsager recombination should be negligible . With an external electric field overwhelming the local fields of ions, this is expected. In the case of zero field, the opposite appears true (), but this is empirically supported in argon , with available evidence being ambiguous for xenon [3, 4]. This is one approach, but we were unable to explain as much data as easily with a fully physically-motivated model for electric field, the Jaffé model, , where is the so-called recombination coefficient and is the electric field. However, the Jaffé model came close to matching significant amounts of data as modified by Aprile et al. to take into account delta electrons [7, 49]. Our approach to the electric field dependence has the advantage of adapting the existing formulae for zero field instead of introducing more models.
Our total compilation of comparisons is available in Figure 7. The data points are shown without any error bars for either the yield or the electric field magnitude, for clarity. (There are few cases where they are even available.) A significant amount of recent data is predicted accurately, having begun by reproducing 122 keV points and the low-energy regime in Dahl’s data. Discrepancies exist, but were inevitable, due to contradictory results being present in such a large compilation (likely due to non-obvious and unreported systematic errors), and due to the fact that we never directly fit any of these experimental yield-dependence curves, even at 122 keV. Instead, we reproduce or predict data with simulation by testing differing underlying recombination probability curves. With a few parameters that describe the recombination process for individual recoils and are functions of electric field, we can effectively explain global aspects such as total light and charge yields. These yields combine individual yields from all recoils spawned by one incident particle. This is a clear break from the past work of others who directly fit the total yields with models (for example, [49, 16, 56]).
We omit works ( for example) which define “total” light and charge yields for partial energy deposition “peaks,” as these would require proper modeling of the detector in question in order to reproduce well. We also exclude older works with measurements that cannot be reconciled with more recent ones, including , and Voronova et al. as cited in . More recent experiments have benefited from technological advancements, leading to more accurate measurements. Recent 9.4 keV data is poorly predicted by the Thomas-Imel model. Curiously, agreements can be achieved by approximating the Kr de-excitation as only a gamma ray and applying the Doke model despite the low energy. Without such adjustments these data contradict Dahl  in this energy range. The contradiction may arise from inaccurate decomposition into constituent electrons and gammas in Geant4.
In spite of these few faults, we have successfully reproduced many measurements based on various detectors. Figure 8 shows our prediction for as superposed on the data from Dahl . While the low-energy regime is a reproduction based on these data, the extension above 15 keV is the prediction from our model. Here we have compared simulation with data at only one of Dahl’s electric field choices, as an example. Comparable agreement is seen at other values of electric field as well.
6 Nuclear Recoil Quenching Factor
We have applied the framework motivated in Section 2, and validated in Section 3, to zero-field NR scintillation yield. NR scintillation efficiency has been a subject of debate, with several contradictory experimental results and theoretical underpinnings [58, 59]. This is an issue of great significance for dark matter detectors. Dark matter particles are expected to produce nuclear, not electron, recoils in most theories .
We offer a unique approach here (applied earlier only by Dahl , but at non-zero field) that uses parameters from a model motivated by gamma-induced ER and applies them to NR directly. Tracks of NR below (100 keV) are always within the Thomas-Imel regime . We again use the single free parameter from the Thomas-Imel model () that successfully explained the low-energy yield of ER. We make only one modification to the mathematical framework for NR, such that Equation 1 is modified [10, 13] as follows:
where is the Lindhard factor (or, effective Lindhard factor for modification to the Lindhard theory made to better explain data). It should not be confused with , which represents the ratio of scintillation light produced by an NR to that produced by a 122 keV gamma, at zero electric field. Both are functions of energy, and the latter is the traditional way to report results on NR. Though not equivalent to , knowledge of enables an empirical determination of , which represents the fraction of energy available for excitation and ionization. Unlike ER, NR will lose most of its energy (80%) to heat, which in this case implies interactions with other nuclei instead of electrons. Nuclear interactions with electrons excite or ionize atoms and thus produce usable signals for a detector [7, 13]. The fractional factor is the reason why two energy scales are often defined: keV and keV (or keV) [2, 10]. The former is for ER and the latter for NR. In this paper we use only “keV” as the total energy of a recoil as reported by Geant4. Without knowing or one can not a priori define the NR energy scale in terms of the “electron equivalent” energy in “keV.”
We use the same value of () as for ER to evaluate the model against empirical data. A choice remains in the application of the Lindhard factor with respect to proposed corrections. We elect to investigate two scenarios in order to treat this intricate subject with due diligence: pure Lindhard, and the Hitachi correction to the Lindhard factor . Lindhard theory is more appropriate for solid crystal scintillators, and for xenon nuclei Hitachi has performed a first-principles recalculation . We thus closely follow the methodology of Sorensen and Dahl, but we differ by starting with an electron-recoil-motivated Thomas-Imel parameter at zero electric field instead of one found from fits to NR data . Therefore, we offer a novel approach.
Our two predictions in Figure 9 are consistent with recent measurements by Plante et al.  and Manzur et al. , but in direct conflict with the older claims of Sorensen et al.  and Aprile et al. . Because NR yields change very weakly with field [22, 54], we did not initially investigate field dependence. Dahl’s study suggests an explanation for this weak dependence: 1 for NR at non-zero field and perhaps is field-dependent, not a fixed 0.06 as in ER . (Energy dependence at low energy is also possible ). If a significant fraction of interactions is due to excitations, then recombination probability is less important. For zero field we saw no need to change the or used successfully earlier. To reproduce field dependence, which is of greater importance now as S2 is used to lower the threshold of xenon detectors , we suggest employing a that changes differently with electric field than it does for ER, rather than attempting a change in . Though this is in apparent contradiction with Dahl’s work, we offer it as another means to the same end. The ability to distinguish between microphysics interpretations (changing or changing ) is likely beyond the reach of macroscopic empirical measurements that observe yields.
We have presented a coherent and comprehensive approach towards understanding the scintillation and ionization processes in liquid xenon. Starting with an ansatz-driven approach, aided by physically motivated models, we have been able to not only explain a large majority of the world data, but make predictions in regions where measurements are sparse. Our model is especially applicable to the low recoil energy regimes of interest for direct dark matter searches, as demonstrated by our prediction of . NEST will be made available to the world community to be used as an add-on to Geant4. Changing a few parameters in NEST, not all free (, , , etc.), will make it possible to simulate noble elements other than xenon. We have also suggested a universal nomenclature () in order to unify different definitions in the literature, which would enable the community to compare and contrast results which are often quoted in varying and confusing fashions. Following in the footsteps of other recent works striving for a unified view [10, 13, 27], we suggest the use of a mean work function (simply labeled as ) as the standard for reporting results and as the benchmark for quoting relative yields. We have provided a framework for simulating electron and nuclear recoils at different energies and applied electric fields, and have confirmed the results against a large sample of data from literature.
Acknowledgements.This work was supported by U.S. Department of Energy grant DE-FG02-91ER40674 at the University of California, Davis as well as performed in part under the auspices of the Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. We would like to thank C. E. Dahl for numerous in-depth discussions and paper draft reviews, but especially for laying the foundation for our paper through his seminal work with T. Shutt and J. Kwong at Case Western Reserve University. We also thank L. Baudis, A. Bolozdynya, L. Kastens, K. Ni, K. Masuda, D.N. McKinsey, and M. Yamashita for answering questions about their past work and for forwarding work with which we were not familiar or could not easily acquire, and R. Svoboda for consulting on drafts and reviewing them.
- Linear energy transfer LET, defined as the quotient of and density, is often used in place of .
- The latter effect is well-known to be small, sub-Poissonian in fact [2, 10, 26].
- This work function is also known as  or , conflicting with others’ use of the same nomenclature .
- This efficiency should not be confused with the light collection efficiency of a detector.
- 0.64 0.03 relative yield with respect to eV .
- Just like the other work functions, the notation is inconsistent. This one is variably called simply [2, 7], or . Our is at times denoted , , , or [18, 42].
- R. W. Schnee, Introduction to dark matter experiments, arXiv:1101.5205.
- E. Aprile and T. Doke, Liquid xenon detectors for particle physics and astrophysics, Rev. Mod. Phys. 82 (2010) 2053.
- K. Ni, E. Aprile, K. L. Giboni, P. Majewski, and M. Yamashita, Gamma ray spectroscopy with scintillation light in liquid xenon, J. Inst. 1 (2006) P09004.
- T. Doke, A. Hitachi, J. Kikuchi, K. Masuda, H. Okada, and E. Shibamura, Absolute scintillation yields in liquid argon and xenon for various particles, Jap. J. Appl. Phys. 41 (2002) 1538.
- S. Agostinelli et. al., GEANT4: A simulation toolkit, Nucl. Instrum. Meth. A506 (2003) 250.
- H. Tawara, T. Sasaki, K. Saito, and E. Shibamura, A Monte-Carlo method for determining absolute scintillation-photon yields and energy resolution of scintillators for gamma rays, in Second International Workshop on EGS, KEK, Aug., 2000.
- E. Aprile, A. E. Bolotnikov, A. I. Bolozdynya, and T. Doke, Noble Gas Detectors. Wiley-VCH Verlag GmbH and Co. KGaA, 2006.
- A. Manzur, A. Curioni, L. Kastens, D. N. McKinsey, K. Ni, and T. Wongjirad, Scintillation efficiency and ionization yield of liquid xenon for monoenergetic nuclear recoils down to 4 keV, Phys. Rev. C 81 (2010) 025808.
- T. Doke, K. Masuda, and E. Shibamura, Estimation of absolute photon yields in liquid argon and xenon for relativistic (1 MeV) electrons, Nucl. Instrum. Meth. A291 (1990) 617.
- C. E. Dahl, The physics of background discrimination in liquid xenon, and first results from XENON10 in the hunt for WIMP dark matter. PhD thesis, Princeton University, Jan., 2009.
- T. Doke, A. Hitachi, S. Kubota, A. Nakamoto, and T. Takahashi, Estimation of Fano factors in liquid argon, krypton, xenon and xenon-doped liquid argon, Nucl. Instrum. Meth. 134 (1976) 353.
- J. Lindhard, V. Nielsen, M. Scharff, and P. V. Thomsen, Integral equations governing radiation effects, Mat. Fys. Medd. Dan. Vid. Selsk. 33 (1963) 1.
- P. Sorensen and C. E. Dahl, Nuclear recoil energy scale in liquid xenon with application to the direct detection of dark matter, Phys. Rev. D 83 (2011) 063501, [arXiv:1101.6080].
- T. Takahashi et. al., Average energy expended per ion pair in liquid xenon, Phys. Rev. A 12 (1975) 1771.
- M. Miyajima et. al., Average energy expended per ion pair in liquid argon, Phys. Rev. A 9 (1974) 1438.
- A. Manalaysay et. al., Spatially uniform calibration of a liquid xenon detector at low energies using Kr, Rev. Sci. Instrum. 81 (2010) 073303, [arXiv:0908.0616].
- E. Aprile, K. L. Giboni, P. Majewski, K. Ni, and M. Yamashita, Observation of anticorrelation between scintillation and ionization for mev gamma rays in liquid xenon, Phys. Rev. B 76 (2007) 014115.
- T. Doke et. al., LET dependence of scintillation yields in liquid argon, Nucl. Instrum. Meth. A269 (1988) 291.
- J. B. Birks, The theory and practice of scintillation counting, International Series of Monographs on Electronics and Instrumentation 27 (1964).
- J. Thomas and D. A. Imel, Recombination of electron-ion pairs in liquid argon and liquid xenon, Phys. Rev. A 36 (1987) 614.
- A. I. Bolozdynya, A. W. Bradley, P. P. Brusov, C. E. Dahl, J. Kwong, and T. Shutt, Using a wavelength shifter to enhance the sensitivity of liquid xenon dark matter detectors, IEEE Transactions on Nuclear Science 55 (2008) 1453.
- E. Aprile et. al., Simultaneous measurement of ionization and scintillation from nuclear recoils in liquid xenon for a dark matter experiment, Phys. Rev. Lett. 97 (2006) 081302.
- E. Conti’s liquid xenon age. http://www.pd.infn.it/conti/LXe.html, 2011.
- A. Mozumder, Free-ion yield and electron-ion recombination rate in liquid xenon, Chem. Phys. Lett. 245 (1995) 359.
- C. E. Dahl. Private communications, 2010-11.
- J. Séguinot, J. Tischhauser, and T. Ypsilantis, Liquid xenon scintillation: photon yield and Fano factor measurements, Nucl. Instrum. Meth. A354 (1995) 280.
- T. Shutt, C. E. Dahl, J. Kwong, A. Bolozdynya, and P. Brusov, Performance and fundamental processes at low energy in a two-phase liquid xenon dark matter detector, Nucl. Instrum. Meth. A579 (2007) 451, [astro-ph/0608137].
- ESTAR. http://physics.nist.gov/PhysRefData/Star/Text/ESTAR.html, 2011.
- I. M. Obodovskii and K. T. Ospanov, Scintillation output of liquid xenon for low-energy -quanta, Instruments and Experimental Techniques 37 (1994) 42.
- T. Doke, E. Shibamura, S. Kubota, and K. Terasawa, Maximum scintillation yields in NaI(Tl) and CsI(Tl) crystals estimated from the scintillation model for liquid rare gases, Nucl. Instrum. Meth. B266 (2008) 5063.
- R. B. Murray and A. Meyer, Scintillation response of activated inorganic crystals to various charged particles, Phys. Rev. 122 (1961) 815.
- F. Salvat, J. M. Fernandez-Varea, and J. Sempau, eds., PENELOPE 2001 - A code system for Monte Carlo simulation of electron and photon transport, (Issy-les-Moulineaux, France), PENELOPE, Nuclear Energy Agency, Nov., 2001.
- I. R. Barabanov, V. N. Gavrin, and A. M. Pshukov, Liquid xenon scintillation spectrometer, Nucl. Instrum. Meth. A254 (1987) 355.
- K. Ni, R. Hasty, T. M. Wongjirad, L. Kastens, A. Manzur, and D. N. McKinsey, Preparation of neutron-activated xenon for liquid xenon detector calibration, Nucl. Instrum. Meth. A582 (2007) 569, [arXiv:0708.1976].
- V. Y. Chepel, M. I. Lopes, R. F. Marques, and A. J. P. L. Policarpo, Primary scintillation yield and ; ratio in liquid xenon, in Proceedings of the 1999 IEEE 13th International Conference on Dielectric Liquids, 1999. (ICDL ’99), p. 52, 1999.
- M. Yamashita, T. Doke, K. Kawasaki, J. Kikuchi, and S. Suzuki, Scintillation response of liquid Xe surrounded by PTFE reflector for gamma rays, Nucl. Instrum. Meth. A535 (2004) 692.
- M. Yamashita, T. Doke, J. Kikuchi, and S. Suzuki, Double phase (liquid/gas) xenon scintillation detector for WIMPs direct search, Astropart. Phys. 20 (2003) 79.
- T. Doke, R. Sawada, and H. Tawara, Non-proportionality of the scintillation yield in liquid xenon and its effect on the energy resolution for gamma-rays, in Technique and Application of Xenon Detectors (Xenon 01)., (Tokyo), p. 17, Proceedings of the International Workshop, 2001.
- P. Belli, R. Bernabei, C. Dai, A. Incicchitti, D. Prosperi, and C. Bacci, Liquid xenon detector for dark matter search at Gran Sasso, Nucl. Instrum. Meth. A336 (1993) 336.
- E. Aprile et. al., Design and performance of the XENON10 dark matter experiment, Astropart. Phys. 34 (2011) 679, [arXiv:1001.2834].
- M. Tanaka et. al., LET dependence of scintillation yields in liquid xenon, Nucl. Instrum. Meth. A457 (2001) 454.
- T. Doke and K. Masuda, Present status of liquid rare gas scintillation detectors and their new application to gamma-ray calorimeters, Nucl. Instrum. Meth. A420 (1999) 62.
- M. Yamashita. Private communications, 2011.
- T. Tojo, Nonlinear scintillation response of thin NaI(Tl) crystals, Nucl. Instrum. Meth. A238 (1985) 153.
- I. M. Obodovskii, Liquid xenon ionization spectrometers for the MeV region, Nucl. Instrum. Meth. A327 (1993) 119.
- J. Séguinot, G. Passardi, J. Tischhauser, and T. Ypsilantis, Liquid xenon ionization and scintillation studies for a totally active-vector electromagnetic calorimeter, Nucl. Instrum. Meth. A323 (1992) 583.
- M. Miyajima, S. Sasaki, and E. Shibamura, Comments on “liquid xenon ionization and scintillation studies for a totally active-vector electromagnetic calorimeter”, Nucl. Instrum. Meth. A352 (1995) 548.
- J. Séguinot, G. Passardi, J. Tischhauser, and T. Ypsilantis, Reply to ‘Comments on “Liquid xenon ionization and scintillation studies for a totally active-vector electromagnetic calorimeter” ’, Nucl. Instrum. Meth. A361 (1995) 623.
- E. Aprile, R. Mukherjee, and M. Suzuki, Performance of a liquid xenon ionization chamber irradiated with electrons and gamma-rays, Nucl. Instrum. Meth. A302 (1991) 177.
- S. Kubota, M. Hishida, M. Suzuki, and J. Z. Ruan, Dynamical behavior of free electrons in the recombination process in liquid argon, krypton, and xenon, Phys. Rev. B 20 (1979) 3486.
- D. Gastler et. al., Measurement of scintillation efficiency for nuclear recoils in liquid argon, arXiv:1004.0373.
- K. Palladino, Neutron-argon interactions 20 MeV: verification of ‘Neutron HP’ in Geant4, DUSEL AARM Collaboration Meeting, Feb., 2011.
- P. Amaudruz et. al., Simultaneous reconstruction of scintillation light and ionization charge produced by 511 keV photons in liquid xenon: Potential application to PET, Nucl. Instrum. Meth. A607 (2009) 668, [arXiv:0906.4075].
- E. Aprile et. al., Scintillation response of liquid xenon to low energy nuclear recoils, Phys. Rev. D 72 (2005) 072006, [astro-ph/].
- E. Conti et. al., Correlated fluctuations between luminescence and ionization in liquid xenon, Phys. Rev. B 68 (2003) 054201, [hep-ex/03].
- J. V. Dawson et. al., A study of the scintillation induced by alpha particles and gamma rays in liquid xenon in an electric field, Nucl. Instrum. Meth. A545 (2005) 690, [physics/0].
- A. S. Barabash, A. A. Golubev, O. V. Kazachenko, V. M. Lobashev, B. M. Ovchinnikov, and B. E. Stern, A pulsed liquid ionization chamber filled with Xe, Ar and CH, Nucl. Instrum. Meth. A236 (1985) 69.
- E. Aprile et. al., First dark matter results from the XENON100 experiment, Phys. Rev. Lett. 105 (2010) 131302.
- J. I. Collar, Light WIMP searches: The effect of the uncertainty in recoil energy scale and quenching factor, arXiv:1010.5187.
- G. Plante et. al., New measurement of the scintillation efficiency of low-energy nuclear recoils in liquid xenon, arXiv:1104.2587.
- P. Sorensen, A coherent understanding of low-energy nuclear recoils in liquid xenon, J. Cosm. Astropart. Phys. 2010 (2010) 033.
- A. Hitachi, Properties of liquid xenon scintillation for dark matter searches, Astropart. Phys. 24 (2005) 247.
- P. Sorensen et. al., The scintillation and ionization yield of liquid xenon for nuclear recoils, Nucl. Instrum. Meth. A601 (2009) 339, [arXiv:0807.0459].
- C. Savage, XENON10/100 dark matter constraints: examining the dependence, arXiv:1012.3926.
- P. Sorensen et. al., Lowering the low-energy threshold of xenon detectors, arXiv:1011.6439.