A review of NEST models for liquid xenon and an exhaustive comparison with other approaches

1 Introduction

For the past 15+ years, leading results in dark matter direct detection searches have been obtained from detectors based on the principle of the dual-phase Time Projection Chamber (TPC) using a liquefied noble element as the detection medium (Baudis, 2018). Liquid xenon (LXe) TPCs, in particular, have produced the most stringent cross-section constraints for Spin-Independent (SI) and neutron Spin-Dependent (SD) interactions between Weakly Interacting Massive Particles (WIMPs) and xenon nuclei. More recently, the use of LXe has also led to WIMP limits using different Effective Field Theory (EFT) operators for mass-energies above O(5 GeV) (Akerib et al., 2021a). EFT extends the set of allowable operators beyond the standard SI and SD interactions and includes searches at higher nuclear recoil energies. Unrelated to dark matter, electron recoil searches up to the MeV regime have set strict constraints on 0 ν β β decay (Anton et al., 2019) and led to observations of double e − capture (Aprile et al., 2019a). XENONnT and PandaX have recently illustrated the potential for precision measurements of ⁸B (Aprile, 2024a; Bo, 2024).

To interpret results from past, present, and future experiments, a reliable Monte Carlo (MC) simulation is required. Recent works have demonstrated the utility of NEST, the cross-disciplinary, detector-agnostic MC software reviewed in this study (Akerib et al., 2021b; Yan et al., 2021; Aprile et al., 2021), for a variety of active detector materials: LAr (Caratelli, 2022; Abud et al., 2023; Westerdale, 2024) and GXe, especially LXe. As the multi-tonne-scale TPCs have commenced data collection (Aalbers et al., 2023; Yan et al., 2021; Aprile et al., 2021), improved MC techniques will not only assist in limit setting but also be essential for determining the mass and cross section of dark matter particles in the event of a WIMP discovery. In either scenario or for the design of a new TPC, predictions of performance are needed on key metrics like the fundamental scintillation light and ionization charge yields for LXe, which is the focus of this work. NEST v2.4 is its default model; different versions are specified as needed. This manuscript is a technical overview of updates to NEST, including new models and comparisons. More pedagogical reviews of the models and related physics are available in the studies of Szydagis et al. (2011) and Szydagis et al. (2021a).

Section 2.1 presents the mean scintillation and ionization yields of electronic recoil (ER) backgrounds, along with comparisons to experimental data. These serve as the basis for the ER background (BG) models in Xe-based dark matter detectors. Section 2.2 summarizes the methods for varying these mean yields to model realistic fluctuations, with variations in the total number of quanta (light and charge) produced. Section 2.3 focuses on the yields of nuclear recoils (NRs) and their fluctuations. These form the foundation for the signal model in an LXe-based dark matter search, as well as for NR backgrounds (such as those from fast neutron scattering and coherent elastic neutrino-nucleus scattering, CEνNS). Lastly, Section 3 compares NEST’s modeling of mean yields (Sections 2.1 and 2.3) with past and present approaches in the existing literature, including some based on first-principles methods, before the conclusion. The strengths and weaknesses of the different approaches are summarized, underscoring NEST’s ability to phenomenologically model data across a broad range of energies and electric fields.

2 Microphysics modeling evaluation

The NEST model choices were justified earlier by Szydagis et al. (2021a) and in the references therein, but they are re-evaluated in this study more comprehensively with newer and more extensive datasets. NEST is openly shared, allowing for regular re-evaluation using the latest calibrations (Szydagis, 2020). Although such data often provide relative light and charge yields, these can be converted to absolute yields if the detector gains are calculable, known as g 1 and g 2 for these respective yields. The light yield gain, g 1, is the primary photon detection efficiency, while the charge gain, g 2, is the average signal size per e − escaping the interaction site. Uncertainties in these gains are a significant source of systematic error, but newer data from higher-quality calibrations help mitigate this issue. Combining calibration data ranging from <1 keV to > 1 MeV energy, NEST predicts the shapes of primary scintillation and ionization yields as functions of energy, E, and drift electric field, E, for different particle interaction types (Conti et al., 2003). The status of the NEST modeling of these shapes is shown in Figure 1.

2.1 Electronic recoils (beta, gamma, and X rays)

NEST begins with a model of the total yield, summing the vacuum ultraviolet (VUV) scintillation photons and ionization electrons produced. IR photons are not included as their yield in LXe is lower by a factor of ∼4 (Bressi et al., 2001), and their wavelength is beyond the sensitivity of most photon sensors commonly used in dark matter experiments. The work function, W q, for the production of quanta depends only on the density, determined using a linear fit based on data collected by Aprile et al. (2008) across different phases (see also Supplementary Appendix SA):

W q e V = 21.94 − 2.93 ρ = W i 1 + N e x / N i . (1)

Here, ρ is the mass density in units of g/cm³. LXe TPCs typically operate at temperatures of 165–180 K and pressures of 1.5–2 bar(a), leading to ρ ≈ 2.9 g/c m 3 and resulting in a W q value between 13 and 14 eV [Equation 1, with discrepant values discussed by Szydagis et al. (2021b)]. The exciton–ion ratio or N e x / N i relates W q to the work function for ionization, W i, which was defined for the charge yields. Moreover, N e x / N i determines the pre-recombination (of e −s with ions) split of quanta into light and charge (see Supplementary Appendix SA, where ρ dependence is explained):where E is the deposited energy in keV for a β interaction or Compton scatter and “erf” refers to the error function. Here, the ρ dependence is based again on Aprile et al. (2008), while the E dependence comes from reconciling Doke et al. (2002); Akerib et al. (2016a); and Lin et al. (2015), given the lines of evidence that light yield approaches 0 as energy E decreases, with lower-E data sets favoring both less recombination and smaller N e x / N i. Ionization electrons can recombine with Xe atoms or escape entirely from the interaction site. Therefore, the number of photons N p h is not simply equal to N e x, providing an anti-correlation between the observed light and charge yields; this motivates the use of both charge and light to measure the energy, E = W q ( N p h + N e − ) (Szydagis et al., 2021a):where r is the recombination probability for e −-ion pairs depending on E, E, and ρ, as well as the particle and interaction type, and S1 and S2 are the experimental observables. Typical values for g 1 are ∼0.1 but O(10) for g 2 due to secondary (gas) scintillation (g 2 is 0.5–1 in single-phase TPCs). The light and charge yields per unit energy are traditionally quoted in experiment, defined as L y ≡ N p h / E and Q y ≡ N e − / E, respectively.

N e x / N i = 0.0674 + 0.0397 ρ × e r f 0.05 E , (2)

N p h = N e x + r E , E , ρ N i = S 1 / g 1 a n d N e − = 1 − r E , E , ρ N i = S 2 / g 2 , (3)

Q y is modeled first; L y is set by W q and subtraction:where N q is the total number of quanta. This procedure leverages the greater reliability of S2 measurements compared to S1 for lower E, as explained by Akerib et al. (2017a) and Szydagis et al. (2021a). Q y in the ER model is a sum of two sigmoids:with m 1 serving as the minimum field-dependent charge yield. m 2 determines the low-E behavior, and m 7 controls the field dependence at high energies. The individual m i values are summarized in Supplementary Appendix SB (with Akerib et al. (2020a) providing more details). Although empirical, the first (left, m 1+…) and second (right, m 5−…) sigmoids of Equation 5 capture the qualitative behavior of two first-principles options, respectively: the Thomas–Imel box model at low energies (Thomas and Imel, 1987) and Doke-modified Birks’ law at higher energies (Doke et al., 1988). Between ∼15 keV and the energy of a minimally ionizing particle (MIP) within Xe (approximately 1 MeV), a track shape is described as cylindrical by Doke for modeling the recombination, and d E / d x decreases with increasing E. The recombination probability r decreases as energy E increases, reducing the ratio of L y to Q y (Szydagis et al., 2022; Szydagis et al., 2011; Berger et al., 2005). Below ∼15 keV, deposits are more amorphous, and straight 1-D track lengths become ill-defined: r and L y increase with the 3-D ionization density and the energy as d E / d x increases with E.

N q ≡ N e x + N i = N p h + N e − = E / W q , w h e r e N e − = Q y E , a n d N p h = N q − N e − , (4)

Q y E , E = m 1 E + m 2 − m 1 E 1 + E m 3 E m 4 E m 9 + m 5 E − m 5 E 1 + E m 7 E m 8 m 10 E , (5)

A Thomas–Imel approach historically uses E and models energy deposits within symmetric boxes or spheres, while the Doke/Birks’ law uses d E / d x and assumes long tracks (cylinders). The former will exhibit r (and therefore L y) only increasing with energy, while the latter will usually exhibit it decreasing, with Q y anti-correlated again.

The recombination fraction or probability, r, is found retroactively in recent NEST versions after fitting to Q y per Equation 5, chosen for matching both the box and Birks models. Using Equation 2 as a constraint avoids the degeneracy of this r with N e x / N i, with the sum N e x + N i (also equal to N p h + N e −) already constrained by Equations 1, 4—the former determines W q, and the latter determines total quanta N q based on W q. Any change in W q (one work function averaging over individual work functions for photon and electron production) should change L y and Q y equally, preserving both their shapes in both energy and field (Anton et al., 2020).

Figure 1 summarizes both L y and Q y for βs and Compton scattering ERs from both data and NEST, with NEST using typical LXe operating conditions of ρ = 2.89 g/c m 3 (T = 173 K and P = 1.57 bar). The non-monotonic energy dependence is obvious. Meanwhile, L y decreases from left to right (top), and Q y correspondingly increases (bottom) as the field increases, suppressing recombination at a fixed N q. However, even at E = 0, there exists a “phantom” Q y, likely caused by an extreme delay in recombination, as explained by Doke et al. (2002) and Szydagis et al. (2021a); this is unobservable, except via long S1 integration times, and by noting that L y vs. energy maintains the same shape at all fields, even at 0. This implies a continuous change in L y as E → 0. Non-zero fields standing in for 0 represent residual stray fields in a detector and/or inherent fields of Xe atoms (Szydagis et al., 2013).

The absorption of any high-energy photon, a γ or x-ray, is modeled as β interactions and Compton scatters but with unique m i (Figure 2) to capture sub-position-resolution multiple scatters and distinct d E / d x. L y is mostly lower and Q y is higher for βs, as explained within the Figure 2 caption. Although it might be possible to merge the γ and β models by relying on differences in d E / d x, γs are treated independently at present. Supplementary Appendix SB lists the β and γ model parameters, in addition to those for NR models.

2.2 Yield fluctuations

Energy resolution typically refers to Gaussian spreads (σ or FWHM) of monoenergetic peaks from high-energy γ-ray photoabsorption, but this is also relevant to lower energies in WIMP searches. The smearing of continuous ER spectra can drive an increase in signal-like background events. However, to understand statistical limitations for high-level parameters like monoenergetic-peak σs or background discrimination, we must start with lower-level parameters that underlie all the relevant stochastic processes involved. This modeling is discussed in depth by Szydagis et al. (2021a), but portions relevant to this work are summarized in this section, culminating in a subsection enumerating the practical steps taken within the NEST code on GitHub.

2.2.1 Total quanta: correlated fluctuations

Realistic smearing of mean yields begins with a Fano-like factor, F q, applied to the total quanta, N q, prior to differentiation into N e x and N i. It is labeled as Fano-like as it does not follow the strict sub-Poissonian definition (Doke et al., 1976). F q may exceed 1, but it is still used in the usual definition of the standard deviation of N q, utilized for decades by Xe experiments to fit their data on combined E (N p h and N e −) scale resolution:where F q is defined for light and charge together as

σ q = F q ⟨ N q ⟩ , (6)

F q = 0.13 − 0.030 ρ − 0.0057 ρ 2 + 0.0016 ρ 3 + δ F ⟨ N q ⟩ E . (7)

The first part of Equation 7 is a spline of data (Aprile et al., 2008) from gas, liquid, and solid. The constant 0.13 represents the theoretical value of the Xe Fano factor, following the traditional definition ( F q < 1 ). O(0.1) matches NEXT gas data on N e − (Alvarez et al., 2013) and Biagi’s Degrad work. The second part of Equation 7 is only for liquid and is data-driven, where δ F = 0.0015 for LXe but is identically 0 for gaseous Xe. The ⟨ N q ⟩ term is included in order to match the data at MeV scales (e.g., for 0 ν β β searches). Such results did not achieve the theoretical minimum in energy resolution even when reconstructing N q, utilizing both channels of information (light and charge), instead of only a single channel. This was true even for the cases where the noise was allegedly subtracted or modeled (Delaquis et al., 2018; Aprile et al., 2020a). As Q y increases with E, the combined E resolution improves. However, the improvement is smaller than naïvely predicted, requiring the E term in F q to match the data (Aprile et al., 2007; Aprile et al., 1991).

There are many possible explanations for F q becoming ≫ 1 as E or E changes. W q may need to be replaced with separate W e x and W i for the excitation and ionization processes (both inelastic scattering), respectively, and then further subdivided into different values that depend on the e − energy shell. Lastly, elastic scattering of orbital e −s may play a role. These mechanisms are discussed by Platzman (1961), but explicit Fano-factor variations can be found in Szydagis et al. (2021a). In NEST, a Gaussian smearing, constrained to be non-negative, is applied to N q with a width defined by Equation 6: N q = G [ ⟨ N q ⟩ , σ q ]. A binomial distribution then divides quanta into excitons versus ions.

2.2.2 Anti-correlated excitation and recombination fluctuations

F q drives resolution on a combined E scale, but such a scale is more relevant for monoenergetic peaks than dark matter searches (Dahl, 2009; Szydagis et al., 2021a). “Recombination fluctuations,” however, describe the redistribution of N p h and N e − caused by widths associated with the means of Equations 3, 4. Often conflated with excitation fluctuations (Equation 2), these are all fundamental and do not originate from detector effects (Aprile at al., 2011; Akerib et al., 2017b); they constitute one of the key factors for the characterization of ER discrimination (Dobi, 2014). Moreover, they are not binomial, despite recombination (or escape) appearing to be a binary decision. Potential explanations for this phenomenon include other energy loss mechanisms, or other effects that break the independence of draws, for instance, δ-ray production (as observed at different energies in both Ar and Xe (Amoruso et al., 2004; Thomas et al., 1988)), the statistics of columnar recombination (Nygren, 2013), and short-lived clustering of Xe dimers (Davis et al., 2016).

While it remains unclear which explanation is correct, NEST proceeds with a fully empirical approach to simply model what is observed in data; following the works by Akerib et al. (2017b) and Akerib et al. (2020a) closely, NEST defines recombination variance as follows:

σ r 2 = ⟨ r ⟩ 1 − ⟨ r ⟩ N i + σ p 2 N i 2 , w h e r e σ r ≈ σ N e − ≈ σ N p h f o r s m a l l F q , σ p = A E e − ⟨ y ⟩ − ξ 2 2 ω 2 1 + e r f α p ⟨ y ⟩ − ξ ω 2 , a n d t h e e − f r a c t i o n y = N e − / N q a n d ⟨ y ⟩ = Q y W q . (8)

⟨ r ⟩ ( 1 − ⟨ r ⟩ ) N i in σ r follows the binomial expectation of σ r ∝ N i. The σ p term leads to σ r ∝ N i, as proposed by Dobi (2014). σ p is a skewed Gaussian (on the third line) with field-dependent amplitude, A, varying from 0.05 to 0.1, as needed to simulate the spectral broadening of ER with higher drift electric field (Akerib et al., 2020a; Akerib et al., 2020b). In NEST versions < 2.1, σ p was simulated as a constant, similar in value to A, but this was found to be inadequate for capturing the full behavior of recombination fluctuations (Akerib et al., 2017b).

σ p’s dependent variable was chosen to be the mean electron fraction ⟨ y ⟩ for simplicity as it is closely related to 1− ⟨ r ⟩. Recombination probability, defined within Equation 3, is degenerate with N e x / N i, while y is directly measurable. It can be written in terms of r: y = ( 1 − r ) / ( 1 + N e x / N i ) (Dahl, 2009). Non-binomial fluctuations decrease as y approaches 0 or 1, causing σ p to vanish. ξ, ω, and α p are the centroid, width, and skew of σ p, respectively. Default NEST values determining the width and skewness of σ p are ω = 0.2 and α p = − 0.2, respectively (future work may recast σ r entirely in terms of y not just σ p).

A skew centroid ξ ≈ 0.4–0.5 was found based on β and γ datasets. The types of datasets included continuous spectra and monoenergetic-peak energy resolutions, both at multiple fields and energies (Dahl, 2009; Aprile et al., 2011; Dobi, 2014). ξ’s value depends on which datasets are used and which other parameters are fixed. A ξ near 0.5 leads to a maximum in σ p (within σ r) near y = 0.5, as would occur within a regular binomial distribution. The asymmetric shape σ p is motivated by observations of recombination fluctuations at lower values of y (low field, high energy) compared to higher values of y (high field, low energy) (Rischbieter, 2022; Dobi, 2014; Akerib et al., 2020a).

Longer, less technical descriptions of all the steps in Section 2.2.2 can be found in the studies by Akerib et al. (2020a) and Rischbieter (2022).

2.2.3 Recombination skewness

We note that the skewed Gaussian σ p ( y ) must not be conflated with E and E-dependent skew defined in Section IVB of Akerib et al. (2020b) as α r; the skew in that study represented the observed asymmetry of the resultant charge yields. NEST uses α r from Equation 13 in Akerib et al. (2020b) to smear the mean N e −, while α p controls the variance of recombination fluctuations, σ r, as described in Equation 8.

A positive α r value can lead to better background discrimination than expected for a WIMP search that uses LXe. Weak rejection was expected due to the recombination fluctuations being greater (worse) than binomial, but positive α r will shift ER events preferentially away from NR (more Q y). This has already been observed by Akerib et al. (2020b).

2.2.4 Uncorrelated fluctuations: detector effects (known and unknown)

Lastly, while the simulated σ q widths predict correlated changes in S1 ( L y ) and S2 ( Q y ) and σ r leads to an anti-correlated change, uncorrelated noise also exists, affecting S1 and S2 independently. S1 and S2 gains are understood sources, assuming position-dependent light collection and field non-uniformities are taken into account. Unknown sources are modeled with a Gaussian smearing proportional to the pulse areas (Szydagis et al., 2021b). A quadratic term may be necessary at the MeV scale (Davis et al., 2016). ER and NR are equally affected by any detector effects (known/unknown). The final E resolutions vs. E are observed for ER, NR, or both (Akerib et al., 2021b; Szydagis et al., 2021b), supplementing the validation of means in Figures 1–3 with their vetting of fluctuations. The scale of the unknown detector effects across experiments is 1%–10% (Szydagis et al., 2021b; Szydagis et al., 2021a; Aalbers et al., 2024) (for S2s and non-integer forms of S1s) but effectively 0% for a spike count of S1 photons. For further details, refer to Supplementary Appendix SA.

2.2.5 Computational implementation

NEST is publicly available as a GitHub repository, which includes the source code, interface scripts, and examples. It is C++-based but can be run with dedicated scripts using either C++ or Python, both of which are available in the repository. These can be used to generate expectation values of yields and their fluctuations for different detectors using Xe or Ar. The step-by-step procedure that NEST follows to perform these tasks is summarized below:

• • F q is used to determine σ q for a normal distribution of total (initially undifferentiated) ER quanta, which can be considered “correlated noise” because, in this case, N p h and N e − increase and decrease together [Eq A1 (Aprile, 2024b)]. Two distinct Fs exist for NR N e x and N i, breaking the correlation (Section 2.3).
• • ER quanta are differentiated (N e x and N i) using a binomial distribution [Eq A2 (Aprile, 2024b)], approximated as normal for computational speed, using the same Box-Muller algorithm as in the first step above. Any non-binomial/ non-Gaussian fluctuation at this stage is essentially degenerate with the next step.
• • A normal or skew-normal [Eq 8–12 (Akerib et al., 2020b)] in N e − capped at N q (minimum of 0) enforces the anti-correlated fluctuation of N p h versus N e −. This step was previously mismodeled by uncorrelated Fano factors. The variance σ r 2 has components proportional to both N i (“binomial style”) and N i 2 (data-driven).

Two more lists cover detector specifics for S1 and S2, closely following Supplementary Appendix SC of Aprile (2024b). First, S1 comprises the following:

• • S1.1 A binomial distribution with probability g 1 (3-D spatially varying) determines the fraction of N p h successfully detected by photo-sensors; g 1 represents the product of geometric × quantum efficiencies.
• • S1.2 Single photo-electrons in sensors are modeled by zero-truncated Gaussians of sensor-specific width. Spike counting is emulated using artificially reduced width but non-zero for matching real data.
• • S1.3 An if-else structure determines whether a second photoelectron is produced due to the secondary PE effect. This step and S1.2 are Gaussian-approximated at high E in the “hybrid” mode or any E in the “parametric” mode.
• • S1.4 Geant4 (G4), Chroma, OptiX, or some other ray-tracer, or NEST’s built-in analytic-approximation ability simulates photon arrival times at S1 sensors and dictates whether a sufficient number of photons were detected in MC with above-threshold (experiment DAQ-specific) pulse areas, based upon stages S1.2 and S1.3 above.

The procedure to model the charge signal or S2 is more intricate, especially in a two-phase experiment:

• • S2.1 Electrons (numbered N e −) diffuse both transversely and longitudinally as they drift at a drift speed determined by the liquid’s field but also influenced by factors such as density, temperature, and pressure (the same applies to diffusion “constants”). Data-driven functions exist for all these phenomena in NEST.
• • S2.2 An electron survival fraction is set by an exponential function depending on the originating depth in a detector and a characteristic electron MFP. It is used as the probability in a binomial distribution.
• • S2.3 Another binomial distribution is utilized to find how many electrons survive extraction from the liquid to the gas. The efficiency is a function of the gas field E g between the liquid/gas boundary and gate grid. NEST offers many options of asymptotic (1 at infinite E g) function based on the past data.
• • S2.4 Each extracted electron produces Y e − S2 photons based on the parameterization described by Chepel and Araújo (2013) depending on E g, gas ρ, and the gap between the liquid surface and gate (thus, E g comes into play twice). Y e − is the mean of an integer-rounded Gaussian with a width of F S 2 Y e −. F S 2 is O(1) and captures grid non-uniformity.
• • S2.5 A binomial of probability g 1 g a s (2-D varying) similar to g 1 is step 1 of a process similar to S1.1–4.

More precise S2 simulation is possible in the optional integration of Garfield with NEST, which also possesses an optional G4 integration for simulating E deposits prior to the first step above. More details on the lists here can be found in Section 2.2 of Szydagis et al. (2021a). Section 2.2.4 explains NEST’s last layer. All values for the first list are provided in Supplementary Table S4 (Supplementary Appendix SB), and examples for S1 and S2 are provided by Rischbieter (2022), especially in Figure 4.3 left.

2.3 Nuclear recoils (neutrons and WIMPs and Boron-8)

NR N q ′ (differentiated in this section from ER with a prime) is well-fit by a power law across >3 orders of magnitude in E [Figure 5 in Szydagis et al. (2021a)]. This is a simplification of the Lindhard approach to modeling the reduced quanta compared with ER but also allows for departures from Lindhard at higher Es, lowering N q ′ ( E )’s rate of change with respect to Lindhard. Fewer equations and parameters are involved compared to Lindhard, which is a combination of multiple power laws inside a rational function (Lindhard, 1963); see Equation 8 in Szydagis et al. (2021a) for more justification. NEST uses that simpler formula:

N q ′ = a E b , w h e r e a = 1 1 − 0.5 + 2.0 a n d b = 1.1 ± 0.05 . (9)

The uncertainties here are > 10 × those reported recently for the same fit as only statistical error was included in Equation 6 of Szydagis et al. (2021a). In this study, systematic uncertainties in S1 detection efficiency and S2 gain (including e − extraction efficiency) are included. They can be found inside the individual references in the caption of Figure 3. Individual power laws were found for each dataset prior to the error-weighted combination so that a dataset with more points was not overly weighted. Equation 9 was also cross-checked with L y ′ and Q y ′ individually extracted from data, as displayed in Figure 3, and the raw S1 and S2 data on continuous energy spectrum sources.

Equation 9 can be used to define “quenching,” L, in Equation 10:which is interpreted as the fraction of total NR energy shared with the electron cloud to produce ions and excitons. L permits one to define the electron equivalent energy in units of keV e e for NR as L × E in keV n r), a best average reconstruction of the (combined-)E of recoiling nuclei. This L should be applicable to neutron calibrations, WIMPs, and CEνNS, such as from ⁸B nuclear fusion (Aprile et al., 2021).

L E , ρ = N q ′ E / N q E , ρ = N q ′ E W q ρ / E , (10)

While the previous equation sets the total quanta, the next equation determines the field- and density-dependent division into individual yields (charge or light) in an anti-correlated fashion, reducing r with higher field:

ς E , ρ = γ E δ ρ ρ 0 υ , w h e r e γ = 0.0480 ± 0.0021 a n d δ = − 0.0533 ± 0.0068 , a n d υ = 0.3 . (11)

The reference density is ρ 0 ≡ 2.90 g/c m 3. The value of 2.89 was a specific example using LUX; the differences in yields are negligible. The exponent υ for the density dependence is hypothetical. It is not well-measured at densities significantly deviating from ρ 0 (Dahl, 2009).

We use Equation 11 to produce a Q y ′ equation:

Q y ′ E , E , ρ = N e − p e r k e V = 1 ς E , ρ E + ϵ p 1 − 1 1 + E ζ η , w h e r e ϵ = 12 . 6 − 2.9 + 3.4 k e V , p = 0.5 , ζ = 0.3 ± 0.1 k e V , a n d η = 2 ± 1 . (12)

Energy deposited is again E (in keV), and ϵ is the reshaping parameter for the E dependence. Higher or lower ς decreases or increases the Q y ′ level, respectively, providing the field-dependent shape of Q y ′ ( E ). ε can be assumed to be the characteristic E where Q y ′ changes in its behavior from ∼constant at O(1 keV) to decreasing at O(10 keV) (note that ς has adaptable units of keV 1 − p).

ζ and η are the two sigmoid parameters that control the Q y ′ roll-off at sub-keV energies. They permit a better match to not only the most recent calibrations (Lenardo et al., 2019; Akerib et al., 2017c) but also to NEST versions pre-2.0 and other past models. Combining Thomas–Imel recombination with Lindhard [Equation 8 of Szydagis et al. (2021a)] produces a roll-off in Q y ′, but it is less steep than that observed in data. Here, η controls steepness, allowing for an improved modeling of low-energy NR (Szydagis et al., 2013; Sorensen and Dahl, 2011), while ζ represents a characteristic scale for NR to ionize one e − (Szydagis et al., 2021a; Sorensen, 2015). At high E, p = 0.5 reproduces Q y ′ ∝ 1 / E (Figure 3, bottom row).

Similar to ER, N p h is derived from N q ′ − N e −, but this is only a temporary anti-correlation enforcement; an additional sigmoid permits L y ′’s flexibility (Equation 13). Future calibration data could show a decrease or even flattening, potentially due to additional N p h from the Migdal effect (Akerib, 2016; Aprile et al., 2019c). An increase in L y ′ is possible even as E → 0. This is not unphysical as long as N p h vanishes in that limit, conserving E.

L y ′ ′ = N q ′ E − Q y ′ . N p h = L y ′ ′ E 1 − 1 1 + E θ ι ; L y ′ = N p h E , w h e r e θ = 0.3 ± 0.05 k e V a n d ι = 2 ± 0.5 . N q ′ = N p h + N e − . (13)

The top row of Figure 3, especially when read from right to left, shows the same L y ′ shape at all fields, once again indicative of a zero-field phantom Q y ′. In the L y ′ calculation, L y ′ ′ is a temporary variable (perfect anti-correlation) used within NEST to calculate the final L y ′ and N q ′ values. The best-fit numbers for θ and ι match those of their counterparts ζ and η for Q y ′. In this modular but smooth approach, the sigmoidal terms in L y ′ and Q y ′ approach 1.0 with increasing E. This method allows for separate fitting of the low- and high-E regimes, enabling the possibility of different physics in the sub-keV region, while avoiding the use of higher-E data to over-constrain lower-E yields.

The two sigmoids reduce the predictive power of NEST for extrapolation into newer, lower-E regimes where no calibrations exist. In the case of L y ′, it will be challenging to achieve any with low uncertainty.

θ is a physically motivated characteristic energy for the release of a single (VUV) photon. Like ζ, its value is 300 eV, in agreement with Sorensen (2015) and NEST pre-v2.0.0 (Szydagis et al., 2013). Fundamental physics models for the L governing total quanta, such as Lindhard (1963) and Hitachi (2005) and Aprile et al. (2006), coupled to the Thomas–Imel “box” model for recombination (Thomas and Imel, 1987), predict a similar value. A larger θ value means more E is needed to produce a single photon (as opposed to excitons), and L y ′ is lowered. This may potentially be detectable for an experiment with sufficient light collection efficiency.

Decreasing ι would also lower L y ′, halving L y ′ across all E when ι = 0. On the other hand, in the limit of infinite ι (and/or θ → 0), the effect of the sigmoid is entirely removed, increasing L y ′ at low E. The same is true for η and ζ in the Q y ′ formulation. A hard cut-off for any quanta was implemented in NEST for E < W q ( N q / N q ′ ) ≈ 200 eV. N q represents the quanta that would have been generated for same-E ER. Below this, no quanta are generated. Sub-keV recoils have been observed at 200–400 V/cm (Figure 3).

In contrast to ER, for which the data suggest strict anti-correlation, simulated ⟨ N q ′ ⟩ is not varied with a common Fano factor shared by both types of quanta for simplicity. For NR, there are (nominally) separate Fano factors for excitation and ionization, which can soften the strict anti-correlation at the level of the fundamental quanta. ⟨ N e x ⟩ is smeared using a Gaussian of standard deviation σ e x = F e x ⟨ N e x ⟩. ⟨ N i ⟩ is similarly varied using σ i = F i ⟨ N i ⟩, as is the standard practice for Fano factors (Fano, 1947). Based on the sparse existing reports of NR E resolution (Akerib, 2016; Lenardo et al., 2019; Plante, 2012), both F e x and F i are set to 0.4 in NEST (as of v2.3.11; 1 earlier) although some data imply F e x ≫ 1 (Akerib, 2016; Plante, 2012). N e x = G [ ⟨ N e x ⟩ , σ e x ] and N i = G [ ⟨ N i ⟩ , σ i ] (G=Gauss).

Using the same functional form as in Equation 8 from ER, NEST models fluctuations in recombination for the redistribution of photons and electrons prior to measurable NR S1 and S2. The new parameters are distinguished using a prime symbol superscript again for NR ( ′ ).

Parameter values are similar but not identical to those from ER: A ′ = 0.04 (as of v2.3.11 and fixed for all fields), ξ ′ = 0.50, and ω ′ = 0.19 (α p ′ = 0). Over time, these appear to have been converging upon values similar to ER’s. These set a final recombination width σ r ′. N e − and N p h distributions have that width but are skewed due to NR recombination asymmetry (α r ′ = 2.25). α r ′ may be higher, but it is difficult to disambiguate NR skew (less L y) in data from unresolved multiple scatters, other detector effects (Akerib et al., 2020b), or Migdal effect ER, which can increase Q y and generate a secondary population (Akerib et al., 2019b).

3 Comparisons to first-principles approaches

By smoothly interpolating datasets taken at individual energies and/or electric fields, NEST is now fully empirical, built upon sigmoids and power laws as needed for a continuous model. However, inherent uncertainty is introduced by extrapolating into new energy and/or field regimes. To assess that and further validate an empirical approach, we show agreement with the models closer to “first principles.” Within NEST’s earliest versions, the Thomas–Imel (T-I) box model (Thomas and Imel, 1987) was used for low energy, while Birks’ law of scintillation was adapted for high energy. Both were qualitatively explained in Section 2.1 but are quantified in this section. The latter approach inside NEST was similar to Doke’s modification (Szydagis et al., 2011) for scintillation alone but applied directly to recombination, allowing it to model both L y and Q y:

⟨ r ⟩ = k A d E d x 1 + k B d E d x + k C , w i t h k C = 1 − k A / k B . (14)

This is Birks’ law for other scintillators (Birks, 1964) but with an additional constant k C that accounts for parent-ion recombination (Doke et al., 2002). Its constraint ensures that ⟨ r ⟩ is between 0 and 1 as it is a probability. A best fit to ER ( γ ) data has a non-zero k C only at 0 V/cm; at non-zero E, Equation 14 contains only one Birks’ constant, k A = k B.

k B’s best-fit value (for 180 V/cm) is 0.28 from a fit to only the high-E portion of the NEST β ER model. That model is, in turn, supported by ³H, ¹⁴C, and ²²⁰Rn data from LUX and XENON. Notably, k B in NEST v0.9x and the first NEST paper, 13 years ago, for this E was 0.257, within 10% of the value in Figure 4 (upper right plot pane), which covers many alternative approaches to NEST.

Despite Birks’ great success in explaining data at high E, that model cannot capture the behavior of ER at E ≲ 50 keV. Although lower-E extensions are possible, such as the addition of higher-order terms in d E / d x for that region, we instead consider the T-I model for lower E:

⟨ r ⟩ = 1 − l n 1 + ξ T I ξ T I , w h e r e ξ T I = N i 4 α T I a T I 2 v d . (15)

ξ T I parameterizes the physical principles. α T I describes diffusion, v d is the e − drift velocity, and N i is again the number of ions. Diffusion is modeled using the relation α T I = D e 2 / ( k T ϵ d ), where D combines e − and positive-ion diffusion coefficients, e is the elementary charge, k is the Boltzmann constant not Birks, T is temperature, and ϵ d = 1.85 × ϵ 0 is the dielectric constant. D = 18.3 cm 2/s is the longitudinal diffusion constant for e −s at 180 V/cm, derived from S2 pulse lengths (Sorensen, 2011b). e − diffusion dominates over cation diffusion. Assuming this D (and T = 173 K from earlier), ϵ d as defined above, and taking v d = 1.51 mm/μs at field E = 180 V/cm (Akerib et al., 2016b), we find α T I = 1.20 × 1 0 − 9 m 3/s. From this, the escape probability ( 1 − ⟨ r ⟩ ) for electrons inside a box is found by solving the relevant (Jaffé) differential equations (refer to Section 6.2 of Dahl (2009) for the details).

We interpret a T I, the size of the “box” surrounding ionized atoms, as corresponding to an (E-independent) e −-ion thermalization distance of 4.6 μm, as calculated by Mozumder (1995). This value was used before as a border in NEST for track length to switch from T-I to Birks. The ultimate value of TIB ≡ α T I / ( a T I 2 v d ) for that case is 0.0376.

Dahl found best-fit values of TIB ranging from 0.03 to 0.04 for both ER and NR data at 60–522 V/cm (Dahl, 2009). Our contemporary fits (for NEST and data), the blue lines at low energies in the first two panels at top in Figure 4, used 0.0300. If v d changes with the drift field [it is typically O(2 mm/μs) (Albert et al., 2017)], then the entire ranges described by Dahl, and by Sorensen and Dahl, are covered: 0.02–0.05 (Sorensen and Dahl, 2011).

For NR, Figure 4 (bottom row) presents many different past models, mainly for L y. NEST originally used T-I for NR, as described by Dahl (2009) and Sorensen and Dahl (2011), represented by the blue lines in Figure 5. This follows the same color convention as Figure 4. T-I fixes r, thus partitioning E into L y and Q y; however, the total yield must still be determined. For maximal distinction, we have selected the original Lindhard formula, as laid out by Lindhard (1963); Sorensen and Dahl (2011); Akerib (2016); and Szydagis et al. (2021a), rather than Equation 9. We set the crucial Lindhard parameter of k L to a value of 0.166, the decades-old default for Xe (Lindhard, 1963). Averaging over E, N q ′ / N q ≈ k L. It is observed that 0.166 is consistent with actual data (Akerib, 2016), Lenardo’s meta-analysis (Lenardo et al., 2015), and NEST v2.3+.

We identify ς from Equation 12 with the TIB value, as justified by Equation 11, where the parameters for the E dependence of ς (γ and δ) overlap at the 1σ level with the power-law field dependence of TIB from Lenardo et al. (2015). At 180 V/cm, ς = 0.0362, which is quite close to earlier theoretical calculation and comparable to a best-fit TIB for ER, assuming N e x / N i = 1.0. Although higher than for ER, it is the most common assumption for NR, and best-fit values from data and theory vary from 0.7 to 1.1 (Sorensen and Dahl, 2011).

An additional quenching is applied to just L y ′ (Manzur et al., 2010). We find a common parameterization of this effect (Bezrukov et al., 2011) to be defined in a manner analogous to Birks’ law or Equation 14:where q < 1 is a multiplicative factor on L y ′. ϵ Z is unitless reduced energy, useful for comparison between elements. Equation 16 is similar to Equation 14. The power law can be identified as proportional to NR d E / d x. If we define d E / d x (or LET) as approximately β ′ ϵ λ, then κ = k B β ′ L. Assuming ER k B (defined as 0.28 for 180 V/cm in Figure 4 top), L ∼ 0.15 (11/73) per an energy-independent approximation of Equation 9, justified by the power being close to 1, and with β ′ = 100, we obtain κ = 4.20, < 0.2 σ away from that determined by Lenardo et al. (2015). A fraction of the quanta removed from L y ′ in Equation 16 may be convertible into Q y ′. Figure 5 (right) explores that with the fraction set to 0.1.

q = 1 1 + κ ϵ Z λ , w i t h ϵ Z ≈ 1 0 − 3 E , (16)

Unlike with ER, Birks’ law models NR over the entire E range of interest (Figure 5, red), with k B = 0.28 and d E / d x = β ′ ϵ λ = 100 ϵ. Although there is disagreement about whether λ is 1.0 or 0.5 depending on the E regime (Hitachi, 2005; Aprile et al., 2006), 1.0 only differs by 1.6 σ from the value of 1.14, as determined by Lenardo et al. (2015).

Looking back at alternatives to Lindhard, Figure 4 shows that NEST’s power law models for L y ′ and Q y ′ align well with results from Mu et al. (2015) and Mu and Ji (2015), and Wang and Mei (2017) and Mei et al. (2008). NEST’s lower 1 σ line intersects with L y ′ determined by Sarkis et al. (2020), which is low due to the exclusion of more recent data points (Akerib, 2016; Akerib, 2022). On the higher-E L y ′ end, NEST’s upper uncertainty band encompasses results from neriX (Aprile et al., 2018b). For Q y ′, NEST lies in between higher values of Wang and Mei (2017) and lower values of Mu and Ji (2015) and Sarkis et al. (2020), also fitting between data from LUX D–D (Akerib, 2016) and LLNL (Lenardo et al., 2019).

The good agreement between the fully empirical NEST model and the first-principle models of both NR and ER shown in this study demonstrates that NEST can accurately simulate potential dark matter signals and backgrounds, respectively. This should hold true even for the regimes where data are still lacking, or they exist but have large uncertainties. In the case of NR, the fully empirical approach reproduces all data more accurately while using a comparable number of free parameters, offering much greater flexibility than semi-empirical approaches. For fluctuations, the number of NEST free parameters increased to two Fano factors (excitation and ionization) and four numbers for recombination width and skew to fully model the E resolution. NEST, justified first using data, is not limited to the operating conditions of previous experiments to make predictions relevant for a future experiment, although L y and Q y must pass through a detector simulation to obtain realistic S1 and S2 pulse areas: the processes in Section 2.2.5 in this study and Supplementary Appendix SA of James et al. (2022).

4 Discussion and future work

Beginning with our models of beta ER, gamma-ray ER, and the NR light and charge yields, along with resolution modeling, a coherent picture was built up inside the NEST framework, which enables a good agreement with data. NEST was also shown to have features from multiple first-principles approaches, such as the box and Birks models. NEST already works for LAr (Szydagis et al., 2021a) using the same formulas as LXe but with unique parameter values. However, it still only works best for point-like interactions, like those in dark matter experiments like DarkSide, not tracks, as will be observed by DUNE. The list of NEST collaborators includes TESSERACT (Biekert et al., 2022) members, so the addition of liquid helium (LHe) to NEST is planned.

Looking beyond LHe, short-term future work includes NEST re-writing to account for the lower W q measured by EXO and Anton et al. (2020); Baudis et al. (2021), but this will be easier if NEST can return to approaches closer to first principles. Therefore, a concerted effort will be made to revisit a semi-empirical formulation through the application of a modified T-I model, as pioneered by ArgoNeuT (Acciarri et al., 2013); this approach will incorporate a literal breakup of long tracks into boxes, as described in the thesis of Dahl, allowing higher energies to exhibit lower light yields without hard-coding this behavior, by virtue of being composed of multiple lower-E interaction sites. High-E modeling is thus accomplished by having one model for all Es but treating high-E interactions as a series of many low-E fragments, where L y will continue to be monotonically increasing with E. The main motivation for this is greater confidence in extrapolations to uncalibrated regions of future detectors.

The modified box model of LArTPC-based high-E neutrino experiments should also be useful for LXe NR. We demonstrate, herein, how it represents a more generalized version of the current NR model:where a ′ ≡ 1 in default T-I [but relaxing this constraint to O(1) as per Acciarri et al. (2013) can better fit data], ξ is short for ξ T I, redefined as β d E / d x with β as a constant (not Equation 15), and α x ≡ N e x / N i for conciseness.where we employ, in order, the approximations b ≈ 1, α x ≈ 1, and a ≈ 10 (Equation 9). Fitting to the SRIM line in Figure 5 of Aprile et al. (2006), one finds that for NR, in normalized (dimensionless) units, stopping power iswhich is valid in the range of 0–100 keV. However, near 50 keV, a square root function with an offset also fits SRIM: 3.4 E + ϵ, with ϵ = 12.6 keV (Equation 12). Making the ansatz β ∼ ς (Equation 11),recovering the high-E portion of Equation 12 at ς = 0.036 (200 V/cm) and E = 50 keV, given Equations 18–20. By modifying the power law for N q to be a E b − C (McMonigle, 2024), it may be possible to eliminate the need for the sigmoids for reducing both Q y and L y at the lowest Es, combining C with an additional degree of freedom, a non-unity a ′ in the natural log. By replacing our present Equations 12, 20 with Equation 17, we should be able to find a sufficiently flexible compromise that fits data with the same number of free parameters or fewer even (eliminating the sigmoid roll-offs and the ϵ offset in d E / d x potentially), all motivated from first principles (T-I). The redefinition of ξ T I in terms of d E / d x permits a non-linearity in the dependence of ξ on N i and an incorporation of d E / d x (as in the Doke/ Birks’ law), while α x could be made E and E-dependent as in Eq. B8 of Aprile (2024b), if absolutely necessary, following the similar increase with E for ER in Equation 2 [mimicked by Eq. A4’s exponential in Aprile (2024b)]. Lastly, the replacement of a E b − C with E / W q in Equation 17 could permit usage for ER, as in LAr, from the keV to the GeV scales.

Q y ′ = N e − E = 1 − r N i E = l n a ′ + ξ ξ N i E = l n a ′ + ξ ξ N q / E 1 + α x = l n a ′ + ξ ξ E a E b 1 + α x , (17)

Q y ′ = a E b − 1 1 + α x l n 1 + β d E d x β d E d x ≈ a 1 + α x l n 1 + β d E d x β d E d x ≈ a 2 l n 1 + β d E d x β d E d x ≈ 5 l n 1 + β d E d x β d E d x , (18)

d E / d x = 120 ϵ Z = 120 0.001 E = 120 0.001 E = 3.8 E , (19)

Q y ′ ≈ 5 l n 1 + ς d E d x ς d E d x = 5 l n 1 + 0.036 3.4 E + 12.6 0.036 3.4 E + 12.6 = 5 3.4 0.677 ς E + ϵ = 1 ς E + ϵ , (20)

Improved modeling of the MeV (ERs) scale is important for searches for neutrinoless double-beta ( 0 ν β β ) decay, for which the key discrimination is not NR vs. ER but between two forms of the latter (β vs. γ). EXO-200 (Anton et al., 2019) and KamLAND-Zen (Abe et al., 2023) have produced the two most stringent half-life limits for ¹³⁶Xe and are highly competitive with the Ge-based experiments. In addition to these results, one must evaluate the prospects of nEXO (Albert et al., 2018), LZ (Akerib et al., 2020c), XENONnT (Aprile et al., 2022), and XLZD (Aalbers et al., 2022) for this field of nuclear physics. Dark-matter-focused experiments have greater ER backgrounds than nEXO but superior energy resolution.

Long-term future work on NEST will involve an ab initio MC approach incorporating cross sections for recombination and the other relevant processes (Piazza et al., 2025), and molecular dynamics modeling of Xe atoms with the 12-6 Lennard-Jones potential for van der Waals forces will be explored (Equation 21). The LXe values for the L-J parameters as well as for other, more advanced versions of the model are well-established (Rutkai et al., 2017):

V d = 4 ϵ L J R d 12 − R d 6 , w h e r e ϵ L J = 1.77 k J / m o l a n d R = 4.10 Å . (21)

While these approaches are challenging at high (MeV) energies, they become more feasible at sub-keV scales, where yields are more uncertain; e.g., for ⁸B, fewer interactions are involved, leading to a more computationally tractable problem.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fdest.2024.1480975/full#supplementary-material

Supplementary material Supplementary material at the source