If you use this program for your work, please acknowledge it and cite the following paper:
Ruedas, T. (2017): Radioactive heat production of six geologically important nuclides. Geochem. Geophys. Geosys., 18(9), 3530-3541, doi:10.1002/2017GC006997.
This software is released under the GNU General Public License v.3.
The purpose of this program is to calculate the energy available as heat or for chemical reactions and the specific heat production of six nuclides ($^{26}$Al, $^{40}$K, $^{60}$Fe, $^{232}$Th, $^{235}$U, and $^{238}$U), which are geologically important as heat sources, especially in the context of the formation, evolution, and geodynamics of terrestrial planetary bodies.
All half-lives are converted into years, with the relation between seconds and years defined as in NUBASE2016 (Audi et al., 2017). Energy calculations use energies in keV. Atomic masses of isotopes are taken from AME2016 (Wang et al., 2017), the elemental averages and isotopic abundances from IUPAC 2013 (Meija et al., 2016a,b), and the half-lives are from NUBASE2016.
The nuclide-specific decay data (energies and intensities) were obtained from the NuDat v.2.7$\beta$ website of the National Nuclear Data Center (NNDC) of the Brookhaven National Laboratory, USA, and the Recommended Data site of the Laboratoire National Henri Becquerel, Saclay, France. In general, the newest available data were used.
The decay chain of the nuclide is followed through most of its branches, leaving out only very minor (improbable) ones in some instances. Total decay energies $Q$ are calculated from the mass difference of parent and decay products of the step in the chain. For $\alpha$ decay, this is equivalent to the energy available for heating. For $\beta^-$ and $\beta^+$ decay, the energy carried by (anti)neutrinos is subtracted; for the latter, the annihilation of the positron and an electron are taken into account as an additional heat-effective energy term. For electron capture (ec), only the radiative part is converted into heat. As all decay series happen to start with nuclides whose half-lives $T_{1/2}$ are orders of magnitude longer than that of any of the daughters, it is the parent half-life that is used for power calculations; the lambda_s decay constants written out for each step are for informational purposes only.
In [ ]:
from math import log
# physical and unit conversion constants
clight0=299792458. # m/s, speed of light, CODATA2014
c2=clight0*clight0
amu=1.660539040e-27 # kg, atomic mass unit, CODATA2014
ec0=1.6021766208e-19 # C, elementary charge, CODATA2014
yr2s=31556926 # yr->s conversion factor (1 yr = 365.2422 d), NUBASE2016
yr2d=365.2422
yr2h=yr2d*24
yr2m=yr2h*60
u_He4=4.00260325413 # AME2016
u_el=9.10938356e-31/amu # m_el/amu, CODATA2014
E_ep0=c2*2*u_el*1e-3*amu/ec0 # electron-positron annihilation energy (in keV)
#
rU85=137.794 # ratio of U-238 to U-235 (Goldmann et al., 2015)
rU84=0.992742/0.000054 # ratio of U-238 to U-234 (IUPAC 2013)
rK19=0.067302/0.932581 # ratio of K-41 to K-39 (IUPAC 2013)
At the core of the calculation is the energy balance of the decay. The following helper functions calculate energies; all include conversions from SI units into keV where necessary, as keV is the commonly used energy unit in nuclear physics.
In $\alpha$ decays, the entire decay energy is converted into heat, which follows from the mass difference between the parent nuclide and the daughter nuclide plus the $^4$He atom that eventually forms from the $\alpha$ particle: \begin{equation} E=c_0^2(m_\mathrm{P}-m_\mathrm{D}-m_{^4\mathrm{He}}). \end{equation}
In [ ]:
def Qa(P,D):
"""Q-value for alpha decay in keV"""
global u
return c2*(u[P]-u[D]-u_He4)*1e-3*amu/ec0
The following is a helper function that calculates the total decay energy for use in $\beta$ or $\epsilon$ (ec) decays according to \begin{equation} E=c_0^2(m_\mathrm{P}-m_\mathrm{D}). \end{equation}
In [ ]:
def Qb(P,D):
"""Q-value for beta/ec decay in keV"""
global u
return c2*(u[P]-u[D])*1e-3*amu/ec0
The heat energies are converted into power per mass using the decay constants $\lambda_{1/2}=\frac{\ln 2}{T_{1/2}}$ of the parent nuclide. \begin{equation} H=\frac{E_H\lambda_{1/2}}{m_\mathrm{P}}. \end{equation} The following function includes a conversion back to SI units.
In [ ]:
def EtoH_W(lam,EH,P):
"""Energy to power conversion in W, with conversion keV->J, yr->s"""
global u
return 1e3*lam*EH*ec0/(u[P]*amu*yr2s)
The decay calculation yields the power per mass of the isotope under consideration, but it is often convenient to know the power per mass of the element, which usually has several isotopes with different total decay heat energies (if they are radioactive). This value changes with time as the proportions of the different isotopes in the isotopic composition of the element change due to the different decay rates. The values for extant isotopes calculated here are present-day values based on observed isotopic abundances and masses: \begin{equation} H_\mathrm{elm}=\frac{1}{\bar{u}}\sum H_i X_i u_i, \end{equation} where $\bar{u}=\sum X_i u_i$ indicates the mean value for the element and the sum is over all of its isotopes.
It must be ensured that $\bar{u}$ and $\sum X_i u_i$ are consistent with each other within the accuracy of the values for $\bar{u}$ from the IUPAC 2013 reports and for the $u_i$ from the AME2016 evaluations. This is the case if the values can be taken from one self-consistent database, in particular if one is interested in present-day values and chooses the isotopic abundances from Meija et al. (2016b); for Th, ${}^{232}$Th is the only isotope with significant abundance, and so its power per mass of isotope can be used for the bulk element. However, if data for the abundances are from different sources, additional considerations are necessary to guarantee internal consistency.
One such case are the extinct isotopes ${}^{26}$Al and ${}^{60}$Fe, for which abundances must be set. Assuming that they are the only radioactive isotopes of their respective elements, one can take the tabulated mean atomic mass as the mean of the stable isotopes, $\bar{u}_\mathrm{st}$. Then the mean of the element in presence of the radionuclide is constrained by the closure condition for the abundances $X_i$, i.e., $\sum X_i=1$, to be \begin{equation} \bar{u}=(1-X_\mathrm{rad})\bar{u}_\mathrm{st}+X_\mathrm{rad}u_\mathrm{rad}, \end{equation} where the subscript "rad" indicates the radioactive nuclide.
In the case of K, there are the two stable isotopes ${}^{39}$K and ${}^{41}$K in addition to the radioactive isotope ${}^{40}$K. The calculation for present-day proportions, which does not need further assumptions, is described further below. If $X_{40}$ is set to a different value, the default for the ratio $r_{40}=X_{40}/X_{39}$ becomes obsolete, and we make instead the additional assumption that $r_{41}=X_{41}/X_{39}$ is constant and given by the IUPAC 2013 value. Then the mean atomic mass of the stable isotopes is given by \begin{equation} \bar{u}_\mathrm{st}=\frac{u_{39}+r_{41}u_{41}}{1+r_{41}} \end{equation} and is inserted into the previous general equation.
In [ ]:
def uelem(Xiso,iso):
"""Mean element atomic mass for elements with 1 radionuclide"""
global u,Xiso_def
# get element name from isotope name
j=iso.find('-')
el=iso[0:j]
# mean of stable isotopes from (present-day) default values
if el == 'K':
ust=(u['K-39']+rK19*u['K-41'])/(1+rK19)
elif Xiso_def == 0.:
ust=u[el] # Al-26, Fe-60: extinct
return (1-Xiso)*ust+Xiso*u[iso]
For U, we assume the existence of three radioactive isotopes and use the recent reevaluation of the ratio $r_5$ of ${}^{238}$U to ${}^{235}$U by Goldmann et al. (2015), which gave $r_5=137.794$, slightly less than the value of 137.804 derived from the IUPAC abundances, which were derived excluding meteoritic materials. In order to determine the isotopic fractions, an additional assumption is necessary to construct a unique solution. The choice made here is to set the ratio $r_4$ of ${}^{238}$U to ${}^{234}$U to the value given by the IUPAC data, based on the assumption that all ${}^{234}$U, which appears only in minute amounts, always stems from ${}^{238}$U, of whose decay chain it is a part, and that the relative amounts of these two coupled isotopes are consistent with each other; this implies that $r_4$ is always the same, even if a different $r_5$ is used. With $r_4=18384.111$ and the closure condition $\sum X_i=1$, we have \begin{align} X_{238}&=\left(\frac{1}{r_4}+\frac{1}{r_5}+1\right)^{-1}\\ \bar{u}_\mathrm{U}&=X_{234}u_{234}+X_{235}u_{235}+X_{238}u_{238}= \left(\frac{u_{234}}{r_4}+\frac{u_{235}}{r_5}+u_{238}\right)X_{238}. \end{align} For the purpose of calculating bulk element power we treat ${}^{234}$U together with ${}^{238}$U.
In [ ]:
def XisoU():
"""Isotope fractions of U isotopes 234, 235, 238"""
global u
X_U238=1/(1+1/rU84+1/rU85)
X_U235=X_U238/rU85
X_U234=X_U238/rU84
u['U']=(u['U-234']/rU84+u['U-235']/rU85+u['U-238'])*X_U238
return [X_U234,X_U235,X_U238]
A conversion of energy to mass is needed in order to determine the mass of an excited nuclear isomer from the mass of the ground state. Such a conversion is needed in the calculations for the decay chain of $^{238}$U, which includes $^{234}$Pa and its isomer $^{234}$Pa$^\mathrm{m}$.
In [ ]:
def keV2u(E):
"""Conversion of energy in keV to mass in u"""
return 1e3*E*ec0/(c2*amu)
The energy available for heating of a decay to the $k$-th level of the daughter consists of the mean energy of the $\beta^-$ emission and the entire radiative (electromagnetic) energy released in the de-excitation of the $k$-th level. The latter is calculated as the difference between the total decay energy and the maximum $\beta^-$ energy rather than read from tabulated observations. The terms of the sum are weighted by the intensity/probability of decay to the $k$-th level.
In [ ]:
def getEHb(Qs,Eb_endp_i,Eb_mean_i,Ib_i):
"""Beta decay energy available for heating"""
EHb=0.
k=0
while (k < len(Ib_i)):
EHb+=(Eb_mean_i[k]+Qs-Eb_endp_i[k])*Ib_i[k]
k+=1
return EHb
As stated above, in $\alpha$ decays, the entire decay energy is converted to heat. The following function returns the $Q$ value of $\alpha$ decays. The weight factor wa is for use in a multi-branch decay, where it would be $<1$ and represents the probability of the $\alpha$ branch.
In [ ]:
def QHalpha(P,D,lmbd,wa=1):
"""Energies and power of alpha decay"""
# In alpha decays, the entire decay energy is converted to heat. In a
# multi-branch decay, wa < 1 is the probability of the alpha branch.
global Etot,Qgstot,silent
Qgs=Qa(P,D)
if (wa == 1):
Etot+=Qgs
Qgstot+=Qgs
else:
Etot+=wa*Qgs
Qgstot+=wa*Qgs
if silent == 1: return
print ("\t(lambda_s=%.4e 1/yr)\tI=%f%%\n"
"\tQ_gs=%.3f keV (%.3e J)\n" \
"\tE_H (per atom)=%.6f keV (%.3e J)") % \
(lmbd,wa*100,Qgs,1e3*Qgs*ec0,Qgs,1e3*Qgs*ec0)
In a $\beta^-$ decay, we have to determine how much energy is carried off by neutrinos, for each target (daughter isotope) energy level. In a multi-branch decay, wb $<1$ is the probability of the $\beta^-$ branch. This weight factor wb thus corresponds to the sum of all $\beta^-$ intensities in a multi-branch decay and ensures that the energies are weighted correctly in the total sum. The $\beta^-$ intensities from the original sources, which do not always sum exactly to 1 even in single-branch decays, are normalized in this function to ensure that the decays to the different target levels of the daughter have the correct relative proportions.
In [ ]:
def QHbetam(P,D,lmbd,Ib_i,Eb_endp_i,Eb_mean_i,wb=1):
"""Energies and power of beta- decay"""
global Etot,Qgstot,silent
Qgs=Qb(P,D)
if (len(Eb_endp_i) == 0): Eb_endp_i=[Qgs]
# normalize intensities to 1
Ib=sum(Ib_i)
for i in range(len(Ib_i)):
Ib_i[i]/=Ib
Eb_endp=sum(EI[0]*EI[1] for EI in zip(Eb_endp_i,Ib_i))
Eb_mean=sum(EI[0]*EI[1] for EI in zip(Eb_mean_i,Ib_i))
fr_b=Eb_mean/Eb_endp
EHb=getEHb(Qgs,Eb_endp_i,Eb_mean_i,Ib_i)
if (wb == 1):
Qgstot+=Qgs
Etot+=EHb
else:
Qgstot+=wb*Qgs
Etot+=wb*EHb
if silent == 1: return
print ("\t(lambda_s=%.4e 1/yr)\tI=%f%%\tX_b=%.6f\n"
"\tQ_gs=%.3f keV (%.3e J)\n" \
"\tE_H (per atom)=%.6f keV (%.3e J)") % \
(lmbd,wb*100,fr_b,Qgs,1e3*Qgs*ec0,EHb,1e3*EHb*ec0)
$^{26}$Al was probably the most important heat source during the early formation stages of planets and is now extinct in planetary bodies. It decays with a half-life of $\sim 717$ ky into the stable $^{26}$Mg, either by $\beta^+$ decay with a probability $I$ of 0.8173 or by ec with a probability of 0.1827 (NuDat v.2.7$\beta$, Basunia & Hurst, 2016). Almost one quarter of the total decay energy is carried off by neutrinos and does not contribute to heat production. - As the isotope is extinct, the power calculation for bulk Al returns zero by default; this can be changed by assigning a non-zero value to Xiso.
In [ ]:
# nuclides
# Aluminum-26 #################################################################
def calc_Al26(Xiso):
"""Decay energy and power of Al-26"""
global u,Xiso_def
u={'Al-26': 25.986891863,'Al': 26.9815385,'Mg-26': 25.982592971}
Xiso_def=0.
if Xiso == "":
Xiso=Xiso_def # default: extinct
ubar=u['Al']
else:
ubar=uelem(Xiso,'Al-26')
Qgs=Qb('Al-26','Mg-26') # ground-state energy difference
# intensities and energies from NuDat v.2.7b
# beta+ branch
print "beta+ decay: Al-26 -> Mg-26"
Ib=0.8173 # beta+ intensity (NuDat v.2.7b)
thalf_b=717e3 # yr
lambda_b=log(2.)/thalf_b
Eb_endp=1173.42
Eb_mean=543.29
fr_b=Eb_mean/Eb_endp
Egb=Qgs-Eb_endp # total gamma from beta+, Mg-26 exc1->gs + e-e' annihilation
EHb=Ib*(Eb_mean+Egb) # total heat from beta+
# electron capture branch
print "electron capture: Al-26 -> Mg-26"
# ec branches to Mg-26 excited states
Ie_i={'1': 0,'2': 0.0274}
Ie_i['1']=1-(Ib+Ie_i['2']) # applying closure condition
thalf_e=thalf_b
lambda_e=log(2.)/thalf_e
# gamma levels
Eg_i={'2-0': 2938, '1-0': Egb-E_ep0, '2-1': 0}
Eg_i['2-1']=Eg_i['2-0']-Eg_i['1-0']
# total heat from ec
EHe=(Eg_i['2-0']+Eg_i['2-1'])*Ie_i['2']+Eg_i['1-0']*Ie_i['1']
# total heat and power
EH=EHb+EHe
He=EtoH_W(lambda_e,EH,'Al-26')
He_Altot=He*Xiso*u['Al-26']/ubar
print (" lambda=%.4e 1/yr\tX_b=%.6f\n"
" Q_gs=%.3f keV (%.3e J)\n" \
" E_H (per atom)=%.3f keV (%.3e J)") % \
(lambda_e,fr_b,Qgs,1e3*Qgs*ec0,EH,1e3*EH*ec0)
print " H(Al-26)=%.6e W/kg\tH(Al)=%.6e W/kg" % (He,He_Altot)
if Xiso != Xiso_def: print (" Mean u_Al used with Xiso=%.3e: %.7f amu") % (Xiso,ubar)
$^{40}$K is the most abundant of the four still extant nuclides considered here and is of special interest because it is assumed to contribute most of the radiogenic heat of the core; its present-day fraction in K is $1.1668\times 10^{-4}$ (Naumenko et al., 2013). It has two important decay branches, both of which produce stable daughters: the most likely ($I=0.8928$) decay is a $\beta^-$ decay directly to the ground state of $^{40}$Ca, the other ($I=0.1072$) produces $^{40}$Ar by electron capture, as does the third, very low-probability $\beta^+$ branch (NuDat v.2.7$\beta$, Chen, 2017). As the branch leading to $^{40}$Ca decays directly into the ground state, there is no $\gamma$ contribution. Note that up to now, the mean $\beta$ energy in the $\beta^-$ decay to $^{40}$Ca was underestimated in nuclear data evaluations. The value of 583.55 keV adopted here is a recent recalculation that overcomes limitations in older evaluations related to the correct treatment of 3rd forbidden unique transitions such as this one (Mougeot, 2015 and pers.comm. 2017); it is the mean of a calculation with the full calculated $\beta$ spectrum and a calculation based on an experimental shape factor, using $Q$ as derived from AME2016.
The present-day isotopic fraction of ${}^{40}$K and the ratio $r_{40}=X_{40}/X_{39}$ were recently redetermined by Naumenko et al. (2013). The fraction was found to be the same as the four decades old one given in Meija et al. (2016b) within error, but for consistency we recalculate the mean atomic mass using the best estimates from Naumenko et al. (2013), $X_{40}=1.1668\times 10^{-4}$ and $r_{40}=1.25116\times 10^{-4}$, and the tabulated masses of the three relevant isotopes, of which two are stable: \begin{equation} \bar{u}_\mathrm{K}=\left(u_{40}-u_{41}+\frac{u_{39}-u_{41}}{r_{40}}\right)X_{40}+u_{41}. \end{equation} The result agrees with the IUPAC value if rounded to the same accuracy.
In [ ]:
# Potassium-40 ################################################################
def calc_K40(Xiso):
"""Decay energy and power of K-40"""
global u,Etot,Qgstot,Xiso_def
# mean element atomic mass of K is calculated from default isotope ratios
u={'Ar-40': 39.96238312378,\
'K-39': 38.96370648661,'K-40': 39.963998166,\
'K-41': 40.96182525796,'K': 0.,\
'Ca-40': 39.962590865}
# default: Naumenko et al. (2013)
Xiso_def=1.1668e-4 # K-40/K
r09=1.25116e-4 # K-40/K-39
if Xiso == "":
Xiso=Xiso_def
u['K']=(u['K-40']-u['K-41']+(u['K-39']-u['K-41'])/r09)*Xiso_def+u['K-41']
ubar=u['K']
else:
ubar=uelem(Xiso,'K-40')
thalf_tot=1.248e9 # yr
# branching ratio from NuDat v.2.7b
lambda_eff=log(2.)/thalf_tot
xb=0.8928
xe_i=[0.1067,0.0005]
xe=sum(xe_i)
br=xe/xb
print "beta- decay: K-40 -> Ca-40"
lambda_s=lambda_eff/(1+br)
# beta energy from Mougeot (2015 and pers.comm., 2017)
Ib_i=[1.]
Eb_endp_i=[] # same as ground-state Q because of absence of gammas
Eb_mean_i=[583.55]
QHbetam('K-40','Ca-40',lambda_s,Ib_i,Eb_endp_i,Eb_mean_i,wb=xb)
print "electron capture: K-40 -> Ar-40"
# a beta+ component with an intensity of 1e-5 is neglected here
lambda_s=lambda_s*br
Qgs=Qb('K-40','Ar-40')
Qgstot+=xe*Qgs
Ig=xe_i[0]/xe # most but not all decays entail a gamma
Eg=1460.82*Ig # NuDat v.2.7b
Etot+=xe*Eg
fr_g=Eg/Qgs
Htot=EtoH_W(lambda_eff,Etot,'K-40')
print ("\tlambda_s=%.4e 1/yr\tX_g=%.6f\n" \
"\tQ_gs=%.3f keV (%.3e J)\n" \
"\tE_H (per atom)=%.6f keV (%.3e J)") % \
(lambda_s,fr_g,Qgs,1e3*Qgs*ec0,Eg,1e3*Eg*ec0)
print "total K-40 decay (weighted mean)"
print (" Q_gs=%.3f keV (%.3e J)\n E_H (per atom)=%.6f keV (%.3e J)\n" \
" H(K-40)=%.6e W/kg\tH(K)=%.6e W/kg") % \
(Qgstot,1e3*Qgstot*ec0,Etot,1e3*Etot*ec0,Htot,Htot*Xiso*u['K-40']/ubar)
print (" Mean u_K used with Xiso=%.4e: %.4f amu") % (Xiso,ubar)
The other major radiogenic heat source during the earliest part of Solar System history may have been $^{60}$Fe ($T_{1/2}\approx 2.62$ Myr), which transmutes by two $\beta^-$ decays via $^{60}$Co ($T_{1/2}\approx 5.2712$ yr) to $^{60}$Ni (NUBASE2016) and is now extinct in planetary bodies; the initial concentrations of this nuclide are very uncertain, and recent work suggests that older studies have overestimated its abundance and hence its importance as a heat source (Boehnke et al., Trappitsch et al., 2017). As the half-life of the intermediate step is more than six orders of magnitude shorter than that of $^{60}$Fe, the $T_{1/2}$ of the latter is assumed to apply to the entire chain. However, it is the second decay step that contributes more than 95% of the total heat. - As the isotope is extinct, the power calculation for bulk Fe returns zero by default; this can be changed by assigning a non-zero value to Xiso.
In [ ]:
# Iron-60 #####################################################################
def calc_Fe60(Xiso):
"""Decay energy and power of the Fe-60 chain"""
global u,Etot,Qgstot,Xiso_def
u={'Fe-60': 59.934070411,'Fe': 55.845,\
'Co-60': 59.933815667,'Ni-60': 59.930785256}
Xiso_def=0.
if Xiso == "":
Xiso=Xiso_def # default: extinct
ubar=u['Fe']
else:
ubar=uelem(Xiso,'Fe-60')
print " Step 1, beta- decay: Fe-60 -> Co-60"
thalf_b=2.62e6 # yr
lambda_eff=log(2.)/thalf_b
# beta intensities and energies from NuDat v.2.7b
Ib_i=[1.]
Eb_endp_i=[178.]
Eb_mean_i=[50.2]
QHbetam('Fe-60','Co-60',lambda_eff,Ib_i,Eb_endp_i,Eb_mean_i)
print " Step 2, beta- decay: Co-60 -> Ni-60 (T1/2=5.2712 yr)"
thalf_s=5.2712 # yr, NUBASE2016
lambda_s=log(2.)/thalf_s
# T1/2(Co-60) << T1/2(Fe-60), hence lambda_b for Fe-60 is rate-limiting
# an intermediate level with Eb_endp=664.19 keV and I=2e-5 is neglected
# beta intensities and energies from NuDat v.2.7b
Ib_i=[0.9988,0.0012]
Eb_endp_i=[317.05,1490.29]
Eb_mean_i=[95.77,625.87]
QHbetam('Co-60','Ni-60',lambda_s,Ib_i,Eb_endp_i,Eb_mean_i)
Htot=EtoH_W(lambda_eff,Etot,'Fe-60')
print "Total heat energy/power for combined Fe-60 -> Co-60 -> Ni-60 decay"
print " Q_gs=%.3f keV (%.3e J)\n E_H (per atom)=%.6f keV (%.3e J)" % \
(Qgstot,1e3*Qgstot*ec0,Etot,1e3*Etot*ec0)
print " Htot(Fe-60)=%.6e W/kg\tHtot(Fe)=%.6e W/kg" % \
(Htot,Htot*Xiso*u['Fe-60']/ubar)
if Xiso != Xiso_def: print (" Mean u_Fe used with Xiso=%.3e: %.3f amu") % (Xiso,ubar)
Of the six radionuclides considered, $^{232}$Th is the one with the longest half-life ($T_{1/2}\approx 14$ Gy), and it is the only isotope of Th that exists today in appreciable amounts. It has a long decay chain consisting of several $\alpha$ and $\beta^-$ decays, but as the half-lives of all unstable daughters are between 9 and 18 orders of magnitude shorter, $T_{1/2}=14$ Gy is used for all rate calculations. While heat production calculation from the $\alpha$ decays is straightforward and requires only the masses of the parent and daughter, care must be taken with the $\beta^-$ decays of $^{228}$Ra (Luca, 2009 2012), $^{228}$Ac (NuDat v.2.7$\beta$, Abusaleem, 2014), $^{212}$Pb (Nichols, 2011b,c), $^{212}$Bi (Nichols, 2011a,d), and $^{208}$Tl (Nichols, 2010, 2016), all of which feed several different excited energy levels of their respective daughters. Moreover, $^{212}$Bi has two decay branches: it can decay via the $\alpha$ branch ($I=0.3593$) to $^{208}$Tl or via the $\beta^-$ branch ($I=0.6407$) to $^{212}$Po (Nichols, 2011d), and the decay energies of the daughters have to be weighted accordingly. Note that in some cases the $\beta^-$ branch of $^{212}$Bi leads to a $^{212}$Po$^\mathrm{m}$ state that proceeds directly via $\alpha$ decay to $^{208}$Pb, i.e., without passing through the $^{212}$Po ground state by $\gamma$ emission; this branch is not accounted for explicitly in the energy calculation, but it is in fact included in the final sum, because it is energetically equivalent to the path through the ground state.
Apart from the final stable daughter $^{208}$Pb, each decay of $^{232}$Th eventually also produces 6 atoms of $^{4}$He.
In [ ]:
# Thorium-232 #################################################################
def calc_Th232():
"""Decay energy and power of the Th-232 series"""
global u,Etot,Qgstot
u={'Tl-208': 207.982017992,\
'Pb-208': 207.976651918,'Pb-212': 211.991895975,\
'Bi-212': 211.991285016,\
'Po-212': 211.988867896,'Po-216': 216.001913506,\
'Rn-220': 220.011392534,\
'Ra-224': 224.020210453,'Ra-228': 228.031068657,\
'Ac-228': 228.031019767,\
'Th-228': 228.028739835,'Th-232': 232.038053689}
print "decay chain: Th-232 -> Pb-208"
# Th-232 -> Ra-228
print " Step 1, alpha decay: Th-232 -> Ra-228"
thalf_Th232=14e9 # yr, AME2016
lambda_eff=log(2.)/thalf_Th232
QHalpha('Th-232','Ra-228',lambda_eff)
# Ra-228 -> Ac-228
print " Step 2, beta- decay: Ra-228 -> Ac-228"
# T1/2(daughters) << T1/2(Th-232), hence lambda_Th232 is rate-limiting for all
thalf_s=5.75 # yr
lambda_s=log(2.)/thalf_s
# beta energies from NuDat v.2.7b, intensities from Luca (2009a, 2012)
Ib_i=[0.3,0.087,0.49,0.12]
Eb_endp_i=[12.7,25.6,39.1,39.5]
Eb_mean_i=[3.22,6.48,9.94,10.04]
QHbetam('Ra-228','Ac-228',lambda_s,Ib_i,Eb_endp_i,Eb_mean_i)
# Ac-228 -> Th-228
print " Step 3, beta- decay: Ac-228 -> Th-228"
thalf_s=6.15/yr2h # conversion h->yr, NUBASE2016
lambda_s=log(2.)/thalf_s
# beta intensities and energies from NuDat v.2.7b
Ib_i=[4.2e-05,6.2e-05,0.00019,0.00054,2.9e-05,0.0028,0.00038,3.3e-05, \
0.00251,0.00046,0.00056,0.00034,0.00077,0.00117,0.00045,5e-05, \
0.00119,0.00062,0.00389,0.00133,0.0176,0.0243,0.00164,0.0112, \
0.0419,0.03,0.0115,0.00095,0.076,0.006,0.012,0.0006,0.00217, \
0.00058,0.0067,0.0017,2.5e-05,0.0311,0.058,0.059,0.0026,0.0007, \
0.0311,0.0033,0.0014,0.299,0.00041,0.0021,0.,0.0059,0.1165, \
0.00116,0.006,0.07]
Eb_endp_i=[11.,97.,104.,111.,120.,124.,147.,175.,189.,197.,205.,227.,234., \
241.,336.,338.,374.,376.,390.,398.,410.,446.,450.,451., \
488.,491.,496.,516.,603.,684.,702.,718.,790.,837.,907., \
959.,960.,966.,981.,1011.,1043.,1074.,1111.,1118.,1154., \
1165.,1190.,1260.,1615.,1806.,1738.,1756.,1947.,2076.]
Eb_mean_i=[2.73,25.15,26.94,29.04,31.51,32.44,38.72,46.79,50.75,52.94,
55.39,61.79,63.71,65.77,94.51,95.13,106.17,106.77,111.29, \
113.94,117.50,129.00,130.51,130.82,142.79,143.73,145.32, \
152.1,181.0,208.9,215.3,220.9,246.4,260.3,288.8,307.6, \
307.9,310.2,315.7,327.0,338.9,350.6,364.7,367.0,381.0, \
385.0,394.4,421.2,536.9,609.3,609.7,616.9,694.3,747.0]
QHbetam('Ac-228','Th-228',lambda_s,Ib_i,Eb_endp_i,Eb_mean_i)
# Th-228 -> Ra-224
print " Step 4, alpha decay: Th-228 -> Ra-224"
thalf_s=1.9124 # yr, NUBASE2016
lambda_s=log(2.)/thalf_s
QHalpha('Th-228','Ra-224',lambda_s)
# Ra-224 -> Rn-220
print " Step 5, alpha decay: Ra-224 -> Rn-220"
thalf_s=3.6319/yr2d # conversion d->yr, NUBASE2012
lambda_s=log(2.)/thalf_s
QHalpha('Ra-224','Rn-220',lambda_s)
# Rn-220 -> Po-216
print " Step 6, alpha decay: Rn-220 -> Po-216"
thalf_s=55.6/yr2s # conversion s->yr, NUBASE2016
lambda_s=log(2.)/thalf_s
QHalpha('Rn-220','Po-216',lambda_s)
# Po-216 -> Pb-212
print " Step 7, alpha decay: Po-216 -> Pb-212"
thalf_s=0.145/yr2s # conversion s->yr, NUBASE2016
lambda_s=log(2.)/thalf_s
QHalpha('Po-216','Pb-212',lambda_s)
# Pb-212 -> Bi-212
print " Step 8, beta- decay: Pb-212 -> Bi-212"
thalf_s=10.64/yr2h # conversion h->yr, NUBASE2016
lambda_s=log(2.)/thalf_s
# beta intensities and energies from Nichols (2011)
Ib_i=[0.0499,0.817,0.133]
Eb_endp_i=[154.6,331.3,569.9]
Eb_mean_i=[41.1,93.5,171.7]
QHbetam('Pb-212','Bi-212',lambda_s,Ib_i,Eb_endp_i,Eb_mean_i)
# two decay branches for Bi-212
# branching ratio from Nichols (2011)
thalf_tot=60.55/yr2m # min (NUBASE2016)
lambda_tot=log(2.)/thalf_tot
xb=0.6407
xa=0.3593
br=xa/xb
# Bi-212 -> Po-212
# Some Po-212m levels decay directly to Pb-208; these are not treated separately
# but are covered equivalently by the gamma emission to Po-212 and the subsequent alpha
# decay Po-212 -> Pb-208 (step 10a below).
print " Step 9a, beta- decay: Bi-212 -> Po-212"
lambda_s=lambda_tot/(1+br)
# beta intensities and energies from Nichols (2011)
Ib_i=[0.0068,0.00032,0.0021,0.019,0.0144,0.045,0.5531]
Eb_endp_i=[446.1,451.2,572.7,631.4,739.4,1524.8,2252.1]
Eb_mean_i=[130.1,131.7,172.4,192.7,230.8,533.1,834.2]
QHbetam('Bi-212','Po-212',lambda_s,Ib_i,Eb_endp_i,Eb_mean_i,wb=xb)
# Bi-212 -> Tl-208
print " Step 9b, alpha decay: Bi-212 -> Tl-208"
lambda_s=lambda_s*br
QHalpha('Bi-212','Tl-208',lambda_s,wa=xa)
# Po-212 -> Pb-208
print " Step 10a, alpha decay: Po-212 -> Pb-208"
thalf_s=294.7e-9/yr2s # s, NUBASE2016
lambda_s=log(2.)/thalf_s
# weight xb is used b/c this continues the original beta- branch
QHalpha('Po-212','Pb-208',lambda_s,wa=xb)
# Tl-208 -> Pb-208
print " Step 10b, beta- decay: Tl-208 -> Pb-208"
thalf_s=3.053/yr2m # min (NUBASE2016)
lambda_s=log(2.)/thalf_s
# beta intensities and energies from Nichols (2016)
Ib_i=[0.00052,0.00017,0.00045,5e-05,0.00102,2e-05,0.00231,0.00174,7e-05, \
0.0317,0.00048,0.0063,0.241,0.221,0.492]
Eb_endp_i=[518.3,615.7,640.3,675.1,702.4,737.1,818.6,873.7,1003.6,1037.8, \
1052.4,1079.0,1290.5,1523.9,1801.3]
Eb_mean_i=[154.3,187.7,196.4,208.6,218.3,230.8,260.4,280.8,329.7,342.8, \
348.4,358.6,441.5,535.4,649.5]
# weight xa is used b/c this continues the original alpha branch
QHbetam('Tl-208','Pb-208',lambda_s,Ib_i,Eb_endp_i,Eb_mean_i,wb=xa)
print "Total heat energy/power for Th-232 decay series"
print " Q_gs=%.3f keV (%.3e J)\n E_H (per atom)=%.6f keV (%.3e J)" % \
(Qgstot,1e3*Qgstot*ec0,Etot,1e3*Etot*ec0)
print " Htot=%.6e W/kg" % EtoH_W(lambda_eff,Etot,'Th-232')
$^{235}$U is the less abundant and shorter-lived ($T_{1/2}=704$ My; Huang & Wang, 2014) of the two major uranium isotopes. Its decay chain is even longer and more complicated than that of $^{232}$Th and yields $^{207}$Pb as well as 7 atoms of $^{4}$He as final products. The decay branch from $^{223}$Fr to $^{219}$At has a very low probability, and therefore the subsequent steps of that branch are omitted from the calculation; they merge with the main decay path at different points. With regard to the energy balance, the $\beta^-$ decays of $^{231}$Th (NuDat v.2.7$\beta$, Browne & Tuli, 2013a), $^{227}$Ac (NuDat v.2.7$\beta$, Kondev et al., 2016), $^{223}$Fr (Huang & Wang, 2009, 2012), $^{211}$Pb (NuDat v.2.7$\beta$, Singh et al., 2013), and $^{207}$Tl (NuDat v.2.7$\beta$, Kondev & Lalkovski, 2011) require attention because of the neutrino component. A possible, very low-probability $\beta^-$ decay path from $^{207}$Tl$^\mathrm{m}$ to $^{207}$Pb is not treated individually but as if it were identical to the path via the $^{207}$Tl ground state.
In [ ]:
# Uranium-235 #################################################################
def calc_U235():
"""Decay energy and power of the U-235 series"""
global u,Etot,Qgstot,HtotU,silent
# mean element atomic mass of U is calculated from isotope ratios
u={'Tl-207': 206.977418586,\
'Pb-207': 206.975896735,'Pb-211': 210.988735356,\
'Bi-211': 210.987268698,
'Po-211': 210.986653085,'Po-215': 214.999418454,\
'At-219': 219.011160647,\
'Rn-219': 219.009478753,\
'Fr-223': 223.019734313,\
'Ra-223': 223.018500719,\
'Ac-227': 227.027750666,\
'Th-227': 227.027702618,'Th-231': 231.036302853,\
'Pa-231': 231.035882575,\
'U-234': 234.040950370,'U-235': 235.043928190,\
'U-238': 238.050786996,'U': 0}
if silent == 0: print "decay chain: U-235 -> Pb-207"
# U-235 -> Th-231
if silent == 0: print " Step 1, alpha decay: U-235 -> Th-231"
thalf_U235=704e6 # yr (NUBASE2016)
lambda_eff=log(2.)/thalf_U235
QHalpha('U-235','Th-231',lambda_eff)
# Th-231 -> Pa-231
if silent == 0: print " Step 2, beta- decay: Th-231 -> Pa-231"
# T1/2(daughters) << T1/2(U-235), hence lambda_U235 is rate-limiting for all
thalf_s=25.52/yr2h # h (NUBASE2016)
lambda_s=log(2.)/thalf_s
# beta intensities and energies from NuDat v.2.7b
Ib_i=[3e-5,0.00065,7.8e-6,0.026,0.003,0.121,0.013,0.12,0.4,0.32,0.0017,\
0.0017,0.00022]
Eb_endp_i=[39.8,71.4,73.6,144.3,173.4,208.1,217.4,289.3,290.2,307.4,313.9,\
333.0,391.6]
Eb_mean_i=[9.8,18.1,18.7,37.8,46.0,55.9,58.6,79.8,80.1,85.3,87.3,93.1,111.4]
QHbetam('Th-231','Pa-231',lambda_s,Ib_i,Eb_endp_i,Eb_mean_i)
# Pa-231 -> Ac-227
if silent == 0: print " Step 3, alpha decay: Pa-231 -> Ac-227"
thalf_s=3.276e4 # yr (NUBASE2016)
lambda_s=log(2.)/thalf_s
QHalpha('Pa-231','Ac-227',lambda_eff)
# two decay branches for Ac-227 (Kondev et al., 2016)
thalf_tot=21.772 # yr (NUBASE2016)
lambda_tot=log(2.)/thalf_tot
xb=0.9862
xa=0.0138
br=xa/xb
# Ac-227 -> Th-227
if silent == 0: print " Step 4a, beta- decay: Ac-227 -> Th-227"
lambda_s=lambda_tot/(1+br)
# beta intensities and energies from NuDat v.2.7b
Ib_i=[0.002999990814,0.0999996938,0.3499989283,0.53999834652]
Eb_endp_i=[6.9,20.3,35.5,44.8]
Eb_mean_i=[1.73,5.11,8.98,11.37]
QHbetam('Ac-227','Th-227',lambda_s,Ib_i,Eb_endp_i,Eb_mean_i,wb=xb)
# Ac-227 -> Fr-223
if silent == 0: print " Step 4b, alpha decay: Ac-227 -> Fr-223"
lambda_s=lambda_s*br
QHalpha('Ac-227','Fr-223',lambda_eff,wa=xa)
# Th-227 -> Ra-223
if silent == 0: print " Step 5a, alpha decay: Th-227 -> Ra-223"
lambda_s=18.697/yr2d # d, NUBASE2016
# weight xb is used b/c this continues the original beta- branch
QHalpha('Th-227','Ra-223',lambda_eff,wa=xb)
# two decay branches for Fr-223 (Huang & Wang, 2012)
thalf_tot=22./yr2m # min, NUBASE2016
lambda_tot=log(2.)/thalf_tot
xb_a=0.9998
xa_a=0.0002
br=xa_a/xb_a
# Fr-223 -> Ra-223
if silent == 0: print " Step 5b-1, beta- decay: Fr-223 -> Ra-223"
lambda_s=lambda_tot/(1+br)
# beta intensities and energies from Huang & Wang (2012)
Ib_i=[0.000012,0.000004,0.0000046,0.0002,0.000082,0.000051,0.00106,\
0.000011,0.00025,0.00088,0.00035,0.0054,0.00014,0.00004,0.0014,\
0.00019,0.0000111,0.00013,0.000046,0.018,0.00037,0.00042,0.00049,\
0.00032,0.00004,0.091,0.0024,0.15,0.0027,0.67,0.06,0.01]
Eb_endp_i=[120.3,124.6,129.9,191.5,205.9,208.4,222.6,243.3,281.9,302.8,\
306.9,323.3,326.0,343.8,345.4,362.1,366.7,555.3,773.1,\
779.9,806.7,814.9,819.4,863.1,869.0,914.5,1025.5,1069.6,\
1087.8,1099.1,1119.3,1149.2]
Eb_mean_i=[31.5,32.7,34.1,51.5,55.6,56.3,60.5,66.6,78.1,84.4,85.7,90.7,\
91.5,97.0,97.5,102.7,104.1,165.6,241.3,243.7,253.3,256.3,\
257.9,273.8,275.9,292.6,333.9,350.5,357.4,361.7,369.4,380.8]
# continuation of orig. alpha branch, beta- fork
QHbetam('Fr-223','Ra-223',lambda_s,Ib_i,Eb_endp_i,Eb_mean_i,wb=xa*xb_a)
# Fr-223 -> At-219
if silent == 0: print " Step 5b-2, alpha decay: Fr-223 -> At-219"
lambda_s=lambda_s*br
# continuation of orig. alpha branch, alpha fork
QHalpha('Fr-223','At-219',lambda_eff,wa=xa*xa_a)
# Ra-223 -> Rn-219
if silent == 0: print " Step 6, alpha decay: Ra-223 -> Rn-219"
thalf_s=11.4377/yr2d # d (NUBASE2016)
lambda_s=log(2.)/thalf_s
QHalpha('Ra-223','Rn-219',lambda_eff)
# Rn-219 -> Po-215
if silent == 0: print " Step 7, alpha decay: Rn-219 -> Po-215"
thalf_s=3.96/yr2s # s (NUBASE2016)
lambda_s=log(2.)/thalf_s
QHalpha('Rn-219','Po-215',lambda_eff)
# Po-215 -> Pb-211
# a very low-probability beta- decay Po-215 -> At-215 is omitted
if silent == 0: print " Step 8, alpha decay: Po-215 -> Pb-211"
thalf_s=1.781e-3/yr2s # s (NUBASE2016)
lambda_s=log(2.)/thalf_s
QHalpha('Po-215','Pb-211',lambda_eff)
# Pb-211 -> Bi-211
if silent == 0: print " Step 9, beta- decay: Pb-211 -> Bi-211"
thalf_s=36.164/yr2m # min (NUBASE2016)
lambda_s=log(2.)/thalf_s
# beta intensities and energies from NuDat v.2.7b
Ib_i=[0.000186,1.3e-05,0.00017,0.0083,4.6e-05,0.00056,0.00022,0.0628,\
0.0163,0.9132]
Eb_endp_i=[96.,133.,171.,258.,264.,287.,416.,535.,962.,1367.]
Eb_mean_i=[25.3,35.3,45.9,71.3,73.1,80.2,120.7,160.2,313.7,471.3]
QHbetam('Pb-211','Bi-211',lambda_s,Ib_i,Eb_endp_i,Eb_mean_i)
# two decay branches for Bi-211
# branching ratio from NuDat v.2.7b
thalf_tot=2.14/yr2m # min, NUBASE2016
lambda_tot=log(2.)/thalf_tot
xb=0.00276
xa=0.99724
br=xa/xb
# Bi-211 -> Po-211
if silent == 0: print " Step 10a, beta- decay: Bi-211 -> Po-211"
lambda_s=lambda_tot/(1+br)
# beta intensity and energy from NuDat v.2.7b
Ib_i=[0.00276]
Eb_endp_i=[574.]
Eb_mean_i=[172.9]
QHbetam('Bi-211','Po-211',lambda_s,Ib_i,Eb_endp_i,Eb_mean_i,wb=xb)
# Bi-211 -> Tl-207
# An intermediate state Tl-207m may decay with very low probability via beta decay into
# Pb-207; this branch is not treated as a beta decay here but lumped together with the
# decay path through the ground state of Tl-207.
if silent == 0: print " Step 10b, alpha decay: Bi-211 -> Tl-207"
lambda_s=lambda_s*br
QHalpha('Bi-211','Tl-207',lambda_eff,wa=xa)
# Po-211 -> Pb-207
if silent == 0: print " Step 11a, alpha decay: Po-211 -> Pb-207"
thalf_s=516e-3/yr2s # s (NUBASE2016)
lambda_s=log(2.)/thalf_s
# weight xb is used b/c this continues the original beta- branch
QHalpha('Po-211','Pb-207',lambda_eff,wa=xb)
# Tl-207 -> Pb-207
if silent == 0: print " Step 11b, beta- decay: Tl-207 -> Pb-207"
thalf_s=4.77/yr2m # min (NUBASE2016)
lambda_s=log(2.)/thalf_s
# beta intensities and energies from NuDat v.2.7b; one level with very low
# probability omitted
Ib_i=[0.00271,0.99729]
Eb_endp_i=[520.,1418.]
Eb_mean_i=[155.0,492.5]
# weight xa is used b/c this continues the original alpha branch
QHbetam('Tl-207','Pb-207',lambda_s,Ib_i,Eb_endp_i,Eb_mean_i,wb=xa)
Htot=EtoH_W(lambda_eff,Etot,'U-235')
X_Ui=XisoU()
if silent == 0:
print "Total heat energy/power for U-235 decay series"
print " Q_gs=%.3f keV (%.3e J)\n E_H (per atom)=%.6f keV (%.3e J)" % \
(Qgstot,1e3*Qgstot*ec0,Etot,1e3*Etot*ec0)
HtotU+=X_Ui[1]*Htot*u['U-235']/u['U']
print " Htot(U-235)=%.6e W/kg\tHtot(U)=%.6e W/kg" % (Htot,HtotU)
print " Isotope fractions: X_U234=%.5f%% X_U235=%.5f%% X_U238=%.5f%%" % \
(X_Ui[0]*100,X_Ui[1]*100,X_Ui[2]*100)
print " Mean u_U: %.5f amu" % u['U']
else:
HtotU=X_Ui[1]*Htot*u['U-235']/u['U']
The other, much more abundant isotope $^{238}$U has a half-life close to the age of the Solar System, $T_{1/2}=4468$ My, and a similarly long and complex decay chain. It yields $^{206}$Pb as well as 8 atoms of $^{4}$He as final products. The isomeric transition (IT) of $^{234}$Pa was treated separately because of the slightly different $\beta^-$ energies, but the very low-probability $\beta^-$ branch of $^{218}$At and two similar branches of $^{210}$Pb to $^{206}$Hg and of $^{210}$Bi to $^{206}$Tl are not included. Energy is lost by neutrinos in the $\beta^-$ decays of $^{234}$Th (Luca, 2009b, 2010), $^{234}$Pa$^\mathrm{m}$ and $^{234}$Pa (Huang & Wang, 2011a-d), $^{218}$Po (Chisté & Bé, 2007a, 2010c), $^{214}$Pb (Chisté & Bé, 2010a,b), $^{214}$Bi (NuDat v.2.7$\beta$, Wu, 2009), $^{210}$Tl, $^{210}$Pb (NuDat v.2.7$\beta$, Basunia, 2014), and $^{210}$Bi (Chisté et al., 2014). The mass of the excited isomer $^{234}$Pa$^\mathrm{m}$ is determined from the mass of the ground state of $^{234}$Pa and the energy of the $\gamma$ photon emitted upon decay to the ground level, which is determined to be 73.92 keV (Luca, 2010); this calculation is done on the fly with the function keV2u.
In [ ]:
# Uranium-238 #################################################################
def calc_U238():
"""Decay energy and power of the U-238 series"""
# Pa-234m is a long-lived Pa-234 excited isomeric state
# mean element atomic mass of U is calculated from isotope ratios
global u,Etot,Qgstot,HtotU,silent
u={'Tl-210': 209.990072970,\
'Pb-206': 205.974465124,'Pb-210': 209.984188301,\
'Pb-214': 213.999803788,\
'Bi-210': 209.984120156,'Bi-214': 213.998710938,\
'Po-210': 209.982873601,'Po-214': 213.995201208,\
'Po-218': 218.008971502,\
'At-218': 218.008693735,\
'Rn-222': 222.017576286,\
'Ra-226': 226.025408455,\
'Th-230': 230.033132358,'Th-234': 234.043599860,\
'Pa-234': 234.043305615,'Pa-234m': 0,\
'U-234': 234.040950370,'U-235': 235.043928190,\
'U-238': 238.050786996,'U': 0}
if silent == 0: print "decay chain: U-238 -> Pb-206"
# U-238 -> Th-234
if silent == 0: print " Step 1, alpha decay: U-238 -> Th-234"
thalf_U238=4.468e9 # yr (NUBASE2016)
lambda_eff=log(2.)/thalf_U238
QHalpha('U-238','Th-234',lambda_eff)
# Th-234 -> Pa-234m
if silent == 0: print " Step 2, beta- decay: Th-234 -> Pa-234m"
# T1/2(daughters) << T1/2(U-238), hence lambda_U238 is rate-limiting for all
thalf_s=24.1/yr2d # d (NUBASE2016)
lambda_s=log(2.)/thalf_s
# beta intensities and energies from Luca (2010)
Ib_i=[0.016,0.00016,0.065,0.141,0.778]
Eb_endp_i=[85.,95.,105.,106.,198.]
Eb_mean_i=[22.,25.,27.,28.,53.]
# mass of excited isomer from ground state mass and excitation energy
Eexc=73.92 # keV
u['Pa-234m']=u['Pa-234']+keV2u(Eexc)
QHbetam('Th-234','Pa-234m',lambda_s,Ib_i,Eb_endp_i,Eb_mean_i)
# two decay branches for Pa-234m (Huang & Wang, 2011)
thalf_tot=1.159/yr2m # minutes
lambda_tot=log(2.)/thalf_tot
xb=0.9985
xIT=0.0015
br=xIT/xb
# Pa-234m -> U-234
if silent == 0: print " Step 3a, beta- decay: Pa-234m -> U-234"
lambda_s=lambda_tot/(1+br)
# beta intensities and energies from Huang & Wang (2011)
Ib_i=[0.0000389,0.000108,0.000452,0.000258,0.0000311,0.000146,0.000021,\
0.000357,0.000024,0.000061,0.0000127,0.000249,0.0000231,\
0.00032,0.000131,0.000092,0.000121,0.000046,0.01006,0.00945,\
0.00049,0.97599]
Eb_endp_i=[299.,332.,358.,394.,406.,460.,473.,488.,575.,602.,667.,677.,\
698.,715.,768.,834.,1032.,1095.,1224.,1459.,1483.,2269.]
Eb_mean_i=[83.0,93.0,101.0,112.3,116.0,133.3,137.4,142.3,171.2,180.1,\
202.5,205.8,213.3,219.2,237.6,261.1,333.1,356.7,405.6,\
496.0,505.3,820.5]
QHbetam('Pa-234m','U-234',lambda_s,Ib_i,Eb_endp_i,Eb_mean_i,wb=xb)
# Pa-234m -> Pa-234 -> U-234
if silent == 0:
print " Step 3b, isomeric transition, then beta- decay: Pa-234m -> Pa-234 -> U-234"
lambda_s=lambda_s*br
Qgs=Eexc
Qgstot+=xIT*Qgs
Etot+=xIT*Qgs
# beta intensities and energies from Huang & Wang (2011)
Ib_i=[0.0042,0.0021,0.00064,0.004,0.0014,0.00055,0.009,0.00112,0.00122,\
0.0059,0.00044,0.0044,0.0035,0.0022,0.0021,0.0025,0.00029,\
0.0017,0.0143,0.0041,0.00061,0.08,0.00129,0.028,0.0078,0.0116,\
0.36,0.084,0.069,0.0095,0.0018,0.00035,0.007,0.0005,0.196,\
0.00078,0.001,0.009,0.0021,0.0025,0.027,0.0012,0.0011,0.00109,\
0.003,0.015,0.019,0.0069,0.08,0.015,0.05,0.031,0.025,0.004,\
0.008,0.008,0.05]
Eb_endp_i=[51.,79.,94.,126.,129.,158.,161.,175.,195.,214.,226.,236.,254.,\
267.,279.,313.,332.,351.,383.,402.,411.,412.,424.,433.,\
457.,458.,472.,472.,502.,542.,545.,576.,606.,613.,642.,\
647.,651.,658.,662.,693.,699.,709.,747.,883.,980.,1000.,\
1067.,1104.,1126.,1171.,1171.2,1206.,1227.,1232.,1247.,\
1346.,2052.]
Eb_mean_i=[13.0,20.4,24.2,33.1,33.8,41.9,42.9,46.7,52.2,57.8,61.3,64.3,\
69.7,73.5,76.9,87.3,93.0,98.9,108.9,114.8,117.6,118.1,\
121.8,124.7,132.3,132.5,137.1,137.2,146.8,160.1,164.6,\
171.4,181.7,184.1,194.0,195.6,197.1,199.3,200.6,211.3,\
213.5,216.9,230.3,278.7,314.2,312.6,346.5,360.2,368.3,\
385.4,385.4,398.5,406.4,408.7,414.4,452.1,732.2]
QHbetam('Pa-234','U-234',lambda_s,Ib_i,Eb_endp_i,Eb_mean_i,wb=xIT)
# U-234 -> Th-230
if silent == 0: print " Step 4, alpha decay: U-234 -> Th-230"
# save total heat energy released up to here for later use with U-234 power
dEH234=Etot
thalf_s=245.5e3 # yr (NUBASE2016)
lambda_s=log(2.)/thalf_s
QHalpha('U-234','Th-230',lambda_s)
# Th-230 -> Ra-226
if silent == 0: print " Step 5, alpha decay: Th-230 -> Ra-226"
thalf_s=75.4e3 # yr (NUBASE2016)
lambda_s=log(2.)/thalf_s
QHalpha('Th-230','Ra-226',lambda_s)
# Ra-226 -> Rn-222
if silent == 0: print " Step 6, alpha decay: Ra-226 -> Rn-222"
thalf_s=1600. # yr (NUBASE2016)
lambda_s=log(2.)/thalf_s
QHalpha('Ra-226','Rn-222',lambda_s)
# Rn-222 -> Po-218
if silent == 0: print " Step 7, alpha decay: Rn-222 -> Po-218"
thalf_s=3.8215/yr2d # d (NUBASE2016)
lambda_s=log(2.)/thalf_s
QHalpha('Rn-222','Po-218',lambda_s)
# two decay branches for Po-218
# branching ratio from Chisté & Bé (2010)
thalf_tot=3.098/yr2m # min (NUBASE2016)
lambda_tot=log(2.)/thalf_tot
xb=0.00022
xa=0.99978
br=xa/xb
# Po-218 -> At-218
if silent == 0: print " Step 8a, beta- decay: Po-218 -> At-218"
lambda_s=lambda_tot/(1+br)
# beta intensities and energies from Chisté & Bé (2010)
Ib_i=[0.00022]
Eb_endp_i=[] # same as ground-state Q because of absence of gammas
Eb_mean_i=[73.]
QHbetam('Po-218','At-218',lambda_s,Ib_i,Eb_endp_i,Eb_mean_i,wb=xb)
# Po-218 -> Pb-214
if silent == 0: print " Step 8b, alpha decay: Po-218 -> Pb-214"
lambda_s=lambda_s*br
QHalpha('Po-218','Pb-214',lambda_s,wa=xa)
# At-218 -> Bi-214
# a low-probability beta- branch to Rn-218 is neglected
if silent == 0: print " Step 9a, alpha decay: At-218 -> Bi-214"
lambda_s=1.5/yr2s # s (NUBASE2016)
# weight xb is used b/c this continues the original beta- branch
QHalpha('At-218','Bi-214',lambda_s,wa=xb)
# Pb-214 -> Bi-214
if silent == 0: print " Step 9b, beta- decay: Pb-214 -> Bi-214"
lambda_s=27.06/yr2m # min (NUBASE2016)
# beta intensities and energies from Chisté & Bé (2010)
Ib_i=[0.02762,0.000196,0.01047,0.4652,0.4109,0.092]
Eb_endp_i=[180.,222.,485.,667.,724.,1019.]
Eb_mean_i=[50.,62.,145.,207.,227.,337.]
# weight xa is used b/c this continues the original alpha branch
QHbetam('Pb-214','Bi-214',lambda_s,Ib_i,Eb_endp_i,Eb_mean_i,wb=xa)
# two decay branches for Bi-214
# branching ratio from NuDat v.2.7b
thalf_tot=19.9/yr2m # min (NUBASE2016)
lambda_tot=log(2.)/thalf_tot
xb=0.99979
xa=0.00021
br=xa/xb
# Bi-214 -> Po-214
if silent == 0: print " Step 10a, beta- decay: Bi-214 -> Po-214"
lambda_s=lambda_tot/(1+br)
# beta intensities and energies from NuDat v.2.7b
Ib_i=[1e-05,1.4e-06,7.9e-06,9e-07,2.1e-05,7e-06,7.3e-05,0.00038,0.00043,\
7e-05,0.000104,4e-05,0.000168,4e-06,0.00041,2.17e-05,1.1e-05,\
0.0002225,1.4e-05,4.5e-05,0.000108,0.000114,0.00014,0.00032,\
0.00033,0.00047,0.0005,0.00538,0.00267,0.00104,0.00292,0.00127,\
0.00036,0.0015,1.8e-06,1e-06,0.00044,0.00123,0.00173,0.01244,\
0.0278,0.00074,6e-05,0.00557,0.00196,0.056,0.00855,0.00444,\
0.04345,0.00102,0.0245,0.01431,0.01177,0.01588,0.0814,0.1696,\
0.00163,0.00125,0.1757,0.0061,0.0312,0.0089,0.0735,0.0058,0.191]
Eb_endp_i=[86.,97.,110.,121.,127.,176.,188.,216.,256.,267.,270.,284.,291.,\
307.,329.,335.,341.,348.,351.,373.,376.,390.,400.,409.,\
443.,484.,500.,541.,551.,571.,575.,608.,639.,665.,708.,\
717.,725.,762.,765.,788.,822.,847.,909.,977.,1004.,1066.,\
1077.,1122.,1151.,1182.,1253.,1259.,1275.,1380.,1423.,\
1505.,1557.,1527.,1540.,1609.,1727.,1855.,1892.,2661.,3270.]
Eb_mean_i=[22.3,25.1,28.6,31.7,33.5,47.1,50.6,58.7,70.5,73.7,74.7,78.9,\
81.1,86.0,92.8,94.8,96.6,98.7,99.4,106.5,107.6,111.8,115.2,\
118.1,129.1,142.5,147.9,161.8,164.9,171.9,173.3,184.4,\
195.4,204.5,219.5,222.8,225.7,238.9,239.9,248.2,260.9,\
269.8,293.0,318.2,328.4,352.1,356.5,373.7,385.1,396.8,\
424.6,427.1,433.5,474.9,492.0,525.3,529.7,534.0,539.4,\
567.2,615.4,668.1,683.7,1007.5,1268.8]
QHbetam('Bi-214','Po-214',lambda_s,Ib_i,Eb_endp_i,Eb_mean_i,wb=xb)
# Bi-214 -> Tl-210
if silent == 0: print " Step 10b, alpha decay: Bi-214 -> Tl-210"
lambda_s=lambda_s*br
QHalpha('Bi-214','Tl-210',lambda_s,wa=xa)
# Po-214 -> Pb-210
if silent == 0: print " Step 11a, alpha decay: Po-214 -> Pb-210"
lambda_s=163.72e-6/yr2s # s (NUBASE2016)
# weight xb is used b/c this continues the original beta- branch
QHalpha('Po-214','Pb-210',lambda_s,wa=xb)
# Tl-210 -> Pb-210
if silent == 0: print " Step 11b, beta- decay: Tl-210 -> Pb-210"
lambda_s=1.3/yr2m # min (NUBASE2016)
# beta intensities and energies from NuDat v.2.7b
Ib_i=[0.02,0.07,0.24,0.1,0.1,0.3,0.2]
Eb_endp_i=[1380.,1600.,1860.,2020.,2413.,4210.,4386.]
Eb_mean_i=[477.,568.,674.,743.,877.1,1635.,1762.6]
# weight xa is used b/c this continues the original alpha branch
QHbetam('Tl-210','Pb-210',lambda_s,Ib_i,Eb_endp_i,Eb_mean_i,wb=xa)
# Pb-210 -> Bi-210
# a very weak alpha decay branch to Hg-206 is neglected
if silent == 0: print " Step 12, beta- decay: Pb-210 -> Bi-210"
lambda_s=22.2 # yr (NUBASE2016)
# beta intensities and energies from NuDat v.2.7b
Ib_i=[0.84,0.16]
Eb_endp_i=[17.,63.5]
Eb_mean_i=[4.16,16.16]
QHbetam('Pb-210','Bi-210',lambda_s,Ib_i,Eb_endp_i,Eb_mean_i)
# Bi-210 -> Po-210
# a very weak alpha decay branch to Tl-206 is neglected
if silent == 0: print " Step 13, beta- decay: Bi-210 -> Po-210"
lambda_s=5.012/yr2d # d (NUBASE2016)
# beta intensities and energies from Chisté et al. (2014)
Ib_i=[1.]
Eb_endp_i=[] # same as ground-state Q because of absence of gammas
Eb_mean_i=[317.]
QHbetam('Bi-210','Po-210',lambda_s,Ib_i,Eb_endp_i,Eb_mean_i)
# Po-210 -> Pb-206
if silent == 0: print " Step 14, alpha decay: Po-210 -> Pb-206"
thalf_s=138.376/yr2d # d (NUBASE2016)
lambda_s=log(2.)/thalf_s
QHalpha('Po-210','Pb-206',lambda_s)
Htot=EtoH_W(lambda_eff,Etot,'U-238')
# H(U-234) uses lambda_eff, because supply from U-238 is still the rate-limiting factor
H234=EtoH_W(lambda_eff,Etot-dEH234,'U-234')
X_Ui=XisoU()
if silent == 0:
print "Total heat energy/power for U-238 decay series"
print " Q_gs=%.3f keV (%.3e J)\n E_H (per atom)=%.6f keV (%.3e J)" % \
(Qgstot,1e3*Qgstot*ec0,Etot,1e3*Etot*ec0)
# weighted mean with U-234, which releases a bit less heat
HtotU+=(X_Ui[0]*H234*u['U-234']+X_Ui[2]*Htot*u['U-238'])/u['U']
print " Htot(U-238+U-234)=%.6e W/kg\tHtot(U)=%.6e W/kg" % (Htot,HtotU)
print " Isotope fractions: X_U234=%.5f%% X_U235=%.5f%% X_U238=%.5f%%" % \
(X_Ui[0]*100,X_Ui[1]*100,X_Ui[2]*100)
print " Mean u_U: %.5f amu" % u['U']
else:
# weighted mean with U-234, which releases a bit less heat
HtotU=(X_Ui[0]*H234*u['U-234']+X_Ui[2]*Htot*u['U-238'])/u['U']
In [ ]:
Etot=0.; Qgstot=0.
silent=0
nuc=raw_input("Known nuclides: Al-26, K-40, Fe-60, Th-232, U-235, U-238\n" \
"Enter nuclide for calculation: ")
if nuc == "Al-26":
Xiso=raw_input("Fraction of isotope in element (0<=Xiso<=1; hit ENTER for default):")
if Xiso != "": Xiso=float(Xiso)
calc_Al26(Xiso)
elif nuc == "K-40":
Xiso=raw_input("Fraction of isotope in element (0<=Xiso<=1; hit ENTER for default):")
if Xiso != "": Xiso=float(Xiso)
calc_K40(Xiso)
elif nuc == "Fe-60":
Xiso=raw_input("Fraction of isotope in element (0<=Xiso<=1; hit ENTER for default):")
if Xiso != "": Xiso=float(Xiso)
calc_Fe60(Xiso)
elif nuc == "Th-232":
calc_Th232()
elif nuc == "U-235":
buf=raw_input("Ratio r=U-238/U-235; hit ENTER for default (r=%.3f):" % rU85)
if buf != "": rU85=float(buf)
silent=1
calc_U238()
Etot=0.; Qgstot=0.
silent=0
calc_U235()
elif nuc == "U-238":
buf=raw_input("Ratio r=U-238/U-235; hit ENTER for default (r=%.3f):" % rU85)
if buf != "": rU85=float(buf)
silent=1
calc_U235()
Etot=0.; Qgstot=0.
silent=0
calc_U238()
else:
print "Unknown nuclide."
Abusaleem, K. (2014): Nuclear data sheets for $A=228$. Nucl. Data Sheets 116, 163-262, doi:10.1016/j.nds.2014.01.002.
Audi, G.; Kondev, F. G.; Wang, M.; Huang, W. J.; Naimi, S. (2017): The NUBASE2016 evaluation of nuclear properties. Chin. Phys. C 41(3), 030001, doi: 10.1088/1674-1137/41/3/030001; https://www-nds.iaea.org/amdc/
Basunia, M.S. (2014): Nuclear data sheets for $A=210$. Nucl. Data Sheets 121, 561-694, doi:10.1016/j.nds.2014.09.004.
Basunia, M.S., Hurst, A.M. (2016): Nuclear data sheets for $A=26$. Nucl. Data Sheets 132, 1-148, doi:10.1016/j.nds.2016.04.001.
Boehnke, P.; McKeegan, K. D.; Stephan, T.; Steele, R. C. J.; Trappitsch, R.; Davis, A. M. (2017): The rise and fall of iron-60. Meteorit. Planet. Sci. 52(S1), 6243, doi:10.1111/maps.12934
Browne, E., Tuli, J.K. (2013a): Nuclear data sheets for $A=231$. Nucl. Data Sheets 114(6-7), 751-840, doi:10.1016/j.nds.2013.05.002.
Browne, E., Tuli, J.K. (2013b): Nuclear data sheets for $A=60$. Nucl. Data Sheets 114(12), 1849-2022, doi:10.1016/j.nds.2013.11.002.
Chen, J. (2017): Nuclear data sheets for $A=40$. Nucl. Data Sheets 140, 1-376, doi:10.1016/j.nds.2017.02.001
Chisté, V., Bé, M.M. (2007): $^{218}$Po - Comments on evaluation of decay data. In Table de Radionucléides, CEA/LNHB, Gif-sur-Yvette, France.
Chisté, V., Bé, M.M. (2010a): $^{214}$Pb - Comments on evaluation of decay data. In Table de Radionucléides, CEA/LNHB, Gif-sur-Yvette, France.
Chisté, V., Bé, M.M. (2010b): ${}_{82}^{214}\mathrm{Pb}_{132}$. In Table de Radionucléides, CEA/LNHB, Gif-sur-Yvette, France.
Chisté, V., Bé, M.M. (2010c): ${}_{84}^{218}\mathrm{Po}_{134}$. In Table de Radionucléides, CEA/LNHB, Gif-sur-Yvette, France.
Chisté, V., Bé, M.M., Kellett, M.A. (2014a): $^{210}$Bi - Comments on evaluation of decay data. In Table de Radionucléides, CEA/LNHB, Gif-sur-Yvette, France.
Chisté, V., Bé, M.M., Kellett, M.A. (2014b): ${}_{83}^{210}\mathrm{Bi}_{127}$. In Table de Radionucléides, CEA/LNHB, Gif-sur-Yvette, France.
Goldmann, A., Brennecka, G., Noordman, J., Weyer, S., Wadhwa, M. (2015): The uranium isotopic composition of the Earth and the Solar System. Geochim. Cosmochim. Acta 148, 145-158, doi:10.1016/j.gca.2014.09.008.
Huang, X., Wang, B. (2009): $^{223}$Fr - Comments on evaluation of the decay data. In Table de Radionucléides, CEA/LNHB, Gif-sur-Yvette, France.
Huang, X., Wang, B. (2011a): $^{234}$Pa - Comments on evaluation of the decay data. In Table de Radionucléides, CEA/LNHB, Gif-sur-Yvette, France.
Huang, X., Wang, B. (2011b): $^{234}$Pa$^\mathrm{m}$ - Comments on evaluation of the decay data. In Table de Radionucléides, CEA/LNHB, Gif-sur-Yvette, France.
Huang, X., Wang, B. (2011c): ${}_{91}^{234}\mathrm{Pa}_{143}$. In Table de Radionucléides, CEA/LNHB, Gif-sur-Yvette, France.
Huang, X., Wang, B. (2011d): ${}_{91}^{234}\mathrm{Pa}_{143}^\mathrm{m}$. In Table de Radionucléides, CEA/LNHB, Gif-sur-Yvette, France.
Huang, X., Wang, B. (2012): ${}_{87}^{223}\mathrm{Fr}_{136}$. In Table de Radionucléides, CEA/LNHB, Gif-sur-Yvette, France.
Huang, X., Wang, B. (2014): ${}_{92}^{235}\mathrm{U}_{143}$. In Table de Radionucléides, CEA/LNHB, Gif-sur-Yvette, France.
Kondev, F.G., Lalkovski, S. (2011): Nuclear data sheets for $A=207$. Nucl. Data Sheets 112(3), 707-853, doi:10.1016/j.nds.2011.02.002.
Kondev, F., McCutchan, E., Singh, B., Tuli, J. (2016): Nuclear data sheets for $A=227$. Nucl. Data Sheets 132, 257-354, doi:10.1016/j.nds.2016.01.002.
Luca, A. (2009a): $^{228}$Ra - Comments on evaluation of decay data. In Table de Radionucléides, CEA/LNHB, Gif-sur-Yvette, France.
Luca, A. (2009b): $^{234}$Th - Comments on evaluation of decay data. In Table de Radionucléides, CEA/LNHB, Gif-sur-Yvette, France.
Luca, A. (2010): ${}_{90}^{234}\mathrm{Th}_{144}$. In Table de Radionucléides, CEA/LNHB, Gif-sur-Yvette, France.
Luca, A. (2012): ${}_{88}^{228}\mathrm{Ra}_{140}$. In Table de Radionucléides, CEA/LNHB, Gif-sur-Yvette, France.
Meija, J., Coplen, T.B., Berglund, M., Brand, W.A., De Bièvre, P., Gröning, M., Holden, N.E., Irrgeher, J., Loss, R.D., Walczyk, T., Prohaska, T. (2016a): Atomic weights of the elements 2013 (IUPAC Technical Report). Pure Appl. Chem. 88(3), 265-291, doi:10.1515/pac-2015-0305.
Meija, J., Coplen, T.B., Berglund, M., Brand, W.A.,De Bièvre, P., Gröning, M., Holden, N.E., Irrgeher, J., Loss, R.D., Walczyk, T., Prohaska, T. (2016b): Isotopic compositions of the elements 2013 (IUPAC Technical Report), Pure Appl. Chem., 88(3), 293-306, doi:10.1515/pac-2015-0503.
Mohr, P.J., Newell, D.B., Taylor, B.N. (2016): CODATA recommended values of the fundamental physical constants: 2014. J. Phys. Chem. Ref. Data 45(4), 043102, doi:10.1063/1.4954402.
Mougeot, X. (2015): Reliability of usual assumptions in the calculation of $\beta$ and $\nu$ spectra. Phys. Rev. C 91(5), 055504, doi: 10.1103/PhysRevC.91.055504
Naumenko, M.O., Mezger, K., Nägler, T.F., Villa, I.M. (2013): High precision determination of the terrestrial $^{40}$K abundance. Geochim. Cosmochim. Acta 122, 353-362, doi:10.1016/j.gca.2013.08.019.
Nichols, A.L. (2010): $^{208}$Tl - Comments on evaluation of decay data. In Table de Radionucléides, CEA/LNHB, Gif-sur-Yvette, France.
Nichols, A.L. (2011a): $^{212}$Bi - Comments on evaluation of decay data. In Table de Radionucléides, CEA/LNHB, Gif-sur-Yvette, France.
Nichols, A.L. (2011b): $^{212}$Pb - Comments on evaluation of decay data. In Table de Radionucléides, CEA/LNHB, Gif-sur-Yvette, France.
Nichols, A.L. (2011c): ${}_{82}^{212}\mathrm{Pb}_{130}$. In Table de Radionucléides, LNHB/CEA, Gif-sur-Yvette, France.
Nichols, A.L. (2011d): ${}_{83}^{212}\mathrm{Bi}_{129}$. In Table de Radionucléides, LNHB/CEA, Gif-sur-Yvette, France.
Nichols, A.L. (2016): ${}_{81}^{208}\mathrm{Tl}_{127}$. In Table de Radionucléides, LNHB/CEA, Gif-sur-Yvette, France.
Singh, B., Abriola, D., Baglin, C., Demetriou, V., Johnson, T., McCutchan, E., Mukherjee, G., Singh, S., Sonzogni, A., Tuli, J. (2013): Nuclear data sheets for $A=211$. Nucl. Data Sheets 114(6-7), 661-749, doi:10.1016/j.nds.2013.05.001.
Trappitsch, R.; Boehnke, P.; Stephan, T.; Telus, M.; Savina, M. R.; Pardo, O.; Davis, A. M.; Dauphas, N.; Huss, G. R. (2017): The life and death of iron-60. Meteorit. Planet. Sci. 52(S1), 6299, doi:10.1111/maps.12934
Wang, M.; Audi, G.; Kondev, F. G.; Huang, W. J.; Naimi, S.; Xu, X. (2017): The AME2016 atomic mass evaluation. (II). Tables, graphs and references. Chin. Phys. C 41(3), 030003, doi: 10.1088/1674-1137/41/3/030003; https://www-nds.iaea.org/amdc/
Wu, S.-C. (2009): Nuclear data sheets for $A=214$. Nucl. Data Sheets 110(3), 681-748, doi:10.1016/j.nds.2009.02.002.