Long-lived and short-lived radionuclides of the Earth¶
Ondřej Šrámek, Charles University, ondrej.sramek@gmail.com, 18 March 2020
Related publication: W F McDonough, O Šrámek, and S A Wipperfurth: Radiogenic power and geoneutrino luminosity of the Earth and other terrestrial bodies through time. Geochemistry, Geophysics, Geosystems, 2020.
Goal: Evaluate radiogenic heating and neutrino luminosity in radioactive decays of natural radionuclides.
Method: We systematically investigate radioactive decays using available data (NuDat + references, other sources), quantify the energies emitted in α,β,γ radiations and carried by neutrinos. Using estimates of Earth's composition and geochemically inferred ratios of radionuclides on early Earth, we quantify the radiogenic power and geoneutrino luminosity of Earth through time.
Units in calculations: mass and energy in keV
How to run: With Jupyter Notebook installed (e.g., via Anaconda), run "jupyter-notebook radionuclidesEarth.ipynb" in Terminal/Linux shell.
References¶
- NIST: Atomic Weights and Isotopic Compositions for All Elements)
- The NIST Reference on Constants, Units, and Uncertainty (physics.nist.gov/cuu/)
- NNDC Chart of Nuclides (www.nndc.bnl.gov/chart)
- NuDat 2.7 (www.nndc.bnl.gov/nudat2)
- Enomoto: Geoneutrino Spectrum and Luminosity (www.awa.tohoku.ac.jp/~sanshiro/research/geoneutrino/spectrum/)
- Bielajew: Introduction to Special Relativity, Quantum Mechanics and Nuclear Physics for Nuclear Engineers, Chapter 15, PDF
- Chambat et al. (2010): "Flattening of the Earth: further from hydrostaticity than previously estimated" Geophys. J. Int., 183(2), 727-732, doi:10.1111/j.1365-246X.2010.04771.x
- Holden et al. (2011): "IUPAC-IUGS common definition and convention on the use of the year as a derived unit of time (IUPAC Recommendations 2011)" Pure Appl. Chem., 83(5), 1159, doi:10.1351/PAC-REC-09-01-22
- Mougeot (2015): "Reliability of usual assumptions in the calculation of β and ν spectra" Phys. Rev. C, 91(5), 055504, doi:10.1103/PhysRevC.91.055504
- Renne et al. (2011): "Response to the comment by W. H. Schwarz et al. ...." Geochim. Cosmochim. Acta, 75(17), 5097–5100, doi:10.1016/j.gca.2011.06.021
- Yoder (1995): "Astrometric and geodetic properties of earth and the solar system" In Global Earth Physics: A Handbook of Physical Constants, AGU Reference Shelf, vol. 1, edited by Ahrens, pp. 1–31, doi:10.1029/RF001p0001
Additional sources referred to throughout the notebook.
Background¶
β− decay
$$_Z^AP ~~\rightarrow~~ _{Z+1}^{~~~~A}D + e^- + \bar\nu_e$$$$Q = m_P - m_D$$The highest possible value of maximum electron kinetic energy $T_e^{max}$ is $Q$ in the case of β− transition to the ground state of daughter nuclide. $T_e^{max}$ can be smaller when the parent transitions to an excited state of the daughter nucleus. The remaining deexcitation energy $E_\gamma$ is released as electromagnetic radiation. The energy $\langle E_\nu\rangle$ carried away by antineutrino, on average, is the difference between $T_e^{max}$ and the mean electron kinetic energy $\langle T_e\rangle$, that is, $T_e^{max}-\langle T_e\rangle$. The energy available as radiogenic heat is then
$$h = \langle T_e\rangle + E_\gamma = Q - \langle E_\nu\rangle = Q - (T_e^{max}-\langle T_e\rangle),$$which for a transition to ground state (i.e., $E_\gamma=0$ and $T_e^{max}=Q$) is simply equal to $\langle T_e\rangle$.
Electron capture
$$_Z^AP + e^- ~~\rightarrow~~ _{Z-1}^{~~~~A}D + \nu_e$$$$Q = m_P - m_D$$$$h = E_\gamma$$The transition energy available as heat only comes from deexcitation of daughter nucleus to ground state ($E_\gamma$). If parent decays directly to ground state of daughter nucleus, all energy is carried away by neutrino.
β+ decay
$$_Z^AP ~~\rightarrow~~ _{Z-1}^{~~~~A}D + e^+ + \nu_e$$$$Q = m_P - m_D - 2m_e$$$$h = \langle T_e\rangle + E_\gamma + 2m_e$$The energy available as radiogenic heat is the sum of mean positron kinetic energy, deexcitation of daughter nucleus (unless transition to ground state), and the energy released upon anihilation of the emitted positron with an ambient electron.
Some basic relationships
For a massive particle x, the total energy $E_x$ is the sum of rest-mass energy $m_xc^2$ and the kinetic energy $T_x$:
$$E_x = m_xc^2 + T_x$$The relativistic relation between its energy $E_x$ and momentum $p_x$ is:
$$E_x^2 = p_x^2c^2+m_x^2c^4$$Fermi's 1934 quantitative theory of β-decay and beyond
Assuming that
- the $Q$ value of the β transition is much smaller than the rest mass energies of the parent and daughter nuclei,
- neutrinos are mass-less,
the total energy $E_0$ available in the COM system is the sum of energy of the electron $E_e=m_ec^2+T_e$ and of the (anti)neutrino $E_\nu=T_\nu$:
$$E_0 = E_e + E_\nu = m_ec^2 + T_e + T_\nu = m_ec^2 + Q$$Working in units $\hbar=m_e=c=1$, we can write quantities as follows:
- $w_0=1+Q$ ... total energy available in the system (i.e., $E_0 = m_ec^2 + Q$)
- $w=1+E$ ... total energy of β-particle (i.e., $E_e = m_ec^2 + T_e$)
- $p=\sqrt{w^2-1}$ ... momentum of β-particle (i.e., $p_xc = \sqrt{E_x^2-m_x^2c^4}$)
- $q$ ... total energy of neutrino = momentum of neutrino (neutrino mass neglected)
The shape of a β spectrum is proportional to:
- phase space factor, $pwq^2$,
- Fermi function, $F(Z,w)$, which accounts for the Coulombic interaction between the nucleus and outgoing β-particle,
- shape factor, $S(w)$, which accounts for the "forbidden" nature of some transitions.
is the electron spectrum. That is, after normalizing to the number of β particles emitted per transition, it is the probability of a β particle to be created with energy in the $dw$ vicinity of $w$, where $w$ goes from $1$ to $w_0=1+Q$ (i.e., includes the rest-mass energy).
It is also the probability of a ν particle to be created with the energy $q=w_0-w$ (i.e., symmetry of electron- and ν-spectra ... except for radiative corrections).
$S(w)$ is often written in the form $S(p,q)$ (for a given $w$, one easily calculates $p$ and $q$ using the relations above).
The theoretical shape factors for unique forbidden transitions are (Bielajew):
Python imports and function definitions¶
In [1]:
#%reset -f
In [2]:
import numpy as np
import scipy.special
from math import sqrt, pi, exp, log
import matplotlib.pyplot as plt
In [3]:
plt.rcParams['figure.figsize'] = (8,5)
In [4]:
def Emesh(Emax, dE):
'''equidistant discretization of energy from 0 to Emax at midpoints of intervals'''
E = np.arange(0.5*dE, Emax, dE)
return E;
In [5]:
def fermi_func(w, Z):
'''Fermi function, after Enomoto's PhD thesis, eqn. 2.7
Neglecting the scalar factors, as spectrum will be normalized to branching'''
w1 = w/me
gamma = sqrt(1-(alpha*Z)**2)
y = alpha * Z * w1/np.sqrt(w1**2-1)
F = (np.sqrt(w1**2-1))**(2*(gamma-1)) * np.exp(pi*y) * abs( scipy.special.gamma( gamma + y*1j ) )**2
return F;
In [6]:
def int_y(dydx, dx):
'''calculates area under curve, simple trapezoidal rule'''
yint = np.sum(dydx)*dx
return yint;
def mean_x(dydx, x, dx):
'''calculates mean x, simple trapezoidal rule'''
xmean = np.sum(dydx*x)*dx / int_y(dydx, dx)
return xmean;
def normalize_dydx(dydx, dx, norm):
'''normalizes spectrum to norm'''
ynorm = dydx / int_y(dydx, dx) * norm
return ynorm;
In [7]:
# decay-independent inputs
me = 510.99895000 # keV mass of electron [2018 CODATA via http://physics.nist.gov/cgi-bin/cuu/Value?mec2mev]
alpha = 7.2973525693e-3 # fine-structure constant [2018 CODATA via https://physics.nist.gov/cgi-bin/cuu/Value?alph]
MeV_per_u = 931.49410242 # [2018 CODATA via http://physics.nist.gov/cgi-bin/cuu/Value?muc2mev]
keV_per_u = MeV_per_u * 1000.
pJ_in_keV = 1.602176634e-4 # [2018 CODATA via https://physics.nist.gov/cgi-bin/cuu/Value?evj]
J_in_keV = pJ_in_keV*1e-12
mass_4He_u = 4.00260325413 # u mass of 4He [https://www.nist.gov/pml/atomic-weights-and-isotopic-compositions-relative-atomic-masses]
NA = 6.02214076e23 # Avogadro's number [2018 CODATA via https://physics.nist.gov/cgi-bin/cuu/Value?na]
secinyr = 3.1556925445e7 # IUPAC Recommendations 2011 (Holden et al 2011 https://doi.org/10.1351/PAC-REC-09-01-22)
Mearth = 5.97218e24 # mass of Earth in kg (Chambat et al 2010 https://doi.org/10.1111/j.1365-246X.2010.04771.x)
Mcore = 1.93e24 # mass of core in kg (Yoder 1995 https://doi.org/10.1029/RF001p0001)
Mbse = Mearth-Mcore # mass of BSE in kg
print("Mass of BSE in kg = {:.4e} (= Mearth - Mcore, and references just above)".format(Mbse))
dE = 1.0 # keV energy discretization
tcai = 0.
tnow = 4.568e9
In [8]:
class Nuclide:
"""
Possible properties:
halflife in years
massfrac_el
mass_atomic atomic mass of parent element in u (=g/mol)
mass_nuclide nuclide mass of paren nuclide in u (=g/mol)
isotfrac natural isotopic mole fraction of parent nuclide (#atoms_nuclide/#atoms_element)
timeref reference time of input abundace (t_CAI for short-lived, t_now for long-lived)
branch_beta branching ratio of beta- decay
branch_ec branching ratio of EC decay
branch_betaplus branching ratio of beta+ decay
Qalpha Q value of alpha decay
Qbeta Q value of beta- decay
Qec Q value of EC decay
Qbetaplus Q value of beta+ decay
Qmean mean Q value for branched decays (weighted by branching ratios)
anu_mean mean antineutrino energy in beta- branch
anu_away energy carried away by antineutrino per decay (accounting for branching)
el_mean mean electron energy in beta- branch
nu_away energy carried away by neutrino per decay (accounting for branching)
pos_mean mean positron energy in beta+ branch
heat_beta heat per decay in beta- branch (accounting for branching)
heat_ec heat per decay in EC branch (accounting for branching)
heat_betaplus heat per decay in beta+ branch (accounting for branching)
heat_gamma heat per decay due to deexcitation energy (electromagnetic radiation)
heat heat per decay
lanu # of antineutrinos emitted per decay
lnu # of neutrinos emitted per decay
facH just an internal scalar value equal to Joules per keV
Methods:
atoms_per_kgrock()
decay_rate_per_kgrock()
Heat_per_kgrock()
Lanu_per_kgrock()
Lnu_per_kgrock()
fac_rate()
facHeat()
facLanu()
facLnu()
heatJoule()
activity_becq()
"""
def __init__(self, name):
self.name = name
def atoms_per_kgrock(self):
return self.massfrac_el / self.mass_atomic * 1000 * NA * self.isotfrac
def decay_rate_per_kgrock(self):
return self.atoms_per_kgrock() * log(2.)/self.halflife/secinyr
facH = pJ_in_keV*1e-12 # Joules per keV
def Heat_per_kgrock(self):
return self.decay_rate_per_kgrock() * self.heat * self.facH
def Lanu_per_kgrock(self):
return self.decay_rate_per_kgrock() * self.lanu
def Lnu_per_kgrock(self):
return self.decay_rate_per_kgrock() * self.lnu
def fac_rate(self):
return log(2.)/self.halflife/secinyr / self.mass_atomic * 1000 * NA * self.isotfrac
def facHeat(self):
return self.fac_rate() * self.heat * self.facH
def facLanu(self):
return self.fac_rate() * self.lanu
def facLnu(self):
return self.fac_rate() * self.lnu
def heatJoule(self):
return self.heat * J_in_keV
def heat_gammaJoule(self):
return self.heat_gamma * J_in_keV
def activity_becq(self, time):
'''returns activity in becquerels at time where time can be an array'''
decc = log(2.)/(self.halflife)
return decc/secinyr * self.atoms_per_kgrock() * np.exp(decc*(self.timeref-time))
Long-lived radionuclides¶
Example: decay of 40K (going slow here...)¶
Branched β− or EC or β+ decay, half-life ~1.25 Gyr, differences in branching and half-life between NuDat (1.248(3)e9 yr) and other sources: Renne et al 2011 (1.253e9 yr), Naumenko-Dèzes et al 2018 (1.261e9–1.264e9 yr)
β− decay of 40K¶
Third unique forbidden transition, experimental shape factor from Mougeot 2015 who refers to Leutz et al 1965 who refer to Bühring 1965 who then also refers to Bühring 1963b90086-5) and Bühring 1963a90290-6).
In [9]:
K40 = Nuclide("40K")
# decay-specific inputs
mass_parent = 39.963998166 # 40K mass in u [NIST]
mass_daughter = 39.962590863 # 40Ca mass in u [NIST]
Z = 20 # atomic number of daughter
K40.mass_nuclide = mass_parent
K40.mass_atomic = 39.0983
#K40.isotfrac = 0.000117 # [NIST]
K40.isotfrac = 0.00011668 # [Naumenko et al 2013 https://doi.org/10.1016/j.gca.2013.08.019]
K40.timeref = tnow
branch_nudat = 0.8928 # NuDat 2.7
lambda_beta = 4.9548e-10 # [Renne et al 2011 https://doi.org/10.1016/j.gca.2011.06.021]
lambda_ec = 5.757e-11 # [Renne et al 2011]
lambda_renne = lambda_beta + lambda_ec
branch_renne = lambda_beta/lambda_renne
lambda_naumenko = 0.5*(5.484e-10+5.498e-10) # [Naumenko-Dèzes et al 2018 https://doi.org/10.1016/j.gca.2017.09.041]
branch_naumenko = 0.5*(0.8925+0.8963)
print ("branching = {:.4f} [Naumenko-Dèzes et al 2018]".format(branch_naumenko))
print ("branching = {:.4f} [Renne et al 2001]".format(branch_renne))
print ("branching = {:.4f} [NuDat]\n".format(branch_nudat))
t12_nudat = 1.248e9 # yr [NuDat 2.7]
t12_renne = log(2.)/lambda_renne
t12_naumenko = log(2.)/lambda_naumenko
print ("half-life = {:.3e} yr [Naumenko-Dèzes et al 2018]".format(t12_naumenko))
print ("half-life = {:.3e} yr [Renne et al 2011]".format(t12_renne))
print ("half-life = {:.3e} yr [NuDat]\n".format(t12_nudat))
# Choosing Naumenko values
branch_beta = branch_naumenko
K40.halflife = t12_naumenko
K40.branch_beta = branch_beta
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
print ('Q = {:.2f} keV (NuDat: 1310.89 ± 0.06 keV)'.format(Q))
K40.Qbeta = Q
In [10]:
Eendpoint = Q
q = Emesh(Eendpoint,dE)
E = Eendpoint - q
w = me + E
p = np.sqrt(w**2-me**2)
dNdE1 = normalize_dydx(p*w*q**2, dE, branch_beta)
plt.plot(q, dNdE1, color="C0", label='neutrino')
plt.plot(E, dNdE1, color="C2", label='electron')
plt.text(400, 0.0001, 'Only phase space factor, p*w*q^2')
plt.xlabel('E (keV)')
plt.ylabel('dN/dE (1/keV/decay)')
plt.legend()
plt.show()
print ("Mean neutrino energy/total available energy = {:.4f}".format(mean_x(dNdE1, q, dE)/Eendpoint))
print ("Mean electron energy/total available energy = {:.4f}".format(mean_x(dNdE1, E, dE)/Eendpoint))
In [11]:
dNdE2 = normalize_dydx(dNdE1 * fermi_func(w,Z), dE, branch_beta)
plt.plot(q, dNdE1, label='w/out Fermi function')
plt.plot(q, dNdE2, label='w/ Fermi function')
plt.text(400, 0.0001, 'Effect of Fermi function')
plt.xlabel('E (keV)')
plt.ylabel('dN/dE (1/keV/decay)')
plt.legend()
plt.title("Neutrino energy")
plt.show()
print ("Mean neutrino energy/total available energy = {:.4f}".format(mean_x(dNdE2, q, dE)/Eendpoint))
print ("Mean electron energy/total available energy = {:.4f}".format(mean_x(dNdE2, E, dE)/Eendpoint))
The positively charged nucleus Coulombically attracts the escaping electron, therefore shifts the electron spectrum toward lower energies, therefore the antineutrino spectrum is shifted toward higher energies.
In [12]:
shape_the = q**6 + 7*q**4*p**2 + 7*q**2*p**4 + p**6 # symmetrical, from Bielajew
shape_exp = 1.05*q**6 + 6.3*q**4*p**2 + 6.25*q**2*p**4 + 0.95*p**6 # experimental from Leutz et al 1965
dNdE_exp = normalize_dydx(dNdE2 * shape_exp, dE, branch_beta)
dNdE_renne_exp = normalize_dydx(dNdE2 * shape_exp, dE, branch_renne)
dNdE_nudat_exp = normalize_dydx(dNdE2 * shape_exp, dE, branch_nudat)
dNdE_the = normalize_dydx(dNdE2 * shape_the, dE, branch_beta)
dNdE_renne_the = normalize_dydx(dNdE2 * shape_the, dE, branch_renne)
dNdE_nudat_the = normalize_dydx(dNdE2 * shape_the, dE, branch_nudat)
plt.plot(q, dNdE2, label='w/out shape factor')
plt.plot(q, dNdE_exp, label='w/ shape factor, experimental (Leutz et al 1965)')
plt.plot(q, dNdE_the, label='w/ shape factor, theoretical (Bielajew)')
plt.text(500, 0.0004, 'Effect of shape factor')
plt.xlabel('E (keV)')
plt.ylabel('dN/dE (1/keV/decay)')
plt.legend()
plt.title("Neutrino energy")
plt.show()
plt.semilogy(q, dNdE_exp, label='w/ experimental shape factor (Leutz et al 1965)')
plt.xlabel('E (keV)')
plt.ylabel('dN/dE (1/keV/decay)')
plt.legend()
plt.title("Neutrino energy")
plt.show()
In [13]:
print ("Mean antineutrino energy = {:.1f} keV using experimental shape factor".format(mean_x(dNdE_exp, q, dE)))
print ("Mean antineutrino energy = {:.1f} keV using theoretical shape factor\n".format(mean_x(dNdE_the, q, dE)))
print ("Accounting for branching, the β- antineutrino carries away on average, in keV:")
print (" {:.1f} (Naumenko, experimental)".format( mean_x(dNdE_exp, q, dE) * branch_beta))
print (" {:.1f} (Renne, experimental)".format( mean_x(dNdE_renne_exp, q, dE) * branch_renne))
print (" {:.1f} (NUDAT, experimental)".format( mean_x(dNdE_nudat_exp, q, dE) * branch_nudat))
print (" {:.1f} (Naumenko, theoretical)".format( mean_x(dNdE_the, q, dE) * branch_beta))
print (" {:.1f} (Renne, theoretical)".format( mean_x(dNdE_renne_the, q, dE) * branch_renne))
print (" {:.1f} (NUDAT, theoretical)\n".format( mean_x(dNdE_nudat_the, q, dE) * branch_nudat))
print ("Neutrinos per decay =", int_y(dNdE_exp, dE), '(check of normalization)\n')
K40.anu_mean = mean_x(dNdE_exp, q, dE)
K40.el_mean = Eendpoint - K40.anu_mean
K40_el_mean_renne = Eendpoint - mean_x(dNdE_renne_exp, q, dE)
K40.anu_away = K40.anu_mean * K40.branch_beta
K40.heat_beta = (Eendpoint - K40.anu_mean) * K40.branch_beta
print ("Mean β− energy = {:.1f} keV (Renne: {:.1f} keV, NuDat: 560.2 ± 1.5 keV)".format(K40.el_mean, K40_el_mean_renne))
print ("Accounting for branching, the β- antineutrino carries away on average {:.1f} keV (Naumenko & exp)".format(K40.anu_away))
print ("Accounting for branching, the β- heat per decay is {:.1f} keV (Naumenko & exp)".format(K40.heat_beta))
Electron capture on 40K¶
In [14]:
# decay-specific inputs
mass_parent = 39.963998166 # 40K mass in u [NIST]
mass_daughter = 39.9623831237 # 40Ar mass in u [NIST]
branch_nudat_ec_gs = 0.00046
branch_nudat_ec_ex = 0.1067
branch_nudat_betaplus = 1e-5
# assuming NNDC branching for the gs->gs EC transition, beta+
branch_ec_ex = 1 - branch_beta - branch_nudat_ec_gs - branch_nudat_betaplus
K40.branch_ec = branch_ec_ex + branch_nudat_ec_gs
K40.branch_betaplus = branch_nudat_betaplus
Eexc = 1460.8 # energy of excited state above GS [NuDat]
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
print ('Q = {:.2f} keV (NuDat: 1504.40 ± 0.06 keV)\n'.format(Q))
K40.Qec = Q
K40.nu_away = branch_nudat_ec_gs * Q + branch_ec_ex * (Q-Eexc)
K40.heat_ec = branch_ec_ex * Eexc
print ("Accounting for branching, the EC neutrino carries away on average {:.3f} keV".format(K40.nu_away))
print ("Accounting for branching, the EC heat per decay is {:.1f} keV".format(K40.heat_ec))
In [15]:
Q = (mass_parent - mass_daughter) * keV_per_u - 2*me # keV Q-value
K40.Qbetaplus = Q
print ('Q = {:.2f} keV (NuDat: 482.40 ± 0.06 keV)\n'.format(Q))
K40.pos_mean = 197.325
K40.heat_betaplus = K40.branch_betaplus * (K40.pos_mean + 2*me)
print ("Mean positron energy = {:.3f} keV (from NuDat)".format(K40.pos_mean))
print ("Accounting for branching, the β+ heat per decay is {:.3f} keV".format(K40.heat_betaplus))
40K decay summary¶
In [16]:
K40.heat = K40.heat_beta + K40.heat_ec
K40.lanu = K40.branch_beta
K40.lnu = 1-K40.branch_beta
K40.Qmean = K40.Qbeta * K40.branch_beta + K40.Qec * K40.branch_ec + K40.Qbetaplus * K40.branch_betaplus
K40.heat_gamma = Eexc * branch_ec_ex # rel.err. < 1e-4
print ("Overall Q-value = {:.1f} keV = {:.4f} pJ \n".format(K40.Qmean, K40.Qmean*pJ_in_keV))
print ("β- decay: mean electron energy = {:.1f} keV/decay".format(K40.el_mean))
print (" mean nuebar energy = {:.1f} keV/decay".format(K40.anu_mean))
print (" carried away by nuebar = {:.1f} keV/decay".format(K40.anu_away))
print (" heats the Earth = {:.1f} keV/decay\n".format(K40.heat_beta))
print ("EC: carried away by nue = {:.1f} keV/decay".format(K40.nu_away))
print (" heats the Earth = {:.1f} keV/decay\n".format(K40.heat_ec))
print ("β+ decay: heats the Earth = {:.3f} keV/decay\n".format(K40.heat_betaplus))
print ("Overall: heats the Earth = {:.1f} keV/decay = {:.4f} pJ/decay\n".format(K40.heat, K40.heat * pJ_in_keV))
print ("(using Naumenko-Dèzes et al 2018 branching and β- spectrum shape factor from Leutz et al 1965)")
87Rb → 87Sr decay¶
β- decay with half-life 48.1 ± 0.9 Gyr (NuDat)
Third forbidden non-unique transition, shape factor from Mougeot 2015 who refers to Grau Carles & Kossert (2007)
In [17]:
# decay-specific inputs
mass_parent = 86.9091805310 # u [NIST]
mass_daughter = 86.9088775 # u [NIST]
Z = 38 # atomic number of daughter
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
print ('Q = {:.1f} keV (NuDat: 282.2 ± 1.1 keV)\n'.format(Q))
# NuDat inputs
branch = 1.
Q = 282.2
el_mean = 81.67
print ("Using NuDat inputs:")
print (" Mean electron energy = {:.2f} keV (from NuDat)\n".format(el_mean))
# calculating spectrum
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
Eendpoint = Q
q = Emesh(Eendpoint,dE)
E = Eendpoint - q
w = me + E
p = np.sqrt(w**2-me**2)
print ("My calculation:")
shape = 1
dNdE = normalize_dydx(p*w*q**2 * fermi_func(w,Z) * shape, dE, branch)
anu_mean = mean_x(dNdE, q, dE)
el_mean = Eendpoint - anu_mean
heat_use = el_mean
print (" Mean electron energy w/ shape factor equal to 1 = {:.2f} keV".format(el_mean))
#
shape = 0.011*p**4 + 0.305*p**2*q**2 + q**4 # Mougeot 2015 <- Grau Carles & Kossert 2007
dNdE = normalize_dydx(p*w*q**2 * fermi_func(w,Z) * shape, dE, branch)
anu_mean = mean_x(dNdE, q, dE)
el_mean = Eendpoint - anu_mean
print (" Mean electron energy w/ shape from Mougeot = {:.2f} keV".format(el_mean))
anu_away = anu_mean * branch
heat = el_mean * branch
Rb87 = Nuclide("87Rb")
Rb87.halflife = 49.61e9 # Villa et al 2015 https://doi.org/10.1016/j.gca.2015.05.025
Rb87.mass_atomic = 85.4677
Rb87.mass_nuclide = mass_parent
Rb87.isotfrac = 0.2783
Rb87.timeref = tnow
Rb87.Qbeta = Q
Rb87.heat = heat_use
Rb87.lanu = 1
Rb87.lnu = 0
Rb87.heat_gamma = 0.
print ("\nHeats the Earth = {:.1f} keV/decay = {:.5f} pJ/decay".format(Rb87.heat, Rb87.heat * pJ_in_keV))
In [18]:
# decay-specific inputs
mass_parent = 137.9071149 # u [NIST]
mass_daughter = 137.905991 # u [NIST]
Z = 58 # atomic number of daughter
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
print ('Q = {:.1f} keV (NuDat: 1052 ± 4 keV)\n'.format(Q))
# NuDat inputs
branch = 0.345
Q = 1052
el_mean = 98
Eexc = La138_heat_gamma_beta = 788.7
print ("Using NuDat inputs:")
print (" Mean electron energy = {:.2f} keV (from NuDat)".format(el_mean))
La138 = Nuclide("La138")
La138.halflife = 1.03e11 # Sato & Hirose 1981; Tanimizu 2000
La138.mass_atomic = 138.90547
La138.mass_nuclide = mass_parent
La138.isotfrac = 0.0008881
La138.Qbeta = Q
La138.branch_beta = branch
La138.heat_beta = (el_mean + Eexc) * La138.branch_beta
La138_heat_gamma_beta = Eexc * La138.branch_beta
print (" Heats the Earth = {:.2f} keV/decay (from NuDat)".format(La138.heat_beta))
EC decay 138La → 138Ba¶
In [19]:
# decay-specific inputs
mass_parent = 137.9071149 # u [NIST]
mass_daughter = 137.90524700 # u [NIST]
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
print ('Q = {:.1f} keV (NuDat: 1742 ± 3 keV)\n'.format(Q))
# NuDat inputs
branch = 0.655
Q = 1742
Eexc = 1435.8
print ("Using NuDat inputs")
La138.Qec = Q
La138.branch_ec = branch
La138.heat_ec = Eexc * La138.branch_ec
La138_heat_gamma_ec = Eexc * La138.branch_ec
print (" Heats the Earth = {:.2f} keV/decay (from NuDat)".format(La138.heat_ec))
138La decay summary¶
In [20]:
La138.Qmean = La138.Qbeta * La138.branch_beta + La138.Qec * La138.branch_ec
La138.heat = La138.heat_beta + La138.heat_ec
La138.timeref = tnow
La138.lanu = La138.branch_beta
La138.lnu = 1 - La138.branch_beta
La138.heat_gamma = La138_heat_gamma_beta + La138_heat_gamma_ec
print ("Qmean = {:.1f}".format(La138.Qmean))
print ("Heats the Earth = {:.2f} keV/decay = {:.5f} pJ/decay (from NuDat)".format(La138.heat, La138.heat * pJ_in_keV))
In [21]:
mass_parent = 146.9149044 # u [NIST]
mass_daughter = 142.9098200 # u [NIST]
Q = (mass_parent - mass_daughter - mass_4He_u) * keV_per_u # keV Q-value
print ('Q = {:.2f} keV (NuDat: 2311 ± 1 keV)'.format(Q))
Sm147 = Nuclide("147Sm")
Sm147.halflife = 106.3e9 # Begemann et al 2001; Tavares & Terranova 2018
Sm147.mass_atomic = 150.3616
Sm147.mass_nuclide = mass_parent
Sm147.isotfrac = 0.14993
Sm147.timeref = tnow
Sm147.Qalpha = Q
Sm147.heat = Q
Sm147.lanu = 0
Sm147.lnu = 0
Sm147.heat_gamma = 0.
print ("\nHeats the Earth = {:.1f} keV/decay = {:.4f} pJ/decay".format(Sm147.heat, Sm147.heat * pJ_in_keV))
In [22]:
# decay-specific inputs
mass_parent = 175.9426897 # u [NIST]
mass_daughter = 175.9414076 # u [NIST]
Z = 72 # atomic number of daughter
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
print ('Q = {:.1f} keV (NuDat: 1190.2 ± 0.8 keV)\n'.format(Q))
# NuDat inputs
Q = 1190.2
el_mean = 181.8
Eexc = 596.95
print ("Using NuDat inputs:")
print (" Mean electron energy = {:.2f} keV (from NuDat)".format(el_mean))
Lu176 = Nuclide("176Lu")
Lu176.halflife = 3.713e10 # Söderlund et al 2004; Hult et al 2014
Lu176.mass_atomic = 174.9668
Lu176.mass_nuclide = mass_parent
Lu176.isotfrac = 0.02599
Lu176.timeref = tnow
Lu176.Qbeta = Q
Lu176.heat = el_mean + Eexc
Lu176.lanu = 1
Lu176.lnu = 0
Lu176.heat_gamma = Eexc
print ("\nHeats the Earth = {:.1f} keV/decay = {:.5f} pJ/decay".format(Lu176.heat, Lu176.heat * pJ_in_keV))
In [23]:
# decay-specific inputs
mass_parent = 186.9557501 # u [NIST]
mass_daughter = 186.9557474 # u [NIST]
Z = 76 # atomic number of daughter
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
print ('Q = {:.3f} keV (NuDat: 2.469 ± 0.004 keV)\n'.format(Q))
# NuDat inputs
Q = 2.469
el_mean = 0.618
Eexc = 0
print ("Using NuDat inputs:")
print (" Mean electron energy = {:.2f} keV (from NuDat)".format(el_mean))
Re187 = Nuclide("187Re")
Re187.halflife = 4.153e10 # Smoliar et al 1996; Selby et al 2007
Re187.mass_atomic = 186.207
Re187.mass_nuclide = mass_parent
Re187.isotfrac = 0.6260
Re187.timeref = tnow
Re187.Qbeta = Q
Re187.heat = el_mean + Eexc
Re187.lanu = 1
Re187.lnu = 0
Re187.heat_gamma = Eexc
print ("\nHeats the Earth = {:.1f} keV/decay = {:.5f} pJ/decay".format(Re187.heat, Re187.heat * pJ_in_keV))
In [24]:
mass_parent = 189.9599297 # u [NIST]
mass_daughter = 185.9538350 # u [NIST]
Q = (mass_parent - mass_daughter - mass_4He_u) * keV_per_u # keV Q-value
print ('Q = {:.2f} keV (NuDat: 3249 ± 6 keV)'.format(Q))
Pt190 = Nuclide("190Pt")
#Pt190.halflife = 6.5e11 # NNDC
Pt190.halflife = 4.899e11 # Cook et al 2004; Braun et al 2017
Pt190.mass_atomic = 195.084
Pt190.mass_nuclide = mass_parent
Pt190.isotfrac = 0.00012
Pt190.timeref = tnow
Pt190.Qalpha = Q
Pt190.heat = Q
Pt190.lanu = 0
Pt190.lnu = 0
Pt190.heat_gamma = 0.
print ("\nHeats the Earth = {:.1f} keV/decay = {:.4f} pJ/decay".format(Pt190.heat, Pt190.heat * pJ_in_keV))
In [25]:
mass_parent = 232.0380558 # 232Th u [NIST]
mass_daughter = 207.9766525 # 208Pb u [NIST]
nalpha = 6
nbeta = 4
Qtotal = (mass_parent - mass_daughter - nalpha * mass_4He_u)* keV_per_u # keV Q-value
# read-in antineutrino spectrum by Enomoto: https://www.awa.tohoku.ac.jp/~sanshiro/research/geoneutrino/spectrum/
x, y = np.loadtxt('AntineutrinoSpectrum_232Th.knt', usecols=(0, 1), unpack=True)
plt.semilogy(x, y)
plt.xlabel('E (keV)')
plt.ylabel('dn/dE (1/keV/decay)')
plt.show()
# integrate antineutrino spectrum
print ('Qtotal = {:.0f} keV = {:.4f} pJ\n'.format(Qtotal, Qtotal * pJ_in_keV))
dx = 1.
nsum = np.sum(y) * dx
Esum = np.sum(x*y) * dx
print ("Check integrated n_anu = {:.4f} (should be {:})\n".format(nsum, nbeta))
print ("Integrated E_anu = {:.1f} keV/decay = {:.4f} pJ/decay\n".format(Esum, Esum * pJ_in_keV))
heat = Qtotal - Esum
Th232 = Nuclide("232Th")
Th232.halflife = 14.1e9 # Farley 1960 https://doi.org/10.1139/p60-114
Th232.mass_atomic = 232.0381
Th232.isotfrac = 1.
Th232.timeref = tnow
Th232.heat = heat
Th232.lanu = nbeta
Th232.lnu = 0
Th232.heat_gamma = -0.000001
print ("Heats the Earth = {:.1f} keV/decay = {:.4f} pJ/decay".format(heat, heat * pJ_in_keV))
In [26]:
mass_parent = 235.0439301 # 235U u [NIST]
mass_daughter = 206.9758973 # 207Pb u [NIST]
nalpha = 7
nbeta = 4
Qtotal = (mass_parent - mass_daughter - nalpha * mass_4He_u)* keV_per_u # keV Q-value
# read-in antineutrino spectrum by Enomoto: https://www.awa.tohoku.ac.jp/~sanshiro/research/geoneutrino/spectrum/
x, y = np.loadtxt('AntineutrinoSpectrum_235U.knt', usecols=(0, 1), unpack=True)
dx = 1.
plt.semilogy(x, y)
plt.xlabel('E (keV)')
plt.ylabel('dn/dE (1/keV/decay)')
plt.show()
# integrate antineutrino spectrum
print ('Qtotal = {:.0f} keV = {:.4f} pJ\n'.format(Qtotal, Qtotal * pJ_in_keV))
dx = 1.
nsum = np.sum(y) * dx
Esum = np.sum(x*y) * dx
print ("Check integrated n_anu = {:.4f} (should be {:})\n".format(nsum, nbeta))
print ("Integrated E_anu = {:.1f} keV/decay = {:.4f} pJ/decay\n".format(Esum, Esum * pJ_in_keV))
heat = Qtotal - Esum
U235 = Nuclide("235U")
U235.halflife = 0.70348e9 # Villa et al 2016 https://doi.org/10.1016/j.gca.2015.10.011
U235.mass_nuclide = mass_parent
U235.mass_atomic = 238.02891
U235.isotfrac = 0.0072037
U235.timeref = tnow
U235.heat = heat
U235.lanu = nbeta
U235.lnu = 0
U235.heat_gamma = -0.000001
print ("Heats the Earth = {:.1f} keV/decay = {:.4f} pJ/decay".format(U235.heat, U235.heat * pJ_in_keV))
In [27]:
mass_parent = 238.0507884 # 235U u [NIST]
mass_daughter = 205.9744657 # 207Pb u [NIST]
nalpha = 8
nbeta = 6
Qtotal = (mass_parent - mass_daughter - nalpha * mass_4He_u)* keV_per_u # keV Q-value
# read-in antineutrino spectrum by Enomoto: https://www.awa.tohoku.ac.jp/~sanshiro/research/geoneutrino/spectrum/
x, y = np.loadtxt('AntineutrinoSpectrum_238U.knt', usecols=(0, 1), unpack=True)
dx = 1.
plt.semilogy(x, y)
plt.xlabel('E (keV)')
plt.ylabel('dn/dE (1/keV/decay)')
plt.show()
# integrate antineutrino spectrum
print ('Qtotal = {:.0f} keV = {:.4f} pJ\n'.format(Qtotal, Qtotal * pJ_in_keV))
dx = 1.
nsum = np.sum(y) * dx
Esum = np.sum(x*y) * dx
print ("Check integrated n_anu = {:.4f} (should be {:})\n".format(nsum, nbeta))
print ("Integrated E_anu = {:.1f} keV/decay = {:.4f} pJ/decay\n".format(Esum, Esum * pJ_in_keV))
heat = Qtotal - Esum
U238 = Nuclide("238U")
U238.halflife = 4.4683e9 # Villa et al 2016 https://doi.org/10.1016/j.gca.2015.10.011
U238.mass_nuclide = mass_parent
U238.mass_atomic = 238.02891
U238.isotfrac = 0.992796
U238.timeref = tnow
U238.heat = heat
U238.lanu = nbeta
U238.lnu = 0
U238.heat_gamma = -0.000001
print ("Heats the Earth = {:.1f} keV/decay = {:.4f} pJ/decay".format(U238.heat, U238.heat * pJ_in_keV))
Recalculation 235U, 238U isotopic fractions and U atomic mass¶
In [28]:
print("Original [NIST]:", U235.isotfrac, U238.isotfrac, U238.mass_atomic)
u238to235 = 137.786 # [Connelly et al 2017 https://doi.org/10.1016/j.gca.2016.10.044]
u234to238 = 5.4970e-5 # [Villa et al 2016 https://doi.org/10.1016/j.gca.2015.10.011]
U238.isotfrac = 1 / (1 + 1/u238to235 + u234to238)
U235.isotfrac = U238.isotfrac / u238to235
U234_isotfrac = U238.isotfrac * u234to238
U235.mass_atomic = U238.mass_atomic = U238.isotfrac * U238.mass_nuclide + U235.isotfrac * U235.mass_nuclide + U234_isotfrac * 234.0409523
print("\nNew [Villa16+]: ", U235.isotfrac, U238.isotfrac, U238.mass_atomic)
Short-lived radionuclides¶
In [29]:
mass_parent = 40.96227792 # u [NIST]
mass_daughter = 40.9618252579 # u [NIST]
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
print ('Q = {:.2f} keV (NuDat: 421.66 ± 0.14 keV)\n'.format(Q))
Ca41 = Nuclide("41Ca")
Ca41.halflife = 9.94e4
Ca41.mass_atomic = 40.078
Ca41.mass_nuclide = mass_parent
Ca41.isotfrac = 4.6e-9 * 0.96941 # 41Ca/40Ca * 40Ca/Ca
Ca41.timeref = tcai
Ca41.Qec = Q
Ca41.heat = 0.
Ca41.lanu = 0
Ca41.lnu = 1
Ca41.heat_gamma = 0.
print ("Heats the Earth = {:.2f} keV/decay = {:.5f} pJ/decay\n".format(Ca41.heat, Ca41.heat * pJ_in_keV))
99Tc → 99Ru¶
β- decay with half-life 0.2111 ± 0.0012 Myr (NuDat references Browne & Tuli 2017 NDS)
Second forbidden non-unique transition, shape factor from Mougeot 2015 who refers to Reich & Schüpferling (1974)
In [30]:
# decay-specific inputs
mass_parent = 98.9062508 # u [NIST]
mass_daughter = 98.9059341 # u [NIST]
Z = 44 # atomic number of daughter
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
# decay to ground state
branch = 0.999984
Eendpoint = Q
q = Emesh(Eendpoint,dE)
E = Eendpoint - q
w = me + E
p = np.sqrt(w**2-me**2)
shape = q**2 + 0.54*p**2 # Mougeot 2015 <- Reich & Schüpferling 1974
dNdE = normalize_dydx(p*w*q**2 * fermi_func(w,Z) * shape, dE, branch)
anu_mean = mean_x(dNdE, q, dE)
el_mean = Eendpoint - anu_mean
Tc99_anu_away0 = anu_mean * branch
Tc99_heat0 = el_mean * branch
plt.semilogy(q, dNdE)
plt.xlabel('E (keV)')
plt.ylabel('dN/dE (1/keV/decay)')
plt.show()
print ('Q = {:.1f} keV (NuDat: 297.5 ± 1.0 keV)\n'.format(Q))
print ("beta- decay to ground state:")
print (" Neutrinos per decay =", int_y(dNdE, dE), "(check of normalization)")
print (" Mean electron energy = {:.2f} keV/decay (NuDat: 84.6 ± 0.5 keV, Mougeot: 95.28/101.03 keV)".format(el_mean))
print (" Mean nuebar energy = {:.1f} keV/decay".format(anu_mean))
print (" Carried away by nuebar = {:.1f} keV/decay".format(Tc99_anu_away0))
print (" Heats the Earth = {:.3f} keV/decay = {:.5f} pJ/decay\n".format(Tc99_heat0, Tc99_heat0 * pJ_in_keV))
# decay to excited state
branch = 0.000016
Eexc = 89.50
Eendpoint = Q - Eexc
el_mean = 81.7 # NuDat, resp. Browne & Tuli 2017 NDS
Tc99_heat1 = el_mean * branch
print ("beta- decay to excited state {:.2f} kev:".format(Eexc))
print (" Neutrinos per decay =", branch)
print (" Mean electron energy = {:.2f} keV/decay (from NuDat, resp. Browne & Tuli 2017 NDS)".format(el_mean))
print (" Heats the Earth = {:.3f} keV/decay = {:.7f} pJ/decay\n".format(Tc99_heat1, Tc99_heat1 * pJ_in_keV))
Tc99 = Nuclide("99Tc")
Tc99.halflife = 2.111e5
Tc99.mass_atomic = 101.07 # Ru atomic mass
Tc99.mass_nuclide = mass_parent
Tc99.isotfrac = 3.9e-5 * 0.1260 # 99Tc/100Ru * 100Ru/Ru
Tc99.timeref = tcai
Tc99.Qbeta = Q
Tc99.heat = Tc99_heat0 + Tc99_heat1
Tc99.lanu = 1
Tc99.lnu = 0
Tc99.heat_gamma = 89.5 * 6.5e-6
print ("Overall heats the Earth = {:.3f} keV/decay = {:.7f} pJ/decay\n".format(Tc99.heat, Tc99.heat * pJ_in_keV))
81Kr → 81Br¶
ε decay (EC only; no β+) with half-life 0.229 ± 0.011 Myr (NuDat references Baglin 2008 NDS)
In [31]:
# decay-specific inputs
mass_parent = 80.9165912 # u [NIST]
mass_daughter = 80.9162897 # u [NIST]
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
print ('Q = {:.1f} keV (NuDat: 280.8 ± 0.5 keV)\n'.format(Q))
# EC to excited state
branch_exc = 0.0030
Eexc = 275.991
Kr81 = Nuclide("81Kr")
Kr81.halflife = 2.29e5
Kr81.mass_nuclide = mass_parent
Kr81.mass_atomic = -999
Kr81.isotfrac = -0
Kr81.timeref = tcai
Kr81.Qec = Q
Kr81.heat = branch_exc * Eexc
Kr81.lanu = 0
Kr81.lnu = 1
#Kr81.heat_gamma = 275.990 * 0.00298
gam_lev = np.array([ 1.48, 11.878, 11.924, 13.284, 13.292, 13.469, 275.990 ])
gam_int = np.array([ 2.18, 15.5, 30.1, 2.11, 4.09, 0.409, 0.298 ]) / 100
Kr81.heat_gamma = sum(gam_lev*gam_int)
print(Kr81.heat_gamma)
print ("Heats the Earth = {:.3f} keV/decay = {:.6f} pJ/decay".format(Kr81.heat, Kr81.heat * pJ_in_keV))
126Sn → 126Sb¶
β- decay with half-life 0.230 ± 0.014 Myr (NuDat referencing Katakura & Kitao 2002 NDS)
Second forbidden non-unique transition. Shape factor?? ... simply using shape factor = 1
In [32]:
# decay-specific inputs
mass_parent = 125.907659 # u [NIST]
mass_daughter = 125.907253 # u [NIST]
Z = 51 # atomic number of daughter
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
print ('Q = {:.1f} keV (NuDat: 378 ± 30 keV)\n'.format(Q))
# NuDat inputs
branch = 1.0
Eexc = 127.9
Eendpoint = Q - Eexc
print ('Endpoint energy = {:.1f} keV (NuDat: 250 ± 30 keV)\n'.format(Eendpoint))
q = Emesh(Eendpoint,dE)
E = Eendpoint - q
w = me + E
p = np.sqrt(w**2-me**2)
shape = 1
dNdE = normalize_dydx(p*w*q**2 * fermi_func(w,Z) * shape, dE, branch)
anu_mean = mean_x(dNdE, q, dE)
el_mean = Eendpoint - anu_mean
print ("Mean electron energy = {:.2f} keV (NuDat: 70 ± 14 keV)".format(el_mean))
anu_away = anu_mean * branch
# there is branching at 17.7 keV level...
Eexc1 = 17.7
branch_it = 0.14
branch_beta = 0.86
deexcite = branch_it * Eexc + branch_beta * (Eexc-Eexc1)
heat = (el_mean + deexcite) * branch
Sn126_heat1 = heat
print ("Heats the Earth = {:.1f} keV/decay = {:.5f} pJ/decay".format(Sn126_heat1, Sn126_heat1 * pJ_in_keV))
#gam_lev = np.array([ 17.7, 21.646, 22.70, 23.280, 42.641, 64.281, 86.938, 87.567 ])
#gam_int = np.array([ 4.6E-5, 1.26, 0.100, 6.4, 0.50, 9.6, 8.9, 37 ]) / 100
#Sn126_heat_gamma1 = sum(gam_lev*gam_int)
gam_lev = np.array([ 3.6, 17.7, 21.646, 22.70, 23.280, 26.111, 26.359, 29.679, 29.726, 30.393, 42.641, 64.281, 86.938, 87.567 ])
gam_int = np.array([ 11.4, 4.6E-5, 1.26, 0.100, 6.4, 8.4, 15.5, 1.42, 2.73, 0.77, 0.50, 9.6, 8.9, 37 ]) / 100
Sn126_heat_gamma1 = sum(gam_lev*gam_int)
In [33]:
# decay-specific inputs
mass_parent = 125.907253 # u [NIST]
mass_daughter = 125.9033109 # u [NIST]
Z = 51 # atomic number of daughter
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
print ('Q = {:.1f} keV (NuDat: 3673 ± 33 keV)\n'.format(Q))
# NuDat inputs
branch = 1.
el_mean_lev = np.array([ 55, 62, 146, 154, 191, 221, 244, 290, 295, 309, 419, 462, 539, 735 ])
el_endp_lev = np.array([ 2.0E+2, 2.2E+2, 4.8E+2, 5.0E+2, 6.0E+2, 6.8E+2, 7.0E+2, 8.6E+2, 8.3E+2, 9.1E+2, 1.18E+3, 1.28E+3, 1.45E+3, 1.90E+3 ])
branch_lev = np.array([ 0.50, 2.09, 29, 5.9, 8.4, 4.2, 0.50, 8.1, 0.8, 4.9, 16, 0.9, 3.0, 20 ]) / 100
print ("Sum of branching ratios = {:.4f}".format(sum(branch_lev)))
print (" renormalizing to branching {:.4f} (NuDat)".format(branch))
norm = sum(branch_lev)
branch_lev *= branch/norm
print ("Sum of branching ratios = {:.4f}".format(sum(branch_lev)))
el_mean = sum(el_mean_lev*branch_lev)
print ("Mean electron energy = {:.1f} keV (NuDat: 340 ± 60 keV)".format(el_mean))
deexcite = sum((Q - el_endp_lev) * branch_lev)
print ("Deexcitation energy per decay = {:.1f} keV".format(deexcite))
heat = el_mean + deexcite
Sn126_heat2 = heat * branch_it
print ("Heats the Earth = {:.1f} keV/decay = {:.5f} pJ/decay".format(Sn126_heat2, Sn126_heat2 * pJ_in_keV))
#gam_lev = np.array([ 148.7, 208.6, 223.9, 278.2, 296.5, 297.1, 414.7, 415.3, 556.3, 573.9, 593.2, 619.9, 638.8, 656.3, 666.5, 674.8, 695.0, 697.0, 720.7, 856.8, 953.7, 958.3, 989.6, 1036.2, 1064.4, 1213.3, 1476.9 ])
#gam_int = np.array([ 0.40, 0.50, 1.39, 2.4, 4.5, 0.50, 83.3, 1.0, 1.69, 6.7, 7.5, 0.90, 0.90, 2.19, 99.60, 3.7, 99.60, 29, 53.8, 17.6, 1.20, 0.50, 6.8, 1.00, 0.9, 2.39, 0.28 ]) / 100
#Sn126_heat_gamma2 = sum(gam_lev*gam_int) * branch_it
gam_lev = np.array([ 3.77, 27.202, 27.472, 30.944, 30.995, 31.704, 148.7, 208.6, 223.9, 278.2, 296.5, 297.1, 414.7, 415.3, 556.3, 573.9, 593.2, 619.9, 638.8, 656.3, 666.5, 674.8, 695.0, 697.0, 720.7, 856.8, 953.7, 958.3, 989.6, 1036.2, 1064.4, 1213.3, 1476.9 ])
gam_int = np.array([ 0.176, 0.482, 0.89, 0.082, 0.158, 0.0457, 0.40, 0.50, 1.39, 2.4, 4.5, 0.50, 83.3, 1.0, 1.69, 6.7, 7.5, 0.90, 0.90, 2.19, 99.60, 3.7, 99.60, 29, 53.8, 17.6, 1.20, 0.50, 6.8, 1.00, 0.9, 2.39, 0.28 ]) / 100
Sn126_heat_gamma2 = sum(gam_lev*gam_int) * branch_it
In [34]:
# decay-specific inputs
mass_parent = 125.907253 # u [NIST]
mass_daughter = 125.9033109 # u [NIST]
Z = 51 # atomic number of daughter
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
print ('Q = {:.1f} keV (NuDat: 3673 ± 33 keV)\n'.format(Q))
# NuDat inputs
branch = branch_beta
Estart = Eexc1
el_mean_lev = np.array([ 287, 341, 470, 743 ])
el_endp_lev = np.array([ 8.5e2, 9.9e2, 1.30e3, 1.92e3 ])
branch_lev = np.array([ 0.86, 1.3, 3.4, 83 ]) / 100
print ("Sum of branching ratios = {:.4f}".format(sum(branch_lev)))
#print (" renormalizing to branching {:.4f} (NuDat)".format(branch))
#norm = sum(branch_lev)
#branch_lev *= branch/norm
#print ("Sum of branching ratios = {:.4f}".format(sum(branch_lev)))
el_mean = sum(el_mean_lev*branch_lev)
print ("Mean electron energy = {:.1f} keV (NuDat: 720 ± 90 keV)".format(el_mean))
deexcite = sum((Q + Estart - el_endp_lev) * branch_lev)
print ("Deexcitation energy per decay = {:.1f} keV".format(deexcite))
heat = el_mean + deexcite
Sn126_heat3 = heat
print ("Heats the Earth = {:.1f} keV/decay = {:.5f} pJ/decay".format(Sn126_heat3, Sn126_heat3 * pJ_in_keV))
#gam_lev = np.array([ 414.5, 620.0, 666.1, 694.8, 928.2, 1034.9, 1061.6, 1476.1 ])
#gam_int = np.array([ 86, 1.54, 86, 82, 1.3, 1.80, 0.51, 0.34 ]) / 100
#Sn126_heat_gamma3 = sum(gam_lev*gam_int)
gam_lev = np.array([ 3.77, 27.202, 27.472, 30.944, 30.995, 31.704, 414.5, 620.0, 666.1, 694.8, 928.2, 1034.9, 1061.6, 1476.1 ])
gam_int = np.array([ 0.139, 0.386, 0.71, 0.066, 0.127, 0.0366, 86, 1.54, 86, 82, 1.3, 1.80, 0.51, 0.34 ]) / 100
Sn126_heat_gamma3 = sum(gam_lev*gam_int)
126Sn → 126Sb → 126Te decay summary¶
In [35]:
Sn126 = Nuclide("126Sn")
Sn126.halflife = 2.35e5
Sn126.mass_atomic = 118.710 # Sn atomic mass
Sn126.isotfrac = 3e-6 * 0.0579 # 126Sn/124Sn * 124Sn/Sn
Sn126.timeref = tcai
Sn126.heat = Sn126_heat1 + Sn126_heat2 + Sn126_heat3
Sn126.lanu = 2
Sn126.lnu = 0
Sn126.heat_gamma = Sn126_heat_gamma1 + Sn126_heat_gamma2 + Sn126_heat_gamma3
print ("126Sn → 126Sb rad. heat = {:6.1f} keV/decay = {:.5f} pJ/decay".format(Sn126_heat1, Sn126_heat1 * pJ_in_keV))
print ("126Sb(GS) → 126Te rad. heat = {:6.1f} keV/decay = {:.5f} pJ/decay".format(Sn126_heat2, Sn126_heat2 * pJ_in_keV))
print ("126Sb(17.7) → 126Te rad. heat = {:6.1f} keV/decay = {:.5f} pJ/decay\n".format(Sn126_heat3, Sn126_heat3 * pJ_in_keV))
print ("Overall heats the Earth = {:.1f} keV/decay = {:.5f} pJ/decay".format(Sn126.heat, Sn126.heat * pJ_in_keV))
36Cl decay¶
β– decay to 36Ar (98.1%) or ε decay to 36S (1.9%) with half-life of 0.301 ± 0.002 Myr (NuDat referencing Nica et al 2012 NDS)
36Cl → 36Ar β– decay¶
Ground state to ground state decay. Second forbidden nonunique transition, shape factor from Mougeot 2015 who refers to Rotzinger et al. 2008
In [36]:
mass_parent = 35.968306809 # u [NIST]
mass_daughter = 35.967545105 # u [NIST]
Z = 18 # atomic number of daughter
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
print ('Q = {:.3f} keV (NuDat: 709.547 ± 0.046 keV)'.format(Q))
branch = 0.9810
Eendpoint = Q
q = Emesh(Eendpoint,dE)
E = Eendpoint - q
w = me + E
p = np.sqrt(w**2-me**2)
#
Cl36 = Nuclide("36Cl")
Cl36.halflife = 3.01e5
Cl36.mass_nuclide = mass_parent
Cl36.mass_atomic = (35.446+35.457)/2. # Cl atomic mass
Cl36.isotfrac = 1.9e-8 * 0.7576 # 36Cl/35Cl * 35Cl/Cl
Cl36.timeref = tcai
Cl36.Qbeta = Q
Cl36.branch_beta = branch
#
print ("\nUniform (no) shape factor:")
shape = 1
dNdE = normalize_dydx(p*w*q**2 * fermi_func(w,Z) * shape, dE, branch)
anu_mean = mean_x(dNdE, q, dE)
el_mean = Eendpoint - anu_mean
print (" Mean electron energy = {:.2f} keV (NuDat: 251.33 keV)".format(el_mean))
Cl36.heat_beta = el_mean * branch
Cl36.lanu = branch
Cl36.lnu = 1-branch
print (" Heats the Earth = {:.1f} keV/decay = {:.5f} pJ/decay".format(Cl36.heat_beta, Cl36.heat_beta * pJ_in_keV))
#
print ("\nExperimental shape factor from Mougeot:")
shape = 1. - 0.282715*w - 0.045988/w - 0.491739*w**2 + 0.438731*w**3
dNdE = normalize_dydx(p*w*q**2 * fermi_func(w,Z) * shape, dE, branch)
anu_mean = mean_x(dNdE, q, dE)
el_mean = Eendpoint - anu_mean
print (" Mean electron energy = {:.2f} keV (NuDat: 251.33 keV)".format(el_mean))
Cl36.heat_beta = el_mean * branch
Cl36.lanu = branch
Cl36.lnu = 1-branch
Cl36_heat_gamma1 = 0.
print (" Heats the Earth = {:.1f} keV/decay = {:.5f} pJ/decay".format(Cl36.heat_beta, Cl36.heat_beta * pJ_in_keV))
In [37]:
mass_parent = 35.968306809 # u [NIST]
mass_daughter = 35.96708071 # u [NIST]
Z = -16 # atomic number of daughter
# electron capture
branch = 0.0189
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
print ('Q(EC) = {:.2f} keV (NuDat: 1142.14 ± 0.19 keV)'.format(Q))
Cl36.Qec = Q
Cl36.heat_ec = 0.
Cl36.branch_ec = branch
# beta+
branch = 0.00014
Q = Q - 2*me
print ('Q(beta+) = {:.2f} keV (NuDat: 120.14 ± 0.19 keV)'.format(Q))
Cl36.Qbetaplus = Q
Cl36.pos_mean = 50.24
Cl36.heat_betaplus = (Cl36.pos_mean +2*me) * branch
Cl36.branch_betaplus = branch
print ("Heats the Earth (beta+) = {:.4f} keV/decay = {:.7f} pJ/decay".format(Cl36.heat_betaplus, Cl36.heat_betaplus * pJ_in_keV))
gam_lev = np.array([ 2.307, 2.308, 2.464, 2.464, 511.0 ])
gam_int = np.array([ 0.042, 0.085, 0.0039, 0.0020, 0.0280 ]) / 100
Cl36_heat_gamma2 = sum(gam_lev*gam_int)
36Cl decay summary¶
In [38]:
Cl36.heat = Cl36.heat_beta + Cl36.heat_ec + Cl36.heat_betaplus
Cl36.heat_gamma = Cl36_heat_gamma1 + Cl36_heat_gamma2
print ("Heats the Earth = {:.2f} keV/decay = {:.5f} pJ/decay".format(Cl36.heat, Cl36.heat * pJ_in_keV))
79Se → 79Br¶
β– decay to ground state of 79Br with half-life of 0.326 ± 0.028 Myr (NuDat referencing Singh 2016 NDS)
First forbidden unique transition. Shape factor?? ... using a symmetrical one from Bielajew.
In [39]:
mass_parent = 78.91849929 # u [NIST]
mass_daughter = 78.9183376 # u [NIST]
Z = 35 # atomic number of daughter
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
print ('Q = {:.3f} keV (NuDat: 150.6 ± 1.3 keV)'.format(Q))
branch = 1.
Eendpoint = Q
q = Emesh(Eendpoint,dE)
E = Eendpoint - q
w = me + E
p = np.sqrt(w**2-me**2)
#
Se79 = Nuclide("79Se")
Se79.halflife = 3.26e5
Se79.mass_atomic = -999
Se79.isotfrac = -0
Se79.timeref = tcai
#
print ("\nUniform (no) shape factor:")
shape = 1
dNdE = normalize_dydx(p*w*q**2 * fermi_func(w,Z) * shape, dE, branch)
anu_mean = mean_x(dNdE, q, dE)
el_mean = Eendpoint - anu_mean
print (" Mean electron energy = {:.2f} keV (NuDat: 52.8 keV)".format(el_mean))
Se79.anu_away = anu_mean * branch
Se79.heat = el_mean * branch
print (" Heats the Earth = {:.1f} keV/decay = {:.5f} pJ/decay".format(Se79.heat, Se79.heat * pJ_in_keV))
#
print ("\nSymmetrical shape factor from Bielajew")
shape = p**2 + q**2 # symmetrical, from Bielajew
dNdE = normalize_dydx(p*w*q**2 * fermi_func(w,Z) * shape, dE, branch)
anu_mean = mean_x(dNdE, q, dE)
el_mean = Eendpoint - anu_mean
print (" Mean electron energy = {:.2f} keV (NuDat: 52.8 keV)".format(el_mean))
Se79.anu_away = anu_mean * branch
Se79.heat = el_mean * branch
Se79.lanu = 1
Se79.lnu = 0
Se79.heat_gamma = 0.
print (" Heats the Earth = {:.1f} keV/decay = {:.5f} pJ/decay".format(Se79.heat, Se79.heat * pJ_in_keV))
26Al → 26Mg¶
ε decay with half-life 0.717 ± 0.024 Myr (NuDat references Basunia & Hurst 2016 NDS)
β+ to 1808.72 level (81.73%) and EC to 2938.41 (2.74%) and 1808.72 (15.51%) levels (total 99.98%)
Using theoretical formula for shape factor, 2nd forbidden unique transition
In [40]:
# decay-specific inputs
mass_parent = 25.986891904 # u [NIST]
mass_daughter = 25.982592968 # u [NIST]
Z = -12 # atomic number of daughter
Al26 = Nuclide("26Al")
Al26.halflife = 7.17e5
Al26.mass_nuclide = mass_parent
Al26.mass_atomic = 26.9815385 # Al atomic mass
Al26.isotfrac = 5.25e-5 * 1. # 26Al/27Al * 27Al/Al
Al26.timeref = tcai
# electron capture
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
print ('Q(EC) = {:.2f} keV (NuDat: 4004.43 ± 0.06 keV)'.format(Q))
branch_exc = np.array([0.1551, 0.0274])
E_exc = np.array([1808.72, 2938.41])
Al26.heat_ec = sum(branch_exc * E_exc)
print ("Heats the Earth = {:.2f} keV/decay = {:.5f} pJ/decay\n".format(Al26.heat_ec, Al26.heat_ec * pJ_in_keV))
# beta+
branch = 0.8173
Q = Q - 2*me # keV Q-value
print ('Q(beta+) = {:.2f} keV (NuDat: n/a)'.format(Q))
E_exc = 1808.72
Eendpoint = Q - E_exc
print ('Eendpoint = {:.2f} keV (NuDat: 1173.42 ± 0.09 keV)'.format(Eendpoint))
q = Emesh(Eendpoint,dE)
E = Eendpoint - q
w = me + E
p = np.sqrt(w**2-me**2)
shape = p**4 + 10/3.*p**2*q**2 + q**4 # Bielajew 2nd forbidden unique
dNdE = normalize_dydx(p*w*q**2 * fermi_func(w,Z) * shape, dE, branch)
nu_mean = mean_x(dNdE, q, dE)
print (nu_mean*branch)
pos_mean = Eendpoint - nu_mean
print ("Mean positron energy = {:.2f} keV (NuDat: 543.29 keV)".format(pos_mean))
Al26.heat_beta = (pos_mean + E_exc + 2*me) * branch
print ("Heats the Earth = {:.2f} keV/decay = {:.5f} pJ/decay\n".format(Al26.heat_beta, Al26.heat_beta * pJ_in_keV))
Al26.heat = Al26.heat_ec + Al26.heat_beta
Al26.lanu = 0
Al26.lnu = 1
#Al26.heat_gamma = 1808.65 * 0.9976 + 1129.67 * 0.0250 + 2938 * 0.0024
gam_lev = np.array([ 1.254, 1.254, 511.0, 1129.67, 1808.65, 2938 ])
gam_int = np.array([ 0.33, 0.167, 163.5, 2.50, 99.76, 0.24 ]) / 100
Al26.heat_gamma = sum(gam_lev*gam_int)
print ("Overall heats the Earth = {:.2f} keV/decay = {:.5f} pJ/decay".format(Al26.heat, Al26.heat * pJ_in_keV))
In [41]:
print ("Looking what happens when I switch sign of Z_daughter in Fermi function...")
plt.plot(E, normalize_dydx(p*w*q**2 * fermi_func(w,0) * shape, dE, branch), label='Z=0')
plt.plot(E, normalize_dydx(p*w*q**2 * fermi_func(w,12) * shape, dE, branch), label='Z=+12 (beta-)')
plt.plot(E, normalize_dydx(p*w*q**2 * fermi_func(w,-12) * shape, dE, branch), label='Z=-12 (beta+)')
plt.xlabel('E (keV)')
plt.ylabel('Fermi function')
plt.legend()
plt.show()
10Be → 10B¶
β- decay with half-life 1.51 ± 0.06 Myr (NuDat referencing Tilley et al 2004 Nucl.Phys.A) or half life 1.387 ± 0.012 Myr according to Bill. Using the latter.
Mougeot 2015: second forbidden unique transition, shape factor referred from Feldman & Wu 1952; in fact a symmetrical theoretical shape factor
In [42]:
# decay-specific inputs
mass_parent = 10.013534695 # u [NIST]
mass_daughter = 10.01293695 # u [NIST]
Z = 5 # atomic number of daughter
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
print ('Q = {:.1f} keV (NuDat: 556.0 ± 0.6 keV)\n'.format(Q))
Be10 = Nuclide("10Be")
Be10.halflife = 1.387e6
Be10.mass_nuclide = mass_parent
Be10.mass_atomic = 9.0121831 # Be atomic mass
Be10.isotfrac = 5.3e-4 * 1. # 10Be/9Be * 9Be/Be
Be10.timeref = tcai
# NuDat inputs
branch = 1.
Q = 556.0
el_mean = 202.56
print ("Using NuDat inputs:")
print (" Mean electron energy = {:.2f} keV (from NuDat)\n".format(el_mean))
# calculating spectrum
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
Eendpoint = Q
q = Emesh(Eendpoint,dE)
E = Eendpoint - q
w = me + E
p = np.sqrt(w**2-me**2)
print ("My calculation:")
shape = 1
dNdE = normalize_dydx(p*w*q**2 * fermi_func(w,Z) * shape, dE, branch)
anu_mean = mean_x(dNdE, q, dE)
el_mean = Eendpoint - anu_mean
print (" Mean electron energy w/ shape factor equal to 1 = {:.2f} keV".format(el_mean))
shape = p**4 + 10/3.*p**2*q**2 + q**4 # Mougeot 2015 <- Feldman & Wu 1952
dNdE = normalize_dydx(p*w*q**2 * fermi_func(w,Z) * shape, dE, branch)
anu_mean = mean_x(dNdE, q, dE)
el_mean = Eendpoint - anu_mean
print (" Mean electron energy w/ shape from Mougeot = {:.2f} keV".format(el_mean))
anu_away = anu_mean * branch
heat = el_mean * branch
**There is a 50 keV difference between *NuDat* and *Mougeot 2015* !!!**
I will use Mougeot 2015, resp. my calculation:
In [43]:
Be10.Qbeta = Q
Be10.heat = heat
Be10.lanu = 1
Be10.lnu = 0
Be10.heat_gamma = 0.
print ("Heats the Earth = {:.2f} keV/decay = {:.5f} pJ/decay\n".format(Be10.heat, Be10.heat * pJ_in_keV))
93Zr → 93Nb¶
β- decay with half-life 1.61 ± 0.05 Myr (NuDat referencing Baglin 2011), 73% to excited level 1st forbidden unique, 27% to ground state 2nd forbidden non-unique
In [44]:
# decay-specific inputs
mass_parent = 92.9064699 # u [NIST]
mass_daughter = 92.9063730 # u [NIST]
Z = 41 # atomic number of daughter
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
print ('Q = {:.1f} keV (NuDat: 90.8 ± 1.6 keV)\n'.format(Q))
Zr93 = Nuclide("93Zr")
Zr93.halflife = 1.61e6
Zr93.mass_atomic = -999
Zr93.isotfrac = -0
Zr93.timeref = tcai
Zr93.Qbeta = Q
print ("Mean electron energy from Nudat = 20 ± 3 keV (from NuDat)\n".format(el_mean))
# calculating spectrum
branch1 = 0.73
Elev1 = 30.8
Eendpoint = Q - Elev1
q = Emesh(Eendpoint,dE)
E = Eendpoint - q
w = me + E
p = np.sqrt(w**2-me**2)
shape = p**2 + q**2
dNdE = normalize_dydx(p*w*q**2 * fermi_func(w,Z) * shape, dE, branch)
anu_mean = mean_x(dNdE, q, dE)
el_mean1 = Eendpoint - anu_mean
#
branch0 = 0.27
Elev0 = 0.
Eendpoint = Q - Elev0
q = Emesh(Eendpoint,dE)
E = Eendpoint - q
w = me + E
p = np.sqrt(w**2-me**2)
shape = 1.
dNdE = normalize_dydx(p*w*q**2 * fermi_func(w,Z) * shape, dE, branch)
anu_mean = mean_x(dNdE, q, dE)
el_mean0 = Eendpoint - anu_mean
#
el_mean = el_mean0 * branch0 + el_mean1 * branch1
heat = (el_mean0 + Elev0) * branch0 + (el_mean1 + Elev1) * branch1
Zr93.heat = heat
Zr93.lanu = 1
Zr93.lnu = 0
#Zr93.heat_gamma = 30.77 * 4.3E-6
gam_lev = np.array([ 2.17, 16.521, 16.615, 18.607, 18.623, 18.952, 30.77 ])
gam_int = np.array([ 2.13, 2.45, 4.7, 0.37, 0.71, 0.161, 4.3E-4 ]) / 100
Zr93.heat_gamma = sum(gam_lev*gam_int)
print ("Mean electron energy with both shape factors equal to one = 17.88 keV".format(el_mean))
print ("Mean electron energy calculated theoretical excited level SF = {:.2f} keV (using this)\n".format(el_mean))
print ("Heats the Earth = {:.1f} keV/decay = {:.4f} pJ/decay".format(Zr93.heat, Zr93.heat * pJ_in_keV))
In [45]:
mass_parent = 153.9244293 # u [NIST]
mass_daughter = 149.9186644 # u [NIST]
Q = (mass_parent - mass_daughter - mass_4He_u) * keV_per_u # keV Q-value
print ('Q = {:.2f} keV (NuDat: 2945 ± 5 keV)\n'.format(Q))
Dy154 = Nuclide("154Dy")
Dy154.halflife = 3.0e6
Dy154.mass_atomic = -999
Dy154.isotfrac = -0
Dy154.timeref = tcai
Dy154.Qalpha = Q
Dy154.heat = Q
Dy154.lanu = 0
Dy154.lnu = 0
Dy154.heat_gamma = 0.
print ("Heats the Earth = {:.1f} keV/decay = {:.4f} pJ/decay".format(Dy154.heat, Dy154.heat * pJ_in_keV))
In [46]:
mass_parent = 149.9186644 # u [NIST]
mass_daughter = 145.9130470 # u [NIST]
Q = (mass_parent - mass_daughter - mass_4He_u) * keV_per_u # keV Q-value
print ('Q = {:.2f} keV (NuDat: 2808 ± 6 keV)\n'.format(Q))
Gd150 = Nuclide("150Gd")
Gd150.halflife = 1.79e6
Gd150.mass_atomic = -999
Gd150.isotfrac = -0
Gd150.timeref = tcai
Gd150.Qalpha = Q
Gd150.heat = Q
Gd150.lanu = 0
Gd150.lnu = 0
Gd150.heat_gamma = 0.
print ("Heats the Earth = {:.1f} keV/decay = {:.4f} pJ/decay".format(Gd150.heat, Gd150.heat * pJ_in_keV))
In [47]:
# decay-specific inputs
mass_parent = 145.9130470 # u [NIST]
mass_daughter = 141.9077290 # u [NIST]
Q = (mass_parent - mass_daughter - mass_4He_u) * keV_per_u # keV Q-value
print ('Q = {:.1f} keV (NuDat: 2528 ± 3 keV)\n'.format(Q))
Sm146 = Nuclide("146Sm")
Sm146.halflife = 1.03e8
Sm146.mass_atomic = 150.36 # Sm atomic mass
Sm146.isotfrac = 8.28e-3 * 0.0307 # 146Sm/144Sm * 144Sm/Sm
Sm146.timeref = tcai
Sm146.heat = Q
Sm146.lanu = 0
Sm146.lnu = 0
Sm146.heat_gamma = 0.
print ("Heats the Earth = {:.1f} keV/decay = {:.4f} pJ/decay".format(Sm146.heat, Sm146.heat * pJ_in_keV))
135Cs → 135Ba¶
β- decay with half-life 2.3 ± 0.3 Myr (NuDat refering to Singh et al 2008 NDS)
Mougeot 2015: second forbidden nonunique transition, shape factor from Lidofsky et al 1953 Phys.Rev.90, p.387 (I am unable to find the PDF)
In [48]:
# decay-specific inputs
mass_parent = 134.9059770 # u [NIST]
mass_daughter = 134.90568838 # u [NIST]
Z = 56 # atomic number of daughter
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
print ('Q = {:.1f} keV (NuDat: 268.7 ± 1.1 keV)\n'.format(Q))
# NuDat inputs
branch = 1.
Q = 268.7
el_mean = 75.7
print ("Using NuDat inputs:")
print (" Mean electron energy = {:.1f} keV (from NuDat)\n".format(el_mean))
# calculating spectrum
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
Eendpoint = Q
q = Emesh(Eendpoint,dE)
E = Eendpoint - q
w = me + E
p = np.sqrt(w**2-me**2)
print ("My calculation:")
shape = 1
dNdE = normalize_dydx(p*w*q**2 * fermi_func(w,Z) * shape, dE, branch)
anu_mean = mean_x(dNdE, q, dE)
el_mean = Eendpoint - anu_mean
print (" Mean electron energy w/ shape factor equal to 1 = {:.2f} keV".format(el_mean))
shape = q**2 + 0.10*p**2 # Mougeot 2015 <- Lidofsky et al 1953
dNdE = normalize_dydx(p*w*q**2 * fermi_func(w,Z) * shape, dE, branch)
anu_mean = mean_x(dNdE, q, dE)
el_mean = Eendpoint - anu_mean
print (" Mean electron energy w/ experimental shape factor from Mougeot = {:.2f} keV (Mougeot: 61.57 keV)".format(el_mean))
anu_away = anu_mean * branch
heat = el_mean * branch
print (" Mean electron energy w/ calc. by Mougeot using ξ approximation = 88.44 keV\n".format(el_mean))
print ("I will use my (~Mougeot) calculation with experimental shape factor:\n")
Cs135 = Nuclide("135Cs")
Cs135.halflife = 2.3e6
Cs135.mass_atomic = 132.90545196 # Cs atomic mass
Cs135.isotfrac = 2.8e-4 * 1. # 135Cs/133Cs * 133Cs/Cs
Cs135.timeref = tcai
Cs135.heat = heat
Cs135.lanu = 1
Cs135.lnu = 0
Cs135.heat_gamma = 0.
print ("Heats the Earth = {:.2f} keV/decay = {:.6f} pJ/decay".format(Cs135.heat, Cs135.heat * pJ_in_keV))
60Fe → 60Co → 60Ni¶
Two-step β- decay with half-lives 2.62 ± 0.04 Myr, then 10.467 ± 0.006 min or 1925.28 ± 0.14 days (NuDat 60Fe, 60Co references Browne & Tuli 2013 NDS)
60Fe(0+) 100% β- decays to 58.603 keV level of 60Co(2+). From then on 99.75% ITs to ground state of 60Co(5+), the remaining 0.25% β- decays to 1332.5 keV (0.24%) or 2158.6 keV (0.0086%) level of 60Ni. 60Co(5+) ground state β- decays to 2505.7 keV (99.88%) or 1332.5 keV (0.12%) level of 60Ni.
First step is second forbidden non-unique transition. Shape factor??
I will use the NuDat (= Browne & Tuli 2013 NDS) inputs.
In [49]:
# decay-specific inputs
mass_parent = 59.9340711 # u [NIST]
mass_daughter = 59.93381630 # u [NIST]
Z = 27 # atomic number of daughter
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
print ('Q = {:.1f} keV (NuDat: 237 ± 3 keV)\n'.format(Q))
Eexc = 58.6
Eendpoint = Q - Eexc
q = Emesh(Eendpoint,dE)
E = Eendpoint - q
w = me + E
p = np.sqrt(w**2-me**2)
shape = 1
dNdE = normalize_dydx(p*w*q**2 * fermi_func(w,Z) * shape, dE, branch)
anu_mean = mean_x(dNdE, q, dE)
el_mean_60fe = Eendpoint - anu_mean
print ("Mean electron energy in 60Fe decay w/ shape factor equal to 1 = {:.2f} keV (NuDat: 50.2 keV)\n".format(el_mean_60fe))
# using NuDat input
E_60co_1 = 58.603
E_60ni_1 = 1332.508
E_60ni_2 = 2158.612
E_60ni_3 = 2505.748
#el_mean_60fe = 50.2
el_mean_60co_0 = 96.41
el_mean_60co_1 = 593.621
branch_it = 0.9975
heat = el_mean_60fe
heat += branch_it * E_60co_1
heat += branch_it * el_mean_60co_0
heat += branch_it * (0.9988 * E_60ni_3 + 0.0012 * E_60ni_1)
heat += (1-branch_it) * el_mean_60co_1
heat += 0.0024 * E_60ni_1 + 0.000086 * E_60ni_2
Fe60 = Nuclide("60Fe")
Fe60.halflife = 2.62e6
Fe60.mass_atomic = 55.845 # Fe atomic mass
Fe60.isotfrac = 1.01e-8 * 0.91754 # 60Fe/56Fe * 56Fe/Fe
Fe60.timeref = tcai
Fe60.heat = heat
Fe60.lanu = 2
Fe60.lnu = 0
#
#Fe60.heat_gamma = 58.603 * 0.020652 + 1173.228 * 0.9985 + 1332.492 * 0.999826
gam_lev = np.array([ 0.78, 6.915, 6.93, 7.649, 7.649, 58.603 ])
gam_int = np.array([ 0.93, 9.1, 18.0, 2.16, 1.11, 2.0652 ]) / 100
Fe60.heat_gamma = sum(gam_lev*gam_int)
gam_lev = np.array([ 0.85, 7.461, 7.478, 8.265, 8.265, 347.14, 826.10, 1173.228, 1332.492, 2158.57, 2505.692 ])
gam_int = np.array([ 3.29E-4, 0.00322, 0.0063, 7.6E-4, 3.91E-4, 0.0075, 0.0076, 99.85, 99.9826, 0.00120, 2.0E-6 ]) / 100
Fe60.heat_gamma += sum(gam_lev*gam_int)
#
print ("Overall heats the Earth = {:.1f} keV/decay = {:.5f} pJ/decay".format(Fe60.heat, Fe60.heat * pJ_in_keV))
Castillo-Rogez et al. 2009 Icarus say 2712 keV per decay. However, they use an earlier inputs (Tuli 2003 NDS) compared to us (Browne & Tuli 2013 NDS). That said, I have not compared the 2003 vs. 2013 versions of nuclear data sheets for the 60Fe → 60Co → 60Ni decay in question.
Mougeot 2015 adresses some of the branches (β- decays of 60Co with endpoint energies of 318 keV and 1492 keV). I have not compared...
53Mn → 53Cr¶
ε decay (EC only; no β+) with half-life 3.7 ± 0.4 Myr, to ground state [NuDat references Junde 2009 NDS]
In [50]:
mass_parent = 52.94128889 # u [NIST]
mass_daughter = 52.94064815 # u [NIST]
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
print ('Q = {:.1f} keV (NuDat: 596.8 ± 0.4 keV)\n'.format(Q))
Mn53 = Nuclide("53Mn")
Mn53.halflife = 3.74e6
Mn53.mass_atomic = 54.938044 # Mn atomic mass
Mn53.isotfrac = 6.28e-6 * 1. # 53Mn/55Mn * 55Mn/Mn
Mn53.timeref = tcai
Mn53.heat = 0.
Mn53.lanu = 0
Mn53.lnu = 1
#Mn53.heat_gamma = 0.
gam_lev = np.array([ 0.57, 5.405, 5.415, 5.947, 5.947 ])
gam_int = np.array([ 0.37, 7.4, 14.6, 1.64, 0.84 ]) / 100
Mn53.heat_gamma = sum(gam_lev*gam_int)
print ("Heats the Earth = {:.2f} keV/decay = {:.5f} pJ/decay".format(Mn53.heat, Mn53.heat * pJ_in_keV))
98Tc → 98Ru¶
β- decay with half-life 4.2 ± 0.3 Myr to 745.35 keV level of 98Ru (NuDat references Singh & Hu 2003 NDS)
Second forbidden non-unique transition. Shape factor??
In [51]:
# decay-specific inputs
mass_parent = 97.9072124 # u [NIST]
mass_daughter = 97.9052868 # u [NIST]
Z = 44 # atomic number of daughter
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
print ('Q = {:.0f} keV (NuDat: 1796 ± 7 keV)\n'.format(Q))
# NuDat inputs
branch = 1.
Q = 1796
el_mean = 118
print ("Using NuDat inputs:")
print (" Mean electron energy = {:.0f} keV (from NuDat)\n".format(el_mean))
# calculating spectrum
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
Eexc = 1397.77
Eendpoint = Q - Eexc
q = Emesh(Eendpoint,dE)
E = Eendpoint - q
w = me + E
p = np.sqrt(w**2-me**2)
print ("My calculation:")
shape = 1
dNdE = normalize_dydx(p*w*q**2 * fermi_func(w,Z) * shape, dE, branch)
anu_mean = mean_x(dNdE, q, dE)
el_mean = Eendpoint - anu_mean
heat = (el_mean + Eexc) * branch
print (" Electron endpoint energy = {:.2f} keV (NuDat: 397 ± 22 keV)".format(Eendpoint))
print (" Mean electron energy w/ shape factor equal to 1 = {:.2f} keV\n".format(el_mean))
#shape = ? # ???
#dNdE = normalize_dydx(p*w*q**2 * fermi_func(w,Z) * shape, dE, branch)
#anu_mean = mean_x(dNdE, q, dE)
#el_mean = Eendpoint - anu_mean
#print (" Mean electron energy w/ shape factor from ??? = {:.2f} keV".format(el_mean))
#anu_away = anu_mean * branch
#heat = (el_mean + Eexc) * branch
print ("Using my calculation with shape factor equal to one (...can we do better??):\n")
Tc98 = Nuclide("98Tc")
Tc98.halflife = 4.2e6
Tc98.mass_atomic = 101.07 # Ru atomic mass
Tc98.isotfrac = 2e-5 * 0.0554 # 98Tc/96Ru * 96Ru/Ru
Tc98.timeref = tcai
Tc98.heat = heat
Tc98.lanu = 1
Tc98.lnu = 0
Tc98.heat_gamma = 652.41 * 1.00 + 745.35 * 1.02
print ("Heats the Earth = {:.2f} keV/decay = {:.6f} pJ/decay".format(Tc98.heat, Tc98.heat * pJ_in_keV))
97Tc → 97Mo¶
ε decay (EC only; no β+) with half-life 4.21 ± 0.16 Myr, to ground state (NuDat references Nica 2010 NDS)
In [52]:
mass_parent = 96.9063667 # u [NIST]
mass_daughter = 96.90601812 # u [NIST]
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
print ('Q = {:.1f} keV (NuDat: 320 ± 4 keV)\n'.format(Q))
Tc97 = Nuclide("97Tc")
Tc97.halflife = 4.21e6
Tc97.mass_atomic = 101.07 # Ru atomic mass
Tc97.isotfrac = 4e-4 * 0.0187 # 97Tc/98Ru * 98Ru/Ru
Tc97.timeref = tcai
Tc97.heat = 0.
Tc97.lanu = 0
Tc97.lnu = 1
#Tc97.heat_gamma = 0.
gam_lev = np.array([ 2.29, 17.374, 17.479, 19.59, 19.607, 19.965 ])
gam_int = np.array([ 3.80, 19.2, 36.6, 2.93, 5.67, 1.24 ]) / 100
Tc97.heat_gamma = sum(gam_lev*gam_int)
print ("Heats the Earth = {:.2f} keV/decay = {:.5f} pJ/decay".format(Tc97.heat, Tc97.heat * pJ_in_keV))
107Pd → 107Ag¶
β- decay with half-life 6.5 ± 0.3 Myr to 745.35 keV level of 98Ru (NuDat references Blachot 2008 NDS)
First forbidden unique transition. Shape factor?? (...check out Kostensalo et al 2007?)
In [53]:
# decay-specific inputs
mass_parent = 106.9051282 # u [NIST]
mass_daughter = 106.9050916 # u [NIST]
Z = 47 # atomic number of daughter
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
print ('Q = {:.1f} keV (NuDat: 34.1 ± 2.7 keV)\n'.format(Q))
# NuDat inputs
branch = 1.
Q = 34.1
el_mean = 9.3
print ("Mean electron energy from NuDat = {:.0f} keV\n".format(el_mean))
# calculating spectrum
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
Eendpoint = Q
print ("Electron endpoint energy = {:.2f} keV (NuDat: 34 ± 3 keV)\n".format(Eendpoint))
q = Emesh(Eendpoint,dE)
E = Eendpoint - q
w = me + E
p = np.sqrt(w**2-me**2)
#
print ("My calculation:")
#
shape = 1
dNdE = normalize_dydx(p*w*q**2 * fermi_func(w,Z) * shape, dE, branch)
anu_mean = mean_x(dNdE, q, dE)
el_mean = Eendpoint - anu_mean
print (" Mean electron energy w/ shape factor equal to one = {:.2f} keV".format(el_mean))
#
shape = p**2 + q**2
dNdE = normalize_dydx(p*w*q**2 * fermi_func(w,Z) * shape, dE, branch)
anu_mean = mean_x(dNdE, q, dE)
el_mean = Eendpoint - anu_mean
print (" Mean electron energy w/ theoretical shape factor = {:.2f} keV".format(el_mean))
print ("\nUsing my calculation with theoretical shape factor:\n")
anu_away = anu_mean * branch
heat = el_mean * branch
Pd107 = Nuclide("107Pd")
Pd107.halflife = 6.5e6
Pd107.mass_atomic = 106.42 # Pd atomic mass
Pd107.isotfrac = 3.5e-5 * 0.2646 # 107Pd/108Pd * 108Pd/Pd
Pd107.timeref = tcai
Pd107.heat = heat
Pd107.lanu = 1
Pd107.lnu = 0
Pd107.heat_gamma = 0.
print ("Heats the Earth = {:.2f} keV/decay = {:.6f} pJ/decay".format(Pd107.heat, Pd107.heat * pJ_in_keV))
182Hf → 182Ta → 182W¶
Two-step β- decay
182Hf → 182Ta¶
β- decay with half-life 8.90 ± 0.09 Myr (NuDat referencing Singh 2015 NDS and Singh & Roediger 2005 NDS)
Unique first forbidden transition?? Using theoretical shape factor.</font>
In [54]:
# decay-specific inputs
mass_parent = 181.9505612 # u [NIST]
mass_daughter = 181.9501519 # u [NIST]
Z = 73 # atomic number of daughter
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
print ('Q = {:.1f} keV (NuDat: 381 ± 6 keV)\n'.format(Q))
# NuDat inputs
branch = 1.
#Q = 381
Eexc = 270.408
el_mean = 32.9
heat_nudat = (el_mean + Eexc) * branch
print ("Using NuDat inputs:")
print (" Mean electron energy = {:.1f} keV".format(el_mean))
print (" Heats the Earth = {:.1f} keV\n".format(heat_nudat))
# calculating spectrum
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
Eexc = 270.408
Eendpoint = Q - Eexc
q = Emesh(Eendpoint,dE)
E = Eendpoint - q
w = me + E
p = np.sqrt(w**2-me**2)
print ("My calculation:")
shape = p**2 + q**2
dNdE = normalize_dydx(p*w*q**2 * fermi_func(w,Z) * shape, dE, branch)
anu_mean = mean_x(dNdE, q, dE)
el_mean = Eendpoint - anu_mean
heat_me = (el_mean + Eexc) * branch
print (" Electron endpoint energy = {:.2f} keV".format(Eendpoint))
print (" Mean electron energy w/ shape factor equal to 1 = {:.2f} keV".format(el_mean))
print (" Heats the Earth = {:.1f} keV\n".format(heat_me))
print ("Using NuDat-derived value:\n")
Hf182_heat1 = heat_me
print ("Heats the Earth = {:.1f} keV/decay = {:.5f} pJ/decay".format(Hf182_heat1, Hf182_heat1 * pJ_in_keV))
182Ta → 182W¶
β- decay with half-life 114.74 ± 0.12 days (NuDat referencing Singh 2015 NDS and Singh & Roediger 2005 NDS)
Third forbidden non-unique transition. Shape factor??
In [55]:
# decay-specific inputs
mass_parent = 181.9501519 # u [NIST]
mass_daughter = 181.94820394 # u [NIST]
Z = 74 # atomic number of daughter
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
print ('Q = {:.1f} keV (NuDat: 1814.5 ± 1.7 keV)\n'.format(Q))
# NuDat inputs
el_mean_lev = np.array([ 72.54, 85.74, 92.86, 107.08, 129.64, 143.95, 158.27, 169.22, 181.82, 529.12, 625.31 ])
el_endp_lev = np.array([ 261.3, 304.3, 327.0, 371.7, 440.7, 483.4, 525.4, 557.1, 593.1, 1485.1, 1714.4 ])
branch_lev = np.array([ 29.23, 0.142, 1.3, 0.565, 20.1, 2.37, 43.2, 0.05, 3.2, 0.096, 0.058 ]) / 100
print ("Sum of branching ratios = {:.4f}".format(sum(branch_lev)))
el_mean = sum(el_mean_lev*branch_lev)
print ("Mean electron energy = {:.1f} keV (NuDat: 127 ± 3 keV)".format(el_mean))
deexcite = sum((Q - el_endp_lev) * branch_lev)
print ("Deexcitation energy per decay = {:.1f} keV".format(deexcite))
heat = el_mean + deexcite
Hf182_heat2 = heat
print ("Heats the Earth = {:.1f} keV/decay = {:.5f} pJ/decay".format(Hf182_heat2, Hf182_heat2 * pJ_in_keV))
#heat_nudat = (el_mean + Eexc) * branch
#print ("Using NuDat inputs:")
#print (" Mean electron energy = {:.1f} keV".format(el_mean))
#print (" Heats the Earth = {:.1f} keV\n".format(heat_nudat))
182Hf → 182Ta → 182W decay summary¶
In [56]:
Hf182 = Nuclide("182Hf")
Hf182.halflife = 8.90e6
Hf182.mass_atomic = 178.49 # Hf atomic mass
Hf182.isotfrac = 1.018e-4 * 0.3508 # 182Hf/180Hf * 180Hf/Hf
Hf182.timeref = tcai
Hf182.heat = Hf182_heat1 + Hf182_heat2
Hf182.lanu = 2
Hf182.lnu = 0
#
#gam_lev = np.array ([ 97.85, 114.32, 156.09, 172.55, 270.408 ])
#gam_int = np.array ([ 0.110, 3.00, 7.00, 0.200, 79.0 ]) / 100
#Hf182.heat_gamma = sum(gam_lev*gam_int)
#gam_lev = np.array ([ 31.7377, 42.7148, 44.66, 65.72215, 67.74970, 84.68024, 100.10595, 110.393, 113.67170, 116.4179, 121.5, 152.42991, 156.3864, 179.39381, 198.35187, 222.1085, 229.3207, 264.0740, 351.02, 829.9, 891.70, 928.00, 959.73, 1001.700, 1035.7, 1044.42, 1113.410, 1121.290, 1157.302, 1158.1, 1180.85, 1189.040, 1221.395, 1223.60, 1231.004, 1257.407, 1273.719, 1289.145, 1342.730, 1373.824, 1387.390, 1410.13, 1453.120 ])
#gam_int = np.array ([ 0.874, 0.268, 0.030, 3.01, 42.9, 2.654, 14.20, 0.107, 1.871, 0.444, 0.0024, 7.02, 2.671, 3.119, 1.465, 7.57, 3.644, 3.612, 0.0113, 0.014, 0.0574, 0.614, 0.350, 2.086, 0.0067, 0.239, 0.445, 35.24, 0.73, 0.29, 0.087, 16.49, 27.23, 0.24, 11.62, 1.509, 0.660, 1.372, 0.2565, 0.2224, 0.0729, 0.0396, 0.0307 ]) / 100
#Hf182.heat_gamma += sum(gam_lev*gam_int)
#
gam_lev = np.array ([ 8.15, 56.28, 57.535, 64.948, 65.222, 66.982, 97.85, 114.32, 156.09, 172.55, 270.408 ])
gam_int = np.array ([ 5.3, 4.52, 7.8, 0.88, 1.71, 0.59, 0.110, 3.00, 7.00, 0.200, 79.0 ]) / 100
Hf182.heat_gamma = sum(gam_lev*gam_int)
gam_lev = np.array ([ 8.4, 31.7377, 42.7148, 44.66, 57.981, 59.318, 65.72215, 66.95, 67.244, 67.74970, 69.067, 84.68024, 100.10595, 110.393, 113.67170, 116.4179, 121.5, 152.42991, 156.3864, 179.39381, 198.35187, 222.1085, 229.3207, 264.0740, 351.02, 829.9, 891.70, 928.00, 959.73, 1001.700, 1035.7, 1044.42, 1113.410, 1121.290, 1157.302, 1158.1, 1180.85, 1189.040, 1221.395, 1223.60, 1231.004, 1257.407, 1273.719, 1289.145, 1342.730, 1373.824, 1387.390, 1410.13, 1453.120 ])
gam_int = np.array ([ 23.6, 0.874, 0.268, 0.030, 10.01, 17.2, 3.01, 1.96, 3.76, 42.9, 1.31, 2.654, 14.20, 0.107, 1.871, 0.444, 0.0024, 7.02, 2.671, 3.119, 1.465, 7.57, 3.644, 3.612, 0.0113, 0.014, 0.0574, 0.614, 0.350, 2.086, 0.0067, 0.239, 0.445, 35.24, 0.73, 0.29, 0.087, 16.49, 27.23, 0.24, 11.62, 1.509, 0.660, 1.372, 0.2565, 0.2224, 0.0729, 0.0396, 0.0307 ]) / 100
Hf182.heat_gamma += sum(gam_lev*gam_int)
#
print ("182Hf → 182Ta rad. heat = {:.1f} keV/decay = {:.5f} pJ/decay".format(Hf182_heat1, Hf182_heat1 * pJ_in_keV))
print ("182Ta → 182W rad. heat = {:.1f} keV/decay = {:.5f} pJ/decay\n".format(Hf182_heat2, Hf182_heat2 * pJ_in_keV))
print ("Overall heats the Earth = {:.1f} keV/decay = {:.5f} pJ/decay".format(Hf182.heat, Hf182.heat * pJ_in_keV))
129I → 129Xe¶
β- decay with half-life 15.7 ± 0.4 Myr (NuDat referencing Timar et al 2014 NDS)
Mougeot 2015: second forbidden non-unique transition, experimental shape factor from der Mateosian & Wu 1954 Phys.Rev.
In [57]:
mass_parent = 128.9049837 # u [NIST]
mass_daughter = 128.9047808611 # u [NIST]
Z = 54 # atomic number of daughter
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
print ('Q = {:.1f} keV (NuDat: 189 ± 3 keV)\n'.format(Q))
# NuDat inputs
branch = 1.
#Q = 189
Eexc = 39.578
el_mean = 40.03
heat_nudat = (el_mean + Eexc) * branch
print ("Using NuDat inputs:")
print (" Mean electron energy = {:.1f} keV\n".format(el_mean))
# calculating spectrum
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
Eexc = 39.578
Eendpoint = Q - Eexc
q = Emesh(Eendpoint,dE)
E = Eendpoint - q
w = me + E
p = np.sqrt(w**2-me**2)
print ("My calculation:")
shape = 1
dNdE = normalize_dydx(p*w*q**2 * fermi_func(w,Z) * shape, dE, branch)
anu_mean = mean_x(dNdE, q, dE)
el_mean = Eendpoint - anu_mean
heat_me1 = (el_mean + Eexc) * branch
print (" Electron endpoint energy = {:.2f} keV".format(Eendpoint))
print (" Mean electron energy w/ shape factor equal to 1 = {:.2f} keV".format(el_mean))
shape = 3.16*q**2 + p**2 # Mougeot 2015 <- der Mateosian & Wu 1954 Phys.Rev.
dNdE = normalize_dydx(p*w*q**2 * fermi_func(w,Z) * shape, dE, branch)
anu_mean = mean_x(dNdE, q, dE)
el_mean = Eendpoint - anu_mean
heat_me_shape = (el_mean + Eexc) * branch
print (" Mean electron energy w/ exp. shape factor from Mougeot = {:.2f} keV (Mougeot: 45.98 keV)".format(el_mean))
print (" Mean electron energy w/ calc. by Mougeot using ξ approx. = 48.62 keV\n".format(el_mean))
print ("I will use my (~Mougeot) calculation with experimental shape factor:\n")
I129 = Nuclide("129I")
I129.halflife = 1.57e7
I129.mass_atomic = 126.90447
I129.isotfrac = 1.4e-4 * 1. # 129I/127I * 127I/I
I129.timeref = tcai
I129.heat = heat_me_shape
I129.lanu = 1
I129.lnu = 0
#I129.heat_gamma = 39.578 * 0.0751
gam_lev = np.array ([ 4.11, 29.461, 29.782, 33.562, 33.624, 34.419, 39.578 ])
gam_int = np.array ([ 7.8, 19.5, 35.9, 3.36, 6.5, 1.97, 7.51 ]) / 100
I129.heat_gamma = sum(gam_lev*gam_int)
print ("Heats the Earth = {:.1f} keV/decay = {:.5f} pJ/decay".format(I129.heat, I129.heat * pJ_in_keV))
205Pb → 205Tl¶
ε decay (EC only; no β+) with half-life 17.3 ± 0.7 Myr, to ground state [NuDat references Kondev 2004 NDS]
In [58]:
mass_parent = 204.9744822 # u [NIST]
mass_daughter = 204.9744278 # u [NIST]
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
print ('Q = {:.1f} keV (NuDat: 50.5 ± 0.5 keV)\n'.format(Q))
Pb205 = Nuclide("205Pb")
Pb205.halflife = 1.73e7
Pb205.mass_atomic = 207.2
Pb205.isotfrac = 1e-3 * 0.014 # 205Pb/204Pb * 204Pb/Pb
Pb205.timeref = tcai
Pb205.heat = 0.
Pb205.lanu = 0
Pb205.lnu = 1
#Pb205.heat_gamma = 0.
Pb205.heat_gamma = 10.3 * 0.220
print ("Heats the Earth = {:.1f} keV/decay = {:.5f} pJ/decay".format(Pb205.heat, Pb205.heat * pJ_in_keV))
92Nb → 92Zr¶
ε decay (EC only; no β+) with half-life 34.7 ± 2.4 Myr, to 1495.6 keV level of 92Zr [NuDat references Baglin 2012 NDS]
In [59]:
mass_parent = 91.9071881 # u [NIST]
mass_daughter = 91.9050347 # u [NIST]
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
print ('Q = {:.1f} keV (NuDat: 2005.9 ± 1.8 keV)\n'.format(Q))
branch = 1.
Eexc = 1495.6
Nb92 = Nuclide("92Nb")
Nb92.halflife = 3.47e7
Nb92.mass_atomic = 92.90637 # Hf atomic mass
Nb92.isotfrac = 1.7e-5 * 1. # 92Nb/93Nb * 93Nb/Nb
Nb92.timeref = tcai
Nb92.heat = Eexc * branch
Nb92.lanu = 0
Nb92.lnu = 1
#Nb92.heat_gamma = 561.1 * 1.00 + 934.5 * 0.74
gam_lev = np.array ([ 2.04, 15.691, 15.775, 17.653, 17.667, 17.969, 561.1, 934.5 ])
gam_int = np.array ([ 3.25, 17.9, 34.3, 2.65, 5.15, 1.05, 100, 74 ]) / 100
Nb92.heat_gamma = sum(gam_lev*gam_int)
print ("Heats the Earth = {:.1f} keV/decay = {:.5f} pJ/decay".format(Nb92.heat, Nb92.heat * pJ_in_keV))
244Pu → 240U → 240Np → 240Pu → 236U → 232Th¶
α → β- → β- → α → α
In [60]:
mass_parent = 244.0642053 # u [NIST]
mass_daughter = 240.0565934 # u [NIST]
Q = (mass_parent - mass_daughter - mass_4He_u) * keV_per_u # keV Q-value
Qtot = Q
print ('Q = {:.1f} keV (NuDat: 4665.5 ± 1.0 keV)\n'.format(Q))
Pu244_heat1 = Q
#Pu244_heat_gamma1 = 44 * 0.00029
Pu244_heat_gamma1 = 44 * 0.00029 + 44 * 0.00029
print ("Heats the Earth = {:.1f} keV/decay = {:.4f} pJ/decay".format(Pu244_heat1, Pu244_heat1 * pJ_in_keV))
In [61]:
mass_parent = 240.0565934 # u [NIST]
mass_daughter = 240.056165 # u [NIST]
Z = 93 # atomic number of daughter
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
Qtot = Qtot + Q
print ('Q = {:.1f} keV (NuDat: 400 ± 16 keV)\n'.format(Q))
# NuDat inputs
el_mean_lev = np.array([ 20.6, 22.2, 25.9, 50.9, 73.6, 73.9, 82.3, 94.1, 107.8 ])
el_endp_lev = np.array([ 100, 105, 120, 210, 288, 289, 317, 356, 400 ])
branch_lev = np.array([ 0.015, 0.005, 0.5, 1.8, 0.09, 1.3, 0.07, 25.0, 75 ]) / 100
print ("Sum of branching ratios = {:.4f}".format(sum(branch_lev)))
print (" renormalizing to 1.000...")
norm = sum(branch_lev)
branch_lev /= norm
el_mean = sum(el_mean_lev*branch_lev)
print ("Mean electron energy = {:.1f} keV (NuDat: 103 ± 13 keV)".format(el_mean))
deexcite = sum((Q - el_endp_lev) * branch_lev)
print ("Deexcitation energy per decay = {:.1f} keV".format(deexcite))
heat = el_mean + deexcite
Pu244_heat2 = heat
#gam_lev = np.array([ 44.10, 49.1, 50.3, 66.5, 78.1, 82.6, 128.3, 145.4, 169.2, 189.7, 212.3, 255.6, 280.1, 294.8, 299.8 ])
#gam_int = np.array([ 1.05, 0.0070, 0.0050, 0.154, 0.0040, 0.0140, 0.0870, 0.0810, 0.115, 0.240, 0.0015, 0.0040, 0.0160, 0.0019, 0.0130 ]) / 100
#Pu244_heat_gamma2 = sum(gam_lev*gam_int)
gam_lev = np.array([ 13.9, 44.10, 49.1, 50.3, 66.5, 78.1, 82.6, 128.3, 145.4, 169.2, 189.7, 212.3, 255.6, 280.1, 294.8, 299.8 ])
gam_int = np.array([ 24, 1.05, 0.0070, 0.0050, 0.154, 0.0040, 0.0140, 0.0870, 0.0810, 0.115, 0.240, 0.0015, 0.0040, 0.0160, 0.0019, 0.0130 ]) / 100
Pu244_heat_gamma2 = sum(gam_lev*gam_int)
print ("Heats the Earth = {:.1f} keV/decay = {:.5f} pJ/decay".format(Pu244_heat2, Pu244_heat2 * pJ_in_keV))
In [62]:
mass_parent = 240.056165 # u [NIST]
mass_daughter = 240.0538138 # u [NIST]
Z = 94 # atomic number of daughter
Q = (mass_parent - mass_daughter) * keV_per_u # keV Q-value
Qtot = Qtot + Q
print ('Q = {:.1f} keV (NuDat: 2188 ± 15 keV)\n'.format(Q))
# NuDat inputs
branch = 0.99880
el_mean_lev = np.array([ 23.4, 57.1, 69.3, 80.1, 113.9, 117.6, 124.2, 145.1, 170.5, 179.1, 195.6, 202.2, 206.9, 220.0, 237.3, 247.0, 278.9, 308.0, 314.5, 346.1, 348.3, 363.8, 412.8, 420.7, 435.1, 450.2, 505.8, 552.2, 733.2, 772.9, 790.2 ])
el_endp_lev = np.array([ 70, 192, 233, 270, 380, 392, 413, 478, 555, 580, 629, 648, 662, 700, 750, 777, 867, 947, 965, 1051, 1057, 1099, 1229, 1250, 1288, 1327, 1539, 1591, 2046, 2145, 2188 ])
branch_lev = np.array([ 0.0038, 0.0080, 0.055, 0.0130, 0.030, 0.041, 0.0050, 0.084, 0.260, 0.070, 0.31, 2.30, 0.160, 0.579, 0.36, 0.20, 0.034, 0.0100, 0.028, 0.014, 0.180, 0.100, 1.26, 1.40, 3.80, 2.50, 0.05, 31.0, 0.7, 42, 10.0 ]) / 100
print ("Sum of branching ratios = {:.4f}".format(sum(branch_lev)))
print (" renormalizing to branching {:.4f} (NuDat)".format(branch))
norm = sum(branch_lev)
branch_lev *= branch/norm
print ("Sum of branching ratios = {:.4f}".format(sum(branch_lev)))
el_mean = sum(el_mean_lev*branch_lev)
print ("Mean electron energy = {:.1f} keV (NuDat: 650 ± 50 keV)".format(el_mean))
deexcite = sum((Q - el_endp_lev) * branch_lev)
print ("Deexcitation energy per decay = {:.1f} keV".format(deexcite))
heat = el_mean + deexcite
Pu244_heat3 = heat
#gam_lev = np.array([ 42.824, 98.860, 152.630, 251.47, 263.37, 289.21, 302.98, 309.99, 340.70, 361.55, 475.0, 496.7, 507.2, 518.2, 554.60, 573.4, 573.4, 580.7, 597.40, 606.10, 658.5, 758.61, 789.59, 796.2, 813.41, 817.2, 817.89, 837.6, 841.11, 857.48, 890.6, 895.3, 900.37, 910.1, 915.98, 928.55, 938.02, 942.39, 959.0, 961.62, 989.2, 1036.5, 1046.62, 1061.6, 1088.3, 1094.2, 1113.2, 1131.0, 1137.0, 1159.2, 1180.1, 1198.0, 1210.5, 1223.0, 1305.8, 1321.1, 1357.2, 1398.5, 1417.2, 1438.5, 1445.3, 1483.0, 1488.2, 1496.9, 1515.9, 1539.62, 1558.8, 1568.6, 1584.1, 1590.5, 1607.6, 1626.6, 1633.33, 1667.6, 1711.0, 1732.4, 1752.9, 1765.2, 1775.3, 1796.2, 1807.9, 1812.8, 1874.9, 1911.4, 1918.0, 1953.6, 1996.7, 2074.8, 2117.5 ])
#gam_int = np.array([ 0.074, 0.17, 0.010, 0.86, 1.139, 0.017, 1.00, 0.044, 0.060, 0.036, 0.011, 0.0100, 0.70, 0.0060, 20.9, 0.0080, 0.0080, 0.0070, 11.7, 0.67, 0.009, 1.18, 0.18, 5E-4, 0.18, 0.05, 1.28, 0.008, 0.150, 0.489, 0.0170, 0.061, 0.160, 0.140, 1.02, 0.150, 1.21, 0.096, 0.0073, 0.130, 0.081, 0.0030, 0.097, 0.029, 0.032, 0.020, 0.018, 0.062, 0.0140, 0.0060, 0.020, 0.0090, 0.015, 0.018, 0.023, 0.034, 0.013, 0.0050, 0.023, 5E-4, 0.380, 0.027, 0.200, 1.33, 0.015, 0.839, 0.0060, 0.0060, 0.0170, 0.097, 0.055, 0.0050, 0.154, 0.019, 0.0020, 0.0020, 0.0054, 0.0070, 0.0030, 0.0030, 0.0020, 0.0050, 0.0120, 0.0140, 8E-4, 0.0023, 1.0E-3, 0.0031, 7E-4 ]) / 100
#Pu244_heat_gamma3 = sum(gam_lev*gam_int)
gam_lev = np.array([ 14.3, 42.824, 98.860, 99.525, 103.374, 116.244, 117.228, 120.54, 152.630, 251.47, 263.37, 289.21, 302.98, 309.99, 340.70, 361.55, 475.0, 496.7, 507.2, 518.2, 554.60, 573.4, 573.4, 580.7, 597.40, 606.10, 658.5, 758.61, 789.59, 796.2, 813.41, 817.2, 817.89, 837.6, 841.11, 857.48, 890.6, 895.3, 900.37, 910.1, 915.98, 928.55, 938.02, 942.39, 959.0, 961.62, 989.2, 1036.5, 1046.62, 1061.6, 1088.3, 1094.2, 1113.2, 1131.0, 1137.0, 1159.2, 1180.1, 1198.0, 1210.5, 1223.0, 1305.8, 1321.1, 1357.2, 1398.5, 1417.2, 1438.5, 1445.3, 1483.0, 1488.2, 1496.9, 1515.9, 1539.62, 1558.8, 1568.6, 1584.1, 1590.5, 1607.6, 1626.6, 1633.33, 1667.6, 1711.0, 1732.4, 1752.9, 1765.2, 1775.3, 1796.2, 1807.9, 1812.8, 1874.9, 1911.4, 1918.0, 1953.6, 1996.7, 2074.8, 2117.5 ])
gam_int = np.array([ 27, 0.074, 0.17, 0.194, 0.308, 0.0372, 0.0731, 0.0286, 0.010, 0.86, 1.139, 0.017, 1.00, 0.044, 0.060, 0.036, 0.011, 0.0100, 0.70, 0.0060, 20.9, 0.0080, 0.0080, 0.0070, 11.7, 0.67, 0.009, 1.18, 0.18, 5E-4, 0.18, 0.05, 1.28, 0.008, 0.150, 0.489, 0.0170, 0.061, 0.160, 0.140, 1.02, 0.150, 1.21, 0.096, 0.0073, 0.130, 0.081, 0.0030, 0.097, 0.029, 0.032, 0.020, 0.018, 0.062, 0.0140, 0.0060, 0.020, 0.0090, 0.015, 0.018, 0.023, 0.034, 0.013, 0.0050, 0.023, 5E-4, 0.380, 0.027, 0.200, 1.33, 0.015, 0.839, 0.0060, 0.0060, 0.0170, 0.097, 0.055, 0.0050, 0.154, 0.019, 0.0020, 0.0020, 0.0054, 0.0070, 0.0030, 0.0030, 0.0020, 0.0050, 0.0120, 0.0140, 8E-4, 0.0023, 1.0E-3, 0.0031, 7E-4 ]) / 100
Pu244_heat_gamma3 = sum(gam_lev*gam_int)
print ("Heats the Earth = {:.1f} keV/decay = {:.5f} pJ/decay".format(Pu244_heat3, Pu244_heat3 * pJ_in_keV))
In [63]:
mass_parent = 240.0538138 # u [NIST]
mass_daughter = 236.0455682 # u [NIST]
Q = (mass_parent - mass_daughter - mass_4He_u) * keV_per_u # keV Q-value
Qtot = Qtot + Q
print ('Q = {:.2f} keV (NuDat: 5255.75 ± 1.4 keV)\n'.format(Q))
Pu244_heat4 = Q
#gam_lev = np.array([ 45.244, 104.234, 160.308, 212.46, 538.09, 642.35, 687.57, 699, 873.92, 958, 960, 967 ])
#gam_int = np.array([ 0.0447, 0.00714, 4.02E-4, 2.9E-5, 1.47E-7, 1.30E-5, 3.50E-6, 1.3E-8, 5.8E-7, 5E-8, 3E-8, 3E-8 ]) / 100
#Pu244_heat_gamma4 = sum(gam_lev*gam_int)
gam_lev = np.array([ 13.6, 45.244, 94.654, 98.434, 104.234, 110.421, 111.298, 114.445, 160.308, 212.46, 538.09, 642.35, 687.57, 699, 873.92, 958, 960, 967 ])
gam_int = np.array([ 9.6, 0.0447, 2.54E-5, 4.05E-5, 0.00714, 5.08E-6, 9.6E-6, 3.73E-6, 4.02E-4, 2.9E-5, 1.47E-7, 1.30E-5, 3.50E-6, 1.3E-8, 5.8E-7, 5E-8, 3E-8, 3E-8 ]) / 100
Pu244_heat_gamma4 = sum(gam_lev*gam_int)
print ("Heats the Earth = {:.1f} keV/decay = {:.4f} pJ/decay".format(Pu244_heat4, Pu244_heat4 * pJ_in_keV))
In [64]:
mass_parent = 236.0455682 # u [NIST]
mass_daughter = 232.0380558 # u [NIST]
Q = (mass_parent - mass_daughter - mass_4He_u) * keV_per_u # keV Q-value
Qtot = Qtot + Q
print ('Q = {:.2f} keV (NuDat: 4573.1 ± 0.9 keV)\n'.format(Q))
Pu244_heat5 = Q
#gam_lev = np.array([ 49.46, 112.79, 171.15 ])
#gam_int = np.array([ 0.078, 0.019, 6.2E-5 ]) / 100
#Pu244_heat_gamma5 = sum(gam_lev*gam_int)
gam_lev = np.array([ 13.0, 49.46, 89.957, 93.35, 104.819, 105.604, 108.582, 112.79, 171.15 ])
gam_int = np.array([ 9.0, 0.078, 0.00124, 0.0020, 2.5E-4, 4.7E-4, 1.8E-4, 0.019, 6.2E-5 ]) / 100
Pu244_heat_gamma5 = sum(gam_lev*gam_int)
print ("Heats the Earth = {:.1f} keV/decay = {:.4f} pJ/decay".format(Pu244_heat5, Pu244_heat5 * pJ_in_keV))
244Pu → 240U → 240Np → 240Pu → 236U → 232Th decay summary¶
In [65]:
Pu244 = Nuclide("244Pu")
Pu244.halflife = 8.11e7
Pu244.mass_atomic = U238.mass_atomic # U atomic mass
Pu244.isotfrac = 8e-3 * U238.isotfrac # 244Pu/238U * 238U/U
Pu244.timeref = tcai
Pu244.heat = Pu244_heat1 + Pu244_heat2 + Pu244_heat3 + Pu244_heat4 + Pu244_heat5
Pu244.lanu = 2
Pu244.lnu = 0
Pu244.heat_gamma = Pu244_heat_gamma1 + Pu244_heat_gamma2 + Pu244_heat_gamma3 + Pu244_heat_gamma4 + Pu244_heat_gamma5
print ("244Pu → 240U rad. heat = {:6.1f} keV/decay = {:.5f} pJ/decay".format(Pu244_heat1, Pu244_heat1 * pJ_in_keV))
print ("240U → 240Np rad. heat = {:6.1f} keV/decay = {:.5f} pJ/decay".format(Pu244_heat2, Pu244_heat2 * pJ_in_keV))
print ("240Np → 240Pu rad. heat = {:6.1f} keV/decay = {:.5f} pJ/decay".format(Pu244_heat3, Pu244_heat3 * pJ_in_keV))
print ("240Pu → 236U rad. heat = {:6.1f} keV/decay = {:.5f} pJ/decay".format(Pu244_heat4, Pu244_heat4 * pJ_in_keV))
print ("236U → 232Th rad. heat = {:6.1f} keV/decay = {:.5f} pJ/decay\n".format(Pu244_heat5, Pu244_heat5 * pJ_in_keV))
print ("Overall heats the Earth = {:.1f} keV/decay = {:.5f} pJ/decay".format(Pu244.heat, Pu244.heat * pJ_in_keV))
print ("Overall Q = {:.1f} keV/decay = {:.5f} pJ/decay".format(Qtot, Qtot * pJ_in_keV))
Fill nuclide arrays¶
In [66]:
# lists of nuclides
Nuclides_long = [K40, Rb87, La138, Sm147, Lu176, Re187, Pt190, Th232, U235, U238]
Nuclides_short = [Ca41, Tc99, Sn126, Cl36, Al26, Be10, Cs135, Fe60, Mn53,
Tc98, Tc97, Pd107, Hf182, I129, Pb205, Nb92, Pu244, Sm146]
Nuclides_special = [Al26, Fe60]
Nuclides_short2 = [Ca41, Tc99, Sn126, Cl36, Be10, Cs135, Mn53, Tc98, Tc97,
Pd107, Hf182, I129, Pb205, Nb92, Pu244, Sm146] # short without Al26, Fe60
Nuclides_all = Nuclides_long + Nuclides_short # all nuclides considered
Nuclides_noabundance = [Kr81, Se79, Zr93, Gd150, Dy154]
Overall summary of heat, (anti)neutrinos per decay¶
In [67]:
def print_heatanunu(nucl):
if nucl.heat<0:
print ("{:s} ...initial concentration NOT AVAILABLE...".format(nucl.name))
else:
if nucl.heat>0:
print ("{:5s} {:10.1f}{:10.4f}{:10.4f}{:10.4f}{:10.4f}".format(nucl.name, nucl.heat, nucl.heat*pJ_in_keV, nucl.lanu, nucl.lnu, nucl.heat_gamma/nucl.heat))
else:
print ("{:5s} {:10.1f}{:10.4f}{:10.4f}{:10.4f} -".format(nucl.name, nucl.heat, nucl.heat*pJ_in_keV, nucl.lanu, nucl.lnu))
return;
print ("Radiogenic heat (keV, pJ), antineutrinos, neutrinos per decay")
print (" keV pJ anu nu")
print ("Long-lived")
for N in Nuclides_long:
print_heatanunu(N)
print ("Short-lived")
for N in Nuclides_short:
print_heatanunu(N)
print ("Short-lived with unknown initial abundance")
for N in Nuclides_noabundance:
print_heatanunu(N)
Rock composition inputs¶
Composition of Bulk Earth taken from McDonough & Sun (1995) and McDonough (2014, Treatise Geochem 3.16)
In [68]:
### BULK EARTH
# McDonough & Sun 1995 [https://doi.org/10.1016/0009-2541(94)00140-4] and
# McDonough 2014 Treatise Geochem 3.16 [https://doi.org/10.1016/B978-0-08-095975-7.00215-1],
# Th/K from Wipperfurth et al 2018 [https://doi.org/10.1016/j.epsl.2018.06.029]
## Long-lived
K40.massfrac_el_BE = 280e-6 * Mbse / Mearth # no K in core, Arevalo & McDonough 2009
Rb87.massfrac_el_BE = Mbse * 0.6e-6 / Mearth # all Rb in BSE
La138.massfrac_el_BE = Mbse * 0.65e-6 / Mearth # all La in BSE
Sm147.massfrac_el_BE = Mbse * 406e-9 / Mearth # all Sm in BSE
Lu176.massfrac_el_BE = Mbse * 68e-9 / Mearth # all Lu in BSE
Re187.massfrac_el_BE = 75e-9 # and Re mostly in the core
Pt190.massfrac_el_BE = 1.9e-6 # and Pt mostly in the core
Th232.massfrac_el_BE = Mbse * 20e-9 * 3.77 / Mearth # all Th in BSE
U235.massfrac_el_BE = Mbse * 20e-9 / Mearth # all U in BSE
## Short-lived
Ca41.massfrac_el_BE = Mbse * 2.53e-2 / Mearth # all Ca in BSE
Tc99.massfrac_el_BE = (5e-9*Mbse + 4e-6*Mcore) / Mearth # Ru mostly in the core
# 81Kr ... initial abundance unknown
Sn126.massfrac_el_BE = (130e-9*Mbse + 500e-9*Mcore) / Mearth # Sn in BSE and core
Cl36.massfrac_el_BE = (17e-6*Mbse + 200e-6*Mcore) / Mearth # Cl in BSE and core
# 79Se ... initial abundance unknown
Al26.massfrac_el_BE = Mbse * 2.35e-2 / Mearth # all Al in BSE
Be10.massfrac_el_BE = Mbse * 68e-9 / Mearth # all Be in BSE
# 93Zr ... initial abundance unknown
# 150Gd ... initial abundance unknown
Cs135.massfrac_el_BE = (Mbse * 21e-9 + Mcore * 65e-9) / Mearth # Cs in BSE and core
Fe60.massfrac_el_BE = (Mbse * 0.0626 + Mcore * 0.85) / Mearth # Fe in BSE and core
# 154Dy ... initial abundance unknown
Mn53.massfrac_el_BE = (Mbse * 1.045e-3 + Mcore * 0.3e-3) / Mearth # Mn in BSE and core
Pd107.massfrac_el_BE = (Mbse * 3.9e-9 + Mcore * 3.1e-6) / Mearth # Pd mostly in the core
Hf182.massfrac_el_BE = Mbse * 283e-9 / Mearth # all Hf in BSE
I129.massfrac_el_BE = (Mbse * 10e-9 + Mcore * 130e-9) / Mearth # I mostly in the core
Pb205.massfrac_el_BE = (Mbse * 0.15e-6 + Mcore * 0.4e-6) / Mearth # Pb mostly BSE and core
Nb92.massfrac_el_BE = Mbse * 658e-9 / Mearth # all Nb in BSE
U238.massfrac_el_BE = -1
Tc98.massfrac_el_BE = -1
Tc97.massfrac_el_BE = -1
Pu244.massfrac_el_BE = -1
Sm146.massfrac_el_BE = -1
### ARCHEAN TTG CRUST
### Johnson et al 2019 EPSL doi:10.1016/j.epsl.2018.10.022
K40.massfrac_el_TTG = 1.34E-02
Rb87.massfrac_el_TTG = 4.43E-05
La138.massfrac_el_TTG = 2.18E-05
Sm147.massfrac_el_TTG = 2.80E-06
Lu176.massfrac_el_TTG = 1.18E-07
Re187.massfrac_el_TTG = 0. # lack of data -> assumption (as siderophile)
Pt190.massfrac_el_TTG = 0. # lack of data -> assumption (as siderophile)
Th232.massfrac_el_TTG = 4.81E-06
U235.massfrac_el_TTG = 9.34E-07
Ca41.massfrac_el_TTG = 2.04E-02
Tc99.massfrac_el_TTG = 3.40E-10 # Ru
Sn126.massfrac_el_TTG = 2.10E-06
Cl36.massfrac_el_TTG = 3.70E-04
Al26.massfrac_el_TTG = 8.33E-02
Be10.massfrac_el_TTG = 2.10E-06
Cs135.massfrac_el_TTG = 1.77E-06
Fe60.massfrac_el_TTG = 1.92E-02
Mn53.massfrac_el_TTG = 3.27E-04
Pd107.massfrac_el_TTG = 5.20E-10
Hf182.massfrac_el_TTG = 4.20E-06
I129.massfrac_el_TTG = 1.40E-06
Pb205.massfrac_el_TTG = 1.16E-05
Nb92.massfrac_el_TTG = 4.45E-06
U238.massfrac_el_TTG = -1
Tc98.massfrac_el_TTG = -1
Tc97.massfrac_el_TTG = -1
Pu244.massfrac_el_TTG = -1
Sm146.massfrac_el_TTG = -1
for N in Nuclides_all:
aBE = N.massfrac_el_BE
aTTG = N.massfrac_el_TTG
if N.massfrac_el_BE > 0: print("{:5s} {:.2e} {:.2e} {:f}".format(N.name, aBE, aTTG, aTTG/aBE))
In [69]:
select_composition = 1 # 1=Bulk Earth, 2=TTG crust
if select_composition==1: # Bulk Earth
for N in Nuclides_all:
N.massfrac_el = N.massfrac_el_BE
elif select_composition==2: # TTG crust
for N in Nuclides_all:
N.massfrac_el = N.massfrac_el_TTG
U238.massfrac_el = U235.massfrac_el # U
Tc98.massfrac_el = Tc99.massfrac_el # Ru
Tc97.massfrac_el = Tc99.massfrac_el # Ru
Pu244.massfrac_el = U235.massfrac_el # U
Sm146.massfrac_el = Sm147.massfrac_el # Sm
In [70]:
print ("Al mass fraction (present-day) =", Al26.massfrac_el)
print ("26Al initial W/kg =", Al26.Heat_per_kgrock())
print ("\nFe mass fraction =", Fe60.massfrac_el)
print ("60Fe initial W/kg =", Fe60.Heat_per_kgrock())
Time series¶
In [71]:
# time array
tstart = 1e5 # in yr
tend = tnow
time = np.exp(np.linspace(log(tstart), log(tend), num=100))
Activity¶
In [72]:
# account for decay steps in multi-step decays
K40.dsteps = 1
Rb87.dsteps = 1
La138.dsteps = 1
Sm147.dsteps = 1
Lu176.dsteps = 1
Re187.dsteps = 1
Pt190.dsteps = 1
Th232.dsteps = 10
U235.dsteps = 11
U238.dsteps = 14
Ca41.dsteps = 1
Tc99.dsteps = 1
Sn126.dsteps = 2
Cl36.dsteps = 1
Al26.dsteps = 1
Be10.dsteps = 1
Cs135.dsteps = 1
Fe60.dsteps = 2
Mn53.dsteps = 1
Tc98.dsteps = 1
Tc97.dsteps = 1
Pd107.dsteps = 1
Hf182.dsteps = 2
I129.dsteps = 1
Pb205.dsteps = 1
Nb92.dsteps = 1
Pu244.dsteps = 5
Sm146.dsteps = 1
In [73]:
total_activity_becq = np.zeros(np.size(time))
total_activity_becq0 = 0.
for N in Nuclides_all:
print(N.name)
total_activity_becq += N.activity_becq(time) * N.dsteps
total_activity_becq0 += N.activity_becq(tcai) * N.dsteps
def figure_activity():
facx = 1e-6
facy = 1
fig, ax = plt.subplots(figsize=(12,8))
#
for N in Nuclides_long:
ax.loglog(time * facx, N.activity_becq(time) * N.dsteps * facy, '--', label=N.name)
for N in Nuclides_special:
ax.loglog(time * facx, N.activity_becq(time) * N.dsteps * facy, '-.', label=N.name)
for N in Nuclides_short2:
ax.loglog(time * facx, N.activity_becq(time) * N.dsteps * facy, label=N.name)
ax.loglog(time * facx, total_activity_becq * facy, 'r-', linewidth=3, label='TOTAL')
#
ax.set_xlim(time[0]*facx,time[-1]*facx)
ax.set_xlabel('Time (million years)')
ax.set_ylabel('Activity (Bq/kg)')
ax.set_title('Activity (Bq/kg) from natural radionuclides')
ax.yaxis.set_ticks_position('both')
ax.tick_params(labeltop=False, labelright=True)
ax.legend(bbox_to_anchor=(1.06, 1), loc=2, borderaxespad=0.)
ax.set_position([0.1,0.1,0.75,0.8])
if select_composition==1:
ax.set_ylim(1e-3,1e6)
fig.savefig('Fig_BE_activity_nuclides.pdf')
elif select_composition==2:
ax.set_ylim(1e-2,1e7)
fig.savefig('Fig_TTG_activity_nuclides.pdf')
figure_activity()
In [74]:
# write ASCII file
if select_composition==1:
f = open('Tab_BE_activity_nuclides.txt', 'w')
elif select_composition==2:
f = open('Tab_TTG_activity_nuclides.txt', 'w')
# header
f.write(' {:>10}'.format('years'))
for N in Nuclides_all:
f.write(' {:>10}'.format(N.name))
f.write('\n')
# data
for i in np.arange(time.size):
f.write(' {:10.3e}'.format(time[i]))
for N in Nuclides_all:
f.write(' {:10.3e}'.format(N.activity_becq(time[i]) * N.dsteps))
f.write('\n')
f.close()
In [75]:
def fig_activity_total(case):
facx = 1e-6
facy = 1
fig, ax = plt.subplots(figsize=(12,8))
#
ax.loglog(time * facx, total_activity_becq * facy)
#
xbars = tend - np.array([4.2, 3.8, 3.4])*1e9
ymin=0
if select_composition==1:
ymax=200 if case==0 else 80 if case==1 else 0
elif select_composition==2:
ymax=1e4 if case==0 else 7e3 if case==1 else 0
ax.vlines(xbars[0]*facx, ymin, ymax, colors='k', linestyles=':')
ax.vlines(xbars[1]*facx, ymin, ymax, colors='k', linestyles=':')
ax.vlines(xbars[2]*facx, ymin, ymax, colors='k', linestyles=':')
ax.text(xbars[0]*facx, ymax, '4.2 Ga', horizontalalignment='center')
ax.text(xbars[1]*facx, ymax, '3.8 Ga', horizontalalignment='center')
ax.text(xbars[2]*facx, ymax, '3.4 Ga', horizontalalignment='left')
#
ax.set_xlabel('Time (million years)')
ax.set_ylabel('Activity (Bq/kg)')
ax.set_title('Total activity (Bq/kg) from natural radionuclides')
ax.yaxis.set_ticks_position('both')
ax.tick_params(labeltop=False, labelright=True)
#ax.legend(bbox_to_anchor=(1.06, 1), loc=2, borderaxespad=0.)
ax.set_position([0.1,0.1,0.75,0.8])
if case==0:
ax.set_xlim(time[0]*facx, time[-1]*facx)
if select_composition==1:
fig.savefig('Fig_BE_activity_total.pdf')
elif select_composition==2:
fig.savefig('Fig_TTG_activity_total.pdf')
elif case==1:
ax.set_xlim(1e8*facx,time[-1]*facx)
if select_composition==1:
ax.set_ylim(1e1,1e2)
fig.savefig('Fig_BE_activity_total_zoom.pdf')
elif select_composition==2:
ax.set_ylim(6e2,1e4)
fig.savefig('Fig_TTG_activity_total_zoom.pdf')
return
fig_activity_total(0)
fig_activity_total(1)
In [76]:
fact = 1e-0
print("TOTAL ACTIVITY")
print("time(yr) Activity(Bq/kg)")
print("{:.1e} {:.1e}".format(tcai*fact, total_activity_becq0))
print("{:.1e} {:.1e}".format(time[0]*fact, total_activity_becq[0]))
print("{:.1e} {:.1e}".format(time[22]*fact, total_activity_becq[22]))
print("{:.1e} {:.1e}".format(time[43]*fact, total_activity_becq[43]))
print("{:.1e} {:.1e}".format(time[64]*fact, total_activity_becq[64]))
print("{:.1e} {:.1e}".format(time[79]*fact, total_activity_becq[79]))
print("{:.1e} {:.1e}".format(time[85]*fact, total_activity_becq[85]))
print("{:.1e} {:.1e}".format(time[-1]*fact, total_activity_becq[-1]))
Power of radiogenic heat sources¶
In [77]:
total_watts = np.zeros(np.size(time))
for N in Nuclides_all:
total_watts += N.activity_becq(time) * N.heatJoule() * Mearth
def figure_power(case = 0):
facx = 1e-6
facy = 1e-12 * Mearth
fig, ax = plt.subplots(figsize=(12,8))
#
for N in Nuclides_long:
if N.heat > 0: ax.loglog(time * facx, N.activity_becq(time) * N.heatJoule() * facy, '--', label=N.name)
for N in Nuclides_special:
if N.heat > 0: ax.loglog(time * facx, N.activity_becq(time) * N.heatJoule() * facy, '-.', label=N.name)
for N in Nuclides_short2:
if N.heat > 0: ax.loglog(time * facx, N.activity_becq(time) * N.heatJoule() * facy, label=N.name)
ax.loglog(time * facx, total_watts * facy/Mearth, 'r-', linewidth=3, label='TOTAL')
#
ax.set_ylim(0.001,2000000)
ax.set_xlabel('Time (million years)')
ax.set_ylabel('Power (TW)')
ax.set_title('Radiogenic power (TW) from natural radionuclides')
ax.yaxis.set_ticks_position('both')
ax.tick_params(labeltop=False, labelright=True)
ax.legend(bbox_to_anchor=(1.06, 1), loc=2, borderaxespad=0.)
ax.set_position([0.1,0.1,0.75,0.8])
if case==0:
ax.set_xlim(time[0]*facx,time[-1]*facx)
if select_composition==1: fig.savefig('Fig_BE_radiogenic_power_nuclides.pdf')
elif case==1:
ax.set_xlim(0.1,100)
if select_composition==1: fig.savefig('Fig_BE_radiogenic_power_nuclides_zoom.pdf')
return
if select_composition==1: figure_power(0)
#if select_composition==1: figure_power(1)
In [78]:
def figure_power_select():
facx = 1e-6
facy = 1e-12 * Mearth
fig, ax = plt.subplots(figsize=(12,8))
select_watts = np.zeros(np.size(time))
#
for N in [Al26, Fe60, Be10, Pu244, Sm146]:
ax.loglog(time * facx, N.activity_becq(time) * N.heatJoule() * facy, label=N.name)
select_watts += N.activity_becq(time) * N.heatJoule() * Mearth
#
for N in [K40, Th232, U238, U235]:
ax.loglog(time * facx, N.activity_becq(time) * N.heatJoule() * facy, '--', label=N.name)
select_watts += N.activity_becq(time) * N.heatJoule() * Mearth
#
rest_watts = total_watts - select_watts
ax.loglog(time * facx, rest_watts * facy/Mearth, ':', label='All other')
#
ax.loglog(time * facx, total_watts * facy/Mearth, 'r-', linewidth=3, label='TOTAL')
#
#ax.set_xlim(time[0]*facx,time[-1]*facx)
ax.set_xlim(1,time[-1]*facx)
ax.set_ylim(0.1,1000000)
ax.set_xlabel('Time (million years)')
ax.set_ylabel('Power (TW)')
ax.set_title('Radiogenic power (TW) from natural radionuclides')
ax.yaxis.set_ticks_position('both')
ax.tick_params(labeltop=False, labelright=True)
ax.legend(bbox_to_anchor=(1.06, 1), loc=2, borderaxespad=0.)
ax.set_position([0.1,0.1,0.75,0.8])
if select_composition==1: fig.savefig('Fig_BE_radiogenic_power_select.pdf')
if select_composition==1: figure_power_select()
Specific power of radiogenic heat sources¶
In [79]:
def figure_specpower(case = 0):
facx = 1e-6
facy = 1 #1e-12 * Mearth
fig, ax = plt.subplots(figsize=(12,8))
#
for N in Nuclides_long:
if N.heat > 0: ax.loglog(time * facx, N.activity_becq(time) * N.heatJoule() * facy, '--', label=N.name)
for N in Nuclides_special:
if N.heat > 0: ax.loglog(time * facx, N.activity_becq(time) * N.heatJoule() * facy, '-.', label=N.name)
for N in Nuclides_short2:
if N.heat > 0: ax.loglog(time * facx, N.activity_becq(time) * N.heatJoule() * facy, label=N.name)
ax.loglog(time * facx, total_watts * facy/Mearth, 'r-', linewidth=3, label='TOTAL')
#
ax.set_xlabel('Time (million years)')
ax.set_ylabel('Power (W/kg)')
ax.set_title('Radiogenic power (W/kg) from natural radionuclides')
ax.yaxis.set_ticks_position('both')
ax.tick_params(labeltop=False, labelright=True)
ax.legend(bbox_to_anchor=(1.06, 1), loc=2, borderaxespad=0.)
ax.set_position([0.1,0.1,0.75,0.8])
if case==0:
ax.set_xlim(time[0]*facx,time[-1]*facx)
if select_composition==1:
ax.set_ylim(1e-14,1e-6)
fig.savefig('Fig_BE_specpow_nuclides.pdf')
if select_composition==2:
ax.set_ylim(1e-14,2e-6)
fig.savefig('Fig_TTG_specpow_nuclides.pdf')
elif case==1:
ax.set_xlim(0.1,100)
if select_composition==1:
ax.set_ylim(1e-14,1e-6)
fig.savefig('Fig_BE_specpow_nuclides_zoom.pdf')
if select_composition==2:
ax.set_ylim(1e-14,2e-6)
fig.savefig('Fig_TTG_specpow_nuclides_zoom.pdf')
return
figure_specpower(0)
#figure_specpower(1)
Specific power of gamma radiation¶
In [80]:
total_watts_gamma = np.zeros(np.size(time))
for N in Nuclides_all:
total_watts_gamma += N.activity_becq(time) * N.heat_gammaJoule() * Mearth
def figure_power_gamma(case = 0):
facx = 1e-6
facy = 1 #1e-12 * Mearth
fig, ax = plt.subplots(figsize=(12,8))
#
for N in Nuclides_long:
if N.heat_gamma > 0: ax.loglog(time * facx, N.activity_becq(time) * N.heat_gammaJoule() * facy, '--', label=N.name)
for N in Nuclides_special:
if N.heat_gamma > 0: ax.loglog(time * facx, N.activity_becq(time) * N.heat_gammaJoule() * facy, '-.', label=N.name)
for N in Nuclides_short2:
if N.heat_gamma > 0: ax.loglog(time * facx, N.activity_becq(time) * N.heat_gammaJoule() * facy, label=N.name)
ax.loglog(time * facx, total_watts_gamma * facy/Mearth, 'r-', linewidth=3, label='TOTAL')
#
ax.set_xlabel('Time (million years)')
ax.set_ylabel('Power (W/kg)')
ax.set_title('Power of gamma radioactivity (W/kg) from natural radionuclides')
ax.yaxis.set_ticks_position('both')
ax.tick_params(labeltop=False, labelright=True)
ax.legend(bbox_to_anchor=(1.06, 1), loc=2, borderaxespad=0.)
ax.set_position([0.1,0.1,0.75,0.8])
if case==0:
ax.set_xlim(time[0]*facx,time[-1]*facx)
if select_composition==1:
ax.set_ylim(1e-14,1e-6)
fig.savefig('Fig_BE_specpow_gamma_nuclides.pdf')
elif select_composition==2:
ax.set_ylim(1e-15,2e-6)
fig.savefig('Fig_TTG_specpow_gamma_nuclides.pdf')
elif case==1:
ax.set_xlim(0.1,100)
if select_composition==1:
ax.set_ylim(1e-14,1e-6)
fig.savefig('Fig_BE_specpow_gamma_nuclides_zoom.pdf')
elif select_composition==2:
ax.set_ylim(1e-15,2e-6)
fig.savefig('Fig_TTG_specpow_gamma_nuclides_zoom.pdf')
return
figure_power_gamma(0)
#figure_power_gamma(1)
(Anti)neutrino luminosity¶
In [81]:
total_lumanu = np.zeros(np.size(time))
total_lumnu = np.zeros(np.size(time))
for N in Nuclides_all:
total_lumanu += N.activity_becq(time) * N.lanu * Mearth
total_lumnu += N.activity_becq(time) * N.lnu * Mearth
total_lum = total_lumanu + total_lumnu
facx = 1e-6
facy = 1 * Mearth
def figure_lumanu():
fig, ax = plt.subplots(figsize=(12,8))
#
for N in Nuclides_long:
if N.lanu > 0: ax.loglog(time * facx, N.activity_becq(time) * N.lanu * facy, '--', label=N.name)
for N in Nuclides_special:
if N.lanu > 0: ax.loglog(time * facx, N.activity_becq(time) * N.lanu * facy, '-.', label=N.name)
for N in Nuclides_short2:
if N.lanu > 0: ax.loglog(time * facx, N.activity_becq(time) * N.lanu * facy, label=N.name)
ax.loglog(time * facx, total_lumanu * facy/Mearth, 'r-', linewidth=3, label='TOTAL')
#
ax.set_xlim(time[0]*facx,time[-1]*facx)
ax.set_ylim(3e24,4e30)
ax.set_xlabel('Time (million years)')
ax.set_ylabel('Luminosity (particles/sec)')
ax.set_title('Antineutrino luminosity (particles/sec) from natural radionuclides (antineutrinos only)')
ax.yaxis.set_ticks_position('both')
ax.tick_params(labeltop=False, labelright=True)
ax.legend(bbox_to_anchor=(1.06, 1), loc=2, borderaxespad=0.)
ax.set_position([0.1,0.1,0.75,0.8])
if select_composition==1: fig.savefig('Fig_BE_luminosity_anu_nuclides.pdf')
def figure_lumnu():
fig, ax = plt.subplots(figsize=(12,8))
#
for N in Nuclides_long:
if N.lnu > 0: ax.loglog(time * facx, N.activity_becq(time) * N.lnu * facy, '--', label=N.name)
for N in Nuclides_special:
if N.lnu > 0: ax.loglog(time * facx, N.activity_becq(time) * N.lnu * facy, '-.', label=N.name)
for N in Nuclides_short2:
if N.lnu > 0: ax.loglog(time * facx, N.activity_becq(time) * N.lnu * facy, label=N.name)
ax.loglog(time * facx, total_lumnu * facy/Mearth, 'r-', linewidth=3, label='TOTAL')
#
ax.set_xlim(time[0]*facx,time[-1]*facx)
ax.set_ylim(3e24,4e30)
ax.set_xlabel('Time (million years)')
ax.set_ylabel('Luminosity (particles/sec)')
ax.set_title('Neutrino luminosity (particles/sec) from natural radionuclides (antineutrinos only)')
ax.yaxis.set_ticks_position('both')
ax.tick_params(labeltop=False, labelright=True)
ax.legend(bbox_to_anchor=(1.06, 1), loc=2, borderaxespad=0.)
ax.set_position([0.1,0.1,0.75,0.8])
if select_composition==1: fig.savefig('Fig_BE_luminosity_nu_nuclides.pdf')
def figure_lum_total():
fig, ax = plt.subplots(figsize=(12,8))
#
ax.loglog(time*facx, total_lumanu * facy/Mearth, label='Antineutrinos')
ax.loglog(time*facx, total_lumnu * facy/Mearth, '--', label='Neutrinos')
#
ax.set_xlim(time[0]*facx,time[-1]*facx)
ax.set_ylim(3e24,4e30)
ax.set_xlabel('Time (million years)')
ax.set_ylabel('Luminosity (particles/sec)')
ax.set_title('Total (anti)neutrino luminosity (particles/sec) from natural radionuclides')
ax.yaxis.set_ticks_position('both')
ax.tick_params(labeltop=False, labelright=True)
ax.legend(bbox_to_anchor=(1.06, 1), loc=2, borderaxespad=0.)
ax.set_position([0.1,0.1,0.75,0.8])
if select_composition==1: fig.savefig('Fig_BE_luminosity_total.pdf')
def figure_lum_select():
select_lum = np.zeros(np.size(time))
fig, ax = plt.subplots(figsize=(12,8))
#
for N in [Al26, Mn53, K40]:
ax.loglog(time * facx, N.activity_becq(time) * N.lnu * facy, '--', label=N.name)
select_lum += N.activity_becq(time) * N.lnu * facy
for N in [Fe60, K40, Be10, Tc98, Sn126, Cs135, Pu244]:
ax.loglog(time * facx, N.activity_becq(time) * N.lanu * facy, '-', label=N.name)
select_lum += N.activity_becq(time) * N.lanu * facy
for N in [Th232, U235, U238, Rb87]:
ax.loglog(time * facx, N.activity_becq(time) * N.lanu * facy, '.', label=N.name)
select_lum += N.activity_becq(time) * N.lanu * facy
#
rest_lum = total_lum - select_lum
ax.loglog(time * facx, rest_lum, ':', label='All other')
#
ax.loglog(time * facx, total_lum * facy/Mearth, 'r-', linewidth=3, label='TOTAL')
#
ax.set_xlim(1,time[-1]*facx)
ax.set_ylim(1e23,3e30)
ax.set_xlabel('Time (million years)')
ax.set_ylabel('Luminosity (particles/sec)')
ax.set_title('Neutrino (dashed) and antineutrino luminosity (particles/sec) from natural radionuclides')
ax.yaxis.set_ticks_position('both')
ax.tick_params(labeltop=False, labelright=True)
ax.legend(bbox_to_anchor=(1.06, 1), loc=2, borderaxespad=0.)
ax.set_position([0.1,0.1,0.75,0.8])
if select_composition==1: fig.savefig('Fig_BE_luminosity_both_select.pdf')
if select_composition==1:
figure_lumanu()
figure_lumnu()
figure_lum_total()
figure_lum_select()
In [82]:
print ("Present-day radiogenic power (TW) = {:.2f}".format(total_watts[-1]*1e-12))
print ("Present-day antineutrino luminosity (1/sec) = {:.2e}".format(total_lumanu[-1]))
print ("Present-day neutrino luminosity (1/sec) = {:.2e}".format(total_lumnu[-1]))
print ("Present-day tot (anu+nu) luminosity (1/sec) = {:.2e}".format(total_lum[-1]))
By nuclide¶
In [83]:
print ("{:8s} {:7s} {:9s} {:s}".format("", "TW", "anu/sec", "nu/sec"))
for N in Nuclides_long:
tw_now = N.activity_becq(time[-1]) * N.heatJoule() * Mearth * 1e-12 # in TW
anu_now = N.activity_becq(time[-1]) * N.lanu * Mearth # in sec-1
nu_now = N.activity_becq(time[-1]) * N.lnu * Mearth # in sec-1
print ("{:8s} {:.4f} {:.2e} {:.2e}".format(N.name, tw_now, anu_now, nu_now))
Values at t = 1 Ma¶
In [84]:
t1 = 1e6
print ("{:5s} {:9s} {:9s} {:s}".format("", "TW", "anu/sec", "nu/sec"))
for N in Nuclides_all:
tw1 = N.activity_becq(t1) * N.heatJoule() * Mearth * 1e-12 # in TW
anu1 = N.activity_becq(t1) * N.lanu * Mearth # in sec-1
nu1 = N.activity_becq(t1) * N.lnu * Mearth # in sec-1
print ("{:5s} {:.3e} {:.3e} {:.3e}".format(N.name, tw1, anu1, nu1))
Total radiogenic energy delivered by 26Al and 60Fe¶
In [85]:
for N in Al26, Fe60:
deccsec = log(2)/N.halflife/secinyr
totErad = N.activity_becq(0) * N.heatJoule() * Mearth / deccsec
print ("{:5s} {:.3e} TW".format(N.name, totErad))
Coefficients for plug-in formulae for radiogenic power, neutrino luminosity¶
We want to write
$h \text{ [W kg-rock}^{-1}] = \eta_\text{K}\text{[K]} + \eta_\text{Rb}\text{[Rb]} + \eta_\text{Sm}\text{[Sm]} + \eta_\text{Th}\text{[Th]} + \eta_\text{U}\text{[U]}$
where [K] is mass fraction of K (in kg-elem/kg-rock) and therefore
$\eta \text{ [W kg-elem}^{-1}] = \frac{N_A X \lambda Q_h}{\mu}$
for each radionuclide, where $N_A$ is Avogadro's number, $X$ is the nuclide's natural molar isotopic fraction, $\mu$ is elemental molar mass, $\lambda$ is decay constant, and $Q_h$ is radiogenic heat released per decay. Below we evaluate the coefficients $\eta$ for various radionuclides.
Similarly for a plugin formula to calculate antineutrino and neutrino luminosities.
In [86]:
print("COEFFICIENTS FOR LONG-LIVED RADIONUCLIDES\n")
fac = 1e9
print("h-fac_K = {:5.3f} x 1e-9 watts/kg-elem".format(K40.facHeat()*fac))
print("h-fac_Rb = {:5.2f} x 1e-9".format(Rb87.facHeat()*fac))
print("h-fac_Sm = {:5.2f} x 1e-9".format(Sm147.facHeat()*fac))
print("h-fac_Th = {:.0f} x 1e-9".format(Th232.facHeat()*fac))
U_facHeat = U235.facHeat() + U238.facHeat()
print("h-fac_U = {:.0f} x 1e-9".format(U_facHeat*fac))
print("\nCheck W/kg in BSE: {:.3e}".format((K40.facHeat()*280e-6+Rb87.facHeat()*0.6e-6+Sm147.facHeat()*406e-9+Th232.facHeat()*20e-9*3.77+U_facHeat*20e-9)))
print("\nCheck TW in BSE: {:.2f}\n".format((K40.facHeat()*280e-6+Rb87.facHeat()*0.6e-6+Sm147.facHeat()*406e-9+Th232.facHeat()*20e-9*3.77+U_facHeat*20e-9)*Mbse*1e-12))
fac = 1e-4
print("lanu-fac_K = {:5.3f} x 1e4 antineutrinos/second/kg-elem".format(K40.facLanu()*fac))
print("lanu-fac_Rb = {:5.2f} x 1e4".format(Rb87.facLanu()*fac))
print("lanu-fac_Th = {:5.0f} x 1e4".format(Th232.facLanu()*fac))
U_facLanu = U235.facLanu() + U238.facLanu()
print("lanu-fac_U = {:5.0f} x 1e4".format(U_facLanu*fac))
print("\nCheck Lanu in BSE: {:.2e} per second\n".format((K40.facLanu()*280e-6+Rb87.facLanu()*0.6e-6+Th232.facLanu()*20e-9*3.77+U_facLanu*20e-9)*Mbse))
fac = 1e-4
print("lnu-fac_K = {:.4f} x 1e4 neutrinos/second/kg-elem".format(K40.facLnu()*fac))
print("\nCheck Lnu in BSE: {:.2e} per second".format((K40.facLnu()*280e-6)*Mbse))
In [87]:
def print_facHeat(Nucl = None):
print("{:14s}{:.3e}".format(Nucl.name, Nucl.facHeat()*fac))
fac = 1e9
print("COEFFICIENTS FOR SHORT-LIVED RADIONUCLIDES\n")
print("coefficients (h/A) in nanowatts per kg-elem")
for N in Nuclides_short:
print_facHeat(N)
print("\n 36Cl beta- {:.3e}".format(Cl36.facHeat()*fac * Cl36.heat_beta/Cl36.heat))
print(" 36Cl eps {:.3e}".format(Cl36.facHeat()*fac * (Cl36.heat_ec+Cl36.heat_betaplus)/Cl36.heat))
print("36Cl {:.3e}".format(Cl36.facHeat()*fac))
print("\n 26Al EC {:.3e}".format(Al26.facHeat()*fac * Al26.heat_ec/Al26.heat))
print(" 26Al beta+ {:.3e}".format(Al26.facHeat()*fac * Al26.heat_beta/Al26.heat))
print("26Al {:.3e}".format(Al26.facHeat()*fac))
In [ ]: