Nuclear decay anomalies as a signature of axion dark matter

A number of nuclear decay anomalies have been reported in the literature, which purport to show periodic variations in the decay rates of certain radioisotopes. If these reports reflect reality, they would necessitate a seismic shift in our understanding of fundamental physics. We provide the first mechanism to explain these findings, via the misalignment mechanism of QCD axion dark matter, wherein oscillations of the effective $\theta$ angle induce periodic variation in nuclear binding energies and hence decay rates. As we expect this effect to be most pronounced in low-$Q$ systems, we analyse 12 years of tritium decay data ($Q\simeq$ 18.6 keV) taken at the European Commission's Joint Research Centre. Finding no statistically significant excess, we exclude axion decay constants below $9.4\times10^{12} - 1.8\times10^{10}$ GeV (95% confidence level) for masses in the $1.7\times{10}^{-23} - 8.7\times 10^{-21}$ eV range.

Introduction.As a cornerstone of modern physics, it is widely accepted that radioactive decay is in general a truly random process, occurring independently of time, space and external influence.This simple fact carries with it a vast array of consequences, across fields as diverse as modern-day nuclear medicine to the Big Bang Nucleosynthesis which occurred in the first few minutes of our cosmic history.
Nonetheless, in recent decades a number of purported nuclear decay anomalies have been reported in the literature, which represent apparent violations of this rule [1,2].Generally these anomalies take the form of periodic variations in the observed decay rate, although notably there are some indications of temporary variations in nuclear decay correlated with solar flares and other astrophysical phenomena [3][4][5][6].
One economical explanation of these conclusions, supported by a body of evidence, is that they are largely the result of random noise, unaccounted-for systematic effects and incomplete uncertainty analysis [7][8][9][10][11][12][13][14][15].For example, the overall preference of some datasets for an annually modulated signal may be simply attributed to the influence of seasonally varying environmental conditions, rather than any exotic deviations from known physics [14].It should also be noted that although radioactive decay is of course a relatively well-explored and understood topic, accurately quantifying the various associated uncertainties is nonetheless nontrivial [15].This perspective may be supported by the fact that such anomalies are furthermore somewhat difficult to explain in the context of ordinary particle physics.Given the annual periodicity of some claimed signals, solar neutrinos are often suggested to be in some way responsible, but there is no known mechanism within the Standard Model which allows this [1].Some somewhat speculative hypotheses have been proposed, as summarised in [9,16], but at present there appears to be no concrete framework within which these anomalies can be understood consistently with other observations.Nonetheless, the absence so far of a compelling explanatory framework does not mean that none exist to be found.Furthermore, regardless of the provenance of these unusual observations, and in particular the key question of if they indeed represent a signature of new physics, if they can be connected to some favoured model of physics beyond the Standard Model, they can provide novel constraints and experimental strategies.
Bearing these motivations in mind we in the following provide the first mechanism to explain this phenomenon, via a class of particles known as axions.Such an approach carries with it the benefit of also addressing certain other problems in fundamental physics, namely the nature of dark matter (DM) and the Strong CP problem in the Standard Model.
Arising originally via a minimal extension of the Standard Model; the Peccei-Quinn solution of the Strong CP problem [17][18][19], axions and axion-like particles occupy a rare focal point in theoretical physics, in that they are simultaneously a generic prediction of the exotic physics of string and M-theory compactifications [20][21][22].Despite the profound differences between these two contexts the resulting axion properties are also largely universal, providing an easily-characterisable theoretical target.As light, long lived pseudoscalar particles they can also influence many aspects of cosmology and astrophysics, creating a wide variety of observational signatures [23].In particular, they are a natural candidate for the mysterious DM comprising much of the mass of our visible universe [24,25], and as such are a topic of intense ongoing investigation [26].
The aforementioned Strong CP problem is the question of why the effective QCD θ parameter, which enters via the Lagrangian term is smaller than 10 −10 in absolute value and not O(1), as would be expected on the grounds of naturalness.Here g is the gauge coupling and G µν the gluon field strength tensor.Axions which enact the Peccei-Quinn solution of the strong CP problem do so by promoting θ to a dynamical variable, which can then relax to zero as required.
In the event that these axions also comprise the DM in our Universe, the misalignment mechanism ensures that the time-dependent θ angle today is where ρ DM is the DM density, ω = m a (1 )), p, v, f a and m a are the axion field momentum, velocity, decay constant and mass respectively, and ϕ is an arbitrary phase.
Since various aspects of nuclear physics are θdependent this has the potential to lead to a variety of observational signatures, from which constraints on the axion parameter space can then be placed.For example, Refs.[27][28][29][30] search for an oscillating neutron electric dipole moment (nEDM), and in Ref. [31] it was also demonstrated that an oscillating θ-angle can lead to underproduction of 4 He during Big Bang Nucleosynthesis (BBN).
For our purposes it is the θ-dependence of nuclear binding energy that is of primary interest.Once nuclear binding energies become time-dependent via (2) we can expect a periodic variation in the nuclear decay rates and hence the possibility to explain the reported decay anomalies.
In particular, we will demonstrate that tritium represents a particularly opportune target to search for these effects.This established, we will then use existing 3 H data to search for the corresponding axioninduced signal.Finding no statistically significant excess, we can then exclude axion decay constants below 9.4 × 10 12 − 1.8 × 10 10 GeV (95 % confidence level) for masses in the 1.7 × 10 −23 − 8.7 × 10 −21 eV range.Conclusions and discussion are presented in closing.
θ-dependence of nuclear decay rates.Our underlying quantity of interest is the fractional change in the beta decay rate Γ(θ) as a function of θ, where E max = (M 2 i +m 2 e −(M f +m ν ) 2 )/2M i , with initial and final state masses M i/f , neutrino and electron masses m ν/e , E e the energy of the emitted electron, and δE(θ) an additional θ-dependent contribution.
In principle this θ-dependence can also enter in other ways, for example modification of the underlying nuclear couplings g A/V .However, in line with the results of Ref. [32] it is reasonable to assume the leading order contributions arise specifically from modifications of the phase space.This point is also specifically analysed in the Supplementary Material [33], in support of this conclusion.
By approximating nuclear beta decay in terms of free neutron decay, it has also been argued in Ref. [2] that for small perturbations to the decay rate, decays with smaller Q ≡ M i − M f − m e resulted in a larger fractional change in the beta decay rate.
However, from the beta decay rates given in Ref. [34][35][36] we can in any case evaluate this integral numerically for small δE(θ) without the free-neutron approximation, finding for 187 Re (Q ≃ 2.6 keV) and 3 H (Q ≃ 18.6 keV) that (4) This result is particularly fortuitous in that if we wish to search for this effect, high quality datasets for 3 H already exist.The use of 3 H decay as a window to new physics more generally has also recently been explored in Ref. [37][38][39].Furthermore, in addition to being more sensitive to the effect of a time-varying θ angle, the θdependence of the properties of lighter nuclei such as 3 H may also be comparatively easier to discern.
Indeed, in Ref. [32] the dependence of the nuclear binding energy upon θ has already been calculated for light nuclei.It was estimated that for three and four-nucleon systems the n-nucleon binding energy B n satisfies Here the bar indicates an average over states which become degenerate in the limit that the approximate Wigner SU (4) symmetry of low energy nuclear physics becomes exact.As noted in Ref. [32] this agrees with the results calculated numerically in Ref. [40].
The averaged 2-nucleon binding energy B 2 (θ) can be found by averaging over the physical deuteron and the spin singlet (dineutron and diproton) channel, with θ dependence parameterised via where B 2 (0) = 2.22, −0.072, −0.787 MeV for the deuteron, dineutron and diproton respectively, and the c i are numerical coefficients given in Ref. [32] (herein we choose the more conservative parameter set II, assuming isospin conservation).
Since this only gives the (averaged) B 3 (θ), rather than the individual initial and final state binding energies B(θ) i/f , we assume in the following that for small θ which is equivalent to assuming that for small θ the individual binding energies scale in proportion to their average.
Knowing the θ-dependence of the binding energy, we can then calculate δM i/f , the θ-dependent part of M i/f , via where the sum runs over nucleons, and the neutron/proton mass difference is given by where , and we assume two degenerate quark flavours, following Ref.[32,41] Expanding about θ = 0 we find which in turn can be calculated via equations ( 7) and ( 9) since we know the θ-dependent parts of (m n − m p ) and (B i − B f ).In turn this allows us to calculate I 0 (θ).However, care is required here.In a Universe where this mechanism is active the measured values of physical quantities such as decay rates are actually not those at θ = 0, but are instead those at ⟨θ 2 ⟩ (since binding energies are O(θ 2 ) at leading order).Therefore when dealing with experimental data, rather than comparing Γ(θ) to Γ(0), we should instead consider where ⟨Γ⟩ is the average value of Γ(θ).Using the known I 0 (θ) = Γ(θ)/Γ(0) − 1, we find after some calculation which oscillates about zero, as expected, leading to alternating excesses and deficits in the number of nuclear decays per unit time.We can also see that the corresponding shift in the decay energy, is also a small correction in regions of the axion parameter space we exclude, justifying our assumptions that δE << Q.Similarly, θ << 1 in these regions.Data analysis and results.These points established, we can now search for evidence of this effect in 3 H.We make use of a dataset, shown in Fig. 1, provided by the European Commission's Joint Research Centre (JRC) at the Directorate for Nuclear Safety and Security in Belgium, which spans approximately 12 years of liquid scintillation counter observations of the decay of an O(microcurie) 3 H source [9].
Here Ṅ is the measured value of decays per second, whilst ⟨ Ṅ ⟩ is its expected value due to the exponential decay law.
It is worth emphasising that the statistical analysis of nuclear decay anomalies has been subject to a number of treatments [42], with no overall consensus on which approach is optimal.This being the case we broadly follow the approach of Refs.[29,30,43], where searches for oscillatory signals in time-series data were also used to place limits on ultralight axion DM.
Since the data points are unevenly spaced in time we will estimate their power spectrum using the Least Squares Spectral Analysis (LSSA) method to construct periodograms [44,45].We compute the power spectrum using the astropy.timeseries.LombScargle class provided by the Python astropy.timeseriespackage [46], evaluated at a set of 8113 evenly spaced trial frequencies.
As the lowest 15 trial frequencies appear to show evidence of uncontrolled systematic effects, possibly longterm drift of the experimental apparatus, we exclude them from our analysis.The largest frequency in the remaining data, 4.2 × 10 −6 Hz, corresponds to a period of ∼ 2.8 days, whilst the smallest, 8.0 × 10 −9 Hz, corresponds to a period of ∼ 4.0 years.The resulting periodogram is then given by the blue curve in Fig. 2.
Under the null hypothesis that the dataset contains no axion-induced signal, the time-series datapoints should follow a Gaussian distribution about I = 0. To place limits on the corresponding power spectrum we there-fore perform Monte Carlo (MC) simulations by generating N = 50000 time-series MC datasets, with the same time spacing as the original data.The MC data points themselves are drawn from a zero-mean Gaussian distribution, with width set by the standard deviation of the original dataset.
For each MC dataset we can calculate a corresponding periodogram, and then use the statistics of these periodograms to construct the Cumulative Distribution Function (CDF) for the power at each frequency.From these CDFs we then determine the false positive (or false alarm) power at 95% confidence level for each frequency, as shown in Fig. 2 in orange.Here and in the following we account for the 'look elsewhere effect' by defining these limits with respect to the global trials factor , where p local is the corresponding local p-value and N f the number of trial frequencies.As can be seen, nowhere does the original dataset exceed the 95 % confidence level threshold, indicating the data are compatible with the null hypothesis.As can be easily seen, the data are compatible with the null hypothesis.
To determine the corresponding limit on the axion parameters we follow a similar approach.For each trial frequency we construct N = 50000 MC datasets containing Gaussian background and injected axion-derived signals and calculate their periodograms.The injected signals are of fixed amplitude and frequency, with the unknown axion phase ϕ drawn from a uniform distribution.As the injected signals are constructed in time series form, we directly match the axion velocity to the lab-frame DM velocity, incorporating modulation effects due to the Earth's passage through the DM halo [47,48].For a given choice of parameters, we can then construct the corresponding CDF for the power at each trial frequency.With the mass fixed by the frequency under consideration, the threshold value of f a can then be determined following a standard frequentist approach, well illustrated in Ref. [49].Having determined the threshold value of f a within this framework, we then subsequently account for the 'stochastic vs deterministic' correction factor occurring in the regime where the measurement time is much less than axion field coherence time, following Ref.[49].
Specifically, from the background-only CDF we first find the false positive threshold power at a desired confidence level.From the background plus signal CDF we can also find the false negative threshold power, for a given choice of f a .The threshold value of f a is then determined by the condition that the false positive threshold from the background-only CDF coincides with the false negative threshold from the background plus signal CDF, at the desired confidence level.Equivalently we can say that the threshold value of f a occurs when the false positive rate α is equal to the false negative rate 1 − β, where CL = 1 − α.The resulting exclusion curve is given in Fig. 3 [27][28][29][30], BBN [31], the spectroscopy of radio-frequency atomic transitions [51], pulsars, gravitational waves and black hole superradiance [52][53][54][55][56], solar/white dwarf observations [57,58], so-called 'axinovae' [59], and the parameter space occupied by canonical QCD axion models (yellow).Figure producing using the AxionLimits code [60].
We simulate future experimental possibilities, via two fiducial cases also shown in Fig. 3.The first is a longterm experiment designed to search lower frequencies, with 1 measurement per hour over 12 years.The second is a short-term experiment to search higher frequencies, with approximately 1 measurement per second over 1 day.Both schemes increase the initial quantity of tritium by a factor of 100 relative to the JRC dataset, up to a presumed limit imposed by detector pileup [50].In the first case we can exclude axion decay constants up to f a < 1.5 × 10 16 GeV (95 % confidence level) and cover masses in the 5.4 × 10 −25 − 1.4 × 10 −17 eV range, while in the latter we can exclude axion decay constants up to f a < 4.5 × 10 11 GeV (95 % confidence level) and cover masses in the 2.4 × 10 −21 − 6.3 × 10 −14 eV range.
Conclusions/discussion. We have explored purported variations in nuclear decay rates, providing a novel and well-motivated explanatory mechanism within the framework of axion physics.The oscillating QCD θangle created by the misalignment mechanism induces a time-varying nuclear binding energy, which then leads to the time-dependence of nuclear decays.
One advantage of this approach is that we can probe regions of axion parameter space which normally are only accessible experimentally via much more sophisticated methods, such as oscillating nEDM searches [27][28][29][30].In contrast, the data we have analysed here were taken using a single O(microcurie) tritium source and a commercially available laboratory liquid scintillation counter.Of course, there are also astrophysical and cosmological bounds on the mass of axion dark matter, such as m a > 3 × 10 −19 eV (99 % confidence level), coming from observations of ultra-faint dwarf galaxies [61] and m a < 2 × 10 −20 eV (95 % confidence level)bounds from the Lyman-alpha forest [62].However, our work provides a new kind of laboratory-based scheme to search for QCD axion dark matter, which is complementary to such constraints without relying on cosmological or astrophysical assumptions and modelling.
However, relying on nuclear decays in this way also presents certain challenges.In particular, issues such as detector pileup may ultimately constrain sensitivity more so than in other experimental approaches.
Another overall consideration here is the validity of this and other approaches in searching for QCD axionderived phenomena far away from the parameter space where these effects are expected to occur: the so-called QCD axion band.It is first of all important to note that (admittedly non-minimal) models of the QCD axion do exist in the parameter space we have probed here, as given in Refs.[57,[63][64][65].Furthermore, we can also understand our efforts here as being a necessary step to ultimately reach the sensitivity to probe standard QCD axion models via methods such as these.
There are several avenues by which these findings could be improved.On the theoretical side, more accurate estimation of the dependence of nuclear binding energies on θ would be of use.Experimentally, dedicated nuclear decay experiments which are optimised to search for this particular effect should also yield stronger constraints, and may be able to probe previously unexplored regions of the axion parameter space.Analysis of other preexisting nuclear decay datasets may also similarly bear fruit.
Nuclear decay anomalies as a signature of axion dark matter Supplementary Material Xin Zhang, Nick Houston, Tianjun Li In this supplementary material we derive in detail the θ-dependence of the nuclear half life of 3 H.This radioisotope is of particular interest due to the low Q-value (18.6 keV), which may result in greater sensitivity to the extremely small axion-induced energy shifts we are searching for.Conveniently 3 H beta decay has for this reason also been studied intensively, due to the relevance it may have for cosmic neutrino background searches.
Our starting point is the Hamiltonian density where is the Cabbibo angle, and e/ν e are respectively the electron and neutrino fields.The strangeness-conserving free nucleon charged current is where p/n are the proton and neutron fields, and g V /g A are the nuclear vector and axial-vector couplings respectively.Following the detailed calculations in Ref. [34], the total β-decay rate is where the Fermi and Gamow-Teller beta strengths are respectively and the phase space integral is where E max is the maximum possible electron energy, p e = |⃗ p e | is the magnitude of the electron 3-momentum, and F 0 is the Fermi function accounting for modification of electron wavefunctions due to the Coulomb field of the nucleus.We will in the following assume isospin conservation, since isospin-violating effects and the θ-dependence we seek to understand are both individually small corrections in our regime of interest.These points established, we can assess where θ-dependence may conceivably enter.Non-QCD parameters.It is evident first of all that parameters which belong to non-QCD sectors, such as m e , G F , and θ C are θ-independent.The PMNS matrix elements responsible for neutrino oscillations and the neutrino masses themselves should also presumably have no such dependence.Nucleon masses.Conversely, it is already known that for free nucleons (m n − m p ) increases with increasing θ, where (m n − m p ) QED = −(0.58± 0.16)MeV is the QED correction to the mass difference [41].We can then write where c 5 = (−0.074± 0.006)GeV −1 is a low energy constant, B 0 = m 2 π /(m u + m d ), and we assume two degenerate quark flavours, following Ref.[32].Nuclear binding energies.We can also expect nuclear binding energies to be θ-dependent.To calculate this dependence we can make use of the results in Ref. [32], where it was empirically estimated that for n = 3, 4 nucleon systems the binding energies B n satisfy Here the bar indicates an average over states which become degenerate in the limit that the approximate Wigner SU (4) symmetry of low energy nuclear physics becomes exact.
Rearranging for n = 3 yields where averaging over the physical 3 H and 3 He states gives B 3 (0) ≃ 8.1 MeV, and B 4 (0) ≃ 28.3 MeV is given by the physical 4 He binding energy [32].
To find B 2 (θ), we average over the physical deuteron and the spin singlet (dineutron and diproton) channel, with θ dependence parameterised via with B 2 (0) ≃ 2.22, −0.072, −0.787 MeV for the deuteron, dineutron and diproton respectively, following Ref.[32], which then yields B 2 (0) ≃ 1.03 MeV.For the c i numerical coefficients we use the (more conservative) choice of parameter set II in Ref. [32], assuming isospin conservation, (c 1 = 3.25, c 2 = 2.55, c 3 = 0.47).Inserting these values we find Since we only know the (averaged) B 3 (θ), rather than the individual initial and final state binding energies B i (θ) and B f (θ), we will assume in the following that for small θ B(θ which is equivalent to assuming that for small θ the individual binding energies scale in proportion to their average.Initial/final state energies.Knowing the dependence of free nucleon masses and nuclear binding energies then allows us to calculate M i and M f , the masses of the initial and final nuclear states.This is defined in terms of the sum of the masses of their constituent (free) nucleons N and their associated binding energies, Combining the previous results we then find Integral endpoint.Modifications to M i/f enter into the beta decay rate via the phase space integral in Eq. ( 18), through the upper limit E max , where where m ν is the neutrino mass.To calculate the leading order modification to E max we can perturb M i/f by δM i/f (representing a θ-dependent correction) and expand to lowest order, finding Since M i/f >> m e , m ν we can further neglect terms proportional to products of δM i/f and m e /m ν to then find The latter terms can be calculated exactly here, but we also can notice that since Q = M i − M f − m e ≃ 18.6 keV is also much less than M i/f , we then have which then leads simply to From the previous section we can then calculate: Fermi function.Before calculating the phase space integral we need to examine the relativistic Fermi function [34][35][36], given by where k = 1, 2, 3 . . ., Z is the atomic number, γ k = k 2 − (αZ) 2 , R 0 = 1.2 × 10 −15 m is the nuclear radius, and z = αZE e /p e .In principle θ-dependence may enter via R 0 : as nucleons become more tightly bound we can expect the effective nuclear radius to decrease.However, given that it is straightforward to see that modification of the Fermi function need not be considered.Phase space integral.These points established we can now add the small perturbation δE to E max and evaluate the phase space integral for a given choice of Fermi function: 2 dE e Emax(0) me where: R is the average radius of a nucleus, a good approximation for the average radius of a nucleus with A nucleons is In the Primakoff-Rosen (non-relativistic) approximation [36] where F PR 0 (Z, E e ) = 2πy/[1 − exp(−2παZ)], y = αZE e /p e .We can find the result exactly, expanding for small δE/E max to give Numerical integration and expanding around δE = 0 then yields which agrees to within 12% with the non-relativistic case.Beta strengths.Following Ref. [34][35][36], on the basis on isospin symmetry and selection rules it can also be shown that If we assume isospin symmetry, these results should carry over to the present context, as selection rules are presumably not θ-dependent for small values of θ.The couplings g A and g V may also have some θ-dependence, however following Ref.[32] it can also be assumed that their θ-dependence is subleading to phase space modifications.
Fractional decay rate.These points established we can now combine these results to establish our quantity of interest for comparison to experimental data, the fractional change in the beta decay rate Γ as a function of θ, where we have used the fact that θ-dependence primarily enters through modification of the phase space integral, and dependence on the various other quantities present in Γ cancels.
In the axion dark matter scenario the form of θ is where ρ DM is the DM density, ω = m a (1 + v 2 + O(v 4 )), p, v, f a and m a are the axion field momentum, velocity, decay constant and mass respectively, and ϕ is an arbitrary phase.This then yields

FIG. 2 .
FIG.2.Periodogram corresponding to the dataset shown in Fig.1(blue), along with the MC-derived 95% confidence level threshold (orange).As can be easily seen, the data are compatible with the null hypothesis.