{'fraction_non_alphanumeric': 0.029866338092304396, 'fraction_numerical': 0.012080465004151823, 'mean_word_length': 4.709435693531121, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A semi-supervised approach to dark matter searches in direct detection data with machine learning

The dark matter sector remains completely unknown. It is therefore crucial to keep an open mind regarding its nature and possible interactions. Focusing on the case of Weakly Interacting Massive Particles, in this work we make this general philosophy more concrete by applying modern machine learning techniques to dark matter direct detection. We do this by encoding and decoding the graphical representation of background events in the XENONnT experiment with a convolutional variational autoencoder. We describe a methodology that utilizes the `anomaly score' derived from the reconstruction loss of the convolutional variational autoencoder as well as a pre-trained standard convolutional neural network, in a semi-supervised fashion. Indeed, we observe that optimum results are obtained only when both unsupervised and supervised anomaly scores are considered together. A data set that has a higher proportion of anomaly score is deemed anomalous and deserves further investigation. Contrary to classical analyses, in principle all information about the events is used, preventing unnecessary information loss. Lastly, we demonstrate the reach of learning-focused anomaly detection in this context by comparing results with classical inference, observing that, if tuned properly, these techniques have the potential to outperform likelihood-based methods.


I. INTRODUCTION
The origin of dark matter (DM) of the Universe is an unresolved mystery, despite the unprecedented experimental effort of the last decades, in particular in the context of Weakly Interacting Massive Particles (WIMPs). One of the most relevant type of searches is looking for nuclear recoils in underground detectors, known as direct detection (DD) [1]. Already now, and in the coming years, several xenon (PandaX [2], XENONnT [3], LZ [4] and ultimately DARWIN [5]) and argon (DEAP-3600 [6], DarkSide and Argo [7], ArDM [8]) ton-scale experiments will improve significantly the sensitivity, see also Refs. [9,10] for the experimental status of DD experiments. In addition, the DAMA annual modulation signal (see Ref. [11] for a review) is being tested in a modelindependent with several experiments using sodium iodine (COSINE [12], ANAIS [13], COSINUS [14] and SABRE North/South [15,16]).
Machine learning (ML) approaches are a growing area of research in the physics community. Indeed, their application in DM phenomenology is starting to become more pronounced. In Ref. [17] a two-phase CNN architecture undertaking first classification and then regression was used to map 3D galaxy distributions to its underlying DM distribution. In Ref. [18] one can see gradient boosted trees used to model DM halo formation.
Similarly, in Ref. [19] a deep CNN architecture is used to model halo formation. In the collider physics context many studies have been published in recent years that attempt to parametrize the concept of anomalousness in data using unsupervised machine learning -see Refs. [20][21][22][23][24][25][26][27][28][29][30][31][32][33]. These studies employ a variety of different methods in a number of different ways; of particular relevance to this paper are Refs. [24] and [33] which employ a variational autoencoder (VAE) on simulated events of the CMS and ATLAS detectors respectively, as well as Ref. [32] where an algorithm dubbed "anomaly awareness" is proposed that utilizes a single supervised CNN architecture to identify anomalous events in boosted hadronic jet structure. In this work we show how using more than one of these such architectures can yield optimal results. In Ref. [34] a convolution neural network was applied to XENON1T TPC detector response images of WIMP signal events (nuclear recoils) and background events (electron recoils) achieving a classification accuracy of above ∼ 95%. In this study, we exploit a pre-trained convolutional neural network with a similar architecture to boost the overall anomaly awareness of our unsupervised model. In anomaly-based unsupervised analyses, the key is to train a model to successfully capture the information contained within a background data set such that data points that don't conform to the properties of the background can be identified as 'anomalous'. In this work we apply this concept to DM direct detection. For this goal we train a convolutional variational autoencoder (CVAE) on a data set of composite images containing the detector response and S1/S2 signals of simulated electron recoil events from the XENONnT DM detector. We then apply the trained CVAE to simulated data sets with a WIMP signal present, so that the reconstruction loss for these events is on average higher, serving as a measure of anomalousness. We then utilize a previously trained convolutional neural network (a slightly tuned version of that from Ref. [35]), capable of classifying background and signal at ∼ 98% accuracy to construct an even more optimized 'anomaly score', observing that its asymptotic properties give the best anomaly identifying power.
While purely unsupervised methods will likely never outperform supervised classification methods they are still powerfully motivated by their ability to identify events which are not background-like in a way that is agnostic to all signal models. In the DM context where the parameter space of possible models is extremely large, the ability to find any anomalous events is helpful in motivating future studies, as well as identifying unknown systematic errors. This paper is structured as follows. In section II we briefly review the ingredients involved in the DM direct detection event rate. In Section III we discuss the XENON experimental TPC. We outline the generation of the data for the analysis in Section IV. Section V provides background on the concept of a variational autoencoder and convolution. In Section VI we present our CVAE architecture and discuss the training. In Section VII we present our analysis pipeline, including an overview of the anomaly score functions utilized for anomaly detection. Here we show the results demonstrating our methods capability of significantly detecting anomalies and forecasting sensitivity, as well as a caparison of projected sensitivity compared to a traditional likelihood based technique. We then conclude in Section VIII.

II. DARK MATTER DIRECT DETECTION
In this section we briefly review the direct detection of DM by measuring the nuclear recoil due to its scattering in underground detectors, see for instance Ref. [36]. We focus on elastic scattering of DM with mass m χ , with a detector consisting of one type of nucleus of mass m A , mass number A and atomic number Z. The differential event rate per unit detector mass is given by where ρ χ is the local WIMP density, f det (v) is the WIMP velocity distribution in the detector rest-frame and v min is the minimum WIMP velocity required to produce a recoil, above the experiment's threshold, of energy E R , with µ χA = m χ m A /(m χ + m A ) being the DM-nucleus reduced mass. For contact interactions, the differential scattering cross section reads where F SI (SD) are the spin-independent (SI) and spindependent (SD) form factors, and σ SI,SD are the corresponding WIMP-nucleus cross sections at zero momentum transfer [37], Here G F is Fermi's constant, f p (n) is the WIMP-proton (neutron) SI coupling, a p (n) is the WIMP-proton (neutron) SD coupling, J the angular momentum of the target nucleus, and S p (n) the spin content of the proton (neutron). Assuming only SI interactions, the differential event rate simplifies to In this paper, we consider SI interactions with xenon nuclei (in the following we drop the SI label). We consider a Standard Halo Model. We will make use of wimprates [38] to generate the DM event rate (see this reference for the values used in this analysis).

III. TPC EVENT RECONSTRUCTION IN THE XENON EXPERIMENT
In this work we simulate background and signal events for the XENONnT experiment located in the Gran Sasso National Laboratory. XENON1T has been running and it is being superseded by XENONnT. The latter uses a dual-phase time projection chamber with 8 tons of natural xenon, of which around 6 tons constitute its fiducial volume. It consists of a liquid xenon target with a gaseous phase on top, with an applied electric field throughout the detector. The detector has two arrays of around 250 photomultipliers (PMTs) in each of the top and bottom layers. It is embedded in a water tank, with muon and neutron vetos [39][40][41].
The main electronic background comes from gamma rays, being generated in the cryostat and in the PMTs, as well as the radioactive decay of 222 Rn and 85 Kr from the detector material. The detection technique makes use of both the prompt scintillation signal or S1, and the scintillation created from the drifted electrons due to the electric field (ionisation signal) or S2. Fiducialisation of the detector is possible thanks to a full 3D position reconstruction by using both the time delay between S1 and S2 signals, as well as the number of photons seen by each PMT (e.g. the hitpatterns).

A. Model independence and supervised classification
The key discrimination parameter to distinguish background (electronic recoils in this context 1 ) from signal (nuclear recoils) is the S2/S1 ratio, which is larger for electron recoils (ER) than for nuclear recoils (NR). For our purposes, it is important to emphasize that this constitutes the main discriminating power between NR and ER, and therefore the classification accuracy achievable by the neural network should be independent of the DM properties (i.e. the DM mass). We note, however, that on an event-by-event basis, the fundamental detector response to a NR vs ER is different, and therefore supervised classification techniques can be used to discriminate between signal and background, as was observed in Ref. [34]. A key result we found whilst validating the results of this study was that such supervised methods work exceptionally well for background rejection regardless of the DM mass and cross-section as well, since the fundamental property being learned by the model is the NR and ER detector response. In this sense, the supervised classifier discriminates between ER and NR completely model-independently 2 . In this study, we focus on the generalized implementation of semi-supervised machine learning for this task, focusing on training the model on ER background-only. For any given data-set, any non-ER event will be labelled as 'anomalous' when passed through the network. It should be noted that in a real experiment there are other sources of NR background (radiogenic neutrons, wall leakage etc.) [34]. We do not consider these in this study since in low background detectors these backgrounds are thoroughly suppressed relative to WIMP signal in the XENONnT experiment [34].

IV. DATA GENERATION
To generate the image event samples used for this analysis, we primarily follow the methodology of Ref. [34] and point readers to their github repository [44] for detailed information about event generation. The ER background spectra are contained within wimprates and are empirically determined using the polynomial fit detailed in Ref. [45]. 3 These spectra are then parsed to laidbax [46], 1 Recently the XENON1T Collaboration reported an excess in electron recoils over known backgrounds at low energies [42,43].
In this work we consider the standard case of a signal produced by WIMPs scattering on nuclei. 2 To clarify, here 'model-independent' refers to the DM model (eg. mass, couplings, etc.). On an event-by-event basis, the supervised classifier discriminates solely among ER and NR, that is, anything that does not classify as ER should then be considered anomalous (non-background) physics. 3 It should be noted that the ER model used in laidbax is not the official model XENON collaboration, which includes several more systematic uncertainties.
which performs the MC generation of events, effectively sampling from the given energy spectrum and producing detector observables for each event. To generate the NR samples that would correspond to WIMP signal, we utilize wimprates [38] to generate energy spectra for an ensemble of WIMP masses in the range 10-1000 GeV. We use a spin-independent interaction cross-section of σ = 10 −45 cm −2 . As discussed in Sec. III A, given that the classification analysis is done on an event-by-event basis, the only relevant property of these images is that they correspond to NR. The coupling strength of the WIMPnucleon interaction as well as the mass of the particle is irrelevant and only serves to determine the overall expected number of events given some exposure. The crosssection therefore simply determines the weights for the Monte Carlo (MC) event sampling. This weighting will effect the projected sensitivity as discussed in Sec. VII A. The physically relevant variables taken for the analysis are recoil type (NR or ER), and corresponding S1, S2 signals, as well several auxiliary parameters. The model output from laidbax is then parsed into PaX (Processor for Analyzing XENON) [47], which generates images of the largest S1 and S2 peaks as well as the hitpatterns for the top and bottom PMTs. Ref. [34] concluded that classification accuracy was maximized when both the S1 and S2 hitpatterns and the largest S1 and S2 peak graphs where displayed in the images. Generating images in this manner meant that only information unique to each WIMP or ER event were present, effectively increasing the robustness of classification. We assume that, irrespective of the architecture used, the composite hitpattern + peak images yield the best classification accuracy. Examples of an ER background event (top) and a WIMP event (bottom) are shown in Fig. 1.

A. Image pre-processing
For the analysis we generate an ensemble of 20000 NR and 20000 ER images. To make training manageable for the architecture presented in this work, the images are scaled down from their original 800 × 800 resolution to 75 × 75 pixels. After validating the results from Ref. [34] with the supervised convolutional neural network, we find that this resolution gives optimal classification accuracy as opposed to 64 × 64 which is often adopted. We therefore also take an input image resolution of 75 × 75 for our CVAE architecture. In Fig. 1 we show examples of an ER and NR event both before (left) and after (right) lowering the image resolution.

V. CONVOLUTIONAL VARIATIONAL AUTOENCODERS
Autoencoders are a type of neural network that attempt to return output that is as similar to their input as possible. This similarity is parameterized by a 'recon- struction loss'. The smaller the loss, the more similar the output is to the input. The aim is to produce an encoding of the input data for the purposes of dimensional reduction, feature learning and produce generative models [48,49].
Autoencoders consist of three main parts: an encoder, a latent space, and a decoder. The data inputs enter into the encoder where they are condensed down into the latent space. This can be thought of as the encoding of the input information to a lower dimensional representation. The decoder then expands this lower dimensional representation of the data back out to the same dimensions as the inputs. For a sample of N inputs, the difference between the input image and the output image is defined as a reconstruction loss, commonly a sum of squared errors, where x is the input vector and y is the output vector. Variational Autoencoders (VAEs) differ from basic autoencoders in that before the latent space vector in an autoencoder we insert two vectors, a vector of means µ and a vector of variances σ 2 of normal distributions N (µ, σ 2 ). Rather than attempting to encode the input vector onto a lower dimensional vector of real numbers, we are now attempting to encode the input information onto a set of Gaussian distributions which when sampled from provide an output vector that is as similar to the input vector as possible. These distributions serve to "cover" the latent space with smooth distributions which is a more faithful representation of smooth and continuous data [50].
Theoretically when a VAE with a standard reconstruction loss from Eq. (7) is trained, it could simply learn Gaussians with extremely small variances (delta functions) or extremely disparate means (poor coverage of the latent space) and thus the VAE would effectively be a standard autoencoder and all benefits would be lost. Thus we regularize the reconstruction loss by penalizing distributions that are different from standard normal distributions. This is done by addition of the Kullback-Leibler divergence to the reconstruction loss, where, for a latent space of K standard normal distributions we have the special case of: which is clearly minimized for σ j = 1 and µ j = 0. The β factor is a parameter of O(1) that scales the overall importance of the KL divergence to the reconstruction loss [51]. The default β = 1 is the definition for a standard VAE, however we find optimal performance for β = 0.5 in this work. A Convolutional Variational Autoencoder (CVAE) is a VAE architecture which utilizes convolutional layers to optimize feature extraction from images and reduce the dimensionality of inputs. For a full discussion of image convolution in the context of machine learning see Ref. [52].

A. Architecture
The CVAE architecture can be seen at the top of Fig. 2 and was implemented in Tensorflow 2 [53]. The network consists of a 75 × 75 × 3 input layer, followed by three convolutional layers each with 128 filters, 3 × 3 kernels and stride lengths of 3. Each convolutional layer uses LeakyRelu activation with α = 0.05, motivated by the hyper-parameter tuning undertaken in Ref [35]. The first two convolutional LeakyRelu layers are followed by pooling layers with dimensions 2×2 and 1×1 respectively. The last convolutional layer is flattened and encoded onto latent space which consists of 512 nodes, corresponding to 256 means and 256 variances, before being decoded via a reflection of the encoder using convolution transpose with "same" padding. The encoder is comprised of three convolutional layers with LeakyRelu activation (yellow), with the first two followed by pooling layers (red). These are then flattened onto a latent space consisting of 256 normal distributions each parameterized by their own means and variances. Samples z ∼ N (µ, σ 2 ) are then drawn from these distributions into a dense layer and then re-broadcast through the decoder using three transposed convolutional layers back into the original image dimensions on which the reconstruction loss (shown in Eq. 8) is calculated between the input and output images (depicted via an arrow). The input and output images shown are actual realizations from our trained CVAE. Bottom: The simple CNN architecture that we used was primarily motivated by the results of Ref. [34] and consists of two convolutional and pooling layers which are then flattened into a single dense output layer with a sigmoid activation function. This is used to calculate the binary cross-entropy in Eq. 11, which will be larger for NR-like images and smaller for ER-like images. Both the reconstruction loss of the CVAE and the binary cross entropy from the CNN are then fed into the anomaly function (Eq. 10), which is then used to do inference.

B. Training and testing on background data
We train the network for 200 epochs on 16000 ER event images in mini-batches of 100 and test on 4000. The training and testing reconstruction loss per epoch can be seen in Fig. 3. Both training and testing appears to converge with similar performance and overtraining does not appear to be an issue.
Once trained, the CVAE should produce generative examples of ER images. This is demonstrated in Fig. 4, where on the left we show a real ER image taken from the testing set and on the right we show the trained CVAE's reconstruction. One can note all primary features of the image are preserved, for example the general shape of the S1/S2 distributions in the bottom half, as well as the colour palate and general topology of the PMT hit patterns.
Before continuing, we briefly draw attention to another well known use for variational autoencoders -that of latent space sampling. As the latent space is a lower dimensional representation of the properties of the input dataset by sampling from the latent space randomly, one may generate 'new' events. This is an on-going and promising area of research, see for example Ref. [31], as it is far more computationally efficient to sample from a latent space of a VAE than to perform costly Monte Carlo. In the case of this study the input data are images, and the pixel information of these images presents an extremely large input space which makes latent space sampling quite difficult. However, using generative models in place of Monte Carlo when the model utilizes raw detector data as inputs is an interesting avenue which we leave for future research.

A. Anomaly detection
The analysis pipeline is summarized in Fig. 2. Anomaly detection refers to the practice of calculating a well defined 'anomaly score' in order to determine if an event or group of events are significantly discrepant from elements of the background distribution [54]. In statistical language, this process involves studying the asymptotic properties of the anomaly function T S(x, y), which is generally a function of the image data inputs x and CVAE outputs y. After data taking, one would calculate T S for all data in the run. The resulting distribution can then be inspected for excesses or 'anomalies' with respect to the simulated background-only distribution.
The function we adopt to perform the anomaly detection analyses is the following function that returns larger values when it is fed data that is more signal-like: The term in parentheses is the CVAE loss function defined in Eq. 8 and the last term is a binary cross-entropy between a vector of zeros and the output of the pretrained CNN. The CNN was trained to classify ER background as 0 and NR as 1 (see Ref. [34] and Fig. 2 for architecture details). Hence, H B is larger in the presence of signal in the data set, e.g. when p(x) → 1: where p(x) is understood as the output of the classifier given the input image x i 4 . The parameter R in Eq. (10) scales the contribution of the cross-entropy term. In the limit R → 0, we recover the reconstruction loss from Eq. (8) for a purely unsupervised analysis, while for R → ∞ we are effectively doing an supervised study 5 . We note that it may seem natural to adopt the reconstruction loss as the anomaly function. However, we found that by generalizing it in this way, we utilize the power of a pre-trained CNN in a semi-supervised fashion which dramatically improves the anomaly awareness, as will be discussed in Sec. VII D.

B. Results for pure signal/background data
We apply the now trained CVAE model to test data sets (data sets the networks were not trained on) consisting of 10000 ER and NR (WIMP) events. To assess the discrimination power of our technique, we histogram the anomaly function defined in Eq. (10) for pure signal injection (NR) vs background-only (ER) (all 10000 events). In other words, we bin TS without re-weighting the distributions based on cross-section and WIMP mass. We do this for two cases: R = 0 (bottom) and R = 170 (top), the latter value being observed to provide optimal results, as will be discussed in Sec. VII D. The results are normalized to unit area and are shown in Fig 5. It can be seen that the distribution of the WIMP events has more density at larger TS on average than the ER events for both values of R, implying that the model is learning something fundamental about the electron recoil events that is not present in the WIMP events. It is directly obvious however, that the case of R = 0 (purely unsupervised) yields significantly worse anomaly awareness than the case of R = 170 (semi-supervised).

C. Results for background plus signal data
The main analysis involves applying the network to a sample of "pseudo-data", created by injecting the WIMP signal NR events into a dataset of ER events scaled to the appropriate normalization. The ER and WIMP events are weighted to an exposure of 5 t y. Furthermore the WIMP signal is weighted to a spin-independent scattering cross-section of σ SI = 10 −45 cm 2 for masses 10, 500 and 1000 GeV. In re-weighting the distributions in such a fashion we are effectively assuming we have used the Asimov data [55]. One may note that the cross-section used is currently excluded for some DM mass values [56]. For the time being, we are simply demonstrating the novelty of the procedure which is exemplified for a relatively large cross-section. An analysis considering non-excluded regions of the DM parameter space is considered later in Sec. VII E.
We expect to see a pseudo-data distribution that matches the electron recoil distribution to a high degree in the low TS areas with excesses seen in tails of the dis-FIG. 7: p-value to reject the background-only (ER) hypothesis as a function of the scale parameter R using an exposure of 5 t y and a spin-independent scattering cross-section of σ = 10 −45 cm 2 . R → 0 implies TS is calculated only using the CVAE (purely unsupervised), while R → ∞ implies TS is dominated by the output of the CNN (purely supervised). The optimal value for R is observed to be in the range R ∈ [150 − 300], indicating the need for a semi-supervised approach for optimal anomaly awareness. Calculation of the p-value is detailed in Sec. VII E.
tributions. This is presented in Fig. 6 where we show the background-only as well as pseudo-data TS distributions for the three WIMP mass benchmarks. The excess is least visible for 10 GeV since this mass gives the lowest number of expected events (13) compared to 500 GeV (187) and 1000 GeV (95). The significance of the excess as a function of exposure as well as a comparison to the XENONnT forecasted 5σ discovery is presented in the next section.

D. Optimizing anomaly awareness
The probability with which we can reject the background-only hypothesis was calculated with a simple Pearson χ 2 test: We assume Wilk's theorem holds such that χ 2 p follows a χ 2 -distribution with one degree of freedom. We expect this to hold given inspection of Fig. 5, clearly displaying asymptotic forms of a nested hypothesis, with the null hypothesis (background-only) being realized explicitly when R = 0. The p-value to reject the null (background-only) hypothesis is then given by where CDF is the cumulative density function of a χ 2distribution with one degree of freedom.
In Fig. 7 we calculated the p-value for rejecting the background only hypothesis while scanning over values of R. We observe that for any given DM mass, a purely supervised anomaly function (large R) has better discriminating power than a purely unsupervised one (R = 0), as expected, with R 10 and R 10 3 yielding a constant p-value. Interestingly, however, we find that there is a range for R centred at about R ∼ 170 for the cross section used, where the p-value is minimized. This appears to be the case regardless of the DM mass (as long as they are above tens of GeV). This does not necessarily apply to varying the cross-section, however, as we will discuss in Sec. VII E. Having the best background hypothesis rejection occur before the plateau at large R indicates the need for both the CVAE and CNN to be utilized in the anomaly function. This is due to the fact that in a purely supervised scenario, constructing the anomaly test in this way would yield significant loss of information due to (nearly) all ER events being binned at T S = 0 and NR events being binned at some arbitrary, large T S value. The inclusion of the CVAE loss into the anomaly function serves to provide a distribution that captures all available information, effectively getting the best of both unsupervised and supervised approaches, qualitatively increasing the statistical power of any observed excess in T S. This demonstrates the overall improvement of our analysis by implementing a semi-supervised approach.
An additional thing to note is that, whilst we see that the anomaly function in Eq. (10) gives good results, one is technically free to choose any such function, as long as the same operation has been consistently done to the background-only and pseudo-data sets. A comprehensive study of the different anomaly functions and their power in the context of semi-supervised learning in different beyond the standard model (BSM) searches is outside the scope of this work and is left for future study.

E. Forecasting sensitivity
When the data being used are a representation of how they would look after some time, the p-value in Eq. (13) can be referred to as the 'projected sensitivity'. We note that a more formal statistical inference may be done, however, to demonstrate the novelty and generalizability of this method, we decided to keep things simple. It is often useful to project or 'forecast' what the p-value or 'sensitivity' of discovery would be (assuming the WIMP exists) after some exposure. In Fig. 8 we plot the p-value calculated as in Eq. 13 as a function of exposure for three benchmark masses and a fixed value of R = 170. For demonstration, we take a currently non-excluded cross section of σ = 10 −45 cm 2 (top) as well as σ = 10 −46 cm 2 (bottom). We additionally include dashed lines representing p-values corresponding to a number of unit Gaussian standard deviations of 2σ and 5σ. We can see that, as expected, the relatively low number of events for a 10 GeV makes it difficult to significantly detect an excess over the background-only distribution, however 500 GeV and 1000 GeV remain reasonably detectable.
Finally, we plot the 5σ sensitivity contour in the crosssection/exposure plane for a variety of R values in the range R ∈ [100, 250]. We present the results as a band to illustrate how the sensitivity can maximally and minimally vary depending on the choice of R. In order to observe the strongest sensitivity, we do this for a 50 GeV WIMP which is the mass which approximately gives the greatest expected number of events for any given crosssection. We show the result in Fig. 9. The variation of sensitivity given a choice of the R parameter exemplifies the need to use a range of R values (or indeed, anomaly functions in general) in a real analysis, since the R value did not enter the training process whatsoever and we are free to choose the one which optimizes anomaly detection (rejecting the background-only hypothesis). In order to qualitatively gauge the performance of our semisupervised anomaly detector relative to a traditional statistical approach, we estimate the 5σ anomaly sensitivity using a traditional frequentist approach by bootstrapping samples of the corrected S1/S2 (cS1/cS2) distribu- tion from laidbax in order to calculate the asymptotic form of a -2ln(Poisson) distribution corresponding to the background-only and background+anomaly hypotheses (see Ref. [55] for details on the standard procedure). We present the forecasted p−value to reject the backgroundonly hypothesis at 5σ in Fig. 9 as a black contour. We observe that our semi-supervised anomaly detector indeed performs better than the likelihood based method by about almost an order of magnitude. We expect this is the case since for classical statistical tests, the choice of likelihood function is paramount in determining the strength of the statistical test. Hence, if one could find the 'correct' likelihood function that exactly describes how nature interacts with the experiment, this would be the optimal choice. Of course, usually we do not have access to such a function in practice, and so we usually take an educated guess (Poisson likelihoods in particle physics for example). What our machine learning method has essentially done is, based on simulations that we consider accurate representations of the physics, 'learn' the underlying properties of the data that would give the 'best estimate' of the likelihood function, even if such an analytical function doesn't exist. This is the general premise of simulation based inference (SBI) (see Ref. [57] for a review). A more detailed study of the general improvement to modern SBI and anomaly detection methods by utilizing multiple, separately trained supervised/unsupervised architectures is beyond the scope and is left for future work.
For completeness, as a dashed black line in Fig. 9 we also include the official XENONnT projected 5σ discovery sensitivity from Ref. [39], which denotes the smallest cross-sections at which the experiment would have a 50% chance of observing an excess with significance greater than 5σ corresponding to a WIMP with mass 50 GeV. It should be noted that this contour should not be directly compared to our analyses, since it explicitly refers to the discovery potential of a 50 GeV WIMP, which comprehensively involves the profiling and minimization over model parameters in the likelihood analysis, as well as the inclusion of several complicated detector dependent terms. The result is that such a sensitivity will be much stronger than that of any anomaly search. It is also important to note that the CVAE was trained on background events only. That is, any events that do not look like background will show up as having an anomalous TS distribution. Whilst the projected sensitivities we present in this section implicitly depend on the WIMP model (mass/cross-section) as we needed to re-weight the T S distributions, the general principle is completely model independent in that one may always run a real data set through the CVAE-anomaly detector to search for anomalous physics of any kind. We thus encourage these methods to be adapted to other low background searches for new physics, especially with detectors for which more than one DM model may be probed. Regardless, it is qualitatively useful to see how the semi-supervised learning approach compares to statistical approaches, even though it is given no direct information about the WIMP model parameters. With further optimization, as well as access to other auxiliary data, it may indeed be very possible to train such an anomaly detector to surpass the currently used state of the art statistical forecasting methods.

VIII. CONCLUSION
In this study we present an example of a modelindependent search for anomalous DM direct detection events using a convolutional variational autoencoder (CVAE). The XENONnT experiment was used as a testbed for this methodology wherein a set of composite images including the detector response and S1/S2 channels were used as input data. The CVAE was trained on a set of electron recoil background events, learning to encode the input images onto a lower dimensional latent space, then decoding from this latent space back to the original image. The reconstruction loss is defined as the difference between the pixel information of the original image and the decoded image. The loss is thus a proxy for anomalousness as images that do not conform to the general behaviour and distribution of background events will inherently produce higher reconstruction losses as the net-work has not learned the properties of these events.
Additionally, in this work we exploit a pre-trained supervised convolution neural network which boosts the anomaly detection power of the CVAE. Interestingly, this means that the best performance for anomaly detection is found for a combination of both unsupervised and supervised networks, in a semi-supervised fashion. We parameterize the influence of the supervised part with a hyperparameter R. We find that the optimal R-value to use is invariant with respect to the DM mass, however, we emphasize that in a real analysis, detailed consideration as well as optimization of any hyperparameters (e.g. R values) should be undertaken.
With respect to future sensitivity, we find that the projected 5σ anomaly sensitivity corresponding to a WIMP mass of 50 GeV outperforms that of a simple classical likelihood based approach. We also find that our anomaly detector performs similarly, although as expected, slightly weaker (e.g. ∼ a factor of a few), than the official prediction of the XENONnT Collaboration for the discovery of a 50 GeV WIMP. Even though these results are not exactly statistically comparable, this is an encouraging result which may render interesting for the Collaboration to run this type of analysis given their access to additional data and more sophisticated experimental techniques.
Finally, we would like to point out that, although a given particle physics interaction was used (point-like, spin-independent, elastic) to train the supervised network, the analysis is quite agnostic regarding the underlying DM model which is generating the excess of anomalous events, whose main characteristic is consisting of nuclear recoils. Given the absence of any positive signals in the last decades, trying to be as agnostic as possible regarding which signals we expect to detect seems the right approach.