You can interact with this notebook online: Launch notebook

Radioactive Decay Energy

Within the ejecta of a supernova, the \(\gamma\)-rays largely come from the decay of \(^{56}Ni\) into \(^{56}Co\), which releases a significant amount of energy.

When \(^{56}Ni\) decays into \(^{56}Co\) it can release a \(\gamma\)-ray at several different transition levels. Each transition level has an energy and an associated probability out of 100 decays. For example, the transition from Energy level 9 to Energy level 7 has an energy of 0.270 Mev and a probability of 36.5 out of 100 decays. To find the total energy per decay you multipliy each energy with its associated probability and add them all up.

import numpy as np
import pandas as pd
from astropy import units as u
from astropy import constants as const

from tardis.energy_input.gamma_ray_channel import (
    create_isotope_dicts,
    create_inventories_dict,
    calculate_total_decays,
    create_isotope_decay_df
)
from tardis.energy_input.gamma_ray_transport import get_taus
from tardis.energy_input.energy_source import get_nuclear_lines_database
from tardis.io.atom_data import (download_atom_data, AtomData)
from tardis.energy_input.tests.test_gamma_ray_channel import test_total_energy_production

/home/ceceliapowers/software/tardis/tardis/__init__.py:23: UserWarning: Astropy is already imported externally. Astropy should be imported after TARDIS.
  warnings.warn(

# energies of each transition
t_energies = np.array([0.270, 0.750, 0.480, 1.56, 0.812, 0.158]) * u.MeV
# probabilities of each transition
t_prob = np.array([.365, .495, .366, .140, .860, 1.00])

energy_per_decay = sum(t_energies * t_prob)
energy_per_decay

$1.7202 \; \mathrm{MeV}$

From the above cell, we get the energy per transition of 1.72 MeV. Note that this comes from a simplified scheme of energies and a more accurate total energy per \(^{56}Ni\) decay is 1.75 MeV. [Nadyozhin94]

The \(^{56}Co\) produced from the decay of \(^{56}Ni\) is also radioactive and will decay into \(^{56}Fe\) and release more \(\gamma\)-rays, however this decay is more complicated than the decay of \(^{56}Ni\). Whereas \(^{56}Ni\) only decays through electron capture, \(^{56}Co\) can decay either by electron capture, which occurs for 81 out of 100 cases, or through positron decay, which occurs for 19 out of 100 cases.

Positron decay produces positrons with a given kinetic energy, that will eventually annihilate with electrons to produce two 0.511 MeV \(\gamma\)-rays. The scheme of decays for \(^{56}Co\) is slightly more complicated than the \(^{56}Ni\) scheme, but to find the total energy per decay, you follow the same process. The total energy per decay from \(\gamma\)-rays is 3.61 MeV and the total kinetic energy of positrons is 0.12 MeV

The total rate of energy production for a mass of \(^{56}Ni\) at a given time is given by the following equation:

\[\epsilon = \frac{M_\odot}{56m_{u}}\frac{1}{\tau_{\text{Co}}-\tau_{\text{Ni}}}[[Q_{\text{Ni}}(\frac{\tau_{\text{Co}}}{\tau_{\text{Ni}}}-1)-Q_{\text{Co}}]\exp(-t/\tau_{\text{Ni}})+Q_{\text{Co}}\exp(-t/\tau_{\text{Co}})]\frac{M_{Ni0}}{M_\odot}\]

\(M_\odot\) is a solar mass. \(56_{u}\) is 56 atomic mass units.

\(\tau_{Ni}\) is the lifetime of \(^{56}Ni\) which is 8.80 days and \(\tau_{Co}\) is the lifetime of \(^{56}Co\) which is 111.3 days.

\(Q_{\text{Ni}}\) is the energy per decay of \(^{56}Ni\) which is 1.75 MeV and \(Q_{\text{Co}}\) is the sum of the energy per decay of \(^{56}Co\) from \(\gamma\)-rays and the kinetic energy from positrons which is 3.73 MeV

If we plug these values into the equation we get the equation:

\[\epsilon = (6.45e43\exp(-t/8.8)+1.45e43\exp(-t/111.3))\frac{M_{Ni0}}{M_\odot} erg/s\]

#time in days
time = 10 * u.day
#mass of Ni56 in solar masses
m_Ni56 = 1 * const.M_sun

energy_production_rate = (6.45e43 * np.exp(-time/(8.8*u.day)) + 1.45e43 * np.exp(-time/(111.3*u.day))) * (m_Ni56/ const.M_sun) * u.erg /u.s

print(f"The total energy production rate for 1 solar mass of 56Ni after 10 days is: {energy_production_rate:.2e}")

The total energy production rate for 1 solar mass of 56Ni after 10 days is: 3.40e+43 erg / s

The total total energy production integrated over time is given by:

\[\mathcal{E} = \mathcal{E}_{Ni} + \mathcal{E}_{Co} = 1.885e50 \frac{M_{Ni0}}{M_\odot} ergs\]

where

\[\mathcal{E}_{Ni} = 6.22e49 \frac{M_{Ni0}}{M_\odot} ergs\]

\[\mathcal{E}_{Co} = 1.26e50 \frac{M_{Ni0}}{M_\odot} erg\]

total_energy_production = 1.885e50 * (m_Ni56/const.M_sun) * u.erg
print(f"The total energy production for 1 solar mass of 56Ni is: {total_energy_production}")

The total energy production for 1 solar mass of 56Ni is: 1.885e+50 erg

Next we show how we can find the same result from functions in TARDIS-HE

First we create a DataFrame with isotope abundances. This DataFrame typically comes from a model, but for our purposes we create an ejecta with 5 shells, made entirely of \(^{56}Ni\).

mass_fractions = {0:1,
                  1:1,
                  2:1,
                  3:1,
                  4:1}

raw_isotope_abundance = pd.DataFrame(mass_fractions, index=pd.MultiIndex.from_tuples([(28, 56)], names=["atomic_number", "mass_number"]))
raw_isotope_abundance

		0	1	2	3	4
atomic_number	mass_number
28	56	1	1	1	1	1

Then we create an array of shell masses. We want to find the total energy production of 1 solar mass of \(^{56}Ni\), so we give each of the five shells a mass of \(\frac{1}{5}M_\odot\)

shell_masses = np.ones(5)*0.2*const.M_sun.to(u.g)
shell_masses

$[3.9768197 \times 10^{32},~3.9768197 \times 10^{32},~3.9768197 \times 10^{32},~3.9768197 \times 10^{32},~3.9768197 \times 10^{32}] \; \mathrm{g}$

Now we can use TARDIS-HE functions to get a DataFrame with information on the decays of this model. The first step is to create a dictionary which has the mass of each isotope in each shell.

isotope_dict = create_isotope_dicts(raw_isotope_abundance, shell_masses)
isotope_dict

{0: {'Ni56': 3.976819741396102e+32},
 1: {'Ni56': 3.976819741396102e+32},
 2: {'Ni56': 3.976819741396102e+32},
 3: {'Ni56': 3.976819741396102e+32},
 4: {'Ni56': 3.976819741396102e+32}}

Next, we create a dictionary that holds the inventory of each isotope given in the `radioactivedecay <https://radioactivedecay.github.io/>`__ package.

inventories = create_inventories_dict(isotope_dict)
inventories

{0: Inventory activities (Bq): {'Ni-56': 5.653446204017355e+48}, decay dataset: icrp107_ame2020_nubase2020,
 1: Inventory activities (Bq): {'Ni-56': 5.653446204017355e+48}, decay dataset: icrp107_ame2020_nubase2020,
 2: Inventory activities (Bq): {'Ni-56': 5.653446204017355e+48}, decay dataset: icrp107_ame2020_nubase2020,
 3: Inventory activities (Bq): {'Ni-56': 5.653446204017355e+48}, decay dataset: icrp107_ame2020_nubase2020,
 4: Inventory activities (Bq): {'Ni-56': 5.653446204017355e+48}, decay dataset: icrp107_ame2020_nubase2020}

Now we calculate the number of decays for each isotope in each shell and create a DataFrame to hold this information.

cumulative_decay_df = calculate_total_decays(inventories, time_delta= np.inf)
print(cumulative_decay_df)

                      number_of_decays
shell_number isotope
0            Ni56         4.281026e+54
             Co56         4.281026e+54
1            Ni56         4.281026e+54
             Co56         4.281026e+54
2            Ni56         4.281026e+54
             Co56         4.281026e+54
3            Ni56         4.281026e+54
             Co56         4.281026e+54
4            Ni56         4.281026e+54
             Co56         4.281026e+54

Next we use the kurucz_cd23_chianti_H_He.h5 dataset to get the decay radiation data for each isotope.

download_atom_data('kurucz_cd23_chianti_H_He')
atom_data_file = 'kurucz_cd23_chianti_H_He.h5'
atom_data = AtomData.from_hdf(atom_data_file)
gamma_ray_lines = atom_data.decay_radiation_data
gamma_ray_lines

WARNING:tardis.io.atom_data.atom_web_download:Atomic Data kurucz_cd23_chianti_H_He already exists in /home/ceceliapowers/Downloads/tardis-data/kurucz_cd23_chianti_H_He.h5. Will not download - override with force_download=True.
INFO:tardis.io.atom_data.util:
        Atom Data kurucz_cd23_chianti_H_He.h5 not found in local path.
        Exists in TARDIS Data repo /home/ceceliapowers/Downloads/tardis-data/kurucz_cd23_chianti_H_He.h5
INFO:tardis.io.atom_data.base:Reading Atom Data with: UUID = 6f7b09e887a311e7a06b246e96350010 MD5  = 864f1753714343c41f99cb065710cace
INFO:tardis.io.atom_data.base:Non provided Atomic Data: synpp_refs, photoionization_data, yg_data, two_photon_data, linelist_atoms, linelist_molecules

	A	Element	Z	N	Parent E(level)	Uncertainty	JPi	Metastable	Decay Mode	Decay Mode Value	...	Radiation	Rad subtype	Rad Energy	Uncertainty.1	EP Energy	Uncertainty.2	Rad Intensity	Uncertainty.3	Dose	Uncertainty.4
Isotope
Nn1	1	Nn	0	1	0.0	NaN	1/2+	False	B-	100.0	...	bm	NaN	301.3700	NaN	782.347	10.0	100.0	NaN	0.301000	NaN
Nn1	1	Nn	0	1	0.0	NaN	1/2+	False	B-	100.0	...	bm av	NaN	301.3700	NaN	NaN	NaN	100.0	NaN	0.301000	NaN
H3	3	H	1	2	0.0	NaN	1/2+	False	B-	100.0	...	bm	NaN	5.6817	12.0	18.591	3.0	100.0	NaN	0.005682	NaN
H3	3	H	1	2	0.0	NaN	1/2+	False	B-	100.0	...	bm av	NaN	5.6820	NaN	NaN	NaN	100.0	NaN	0.005680	NaN
He6	6	He	2	4	0	NaN	0+	False	B-	100.0	...	bm	NaN	1567.6200	54.0	3507.800	11.0	100.0	NaN	1.567600	NaN
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
Lv292	292	Lv	116	176	0	NaN	0+	False	A	100.0	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
Lv293	293	Lv	116	177	0	NaN	NaN	False	A	100.0	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
Ts293	293	Ts	117	176	0	NaN	NaN	False	A	100.0	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
Og294	294	Og	118	176	0	NaN	0+	False	A	100.0	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
Ts294	294	Ts	117	177	0	NaN	NaN	False	A	100.0	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

242295 rows × 41 columns

Finally, we can create a DataFrame with that holds all the information about the decays of each isotope in each shell.

isotope_decay_df = create_isotope_decay_df(cumulative_decay_df, gamma_ray_lines)
isotope_decay_df

		number_of_decays	decay_mode	radiation	radiation_energy_keV	radiation_intensity	energy_per_channel_keV	decay_energy_keV	decay_energy_erg
shell_number	isotope
0	Ni56	4.281026e+54	EC	bp	65.000	0.00060	0.000390	1.669600e+51	2.674994e+42
	Ni56	4.281026e+54	EC	bp	408.000	0.00003	0.000122	5.239975e+50	8.395366e+41
	Ni56	4.281026e+54	EC	bp	478.000	0.00003	0.000143	6.138991e+50	9.835748e+41
	Ni56	4.281026e+54	EC	bp av	100.000	0.00070	0.000700	2.996718e+51	4.801271e+42
	Ni56	4.281026e+54	EC	e	6.070	54.50000	3.308150	1.416227e+55	2.269047e+46
...	...	...	...	...	...	...	...	...	...
4	Co56	4.281026e+54	EC	g	3369.860	0.01010	0.340356	1.457072e+54	2.334487e+45
	Co56	4.281026e+54	EC	g	3451.232	0.94900	32.752192	1.402130e+56	2.246459e+47
	Co56	4.281026e+54	EC	g	3548.050	0.19550	6.936438	2.969507e+55	4.757674e+46
	Co56	4.281026e+54	EC	g	3600.800	0.01670	0.601334	2.574325e+54	4.124523e+45
	Co56	4.281026e+54	EC	g	3611.530	0.00860	0.310592	1.329651e+54	2.130335e+45

485 rows × 8 columns

Now we can use the data from the isotope decay DataFrame to find the energy production rate and total decay energy for the model and compare our results to the simplified analytical solution above.

def percent_error(observed, expected):
    return np.abs(expected - observed)/expected * 100

First we need to find the lifetimes of each isotope using the get_taus function.

tau_dict = get_taus(raw_isotope_abundance)

tau_co56 = tau_dict[0]['Co-56'] * u.s
tau_ni56 = tau_dict[0]['Ni-56'] * u.s
tau_ni56_days = tau_ni56.to(u.day)
tau_co56_days = tau_co56.to(u.day)
print(f"The lifetime of 56Ni is: {tau_ni56_days}")
print(f"The lifetime of 56Co is: {tau_co56_days}")
print(f"The percent error for 56Ni lifetime is: {percent_error(tau_ni56_days.value, 8.8)} %")
print(f"The percent error for 56Co lifetime is: {percent_error(tau_co56_days.value, 111.3)} %")

The lifetime of 56Ni is: 8.764372373400452 d
The lifetime of 56Co is: 111.41933800785463 d
The percent error for 56Ni lifetime is: 0.40485939317669306 %
The percent error for 56Co lifetime is: 0.1072219297885327 %

Next we find the energy per decay or Q value for both \(^{56}Ni\) and \(^{56}Co\) using the data in our isotope decay DataFrame, following a similar method as above. Here the “probabilities” are found in the radiation_intensity column and are in the form of number of decays per 100 decays so we have to divide by 100. The energy for each transition is found in the radiation_energy_keV column.

# For Ni56
ni56 = isotope_decay_df.xs('Ni56', level='isotope')
ni56_energy_per_decay_MeV = (ni56['radiation_intensity']/100 * (ni56['radiation_energy_keV']*u.keV.to("MeV"))).sum()/5 * u.MeV

# For Co56
co56 = isotope_decay_df.xs('Co56', level='isotope')
co56_energy_per_decay_MeV = (co56['radiation_intensity']/100 * (co56['radiation_energy_keV']* u.keV.to("MeV"))).sum()/5 * u.MeV

print(f"The energy per decay for Ni56 is: {ni56_energy_per_decay_MeV}")
print(f"The energy per decay for Co56 is: {co56_energy_per_decay_MeV}")
print(f"The percent error for Ni56 energy per decay is: {percent_error(ni56_energy_per_decay_MeV.value, 1.75)} %")
print(f"The percent error for Co56 energy per decay is: {percent_error(co56_energy_per_decay_MeV.value, 3.73)} %")

The energy per decay for Ni56 is: 1.7249433621 MeV
The energy per decay for Co56 is: 3.8521702007573437 MeV
The percent error for Ni56 energy per decay is: 1.4318078799999943 %
The percent error for Co56 energy per decay is: 3.2753405028778486 %

Next we can plug the \(\tau\) values and Q values into the energy production rate equation to get the coefficients in the second simplified energy production rate equation.

outside_factor = (1/ (tau_co56 - tau_ni56)) * ((const.M_sun.to('g'))/(56 * 1.67e-24 * u.g))
ni56_coeff = outside_factor * (ni56_energy_per_decay_MeV.to("erg") * ((tau_co56/tau_ni56) - 1) - co56_energy_per_decay_MeV.to("erg"))
co56_coeff = outside_factor * co56_energy_per_decay_MeV.to("erg")

print(f"The coefficient for Ni56 is: {ni56_coeff}")
print(f"The coefficient for Co56 is: {co56_coeff}")
print(f"The percent error for Ni56 coefficient is: {percent_error(ni56_coeff.value, 6.45e43)} %")
print(f"The percent error for Co56 coefficient is: {percent_error(co56_coeff.value, 1.45e43)} %")

The coefficient for Ni56 is: 6.280298471107531e+43 erg / s
The coefficient for Co56 is: 1.4795292955420012e+43 erg / s
The percent error for Ni56 coefficient is: 2.631031455697185 %
The percent error for Co56 coefficient is: 2.0365031408276604 %

Then we can use these coefficients to find the energy production rate after 10 days.

energy_production_rate_from_HE = ni56_coeff * np.exp(-time/tau_ni56_days) + co56_coeff * np.exp(-time/tau_co56_days) *  (m_Ni56/ const.M_sun)

print(f"The energy production rate from TARDIS-HE is: {energy_production_rate_from_HE}")
print(f"The percent error for energy production rate from TARDIS-HE is: {percent_error(energy_production_rate_from_HE.value, energy_production_rate.value)} %")

The energy production rate from TARDIS-HE is: 3.3591101069090006e+43 erg / s
The percent error for energy production rate from TARDIS-HE is: 1.0789985177098882 %

We can also sum the decay energy column of the DataFrame to find the total energy produced.

total_energy_production_HE = sum(isotope_decay_df["decay_energy_erg"]) * u.erg

print(f"The total energy production from TARDIS-HE is: {total_energy_production_HE}")
print(f"The percent error for total energy production from TARDIS-HE is: {percent_error(total_energy_production_HE.value, total_energy_production.value)} %")

The total energy production from TARDIS-HE is: 1.9126597273910088e+50 erg
The percent error for total energy production from TARDIS-HE is: 1.4673595432896005 %