New Projections for Dark Matter Searches with Paleo-Detectors

Paleo-detectors are a proposed experimental technique to search for dark matter (DM). In lieu of the conventional approach of operating a tonne-scale real-time detector to search for DM-induced nuclear recoils, paleo-detectors take advantage of small samples of naturally occurring rocks on Earth that have been deep underground ($\gtrsim 5$ km), accumulating nuclear damage tracks from recoiling nuclei for $\mathcal{O}(1)$ Gyr. Modern microscopy techniques promise the capability to read out nuclear damage tracks with nanometer resolution in macroscopic samples. Thanks to their $\mathcal{O}(1)$ Gyr integration times, paleo-detectors could constitute nuclear recoil detectors with keV recoil energy thresholds and 100 kilotonne-yr exposures. This combination would allow paleo-detectors to probe DM-nucleon cross sections orders of magnitude below existing upper limits from conventional direct detection experiments. In this article, we use improved background modeling and a new spectral analysis technique to update the sensitivity forecast for paleo-detectors. We demonstrate the robustness of the sensitivity forecast to the (lack of) ancillary measurements of the age of the samples and the parameters controlling the backgrounds, systematic mismodeling of the spectral shape of the backgrounds, and the radiopurity of the mineral samples. Specifically, we demonstrate that even if the uranium concentration in paleo-detector samples is $10^{-8}$ (per weight), many orders of magnitude larger than what we expect in the most radiopure samples obtained from ultra basic rock or marine evaporite deposits, paleo-detectors could still probe DM-nucleon cross sections below current limits. For DM masses $\lesssim 10$ GeV/$c^2$, the sensitivity of paleo-detectors could still reach down all the way to the conventional neutrino floor in a Xe-based direct detection experiment.


Introduction
More than a century after evidence for the existence of Dark Matter (DM) in our Universe started to emerge [1], we still do not know what DM is made out of. Weakly Interacting Massive Particles (WIMPs) remain one of the best-motivated DM candidates to date [2][3][4]. Among the most sensitive probes of the WIMP DM paradigm are direct detection experiments [5][6][7][8][9], where one searches for interactions of DM with atomic nuclei. Over the past few decades, direct detection experiments have made impressive progress: liquid noble gas detectors have reached O(tonne-yr) exposures with nuclear recoil energy thresholds as small as a few keV [10][11][12][13], while cryogenic bolometric detectors have pushed the nuclear recoil energy threshold below 100 eV, albeit as yet only with detectors accumulating O(kg days) of exposure [14][15][16][17][18]. A large experimental program is underway to extend the sensitivity of direct detection experiments both to lower WIMP masses and smaller WIMP-nucleon cross sections [19][20][21][22][23], as well as to develop directional detectors [24][25][26][27][28][29][30][31][32][33][34]. However, these detectors are becoming ever more costly and more challenging to operate.
Paleo-detectors [35][36][37] are a proposed alternative technique to search for interactions of DM with nuclei: instead of operating a nuclear recoil detector in real time, as in a conventional direct detection experiment, one would search for the persistent traces of WIMPnucleon interactions recorded in natural minerals over timescales as long as O(1) Gyr. Many minerals are excellent solid state nuclear track detectors [38][39][40][41] and retain these traces (i.e. damage tracks) caused by a recoiling atomic nucleus for timescales which can exceed the age of the Earth by many orders of magnitude (i.e. 4.5 Gyr). Natural minerals commonly found on Earth are up to O(1) Gyr old and can be dated to a few-percent accuracy [42][43][44]. Modern microscopy technologies, including electron microscopy [40,45], Helium Ion Beam Microscopy (HIBM) [46,47], Small Angle X-ray scattering (SAXs) [48][49][50], and optical microscopy [41,51,52] potentially allow for the readout of such nuclear damage tracks in macroscopic samples of target material. For example, HIBM may allow one to read out O(10) mg of material with O(1) nm precision when combined with pulsedlaser and fast-ion-beam ablation techniques [53][54][55][56]. On the other hand, SAXs promises spatial resolutions of O(15) nm in O(100) g of target material [49,50] (see the discussion in Refs. [35,36]). These advanced microscopy techniques enable paleo-detectors to reach much larger exposures with better track length resolutions than previous ideas to probe rare events using natural minerals; see Refs. [57][58][59][60][61][62][63][64][65][66][67][68][69][70][71] and, in particular, the work by Snowden-Ifft et al. searching for damage tracks from WIMP-nucleus interactions in ancient phlogopite mica samples [72][73][74][75]. While not the focus of this article, paleo-detectors have also been proposed as a probe of astrophysical neutrino fluxes from Galactic supernovae [76], cosmic ray interactions in Earth's atmosphere [77], and from our Sun [78]. There has been recent related work on using crystal-defects created in certain minerals by WIMP-nucleus interactions to search for light DM [79][80][81] or to monitor nuclear reactors [82], and to use very long damage tracks in ancient quartz samples to search for ultra-heavy DM [83].
Using natural minerals found on Earth to record DM-nucleus interactions over O(100) Myr− O(1) Gyr timescales provides paleo-detectors with two crucial advantages over conventional direct detection experiments. First, the enormous integration times promise exposures many orders of magnitude larger than what is feasible in laboratory experiments. Reading out 100 g of material that has been recording nuclear damage tracks for 1 Gyr yields an exposure of ε = 100 g Gyr, or, in units more appropriate for a real-time experiment, ε = 100 kilotonne yr. The O(15) nm spatial resolution achievable with SAXs readout corresponds to a nuclear recoil energy threshold of O(1) keV. Thus, paleo-detectors would combine the recoil energy threshold achievable in direct detection experiments utilizing liquid noble gases as target materials with the exposures achievable in large (planned) neutrino detectors such as Super-/Hyper-Kamiokande or DUNE. Second, while conventional direct detection experiments measure events in real time, paleo-detectors integrate over timescales up to O(1) Gyr, comparable to the age of the Solar System, t ∼ 4.5 Gyr, and the orbital period of the Solar System around the Milky Way, T ∼ 250 Myr. Using a series of samples with different ages, paleo-detectors offer a unique pathway to gaining information about a signal's time-evolution over Gyr timescales [76][77][78]84].
Of course, paleo-detectors must pay a price for these advantages: compared to a conventional detector operated in a controlled laboratory environment, background mitigation might be much more challenging. The background sources for DM searches in paleodetectors are qualitatively the same as in conventional direct detection experiments since both approaches are, at their heart, nuclear recoil searches. Extensive background studies for paleo-detectors have been performed in previous work [35,36,76] and we will discuss the different sources further below. As we will see, the most important background sources are cosmic rays, solar neutrinos, and radioactivity. Regarding the first, the fundamental strategy for paleo-detectors is to use natural minerals extracted from deep within the Earth where the rock-overburden has shielded them from cosmogenic backgrounds. Conventional direct detection experiments are operated in deep underground laboratories for the same reason. However, paleo-detectors require only (at most) a few kg of target material compared to actively operating a tonne-scale experiment; samples of such modest sizes can be extracted from much greater depths than those of existing underground laboratories (for example from existing boreholes drilled for geological R&D and oil exploration). For the backgrounds from solar neutrinos and radioactivity it is important to note that, due to their large exposures, paleo-detectors would contain a sizable number of background events. A rare event search, such as for DM-induced nuclear recoils, would thus not proceed as in conventional direct detection experiments where one typically attempts to find a background-free signal region. Instead, a paleo-detector analysis is more similar to a collider physics search: one would model the distribution of background events (for example in track-length space) and search for deviations of the observed data from the backgroundonly prediction using a signal template. One important consequence of the large number of background (and possible signal) events is that one can characterize many aspects of the backgrounds from the data themselves. Just as in collider searches, paleo-detectors can use control regions, where the backgrounds dominate, to better characterize the background contribution in the signal region of the search.
The remainder of this article is organized as follows: in Section 2 we discuss the calculation of DM signals in paleo-detectors followed by the most important background sources in Section 3. In Section 4 we show projections for the reach of paleo-detectors in the conventional WIMP mass-cross section plane. Compared to previous estimates of the sensitivity of paleo-detectors to DM signals [35][36][37], in this work we make use of both improved background modeling and new statistical techniques. In this sense, the results presented here supersede the sensitivity projections shown in Refs. [35][36][37]. In Section 4 we also explore the robustness of the sensitivity projections to changes in the assumptions we make in our analysis; for example, we remove external constraints on the age of the minerals and the parameters controlling the background normalizations, and explore the effect of systematic mismodeling of the spectral shape of the backgrounds. Furthermore, we show how the projected sensitivity changes with the radiopurity of the samples, and demonstrate that, even if the uranium concentration in paleo-detector samples is 10 −8 (per weight, denoted as "g/g" throughout), many orders of magnitude larger than what we expect in the most radiopure samples obtained from ultra basic rock or marine evaporite deposits, paleo-detectors could still probe DM-nucleon cross sections below current limits. Specifically, for WIMP masses m χ 10 GeV/c 2 , the sensitivity of paleo-detectors could still reach down all the way to the conventional neutrino floor in a Xe-based direct detection experiment. We summarize and conclude in Section 5. We make the code used in this work available: paleoSpec [85] for the computation of the signal and background spectra, and paleoSens [86] for the sensitivity forecasts.

Dark Matter Signals in Paleo-Detectors
The signature of WIMP DM in a paleo-detector is the spectrum of damage tracks left by nuclei in the mineral sample after receiving a "kick" from WIMP-nucleus scattering. In this section, we discuss how this signal is computed. The differential rate of nuclear recoils per unit target mass for a WIMP with mass m χ scattering off target nuclei N with recoil energy E R is given by [87,88] dR 1) where A N is the atomic mass number of N , µ p = m p m χ /(m p +m χ )is the reduced mass of the WIMP-proton system with the proton mass m p , and we have assumed isospin-conserving spin-independent WIMP-nucleon interactions with the (zero momentum transfer) WIMPproton cross section σ SI p . Expressions analogous to Eq. (2.1) for isospin-violating SI interaction and spin-dependent WIMP-nucleon interactions can be found in Refs. [87][88][89][90][91], and for more general WIMP-nucleus interactions in Refs. [92,93] . Note that while we focus on spin-independent interactions in this work, paleo-detectors could also probe spin-dependent WIMP-nucleus interactions, see Ref. [36].
In Eq. (2.1), the nuclear form factor F N (E R ) accounts for the finite size of the nucleus; we will use the Helm parameterization [94][95][96] in our numerical calculations. 1 The distribution of WIMPs is described by the local WIMP (mass) density, ρ χ , and the mean inverse speed, of the WIMP relative to the detector. In Eq. (2.2), v min = m N E R /2µ 2 N accounts for the scattering kinematics; v min is the minimal speed, with respect to a nucleus of mass m N , that a WIMP must have to give rise to a recoil with energy E R (where µ N is the reduced mass of the WIMP-nucleus system). In our numerical calculations we set ρ χ = 0.3 GeV/cm 3 while for f (v) we use a Maxwell-Boltzmann velocity distribution in the Galactic rest frame with velocity dispersion σ v = 166 km/s [101], truncated at the escape speed v esc = 550 km/s [102], and boosted to the Solar System frame by v = 248 km/s [103] as in the Standard Halo Model (SHM) [7,95,104].
In order to obtain the differential track length spectrum (per unit target mass), dR/dx, we use the range of a nucleus N with energy E R , where (dE/dx) N is the stopping power of N in the target material. The true length of a damage track in a given target may differ from the range if, for example, the nucleus creates a lasting damage track only along some portion of the length it travels through the material, or if the nucleus' trajectory deviates significantly from a straight line. Our previous studies suggest that such effects should be small (see the discussion in Ref. [36]), making Eq. (2.3) a good proxy for the track length. Furthermore, the effects of thermal annealing could potentially be significant over geological timescales; fortunately, any associated modifications to the track lengths would be similar for both the signal and background recoils. While detailed experimental studies for each target material will be required to account for these possible corrections, for the purposes of this work we will use Eq. (2.3) for the track length and SRIM [106,107] to calculate stopping powers. The left panel of Figure 1 shows Right: Differential rate of tracks per track length and unit target mass in gypsum from a 5 GeV/c 2 and a 500 GeV/c 2 WIMP compared to background spectra induced by neutrinos (ν), radiogenic neutrons (n), and 238 U → 234 Th + α recoils ( 234 Th ), see the discussion in Section 3. For the normalization of the WIMP DM spectra we set the spin-independent WIMP-nucleon cross sections to σ SI p = 10 −43 cm 2 for the 5 GeV/c 2 case and σ SI p = 10 −46 cm 2 for the 500 GeV/c 2 case, compatible with current upper bounds [12,105]. the range of the different constituent nuclei 2 in gypsum [Ca(SO 4 )·2(H 2 O)] as an example; we see that a 1 keV recoil gives rise to a few-nm long track, while a 100 keV recoil produces a track with a length of a few hundred nm.
In a composite target material composed of different isotopes with mass fraction ξ N , the differential rate at which tracks are produced (per unit target mass) is then given by with dR/dE R given in Eq. (2.1). In the right panel of Figure 1 we show dR/dx for a 5 GeV/c 2 (with cross section σ SI p = 10 −43 cm 2 ) and a 500 GeV/c 2 (with σ SI p = 10 −46 cm 2 ) WIMP in gypsum as an example together with various background spectra which we will discuss in Section 3. Note that for the 5 GeV/c 2 case, the largest contribution to the signal is at track lengths of a few nm, while heavier WIMPs, such as the 500 GeV/c 2 case shown in Figure 1, give rise to longer tracks with typical lengths of order 100 nm.
The nuclear recoil tracks in a paleo-detector sample can be read out with a variety of readout techniques as briefly discussed in the introduction (see Ref. [36] for a more detailed discussion). Here, we will consider two scenarios: • High-resolution readout scenario: We assume that 10 mg of material can be read out with a track length resolution of σ x = 1 nm. This scenario may be achievable with Helium Ion Beam Microscopy (HIBM) combined with pulsed-laser and fast-ion-beam ablation techniques [46,47,[53][54][55][56]. As we will see, this scenario is advantageous for low-mass (m χ 10 GeV/c 2 ) WIMP searches.
• High-exposure readout scenario: We assume that 100 g of material can be read out with σ x = 15 nm. This scenario could be realized via Small Angle X-ray scattering (SAXs) tomography at a synchrotron facility [48][49][50]; it is better suited for heavier (m χ 10 GeV/c 2 ) WIMP searches.
These are the same readout scenarios as what we have used in previous paleo-detector studies [35-37, 76, 77]. We note that these scenarios represent challenging applications of the respective techniques; demonstrating their feasibility and especially their scalability to the proposed sample masses would be an important step towards the realization of a paleo-detector experiment. There are ongoing efforts to demonstrate the readout of few-keV recoils using a variety of microscopy techniques at JAMSTEC [108] and at Rensselaer Polytechnic Institute [109].
In order to assess the sensitivity of paleo-detectors, we will consider the binned track length spectrum as the experimental observable. Accounting for the finite resolution of the readout method, the number of tracks in the i-th bin with reconstructed length x ∈ x min i , x max i in a sample of mass M which has been recording tracks for a time t age is where dR/dx is given by Eq. (2.4) and we have assumed that dR/dx is independent of time. 3 The window function, W , accounts for the finite resolution effects of the readout. In our numerical calculations, we will assume that the probability of measuring a track length x from a track with true length x is Gaussian-distributed with variance σ 2 x . The corresponding window function is Note that in order to avoid artificial sensitivity to tracks which are much shorter than the spatial resolution of the readout, we remove all tracks with true track length x < σ x /2 from the spectra (dR/dx) before applying Eq. (2.5). For all numerical results in this work, we use 100 log-spaced bins between σ x /2 and 1000 nm. Figure 2 shows the same 5 GeV/c 2 and 500 GeV/c 2 WIMP DM spectra as the right panel of Figure 1 together with the various background sources after taking into account finite resolution effects. The left panel is for the high-resolution scenario, while the right panel is for the high-exposure scenario. We will discuss this figure further at the end of Section 3 after introducing the relevant background sources. Let us here only note that the high-resolution scenario is better suited for relatively light WIMPs with masses m χ 10 GeV/c 2 , while the high-exposure scenario is aimed at heavier WIMPs, m χ 10 GeV/c 2 . This is for two reasons: First, lighter WIMPs give rise to softer nuclear recoils, resulting in (on average) shorter signal tracks, making excellent track length resolution crucial to  Figure 2. Track length spectra for a 5 GeV/c 2 and a 500 GeV/c 2 WIMP compared to the background spectra induced by neutrinos (ν), radiogenic neutrons (n), and 238 U → 234 Th + α recoils ( 234 Th ) after taking finite resolution effects into account. The left panels are for our high resolution scenario (M = 10 mg of sample material read out with track length resolution of σ x = 1 nm) while the right panels are for the high exposure scenario (M = 100 g, σ x = 15 nm). We assume that the samples have been recording tracks for t age = 1 Gyr. In all panels, we are using 100 logarithmically spaced bins between σ x /2 ≤ x ≤ 10 3 nm; note that these bins and the scales of the axes differ between the left and the right panels. The upper panels show the number of tracks per bin from the respective sources. For the backgrounds, the shaded bands around the respective spectra show the associated error (systematic and statistical added in quadrature). The bottom panels show the ratio of the number of signal (S i ) to the (summed) number of background (B i ) events per bin, and the sand-colored shade shows the relative uncertainty of the total number of background events per bin. The signal-to-noise ratio per bin can be read off from the lower panels by dividing S i /B i for the respective signal with the relative uncertainty of the background. For the normalization of the WIMP DM spectra we set the spin-independent WIMP-nucleon cross sections to σ SI p = 10 −43 cm 2 for the 5 GeV/c 2 case and σ SI p = 10 −46 cm 2 for the 500 GeV/c 2 case, as in Figure 1.
resolving the details of the signal and background spectra in the region where the ratio of signal to background events is largest for such light WIMPs. Second, current bounds from direct detection experiments rule out much smaller cross sections for heavier WIMPs with m χ 10 GeV/c 2 , hence, searches for WIMPs with m χ 10 GeV/c 2 and cross sections smaller than current upper limits require very large exposures.

Backgrounds
The background sources for DM searches with paleo-detectors are the same as for conventional direct detection experiments. However, the relative importance of the respective sources is different because, compared to conventional direct detection experiments, paleo-detectors use much smaller target masses ( 1 kg) and integrate over much longer timescales (t age ∼ 0.1 − 1 Gyr) without information about the timing of individual signal or background events. The experimental observables in paleo-detectors are nuclear damage tracks with, for all practical purposes, perfect rejection of electronic recoils. Furthermore, natural defects in minerals are either single-site or span across the entire (mono)-crystalline volume and do not resemble nuclear damage tracks. Hence, the only relevant backgrounds are nuclear recoils. In the remainder of this section we will give a brief overview of the most important background sources; we refer the reader to Refs. [36,76] for more detailed discussions.
• Cosmogenics: In order to mitigate backgrounds from cosmic rays, minerals to be used as paleo-detectors must have been shielded by a large overburden for as long as they have been recording nuclear damage tracks (just as conventional direct detection experiments are operated in deep underground laboratories). Paleo-detectors require only (at most) a few kg of target material rather than actively operating a tonnescale experiment; such modest amounts of materials can be sourced from even greater depths than those of existing underground laboratories. One promising source of samples are (existing) boreholes drilled for geological R&D or oil exploration. For example, for an overburden of 5 km rock, the cosmogenic-muon-induced neutron flux is O(10 2 ) cm −2 Gyr −1 [110], leading to negligible backgrounds for the purposes of a paleo-detector. We note that at depths 6 km rock, backgrounds from atmospheric neutrinos producing neutrons in the vicinity of the target become comparable to the backgrounds from cosmogenic muons [111]. We also note that paleo-detector samples can be stored near the surface after extraction from deep in the Earth before readout. For example, the cosmogenic-muon-induced neutron flux in a ∼50 m deep storage facility is 0.2 cm −2 yr −1 .
• Astrophysical Neutrinos: Neutrinos scattering off nuclei in a paleo-detector sample give rise to nuclear damage tracks. We include neutrinos from the Sun, supernovae, and the interactions of cosmic rays with the Earth's atmosphere in our background modeling. We take the solar and atmospheric neutrino fluxes from Ref. [112]. Because of the long integration times, paleo-detectors are sensitive not only to neutrinos from supernovae in far-away galaxies throughout our Universe (the Diffuse Supernova Neutrino Background, DSNB), but also to neutrinos from local supernovae -the supernova rate in the Milky Way is estimated to be 2-3 per 100 years [113][114][115][116][117][118]. We compute the DSNB spectrum and the contribution from Galactic supernovae as in Ref. [76], see also Refs. [118][119][120][121][122]. The neutrino-induced background is dominated by solar neutrinos at track lengths below ∼100 nm, with supernova neutrinos giving the largest contributions for tracks of a few 100 nm, before atmospheric neutrinos dominate the neutrino-induced backgrounds at even longer track lengths. Note that while neutrino-induced nuclear recoils are a background for DM searches, they can also be an interesting signal for paleo-detectors, see Refs. [76][77][78].
• Radiogenics: Any natural mineral used as a paleo-detector will contain trace amounts of radioactive materials. In order to mitigate radiogenic backgrounds, it is crucial to use minerals with as low a concentration of radioactive elements as possible. The most important radioactive isotope for paleo-detectors is 238 U . Typical minerals formed in the Earth's crust have 238 U concentrations of order C 238 ∼ 10 −6 g/g, leading to prohibitively large backgrounds for DM searches. However, minerals which constitute so-called Ultra Basic Rocks (UBRs), formed from the material of the Earth's mantle, and Marine Evaporites (MEs), formed from evaporated sea water, have typical 238 U concentrations orders of magnitude lower [123][124][125][126][127][128][129], (see also Ref. [36] and especially the discussion in the appendix of Ref. [76]). As in previous works on paleo-detectors [35-37, 76, 77] we will assume benchmark 238 U concentrations of C 238 = 10 −10 g/g for UBRs and C 238 = 10 −11 g/g for MEs.
Since paleo-detectors record only nuclear recoils, the most relevant radiogenic background arises from α-decays and spontaneous fission. Note that it is not the αparticles themselves or the fission fragments which lead to backgrounds for DM searches: α-particles (i.e. He nuclei) do not give rise to lasting damage tracks in most materials and, even if they would, the range of the α's is of order 10 µm in natural minerals, orders of magnitude larger than the tracks from DM-induced nuclear recoils. Similarly, the fission fragments from spontaneous fission of heavy elements have typical energies of tens of MeV, producing much longer tracks than DM-induced recoils. However, in an α-decay, the decaying nucleus recoils with typical energies of some tens of keV against the α-particle. Due to the long integration time of paleo-detectors, almost all of the 238 U nuclei in the mineral which undergo the first ( 238 U → 234 Th + α) decay in the 238 U decay chain will have continued on to decay to the stable 206 Pb. There are eight α-decays in the 238 U decay chain, leading to a characteristic pattern of connected nuclear recoil tracks. We expect such patterns of connected tracks to be easily distinguishable from the isolated recoil tracks DM-nucleus interactions would cause and, hence, do not consider them a relevant background for DM searches. However, the second α-decay in the 238 U decay chain ( 234 U → 230 Th + α) has a relatively long half-life, T 1/2 ( 234 U) = 0.25 Myr. This leads to a population of isolated tracks from 238 U atoms which have undergone the initial ( 238 U → 234 Th + α) decay, but not the second α-decay ( 234 U → 230 Th + α) in the uranium series [130,131]. The 234 Th nucleus receives 72 keV of recoil energy in the ( 238 U → 234 Th + α) decay, leading to a population of monochromatic isolated tracks with approximate number density n( 234 Th) 10 6 g −1 × C 238 /10 −11 g/g [36].
The second important source of radiogenic backgrounds are neutrons: spontaneous fission of heavy isotopes as well as (α, n)-reactions between the α-particles produced in α-decays and nuclei comprising the target material give rise to neutrons with MeV energies. These neutrons in turn scatter off the nuclei in the paleo-detector sample, giving rise to a broad spectrum of nuclear recoil tracks. In order to model this background, we use SOURCES-4A [132] to compute the neutron spectra from spontaneous fission and (α, n)-reactions of α's from all α-decaying nuclei in the 238 U decay chain; note that the contributions from spontaneous fission and (α, n)-reactions to the neutron flux are typically of the same order of magnitude, the precise balance depends on the chemical composition of the target. From the neutron spectra, we calculate the induced nuclear recoil spectra using TENDL-2017 [133][134][135][136] neutron-nucleus cross sections tabulated in the JANIS4.0 [137] database. 4 The neutron-induced backgrounds are strongly suppressed in target materials containing hydrogen: since neutrons and H nuclei have approximately the same mass, neutrons lose a large fraction of their energy in a single interaction with hydrogen, leading to efficient moderation of the radiogenic neutrons even if hydrogen comprises only a small fraction of the atoms in the target material [35,36].
In Figures 1 and 2 we show the track-length spectra produced by these various background sources together with two benchmark DM signals in gypsum [Ca(SO 4 )·2(H 2 O)]. At short track lengths, the background budget is dominated by (solar) neutrino induced nuclear recoils. At intermediate ranges between a few tens and a hundred nm (the precise values depend on the particular mineral), the isolated 234 Th tracks from 238 U → 234 Th+α decays contribute, before the background budget becomes dominated by radiogenic-neutroninduced backgrounds at larger track lengths. Comparing to the benchmark WIMP spectra, we see that solar neutrinos are the most relevant background source for WIMPs with masses in the few GeV/c 2 range, while radiogenic neutrons are the limiting background in searches for heavier WIMPs with masses larger than ∼10 GeV/c 2 .
In Figure 2 we show the (binned) number of tracks one would observe in the highresolution (left panel) and high-exposure (right panel) scenarios after accounting for the effects of finite readout resolution. Because of the large exposure of paleo-detectors, a sizable number of background events will be present. A DM search with paleo-detectors would thus proceed quite differently from a traditional direct detection experiment, in which one typically attempts to construct a signal region with very few (or, ideally, zero) background events, and then searches for an O(1) number of DM-induced events. Given the large number of background events, in a paleo-detector analysis one would model the distribution of background events (for example, in track-length space) and search for deviations of the observed data from the background-only prediction. Thus, the uncertainty in the background prediction is much more important than the absolute number of background events.
For illustration, we include error bands around the respective background components in Figure 2. These bands include both the (statistical) Poisson noise and a systematic error in the background prediction. For the latter, we make well-motivated assumptions in order to illustrate how the sensitivity of paleo-detectors to signals from DM of various masses would depend on different background components in a cut-and-count analysis, see previous studies of paleo-detectors [35-37, 76, 77]. However, we note that in the spectral analysis discussed in Section 4, the sensitivity of paleo-detectors to a WIMP DM signal is virtually independent of any prior knowledge of the background normalizations. Let us briefly describe our choices for the systematic errors in Figure 2. Neutrino-induced backgrounds are controlled by sources which are not associated with the environment of the paleodetector target sample: the fluxes of solar, supernova, or atmospheric neutrinos. These fluxes may vary by an O(1) factor over the ∼ Gyr integration time of paleo-detectors [76][77][78], hence, we assign a (relative) systematic uncertainty of ±100 % to the neutrino-induced backgrounds. 5 The radiogenic backgrounds, on the other hand, are controlled by the 238 U concentrations which do not change over time. Furthermore, the shape of the track length spectra induced by radiogenic backgrounds can be calibrated by studying the track length distributions in mineral samples with large concentrations of 238 U . Finally, the 238 U concentration in a given sample can be measured using, for example, mass spectrometry or gamma ray spectroscopy [141,142]. Hence, we expect the uncertainty on the radiogenic background predictions to be much smaller and assign them a (relative) systematic error of ±1 %.
From the bottom panels in Figure 2 we can directly read off the signal-to-noise ratio (per bin), i.e. the ratio of the number of signal events to the uncertainty of the background prediction: in these panels, the lines show the ratio of the number of signal to the (total) number of background events (S/B), while the sand-colored shaded region shows the relative uncertainty of the background prediction. For a given signal, the ratio of S/B to the (relative) background uncertainty directly corresponds to the significance one would obtain in a simple cut-and-count analysis. As discussed in Refs. [35,36], any signal for which S/B is larger than the shaded region indicating the background uncertainty in the bottom panels of Figure 2 in at least one bin could be probed with a simple cut-and-count analysis.
Given the large number of events expected in a paleo-detector, one can go beyond the simple cut-and-count analysis and perform a spectral analysis in track-length space: Armed with predictions for the spectral shape of the various background components and of a possible signal, one can fit all of these spectral templates to the data and ask if the data is better described by either the backgrounds-only or the backgrounds+signal hypothesis. Heuristically, such a search uses control regions, where some background dominates, to measure that background component, and uses this information to reduce the error of the background prediction in the signal region, where the signal-to-background ratio is largest. For example, short track lengths can serve as a control region for the solarneutrino-induced backgrounds (as well as the 234 Th tracks in the high-exposure scenario), while track lengths longer than ∼500 nm can provide an excellent control region for the neutron-induced backgrounds. We have first used such a spectral analysis in Ref. [37] and will discuss an alternative approach for the spectral analysis in the next section.

Sensitivity Forecasts
In Section 2 we described how to compute the DM signal in paleo-detectors. In particular, we showed how to obtain the track length spectrum one would observe in a paleo-detector from the differential recoil rate per recoil energy and how to account for finite track-length resolution effects. In Section 3 we discussed the most relevant background sources and briefly sketched how one would perform a DM search with a paleo-detector. Due to the potentially large exposure, a sizable number of background (and, possibly, signal) events will be present in any paleo-detector sample. In order to search for DM, one would perform a spectral analysis of the observed recoil spectrum in track-length space and ask if that spectrum is better fit by a background-only model or by a background+signal hypothesis. In Refs. [35,36] we have used a much simpler cut-and-count analysis to project the sensitivity of paleo-detectors to DM, while in Ref. [37] we have estimated the sensitivity using a spectral analysis based on the Fisher-matrix formalism using the swordfish package [143][144][145]. In the remainder of this section, we describe a profile likelihood ratio based approach to calculate the projected sensitivity of paleo-detectors to a DM signal with a spectral analysis. As we will see, one advantage of using this statistical approach is that we can easily incorporate the effect of background mismodeling into our sensitivity forecast. We show results for the sensitivity forecasts in Figures 3 -5.
The predicted track length spectrum in a paleo-detector is To easily compare to existing upper limits and the projected sensitivities of conventional direct detection experiments, we calculate the projected upper (exclusion) limit on σ SI p at 90 % confidence level as a function of m χ for a paleo-detector. The set of nuisance parameters in Eq. (4.1) is then Since paleo-detectors are ultimately a counting experiment, we expect the data to be Poisson-distributed. The appropriate log-likelihood to observe the data D given the parameters θ; σ SI p , m χ is thus 6 The nuisance parameters will be constrained by other measurements. For example, the neutrino fluxes will be constrained by existing neutrino experiments and theoretical modeling. The mass, age, and uranium concentration of the sample will be constrained by direct measurements of the sample mineral. In a frequentist analysis, these constraints can be taken into account by jointly analyzing these ancillary measurements with the track length spectrum data. To mimick such a joint analysis, we include a set of external (Gaussian) constraints on the nuisance parameters, Here, θ j is the central value for the j-th nuisance parameter supplied by the ancillary measurement(s), and c j is the relative uncertainty of the measurement(s). The likelihood is then ln L D θ; σ SI p , m χ = ln L Poisson D θ; σ SI p , m χ + ln L ext. const. (θ) . To calculate the projected sensitivity of paleo-detectors we consider the maximum log-likelihood ratio [146][147][148] In the numerator,θ is the set of nuisance parameters that maximizes L for given values of σ SI p (and m χ ), while in the denominator, L is maximized over both θ and σ SI p (for a given value of m χ ), with best-fit valuesθ andσ SI p . In order to calculate the projected exclusion limit, we use the Asimov data set [146] under the assumption of no signal,  Table 1. List of minerals considered in this work with their chemical composition and our fiducial assumption for the 238 U concentration (C 238 ) in radiopure samples of these minerals.
where θ is a set of fiducial choices for the nuisance parameters. The 90 % confidence level exclusion limit for a given value of m χ is obtained by finding the smallest value of σ SI p for which q(σ SI p ; m χ ) ≥ q crit = 2.71. 7 We describe the the fiducial values for the nuisance parameters listed in Eq. (4.2) and the relative uncertainty of the external constraints on those parameters, c j , defined by Eq. (4.4). For the high-resolution (high-exposure) scenario, we choose a fiducial value of M = 10 mg (M = 100 g) for the sample mass and, for both scenarios, t age = 1 Gyr for the time the sample has been recording damage tracks. We assume that these parameters can be constrained by ancillary measurements to c M = 0.01 % and c tage = 5 %; while measuring the mass of a target sample to high precision is trivial, its age can be inferred with few-percent precision using standard geological dating techniques [42][43][44]. For the normalization of the neutrino fluxes, Φ i ν , and the 238 U concentrations in our target minerals, C 238 , we take fiducial values described in Section 3 (the latter are also summarized in Table 1). As a benchmark, we assume that the associated nuisance parameters can be constrained to c Φ i ν = 100 % and c C 238 = 1 %; these choices are equivalent to the respective systematic uncertainties of the backgrounds assumed in previous studies [35-37, 76, 77]. However, as we will discuss below, the projected sensitivity of paleo-detectors to a DM signal is largely independent of external constraints on the nuisance parameters.
In Figure 3 we show the projected 90 % confidence level upper limit under these fiducial assumptions in both the high-resolution (sample mass M = 10 mg, track-length resolution σ x = 1 nm; left panel) and the high-exposure (M = 100 g, σ x = 15 nm; right panel) readout scenario. We show results for a selection of minerals, see Table 1: three examples of Ultra-Basic Rocks (UBRs), the very common olivine and two less common UBRs which contain hydrogen, phlogopite mica and nchwaningite. We also show the projected sensitivity for three examples of Marine Evaporites (MEs): halite, gypsum, and sinjarite where the latter two contain hydrogen. As discussed above, the high-resolution scenario has better sensitivity for WIMP masses m χ 10 GeV/c 2 where the dominant background   Table 1. The gray-shaded region of parameter space is disfavored by current upper limits from direct detection experiments [12,14,17,105,150], while the sand-colored region indicates the neutrino floor for a Xe-based experiment [151]. Colors and linestyles are the same in both panels. stems from solar neutrinos; here, excellent track length resolution is crucial to resolving the differences between the shapes of the signal and background track length spectra. For m χ 10 GeV/c 2 , the high-exposure scenario has better sensitivity since for such WIMP masses radiogenic backgrounds are dominant and maximizing the exposure is crucial in order to achieve the best sensitivity. Comparing the sensitivity for the different target materials we can note that materials containing hydrogen generally have better sensitivity than those without hydrogen. This is due to the effective suppression of the radiogenic neutron backgrounds by hydrogen, see the discussion in Section 3. Furthermore, we note that MEs typically have better sensitivity than UBRs; this is because we assume a lower 238 U concentration (C 238 = 10 −11 g/g) for MEs than for UBRs (C 238 = 10 −10 g/g). Out of the minerals considered in Figure 3, throughout the shown WIMP mass range and for both the high-resolution and the high-exposure readout scenarios, sinjarite promises the best sensitivity. Primarily, this is due to sinjarite being an ME containing hydrogen, and hence having relatively small 238 U -induced backgrounds. However, the sensitivity is also enhanced by its particular chemical composition: the other example of an ME containing hydrogen shown in Figure 3, gypsum, is very oxygen-rich, and hence, both the DM-induced and the background spectra are dominated by the contribution of oxygen. Sinjarite on the other hand has one calcium, two oxygen, and two chlorine atoms; each of these elements can therefore contribute comparably to the recoil spectra. This adds additional structure to the track length spectra which is advantageous in a spectral analysis.
In Figure 3 we also show existing limits from conventional direct detection experiments [12,14,17,105,150] as well as the neutrino floor for a Xe-based experiment [151] with the shaded regions. At small WIMP masses, m χ 10 GeV/c 2 , paleo-detectors could probe WIMPs with WIMP-nucleon cross sections many orders of magnitude smaller than those currently probed by conventional experiments; for example, at m χ = 1 GeV/c 2 the best existing upper limit is σ SI p 10 −38 cm 2 [14,17]. Moreover, paleo-detectors offer sensitivity deep into the conventional (solar) neutrino floor in the m χ 10 GeV/c 2 mass range. Due to the large exposures, a paleo-detector would observe O(100 − 1000) recoil tracks from solar neutrinos in the high-resolution readout scenario, see Figure 2. In a spectral analysis, this large number of neutrino-induced events combined with the excellent tracklength (and, in turn, recoil energy) resolution would allow one to characterize the solar neutrino background in the "signal region" and hence probe WIMP-nucleon cross sections deep into the conventional neutrino floor. For heavier WIMPs (m χ 10 GeV/c 2 ), we see that paleo-detectors could probe WIMP-nucleon cross sections roughly one and a half orders of magnitude below current upper limits, but somewhat larger than the conventional neutrino floor in a Xe-based direct detection experiment. In this WIMP mass region, the sensitivity of paleo-detectors is limited by the radiogenic neutron backgrounds, as we will discuss further below.
In the remainder of this section, we will explore the robustness of the sensitivity projections shown in Figure 3 to the assumptions we have made. Let us first note that the projected exclusion limits are virtually independent of the ancillary measurements. For fixed normalization of the backgrounds in the Asimov data, removing the Gaussian constraints [defined by Eq. (4.4)] on the parameters controlling the normalization of the neutrinoinduced backgrounds (Φ i ν ) and the radiogenic backgrounds (C 238 ) as well as removing the external constraint on the age of the samples (t age ) leads to a depreciation of the projected upper limits by less than a factor two across the WIMP mass range shown in Figure 3. We find the most pronounced depreciation of the upper limits in the m χ 1 GeV/c 2 WIMP range when the constraints on the Φ i ν are removed. Hence, the data from the track-length spectrum alone suffices to search for DM with paleo-detectors, even in the absence of any ancillary information about the normalization of the backgrounds.
In order to explore how robust the sensitivity forecasts are with respect to mismodeling of the shape of each background component, we replace the Poisson log-likelihood, Eq. (4.3), with a χ 2 distribution, where σ 2 i is the variance in the i-th bin. If we include only the Poisson error of the data, σ 2 i → D i , the Poisson log-likelihood and the χ 2 likelihood become asymptotically equivalent in the large-N i limit. Given the number of signal events one would observe in a paleo-detector, see Figure 2, we thus expect to find the same sensitivity for both forms of the likelihood. The χ 2 distribution allows us to include additional uncertainties, which we model as  where ∆ sys Bkg parameterizes the (relative) systematic modeling error. We note that it is straightforward to extend this to σ 2 , allowing for the assignment of a systematic error to the j-th contribution to the data in the i-th bin. In this work, we will only employ a simple overall systematic error ∆ to explore the robustness of the sensitivity projections.
In Figure 4 we show the projected upper limit for a gypsum paleo-detector under different assumptions for the background-shape mismodeling. First, we note that for ∆ sys Bkg = 0, i.e. accounting only for the Poisson error in the variance, σ 2 i = D i , both the Poisson loglikelihood and the χ 2 distribution yield the same projected limit. In the high-resolution scenario (left panel of Figure 4) we find that if we include an error of ∆ sys Bkg = 0.1, corresponding to 10 % bin-to-bin mismodeling of the shape of the backgrounds, the sensitivity does not depreciate significantly. For ∆ sys Bkg = 1, the projected upper limit is up to an order of magnitude worse than in the ∆ sys Bkg = 0 case for m χ 30 GeV/c 2 , with smaller changes at larger WIMP masses. In the high-exposure scenario (right panel of Figure 4), the sensitivity is not quite as robust against background-shape mismodeling. For small systematic errors, ∆ sys Bkg = 0.01, we find that the sensitivity becomes only slightly worse than for ∆ sys Bkg = 0, while for ∆ sys Bkg = 0.1 the projected upper limit is a factor of ∼ 5 larger than for ∆ sys Bkg = 0. The high-exposure-readout scenario is most relevant for the mass range m χ 10 GeV/c 2 , where the dominant background stems from radiogenic neutrons. The corresponding track length spectrum is straightforward to calibrate in paleo-detectors by measuring the tracks in samples of the same mineral but with larger 238 U concentration. Thus, it seems unlikely that the systematic (bin-to-bin) shape errors would be as large as ∆ sys Bkg = 0.1 in an actual paleo-detector experiment. However, we note from Figure 4 that even for errors as large as ∆ sys Bkg = 0.1, a gypsum paleo-detector could still probe WIMPnucleon cross sections almost an order of magnitude below current experimental limits for m χ 100 GeV/c 2 .
In Figure 5 we show the effect of assuming different 238 U concentrations on the projected upper limits in a gypsum paleo-detector. Trivially, we find that paleo-detectors could probe smaller WIMP cross sections if the samples are more radiopure. For most of the range of 238 U concentrations, 10 −12 g/g ≤ C 238 ≤ 10 −7 g/g shown in Figure 5 and most of the shown WIMP mass range, the projected upper limit scales like σ SI p ∝ √ C 238 . This scaling arises from the uncertainty of the radiogenic backgrounds in the "signal region" being given by the associated Poisson error, with the normalization of the radiogenic backgrounds inferred from the relevant "control regions". The ∝ √ C 238 scaling of the sensitivity is broken if other backgrounds provide significant contributions in the "signal region." For example, the radiogenic background contribution can be subdominant to that of solar neutrinos for m χ 10 GeV/c 2 or, for very small 238 U concentrations, atmospheric neutrinos for m χ 100 GeV/c 2 . Note that in the high-resolution scenario, paleo-detectors could probe WIMP-nucleon cross sections down to the conventional (solar) neutrino floor even for 238 U concentrations as large as C 238 = 10 −8 g/g in the m χ 10 GeV/c 2 mass range. For heavier WIMPs, we find that a gypsum paleo-detector in the high-exposure readout scenario could still probe WIMP-nucleus cross sections below current experimen-tal limits even if the 238 U concentration is orders of magnitude larger than our fiducial assumption of C 238 = 10 −11 g/g.
Before moving on, let us note that the statistical framework used here to project the sensitivity of paleo-detectors can easily be adapted to answer different statistical questions or to incorporate additional uncertainties. For example, by using the Asimov data as defined above (backgrounds only) but also including a DM signal contribution, one can calculate the projected discovery reach rather than the exclusion limit. Using the χ 2 likelihood, it is straightforward to include arbitrarily complicated uncertainties; another approach to including (modeling) uncertainties would be to use random realizations of noisy mock data instead of the Asimov data set. It is also straightforward to extend the framework to simultaneously analyze a series of paleo-detectors. Samples of different ages can be used to explore the sensitivity of paleo-detectors to time-dependent signals. For instance, the Solar System could have traversed some denser region of the DM halo during its past few revolutions around our Galaxy [84]. Time-dependent neutrino signals may also arise from the evolution of our Sun [78], the Galactic supernova rate [76], or the flux of cosmic rays hitting Earth's atmosphere [77]. The code used in this work is available here [86]. Some of the options mentioned in this paragraph are already implemented in the code, and we invite the reader to add the capabilities in order to answer their own analysis questions.

Summary and Conclusions
This work surveys the prospects of searching for DM with paleo-detectors, natural minerals which can record and retain nuclear damage tracks for geological timescales. The exposure times of paleo-detectors can be of the order of a billion years, and modern microscopy techniques may allow for the readout of these damage tracks with nanometer spatial resolution in macroscopic samples. Thus, paleo-detectors could constitute nuclear recoil detectors with ∼ keV recoil energy thresholds but with exposures rivaling those of large (planned) neutrino detectors such as Super/Hyper-Kamiokande or DUNE. Building on our previous work [35-37, 76, 77], we have discussed the relevant backgrounds for DM searches with paleo-detectors and presented updated sensitivity forecasts using a full spectral analysis. For the first time, the background induced by neutrinos from Galactic supernovae is included in our projections and, in this sense, the results presented here supersede the forecasts in Refs. [35][36][37]. We have also put particular emphasis on exploring the robustness of our sensitivity forecasts: we have demonstrated that DM-searches with paleo-detectors are possible even in the absence of any external information about the age of the samples, or the neutrino fluxes and the 238 U concentration, which control the normalization of the dominant backgrounds in paleo-detectors; the track length spectra alone contain sufficient information to infer these parameters. Furthermore, we have shown that the projected sensitivity is robust against uncertainties in the modeling of the spectral shape of the backgrounds. Finally, we have demonstrated that even if the uranium concentration in paleo-detector samples is C 238 = 10 −8 g/g, many orders of magnitude larger than what we expect in the most radiopure samples obtained from ultra basic rock or marine evaporite deposits, paleo-detectors could still probe WIMP-nucleon cross sections below current limits. In particular, for WIMP masses m χ 10 GeV/c 2 , the sensitivity of paleo-detectors could still reach down all the way to the conventional neutrino floor in a Xe-based direct detection experiment.
We hope that this article serves as further motivation for experimental efforts towards paleo-detectors. Ongoing efforts to demonstrate the feasibility of reading out tracks from keV-scale nuclear recoils [108,109] are exciting first steps in this direction. Conventional direct detection experiments have been making impressive progress in their sensitivity to DM during the last few decades. However, despite detectors becoming ever larger and more sophisticated, they have failed to deliver (conclusive) evidence of WIMP-nucleon interactions to date. Future experiments are becoming increasingly expensive and challenging to operate. Thus, paleo-detectors are a timely proposal to probe the remaining WIMP parameter space. Moreover, paleo-detectors have applications to physics beyond DM searches, such as measuring solar, supernova, or atmospheric neutrinos, and offer a unique tool to explore the time-dependence of any signal over gigayear timescales.