New Projections for Dark Matter Searches with Paleo-Detectors

Abstract

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 (≳5 km), accumulating nuclear damage tracks from recoiling nuclei for $\mathcal{O}\left(1\right)$Gyr. Modern microscopy techniques promise the capability to read out nuclear damage tracks with nanometer resolution in macroscopic samples. Thanks to their $\mathcal{O}\left(1\right)$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 ≲ 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.

1. 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 $\mathcal{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 of yet only with detectors accumulating $\mathcal{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 WIMP–nucleon interactions recorded in natural minerals over timescales as long as $\mathcal{O}\left(1\right)$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 that 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 $\mathcal{O}\left(1\right)$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 $\mathcal{O}\left(10\right)$mg of material with $\mathcal{O}\left(1\right)$nm precision when combined with pulsed-laser and fast-ion-beam ablation techniques [53,54,55,56]. On the other hand, SAXs promises spatial resolutions of $\mathcal{O}\left(15\right)$nm in $\mathcal{O}\left(100\right)$g of target material [49,50] (see the discussion in References [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 References [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 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 $\mathcal{O}\left(100\right)\mathrm{Myr}-\mathcal{O}\left(1\right)$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 $\epsilon =100$g Gyr, or, in units more appropriate for a real-time experiment, $\epsilon =100$kilotonne yr. The $\mathcal{O}\left(15\right)$nm spatial resolution achievable with SAXs readout corresponds to a nuclear recoil energy threshold of $\mathcal{O}\left(1\right)$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 $\mathcal{O}\left(1\right)$Gyr, comparable to the age of the Solar System, $t_{\odot }$ ∼ $4.5$Gyr, and the orbital period of the Solar System around the Milky Way, $T_{\odot }$ ∼ $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 paleo-detectors 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 background-only 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 References [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_{\chi }\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. 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.

2. 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_{\chi }$ scattering off target nuclei N with recoil energy $E_R$ is given by [87,88]

$${\left(\frac{dR}{d{E}_R}\right)}_N=\frac{A_N^2{\sigma }_p^{\mathrm{SI}}{F}_N^2\left(E_R\right)}{2{\mu }_p^2{m}_{\chi }}{\rho }_{\chi }{\eta }_{\chi }\left(v_{\min }\right),$$

(1)

where $A_N$ is the atomic mass number of N, ${\mu }_p=m_p{m}_{\chi }/(m_p+m_{\chi })$ 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) WIMP-proton cross section ${\sigma }_p^{\mathrm{SI}}$. Expressions analogous to Equation (1) for isospin-violating SI interaction and spin-dependent WIMP–nucleon interactions can be found in References [87,88,89,90,91], and for more general WIMP–nucleus interactions in References [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 Reference [36].

In Equation (1), the nuclear form factor $F_N\left(E_R\right)$ accounts for the finite size of the nucleus; we will use the Helm parameterization [94,95,96] in our numerical calculations. (Note that more refined calculations of the form factors are available, although only for a few isotopes; see e.g., References [97,98,99,100]). The distribution of WIMPs is described by the local WIMP (mass) density, ${\rho }_{\chi }$, and the mean inverse speed,

$${\eta }_{\chi }\left(v_{\min }\right)={\int }_{v\ge v_{\min }}{d}^3\mathit{v}\frac{f\left(\mathit{v}\right)}{v},$$

(2)

where $f\left(\mathit{v}\right)$ is the velocity distribution of WIMPs with $\mathit{v}$ ($v=\left|\mathit{v}\right|$) the velocity (speed) of the WIMP relative to the detector. In Equation (2), $v_{\min }=\sqrt{m_N{E}_R/2{\mu }_N^2}$ 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 ${\mu }_N$ is the reduced mass of the WIMP–nucleus system). In our numerical calculations, we set ${\rho }_{\chi }=0.3$GeV/cm${}^3$, while for $f\left(\mathit{v}\right)$, we use a Maxwell–Boltzmann velocity distribution in the Galactic rest frame with velocity dispersion ${\sigma }_v=166$km/s [101], truncated at the escape speed $v_{\mathrm{esc}}=550$km/s [102], and boosted to the Solar System frame by $v_{\odot }=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$,

$$x_N\left(E_R\right)={\int }_0^{E_R}dE{\left|\frac{dE}{dx}\right|}_N^{-1},$$

(3)

where ${\left(dE/dx\right)}_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 Reference [36]), making Equation (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 Equation (3) for the track length and SRIM [105,106] to calculate stopping powers. The left panel of Figure 1 shows the range of the different constituent nuclei in gypsum [Ca(SO${}_4)·2$(H${}_2$O)] as an example. (Note that we do not show the range of H in Figure 1 since we do not expect ions with charge $Z\le 2$ (in particular, H and He) to produce lasting damage tracks in most materials; see the discussion in [36]). We see that a 1 keV recoil gives rise to a track a few nm long, 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 ${\xi }_N$, the differential rate at which tracks are produced (per unit target mass) is then given by

$$\frac{dR}{dx}=\sum _N{\xi }_N{\left(\frac{dR}{d{E}_R}\right)}_N{\left(\frac{dE}{dx}\right)}_N,$$

(4)

with $dR/d{E}_R$ given in Equation (1). In the right panel of Figure 1 we show $dR/dx$ for a 5 GeV/c${}^2$ (with cross section ${\sigma }_p^{\mathrm{SI}}={10}^{-43}{\mathrm{cm}}^2$) and a 500 GeV/c${}^2$ (with ${\sigma }_p^{\mathrm{SI}}={10}^{-46}{\mathrm{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 Reference [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 ${\sigma }_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_{\chi }\lesssim 10$GeV$/c^2$) WIMP searches.
• High-exposure readout scenario: We assume that 100 g of material can be read out with ${\sigma }_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_{\chi }\gtrsim 10$GeV$/c^2$) WIMP searches.

These are the same readout scenarios as what we have used in previous paleo-detector studies [35,36,37,76,77]. We note that these scenarios represent challenging applications of the respective techniques; demonstrating that 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\in \left[x_i^{\min },x_i^{\max }\right]$ in a sample of mass M, which has been recording tracks for a time $t_{\mathrm{age}}$, is

$$N_i=M\times t_{\mathrm{age}}\int d{x}^{\prime }W(x^{\prime };x_i^{\min },x_i^{\max })\frac{dR}{dx}\left(x^{\prime }\right),$$

(5)

where $dR/dx$ is given by Equation (4) and we have assumed that $dR/dx$ is independent of time. (For time-dependent signals, replace $t_{\mathrm{age}}\to {\int }_0^{t_{\mathrm{age}}}dt$ in Equation (5)). 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^{\prime }$ is Gaussian-distributed with variance ${\sigma }_x^2$. The corresponding window function is

$$W(x^{\prime };x_i^{\min },x_i^{\max })=\frac{1}{2}\left[\mathrm{erf}\left(\frac{x^{\prime }-x_i^{\min }}{\sqrt{2}{\sigma }_x}\right)-\mathrm{erf}\left(\frac{x^{\prime }-x_i^{\max }}{\sqrt{2}{\sigma }_x}\right)\right].$$

(6)

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^{\prime }<{\sigma }_x/2$ from the spectra ($dR/dx$) before applying Equation (5). For all numerical results in this work, we use 100 log-spaced bins between ${\sigma }_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_{\chi }\lesssim 10$GeV/c${}^2$, while the high-exposure scenario is aimed at heavier WIMPs, $m_{\chi }\gtrsim 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 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_{\chi }\gtrsim 10$GeV/c${}^2$; hence, searches for WIMPs with $m_{\chi }\gtrsim 10$GeV/c${}^2$ and cross sections smaller than current upper limits require very large exposures.

3. 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_{\mathrm{age}}\sim$  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 References [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 tonne-scale experiment; such modest amounts of materials can be sourced from even greater depths than those of existing underground laboratories. One promising source of samples is (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 $\mathcal{O}\left({10}^2\right){\mathrm{cm}}^{-2}{\mathrm{Gyr}}^{-1}$ [110], leading to negligible backgrounds for the purposes of a paleo-detector. We note that at depths $\gtrsim 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{\mathrm{cm}}^{-2}{\mathrm{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 Reference [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 Reference [76]; see also References [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 References [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}\mathrm{U}$. Typical minerals formed in the Earth’s crust have ${}^{238}\mathrm{U}$ concentrations of order $C^{238}$ ∼ ${10}^{-6}$g/g, leading to prohibitively large backgrounds for DM searches. However, minerals that 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}\mathrm{U}$ concentrations orders of magnitude lower [123,124,125,126,127,128,129], (see also Reference [36] and especially the discussion in the appendix of Reference [76]). As in previous works on paleo-detectors [35,36,37,76,77], we will assume benchmark ${}^{238}\mathrm{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 $\alpha$-decays and spontaneous fission. Note that it is not the $\alpha$-particles themselves or the fission fragments that lead to backgrounds for DM searches: $\alpha$-particles (i.e., He nuclei) do not give rise to lasting damage tracks in most materials, and even if they do, the range of the $\alpha$’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 $\alpha$-decay, the decaying nucleus recoils with typical energies of some tens of keV against the $\alpha$-particle. Due to the long integration time of paleo-detectors, almost all of the ${}^{238}\mathrm{U}$ nuclei in the minerals that undergo the first (${}^{238}\mathrm{U}\to {}^{234}\mathrm{Th}+\alpha$) decay in the ${}^{238}\mathrm{U}$ decay chain will have continued on to decay to the stable ${}^{206}$Pb. There are eight $\alpha$-decays in the ${}^{238}\mathrm{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 $\alpha$-decay in the ${}^{238}\mathrm{U}$ decay chain (${}^{234}\mathrm{U}\to {}^{230}\mathrm{Th}+\alpha$) has a relatively long half-life, $T_{1/2}{(}^{234}\mathrm{U})=0.25$Myr. This leads to a population of isolated tracks from ${}^{238}\mathrm{U}$ atoms that have undergone the initial (${}^{238}\mathrm{U}\to {}^{234}\mathrm{Th}+\alpha$) decay, but not the second $\alpha$-decay (${}^{234}\mathrm{U}\to {}^{230}\mathrm{Th}+\alpha$) in the uranium series [130,131]. The ${}^{234}\mathrm{Th}$ nucleus receives $72$keV of recoil energy in the (${}^{238}\mathrm{U}\to {}^{234}\mathrm{Th}+\alpha$) decay, leading to a population of monochromatic isolated tracks with approximate number density $n\left({}^{234}\mathrm{Th}\right)\simeq {10}^6{\mathrm{g}}^{-1}\times \left(C^{238}/{10}^{-11}\mathrm{g}/\mathrm{g}\right)$ [36].

The second important source of radiogenic backgrounds is neutrons: spontaneous fission of heavy isotopes as well as ($\alpha ,n$)-reactions between the $\alpha$-particles produced in $\alpha$-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 $(\alpha ,n)$-reactions of $\alpha$’s from all $\alpha$-decaying nuclei in the ${}^{238}\mathrm{U}$ decay chain; note that the contributions from spontaneous fission and ($\alpha ,n$)-reactions to the neutron flux are typically of the same order of magnitude, and 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. (We take only elastic neutron-nucleus scattering into account; this yields a conservative estimate of the background because neutrons typically lose a larger fraction of their energies through inelastic processes than in elastic scattering. Note also that our Monte Carlo simulation of the nuclear recoils induced by radiogenic neutrons has been validated with a calculation of the nuclear recoils induced by the same neutron spectra with FLUKA [138,139,140] in Reference [77] for the particular case of a halite target.) 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 Figure 1 and Figure 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}\mathrm{Th}$ tracks from ${}^{238}\mathrm{U}\to {}^{234}\mathrm{Th}+\alpha$ decays contribute, before the background budget becomes dominated by radiogenic-neutron-induced 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 high-resolution (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 $\mathcal{O}\left(1\right)$ 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,36,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 that are not associated with the environment of the paleo-detector target sample: the fluxes of solar, supernova, or atmospheric neutrinos. These fluxes may vary by an $\mathcal{O}\left(1\right)$ factor over the ∼Gyr integration time of paleo-detectors [76,77,78]; hence, we assign a (relative) systematic uncertainty of $\pm 100\%$ to the neutrino-induced backgrounds. (In Figure 2, we plot all neutrino-induced backgrounds together to avoid excessive clutter. In our sensitivity projections, we treat the solar, DSNB, Galactic supernova, and atmospheric-neutrino-induced backgrounds as four independent background contributions.) The radiogenic backgrounds, on the other hand, are controlled by the ${}^{238}\mathrm{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}\mathrm{U}$. Finally, the ${}^{238}\mathrm{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 we assign them a (relative) systematic error of $\pm 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 References [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 solar-neutrino-induced backgrounds (as well as the ${}^{234}\mathrm{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 Reference [37] and will discuss an alternative approach for the spectral analysis in the next section.

4. 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 References [35,36] we used a much simpler cut-and-count analysis to project the sensitivity of paleo-detectors to DM, while in Reference [37], we 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 Figure 3, Figure 4 and Figure 5.

The predicted track length spectrum in a paleo-detector is

$$\begin{array}{cc}\hfill \mathit{N}(\mathit{\theta };{\sigma }_p^{\mathrm{SI}},m_{\chi })& ={\mathit{N}}_{\nu }^{\mathrm{sol}}\left({\mathsf{\Phi }}_{\nu }^{\mathrm{sol}}\right)+{\mathit{N}}_{\nu }^{\mathrm{DSNB}}\left({\mathsf{\Phi }}_{\nu }^{\mathrm{DSNB}}\right)+{\mathit{N}}_{\nu }^{\mathrm{GSNB}}\left({\mathsf{\Phi }}_{\nu }^{\mathrm{GSNB}}\right)+{\mathit{N}}_{\nu }^{\mathrm{atm}}\left({\mathsf{\Phi }}_{\nu }^{\mathrm{atm}}\right)\hfill \\ \hfill & +{\mathit{N}}_{\mathrm{rad}}^{{}^{234}\mathrm{Th}}\left(C^{238}\right)+{\mathit{N}}_{\mathrm{rad}}^n\left(C^{238}\right)\hfill \\ \hfill & +{\mathit{N}}^{\mathrm{DM}}({\sigma }_p^{\mathrm{SI}};m_{\chi }),\hfill \end{array}$$

(7)

where the $\mathit{N}$s are the track length spectra described in Equation (5). Each $\mathit{N}$ is computed from the appropriate recoil energy distribution, given by Equation (1), for the DM signal and the corresponding expressions for the different background components discussed in Section 3. In Equation (7), the terms in the first line describe the contributions induced by solar neutrinos (${\mathit{N}}_{\nu }^{\mathrm{sol}}$), the DSNB (${\mathit{N}}_{\nu }^{\mathrm{DSNB}}$), the Galactic Supernova Neutrino Background (${\mathit{N}}_{\nu }^{\mathrm{GSNB}}$), and atmospheric neutrinos (${\mathit{N}}_{\nu }^{\mathrm{atm}}$). The terms in the second line describe the radiogenic backgrounds, separated into: the isolated ${}^{234}\mathrm{Th}$ decays from ${}^{238}\mathrm{U}$ nuclei, which have undergone only the inital $\alpha$-decay during the time the sample has been recording nuclear damage tracks (${\mathit{N}}_{\mathrm{rad}}^{{}^{234}\mathrm{Th}}$), and the background induced by radiogenic neutrons (${\mathit{N}}_{\mathrm{rad}}^n$). Finally, the term in the third line is the DM signal. For each $\mathit{N}$, we indicate in parantheses the parameter controlling its normalization: the flux (${\mathsf{\Phi }}_{\nu }^i$) for each neutrino contribution, the ${}^{238}\mathrm{U}$ concentration in the sample ($C^{238}$) for the radiogenic backgrounds, and the WIMP-nucleon cross section (${\sigma }_p^{\mathrm{SI}}$) for the DM signal. Note that all contributions are also proportional to the mass of the sample, M, and all contributions except for ${\mathit{N}}_{\mathrm{rad}}^{{}^{234}\mathrm{Th}}$ are proportional to the length of time the sample has been recording damage tracks, $t_{\mathrm{age}}$; we have suppressed these dependencies in Equation (7) for compactness. Finally, the DM signal also depends on the mass of the WIMP, $m_{\chi }$.

To easily compare to existing upper limits and the projected sensitivities of conventional direct detection experiments, we calculate the projected upper (exclusion) limit on ${\sigma }_p^{\mathrm{SI}}$ at 90% confidence level as a function of $m_{\chi }$ for a paleo-detector. The set of nuisance parameters in Equation (7) is then

$$\mathit{\theta }=\left\{M,t_{\mathrm{age}},{\mathsf{\Phi }}_{\nu }^{\mathrm{sol}},{\mathsf{\Phi }}_{\nu }^{\mathrm{DSNB}},{\mathsf{\Phi }}_{\nu }^{\mathrm{GSNB}},{\mathsf{\Phi }}_{\nu }^{\mathrm{atm}},C^{238}\right\}.$$

(8)

Since paleo-detectors are ultimately a counting experiment, we expect the data to be Poisson-distributed. The appropriate log-likelihood to observe the data $\mathit{D}$ given the parameters $\left\{\mathit{\theta };{\sigma }_p^{\mathrm{SI}},m_{\chi }\right\}$ is thus (here and in the following, we drop constant factors in the expression of the likelihood, which cancel in the likelihood ratio we are ultimately interested in):

$$ln{\mathcal{L}}_{\mathrm{Poisson}}\left(\mathit{D}|\mathit{\theta };{\sigma }_p^{\mathrm{SI}},m_{\chi }\right)=\sum _i\left[D_{i}ln{N}_i(\mathit{\theta };{\sigma }_p^{\mathrm{SI}},m_{\chi })-N_i(\mathit{\theta };{\sigma }_p^{\mathrm{SI}},m_{\chi })\right].$$

(9)

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:

$$ln{\mathcal{L}}_{\mathrm{ext}.\mathrm{const}.}\left(\mathit{\theta }\right)=\sum _j\left[-\frac{1}{2}{\left(\frac{{\theta }_j-{\overline{\theta }}_j}{c_j{\overline{\theta }}_j}\right)}^2\right].$$

(10)

Here, ${\overline{\theta }}_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\mathcal{L}\left(\mathit{D}|\mathit{\theta };{\sigma }_p^{\mathrm{SI}},m_{\chi }\right)=ln{\mathcal{L}}_{\mathrm{Poisson}}\left(\mathit{D}|\mathit{\theta };{\sigma }_p^{\mathrm{SI}},m_{\chi }\right)+ln{\mathcal{L}}_{\mathrm{ext}.\mathrm{const}.}\left(\mathit{\theta }\right).$$

(11)

To calculate the projected sensitivity of paleo-detectors, we consider the maximum log-likelihood ratio [146,147,148]:

$$q({\sigma }_p^{\mathrm{SI}};m_{\chi })=-2ln\left[\frac{\mathcal{L}\left(\mathit{D}|\widehat{\widehat{\mathit{\theta }}};{\sigma }_p^{\mathrm{SI}},m_{\chi }\right)}{\mathcal{L}\left(\mathit{D}|\widehat{\mathit{\theta }};{\widehat{\sigma }}_p^{\mathrm{SI}},m_{\chi }\right)}\right].$$

(12)

In the numerator, $\widehat{\widehat{\mathit{\theta }}}$ is the set of nuisance parameters that maximizes $\mathcal{L}$ for given values of ${\sigma }_p^{\mathrm{SI}}$ (and $m_{\chi }$), while in the denominator, $\mathcal{L}$ is maximized over both $\mathit{\theta }$ and ${\sigma }_p^{\mathrm{SI}}$ (for a given value of $m_{\chi }$), with best-fit values $\widehat{\mathit{\theta }}$ and ${\widehat{\sigma }}_p^{\mathrm{SI}}$. In order to calculate the projected exclusion limit, we use the Asimov data set [146] under the assumption of no signal:

$$\mathit{D}=\mathit{N}(\overline{\mathit{\theta }};{\sigma }_p^{\mathrm{SI}}=0),$$

(13)

where $\overline{\mathit{\theta }}$ is a set of fiducial choices for the nuisance parameters. The 90% confidence level exclusion limit for a given value of $m_{\chi }$ is obtained by finding the smallest value of ${\sigma }_p^{\mathrm{SI}}$ for which $q({\sigma }_p^{\mathrm{SI}};m_{\chi })\ge q_{\mathrm{crit}}=2.71$. (By Wilks’ theorem [149], the distribution of q is (asymptotically) ${\chi }^2$ distributed, and for a one-dimensional ${\chi }^2$-distribution, ${\chi }^2=2.71$ for a p-value of 0.1.)

We describe the the fiducial values for the nuisance parameters listed in Equation (8) and the relative uncertainty of the external constraints on those parameters, $c_j$, defined by Equation (10). For the high-resolution (high-exposure) scenario, we choose a fiducial value of $\overline{M}=10$mg ($\overline{M}=100$g) for the sample mass, and for both scenarios, ${\overline{t}}_{\mathrm{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_{t_{\mathrm{age}}}=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, ${\overline{\mathsf{\Phi }}}_{\nu }^i$, and the ${}^{238}\mathrm{U}$ concentrations in our target minerals, ${\overline{C}}^{238}$, we take fiducial values described in Section 3 (the latter are also summarized in

Table 1. List of minerals considered in this work with their chemical composition and our fiducial assumption for the ${}^{238}\mathrm{U}$ concentration ($C^{238}$) in radiopure samples of these minerals.

Mineral	Composition	Fiducial ${}^{238}$U Concentration [per Weight, g/g]
Halite	NaCl	${10}^{-11}$
Gypsum	Ca(SO${}_4)·2$(H${}_2$O)	${10}^{-11}$
Sinjarite	CaCl${}_2·2($H${}_2$O)	${10}^{-11}$
Olivine	Mg${}_{1.6}$Fe${}_{0.4}^{2+}$(SiO${}_4$)	${10}^{-10}$
Phlogopite	KMg${}_3$AlSi${}_3$O${}_{10}$F(OH)	${10}^{-10}$
Nchwaningite	Mn${}_2^{2+}$SiO${}_3$(OH)${}_2·($H${}_2$O)	${10}^{-10}$

). As a benchmark, we assume that the associated nuisance parameters can be constrained to $c_{{\mathsf{\Phi }}_{\nu }^i}=100\%$ and $c_{C^{238}}=1\%$; these choices are equivalent to the respective systematic uncertainties of the backgrounds assumed in previous studies [35,36,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 ${\sigma }_x=1$nm; left panel) and the high-exposure ($M=100$g, ${\sigma }_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_{\chi }\lesssim 10$GeV/c${}^2$ where the dominant background 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_{\chi }\gtrsim 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}\mathrm{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}\mathrm{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,107,150] as well as the neutrino floor for a Xe-based experiment [151] with the shaded regions. At small WIMP masses, $m_{\chi }\lesssim 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_{\chi }=1$GeV/c${}^2$ the best existing upper limit is ${\sigma }_p^{\mathrm{SI}}\lesssim {10}^{-38}{\mathrm{cm}}^2$ [14,17]. Moreover, paleo-detectors offer sensitivity deep into the conventional (solar) neutrino floor in the $m_{\chi }\lesssim 10$GeV/c${}^2$ mass range. Due to the large exposures, a paleo-detector would observe $\mathcal{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 track-length (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_{\chi }\gtrsim 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 Equation (10)) on the parameters controlling the normalization of the neutrino-induced backgrounds (${\mathsf{\Phi }}_{\nu }^i$) and the radiogenic backgrounds ($C^{238}$). as well as removing the external constraint on the age of the samples ($t_{\mathrm{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_{\chi }\lesssim 1$GeV/c${}^2$ WIMP range when the constraints on the ${\mathsf{\Phi }}_{\nu }^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, Equation (9), with a ${\chi }^2$ distribution,

$$ln{\mathcal{L}}_{{\chi }^2}\left(\mathit{D}|\mathit{\theta };{\sigma }_p^{\mathrm{SI}},m_{\chi }\right)=\sum _i\left[-\frac{1}{2}{\left(\frac{D_i-N_i(\mathit{\theta };{\sigma }_p^{\mathrm{SI}},m_{\chi })}{{\sigma }_i}\right)}^2\right],$$

(14)

where ${\sigma }_i^2$ is the variance in the i-th bin. If we include only the Poisson error of the data, ${\sigma }_i^2\to D_i$, the Poisson log-likelihood and the ${\chi }^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 ${\chi }^2$ distribution allows us to include additional uncertainties, which we model as

$${\sigma }_i^2=D_i+{\left({\Delta }_{\mathrm{Bkg}}^{\mathrm{sys}}{D}_i\right)}^2,$$

(15)

where ${\Delta }_{\mathrm{Bkg}}^{\mathrm{sys}}$ parameterizes the (relative) systematic modeling error. We note that it is straightforward to extend this to ${\sigma }_i^2=D_i+{\sum }_j{\left({\Delta }_i^j{D}_i^j\right)}^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 $\Delta$ 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 ${\Delta }_{\mathrm{Bkg}}^{\mathrm{sys}}=0$, i.e., accounting only for the Poisson error in the variance, ${\sigma }_i^2=D_i$, both the Poisson log-likelihood and the ${\chi }^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 ${\Delta }_{\mathrm{Bkg}}^{\mathrm{sys}}=0.1$, corresponding to 10% bin-to-bin mismodeling of the shape of the backgrounds, the sensitivity does not depreciate significantly. For ${\Delta }_{\mathrm{Bkg}}^{\mathrm{sys}}=1$, the projected upper limit is up to an order of magnitude worse than in the ${\Delta }_{\mathrm{Bkg}}^{\mathrm{sys}}=0$ case for $m_{\chi }\lesssim 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, ${\Delta }_{\mathrm{Bkg}}^{\mathrm{sys}}=0.01$, we find that the sensitivity becomes only slightly worse than for ${\Delta }_{\mathrm{Bkg}}^{\mathrm{sys}}=0$, while for ${\Delta }_{\mathrm{Bkg}}^{\mathrm{sys}}=0.1$, the projected upper limit is a factor of ∼5 larger than for ${\Delta }_{\mathrm{Bkg}}^{\mathrm{sys}}=0$. The high-exposure-readout scenario is most relevant for the mass range $m_{\chi }\gtrsim 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}\mathrm{U}$ concentration. Thus, it seems unlikely that the systematic (bin-to-bin) shape errors would be as large as ${\Delta }_{\mathrm{Bkg}}^{\mathrm{sys}}=0.1$ in an actual paleo-detector experiment. However, we note from Figure 4 that even for errors as large as ${\Delta }_{\mathrm{Bkg}}^{\mathrm{sys}}=0.1$, a gypsum paleo-detector could still probe WIMP-nucleon cross sections almost an order of magnitude below current experimental limits for $m_{\chi }\gtrsim 100$GeV/c${}^2$.

In Figure 5, we show the effect of assuming different ${}^{238}\mathrm{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}\mathrm{U}$ concentrations, ${10}^{-12}\mathrm{g}/\mathrm{g}\le C^{238}\le {10}^{-7}$g/g shown in Figure 5 and most of the shown WIMP mass range, the projected upper limit scales like ${\sigma }_p^{\mathrm{SI}}\propto \sqrt{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 $\propto \sqrt{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_{\chi }\lesssim 10$GeV/c${}^2$ or, for very small ${}^{238}\mathrm{U}$ concentrations, atmospheric neutrinos for $m_{\chi }\gtrsim 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}\mathrm{U}$ concentrations as large as $C^{238}={10}^{-8}$g/g in the $m_{\chi }\lesssim 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 experimental limits even if the ${}^{238}\mathrm{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 ${\chi }^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.

5. Summary and Conclusions

This work surveys the prospects of searching for DM with paleo-detectors, natural minerals that 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,36,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 References [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}\mathrm{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_{\chi }\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.

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.