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.

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.

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: q = 1 1 + κ ϵ Z    λ ,   w i t h   ϵ Z ≈ 1 0 − 3 E , (16) 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.

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).

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 E s 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: 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) 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. 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) 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 is d E / d x = 120 ϵ Z = 120 0.001 E = 120 0.001 E = 3.8 E , (19) which 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), 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) 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 E s, 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.

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.