You can interact with this notebook online: Launch notebook

Basic Spectrum Generation

The first and most simple way in which TARDIS calculates spectra is by calculating it directly from the Monte Carlo packets after the final Monte Carlo Iteration. This simply requires knowledge of each packet’s energy and frequency in the lab frame (see Reference Frames) at the end of the iteration. The only other quantity needed is the time duration of the simulation \(\Delta t\), which is calculated based off of the luminosity of the supernova’s photosphere (see Energy Packet Initialization).

Note

The only packets which are used for this calculation are the packets which escape the outer boundary of the computational domain – those reabsorbed into the photosphere are not included (see Packet Propagation).

The spectrum calculation is very straightforward. A packet of energy \(E_\mathrm{packet}\) contributes a luminosity

\[L_\mathrm{packet} = \frac{E_\mathrm{packet}}{\Delta t}\]

to the spectrum at its frequency.

In the code below, we will see an issue with merely relying on luminosity to give us a spectrum, which will allow us to develop the concept of luminosity density, and then correctly plot the TARDIS spectrum.

We start by importing the necessary packages, loading a configuration, and setting up a simulation object (see Setting up the Simulation):

from tardis.io.configuration.config_reader import Configuration
from tardis.simulation import Simulation
from tardis.montecarlo.spectrum import TARDISSpectrum
from tardis.io.atom_data.util import download_atom_data
from astropy import units as u
import numpy as np
import matplotlib.pyplot as plt

/home/runner/micromamba/envs/tardis/lib/python3.11/site-packages/setuptools_scm/git.py:162: UserWarning: "/home/runner/work/tardis/tardis" is shallow and may cause errors
  warnings.warn(f'"{wd.path}" is shallow and may cause errors')

# We download the atomic data needed to run the simulation
download_atom_data('kurucz_cd23_chianti_H_He')

Atomic Data kurucz_cd23_chianti_H_He already exists in /home/runner/Downloads/tardis-data/kurucz_cd23_chianti_H_He.h5. Will not download - override with force_download=True.

tardis_config = Configuration.from_yaml('tardis_example.yml')

sim = Simulation.from_config(tardis_config)

Number of density points larger than number of shells. Assuming inner point irrelevant
/home/runner/work/tardis/tardis/tardis/plasma/properties/radiative_properties.py:93: RuntimeWarning: invalid value encountered in divide
  (g_lower * n_upper) / (g_upper * n_lower)
/home/runner/work/tardis/tardis/tardis/plasma/properties/radiative_properties.py:93: RuntimeWarning: invalid value encountered in divide
  (g_lower * n_upper) / (g_upper * n_lower)

We now select a number of packets to run through the Monte Carlo simulation, and then run one Monte Carlo iteration (see Monte Carlo Iteration):

N_packets = 5000

# Using the commented out code below, we can also get the number of packets
# from the configuration -- try it out:
#N_packets = tardis_config.no_of_packets

sim.iterate(N_packets)

We call the arrays of each packet’s frequency and how much energy the packet has:

nus = sim.transport.transport_state.output_nu
nus

$[9.1172974 \times 10^{14},~1.6448802 \times 10^{14},~4.221768 \times 10^{14},~\dots,~1.1955602 \times 10^{15},~9.1111062 \times 10^{14},~4.662803 \times 10^{14}] \; \mathrm{Hz}$

energies = sim.transport.transport_state.output_energy
energies

$[0.00019893443,~-0.00019294683,~0.00020609497,~\dots,~0.00020733084,~0.00020628066,~0.000190025] \; \mathrm{erg}$

Notice that some energies are negative. This means that the packet ended up being reabsorbed into the photosphere (we will separate out these packets later). Also note that the number of elements of our arrays of frequencies and energies is the same as the number of packets that we ran (as it should be):

len(nus), len(energies)

(5000, 5000)

Note

The energies used for the spectrum are in the lab frame (see Reference Frames). Recall that while packets are all initialized with the same energy in the lab frame, throughout the simulation energy is conserved in the co-moving frame, allowing the lab frame energies of the packets to change though interactions (see Performing an Interaction), creating the varied lab-frame energies that we see at the end of the simulation.

TARDIS will then calculate a list of the packet luminosities by dividing each element in the array of energies by the time of the simulation \(\Delta t\):

luminosities = energies / sim.transport.transport_state.time_of_simulation
luminosities

$[2.1072654 \times 10^{39},~-2.0438401 \times 10^{39},~2.1831153 \times 10^{39},~\dots,~2.1962066 \times 10^{39},~2.1850823 \times 10^{39},~2.0128899 \times 10^{39}] \; \mathrm{\frac{erg}{s}}$

Now, as mentioned before, we only want to include the packets that make it through to the outer boundary. To do this, TARDIS creates an array of booleans (either True or False) called a mask that tells us if the packet should be counted in the spectrum. We then can use that mask to get an array of the frequencies and energies of only the packets which we are interested in:

emitted_mask = sim.transport.transport_state.emitted_packet_mask
emitted_mask

array([ True, False,  True, ...,  True,  True,  True])

emitted_nus = nus[emitted_mask]
emitted_nus

$[9.1172974 \times 10^{14},~4.221768 \times 10^{14},~5.6912554 \times 10^{14},~\dots,~1.1955602 \times 10^{15},~9.1111062 \times 10^{14},~4.662803 \times 10^{14}] \; \mathrm{Hz}$

emitted_luminosities = luminosities[emitted_mask]
emitted_luminosities

$[2.1072654 \times 10^{39},~2.1831153 \times 10^{39},~2.1648462 \times 10^{39},~\dots,~2.1962066 \times 10^{39},~2.1850823 \times 10^{39},~2.0128899 \times 10^{39}] \; \mathrm{\frac{erg}{s}}$

The length of these lists is the number of packets that made it out of the supernova:

len(emitted_nus), len(emitted_luminosities)

(3680, 3680)

Now, let’s plot frequency versus luminosity. We will see a very strange graph, which will lead us into developing a new strategy for plotting the spectrum:

plt.scatter(emitted_nus, emitted_luminosities)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Luminosity (erg/s)');

This is not the type of plot that we are looking for. We also cannot solve this problem by adding up the luminosities of all packets with the same frequency – in fact, since frequency is a continuum, it would be extremely unlikely that any two packets have the same exact frequency. To solve this problem, we will need the concept of luminosity density.

Luminosity Density

What we will have to do instead is plot a histogram where we bin up different frequencies that are close to each other into a certain number of bins. We then add up the luminosities of each packet in the bin and divide by the “width” of the bins. For example, if we are plotting between 0 Hz and 50 Hz in 5 bins, we would add up the luminosities of the packets between 0 Hz and 10 Hz, between 10 Hz and 20 Hz, between 20 Hz and 30 Hz, etc., and then divide each value by 10 Hz, which is the width of each bin (note that the bin widths need not be uniform). This will give us the luminosity density with respect to frequency, denoted by \(L_\nu\), and measured in ergs per second per Hertz. This can be interpreted as the luminosity per unit Hertz, i.e. how much luminosity will be in an interval with a width of 1 Hz.

The division step here is crucial, as otherwise the values on the y-axis will, for example, approximately double if we double the widths of our bins (as we would be adding about double the luminosity contributions into that bin). We clearly want the values on the y-axis to be independent of the number of bins we break the luminosity into, thus making luminosity density the best way to plot a spectrum.

Note that we can also have luminosity density with respect to wavelength, \(L_\lambda\), which we get by binning the packets by wavelength instead of frequency. Since the width of the bins would now have the dimensions of length, \(L_\lambda\) will have units of ergs per second per Angstrom.

Now, to generate our spectrum, we select the bounds for our spectrum and the number of bins that we group the packets into. Feel free to change the number of bins to see how it affects the spectrum!

# The lowest frequency we plot
freq_start = 1.5e14 * u.Hz

# The highest frequency we plot
freq_stop = 3e15 * u.Hz

# The number of bins
N = 500

The above information can also be retrieved from the configuration using the commented-out code below. Note that the configuration has the bounds specified in terms of wavelengths. Therefore we must convert these to frequencies, noting that the frequency is the speed of light divided by the wavelength (Astropy has a built-in way to do this, which we shall use). Additionally, since wavelength and frequency are inversely related, the lower bound for the wavelengths is the upper bound for the frequencies and vice versa.

#freq_start = tardis_config.spectrum.stop.to('Hz', u.spectral())

#freq_stop = tardis_config.spectrum.start.to('Hz', u.spectral())

#N = tardis_config.spectrum.num

Next, TARDIS generates the list of frequency bins. The array shown contain the boundaries between successive bins as well as the lower and upper bounds of the spectrum:

spectrum_frequency = np.linspace(freq_start, freq_stop, N+1)
spectrum_frequency

$[1.5 \times 10^{14},~1.557 \times 10^{14},~1.614 \times 10^{14},~\dots,~2.9886 \times 10^{15},~2.9943 \times 10^{15},~3 \times 10^{15}] \; \mathrm{Hz}$

Then, TARDIS creates a histogram where we add up the luminosity in each bin:

emitted_luminosity_hist = u.Quantity(np.histogram(emitted_nus,
                                                  weights=emitted_luminosities,
                                                  bins=spectrum_frequency,
                                                 )[0], "erg / s",)
emitted_luminosity_hist

$[8.6712772 \times 10^{39},~1.2946958 \times 10^{40},~1.5022143 \times 10^{40},~\dots,~0,~0,~0] \; \mathrm{\frac{erg}{s}}$

Finally, we input this information into the TARDISSpectrum class, which will generate the luminosity density with respect to both wavelength and frequency and allow us to plot both of these.

spectrum = TARDISSpectrum(spectrum_frequency, emitted_luminosity_hist)

spectrum.plot(mode='frequency')

spectrum.plot(mode='wavelength')

Note

Most of this process is done internally by TARDIS. Given a simulation object sim that has been run, calling sim.transport.spectrum will give you a TARDISSpectrum object that can then be plotted using the .plot() method. See, for example, our Quickstart Guide. This notebook just demonstrates how TARDIS generates this spectrum when sim.transport.spectrum is called.

You may notice that the bins are not uniformly spaced with respect to wavelength. This is because we use the same bins for both wavelength and frequency, and thus the bins which are uniformly spaced with respect to frequency are not uniformly spaced with respect to wavelength. This is okay though, since luminosity density allows us to alter the bins without significantly altering the graph!

Another thing you may notice in these graphs is the lack of a smooth curve. This is due to noise: the effects of the random nature of Monte Carlo simulations. This makes it very difficult to get a precise spectrum without drastically increasing the number of Monte Carlo packets. To solve this problem, TARDIS uses virtual packets and the formal integral method to generate a spectrum with less noise.