Measuring Changes in the Atmospheric Neutrino Rate Over Gigayear Timescales

Measuring the cosmic ray flux over timescales comparable to the age of the solar system, $\sim 4.5\,$Gyr, could provide a new window on the history of the Earth, the solar system, and even our galaxy. We present a technique to indirectly measure the rate of cosmic rays as a function of time using the imprints of atmospheric neutrinos in paleo-detectors, natural minerals which record damage tracks from nuclear recoils. Minerals commonly found on Earth are $\lesssim 1\,$Gyr old, providing the ability to look back across cosmic ray history on timescales of the same order as the age of the solar system. Given a collection of differently aged samples dated with reasonable accuracy, this technique is particularly well-suited to measuring historical changes in the cosmic ray flux at Earth and is broadly applicable in astrophysics and geophysics.

The sources, composition, and propagation of cosmic rays through our galaxy remain some of the biggest mysteries in astrophysics despite a century of experimental effort using a variety of techniques. Today, balloonand space-based experiments measure the flux of cosmic rays directly, while terrestrial experiments probe the nature of cosmic rays indirectly via measurements of secondary particles created by their interactions with the Earth's atmosphere. Notably, atmospheric neutrinos arising from such cascades can be detected in terrestrial experiments [1][2][3][4][5].
Characterizing atmospheric neutrinos provides a window on a variety of physical processes, many of which cannot easily be studied by other means. The flux of primary cosmic rays incident on the atmosphere depends on the sources and composition of the cosmic rays themselves as well as their subsequent propagation through the cosmos. Upon reaching our solar system, the cosmic ray trajectories are further altered by the magnetic fields of the Sun and Earth [6][7][8][9]. Finally, the composition and structure of the Earth's atmosphere impact the observed neutrino flux at the surface [7,10,11].
Compared to other tracers of the primary cosmic ray flux, e.g. cosmogenic muons, atmospheric neutrinos have the advantage that they are not screened by the Earth. The cosmogenic muon flux is exponentially sensitive to the height of the overburden; thus, using muon-induced, rather than neutrino-induced recoils to infer the cosmic ray flux with paleo-detectors would require exquisite knowledge of the sample's geological history. Using atmospheric neutrinos as tracers, the required geological knowledge is only that samples have been buried at depths large enough for cosmogenic muons to be well shielded. Paleo-detector measurements of the atmospheric neutrino flux can thus be sensitive to the evolution of the cosmic ray flux over geological timescales.
To date, the evolution has only been measured by studying rare isotopes produced in cosmic ray interactions. For example, the deposition of cosmogenic nuclides, such as 10 Be , in the Earth's crust allow us to track the evolution of the cosmic ray flux and Earth's atmosphere going back as far as ∼ 10 Myr [42,43]. Cosmogenic nuclides in meteorites can include radiogenic isotopes with much longer half-lives, (e.g. 40 K with T 1/2 ∼ 1.2 Gyr), providing information about the galactic cosmic ray flux on gigayear timescales (for a review, see [43]). However, such studies infer the flux of cosmic rays at a given meteorite's location and are hence not applicable to studying local influences on the cosmic ray flux at Earth.
Here, we consider the capability of paleo-detectors to differentiate several generic scenarios of cosmic ray flux evolution without regard for the specific underlying physics -such changes could be gradual (e.g. originating from slowly evolving star formation or supernova rates) or transient (e.g. caused by a nearby supernova or merger of neutron stars). For definiteness, we will consider halite (NaCl) samples, but our conclusions hold for other target materials. With the present-day flux, approximately 6×10 4 tracks from nuclear recoils created by atmospheric neutrino interactions are expected in the background-free signal region of a 1 Gyr old, 100 g sample. Given this large number of signal tracks and the ability to control systematics with multiple mineral samples, we anticipate that paleo-detectors could be highly sensitive to changes in the cosmic ray flux over ∼ 1 Gyr.
Atmospheric Neutrino Signal.-The atmospheric neutrino flux depends on a variety of physical systems which all evolve independently over geological timescales. Cosmogenic nuclide studies of meteorites suggest that the primary cosmic ray intensity in our galaxy has increased by a factor of ∼ 1.5 over the last ∼ 1 Gyr [43]. While the propagation of cosmic rays through our solar system depends on the strength and orientation of the heliomagnetic field, such effects are modulated much faster than the ∼ 100 Myr geological timescales relevant for paleodetectors (see [9] and references within).
Cosmic rays are also deflected by the geomagnetic field before interacting with the Earth's atmosphere. The geomagnetic field is well-approximated by a dipole and studies of the geodynamo and paleomagnetic records (for example, see [8]) suggest that the dipole moment could have varied by an O(1) factor over Gyr timescales. The main effect of the geomagnetic field on the primary cosmic ray spectrum hitting the atmosphere is a rigidity cutoff, providing a lower-energy limit on primary cosmic rays, and in turn, the associated atmospheric neutrino spectra. The rigidity cutoff depends on the location of the target mineral relative to the orientation of the geomagnetic fields and is directly proportional to the dipole moment (for example, see [6]); we will discuss the implications for paleo-detectors below.
Finally, the atmospheric neutrino flux depends on the composition and density of the atmosphere. The change in atmospheric composition over the relevant timescales 1 Gyr considered here is primarily the replacement of a fraction of the N 2 with O 2 (for example, see [44]), which does not significantly alter the cascades that yield atmospheric neutrinos. A variety of studies (see [11] and references within) suggest that the density of the atmosphere could have varied by as much as a factor of two downwards or an order of magnitude upwards going back O(1) Gyr. Since the atmosphere today is thicker than ten interaction lengths of a typical cosmic ray proton [7,10], we expect no significant modification of the atmospheric neutrino flux from such density changes.
Before assessing the sensitivity of paleo-detectors to time-varying effects, we discuss the modeling of the atmospheric neutrino signal. For simplicity, we assume that the atmospheric neutrino spectrum is constant in time. Changes in the spectral shape would affect the expected signal, but for the simple counting experiment we propose, these changes just alter the observed rate in the signal regions and hence would be largely indistinguishable from a change in overall normalization. A spectral analysis could exploit the track length distribution's dependence on the shape of the atmospheric neutrino flux; however, the nuclear recoil spectrum is only weakly dependent on the primary cosmic ray spectrum because of the chain of intermediate interactions.
We calculate the neutrino flux with the multi-particle transport and interaction code FLUKA [45][46][47] as in [10] for minimum solar modulation, at the location of the Gran Sasso underground laboratory (42.44 We consider all nuclear recoil tracks stemming from interactions of atmospheric neutrinos with the target nuclei to be signal. Signal tracks arise from the constituent nuclei of halite (NaCl) as well as lighter nuclei produced by inelastic neutrino interactions with the target. Secondary nuclear recoils can be induced by products of neutrino interactions (in particular neutrons, protons, pions, muons, and kaons) with the target. We model the signal with FLUKA, which implements quasi-elastic, resonant (∆ production only), and deep inelastic scattering (DIS) interactions along with initial-and final-state effects.
Accurate modeling of neutrons in the target material is crucial to predicting both the signal component from neutrons produced in neutrino interactions and the background from radiogenic neutrons (discussed below). Details on the treatment of neutrons in previous versions of FLUKA can be found in [46]. For this work, an unreleased version of FLUKA containing new algorithms for the interactions of low energy neutrons is used. The new treatment ensures energy conservation for each interaction and provides an estimate of the recoil energy. These algorithms are still under development and are thus only available for a subset of isotopes. In particular, the cross sections for Cl are not implemented and so a similar nucleus (P) is substituted.
We generate neutrino interactions at the center of a large (20 m) 3 uniform cube to ensure that all products from atmospheric neutrino interactions are contained in the simulated volume. In the left panel of Fig. 1, we show the distribution of atomic numbers and energies of recoils induced by atmospheric neutrinos in a halite target (including ν e ,ν e , ν µ , andν µ for both neutral and charged current interactions). The track lengths of the recoils associated with both the signal and the backgrounds are calculated using stopping powers from SRIM [48,49]; we show the resulting track length spectrum in the right panel of Fig. 1.
We assume the same read-out benchmark previously used in studies of paleo-detectors for dark matter and core collapse supernova neutrinos. By imaging with small angle X-ray scattering at a synchrotron facility, we estimate that a ∼ 100 g sample could be read out with threedimensional spatial resolution of ∼ 15 nm [50,51]. We note that read-out techniques with lower spatial resolution and higher throughput could potentially be used to image the higher energy recoil tracks investigated in this work (see [17] for discussion).
Radiogenic Backgrounds.-A large number of background tracks from α-particles produced in the decay chains of radioactive contaminants such as 238 U could potentially impact the sensitivity of paleo-detectors to atmospheric neutrinos. However, very light nuclei (A < 5, i.e. H and He) generally have too small stopping powers to leave robust tracks in minerals, although some recent work has shown that it is possible to read out α-particle tracks in some materials [20,21]. Here, we assume that recoiling nuclei with A < 5 do not leave resolvable tracks. For materials and readout technologies where such tracks are visible, the qualitative results of this paper will remain the same.
The most important remaining backgrounds from 238 U come in two forms: 1) neutrons produced in radioactive decays which scatter off of target nuclei and 2) daughter nuclei from 238 U spontaneous fission (SF) which are generally much heavier and more energetic than the nuclei created by atmospheric neutrino interactions with halite. Note that the tracks from SF daughters can be used for fission track dating of the samples.
We calculate the primary neutron spectra from the en-tire 238 U → . . . → 206 Pb decay chain with SOURCES [52]. Our calculation includes the neutrons produced by 238 U SF (and the nuclei in its decay chain) as well as neutrons produced in (α, n) reactions of 238 U decay chain α-particles with Na and Cl. These neutrons lose their energy predominantly via elastic interactions with the nuclei comprising the target material, giving rise to a large number of relatively soft nuclear recoils. We propagate the neutrons through the material and calculate the neutron-induced recoil spectrum with FLUKA; the associated tracks dominate the track length spectrum at lengths 1 µm.
Daughters from 238 U SF are modeled with FREYA [53], which produces correlated fission secondaries using a combination of data and analytical models. The daughter nuclei from each fission event come out approximately back-to-back, so we treat the two daughter tracks as a single longer track. We conservatively assume that the chosen readout technique is insensitive to the distribution of energy loss along each track although this would in principle be useful for distinguishing different recoiling nuclei and identifying the two fission-daughter tracks. The SF daughters produce tracks with characteristic length ∼ 25 − 40 µm in halite.
Sensitivity to Atmospheric Neutrinos.-In previous studies of paleo-detectors [15][16][17][18], which investigated the sensitivity to signals associated with keV-scale recoils, the dominant background was due to radiogenic neutrons. In Fig. 1, we see that atmospheric neutrino interactions are at a significantly higher energy scale and can create tracks much longer than those arising from radio-genic neutrons. Consequently, the level of radiopurity in the sample does not strongly affect the sensitivity of the proposed measurement and the qualitative conclusions presented below should hold for similar minerals which can form and retain tracks on geological timescales.
Mineral samples of different ages should be associated with the same type of host rocks in order to help control for any potential systematics arising from different characteristic geological histories. For instance, the contribution of non-neutrino cosmic rays close to the surface (e.g. muons and neutrons) could vary significantly. Changes to the depth of host rocks with time could raise or lower the cosmogenic muon/neutron flux. As discussed in [15,16], many halite-bearing evaporite deposits feature salt structures which become accessible to bore-hole drilling after intruding into the overlying rock. During this intrusion (called diapirism) the halite typically re-crystallizes, effectively erasing any latent tracks induced by cosmic rays during the initial burial period. We therefore simply consider cosmogenic muons/neutrons as a background which can be avoided by retrieving samples with a variety of ages from evaporite deposits at sufficient depths 5 km in different locations (for example, see [54,55]).
Systematics associated with flux variations as a function of either latitude or time arising from (necessarily) excavating different-age samples at different locations can potentially be constrained or measured in the background-free region of the track length spectrum at shorter track lengths. The tracks in the 2 µm ≤ x ≤ 20 µm signal region primarily arise from nuclei with energies 10 MeV and masses close to or matching those of the nuclei comprising the target mineral. Such recoils typically arise in (quasi)elastic scattering of neutrinos or secondary neutrons off of the target nuclei. The rigidities of cosmic ray particles capable of producing atmospheric neutrinos at energies where (quasi)elastic scattering dominates the interactions with the target are similar to the range of geomagnetic rigidity cutoffs one might expect at different latitudes today.
The signal tracks in the background-free region at longer track lengths, 50 µm ≤ x ≤ 1 mm, arise exclusively from nuclei with masses much lighter than those comprising the target. Such lighter nuclei are produced by DIS of neutrinos with the target nuclei. Cosmic rays giving rise to atmospheric neutrinos sufficiently energetic to interact with the target via DIS have energies well above the largest rigidity cutoff one expects at any latitude, even after accounting for a potential O(1) variation in the geomagnetic field strength over the last Gyr. This signal region is thus largely independent of the systematics associated with the geomagnetic field.
Here, we do not include effects of annealing or track fading. While detailed studies are required to determine a precise model for track formation and retention in the variety of possible paleo-detector minerals, tracks associated with the much higher energy nuclear recoils induced . Error bands represent the 1σ error one would obtain on the average recoil rate from a sample of a given age, assuming that the age of samples can be determined with relative uncertainty of 5% and its mass with 1% uncertainty.
by atmospheric neutrinos should be more robust to annealing compared to tracks induced by dark matter or neutrinos from core collapse supernovae. Such studies have been carried out for recoils induced by SF daughters [56,57], which produce tracks similar to the signal. Any effects of track fading measured in an analysis of the fission tracks could, in principle, provide for a template to model the fading of signal tracks. In any case, signal and background track fading should be similar and, thus, large signal-rich, background-free swaths of the track length spectrum will remain.
Measuring the time evolution of the atmospheric neutrino flux requires counting the number of atmospheric neutrino induced tracks N i in (the background free regions of) multiple mineral samples of different ages t i . A simple quantity one could construct from such measurements is the average recoil rate per unit target mass over the age of the sample, r i = N i /t i M i , where M i is the sample mass. Given the large number of tracks expected in a representative sample (see Fig. 1), a good estimate for the relative error on r i can be obtained by summing the relative systematic errors of the quantities entering the rate in quadrature. One source of systematic errors is the relation of the measured number of tracks N i and the true number of recoils induced by atmospheric neutrino interactions. Such errors could e.g. be induced via mis-modeling of backgrounds, track fading, or read-out efficiencies. Further systematic errors are due to imperfect determination of the samples' ages and masses.
In Fig. 2 we show a cartoon of the average recoil rate r as a function of the age of a sample for three different scenarios: constant atmospheric neutrino flux (equal to the present-day flux), linearly decreasing flux with a 50% drop over 1 Gyr, and a transient event where the flux is elevated by 100% between 450 Myr and 550 Myr ago. To get an idea for the uncertainty on r, we assume that the samples' ages and masses can be determined with relative precision of ∆ t = 5% and ∆ M = 1%, and that such errors dominate over systematic errors on N i . The shaded bands indicate the corresponding error estimate. Quantitative statements about differentiating various scenarios for the atmospheric neutrino flux history are dependent on the specific parameters of the excavated samples (M i , t i ) and the statistical procedure used for a particular analysis, cf. [15]. However, Fig. 2 indicates that with the above error assumptions, a modest number of samples with different ages covering t i ∼ 100 Myr − 1 Gyr would allow us to decisively differentiate scenarios which differ at least as much from each other as those shown in Fig. 2.
Conclusions.-Paleo-detectors are a unique tool to study the evolution of the cosmic ray flux at Earth over gigayear timescales using the traces of nuclear recoils induced by atmospheric neutrino interactions. Due to the large exposures that can be achieved using natural minerals (∼ 0.1 kg Gyr) and the presence of backgroundfree track length regions, paleo-detectors are well-suited to measuring changes in the atmospheric neutrino flux over time. In a simple analysis which counts the recoil tracks in two background-free track length regions (e.g. 2 µm ≤ x ≤ 20 µm and 50 µm ≤ x ≤ 1 mm in halite), the total number of recoils expected in a 100 g, 1 Gyr old sample is ∼ 6 × 10 4 .
We note that while a simple counting experiment is unable to isolate the effects of the different physical mechanisms affecting the flux, progress towards this end is possible with more sophisticated analysis techniques. Combining measurements from a modest sample of ancient minerals of varying ages allows one to distinguish scenarios of cosmic ray flux evolution and constrain systematic uncertainties due to differences in sample composition, age, and location. These measurements offer a novel approach to studying the history of cosmic rays at Earth over gigayear timescales. Paleo-detectors could thus open a new window into the evolution of the Earth, in addition to providing for a more complete picture of astrophysics in our galaxy.