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

  • 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

  1. the $Q$ value of the β transition is much smaller than the rest mass energies of the parent and daughter nuclei,
  2. 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)
$$w_0 = w+q = 1+E+q = 1+Q \quad \text{or} \quad E+q=Q$$

The shape of a β spectrum is proportional to:

  1. phase space factor, $pwq^2$,
  2. Fermi function, $F(Z,w)$, which accounts for the Coulombic interaction between the nucleus and outgoing β-particle,
  3. shape factor, $S(w)$, which accounts for the "forbidden" nature of some transitions.
$$\frac{dN}{dw} \propto pwq^2 F(Z,w) S(w)$$

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):

bielajew_bookchapter15_tab15.1

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
Mass of BSE in kg = 4.0422e+24  (= Mearth - Mcore, and references just above)
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).

40K beta-decay scheme

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
branching = 0.8944 [Naumenko-Dèzes et al 2018]
branching = 0.8959 [Renne et al 2001]
branching = 0.8928 [NuDat]

half-life = 1.262e+09 yr [Naumenko-Dèzes et al 2018]
half-life = 1.253e+09 yr [Renne et al 2011]
half-life = 1.248e+09 yr [NuDat]

Q = 1310.89 keV  (NuDat: 1310.89 ± 0.06 keV)
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))
Mean neutrino energy/total available energy = 0.5924
Mean electron energy/total available energy = 0.4076
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))
Mean neutrino energy/total available energy = 0.6125
Mean electron energy/total available energy = 0.3875

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))
Mean antineutrino energy = 726.8 keV using experimental shape factor
Mean antineutrino energy = 722.9 keV using theoretical  shape factor

Accounting for branching, the β- antineutrino carries away on average, in keV:
    650.0 (Naumenko, experimental)
    651.1 (Renne,    experimental)
    648.8 (NUDAT,    experimental)
    646.5 (Naumenko, theoretical)
    647.6 (Renne,    theoretical)
    645.4 (NUDAT,    theoretical)

Neutrinos per decay = 0.8944000000000001 (check of normalization)

Mean β− energy = 584.1 keV  (Renne: 584.1 keV, NuDat: 560.2 ± 1.5 keV)
Accounting for branching, the β- antineutrino carries away on average 650.0 keV (Naumenko & exp)
Accounting for branching, the β- heat per decay is 522.5 keV (Naumenko & exp)

Electron capture on 40K

40K EC-decay scheme

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))
Q = 1504.40 keV  (NuDat: 1504.40 ± 0.06 keV)

Accounting for branching, the EC neutrino carries away on average 5.276 keV
Accounting for branching, the EC heat per decay is 153.6 keV

β+ decay of 40K

(1.00 ± 0.13) × 10-5 branch with mean β+ energy 200 ± 40 keV (NuDat)

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))
Q = 482.40 keV  (NuDat: 482.40 ± 0.06 keV)

Mean positron energy = 197.325 keV (from NuDat)
Accounting for branching, the β+ heat per decay is 0.012 keV

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)")
Overall Q-value = 1331.3 keV = 0.2133 pJ 

β- decay: mean electron energy   = 584.1 keV/decay
          mean nuebar energy     = 726.8 keV/decay
          carried away by nuebar = 650.0 keV/decay
          heats the Earth        = 522.5 keV/decay

EC:       carried away by nue    =   5.3 keV/decay
          heats the Earth        = 153.6 keV/decay

β+ decay: heats the Earth        = 0.012 keV/decay

Overall:  heats the Earth        = 676.0 keV/decay = 0.1083 pJ/decay

(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)

87RB_beta_decay_scheme

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))
Q = 282.3 keV  (NuDat: 282.2 ± 1.1 keV)

Using NuDat inputs:
    Mean electron energy = 81.67 keV  (from NuDat)

My calculation:
    Mean electron energy w/ shape factor equal to 1 = 81.87 keV
    Mean electron energy w/ shape from Mougeot      = 56.72 keV

Heats the Earth = 81.9 keV/decay = 0.01312 pJ/decay

138La decay

Branched β or EC decay, half-life 102 ± 1 Gyr (NuDat)

β decay 138La → 138Ce

138LA_beta_decay_scheme

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))
Q = 1046.9 keV  (NuDat: 1052 ± 4 keV)

Using NuDat inputs:
    Mean electron energy = 98.00 keV  (from NuDat)
    Heats the Earth = 305.91 keV/decay (from NuDat)

EC decay 138La → 138Ba

138LA_eps_decay_scheme

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))
Q = 1739.9 keV  (NuDat: 1742 ± 3 keV)

Using NuDat inputs
    Heats the Earth = 940.45 keV/decay (from NuDat)

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))
Qmean = 1504.0
Heats the Earth = 1246.36 keV/decay = 0.19969 pJ/decay (from NuDat)

147Sm → 143Nd decay

α decay with half-life 107 ± 1 Gyr (NuDat)

147SM_alpha_decay_scheme

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))
Q = 2311.17 keV  (NuDat: 2311 ± 1 keV)

Heats the Earth = 2311.2 keV/decay = 0.3703 pJ/decay

176Lu → 176Hf

β- decay with half-life 37.6 ± 0.7 Gyr (NuDat)

176LU_beta_decay_scheme

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))
Q = 1194.3 keV  (NuDat: 1190.2 ± 0.8 keV)

Using NuDat inputs:
    Mean electron energy = 181.80 keV  (from NuDat)

Heats the Earth = 778.8 keV/decay = 0.12477 pJ/decay

187Re → 187Os

β- decay with half-life 43.3 ± 0.7 Gyr (NuDat)

187RE_beta_decay_scheme

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))
Q = 2.515 keV  (NuDat: 2.469 ± 0.004 keV)

Using NuDat inputs:
    Mean electron energy = 0.62 keV  (from NuDat)

Heats the Earth = 0.6 keV/decay = 0.00010 pJ/decay

190Pt → 186Os

α decay with half-life 650 ± 3 Gyr (NuDat)

190PT_alpha_decay_scheme

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))
Q = 3252.26 keV  (NuDat: 3249 ± 6 keV)

Heats the Earth = 3252.3 keV/decay = 0.5211 pJ/decay

232Th → ... → 208Pb decay chain

Decay network with 6 α and 4 β decays per one decayed parent atom

Decay_chain_232Th

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))
Qtotal = 42647 keV = 6.8329 pJ

Check integrated n_anu = 3.9714 (should be 4)

Integrated E_anu = 2229.9 keV/decay = 0.3573 pJ/decay

Heats the Earth = 40417.4 keV/decay = 6.4756 pJ/decay

235U → ... → 207Pb decay chain

Decay network with 7 α and 6 β decays per one decayed parent atom

Decay_chain_235U

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))
Qtotal = 46398 keV = 7.4337 pJ

Check integrated n_anu = 4.0226 (should be 4)

Integrated E_anu = 2019.7 keV/decay = 0.3236 pJ/decay

Heats the Earth = 44378.0 keV/decay = 7.1101 pJ/decay

238U → ... → 206Pb decay chain

Decay network with 8 α and 6 β decays per one decayed parent atom

Decay_chain_238U

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))
Qtotal = 51695 keV = 8.2824 pJ

Check integrated n_anu = 6.0077 (should be 6)

Integrated E_anu = 4049.9 keV/decay = 0.6489 pJ/decay

Heats the Earth = 47644.9 keV/decay = 7.6336 pJ/decay

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)
Original [NIST]: 0.0072037 0.992796 238.02891

New [Villa16+]:   0.007204944512094626 0.9927404845434701 238.02890533224976

Short-lived radionuclides

41Ca → 41K

ε decay (EC only; no β+) with half-life 99.4 ± 1.5 kyr (NuDat)

41CA_eps_decay_scheme

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))
Q = 421.65 keV  (NuDat: 421.66 ± 0.14 keV)

Heats the Earth = 0.00 keV/decay = 0.00000 pJ/decay

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)

99TC_beta_decay_scheme

99TC_beta_decay_scheme_NDS2017

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))
Q = 295.0 keV  (NuDat: 297.5 ± 1.0 keV)

beta- decay to ground state:
    Neutrinos per decay = 0.9999840000000002 (check of normalization)
    Mean electron energy   = 95.69 keV/decay (NuDat: 84.6 ± 0.5 keV, Mougeot: 95.28/101.03 keV)
    Mean nuebar energy     = 199.3 keV/decay
    Carried away by nuebar = 199.3 keV/decay
    Heats the Earth        = 95.687 keV/decay = 0.01533 pJ/decay

beta- decay to excited state 89.50 kev:
    Neutrinos per decay = 1.6e-05
    Mean electron energy   = 81.70 keV/decay (from NuDat, resp. Browne & Tuli 2017 NDS)
    Heats the Earth        = 0.001 keV/decay = 0.0000002 pJ/decay

Overall heats the Earth = 95.689 keV/decay = 0.0153310 pJ/decay

81Kr → 81Br

ε decay (EC only; no β+) with half-life 0.229 ± 0.011 Myr (NuDat references Baglin 2008 NDS)

81KR_eps_decay_scheme

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))
Q = 280.8 keV  (NuDat: 280.8 ± 0.5 keV)

7.16395161
Heats the Earth = 0.828 keV/decay = 0.000133 pJ/decay

126Sn → 126Sb → 126Te

Two-step β- decay

126SN_decay_scheme

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

126SN_beta_decay_scheme

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)
Q = 378.2 keV  (NuDat: 378 ± 30 keV)

Endpoint energy = 250.3 keV  (NuDat: 250 ± 30 keV)

Mean electron energy   = 70.47 keV (NuDat: 70 ± 14 keV)
Heats the Earth = 183.1 keV/decay = 0.02934 pJ/decay

126Sb → 126Te from ground state

At 14% branching the β– decay starts at ground state of 126Sb, half-life 12.35 ± 0.06 days (NuDat)

126SB_beta_decay_scheme0

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
Q = 3672.0 keV  (NuDat: 3673 ± 33 keV)

Sum of branching ratios = 1.0429
    renormalizing to branching 1.0000 (NuDat)
Sum of branching ratios = 1.0000
Mean electron energy = 340.3 keV  (NuDat: 340 ± 60 keV)
Deexcitation energy per decay = 2711.8 keV
Heats the Earth = 427.3 keV/decay = 0.06846 pJ/decay

126Sb → 126Te from 17.7 keV level

At 86% branching the β– decay starts at 17.7keV level of 126Sb, half-life 19.15 ± 0.08 mins (NuDat)

126SB_beta_decay_scheme1

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)
Q = 3672.0 keV  (NuDat: 3673 ± 33 keV)

Sum of branching ratios = 0.8856
Mean electron energy = 639.6 keV  (NuDat: 720 ± 90 keV)
Deexcitation energy per decay = 1609.7 keV
Heats the Earth = 2249.2 keV/decay = 0.36037 pJ/decay

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))
126Sn       → 126Sb rad. heat =  183.1 keV/decay = 0.02934 pJ/decay
126Sb(GS)   → 126Te rad. heat =  427.3 keV/decay = 0.06846 pJ/decay
126Sb(17.7) → 126Te rad. heat = 2249.2 keV/decay = 0.36037 pJ/decay

Overall heats the Earth = 2859.7 keV/decay = 0.45817 pJ/decay

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

36CL_beta_decay_scheme

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))
Q = 709.523 keV  (NuDat: 709.547 ± 0.046 keV)

Uniform (no) shape factor:
    Mean electron energy = 250.86 keV (NuDat: 251.33 keV)
    Heats the Earth = 246.1 keV/decay = 0.03943 pJ/decay

Experimental shape factor from Mougeot:
    Mean electron energy = 340.80 keV (NuDat: 251.33 keV)
    Heats the Earth = 334.3 keV/decay = 0.05356 pJ/decay

36Cl → 36S ε decay

EC (1.89%) and β+ (0.014%) decay, both ground state to ground state

36CL_eps_decay_scheme

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)
Q(EC) = 1142.10 keV  (NuDat: 1142.14 ± 0.19 keV)
Q(beta+) = 120.11 keV  (NuDat: 120.14 ± 0.19 keV)
Heats the Earth (beta+) = 0.1501 keV/decay = 0.0000241 pJ/decay

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))
Heats the Earth = 334.47 keV/decay = 0.05359 pJ/decay

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.

79BR_beta_decay_scheme

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))
Q = 150.613 keV  (NuDat: 150.6 ± 1.3 keV)

Uniform (no) shape factor:
    Mean electron energy = 40.75 keV (NuDat: 52.8 keV)
    Heats the Earth = 40.7 keV/decay = 0.00653 pJ/decay

Symmetrical shape factor from Bielajew
    Mean electron energy = 55.95 keV (NuDat: 52.8 keV)
    Heats the Earth = 55.9 keV/decay = 0.00896 pJ/decay

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

26AL_eps_decay_scheme

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))
Q(EC) = 4004.43 keV  (NuDat: 4004.43 ± 0.06 keV)
Heats the Earth = 361.04 keV/decay = 0.05785 pJ/decay

Q(beta+) = 2982.44 keV  (NuDat: n/a)
Eendpoint = 1173.72 keV  (NuDat: 1173.42 ± 0.09 keV)
513.5571983792743
Mean positron energy = 545.36 keV  (NuDat: 543.29 keV)
Heats the Earth = 2759.27 keV/decay = 0.44208 pJ/decay

Overall heats the Earth = 3120.31 keV/decay = 0.49993 pJ/decay
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()
Looking what happens when I switch sign of Z_daughter in Fermi function...

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

10BE_beta_decay_scheme

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
Q = 556.8 keV  (NuDat: 556.0 ± 0.6 keV)

Using NuDat inputs:
    Mean electron energy = 202.56 keV  (from NuDat)

My calculation:
    Mean electron energy w/ shape factor equal to 1 = 202.78 keV
    Mean electron energy w/ shape from Mougeot      = 252.65 keV

**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))
Heats the Earth = 252.65 keV/decay = 0.04048 pJ/decay

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

93ZR_beta_decay_scheme

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))
Q = 90.3 keV  (NuDat: 90.8 ± 1.6 keV)

Mean electron energy from Nudat = 20 ± 3 keV  (from NuDat)

Mean electron energy with both shape factors equal to one = 17.88 keV
Mean electron energy calculated theoretical excited level SF = 23.14 keV (using this)

Heats the Earth = 45.6 keV/decay = 0.0073 pJ/decay

154Dy → 150Gd

α decay with half-life 3.0 ± 1.5 Myr (NuDat)

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))
Q = 2945.05 keV  (NuDat: 2945 ± 5 keV)

Heats the Earth = 2945.1 keV/decay = 0.4718 pJ/decay

150Gd → 146Sm

α decay with half-life 1.79 ± 0.08 Myr (NuDat)

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))
Q = 2807.66 keV  (NuDat: 2808 ± 6 keV)

Heats the Earth = 2807.7 keV/decay = 0.4498 pJ/decay

146Sm → 142Nd

α decay with half life of 103 ± 5 Myr according to NuDat or 68 ± 7 Myr according to Bill. Using the NuDat value.

146SM_alpha_decay_scheme

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))
Q = 2528.8 keV  (NuDat: 2528 ± 3 keV)

Heats the Earth = 2528.8 keV/decay = 0.4052 pJ/decay

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)

135CS_beta_decay_scheme

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))
Q = 268.8 keV  (NuDat: 268.7 ± 1.1 keV)

Using NuDat inputs:
    Mean electron energy = 75.7 keV  (from NuDat)

My calculation:
    Mean electron energy w/ shape factor equal to 1                = 75.64 keV
    Mean electron energy w/ experimental shape factor from Mougeot = 61.51 keV  (Mougeot: 61.57 keV)
    Mean electron energy w/ calc. by Mougeot using ξ approximation = 88.44 keV

I will use my (~Mougeot) calculation with experimental shape factor:

Heats the Earth = 61.51 keV/decay = 0.009855 pJ/decay

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.

60FEdecay_scheme 60COdecay_scheme 60FE_beta_decay_scheme_NDS2013 60NI_IT_decay_scheme_NDS2013 60NI5+_beta_decay_scheme_NDS2013 60NI2+_beta_decay_scheme_NDS2013

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))
Q = 237.3 keV  (NuDat: 237 ± 3 keV)

Mean electron energy in 60Fe decay w/ shape factor equal to 1 = 50.12 keV   (NuDat: 50.2 keV)

Overall heats the Earth = 2707.7 keV/decay = 0.43382 pJ/decay

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]

53MN_eps_decay_scheme

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))
Q = 596.8 keV  (NuDat: 596.8 ± 0.4 keV)

Heats the Earth = 0.00 keV/decay = 0.00000 pJ/decay

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??

98TC_beta_decay_scheme 98TC_beta_decay_scheme_NDS2003

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))
Q = 1794 keV  (NuDat: 1796 ± 7 keV)

Using NuDat inputs:
    Mean electron energy = 118 keV  (from NuDat)

My calculation:
    Electron endpoint energy = 395.92 keV  (NuDat: 397 ± 22 keV)
    Mean electron energy w/ shape factor equal to 1 = 118.69 keV

Using my calculation with shape factor equal to one (...can we do better??):

Heats the Earth = 1516.46 keV/decay = 0.242964 pJ/decay

97Tc → 97Mo

ε decay (EC only; no β+) with half-life 4.21 ± 0.16 Myr, to ground state (NuDat references Nica 2010 NDS)

97TC_eps_decay_scheme

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))
Q = 324.7 keV  (NuDat: 320 ± 4 keV)

Heats the Earth = 0.00 keV/decay = 0.00000 pJ/decay

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?)

107PD_beta_decay_scheme

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))
Q = 34.1 keV  (NuDat: 34.1 ± 2.7 keV)

Mean electron energy from NuDat = 9 keV

Electron endpoint energy = 34.09 keV  (NuDat: 34 ± 3 keV)

My calculation:
    Mean electron energy w/ shape factor equal to one = 8.74 keV
    Mean electron energy w/ theoretical shape factor = 13.27 keV

Using my calculation with theoretical shape factor:

Heats the Earth = 13.27 keV/decay = 0.002126 pJ/decay

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>

182HF_beta_decay_scheme 182HF_beta_decay_scheme_NDS2010

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))
Q = 381.3 keV  (NuDat: 381 ± 6 keV)

Using NuDat inputs:
    Mean electron energy =  32.9 keV
    Heats the Earth =  303.3 keV

My calculation:
    Electron endpoint energy = 110.85 keV
    Mean electron energy w/ shape factor equal to 1 = 41.11 keV
    Heats the Earth =  311.5 keV

Using NuDat-derived value:

Heats the Earth = 311.5 keV/decay = 0.04991 pJ/decay

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??

182TA_beta_decay_scheme

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))
Q = 1814.5 keV  (NuDat: 1814.5 ± 1.7 keV)

Sum of branching ratios = 1.0031
Mean electron energy = 127.8 keV  (NuDat: 127 ± 3 keV)
Deexcitation energy per decay = 1388.3 keV
Heats the Earth = 1516.1 keV/decay = 0.24290 pJ/decay

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))
182Hf → 182Ta rad. heat = 311.5 keV/decay = 0.04991 pJ/decay
182Ta → 182W  rad. heat = 1516.1 keV/decay = 0.24290 pJ/decay

Overall heats the Earth = 1827.6 keV/decay = 0.29281 pJ/decay

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.

129I_beta_decay_scheme 129I_beta_decay_scheme_NDS2014

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))
Q = 188.9 keV  (NuDat: 189 ± 3 keV)

Using NuDat inputs:
    Mean electron energy =  40.0 keV

My calculation:
    Electron endpoint energy = 149.37 keV
    Mean electron energy w/ shape factor equal to 1          = 40.32 keV
    Mean electron energy w/ exp. shape factor from Mougeot   = 45.74 keV  (Mougeot: 45.98 keV)
    Mean electron energy w/ calc. by Mougeot using ξ approx. = 48.62 keV

I will use my (~Mougeot) calculation with experimental shape factor:

Heats the Earth = 85.3 keV/decay = 0.01367 pJ/decay

205Pb → 205Tl

ε decay (EC only; no β+) with half-life 17.3 ± 0.7 Myr, to ground state [NuDat references Kondev 2004 NDS]

205PB_eps_decay_scheme

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))
Q = 50.7 keV  (NuDat: 50.5 ± 0.5 keV)

Heats the Earth = 0.0 keV/decay = 0.00000 pJ/decay

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]

92NB_eps_decay_scheme 92NB_eps_decay_scheme_NDS2012

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))
Q = 2005.9 keV  (NuDat: 2005.9 ± 1.8 keV)

Heats the Earth = 1495.6 keV/decay = 0.23962 pJ/decay

244Pu → 240U → 240Np → 240Pu → 236U → 232Th

α → β- → β- → α → α

244Pu → 240U

α decay with half-life 81.1 ± 0.3 Myr (NuDat)

244PU_alpha_decay_scheme

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))
Q = 4665.5 keV  (NuDat: 4665.5 ± 1.0 keV)

Heats the Earth = 4665.5 keV/decay = 0.7475 pJ/decay

240U → 240Np

β- decay with half-life 14.1 ± 0.1 hours (NuDat)

240U_beta_decay_scheme

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))
Q = 399.1 keV  (NuDat: 400 ± 16 keV)

Sum of branching ratios = 1.0378
    renormalizing to 1.000...
Mean electron energy = 102.6 keV  (NuDat: 103 ± 13 keV)
Deexcitation energy per decay = 15.9 keV
Heats the Earth = 118.5 keV/decay = 0.01899 pJ/decay

240Np → 240Pu

β- decay with half-life 7.22 ± 0.02 mins (NuDat)

240NP_beta_decay_scheme

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))
Q = 2190.1 keV  (NuDat: 2188 ± 15 keV)

Sum of branching ratios = 0.9755
    renormalizing to branching 0.9988 (NuDat)
Sum of branching ratios = 0.9988
Mean electron energy = 644.4 keV  (NuDat: 650 ± 50 keV)
Deexcitation energy per decay = 369.3 keV
Heats the Earth = 1013.7 keV/decay = 0.16241 pJ/decay

240Pu → 236U

α decay with half-life 6561 ± 7 yr (NuDat)

240PU_alpha_decay_scheme

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))
Q = 5255.81 keV  (NuDat: 5255.75 ± 1.4 keV)

Heats the Earth = 5255.8 keV/decay = 0.8421 pJ/decay

236U → 232Th

α decay with half-life 23.42 ± 0.04 Myr (NuDat)

236U_alpha_decay_scheme

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))
Q = 4572.84 keV  (NuDat: 4573.1 ± 0.9 keV)

Heats the Earth = 4572.8 keV/decay = 0.7326 pJ/decay

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))
244Pu → 240U  rad. heat = 4665.5 keV/decay = 0.74750 pJ/decay
240U  → 240Np rad. heat =  118.5 keV/decay = 0.01899 pJ/decay
240Np → 240Pu rad. heat = 1013.7 keV/decay = 0.16241 pJ/decay
240Pu → 236U  rad. heat = 5255.8 keV/decay = 0.84207 pJ/decay
236U  → 232Th rad. heat = 4572.8 keV/decay = 0.73265 pJ/decay

Overall heats the Earth = 15626.4 keV/decay = 2.50362 pJ/decay
Overall Q = 17083.4 keV/decay = 2.73706 pJ/decay

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)
Radiogenic heat (keV, pJ), antineutrinos, neutrinos per decay
            keV        pJ       anu       nu
Long-lived
40K        676.0    0.1083    0.8944    0.1056    0.2272
87Rb        81.9    0.0131    1.0000    0.0000    0.0000
La138     1246.4    0.1997    0.3450    0.6550    0.9729
147Sm     2311.2    0.3703    0.0000    0.0000    0.0000
176Lu      778.8    0.1248    1.0000    0.0000    0.7665
187Re        0.6    0.0001    1.0000    0.0000    0.0000
190Pt     3252.3    0.5211    0.0000    0.0000    0.0000
232Th    40417.4    6.4756    4.0000    0.0000   -0.0000
235U     44378.0    7.1101    4.0000    0.0000   -0.0000
238U     47644.9    7.6336    6.0000    0.0000   -0.0000
Short-lived
41Ca         0.0    0.0000    0.0000    1.0000        -
99Tc        95.7    0.0153    1.0000    0.0000    0.0000
126Sn     2859.7    0.4582    2.0000    0.0000    0.6966
36Cl       334.5    0.0536    0.9810    0.0190    0.0004
26Al      3120.3    0.4999    0.0000    1.0000    0.8573
10Be       252.6    0.0405    1.0000    0.0000    0.0000
135Cs       61.5    0.0099    1.0000    0.0000    0.0000
60Fe      2707.7    0.4338    2.0000    0.0000    0.9259
53Mn         0.0    0.0000    0.0000    1.0000        -
98Tc      1516.5    0.2430    1.0000    0.0000    0.9316
97Tc         0.0    0.0000    0.0000    1.0000        -
107Pd       13.3    0.0021    1.0000    0.0000    0.0000
182Hf     1827.6    0.2928    2.0000    0.0000    0.8447
129I        85.3    0.0137    1.0000    0.0000    0.2780
205Pb        0.0    0.0000    0.0000    1.0000        -
92Nb      1495.6    0.2396    0.0000    1.0000    0.8441
244Pu    15626.4    2.5036    2.0000    0.0000    0.0208
146Sm     2528.8    0.4052    0.0000    0.0000    0.0000
Short-lived with unknown initial abundance
81Kr         0.8    0.0001    0.0000    1.0000    8.6524
79Se        55.9    0.0090    1.0000    0.0000    0.0000
93Zr        45.6    0.0073    1.0000    0.0000    0.0321
150Gd     2807.7    0.4498    0.0000    0.0000    0.0000
154Dy     2945.1    0.4718    0.0000    0.0000    0.0000

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))
40K   1.90e-04 1.34e-02 70.707260
87Rb  4.06e-07 4.43e-05 109.086175
La138 4.40e-07 2.18e-05 49.551907
147Sm 2.75e-07 2.80e-06 10.189415
176Lu 4.60e-08 1.18e-07 2.563837
187Re 7.50e-08 0.00e+00 0.000000
190Pt 1.90e-06 0.00e+00 0.000000
232Th 5.10e-08 4.81e-06 94.252086
235U  1.35e-08 9.34e-07 68.997622
41Ca  1.71e-02 2.04e-02 1.191316
99Tc  1.30e-06 3.40e-10 0.000262
126Sn 2.50e-07 2.10e-06 8.414437
36Cl  7.61e-05 3.70e-04 4.859520
26Al  1.59e-02 8.33e-02 5.237142
10Be  4.60e-08 2.10e-06 45.627600
135Cs 3.52e-08 1.77e-06 50.256588
60Fe  3.17e-01 1.92e-02 0.060556
53Mn  8.04e-04 3.27e-04 0.406594
107Pd 1.00e-06 5.20e-10 0.000518
182Hf 1.92e-07 4.20e-06 21.927044
129I  4.88e-08 1.40e-06 28.700399
205Pb 2.31e-07 1.16e-05 50.261867
92Nb  4.45e-07 4.45e-06 9.991975
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())
Al mass fraction (present-day) = 0.015905620728109335
26Al initial W/kg = 2.8543948443064884e-07

Fe mass fraction = 0.3170601803696473
60Fe initial W/kg = 1.15237683486641e-10

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()
40K
87Rb
La138
147Sm
176Lu
187Re
190Pt
232Th
235U
238U
41Ca
99Tc
126Sn
36Cl
26Al
10Be
135Cs
60Fe
53Mn
98Tc
97Tc
107Pd
182Hf
129I
205Pb
92Nb
244Pu
146Sm