Digging for dark matter: Spectral analysis and discovery potential of paleo-detectors

Paleo-detectors are a recently proposed method for the direct detection of dark matter (DM). In such detectors, one would search for the persistent damage features left by DM – nucleus interactions in ancient minerals. Initial sensitivity projections have shown that paleo-detectors could probe much of the remaining weakly interacting massive particle (WIMP) parameter space. In this paper, we improve upon the cut-and-count approach previously used to estimate the sensitivity by performing a full spectral analysis of the background-and DM-induced signal spectra. We consider two scenarios for the systematic errors on the background spectra: (i) systematic errors on the normalization only, and (ii) systematic errors on the shape of the backgrounds. We find that the projected sensitivity is rather robust to imperfect knowledge of the backgrounds. Finally, we study how well the parameters of the true WIMP model could be reconstructed in the hypothetical case of a WIMP discovery.

Progress in the search for DM heavier than ∼ 10 GeV relies on maximizing the exposure (the product of target mass and integration time) of the experiment.Instead, to probe ever lighter DM, experiments must achieve sensitivity to lower and lower nuclear recoil energies.In addition, both mass regimes require exquisite control of a variety of possible background sources, from cosmic rays to intrinsic radioactivity.A number of experiments with lower energy thresholds, larger exposures or, ideally, both will continue the search for lighter and more weakly-interacting DM in the next 5-10 years [32][33][34][35][36][37].
Recently, Refs.[38,39] proposed paleo-detectors for the direct detection of DM 1 .In paleo-detectors, one would examine ancient minerals extracted from O(10) km below the surface of the Earth for traces of DM interactions with atomic nuclei.In lieu of the multi-ton target masses used in conventional direct detection experiments, paleodetectors take advantage of the fact that DM may have been interacting with the target material for as long as ∼ 1 Gyr.In certain minerals, including those long used as Solid State Track Detectors (SSTDs) [40][41][42][43], DMinduced nuclear recoils would give rise to 1 − 500 nm long damage tracks.In many materials, such damage tracks would be preserved over timescales much longer than 1 Gyr.Paleo-detectors could thus obtain much larger exposures than conventional detectors even if only a small amount of target material can be investigated.For example, reading out 100 g of material which has been recording DM-induced events for 1 Gyr provides an exposure which could only be matched in the laboratory by observing 10 4 tons of target mass for 10 yr.Further, the relatively small target masses required for paleo-detectors can be obtained from depths much greater than those of underground laboratories in which conventional direct detection experiments are usually operated, providing an unprecedented level of shielding from cosmic rays.
In conventional direct detection experiments, nuclear recoils are detected by scintillation, ionization, and phonon signals in the detector [44].In paleo-detectors, the observational signature would be nano-scale defects in the minerals.These may be observed through a variety of read-out methods such as X-ray microscopy [45][46][47] or ion-beam microscopy [48].Mineral-based searches, initially for monopoles and then for DM, have been proposed and performed before [49][50][51][52][53][54][55][56].However, modern high-resolution imaging techniques [57][58][59][60][61][62][63] as well as the availability of rocks from deeper underground may significantly improve the prospects for DM detection.In particular, measurements of nanometre-length tracks could provide sensitivity to recoil energies as low as 100 eV.A more detailed study of backgrounds, target minerals and read-out methods for paleo-detectors was recently presented in Ref. [39].
In this work, we explore the prospects for excluding or discovering DM with paleo-detectors.In Refs.[38,39], a simple cut-and-count analysis with a sliding signal window in track length was employed to estimate the sensitivity of paleo-detectors.Here, we adopt more sophisticated statistical techniques, making use of information contained in the spectral shape of the track length distributions.Using realistic distributions for backgrounds induced by neutrinos [64] and radioactivity (developed in Refs.[38,39]), this approach allows us to explore the impact of different systematic background uncertainties on projected sensitivities.Further, we examine how well the DM properties (mass and cross section) could be reconstructed by paleo-detectors in the case of a discovery.
As in Refs.[38,39] we consider two fiducial read-out scenarios.The first assumes a high track-length resolution but a relatively small exposure (referred to as high resolution in the following).The second assumes worse track-length resolution, which should in principle allow more material to be analyzed, facilitating a larger exposure (referred to as high exposure).The high resolution configuration is particularly well suited to probing DM with masses below ∼ 10 GeV while the high exposure configuration is geared more towards heavier DM.
A wide range of minerals are well suited to be paleodetectors.As described in the discussion of mineral optimization in Ref. [39], minerals can be broadly divided into classes suitable for different applications of paleo-detectors, based on their chemical composition.We consider 4 different minerals in this work, chosen to represent paleo-detectors suitable for probing spinindependent DM-nucleus interactions: ).These particular target materials are also selected for their low levels of natural radioactive contamination, helping to suppress radioactivity-induced backgrounds.Average uranium concentrations in the Earth's crust are a few parts per million (ppm) by weight, which would lead to unacceptably high levels of background due to radioactivity.Minerals formed in ultra-basic rock deposits, stemming from the Earth's mantle, are much more radiopure.Examples of such minerals investigated here are olivine and nchwaningite, for which we assume uranium concentrations of 0.1 parts per billion.Even less contaminated by radioactive elements may be minerals in marine evaporite deposits formed at the bottom of evaporating oceans.We assume uranium concentrations of 0.01 parts per billion for such minerals and use halite and sinjarite as examples.Halite and olivine are very common minerals.In contrast, sinjarite and nchwaningite are less abundant but contain hydrogen.Although hydrogen makes up only a small fraction of these minerals by mass, its presence plays an important role in reducing neutron-induced backgrounds, as we discuss later.
The rest of this paper is organized as follows.In Sec.II, we discuss paleo-detectors in more detail, including the calculation of signal spectra, track lengths, and the most relevant background sources.In Sec.III, we present the projected upper limits and discovery reach for the minerals listed above.In Sec.IV, we use benchmark-free techniques to determine contours in the (DM mass)-(cross section) plane which could be reconstructed with paleodetectors in the hypothetical case of a future discovery.We discuss a number of challenges for paleo-detectors in Sec.V. Our conclusions are presented in Sec.VI.
Code for performing all calculations presented in this paper is publicly available here [65].

A. Signal from WIMP Scattering
For elastic scattering of DM with mass m χ off nuclei with mass m N , the differential recoil rate as a function of recoil energy E R per unit target mass is given by [66]: The integral is over DM velocities v, with v = |v| and We assume standard spinindependent (SI) interactions, with equal couplings to protons and neutrons.In this case, the differential cross section can be re-written in terms of the DM-nucleon cross section at zero momentum transfer σ SI n as: Here, µ χN ≡ m χ m N /(m χ + m N ) is the reduced mass of the DM-nucleus system (and similarly for the DMproton reduced mass µ χp ).The factor of A2 corresponds to the coherent enhancement for a nucleus composed of A nucleons.The internal structure of the nucleus is encoded in the form factor F 2 (E R ), for which we assume the Helm parametrization [67][68][69].The differential recoil rate then takes the standard form: For the DM distribution, we assume the Standard Halo Model (SHM), fixing a benchmark value for the local density of ρ χ = 0.3 GeV cm −3 [70,71] in order to compare directly with other direct detection experiments.However, we note that observational estimates of ρ χ vary substantially [72].In the SHM, the DM velocity distribution f (v) follows a truncated Maxwell-Boltzmann distribution, for which we fix values for the Sun's speed v = 248 km/s [73], the local circular speed v c = 235 km/s [74], and the Galactic escape speed v esc = 550 km/s [75].We do not consider here uncertainties on the speed distribution [71] or more recently suggested refinements to the SHM [76].
Recoil spectra for DM and neutrino scattering are calculated using the publicly available WIMpy code [77].

B. Paleo-Detector Rates
The rate of tracks produced with length x T per unit target mass is given by: where the index i runs over the different target nuclei which make up the mineral.The rate of nuclear recoils (with initial energy E R ) per unit mass is given by dR i /dE R and we weight by the mass fraction ξ i of each nucleus i.The track length as a function of initial recoil energy is calculated as: The stopping power dE/dx T as a function of energy must be calculated for each of the recoiling nuclei in a given target material.We use the publicly available SRIM package (Stopping and Range of Ions in Matter) [78,79] to tabulate the stopping power, although analytic estimates are also possible [39,80].A more detailed discussion of the calculation of track lengths can be found in Ref. [39].
In this paper we assume that recoiling hydrogen nuclei and α-particles do not produce tracks which can be reconstructed.Whether such low-Z tracks are observable will depend on the target material and read-out method and is a question which must be determined empirically.A discussion of this issue as well as a comparison of results with and without low-Z tracks can be found in Ref. [39].
The resolution at which track lengths can be measured depends on the read-out technique used.For a true track length of x , we assume that the measured track length x is Gaussian distributed 2 with track length resolution σ x T : The number of tracks with lengths in the range [x a , x b ] is then: where is the exposure, given by the product of the age of the mineral and the total mass of the sample analyzed.
The window function W captures resolution effects and is given by: We assume that the smallest measurable track length is σ x T /2 and consider tracks as long as 1000 nm.As we will see in Sec.II E, DM-induced recoils do not typically induce tracks longer than this.
We will consider two scenarios for the analysis of paleodetectors, as in Ref. [39]: • High resolution -we assume a track length resolution of σ x T = 1 nm, which may be achievable with helium ion beam microscopy [48] (using a focused ion beam [61,81] and/or pulsed lasers [82][83][84] to remove layers of material which have already been imaged).In this case, we assume an exposure of = 0.01 kg Myr, which would correspond to a few O(1 mm 3 ) mineral samples, each with an age of 1 Gyr.
• High exposure -we assume a track length resolution of σ x T = 15 nm.Small angle X-ray scattering has been demonstrated to achieve such spatial resultions in three dimensions [47,62,63].However, such resolutions have not yet been demonstrated when imaging damage tracks arising from nuclear recoils.Here we assume an exposure of = 100 kg Myr, corresponding to the analysis of larger samples of O(10 cm3 ).This is not an exhaustive list of possible scenarios, see Ref. [39] for a discussion of a variety of read-out techniques.
We note that a number of techniques (both stratigraphic and radiometric) are used for dating rock samples [85].Perhaps most relevant for O(Gyr)-aged rocks is fission track dating, which should allow an age estimate which is accurate to ∼ 10% [86,87], though we will neglect dating uncertainties in our analysis.
Given the target materials we analyze, we note that DM candidates with m χ 500 MeV do not give rise to a significant number of recoil tracks longer than ∼ 1 nm, the best track length resolution assumed in our analysis.Thus, we only consider DM with m χ 500 MeV, though it may be possible to probe lower mass DM with either better track length resolution or with target materials which allow for longer tracks.

C. Backgrounds
Here, we summarize the most problematic backgrounds for DM searches with ancient minerals.While follow-up studies and direct calibration will likely lead to refined background modeling, we expect the estimates presented here to be representative.A more detailed discussion can be found in Ref. [39].
We note that cosmogenic backgrounds should be negligible for materials obtained from a depth below ∼ 5 km, meaning that the dominant backgrounds will be neutrino interactions and intrinsic radioactive backgrounds in the target materials themselves.
a. Neutrinos -Being weakly interacting particles, neutrinos represent an irreducible background to (nondirectional) DM searches [88][89][90][91].Neutrinos with MeVenergies can produce keV-scale nuclear recoils, thereby mimicking DM signatures.We calculate the neutrinonucleus scattering rate following Ref.[91] (and references therein).In our analysis, we include solar, atmospheric and diffuse supernova background (DSNB) neutrino fluxes, with spectra and normalizations as compiled in Ref. [64].Uncertainties on the present day neutrino fluxes vary substantially, from O(1 %) for 8 B and hep neutrinos [92] to as much as 50 % for DSNB neutrinos [93].Paleo-detectors probe neutrino fluxes over the past O(1) Gyr, which may differ from the current values.We therefore conservatively assume a Gaussian systematic uncertainty of 100 % on the normalization of each neutrino component independently 3 , including each of the components from different nuclear processes in the Sun.
b. Backgrounds from α-decays -One possible background is from the 'uranium series' of uranium-238, a decay chain which proceeds via a series of α and β decays.With each α-decay in the series, the child nucleus recoils against the α particle with O(10 − 100) keV energy.The half-life of 238 U is T 1/2 ∼ 4.5 × 10 9 yr, while the subsequent decays occur much more quickly (T 1/2 2.5 × 10 5 yr).Thus, the vast majority of 238 U nuclei which decay will have completed the entire decay chain over the age of the mineral (see Ref. [39] for a more in-depth discussion).Even if the α tracks are not observable, the numerous decays in the chain will lead to a characteristic pattern of tracks.We assume that all such track arrangements can be rejected as background.However, we note that in a 10 mg sample of sinjarite there would be O(10 7 ) completed decay chains and further work is required to estimate whether such large rejection factors will be achievable in a real experiment.
A more problematic background comes from uranium-238 nuclei which have only undergone a single α decay ( 238 U → 234 Th + α).In this case, the thorium-234 child has a characteristic recoil energy of 72 keV and (assuming that the α track is not seen) is indistinguishable from a DM-induced recoil.The number of such 1α-thorium tracks depends on the relative half-lives of the 238 U and 234 U decays [( 234 U → 230 Th + α) is the second α-decay in the uranium-238 decay chain] and is roughly where C 238 is the uranium-238 concentration (by weight).
It should be possible to estimate the normalization of radioactive backgrounds to a high precision, for example by measuring the number of full 238 U decay chains in the sample.We therefore assume a 1 % systematic uncertainty on the normalization of the 1α-thorium background, though as we will see in Sec.III, the projected sensitivity of paleo-detectors is limited more by uncertainties in background shapes.
c. Neutron-Induced Backgrounds -Fast neutrons produced in or around the target minerals will scatter elastically with nuclei.Such neutrons have a mean free path of a few centimetres and typically give rise to 10 − 100 nuclear recoil events with energies comparable to those caused by DM.This large number of neutroninduced tracks is therefore difficult to reject as background.
Fast neutrons may be produced in the spontaneous fission (SF) of 238 U.This accounts for roughly 1 in every 2 × 10 6 decays of 238 U, producing ∼ 2 fast neutrons with MeV energies per SF decay.Neutrons may also be produced in (α, n) reactions, in which nuclei absorb an incident α particles and emit fast neutrons.The neutron-induced recoil spectra are estimated using the SOURCES-4A code [94], including both SF and (α, n) contributions, as described in Ref. [39].As in the case of 1αthorium backgrounds, the normalization of the neutroninduced background scales with uranium-238 contami-nation.For minerals found in ultra-basic rocks (nchwaningite and olivine), we assume a uranium-238 contamination of C 238 = 0.1 ppb by weight, while for those found in marine evaporites (halite and sinjarite) we assume C 238 = 0.01 ppb.
We assume a 1 % systematic uncertainty on the normalization of the neutron-induced backgrounds.We discuss the sensitivity of our results to the assumed background uncertainty and explore more extended shape systematics in Sec.III.

D. Analysis Theory
To estimate projected upper limits and discovery reaches for paleo-detectors, we use the statistical techniques developed in Ref. [95].These extend traditional Fisher forecasting methods to the Poisson regime, which allows one to approximate the median result obtained via a Monte-Carlo simulation with minimal computational overhead.We briefly summarize the technique below.
Traditional Fisher forecasting methods are sufficient to accurately calculate expected exclusion limits as well as the discovery reach when in the Gaussian regime.Typically, in direct detection experiments, however, the small number of signal and background events means that the number of expected counts is not well described by a Gaussian distribution.To remain accurate in this Poissonian regime, we adopt the equivalent counts method, as developed in Ref. [95].The basic procedure for the equivalent counts method is to define a mapping between the expected background counts, their associated uncertainty, and the expected signal such that the full profile log-likelihood is approximated by the Poisson loglikelihood ratio.This mapping is given by, where D A (S 0 ) is the Asimov data set [96] given no expected signal events and S is the expected number of signal events.P (a|a ) represents the Poisson probability mass function, i.e., the probability of a events given a expected events.The equivalent signal and background events are denoted by s eq and b eq , respectively.They are defined such that the log-likelihood ratio of a simple onebin Poisson process approximates the full log-likelihood ratio.We found expressions for s eq and b eq in terms of the Fisher matrix of the full problem, which are given in Eq. ( 6) and Eq. ( 7) of Ref. [95].The procedure leadsper definition -to exact results in the limits where the full problem is a one-bin Poisson process or Gaussian, and approximates very well Monte Carlo results in the intermediate range.
The discovery reach and exclusion limit can then, trivially, be calculated by solving − 2 ln P (s eq + b eq |b eq ) P (s eq + b eq |s eq + b eq ) = Z 2 , ( and respectively, and mapping this back on to the signal parameters of the full model.The significance level α determines the critical value Z, which is derived from the inverse of the standard normal cumulative distribution, denoted F N , as For example, Z(0.05) = 1.64 represents a 95% confidence level.Reference [97] showed that for most cases the equivalent counts method is accurate to the percent level.
In some exceptional cases the derived upper limit exhibits a maximum deviation of 40% from a fully coverage corrected Monte Carlo calculation, for more detailed discussion see Ref. [97].These deviations are most relevant when two distinct regions of the spectrum are present; one in which the signal dominates and the background is in the Poissonian regime, and the other in which the background dominates over the signal.In almost all our cases we expect a significant background so this is unlikely to be a problem.The information flux [95] provides an intuitive illustration of which region of signal space provide most information about the DM-induced signal.It is a generalized signal-to-noise ratio which allows for the inclusion of extended systematics and covariances.Technically, the information flux is obtained by taking functional derivatives of the Fisher information matrix with respect to the exposure at different track lengths.It can be thought of as the rate at which the error bar on the parameter of interest will be reduced from an infinitesimal increase in exposure (for a particular track length).We stress that the information flux is used for illustration only and does not enter directly into the calculation of projected upper limits or the discovery reach.
The statistical techniques outlined above are implemented in the software package swordfish4 [97].This package provides a straightforward interface, allowing the user to input signal and background spectra, and efficiently computes the projected exclusion limits and discovery reach.Although the forecasting techniques developed in Ref. [95] are applicable to unbinned data, the implementation in swordfish requires, for practical reasons, that the data be binned.We use 70 log-spaced bins throughout this work, with the range defined by the resolution as σ x T /2 x T 1000 nm.The number of bins was chosen by incrementally increasing the bin width until the projected constraints on the cross section begin to weaken.This way we use the minimum number of bins required to resolve all features in the track length spectra.We do this both to minimize computation time and to ensure that the systematics study in Sec.III is conservative (using a small number of bins enhances the impact of bin-to-bin variations in the backgrounds).The same number of bins is used in both the high resolution and high exposure cases to allow for an easier comparison between the two read-out methods.

E. Track Length Spectra
We now present the distribution of track lengths expected for DM signals as well as the backgrounds already described in this section.These are shown in Fig. 1, with each panel showing a different mineral.We fix the DM-nucleon cross section to σ SI n = 10 −45 cm 2 and show the DM signal spectra for three different DM masses, m χ = 5, 50 and 500 GeV, as solid lines.We see that at short track lengths, x T 10 nm, where signals of lighter WIMPs would appear, the dominant background often comes from solar neutrinos.At longer track lengths, where the signal from heavier DM typically peaks, the dominant backgrounds are radioactive: mono-energetic 1α-thorium recoils and neutron-induced recoils.
In the upper part of each panel, we also plot the information flux, as described in Sec.II D. The most pronounced feature is a sharp drop for track lengths corresponding to 1α-thorium recoils.Little information about a DM signal can be gained by studying tracks of this length, as signal events are here degenerate with 1αthorium tracks.The various peaks in the information flux indicate the track lengths that provide the most constraining power for the DM signal.These maxima appear either when the signal is large, corresponding to a large signal-to-background ratio, or when one of the backgrounds is prominent.In the latter case, track lengths corresponding to peaks in the information flux allow us to constrain the normalization of a particular background component.This in turn leads to improved constraints on the signal.In between these peaked regions, the information flux is typically suppressed.Note that the detailed shape of the information flux depends significantly on the specific assumptions that are made about the correlation of background systematics.
We discuss the track length spectra, information fluxes and their impact on paleo-detector sensitivity in more detail in the next section.

III. PROJECTED SENSITIVITY
Here we present the results of the sensitivity analysis.First we discuss a simple benchmark case in which we consider systematic errors only on the normalization of the individual background components.We also use this scenario for the mass reconstruction projections in Sec.IV.In addition, we also consider bin-to-bin systematics in order to assess the impact of shape uncertainties of the background spectra on the projected sensitivity.

A. Background Normalization Systematics
The backgrounds described in Sec.II all have an associated uncertainty which must be accounted for within the analysis framework.For our background normalization systematics scenario we assume that the shapes of the signal and backgrounds are fixed with only a systematic uncertainty on the normalization of each background component.The systematic uncertainties we assign to the normalization of different backgrounds are detailed in Sec.II C. We ignore covariances between the signal and background and between individual backgrounds components.With careful calibration, we may be able to understand the shape of the background to a high degree of precision.In practice, it should be straightforward to produce target samples with high levels of radioactivityinduced tracks in the laboratory, though such an approach is more challenging for the neutrino-induced backgrounds.
Uncertainties in the normalization of the backgrounds can be mitigated when a good 'control region' is available, where the signal is sub-dominant and the overall background rate can be well-constrained.This is typically the case for the broad distribution of neutron-induced recoils; even for heavy WIMPs, the signal drops off well below 1000 nm, providing a good control region at large track lengths.Instead, neutrino-induced backgrounds may mimic the DM signal for certain DM masses.If this is the case, no control region is available and limits are severely weakened.
We show in Fig. 2 the projected sensitivity for the four minerals we consider in this work.The top two panels show the projected 90% confidence exclusion limits while the bottom two panels show the projected 5 σ discovery reach [98], which we define as the line above which the paleo-detector setups we consider would have at least 50 % chance of achieving a 5 σ discovery of DM.In gray, we show current bounds from conventional direct detection experiments, coming predominantly from XENON1T [8] and SuperCDMS [11] in this mass and cross section range.
In the left panels, we show the high resolution case, with σ x T = 1 nm and an exposure of = 0.01 kg Myr.The 'bump' in the limits at m χ ∼ 7 GeV is due to the WIMP signal spectrum near this mass being mimicked by the spectra of 8 B and hep neutrinos.This 'Solar neutrino floor' has been studied in detail in, for example, Refs.[64,91,99].The loss in sensitivity is even more pronounced in our case, owing to the larger systematic uncertainty we have assigned to the Solar neutrino flux.In addition, moving away from m χ ∼ 7 GeV, the spectra are no longer degenerate, meaning that control regions are rapidly regained and sensitivity restored.In contrast to conventional direct detection experiments, paleo-detectors have large enough exposures to directly constrain the normalization of the 8 B/hep Solar neutrino fluxes.
Even accounting for these low energy background neu- trinos, in the left panels, the limits at low mass are substantially stronger than limits from any conventional detector.This is due to the ability of modern imaging techniques to measure small track lengths.Tracks of 1 nm in length typically correspond to recoil energies around 100 eV, giving a threshold comparable to the CRESST-2017 Surface Run [12], albeit with a much larger exposure.Thus, for DM lighter than ∼ 10 GeV, paleo-detectors may probe DM-nucleon cross sections much smaller than the projected sensitivity of conventional direct detection experiments.
For larger DM masses, the sensitivity is severely limited by neutron-induced backgrounds.Nchwaningite and sinjarite contain hydrogen, which is an efficient modera-tor of fast neutrons.Thus, the rate of neutron-induced backgrounds is much smaller than in olivine and halite, which do not contain hydrogen.In addition, for halite and sinjarite we assume an intrinsic contamination of C 238 = 0.01 ppb, whereas for olivine and nchwaningite we assume C 238 = 0.1 ppb.For nchwaningite and sinjarite, the minerals containing hydrogen, cross sections as low as 10 −46 cm 2 could be probed for a DM mass of 50 GeV, assuming an exposure of 0.01 kg Myr, as shown in the left panels.At higher masses, the projected sensitivity of paleo-detectors is comparable to current XENON1T bounds.
In the right panels of Fig. 2, we show the projected sensitivity for the lower resolution case, with σ x T = 15 nm

FIG. 2.
Projected paleo-detector sensitivity using a spectral analysis.Top: Projected 90% Upper Limits.Bottom: 5 σ Discovery Reach.In the left column, we assume σx T = 1 nm; = 0.01 kg Myr (high resolution case).In the right column, we assume σx T = 15 nm; = 100 kg Myr (high exposure case).We assume a 100 % systematic uncertainty on the normalization of each individual neutrino component, as well as a 1 % systematic uncertainty on the normalization of both the neutron background and the 1α-thorium background.Gray regions show the parameter space currently excluded by conventional direct detection experiments [8,14].In the upper panels, we also show the projected limits from nchwaningite using a sliding window analysis (dotted purple), as reported recently in Ref. [39].and a larger exposure of = 100 kg Myr.For DM heavier than 100 GeV, it is possible to probe DM-nucleon cross sections a factor of 30 and 100 smaller than current XENON1T bounds using nchwaningite and sinjarite, respectively.For DM lighter than 10 GeV, the sensitivity is marginally weaker than in the high resolution case; at low DM mass, the signal spectrum peaks towards shorter track lengths and is thus challenging to resolve with worse resolution.
We note that the projected limits are comparatively weak between 20 and 100 GeV.This is because the peak due to 1α-thorium recoils (vertical dashed lines in Fig. 1), typically coincides with the peak in the signal spectra in this mass range 5 .We see in the upper panels of Fig. 1 that the broad peak in the information flux for the 50 5 For mχ ∼ 50 GeV, the typical recoil energy for nuclei in the GeV case typically occurs at the same position as the 'dip' caused by the 1α-thorium background.
We now compare our results with those obtained using a simpler 'sliding window' analysis in Refs.[38,39], where a cut-and-count analysis is performed over a window chosen to optimize the signal-to-noise ratio.In the top two panels of Fig. 2, we show the limit obtained in Ref. [39] for nchwaningite.
In the high resolution case (left panels of Fig. 2), we find that the 'sliding window' analysis is only marginally less sensitive at the highest masses.Near m χ ∼ 50 GeV, the full spectral analysis is more sensitive by a factor of a few due to better rejection of the sharply peaked minerals we consider is ∼ 20 keV, much smaller than the 72 keV 1α-thorium recoil.However, the stopping power for lighter nuclei such as Ca, Cl and Na is smaller, leading to similar track lengths.
1α-thorium background.Going to lower masses, the 'bump' corresponding to the 8 B/hep Solar neutrinos becomes more pronounced in a full spectral analysis.For lighter DM, the neutrino-and DM-induced spectra become distinguishable again.The spectral analysis can effectively reduce the error on the normalization of the neutrino-induced background by using information from track lengths where neutrino-induced events dominate.For example, at masses of 1 GeV this allows the projected sensitivity from the spectral analysis to improve nearly an order of magnitude with respect to that found in the 'sliding window' analysis.
In the high exposure case (right panels of Fig. 2), the spectral analysis gains roughly an order of magnitude in sensitivity with respect to the 'sliding window' analysis for DM heavier than ∼ 100 GeV.The lower resolution makes it more difficult to exploit subtle differences in the shape of the signal and background spectra.However, at the longest track lengths considered, the tracks induced by neutrons always dominate over those induced by DM, cf.Fig. 1.Because of the larger exposure, this 'control region' has sufficient statistics to tightly constrain the normalization of the neutron-induced background.Thus, the sensitivity to higher mass WIMPs is improved with respect to the 'sliding window' analysis.

B. Background Shape Systematics
As discussed in the previous section, our sensitivity to DM at high masses is typically limited by neutroninduced backgrounds.Conversely, neutrinos are the dominant background for low DM mass.As we show in Fig. 3, the maximum signal-to-background ratio (as a function of track length) for the high exposure case is typically much smaller than 10 %.This means that the strength of the limits depends in principle on a %-level understanding of the shape of the background distributions.For the high resolution case the maximum signal-to-background is typically closer to 30 %.Therefore, the projected sensitivity should be more robust to shape uncertainties in that case.
In order to explore how the sensitivity of paleodetectors is affected by such background shape uncertainties, we assign a Gaussian systematic error to the normalization of each bin of each background component.Our analysis therefore allows the backgrounds in individual bins to fluctuate independently.Because we no longer assign systematic uncertainties to the overall normalization of the backgrounds, in some situations (e.g. when the bin-to-bin uncertainty is chosen to be smaller than the normalization uncertainty in Sec.III A) the projected limit with shape uncertainties may be stronger than in the normalization systematics case.In such situations, we set the projected limit equal to the normalization systematics case.
In Fig. 4, we plot the limits obtained with various bin-to-bin background systematics for a sinjarite paleo- 3. Maximum signal-to-background ratio at the 90 % confidence limit.Both lines show the maximum signal-to-background ratio (over all track lengths) for sinjarite with σ SI n set to the projected 90 % confidence exclusion limit at each mass.The orange (dashed) line shows the high-resolution case.Here we see the signal-to-background is relatively constant around 10 − 30 %.The blue (solid) line shows the high-exposure case.At low masses the maximum signal-to-background is roughly 10 % whereas at high masses this is reduced to 0.4−0.5 %.This is reflected by the increased sensitivity of the limit to bin-to-bin background shape uncertainties, as shown in Fig. 4.
detector.For comparison, we also show the sensitivity obtained in the background normalization systematics case described in the previous sub-section.These results are mostly illustrative since the actual level of shape uncertainties is hard to anticipate a priori.However, they indicate the level of uncertainty that can be tolerated in practice.
As expected, allowing some variation in the shape of the backgrounds degrades the sensitivity of paleodetectors.In the upper panel of Fig. 4 we show the high resolution case.The limit is unaffected by bin-to-bin systematics until they are increased to 50 %.The more relevant uncertainty is therefore the overall normalization of the background components.
The high exposure case is shown in the lower panel of Fig. 4. We find a much greater dependence on the background shape systematics.For 10 % bin-to-bin systematic uncertainties, the sensitivity is degraded by an order of magnitude compared to the background normalization systematics case for DM heavier than ∼ 40 GeV.Another factor of ∼ 5 is lost when increasing the bin-tobin systematics from 10 % to 50 %.
The high resolution case is more robust to background shape systematics primarily because of its larger signalto-background ratio, as shown in Fig. 3.In addition, spectral features in the high resolution case are more pronounced, allowing for easier distinction between signal and background even when there are significant uncertainties in the background shapes.Projected Limits: LUX-Zeplin SuperCDMS SNOLAB FIG. 4. Projected 90 % confidence limits for sinjarite including background shape systematics.Top and bottom panels show the high resolution and high exposure cases, respectively.Blue line: Background normalization systematics case with systematic normalization uncertainties for each background component.The normalization systematic on neutrinos here is set to 100 % whereas for the radioactive backgrounds we assume 1 % normalization error.Red, Orange, and green lines: Background shape systematics case where we allow the normalization of each background component in each of the 70 log-spaced track-length bins to fluctuate independently.The red, green, and orange lines show results for 50 %, 10 %, and 1 % systematic uncertainty per bin, respectively.Note that where the bin-to-bin systematics produce a limit stronger than that of the normalization systematics case, we set the projected limit assuming normalization systematics.In the top panel the 1% (orange dotted) and 10 % (green dashed) bin-to-bin systematic lines are therefore not distinguishable from the normalization systematics case (blue solid).Shaded regions show projected limits from LUX-Zeplin [37] and SuperCDMS SNOLAB (Ge) [33].
For comparison with our projections, we also show in Fig. 4 the projected 90 % confidence exclusion limits from LUX-Zeplin [37] and SuperCDMS SNOLAB (Ge) [33] (planned for data-taking from 2020 onwards), with respective exposures of 1.5×10 4 kg yr and 44 kg yr.For DM masses below 10 GeV, the high resolution case can improve upon future SuperCDMS SNOLAB constraints by up to one and a half orders of magnitude.For the case of 50 % bin-to-bin shape systematics, sinjarite would still be an order of magnitude more sensitive than SuperCDMS SNOLAB projections at 2 GeV.The high exposure case can achieve the same sensitivity as LUX-Zeplin to higher mass DM only if the background shape uncertainties are kept at the 1 % level.

IV. CONSTRAINING THE DARK MATTER MASS
In this section, we investigate to what extent the properties of a DM candidate, in particular its mass, could be constrained in the hypothetical case of a DM discovery.Thus, we switch from projecting limits, as in Sec.III, to parameter reconstruction.A priori there is no reason for the DM to appear in any particular region of the parameter space.Instead of employing benchmark scenarios as often done in the literature [100][101][102][103], we perform a benchmark-free study using the Euclideanized signal method [97,104].
The Euclideanized signal method maps points in the model parameter space -here, (m χ , σ SI n ) -onto points in a 'Euclideanized signal' space, taking into account systematic uncertainties, covariances and nuisance parameters.Likelihood ratios between points in the model space are mapped to Euclidean distances in the 'Euclideanized signal' space.Comparing the Euclidean distances between large numbers of points is computationally fast (using clustering algorithms), allowing us to efficiently map out the reconstruction prospects over a wide range of the parameter space.Full details of the Euclideanized signal method can be found in Refs.[97,104] and a summary is given in App. A.
For regions where a future signal would provide a closed constraint on the DM mass, we calculate the accuracy to which this is possible by defining the fractional uncertainty as Here, m χ,max and m χ,min are the maximum and minimum edges of the two-dimensional 2 σ confidence contour around a point with a given (m χ , σ SI n ).Note that these contours are typically quite asymmetric, usually extending much further towards masses larger than the true mass than towards masses smaller.Thus, a fractional uncertainty ∆m χ /m χ 1 does not necessarily imply that no information about the DM mass can be obtained.Rather, ∆m χ /m χ 1 typically implies that ∆m χ /m χ ∼ m χ,max /m χ , while usually m χ,min is not much smaller than m χ .
In Fig. 5, we show the ability of a sinjarite paleodetector to constrain the DM mass from a future signal.The color scale shows contours in ∆m χ /m χ ; in the following we refer to the colored regions as those where the mass can be reconstructed.The gray regions indicate points in the parameter space where the 2 σ confidence contours are not closed, i.e. the reconstructed mass would be unbounded from above.Thus, we refer to the gray regions as portions of parameter space where the mass cannot be reconstructed.In Fig. 5, we quantify how well the mass could be reconstructed for DM-nucleon cross sections between the projected 90 % confidence exclusion limits, cf.Sec.III A, and cross sections a factor 100 larger than this (the region bounded by the two black curves).Note that some of this region is already ruled out by XENON1T, cf.Fig. 2.
Direct detection experiments suffer from an almost exact degeneracy between the mass and cross section at large DM masses [105].The degeneracy occurs when m χ m N .This is because, for a given nuclear recoil energy, v min depends only on the reduced mass of the DM-nucleus system.For m χ m N , the reduced mass becomes independent of the DM mass.
For traditional direct detection experiments, the ability to reconstruct the mass of a hypothetical DM particle has been studied extensively, see e.g.[103][104][105][106][107].These studies found that a Xenon experiment would only be able to constrain the DM mass up to ∼ 200 GeV.For paleo-detectors, we instead find that mass reconstruction is possible for DM masses as large as ∼ 1 TeV.
We show results for the high-resolution configuration in the upper panel of Fig. 5.In the case of DM-nucleon cross sections just below the current limits, we find that the largest DM mass for which the mass could be reconstructed is ∼ 250 GeV.In the high exposure configuration (Fig. 5, lower panel) the DM mass could be reconstructed for masses as high as ∼ 1 TeV if the DM-nucleon cross section is just below current limits.For cross sections further below current limits, we see that in the high resolution case the colored region extends to slightly larger masses than in the high exposure case.
For a DM candidate with mass below ∼ 10 GeV, the mass could be reconstructed for all cross sections in reach of a sinjarite paleo-detector in either the high resolution or high exposure configuration.Although the precision of the mass reconstruction depreciates with decreasing DMnucleon cross section, there are large regions of parameter space which are currently unconstrained and where the DM mass could be reconstructed reasonably well by paleo-detectors in case of a discovery.
The potential for paleo-detectors to tightly constraint the mass of a DM candidate stems in part from their large exposure.For example, in the high resolution case, a 6 GeV DM candidate with cross section at the 90 % confidence exclusion limit would give rise to ∼ 10 5 signal events.Such large numbers of events would allow us to accurately map out the track length spectrum.For DM masses of ∼ 1 TeV, paleo-detectors would only measure O( 103 ) events at the exclusion limit, making the reconstruction of the associated track length spectrum more challenging.This should be contrasted with ex- Constraints on the mass of a DM particle from a future signal.Gray shaded regions correspond to parameter points where the DM mass is unconstrained from above at the 2 σ-level.The colored contours indicate the fractional uncertainty on the DM mass obtained by constraining a future signal, as defined in Eq. ( 14).Top and bottom panels show the high resolution and high exposure cases, respectively.The lower black lines in both panels correspond to the projected 90 % confidence limit in Fig. 2. We consider regions between these lower lines and a factor of 100 larger, indicated by the upper black lines.Note that some of these regions are already excluded by current experiments (see shaded regions in Fig. 2).
posures possible in conventional direct detection experiments, where at their exclusion limits only O(1-10) events would be detected.Such a low number of signal events would not provide enough information to resolve minute differences in the recoil spectra required for mass reconstruction.
Further, paleo-detectors could probe track lengths over three orders of magnitude, which corresponds to sensitivity spanning a large range of recoil energies.In particular, a 1000 nm track corresponds roughly to a ∼ 1 MeV nuclear recoil whereas a 1 nm track equates to a ∼ 100 eV recoil 6 .The high energy part of the recoil spectra has a significant dependence on the DM mass [102,103].Unlike the traditional energy window of direct detection experiments, we exploit this information by observing a wide variety of track lengths.
Finally, the target materials we consider here contain a variety of constituent nuclei with different masses.Sinjarite contains nuclei with masses from 1 GeV (H) to 37 GeV (Ca).Since the observed signal is a weighted sum over the contributions from the respective target nuclei, the resulting track length spectrum is richer in features than the simple slope dependence one finds for a single target nucleus.As these features can be exploited efficiently with a spectral analysis, paleo-detectors are particularly well suited for DM mass reconstruction.
We note that we have not considered any uncertainties in the DM velocity distribution in this analysis.For example, the longest tracks (which we leverage to constrain very heavy DM) are produced by the fastest moving DM particles, close to the Galactic escape velocity.The length of these longest tracks is therefore dependent on uncertainties in the escape velocity.More generally, allowing for variations in the DM distribution will widen the constraints on the DM mass.A number of techniques have been developed to incorporate velocity distribution uncertainties in direct detection (see for example Refs.[100,[108][109][110][111][112][113]) but we leave this more detailed analysis to future work.We note however that previous studies have shown that using multiple experiments with different target nuclei greatly reduces the impact of these uncertainties [102,104].We expect the same to be true for paleo-detectors since the minerals investigated in this work contain a variety of target nuclei.

V. CHALLENGES
Throughout this work we have shown that the sensitivity of paleo-detectors may well exceed that of current direct detection experiments.In Sec.III A, we projected the sensitivity assuming systematic errors on the overall normalization of the different background components only.For the neutrino-induced backgrounds, we assumed 100 % systematics, while we assumed 1 % systematics for backgrounds induced by radioactivity.In order to check the robustness of our results, we increased the normalization systematics on the neutron-induced backgrounds to 5 % and found that the sensitivities are unaffected.In the following, we discuss some of the other potential issues moving forward.
In our background normalization systematics case we assumed no covariance between the normalization of the background components.This assumption should hold for many background components, for example we expect no covariance between the spectra induced by solar neutrinos and diffuse supernova neutrinos.For the radiogenic backgrounds there may exist some covariance since they have a common origin.
For the background shape systematics case, we assume no bin-to-bin covariances.This may be an optimistic or pessimistic assumption depending on the covariances one might expect.The most troublesome scenario would be a bin-to-bin covariance that makes signal-like variations in the background more likely.Due to the lack of theoretical guidance, we have chosen not to explore bin-to-bin covariances.Instead we attempt to maximize the error introduced by bin-to-bin systematics by using the minimum number of bins required to resolve all features in the spectra.We leave careful study of covariances to future analyses.
One of the primary assumptions throughout our analysis is that we can reject damage features in the minerals that are not tracks arising from nuclear recoils.Further, we do not consider the background produced by a series of linked α-recoil tracks in the uranium-238 decay chain.
The assumption is that the characteristic track pattern is easily recognizable and therefore rejected with 100% efficiency.In reality there will be > O(10 7 ) tracks within a sample, many of which will exhibit these characteristic track patterns.We therefore require an automated tagging and rejection system.Since the characteristic track pattern is quite distinct from a normal track, it is possible that this system can be very efficient.However, this is yet to be validated for a large data set, a task we leave to future work.
The data produced by scanning significant amounts of material at high precision could present an issue by itself.Naively, scanning 1mm 3 of material at 1nm precision will produce 10 6 terabytes of data.It is a monumental task to analyze such a large data set.Luckily, much of the mineral will contain no track information at all, therefore a suitable compression format can be adopted to make the analysis more tracktable.Analysis of the data will require automated track-recognition, an ideal application of machine learning algorithms.We will also address this in future publications.
Naively, one would expect paleo-detectors to be able to exploit directional information from the orientation of the recoil tracks.However, the target minerals we consider here are O(1 Gyr)-old, which is comparable to the period of the Sun's revolution around the Galactic center.Also, geological processes occur on timescales shorter than O(1 Gyr), further complicating the expected directionality of the DM-induced signal.Reference [55] attempted to quantify the directional dependence of the DM-induced tracks within ancient minerals, showing that there is a preferred direction.Unfortunately for an O(1 Gyr) mineral the effect in ancient mica was calculated to be only O(1%).Because it is unlikely that we will be able to resolve the head/tail orientation of tracks at the nm scale, the induced anisotropy would need to be much larger than O(1%) in order to be statistically observable [114].
Finally, the translation of the range of the nucleus x T to the reconstructed track length after read-out is a source of uncertainty.Quantifying such an uncertainty requires detailed studies for different combinations of minerals and read-out methods [39].However, in the case of a claimed detection, we would be able to confirm a signal using minerals with different constituent nuclei and ages, allowing one to mitigate some of these systematic issues.

VI. CONCLUSIONS
In this work, we have explored the prospects for probing Weakly Interacting Massive Particle (WIMP) Dark Matter (DM) with paleo-detectors.In particular, we have extended previous studies by performing a full spectral analysis, including information about the expected distributions of track lengths left in the minerals by a DM signal as well as by neutrino-induced and radiogenic backgrounds.We further explored how systematic uncertainties on the normalization and shape of backgrounds impact projected limits.Finally, we have studied how well the DM mass could be measured in case of a future discovery.
We .Sinjarite is the most sensitive out of the minerals examined here due to the assumed low levels of radioactive contamination and efficient neutron moderation by hydrogen (Ref.[39] came to a similar conclusion for the mineral epsomite [Mg(SO 4 )•7(H 2 O)]).
For moderate track length resolutions, σ x T = 15 nm, we find that the full spectral analysis extends the projected paleo-detector sensitivity to DM-nucleon cross sections roughly an order of magnitude smaller than the sliding window analysis of Refs.[38,39].This improvement is driven by the fact that a full spectral analysis automatically entails the use of optimal 'control regions', where the signal is sub-dominant, helping to pin down the normalization and shape of the backgrounds.
We find that by analyzing an O(10 cm 3 ) sample of sinjarite using small angle X-ray scattering, it could be possible to probe DM-nucleon cross sections roughly a factor of 100 smaller than current direct detection experiments for DM heavier than ∼ 100 GeV.Including systematic uncertainties in background shapes at the 10 %-level, projected limits remain a factor of 7 − 8 more stringent than current XENON1T bounds [8].The sensitivity depreciates further if systematics larger than 10 % are assumed for the shapes of backgrounds.
Analyzing smaller samples of O(1 mm 3 ) at nmresolution (e.g. using helium ion beam spectroscopy), we find that paleo-detectors may be able to probe DMnucleon scattering cross sections many orders of mag-nitude below current limits, for 500 MeV m χ 10 GeV.Probing O(nm) track lengths corresponds to an O(100 eV) energy threshold, exploring significant regions of the recoil spectra from low-mass WIMPs.With highresolution read-out methods, the limits would be robust to systematic uncertainties in the background shapes as large as ∼ 50 %.
In addition, we have investigated the prospects for paleo-detectors to constrain the DM parameters in the case of a future signal.As an example, we calculate the regions in which the mass and cross section become degenerate for a sinjarite paleo-detector.We find that below m χ 15 GeV, it would be possible to reconstruct the mass of the DM particle with a relative error of less than 10 % if the cross section is large enough for a 5 σ discovery.For m χ 15 GeV the signal becomes increasingly insensitive to changes in m χ , making the mass harder to constrain.In spite of this, paleo-detectors should be able to obtain both a lower and an upper limit on the DM mass for m χ 1 TeV if the cross section is just below current limits.In contrast, conventional direct detection experiments could provide only a lower bound on the DM mass if the true mass is larger than ∼ 200 GeV [104].
Paleo-detectors could also be used to investigate a number of interesting questions beyond searches for DM.For example, in the absence of a DM signal, the analysis outlined here would straightforwardly allow us to measure neutrino-induced events.It could therefore be possible to use mineral samples of different ages as a unique probe of the neutrino history of our galaxy.This will potentially allow us to study both historical neutrino processes in the Sun and the signal from supernovae.
Paleo-detectors represent an excellent opportunity to probe large areas of the WIMP DM parameter space in the near future.The next steps involve assessing the practical challenges of reaching the required exposures to achieve these sensitivities, as well as more detailed modeling of backgrounds.We leave both of these tasks to future work.However, we note that WIMP-nucleon cross sections much smaller than projected in this work may be probed by paleo-detectors if either novel ideas to control the backgrounds emerge (akin to progress made in conventional direct detection experiments in recent decades) or if target materials with significantly lower levels of radioactivity are available.With uranium concentrations of 10 −15 , radioactive backgrounds would no longer dominate at high DM masses.In such a case, paleo detectors could perhaps probe WIMP-nucleon cross sections all the way down to the diffuse supernova and atmospheric neutrino floor.

FIG. 1 .
FIG. 1.Track-length spectra and information flux for three different WIMP masses.Each panel shows the spectrum of track lengths xT expected for a different mineral target.WIMP signals are shown as solid lines, while background distributions are shown as dotted and dashed lines.We fix the signal normalization to σ SI n = 10 −45 cm 2 .The information flux, shown above each set of spectra, is a generalized signal-to-noise ratio discussed in more detail in Sec.II D. The information flux has been calculated for the high-resolution case, with a track-length resolution of σx T = 1 nm.

2 ] 2 ] 2 ] 2 ]
Projected 90% Upper limits; σ xT = 1 nm; = 0.01 kg Myr Projected 90% Upper limits; σ xT = 15 nm; = WIMP-nucleon cross section σ SI n [cm 5σ discovery reach; σ xT = 1 nm; = 0.01 kg Myr WIMP-nucleon cross section σ SI n [cm 5σ discovery reach; σ xT = 15 nm; = 100 kg Myr Nchwaningite Olivine Sinjarite Halite FIG.3.Maximum signal-to-background ratio at the 90 % confidence limit.Both lines show the maximum signal-to-background ratio (over all track lengths) for sinjarite with σ SI n set to the projected 90 % confidence exclusion limit at each mass.The orange (dashed) line shows the high-resolution case.Here we see the signal-to-background is relatively constant around 10 − 30 %.The blue (solid) line shows the high-exposure case.At low masses the maximum signal-to-background is roughly 10 % whereas at high masses this is reduced to 0.4−0.5 %.This is reflected by the increased sensitivity of the limit to bin-to-bin background shape uncertainties, as shown in Fig.4.

2 ] 2 ]
Projected 90% Upper limits; σ xT = 1 nm; = 0.01 kg Myr Sinjarite -CaCl 2 •2(H 2 O) WIMP-nucleon cross section σ SI n [cm Projected 90% Upper limits; σ xT = 15 nm; = 100 kg Myr FIG.5.Constraints on the mass of a DM particle from a future signal.Gray shaded regions correspond to parameter points where the DM mass is unconstrained from above at the 2 σ-level.The colored contours indicate the fractional uncertainty on the DM mass obtained by constraining a future signal, as defined in Eq. (14).Top and bottom panels show the high resolution and high exposure cases, respectively.The lower black lines in both panels correspond to the projected 90 % confidence limit in Fig.2.We consider regions between these lower lines and a factor of 100 larger, indicated by the upper black lines.Note that some of these regions are already excluded by current experiments (see shaded regions in Fig.2).