Worksheet Geoneutrinos (for Encyclopedia of Geology, 2nd edition)¶
Ondřej Šrámek, Charles University, ondrej.sramek@gmail.com, 22 July 2020
References
- PDG 2020 = P.A. Zyla et al. (Particle Data Group), to be published in Prog. Theor. Exp. Phys. 2020, 083C01 (2020) http://pdg.lbl.gov
- Vogel & Beacom (1999): "Angular distribution of neutron inverse beta decay, νe + p → e+ + n." Phys. Rev. D 60, 053003, doi:10.1103/PhysRevD.60.053003
In [1]:
import numpy as np
import matplotlib.pyplot as plt
#print (plt.rcParams['figure.figsize'])
plt.rcParams['font.size'] = 12
In [2]:
me_MeV = 0.5109989461 # electron mass [PDG 2020]
mp_MeV = 938.272081 # proton mass [PDG 2020]
mn_MeV = 939.565413 # neutron mass [PDG 2020]
md_MeV = 1875.61294257 # deuteron mass [2018 CODATA https://physics.nist.gov/cgi-bin/cuu/Value?mdc2mev]
#nu_e mass < 1.1 eV [PDG 2020]
secinyr = 31556925.445 # Holden et al. 2011 https://doi.org/10.1351/PAC-REC-09-01-22
ibd_thresh_lab = ((mn_MeV + me_MeV)**2 - mp_MeV**2) / (2*mp_MeV) # Vogel & Beacom (1999) eqn. A2
print ('IBD threshold lab =', ibd_thresh_lab, 'MeV')
Eprompt = mn_MeV - me_MeV - mp_MeV
print ('E prompt = E_nu +', Eprompt, 'MeV')
Edelayed = mp_MeV + mn_MeV - md_MeV
print ('E delayed =', Edelayed, 'MeV')
print ('mn - mp =', mn_MeV - mp_MeV, 'MeV')
print ('secinyr * 1e32 =', secinyr * 1e32)
In [3]:
# plot pie chart
labels1 = '$^{238}$U', '$^{232}$Th', '$^{235}$U (0.4%)', '$^{40}$K antineutrinos', '$^{40}$K neutrinos', '$^{87}$Rb', '$^{187}$Re (1.0%)'
geonu = [12.20, 10.04, 0.37, 64.48, 7.61, 4.29, 0.98]
explode1 = (0.08, 0.08, 0.08, 0.08, 0.08, 0.08, 0.08)
labels2 = '$^{238}$U', '$^{232}$Th', '$^{235}$U', '$^{40}$K', '$^{147}$Sm + $^{87}$Rb'
heat = [38.36, 40.17, 1.65, 19.30, 0.52]
explode2 = (0.08, 0.08, 0.08, 0.08, 0.08)
def my_autopct(pct):
return ('%.1f%%' % pct) if pct > 2 else ''
fig1, (ax1, ax2) = plt.subplots(1,2,figsize=(11,4))
ax1.pie(geonu, explode=explode1, labels=labels1, autopct=my_autopct, shadow=True, startangle=340)
ax1.axis('equal')
ax1.set_title('Geoneutrino luminosity')
ax2.pie(heat, explode=explode2, labels=labels2, autopct=my_autopct, shadow=True, startangle=220)
ax2.axis('equal') # Equal aspect ratio ensures that pie is drawn as a circle.
ax2.set_title('Heat production')
ax2.text(-1.45, 0.49, '(1.7%)', fontsize=12)
ax2.text(-1.5, -0.93, '(0.5%)', fontsize=12)
#plt.savefig('../manuscript_latex/fig1_pies.pdf', bbox_inches='tight')
plt.show()
In [4]:
# spectra from http://www.awa.tohoku.ac.jp/~sanshiro/research/geoneutrino/spectrum/
class Nuclide:
def __init__(self, name):
self.name = name
# units of input spectra files: keV dn/keV
k40 = Nuclide("40K")
k40.x, k40.y = np.loadtxt('AntineutrinoSpectrum_40K.knt', usecols=(0, 1), unpack=True)
rb87 = Nuclide("87Rb")
rb87.x, rb87.y = np.loadtxt('AntineutrinoSpectrum_87Rb.knt', usecols=(0, 1), unpack=True)
th232 = Nuclide("232Th")
th232.x, th232.y = np.loadtxt('AntineutrinoSpectrum_232Th.knt', usecols=(0, 1), unpack=True)
u235 = Nuclide("235U")
u235.x, u235.y = np.loadtxt('AntineutrinoSpectrum_235U.knt', usecols=(0, 1), unpack=True)
u238 = Nuclide("238U")
u238.x, u238.y = np.loadtxt('AntineutrinoSpectrum_238U.knt', usecols=(0, 1), unpack=True)
facx = 1e-3
facy = 1e3
for N in u238, th232, u235, k40, rb87:
N.x = N.x * facx
N.y = N.y * facy
In [5]:
# plot geoneutrino spectra figure
fig, ax = plt.subplots(1,1,figsize=(8,6))
for N in u238, th232, u235, k40, rb87:
ax.semilogy(N.x, N.y, label=N.name)
ax.set_xlim(0,3.5)
ax.set_ylim(1e-2,6e1)
ax.set_xlabel('E (MeV)')
ax.set_ylabel('dn/dE (1/MeV/decay)')
ax.legend()
plt.show()
In [6]:
# number of antineutrinos above IBD threshold per decay
print(th232.x[1804:])
#print(u238.x[1804:])
dx = 0.001
for N in th232, u238,:
N.ntot = np.sum(N.y) * dx
N.nIBD = np.sum(N.y[1804:-1]) * dx
N.fIBD = N.nIBD / N.ntot
print ("{:5s} ntot = {:.4f} nIBD = {:.4f} fIBD = {:.5f}".format(N.name, N.ntot, N.nIBD, N.fIBD))
In [7]:
x_MeV = th232.x
sigmaIBD = np.zeros(x_MeV.size)
xtmp = x_MeV[1804:]
Delta = mn_MeV - mp_MeV
sigmaIBD[1804:] = 9.52e-44 * (xtmp-Delta)**2 * np.sqrt(1.-(me_MeV/(xtmp-Delta))**2)
fig, ax = plt.subplots(1,1,figsize=(8,6))
plt.semilogy(x_MeV, sigmaIBD)
plt.xlim(0,3.5)
plt.ylim(1e-46,1e-42)
plt.xlabel('E (MeV)')
plt.ylabel('$\sigma_{IBD}$ (cm$^2$)')
plt.show()
In [8]:
fig, ax = plt.subplots(1,1,figsize=(5,6))
for N in u238, th232:
ax.semilogy(N.x, N.y * sigmaIBD, label=N.name)
plt.xlim(0,3.5)
plt.ylim(1e-46,1e-43)
plt.xlabel('E (MeV)')
plt.ylabel('Interaction spectrum (cm$^2$/MeV/decay)')
plt.legend()
plt.show()
In [9]:
fig1, (ax1, ax2) = plt.subplots(1, 2, figsize=(14,6), gridspec_kw={'width_ratios': [1.8, 1], 'wspace': 0.4})
#
for N in u238, th232, u235, k40, rb87:
ax1.semilogy(N.x, N.y, label=N.name)
ax1.set_xlim(0,3.5)
ax1.set_ylim(1e-2,6e1)
ax1.set_title('Geoneutrino production spectrum')
ax1.set_xlabel('E (MeV)')
ax1.set_ylabel('dn/dE (1/MeV/decay)')
ax1.legend(bbox_to_anchor=(0.3, 1))
ax1b = ax1.twinx()
ax1b.semilogy(x_MeV, sigmaIBD, 'k--', label='sigma')
ax1b.text(2.3, 1.3e-43, 'IBD cross section', fontsize=12, rotation=30)
ax1b.set_ylim(1e-45,1e-42)
ax1b.set_ylabel('$\sigma_{IBD}$ (cm$^2$)')
#
for N in u238, th232:
ax2.semilogy(N.x, N.y * sigmaIBD, label=N.name)
ax2.set_xlim(0,3.5)
ax2.set_ylim(1e-46,1e-43)
ax2.set_title('IBD interaction spectrum')
ax2.set_xlabel('E (MeV)')
ax2.set_ylabel('Interaction spectrum (cm$^2$/MeV/decay)')
ax2.legend()
#
#plt.savefig('../manuscript_latex/fig2_spectra.pdf', bbox_inches='tight')
plt.show()
In [10]:
fig1, (ax1, ax2) = plt.subplots(1,2,figsize=(8,5))
N = 3
ind = np.arange(N) # the x locations for the groups
width = 0.5 # the width of the bars: can also be len(x) sequence
y = 32.1
p = 5.0
m = 5.0
means = (23.6, 27.9, 26.1)
plus = (6.8, 7.3, 2.6)
minus = (5.2, 5.8, 2.6)
ax1.plot((-1,N), (y,y), 'C4')
ax1.plot((-1,N), (y+p,y+p), '--C4')
ax1.plot((-1,N), (y-m,y-m), '--C4')
ax1.text(-0.3, 33, 'measurement', color='C4')
ax1.set_xlim(-0.5,2.5)
ax1.set_ylim(0,58)
ax1.set_ylabel('Geoneutrino signal in TNU')
ax1.bar(ind, means, width, yerr=(minus,plus), color=('C0','C2','C1'))
ax1.set_title('KamLAND')
ax1.set_xticks(ind)
ax1.set_xticklabels(('Huang+13', 'Wipp+20', 'KamLAND19'))
y = 47.0
p = 8.6
m = 8.1
means = (31.9, 32.3, 25.9)
plus = (7.3, 7.9, 4.9)
minus = (5.8, 6.4, 4.1)
ax2.plot((-1,N), (y,y), 'C4')
ax2.plot((-1,N), (y+p,y+p), '--C4')
ax2.plot((-1,N), (y-m,y-m), '--C4')
ax2.text(-0.3, 48, 'measurement', color='C4')
ax2.set_xlim(-0.5,2.5)
ax2.set_ylim(0,58)
ax2.bar(ind, means, width, yerr=(minus,plus), color=('C0','C2','C3'))
ax2.set_title('Borexino')
plt.xticks(ind, ('Huang+13', 'Wipp+20', 'Borexino20'))
#plt.savefig('../manuscript_latex/fig7_discrepancy.pdf', bbox_inches='tight')
plt.show()
In [11]:
N = 7
label = ('OBD', 'ANDES', 'KamLAND', 'Borexino', 'JNE', 'SNO+', 'JUNO')
# geonu from Uncle O [CRUST1.0] as in Šrámek et al. 2016 Sci.Rep.
# total, mantle, far-field litho, near-field litho
gtot = np.array([ 11.7, 41.8, 34.8, 41.4, 58.5, 44.2, 38.9 ])
gman = np.array([ 8.7, 8.3, 8.3, 8.2, 8.1, 8.2, 8.2 ])
gffl = np.array([ 2.4, 10.7, 9.4, 15.4, 20.6, 16.5, 13.5 ])
gnlf = gtot - gman - gffl
# reactor in geonu energy range
# Baldoncini et al. 2015; Wan et al. arXiv:1612.00133)
reac = np.array([ 0.9, 0.9, 18.3, 22.2, 6.8, 47.8, 354.5 ])
# total antineutrino
nutot = gtot + reac
# % total signal
relgman = gman / nutot * 100
relgffl = gffl / nutot * 100
relgnlf = gnlf / nutot * 100
relreac = reac / nutot * 100
print ('NFL(TNU) =', gnlf)
print ('Tot(TNU) =', nutot)
print ('Man/gTot =', gman/gtot)
print ('FFL/gTot =', gffl/gtot)
print ('NFL/gTot =', gnlf/gtot)
print ('Man/Tot =', relgman/100.)
ind = np.arange(N) + 0.2 # the x locations for the groups
width = 0.6 # the width of the bars: can also be len(x) sequence
# Two subplots, the axes array is 1-d
f, (ax0, ax1) = plt.subplots(2, sharex=True, figsize=(8, 6))
#ax0.set_title('Antineutrino detectors')
ax0.set_xticks(ind)
ax0.set_xticklabels(label)
# top subplot
ax0.bar(ind, relgman, width, edgecolor='black', color='#8bc42d')
ax0.bar(ind, relgffl, width, edgecolor='black', color='#206eb0', bottom=relgman)
ax0.bar(ind, relgnlf, width, edgecolor='black', color='#e85d1e', bottom=relgman+relgffl)
ax0.bar(ind, relreac, width, edgecolor='black', color='#aaaaaa', bottom=relgman+relgffl+relgnlf)
ax0.set_ylim((0,108))
ax0.set_yticks(np.arange(0, 104, 10))
ax0.set_ylabel('% antineutrino signal')
# annotate mantle/antinu percent
for xp, yp, zp in zip(ind, relgman/2., relgman):
if zp > 8.:
ax0.text(xp, yp, "{:.0f}%".format(zp), horizontalalignment='center', verticalalignment='center')
# bottom subplot
ax1.bar(ind, gman, width, edgecolor='black', color='#8bc42d')
ax1.bar(ind, gffl, width, edgecolor='black', color='#206eb0', bottom=gman)
ax1.bar(ind, gnlf, width, edgecolor='black', color='#e85d1e', bottom=gman+gffl)
ax1.bar(ind, reac, width, edgecolor='black', color='#aaaaaa', bottom=gman+gffl+gnlf)
ax1.set_ylim((0,100))
ax1.set_yticks(np.arange(0, 101, 10))
ax1.set_ylabel('Antineutrino signal in TNU')
# legend
ax1.scatter(0, 90, s=200, color='#aaaaaa', marker='s', edgecolors='black')
ax1.scatter(0, 78, s=200, color='#e85d1e', marker='s', edgecolors='black')
ax1.scatter(0, 66, s=200, color='#206eb0', marker='s', edgecolors='black')
ax1.scatter(0, 54, s=200, color='#8bc42d', marker='s', edgecolors='black')
ax1.text(0.2, 90, 'Reactor background', verticalalignment='center')
ax1.text(0.2, 78, 'Crust: closest 500 km', verticalalignment='center')
ax1.text(0.2, 66, 'Crust: rest of the world', verticalalignment='center')
ax1.text(0.2, 54, 'Mantle', verticalalignment='center')
ap = dict(width=0.2, facecolor='black', headwidth=8, headlength=8 )
xp = ind[N-1]-0.15
ax1.annotate("{:.0f} TNU".format(nutot[N-1]), xy=(xp-0.06, 98), xytext=(xp-0.06, 78), arrowprops=ap)
# save figure
#plt.savefig('../manuscript_latex/fig8_7detectors.pdf', bbox_inches='tight')
plt.show()
In [ ]: