Finding Dark Matter Faster with Explicit Profile Likelihoods
J. Aalbers B. Pelssers V. C. Antochi P. L. Tan J. Conrad
physics.ins-det astro-ph.CO astro-ph.IM hep-ex



[]jelle.aalbers@fysik.su.se[]bart.pelssers@fysik.su.se[]cristian.antochi@fysik.su.se[]pueh-leng.tan@fysik.su.se[]conrad@fysik.su.se


Liquid xenon time-projection chambers are the world's most sensitive detectors for a wide range of dark matter candidates. We show that the statistical analysis of their data can be improved by replacing detector response Monte Carlo simulations with an equivalent deterministic calculation.
This allows the use of high-dimensional undiscretized models, yielding up to 2 times better discrimination of the dominant backgrounds. In turn, this could significantly extend the physics reach of upcoming experiments such as XENONnT and LZ, and bring forward a potential 5 σ dark matter discovery by over a year.



    95.35.+d,  
    14.80.Ly,  
    29.40.-n,  
    95.55.Vj   



Finding Dark Matter Faster with Explicit Profile Likelihoods
    J. Conrad
    December 30, 2023
============================================================






§ INTRODUCTION




 §.§ Motivation

Astrophysical and cosmological measurements have established that about four-fifths of the mass of the universe's matter consists of dark matter <cit.>. Several experiments attempt to detect dark matter particles on earth <cit.>. Liquid xenon time projection chambers (LXe TPCs) lead the field of direct detection for many dark matter models <cit.>, including the prominent category of models with dark matter mass around 100GeV/c^2 known as weakly interacting massive particles (WIMPs).

Along with other rare-event searches <cit.>, LXe TPCs use (profile) likelihood ratio tests to derive their final physics results, since these tests have (nearly) optimal statistical power <cit.>. In this work, we suggest an improvement to how upcoming LXe TPC experiments such as 
XENONnT <cit.>,
LZ <cit.>,
and PandaX-4T <cit.>
calculate their likelihood. We show this leads to a higher sensitivity through improved signal/background discrimination, and more robust results through supporting more simultaneous correlated nuisance parameters. 




 §.§ Background


Particle physics experiments often specify their detector response model implicitly through a Monte Carlo (MC) simulation. 
However, their inference uses (profile) likelihoods, which need the differential expected event rate – the number of events expected in an infinitesimal volume around the event, in the space of observables such as signal amplitude or reconstructed position.
To obtain differential rates from a simulation, experiments use density estimation: usually, they histogram the  of a large simulation run, e.g. 𝒪(10^7) events. Such histogram-derived densities are known as templates.

Parameters of the detector response model are implemented as configuration options of the simulation. During inference, an optimizer
explores many parameter-value combinations. To avoid re-doing simulations at every point during inference, experiments often precompute templates at strategically chosen points in the parameter space, and interpolate between them; linearly, or using more advanced techniques – this is known as template morphing <cit.>.

To maximize sensitivity, an experiment would like to use all relevant simultaneous observables measured in an event. In particular, unless the detector response is homogeneous and constant in time, reconstructed position and time are relevant. This means templates must be high-dimensional histograms, which require a large number of simulations.

To derive robust results, experiments account for uncertainties as `nuisance' (i.e. additional) parameters of the likelihood, which are profiled during inference (or marginalized, in Bayesian methods).
The more nuisance parameters are used, the more templates must be precomputed. In particular, correlated parameters require templates precomputed on a multidimensional grid, causing an exponential scaling of cost with number of parameters. 

These requirements have an escalating price in terms of computation. Experiments therefore often make compromises that decrease their physics potential, such as discretizing the model in a few spatial- or temporal `bins', or reduce robustness, such as considering only a few key or aggregate uncertainties as nuisance parameters.



 §.§ Concept and outline



In this paper, we show how to compute the differential expected event rate for events in LXe TPCs directly, without MC simulations or templates. The result is not an approximation of an MC. Rather, we directly obtain the result of running an MC simulation with infinite statistics to fill an infinitesimally binned histogram, separately for the detector conditions appropriate for each individual event.

We build on the basic idea that every MC simulation approximates an integral or sum. For example, take a trivial simulation such as:

    N    ∼Poisson(λ) 
    
        S    ∼Gauss(μ = N, σ = 0.1 √(N)) .

This might represent an experiment observing signals of size S, obtained by Gaussian smearing from a number of photons N, which is in turn is sampled from a Poisson distribution with mean λ. 
To compute the probability density P(s) at one or more observations, we could run the simulation many times, populate a histogram, and look up the estimated density at each observed s. Alternatively, we could compute the following sum:

    P(s)    = ∑_n P(s|n) P(n) 
       = ∑_n Gauss(s-n, 0.1 √(n)) Poisson(n | λ) ,

where Poisson(n | λ) = λ^n e^-λ/n! and Gauss(x, σ) = exp(-x^2/(2σ^2)/√(2 πσ^2). This method works even if the mean of the Gaussian distribution is a complicated function of the number of detected photons, e.g. involving an externally determined response map.
The sum runs over all non-negative integers, but evaluating only e.g. n ∈{s-0.5√(s), s-0.5√(s) + 1, ..., s+0.5√(s) } will give accurate results except for ≳ 5 σ outliers.

For the LXe emission and detector response model, the summation is more complicated to construct, and sensible bounds are more difficult to estimate – but the principle is the same. The summation becomes a large matrix/tensor multiplication, implemented using TensorFlow <cit.>, to take advantage of GPU acceleration and automatic differentiation. The latter gives access to the gradient and Hessian matrix in the space of nuisance parameters, which makes inference with a high number of nuisance parameters practical.

We present an open-source Python package  <cit.> which integrates the direct differential rate computation in an inference framework. Instead of specializing to one experiment's choice of LXe model, we provide a framework in which functions such as the parametrizations of yields can be easily customized. 
Our model skeleton is inspired by, though not exactly equivalent to, NEST <cit.> and models used in XENON1T <cit.>.
 includes utilities for statistical inference with the MINUIT and SciPy optimizers <cit.>.
Finally, we show that 's high-dimensional likelihood gives up to a factor 2 better discrimination of the dominant electronic recoil background when compared with a classical 2d MC-template approach, and a corresponding increase in physics reach for dark matter searches.

Section <ref> reviews the principle of LXe TPCs and their likelihoods as they are currently used. Section <ref> discusses  's method and models in detail. Section  <ref> discusses 's performance on a test case that aims to resemble future detectors such as XENONnT and LZ. We end with a short summary and outlook in section <ref>.









§ LXE TPCS AND THEIR LIKELIHOODS



In this section, we briefly recapitulate the operational principle of LXe TPCs, and introduce notation used later in the paper.
A thorough exposition of LXe TPCs can be found in <cit.>.




 §.§ LXe TPCs




TPCs used in modern dark matter searches are meter-sized cylinders filled with liquid xenon, with a few centimeter thick layer of gaseous xenon at the top. Figure <ref> illustrates the working principle of such a LXe TPC.
Particles scatter off either the electrons or the nucleus of a xenon atom, resulting in electronic (ER) or nuclear recoils (NR), respectively. 
Recoils create excited xenon dimers, which decay within 𝒪(10ns), emitting UV light observed as a signal (S1) by the two arrays of photomultipliers (PMTs) at the top and bottom of the detector. Electron/ion pairs are also created in the recoil. The electrons are drifted up towards the gas by a 𝒪(100V/cm) electric field. The electrons then produce a larger secondary signal (S2) in the gaseous xenon under the influence of a 𝒪(10kV/cm) electric field. The light pattern of the S2 signal on the top PMT array indicates the horizontal (x, y) position of the event; the drift time of the electrons (measured by the time delay between S1 and S2) indicates the depth (z) of the interaction. The relative size of the S1 and S2 signals depends on the recoil's energy and type (ER or NR). This allows LXe TPCs to discriminate the NR signals expected for dark matter (and e.g. radiogenic neutron backgrounds) from ERs, primarily from radioactive contaminants in the LXe. 




 §.§ Observables and likelihoods


Likelihood tests require two key inputs: reconstructed quantities of observed events, and expected differential rates for signal and background sources.

LXe TPC likelihoods use a vector of observed quantities, s, for each event that passes the selection criteria. The s usually includes the amplitudes S1 and S2 of the main detector signals in photoelectrons (PE), sometimes the reconstructed event position – x, y, z or r, z, ϕ depending on whether Cartesian or cylindrical coordinates are used – and sometimes the measured absolute event time t. The likelihood is always expressed in directly observable quantities, both in this work and all implementations we are aware of in the field. For instance, the event positions are reconstructed positions, not presumed true positions.

Likelihoods also require the differential event rate R^j(s) for different signal or background sources j. For example, if s = [S1, S2, x, y, z, t], R^j relates to μ^j, the total expected number of events after selections of the source j, as

    μ^j = ρ∫ dS1  dS2  dx   dy   dz   dt 
          R^j(S1, S2, x, y, z, t)
        .

Here ρ is the liquid xenon density in the detector, so that R^j has units of events/(tonne×year×PE^2). During inference, μ^j and R^j depend on several parameters, such as the dark matter cross-section, the rates of different backgrounds, parameters of the LXe charge and light yield response functions, or experimental parameters such as electron lifetime.

To compare models against the observed data, LXe TPCs commonly use an extended unbinned likelihood:

    L = 
        Poisson(N_tot | μ) 
        ∏_i^events ∑_j^sourcesR^j(s_i)/μ
        .

Here, Poisson is the Poisson probability mass function, μ = ∑_j μ^j is the total expected number of events from all sources, and N_tot is the actual observed number of events. Eq. <ref> ensures that the sum calculates the probability density function of events in the total model. In practice, the logarithm of eq. <ref> is used:

    log L = - μ + ∑_i^eventslog( ∑_j^sources R^j(s_i) )  + constant .

The constant depends on the dataset but not the model, and cancels in likelihood ratio inference. It is thus omitted.

The likelihood can have many parameters, but one type is worth highlighting: rate multipliers r^j, which are unitless scales of both μ^j and R^j. That is, 

    μ^j(θ)    = r^j × m^j(θ) 
    
        R^j(θ)    = r^j × M^j(θ),

where m^j(θ) and M^j(θ) are the expected number of events and differential rate, respectively, for the source j if the multiplying parameter r^j = 1, and θ encapsulates the other (nuisance) parameters in the likelihood, known as shape parameters. The r^j = 1 represents a reference level, e.g. a cross-section of 1e-47cm^2 for a WIMP, or a particular nominal level for the internal ER background. 



 §.§ Inhomogeneity and corrections


LXe TPCs do not respond in the same way to signals produced at different positions in the detector or at different times during a science run. The S1 photon detection efficiency is generally highest close to the bottom PMT array, and lowest at the top edge of the liquid volume near the liquid-gas interface (which reflects photons impinging on it at a sufficiently large angle). The survival probability of S2 electrons decreases with depth because impurities in the LXe absorb drifting electrons, as quantified by the mean electron lifetime, τ. The electron lifetime during a science run is often time-dependent. There are also other, usually smaller, inhomogeneities in the response of LXe TPCs, such as variations in the drift and extraction fields, PMT gains, PMT quantum efficiencies, etc.

 operates on the full (S1, S2, x, y, z, t) space of observables, and can thus take full account of these effects. The field currently uses template-based methods, for which, as mentioned above, using six dimensions is impractical, except by coarsely discretizing some dimensions. Instead, collaborations define space-and time dependent corrections for their signals, so they can use lower-dimensional likelihoods.
The corrected signals are called (cS1, cS2), or sometimes (S1_c, S2_c). 
Using a template in (cS1, cS2), rather than (S1, S2), gives a much better signal/background discrimination.

Let G_1(x, y, z, t) and G_2(x, y, z, t) denote the mean expected signal in PE per released photon or electron, respectively, at a given position, S1/G_1 and S2/G_2 are asymptotically unbiased estimates of the number of photons and electrons released at the interaction site. These are usually multiplied by constant scale factors g_1 and g_2 to obtain the corrected signals:

    cS1(S1, x, y, z, t)    = S1g_1/G_1(x, y, z, t)
    cS2(S2, x, y, z, t)    = S2g_2/G_2(x, y, z, t) .

The g_1 and g_2 factors are chosen so that, for a homogeneously distributed source, cS1 = S1 on average inside the fiducial volume, and cS2 = S2 at the liquid-gas interface. To be precise:

    g_1    = 1/V δ T∫ dx   dy   dz   dt   G_1(x, y, z, t) 
    
    g_2    = 1/A δ T∫ dx   dy   dt   G_2(x, y, z = 0, t) ,

with the integrals running over the bounds of the fiducial volume V and its projected (x, y) area A respectively, and δ T the duration of the search.

 uses the full space of observables, (S1, S2, x, y, z, t), rather than just (cS1, cS2). Since cS1 and cS2 are direct functions of these,  can still use efficiencies or other functions parametrized in these corrected variables.

Using the full observable space is preferable for two reasons. First, signal and background intensities and spectra often vary with (x, y, z, t), even if the detector response was homogeneous, and even in the fiducial volume of the detector. Second, signal corrections only compensate for a change in the mean response, not for differences in resolution. For example, the electron detection efficiency decreases with depth in a TPC, worsening the S2 resolution further down the detector. For this reason, top regions of TPCs often show the best signal/background discrimination. If the analysis does not distinguish these from the bottom regions, it loses physics potential.




§ METHOD





 §.§ Overview




Figure <ref> sketches the modeling functionality of . Users can specify the desired physics through model functions – e.g. how the electron lifetime and ER charge yield vary with space and time (and, for the latter, energy). Model functions will be highlighted in green in the text below. Each can depend on an arbitrary number of observables such as (x, y, z, t), and some can additionally depend on hidden variables, such as energy. The model functions are used for three main computations in . 

First, given data (provided as a Pandas DataFrame <cit.>),  estimates reasonable ranges of the hidden variables that are summed over in the computation. A parameter  controls the width of the bounds: increasing it yields more accurate results for outlier events at the cost of computational speed. For this paper, we use = 5, to much the same effect as in the example in section <ref>.

Second,  uses the model functions inside the differential rate computation itself. We first convert the data to a series of TensorFlow tensors, segmented in batches of configurable size to control the memory required for the computation. Next, we compute a series of tensors that are multiplied together to yield the differential rate, as explained in section <ref> - <ref>.

Finally,  also contains an ordinary Monte Carlo code for simulating events using the same model functions. Besides its use in creating toy datasets, we will use it in section <ref> to confirm that 's differential rate computation matches the result of an equivalent template-based method.









 §.§ Mean estimation



The simulation code also plays a role in 's inference. Besides differential rates, the likelihood in eq. <ref> uses the expected total number of events μ. Computing this by integrating the differential rate computation over the space of possible observables would be very inefficient. Instead, we create a toy dataset with a configurable number of events (by default 10^5) and count how many yield detectable signals that pass all efficiencies.

Like differential rates, the μ^j depend on shape parameters θ: see eq. <ref>.  is able to use the same method as classical template morphing: precompute m^j at some particular θ values, then interpolate between the results while θ varies during the inference. This is not a significant computational burden, since a single number is much simpler to estimate and interpolate than a multi-dimensional template. For example, as discussed in <cit.>, conventional (`vertical') interpolation performs poorly in modeling shifts in distributions, even in one dimension. For a single number such problems do not occur.

Interestingly, template-like interpolation of m^j is not needed at all in many cases, or a very coarse estimate suffices. If the error on the estimate of m^j(θ) is well below the width of the external constraint on r^j, fits or limit constructions will converge to the correct values of the shape parameters θ̂, and the inaccuracy in m^j(θ) will be absorbed by the fitted rate multiplier r̂^j alone. After the fit, a single simulation can then be used to estimate m^j(θ̂) accurately, and we find the correct r̂^j by dividing the fitted r̂^j by the ratio of the accurate (post-hoc) and inaccurate (as used during the fit) estimates of m^j(θ̂).

Indeed, LXe TPC likelihoods generally have no or weak external constraints on the rate multipliers r^j of the dark matter signal and the dominant internal ER background – since the dark matter rate is the main parameter of interest, and the most sensitive measurement of the internal ER background usually comes from the data itself. For such a source,  pre-fit estimates of m^j(θ) are not needed at all. For some backgrounds, such as radiogenic neutrons, external constraints on r^j are relevant, since their expected rate is too low to be effectively constrained by the data itself. However, for rare backgrounds, small inaccuracies in m^j(θ) will only have a marginal effect.





 §.§ Model structure




In the remainder of this section, we will describe the ER/NR light- and charge emission and detector response model implemented in . More precisely, we describe the structure in which several physics models can be implemented; the exact form of model functions such as light and charge yield can be chosen by collaborations themselves.  also has an interface for users to change the core structure of the model, e.g. to add extra factors or change assumed distributions fundamentally, which is described in our software documentation <cit.>. As mentioned above, we provide reasonable defaults inspired by NEST <cit.> and our previous work on XENON1T <cit.>. 

Throughout this section, an “event” will be taken to mean a single localized energy deposition in the LXe, causing a single S1 and a single S2 signal. At present  only models single-scatter interactions, or multi-site scatters occurring sufficiently close together that they can be regarded as single scatters for all purposes. These comprise the dominant signals and backgrounds in LXe dark matter searches; resolvable multi-site scatters are efficiently rejected by data quality cuts <cit.>.

Figure <ref> illustrates the implementation of the differential rate calculation in . In outline, we start with an input energy spectrum R_0(E). For each event, this is a vector whose elements range over different energies – see section <ref>. The computation is vectorized over a batch of events, indicated as the depth dimension in figure <ref>.
Since each event yields a stochastic number of detectable quanta n_q, we transform the expected spectrum at each event to a vector R_1(n_q) over n_q – see section <ref>. Some quanta manifest as scintillation photons (n^_), others as drifting electrons (n^_) – see section <ref>. Accounting for this, the differential rate becomes a matrix R_2(n^_, n^_). Some n^_ (n^_) of the produced photons (electrons) survive to be detected, as quantified by a matrix  P(n^_ | n^_) (P(n^_ | n^_)) – see section <ref>. Finally, the detected photons (electrons) determine the S1 (S2) amplitude, as discussed in section <ref>. In brief,  P(S2 | n^_) is assumed to be Gaussian, while for S1, two factors 
P(n^_ | n^_) and P(S1 | n^_) are involved, representing, respectively, the probability of detecting n^_ photoelectrons in the PMTs, then having these produce the measured S1 signal. Note that P(S1 | n^_) (P(S2 | n^_)) is only a vector in n^_ (n^_) since the S1 (S2) value is already observed.
We multiply the P's and R_2 together to get R(S1, S2), the differential expected event rate at the observed event.





 §.§ Input energy spectrum



's starting point is an input energy spectrum, represented as a comb of delta functions. This is a natural input format for spectra specified as a finely (but possibly non-homogeneously) binned histogram.
Formally, flamedisx requires a differential rate R_0, where

    μ_0 = ρ∑_E ∫ dx dy dz dt R_0(E, x, y, z, t)

is the total number of expected events in a hypothetical detector with perfect efficiency. Thus R_0 has units of e.g. events/(tonneyear). There is no d in eq. <ref>; this is absorbed in R_0.

For a homogeneous internal background, R_0 is only a function of energy (and perhaps time); for dark matter, R_0 is a function of energy and time (since dark matter signals are expected to have an annual modulation <cit.>); for external backgrounds, R_0 can be a function of all observables.
Like all model functions, the x, y, z, t have no particular special status, and the user can make R_0 depend on an arbitrary number of observables.
Dependencies on hidden variables – such as E for R_0 – are fixed, since they determine the structure of 's computation.




 §.§ Detectable quanta



The differential rate as a function of the number of quanta is:

    R_1(n_q) = ∑_ P(n_q|) R_0() .

We omitted the dependence on (x, y, z, t), and will do so below, since all model functions can depend on arbitrarily many observables.

The probability to generate n_q quanta (scintillation photons and ionization electrons) given some deposited energy  is given by

    P(n_q|)_ER   = δ(n_q = ⌊/ W⌋) 
    
        P(n_q|)_NR   = Poisson(n_q |ℒ() / W),

for ER and NR, respectively. Here δ is the Kronecker delta function, Poisson the Poisson probability mass function, 𝒲≈13.8eV the LXe work function, i.e. the energy required to generate one detectable ER quantum, and ℒ() the mean fraction of NR energy used for making detectable quanta (rather than lost as heat). This is often called the Lindhard factor, though the user can specify any functional form of ℒ, not just the one proposed by Lindhard <cit.>.





 §.§ Splitting over photons and electrons



The differential rate of producing n^_ photons and n^_ electrons is given by:

    R_2(n^_, n^_) = ∑_n_q P(n^_| n_q) R_1(n_q) 
    δ(n_q=n^_ + n^_).


The detected quanta will split binomially over photons or electrons. Thus n^_ = n_q - n^_, and

    P(n^_| n_q) = Binom(n^_| n_q, p_e),

where Binom is the binomial probability mass function and p_e the fraction of detectable quanta that manifest as electrons. ERs show an overdispersion on top of this, commonly modeled as a Gaussian fluctuation on top of p_e, likely due variations in recombination <cit.>. The Gaussian model cannot be physical, since it allows probabilities outside [0,1]. Instead, we assume the fluctuation is described by a beta distribution:

    P(p_e| n_q) = Beta(n^_|α, β) .

We can express the beta distribution's formal parameters α and β in the mean μ_p_e and standard deviation σ_p_e as

    α = μ_p_e( μ_p_e/σ_p_e^2 - μ_p_e^2/σ_p_e^2 - 1 ),
        
          β= α (1/μ_p_e - 1)

Since the beta distribution is a conjugate distribution of the binomial, the process is equivalent to drawing n^_ directly from a beta-binomial distribution. Thus, our model for n^_ for ERs and NRs is:

    P(n^_| n_q)_ER   = BetaBinom(n^_| n_q, μ_p_e(n_q), σ_p_e(n_q)) 
    
        P(n^_| n_q)_NR   = Binom(n^_| n_q,  p_e(n_q)) .

The parameters are functions of n_q rather than . For ER,  and n_q have an identity relation (see eq. <ref>), so this is merely a change of notation. 
For low-energy NRs, it is an open question whether a direct dependence on  instead of n_q would better fit the data. As low-energy NR tracks are much smaller than the electrons' thermalization length <cit.>, at the time of recombination, there should be no memory of  other than through the number of surviving excitons and electron/ion pairs n_q. Thus a direct dependency on n_q may well be favored. On the other hand, post-thermalization recombination is not the only process involved.

Current models often apply a separate factor for suppression of the NR light yield due to Penning quenching <cit.>. At present,  implements this by unphysically making the photon detection efficiency a function of the number of produced photons; a more elegant implementation would be to adjust the forms and parameters for p_e and the Lindhard factor ℒ.













 §.§ Signal quanta losses



Some of the electrons and photons produce detectable signals in the PMTs, while others are lost. These efficiencies are usually position-dependent and include a variety of effects, such as absorption by impurities and lossy extraction through the liquid-gas interface for electrons, and imperfect collection and quantum efficiency for photons. Each signal photon and electron undergoes these processes individually, without interacting with other photons or electrons. Thus the probability of each photon to be detected (given a particular starting position, time, etc.) must be the same and independent of what happens to the other photon, and likewise for electrons. The fraction of detected photons (and likewise electrons) must therefore follow a binomial distribution.

Separately from the per-quanta efficiencies η^ and η^, experiments have per-event detection efficiencies. For example, two-photon S1s are often excluded from analysis, since these are likely to arise from accidental coincidence of PMT dark counts. These efficiencies are usually parametrized as functions of either the raw S1 or S2 signal amplitude (see the next section) or the number of detected photons and electrons. We call the latter ζ^ and ζ^, respectively.

Given n^_ produced photons (n^_ electrons), the probability of detecting n^_ photons (n^_ electrons) is thus:

    P(n^_| n^_) 
               = ζ^(n^_) Binom(n^_| n^_, η^) 
    
        P(n^_| n^_)
               = ζ^(n^_) Binom(n^_| n^_, η^) .





 §.§ Final signal production
 

For sufficiently small S1s, the non-Gaussian response of PMTs to single photons is relevant. PMTs in LXe TPCs have a probability
p_DPE≈ 0.2
of producing two rather than one photoelectron for each detected LXe scintillation photon <cit.>.  accounts for this with a binomial factor:

    P(n^_| n^_) = Binom(
        n^_ - n^_ |  n^_, p_DPE ) ,

where n^_ is the number of detected photoelectrons in response to n^_ detected photons. The factor is set to zero for n^_ < n^_. This binomial cannot be combined with the one in eq. <ref>, because it is evaluated at n^_ - n^_, not n^_. 

We only apply this factor for S1, because S2s signals are large enough (> 100 photoelectrons) that the non-Gaussian nature of the double-photoelectron emission effect is irrelevant. The double-photoelectron effect affects the interpretation of the Gaussian S2 response parameters (described below), but not the functional form of the response.

The observed S1 and S2 signal size given a number of detected quanta – photoelectrons for S1, and extracted electrons for S2 – is assumed to follow a Gaussian distribution. We thus have:

    P(S1| n^_)
               = ξ_S1Normal(S1|μ_S1(n_^) ,
            σ_S1(n_^)) 
    
        P(S2| n^_)
               = ξ_S2Normal(S2|μ_S2(n^_),  σ_S2( n^_) ),

with μ_S1 and σ_S1 the mean and standard deviation S1 area in PE produced by a single photoelectron, ξ^S1 the per-event detection efficiencies parametrized as function of S1 (and thus not included in ζ^), and similar quantities for electrons and S2s.

Assuming a linear detector response to photoelectrons and extracted electrons, the final signal means and standard deviations are fully specified in terms of those of a single quantum:

    μ_S1(n_^)    = μ_1^S1 n_^
    σ_S1(n_^)    = σ_1^S1√(n_^)
    μ_S2(n_^)    = μ_1^S2 n_^
    σ_S2(n_^)    = σ_1^S2√(n_^)

with μ_1^S1 and σ_1^ the mean and standard deviation of the S1 area produced per single detected photoelectron, respectively, and similar quantities for S2.






§ RESULTS








 §.§ Model verification




A first and essential check is to verify that 's direct differential rate computation (figure <ref>) produces the same result as a lookup in finely-binned template created from a high-statistic simulation. A divergence would point to an error in . 

Computing even a single 6-dimensional (, , x, y, z, t) histogram with sufficient statistics for this check would be a significant undertaking. It is much more practical to fix the simulated event positions and time,
so only a two-dimensional (S1, S2) histogram is needed for comparison. We take the (S1, S2) histogram to have 70 bins in both S1 and S2, uniformly increasing in log S1 and log S2 respectively to better capture the quickly-varying region at low energy. We fill the histogram with 10^8 simulated events, in accordance with (or excess over) current practice.

Figure <ref> compares the  model against the high-statistic template, for two sources: a 0-10keV flat ER background and a 30GeV/c^2 WIMP, positioned at the center of a XENON1T-like detector. We see that  agrees with the template deep into the tails of the distribution, until the high-statistics template becomes unreliable due to statistical errors. The result of the two methods should be exactly the same only in the limit of an infinitely large simulation, a histogram with infinitesimally small bins, and →∞.

It is simple to repeat this test, though it should not be necessary often. Users will frequently change
only the model functions, which are shared between the simulator and the direct computation. Changes in these cannot by themselves cause a divergence if the model is implemented correctly.




 §.§ Single-event discrimination


 can boost the physics potential of LXe detectors by making a full, undiscretized (S1, S2, x, y, z, t) likelihood viable, in a robust analysis with several correlated nuisance parameters. In this section, we will show the benefits of this, by comparing 's likelihood with a two-dimensional (cS1, cS2) likelihood.

A fully specified model is necessary for this comparison. 
We used the ER and NR models described in <cit.> (reparametrized to fit the  structure) and the S1 and S2 efficiency maps from <cit.>. The exact model functions can be inspected along with the the  source code at <cit.>. 
We assume a fiducial volume of 5tonnes, an electron lifetime of 500us, a TPC length of 1.5m, a homogeneous drift field and perfectly stable detector conditions.



Using this model, we compute two-dimensional (cS1, cS2) ER and NR templates, using the same binning discussed above, except that it is linear in cS1 to match existing practice. We fill each template with 10^8 events each, then use the templates as well as 's direct computation to estimate the differential rate for simulated events. The ratio of expected differential rates of two sources is the Neyman-Pearson optimal discrimination test statistic to discriminate single events.

Figure <ref> shows several discrimination efficiencies (“ROC plots”) using this statistic. Specifically, it shows, for different backgrounds, the fraction of remaining expected WIMP signal events if a given fraction of the background must be removed. We see that, at 50% signal acceptance,  effectively reduces the dominant homogeneous flat ER background by a factor two. This is almost entirely due to separating different z regions of the detector, which have different S2 resolution (and thus different discrimination power) due to the finite electron lifetime. The neutron background is discriminated by an even larger factor, due to its spatial dependence within the assumed fiducial volume. The WIMP signal's annual modulation contributes almost nothing to this discrimination: a hypothetical homogeneous NR background with the same mean energy spectrum as the signal (but without modulation) can barely be distinguished from the signal.

We stress that the same results would be obtained in a traditional template-based method with densely binned high-dimensional templates, if it would be feasible.




 §.§ Speed





 has a fundamentally different performance profile than classic template-based likelihoods. For templates, filling the histograms with simulated events is usually the most expensive step: the number of events needed grows exponentially with the histogram dimensionality and the number of potentially correlated parameters that shape the distributions.  does not use templates, and thus requires essentially no precomputation – only 𝒪(10s) of TensorFlow graph construction, which does not need to be repeated when a new (toy) dataset is considered.

Figure <ref> shows an estimate of the required template computation time for a two-dimensional (cS1, cS2) and a six-dimensional (S1, S2, x, y, z, t) likelihood. For the central model (bands) we assumed a simulator capable of producing an accurate 2D-template model in 100 (50-200) seconds
[To see cS1/cS2 discrimination at the one in a thousand level, we need 𝒪(10^6) events per cS1 bin. Having 𝒪(100) cS1 bins thus means we need 𝒪(10^8) events to fill a template accurately. Assuming the simulator can create 𝒪(10^6) events per second – flamedisx's built-in simulator achieves about half this speed on currently common CPUs – implies that 𝒪(100) seconds are needed to construct a 2-dimensional template.]
and 7 (5-10) templates per parameter. For the 6D-template model, we assume 7 (5-10) bins each for (x, y, z and t). Clearly, with more than a handful of nontrivial parameters, even two-dimensional templates become unwieldy, and six-dimensional template grids require computation times measured on geological scales.

After precomputation, template-based likelihoods are relatively fast: computing differential rates only requires a few lookups and simple interpolation.  must instead do a new differential rate computation for each event, which is generally the rate-limiting step.

On a single Tesla K80 GPU, 's differential rate computation runs at 3events/ms for events from a 0-10keV ER model. NR models are nearly three times faster to evaluate, mainly because they lack the beta-binomial in equation <ref>. Thus, the likelihood of a 𝒪(1000event) dataset under a model with both ER and NR sources can be computed in about half a second. 's differential rate computation can also run on a CPU, but then it is 𝒪(10^2) times slower.

During inference, the likelihood must be computed at many different points in parameter space.   exploits TensorFlow's automatic differentiation to compute the gradient and (optionally) the Hessian in parameter space. With this information, high-dimensional optimization needs far fewer iterations, and becomes more robust. Figure <ref> shows the mean duration of fitting an ER model to a 1000-event ER calibration dataset. Specifically, we fitted polynomial models of different orders for μ_p_e(n_q) – effectively controlling the mean ER charge yield as a function of energy.

Figure <ref> is clearly not a level comparison of  and a template-based inference framework: it shows an estimated precomputation time on a CPU for the template method, and a per-fit measured time on a GPU for . We juxtapose the two pieces of information to highlight the fundamentally different scaling with the number of nontrivial parameters. A template-based framework could be ported to GPUs, but the only benefit would be to allow, perhaps, two or three more dimensions or parameters within the same computation time frame.




 §.§ Physics reach improvement




The improved background discrimination translates into an improved physics reach. 

One way to measure the physics reach is by the expected discovery significance, i.e. the p-value of the background-only hypothesis, for one particular assumed WIMP signal.
Discovery significances can be computed using the profile likelihood test statistic:

    t_0 = -2 logL(σ = 0, θ̂̂̂)/L(σ̂, θ̂) ,

where σ is the WIMP cross section, θ represents nuisance parameters (here only the rate of the internal ER background), σ̂, θ̂ denote the global best-fit value, and θ̂̂̂ the best fit of θ conditional on σ = 0. For sufficiently large exposures, as well as certain other technical conditions <cit.>, t_0 is asymptotically distributed as 1/2δ(0) + χ^2_ν=1, with δ the Dirac delta function and χ^2_ν=1 the chi-squared distribution with one degree of freedom.
We verified that the distribution of t_0 is well-described by this approximation for 200GeV/c^2 WIMPs and exposures of 10tonne year and higher, using fits to background-only toy MCs. Specifically, a 3 σ asymptotic significance corresponds to a 3.0 ± 0.2 σ true significance, for both  and the (cS1, cS2) template method.

Figure <ref> shows the expected asymptotic discovery significances for a σ = 2e-47cm^2, 200GeV/c^2 WIMP in different exposures of the XENONnT/LZ like test model discussed above. This model is not excluded by current experiments <cit.>, but well within range of the next generation of detectors <cit.>. We assumed a flat ER background level of 75events/(tonne × year) and a radiogenic NR background of 0.04events/(tonne × year).
Clearly,  gives a substantial advantage over using a classic (cS1, cS2) likelihood. With a 5tonne fiducial volume, it cuts the time needed to get a 5 σ discovery for this model by around a year.
More details on this projection can be found in <cit.>.

We stress that our model should be regarded as a test model only. Accurate projections of the performance of future detectors are more complex and rely on information available only to collaborations themselves. For example, we omitted the accidental coincidence and coherent neutrino nucleus-scattering background completely, and our neutron background spectral shape and its spatial rate dependence are chosen to roughly approximate published figures <cit.>.

The physics reach can also be characterized by the median exclusion sensitivity in case dark matter does not exist. In this case, background discrimination is much less relevant, since the dark matter models probed are by definition at the statistical limit of what the experiment is capable of distinguishing. Only in the extreme background-limited case, where the sensitivity scales with the square root of exposure, do discrimination improvements (and ) yield a proportional improvement of the exclusion sensitivity. LXe DM searches so far have usually stopped running to build a bigger detector well before reaching this limit. 



§ DISCUSSION




 §.§ Validity of the model


As discussed in section <ref>, our model structure is not exactly equivalent to models previously used in the field, but designed to have ample expressive power to fit real data. Since LXe collaborations have not released their raw data to the public, we cannot report on our efforts to verify this here. For future  versions, we hope to implement defaults that more closely track NEST <cit.> for a variety of detector conditions.

We do not suggest experiments use the  defaults of  out of the box, but that they use  to fit their calibration data first. Regardless of which modeling framework is used, assessing its goodness of fit is critical. 's integrated simulator can be used to produce lower-dimensional templates for use in non-parametric tests, e.g. chi-squared or Kolmogorov-Smirnov. If the fit is poor, users can re-fit to calibration data, change the model functions, or even change the model structure fundamentally.

A standard way to address (potential) mismodeling is to add nuisance parameters for uncertain aspects of the model.  helps here, because it can handle more correlated nontrivial uncertainties than a MC-template based framework. Finally,  is a single framework capable of both calibration data fitting and the final scientific inference. Calibration data can also be fitted simultaneously with the final inference, eliminating the sometimes difficult choice of which uncertainties to propagate, and thereby another source of potential mismodeling.

Using a `six-dimensional likelihood' might give some readers pause. Would this mean that orders of magnitude more calibration data must be collected to verify that the model fits? This is not the case, for several reasons.
First, we must not confuse the number of observables with the number of parameters in the model.
, like an MC simulation, can use as many or few parameters in its functions as the user wants, which determines whether it will under- or overfit the data.  looks at more observables (six) than a two-dimensional template for the same number of events, so it is a more sensitive instrument to decide between models. If anything, using  would mean that less calibration data is needed by an experiment.

Moreover, a two-dimensional template is equivalent to a six-dimensional template that assumes the response is constant over four of the dimensions. A model that actually looks at six observables only has to improve on this rather low bar to be superior – e.g. by accounting for simple and easily verified effects, such as the reduction in S2 resolution by electron lifetime. Finally, since  obviates computationally expensive template morphing, it becomes simpler to add additional nuisance parameters representing model uncertainties to the inference. Properly used, this will increase the robustness of results.




 §.§ Impact on real experiments

's advantage in physics reach is entirely due to allowing a full undiscretized (S1, S2, x, y, z, t) likelihood. We compared  against a (cS1, cS2) likelihood above (e.g. figure <ref>), but LXe dark matter search likelihoods often incorporate space and time dimensions already in a limited way. For two reasons, we still believe using  would significantly improve the physics reach of the imminent generation of experiments.

First, the space and time dimensions are very sparsely discretized in current likelihoods. That is, the histogram templates may cover more than two dimensions, but the number of bins in the spatio-temporal dimensions is low.
The likelihood itself is still unbinned, but less accurate than it would be with an undiscretized model.
XENON1T initially used a 2-d (cS1, cS2) likelihood <cit.>. In XENON1T's full science run, a sparse spatial discretization was used, mostly in the r dimension <cit.>. LUX's likelihood was discretized in four z segments and four time segments <cit.>. PandaX-II used 18 time segments <cit.>. LZ's sensitivity projection assumes a 2-d (cS1, cS2) likelihood <cit.>.

Second, we assumed constant detector conditions and homogeneous drift and extraction fields in our example experiment. Neither of these were realized in the last generation of detectors <cit.>. For detectors with variable conditions, space and time dimensions in the likelihood are more important, so the benefit of using  is greater.




§ CONCLUSIONS AND OUTLOOK



We described a new framework for computing LXe TPC likelihoods, , which replaces simulation-based templates with a direct differential rate computation implemented in TensorFlow. This enables high-dimensional undiscretized likelihoods, which will increase the physics potential of LXe detectors. It also enables consideration of more correlated nuisance parameters, leading to more robust results.

We hope  will be useful for experimental collaborations. 
The  source code is released at <cit.> under a permissive open-source license, and will continue to be developed.
As members of XENON, we are particularly looking forward to possible application of  in the upcoming XENONnT experiment <cit.>. Eventually, we hope the methods used in  might inspire similar changes in likelihoods used by other types of particle physics experiments.



 §.§ Acknowledgements

The authors would like to thank Knut D. Morå, Fei Gao, and others in the XENON Collaboration for useful discussions. The authors gratefully acknowledge support from the Knut and Alice Wallenberg Foundation and the Swedish Research Council. 


apsrev4-1







