Prepared by: Benoit Côté
This notebook describe the radioactive isotope implementation in NuPyCEE and shows how to run a SYGMA and OMEGA with radioactive yields.
In [1]:
# Import python modules
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
# Import the NuPyCEE codes
import sygma
import omega
The inputs that need to be provided to activate the radioactive option are:
The list of isotopes is declared in the yield_tables/decay_info.txt file and can be modified prior any simulation. The radioactive yields are found (or need to be added) in the yield_tables/ folder. Each stable yield table can have their associated radioactive yield table:
Each enrichment source can be activated independently by providing its input radioactive yield table. The radioactive yield table file format needs to be identical to its stable counterpart.
Warning: Radioactive isotopes will decay into stable isotopes. When using radioactive yields, please make sure that the stable yields do not include the decayed isotopes already.
If the radioactive isotopes you selected have only one decay channel, you can use the default decay option, which uses the following exponential law,
$N_r(t)=N_r(t_0)\,\mathrm{exp}\left[\frac{-(t-t_0)}{\tau}\right],$
$\tau=\frac{T_{1/2}}{\mathrm{ln}(2)},$
where $t_0$ is the reference time where the number of radioactive isotopes was equal to $N_0$. $T_{1/2}$ is the half-life of the isotope, which needs to be specified in yield_tables/decay_info.txt. The decayed product will be added to the corresponding stable isotope, as defined in yield_tables/decay_info.txt.
Below, a SYGMA simulation is ran with no star formation to better isolate the decay process. Here we choose Al-26 as an example, which decays into Mg-26.
In [2]:
# Number of timesteps in the simulaton.
# See https://github.com/NuGrid/NuPyCEE/blob/master/DOC/Capabilities/Timesteps_size_management.ipynb
special_timesteps = -1
nb_dt = 100
tend = 2.0e6
dt = tend / float(nb_dt)
# No star formation.
sfr = ''
starbursts = np.zeros(nb_dt)
# Dummy neutron star merger yields to activate the radioactive option.
nsmerger_table_radio = 'yield_tables/extra_table_radio_dummy.txt'
In [3]:
# Add 1 Msun of radioactive Al-26 in the gas.
# The indexes of this array reflect the order seen in the yield_tables/decay_file.txt file
# Index 0, 1, 2 --> Al-26, K-40, U-238
ism_ini_radio = [1.0, 0.0, 0.0]
In [4]:
# Run SYGMA (or in this case, the decay process)
s = sygma.sygma(iniZ=0.0, sfr=sfr, starbursts=starbursts, ism_ini_radio=ism_ini_radio,\
special_timesteps=special_timesteps, tend=tend, dt=dt,\
decay_file='yield_tables/decay_file.txt', nsmerger_table_radio=nsmerger_table_radio)
In [5]:
# Get the Al-26 (radioactive) and Mg-26 (stable) indexes in the gas arrays
i_Al_26 = s.radio_iso.index('Al-26')
i_Mg_26 = s.history.isotopes.index('Mg-26')
# Extract the evolution of these isotopes as a function of time
Al_26 = np.zeros(s.nb_timesteps+1)
Mg_26 = np.zeros(s.nb_timesteps+1)
for i_t in range(s.nb_timesteps+1):
Al_26[i_t] = s.ymgal_radio[i_t][i_Al_26]
Mg_26[i_t] = s.ymgal[i_t][i_Mg_26]
In [7]:
# Plot the evolution of Al-26 and Mg-26
%matplotlib nbagg
plt.figure(figsize=(8,4.5))
plt.plot( np.array(s.history.age)/1e6, Al_26, '--b', label='Al-26' )
plt.plot( np.array(s.history.age)/1e6, Mg_26, '-r', label='Mg-26' )
plt.plot([0,2.0], [0.5,0.5], ':k')
plt.plot([0.717,0.717], [0,1], ':k')
# Labels and fontsizes
plt.xlabel('Time [Myr]', fontsize=16)
plt.ylabel('Mass of isotope [M$_\odot$]', fontsize=16)
plt.legend(fontsize=14, loc='center left', bbox_to_anchor=(1, 0.5))
plt.subplots_adjust(top=0.96)
plt.subplots_adjust(bottom=0.15)
plt.subplots_adjust(right=0.75)
matplotlib.rcParams.update({'font.size': 14.0})
If the radioactive isotopes you selected have more than one decay channel, you need to use the provided decay module. This option can be activated by adding use_decay_module=True in the list of parameters when creating an instance of SYGMA and OMEGA. When using the decay module, the yield_tables/decay_file.txt file still needs to be provided as an input to define which radioactive isotopes are selected for the calculation.
Below we still run a SYGMA simulation with no star formation to better isolate the decay process. A fraction of K-40 decays into Ca-40, and another fraction decays into Ar-40.
In [8]:
# Add 1 Msun of radioactive K-40 in the gas.
# The indexes of this array reflect the order seen in the yield_tables/decay_file.txt file
# Index 0, 1, 2 --> Al-26, K-40, U-238
ism_ini_radio = [0.0, 1.0, 0.0]
In [9]:
# Number of timesteps in the simulaton.
# See https://github.com/NuGrid/NuPyCEE/blob/master/DOC/Capabilities/Timesteps_size_management.ipynb
special_timesteps = -1
nb_dt = 100
tend = 5.0e9
dt = tend / float(nb_dt)
In [10]:
# Run SYGMA (or in this case, the decay process)
# with the decay module
s = sygma.sygma(iniZ=0.0, sfr=sfr, starbursts=starbursts, ism_ini_radio=ism_ini_radio,\
special_timesteps=special_timesteps, tend=tend, dt=dt,\
decay_file='yield_tables/decay_file.txt', nsmerger_table_radio=nsmerger_table_radio,\
use_decay_module=True, radio_refinement=1)
In [11]:
# Get the K-40 (radioactive) and Ca-40 and Ar-40 (stable) indexes in the gas arrays
i_K_40 = s.radio_iso.index('K-40')
i_Ca_40 = s.history.isotopes.index('Ca-40')
i_Ar_40 = s.history.isotopes.index('Ar-40')
# Extract the evolution of these isotopes as a function of time
K_40 = np.zeros(s.nb_timesteps+1)
Ca_40 = np.zeros(s.nb_timesteps+1)
Ar_40 = np.zeros(s.nb_timesteps+1)
for i_t in range(s.nb_timesteps+1):
K_40[i_t] = s.ymgal_radio[i_t][i_K_40]
Ca_40[i_t] = s.ymgal[i_t][i_Ca_40]
Ar_40[i_t] = s.ymgal[i_t][i_Ar_40]
In [12]:
# Plot the evolution of Al-26 and Mg-26
%matplotlib nbagg
plt.figure(figsize=(8,4.5))
plt.plot( np.array(s.history.age)/1e6, K_40, '--b', label='K-40' )
plt.plot( np.array(s.history.age)/1e6, Ca_40, '-r', label='Ca-40' )
plt.plot( np.array(s.history.age)/1e6, Ar_40, '-g', label='Ar-40' )
# Labels and fontsizes
plt.xlabel('Time [Myr]', fontsize=16)
plt.ylabel('Mass of isotope [M$_\odot$]', fontsize=16)
plt.legend(fontsize=14, loc='center left', bbox_to_anchor=(1, 0.5))
plt.subplots_adjust(top=0.96)
plt.subplots_adjust(bottom=0.15)
plt.subplots_adjust(right=0.75)
matplotlib.rcParams.update({'font.size': 14.0})
In [13]:
# Add 1 Msun of radioactive U-238 in the gas.
# The indexes of this array reflect the order seen in the yield_tables/decay_file.txt file
# Index 0, 1, 2 --> Al-26, K-40, U-238
ism_ini_radio = [0.0, 0.0, 1.0]
In [14]:
# Number of timesteps in the simulaton.
# See https://github.com/NuGrid/NuPyCEE/blob/master/DOC/Capabilities/Timesteps_size_management.ipynb
special_timesteps = -1
nb_dt = 100
tend = 5.0e9
dt = tend / float(nb_dt)
In [15]:
# Run SYGMA (or in this case, the decay process)
# with the decay module
s = sygma.sygma(iniZ=0.0, sfr=sfr, starbursts=starbursts, ism_ini_radio=ism_ini_radio,\
special_timesteps=special_timesteps, tend=tend, dt=dt,\
decay_file='yield_tables/decay_file.txt', nsmerger_table_radio=nsmerger_table_radio,\
use_decay_module=True, radio_refinement=1)
In the case of U-238, there are many isotopes that are resulting from the multiple decay channels. Those new radioactive isotopes are added automatically in the list of isotopes in NuPyCEE.
In [16]:
print(s.radio_iso)
In [ ]: