Geoneutrino flux to TNU conversion¶
Ondřej Šrámek, Charles University, ondrej.sramek@gmail.com, 27 Sep 2018
Question addressed: Given a geoneutrino flux at a point, what is the equivalent TNU signal?
Earth's natural antineutrinos detectable by inverse beta decay (IBD) come from two decay chains: 232Th and 238U.
The antineutrino production spectra $dn_\nu/dE_\nu$ were calculated by Enomoto: http://www.awa.tohoku.ac.jp/~sanshiro/research/geoneutrino/spectrum/; the spectral data files can be downloaded, the site also lists the assumptions made in the calculation. The spectra are normalized to the number of antineutrinos emitted per decay of one atom of parent nuclide (i.e., $n_\nu=4$ for 232Th and $n_\nu=6$ for 238U).
The geoneutrino (oscillated) flux spectrum $d\phi/dE_\nu$ at location $r$ is calculated from (don't need to worry now about the complete notation)
$$\frac{d\phi}{dE_\nu}(E_\nu,r) = \frac{Xl}{n_\nu} \frac{dn_\nu}{dE_\nu}(E_\nu) \int\limits_\Omega \frac{A(r')\rho(r')P_{ee}(E_\nu,|r-r'|)}{4\pi|r-r'|^2} d^3r.$$
Similarly to the Solar Neutrino Unit (SNU), the Terrestrial Neutrino Unit (TNU) is defined with reference to a specific detection mechanism (IBD, i.e., $\bar\nu_e+p\rightarrow n+e^+$): 1 TNU means to 1 fully efficient detection per 1 year per 1e32 free protons.
To get from flux to TNU signal, one needs account for the probability of IBD interaction (expressed as the interaction cross section $\sigma_\text{IBD}$),
$$S_\text{TNU}(r) = N_p T_\text{exp} F_\text{eff} \int \frac{d\phi}{dE_\nu}(E_\nu,r) ~\sigma_\text{IBD}(E_\nu) ~dE_\nu,$$
where $N_p=10^{32}$ (protons), $T_\text{exp}=1$ year, $F_\text{eff}=1$ (100% efficiency).
Now, if one accounts for the full enegy- and distance-dependent neutrino oscillation, the flux spectrum $d\phi/dE_\nu$ will have different shape than the neutrino production spectrum $dn_\nu/dE_\nu$, and will be position dependent. In this case, one needs to use the above equation to calculate TNU signal $S_\text{TNU}$.
If, however, one assumes an average survival probability (which is simply a numerical factor between 0 and 1), geoneutrino flux at each position has same shape as the antineutrino production spectrum and one can write
$$\frac{d\phi}{dE_\nu}(E_\nu,r) = \phi(r) \times \frac{1}{n_\nu} \frac{dn_\nu}{dE_\nu}(E_\nu),$$
where $\phi=\int(d\phi/dE_\nu)dE_\nu$ is the geoneutrino flux (in $cm^{-2}~s^{-1}$), and deduce the TNU to flux ratio
$$\frac{S_\text{TNU}(r)}{\phi(r)} = N_p T_\text{exp} F_\text{eff} \int \frac{1}{n_\nu}\frac{dn_\nu}{dE_\nu}(E_\nu) ~\sigma_\text{IBD}(E_\nu) ~dE_\nu.$$
I performed the last equation's calculation for both 232Th and 238U using these inputs:
- Enomoto's antineutrino productions spectra $dn_\nu/dE_\nu$
- IBD cross section $\sigma_\text{IBD}$ formula from Dye 2012 equation 11
Noting that with my trivial trapezoidal integration of Enomoto's spectra with 1 keV discretization, they do not yield exactly 4 resp 6 antineutrinos, but close enough and I do not worry about it.
I do recover Dye 2012 results for $\phi/S_\text{TNU}$ in units of $cm^{-2}~s^{-1}~TNU^{-1}$. Please see below.
import numpy as np
import matplotlib.pyplot as plt
from math import sqrt
secinyr = 365.25 * 24 * 3600
me = 0.5109989461
Mp = 938.2720813
Mn = 939.5654133
Delta = Mn - Mp
print ('Mn - Mp (MeV) =', Delta)
print ('IBD threshold (MeV) =', Delta+me)
def cross_section_IBD(E):
'''simply Dye 2018 RG, equation 11'''
cs = np.zeros(E.size)
for i in np.arange(E.size):
if E[i] > Delta+me:
cs[i] = 9.52e-44 * (E[i]-Delta)**2 * sqrt(1-(me/(E[i]-Delta))**2)
return cs
232Th¶
# read-in Enomoto's spectrum
E, dndE = np.loadtxt('AntineutrinoSpectrum_232Th.knt', usecols=(0, 1), unpack=True)
# keV to MeV
fac = 0.001
E = E * fac
dndE = dndE / fac
dE = E[1]-E[0]
print('dE =', dE)
nnu = np.sum(dndE) * dE
print('n_nu =', nnu)
sigma = cross_section_IBD(E)
TNUperFLUX = 1e32 * secinyr /nnu * np.sum(dndE*sigma)*dE
FLUXperTNU = 1 / TNUperFLUX
print('\nFLUXperTNU (cm-2 s-1 TNU-1) = {:.3e}'.format(FLUXperTNU))
print('\n[Compare to Dye 2012: 2.5e5]')
dndE_Th = dndE
238U¶
# read-in Enomoto's spectrum
E, dndE = np.loadtxt('AntineutrinoSpectrum_238U.knt', usecols=(0, 1), unpack=True)
# keV to MeV
fac = 0.001
E = E * fac
dndE = dndE / fac
dE = E[1]-E[0]
print('dE =', dE)
nnu = np.sum(dndE) * dE
print('n_nu =', nnu)
sigma = cross_section_IBD(E)
TNUperFLUX = 1e32 * secinyr /nnu * np.sum(dndE*sigma)*dE
FLUXperTNU = 1 / TNUperFLUX
print('\nFLUXperTNU (cm-2 s-1 TNU-1) = {:.3e}'.format(FLUXperTNU))
print('\n[Compare to Dye 2012: 7.6e4]')
dndE_U = dndE
plt.rcParams['figure.figsize'] = (8,6)
# plot of 232Th antineutrino production spectrum
plt.semilogy(E, dndE_Th, label="232Th")
plt.semilogy(E, dndE_U, label="238U")
plt.title("Antineutrino production spectrum (Enomoto)")
plt.xlabel("Neutrino energy (MeV)")
plt.ylabel("dn/dE (1/MeV)")
plt.legend()
plt.show()
# plot of IBD cross section
plt.semilogy(E, sigma)
plt.title("IBD cross section (following Dye 2012)")
plt.xlabel("Neutrino energy (MeV)")
plt.ylabel("Cross section (cm2)")
plt.show()
