This notebook is available at https://github.com/tardis-sn/tardis/tree/master/docs/physics/montecarlo/initialization.ipynb

---

Energy Packet Initialization

While it is instructive to think about tracking the propagation history of individual photons when illustrating the basic idea behind Monte Carlo radiative transfer techniques, there are important numerical reasons for using a different discretization scheme. Instead of thinking in the photon picture, it brings significant advantages to follow the idea of [AbbottLucy85] and [Lucy99a] and consider parcels of radiant energy as the fundamental building blocks of the Monte Carlo calculation. These basic Monte Carlo quanta are commonly referred to as “energy packets” or simply “packets”, and are composed of many photons with the same frequency.

During a Monte Carlo calculation, \(N\) (a large number) packets, all with a certain energy \(\varepsilon\), are created at the inner boundary of the computational domain (which is discussed in Model of Supernova Domain) known as the photosphere. Currently, the photosphere is modeled as a spherical blackbody with a radius \(R_\mathrm{phot}\) and temperature \(T_\mathrm{phot}\). In TARDIS, all packets are assigned identical energies, and the total energy of the packets is 1 erg (and thus each packet has an energy of \(\frac{1}{N}\) ergs).

Note

The indivisible energy packet scheme does not require that all packets have the same energy. This is just a convenient and simple choice adopted in TARDIS.

Since the photosphere is modeled as a blackbody, its total luminosity \(L\) (recall that luminosity is energy emitted divided by the time in which it is emitted) is

\[L=\frac{N\varepsilon}{\Delta t}=4 \pi R_{\mathrm{phot}}^2 \sigma_{\mathrm{R}} T_{\mathrm{phot}}^4\]

where \(\sigma_\mathrm{R}\) is the Stefan-Boltzmann constant and \(\Delta t\) is the physical duration of the simulation. In order to make this relationship hold (remembering that \(N\varepsilon = 1\) erg), we use

\[\Delta t = \frac{1}{L}=\frac{1}{4 \pi R_{\mathrm{phot}}^2 \sigma_{\mathrm{R}} T_{\mathrm{phot}}^4}.\]

During packet initialization, each packet is assigned an initial propagation direction \(\mu\) which is the cosine of the angle \(\theta\) which the packet’s path makes with the radial direction. Using a pseudo-random number generator which generates numbers \(z\) uniformly distributed on the interval \([0,1]\), the propagation direction is determined according to

\[\mu = \sqrt{z}.\]

This sampling is demonstrated in the code below.

Finally, each packet is assigned an initial frequency (or more precisely, the initial frequency of its consitiuent photons). Note that since each packet has the same energy, each packet will represent a different number of real photons. The sampling on packet frequencies is more involved than that of the propagation direction, as it involves sampling the Planck distribution (see below). TARDIS uses the technique described in [CarterCashwell75] and summarized in [BjorkmanWood01] for this purpose.

During the simulation, the energy of the packet remains constant in the local co-moving frame (see Reference Frames). This naturally ensures energy conservation and constitutes the main advantage of this discretization scheme. However, while the energy of the packets is conserved in the co-moving frame, the frequency of the constituent photons (in the local co-moving frame) may vary over the course of the simulation. Thus, a packet may represent several different numbers of real photons throughout their lifetimes.

We now demonstrate the TARDIS packet initialization framework:

import numpy as np
from tardis.montecarlo.packet_source import BlackBodySimpleSource
from astropy import units as u
from tardis import constants as const
import matplotlib.pyplot as plt

#The random number generator that will be used
rng = np.random.default_rng()

/usr/share/miniconda3/envs/tardis/lib/python3.7/importlib/_bootstrap.py:219: QAWarning: pyne.data is not yet QA compliant.
  return f(*args, **kwds)

[py.warnings         ][WARNING]  /usr/share/miniconda3/envs/tardis/lib/python3.7/site-packages/traitlets/traitlets.py:3036: FutureWarning: --rc={'figure.dpi': 96} for dict-traits is deprecated in traitlets 5.0. You can pass --rc <key=value> ... multiple times to add items to a dict.
  FutureWarning,
 (warnings.py:110)

The following cell contains values that you can change to see how it affects the spectrum:

# Seed for the pseudo-random number generator
seed = 1

# Radius of the supernova's photosphere in cm
r_phot = 1e15 * u.cm

# Number of packets generated
n_packets = 40000

You can either set a temperatature of the photosphere, which will determine its luminosity; or you can set the luminosity of the photosphere, which will determine its temperature.

# Temperature in K
temperature = 10000 * u.K

luminosity = 4 * np.pi * (r_phot**2) * const.sigma_sb * (temperature**4)

# Makes sure the luminosity is given in erg/s
luminosity = luminosity.to('erg/s')

print('Luminosity:', luminosity)

Luminosity: 7.125600863965567e+42 erg / s

# Luminosity in erg/s
luminosity = 7e42 * u.erg / u.s

temperature = (luminosity / (4 * np.pi * const.sigma_sb))**0.25 / np.sqrt(r_phot)

# Makes sure the termperature is given in K
temperature = temperature.to('K')

print('Temperature:', temperature)

Temperature: 9955.638923899083 K

We now generate the ensemble of packets. The array of packet energies and radii are also shown.

# We define our packet source
packet_source = BlackBodySimpleSource(seed)

radii, nus, mus, energies = packet_source.create_packets(
    temperature.value,
    n_packets,
    rng,
    r_phot)

# Sets the energies in units of ergs
energies *= u.erg

# Sets the frequencies in units of Hz
nus *= u.Hz

print('Energies:', energies)
print('Radii:', radii)

Energies: [2.5e-05 2.5e-05 2.5e-05 ... 2.5e-05 2.5e-05 2.5e-05] erg
Radii: [1.e+15 1.e+15 1.e+15 ... 1.e+15 1.e+15 1.e+15] cm

We set the timespan of the simulation so that each packet contributes the appropriate luminosity to the spectrum.

# Time of simulation
t_simulation = 1 * u.erg / luminosity
print('Time of simulation:', t_simulation)

# Array of luminosity contribution by each packet
lumin_per_packet = energies / t_simulation
print('Luminosity per packet:', lumin_per_packet)

Time of simulation: 1.4285714285714285e-43 s
Luminosity per packet: [1.75e+38 1.75e+38 1.75e+38 ... 1.75e+38 1.75e+38 1.75e+38] erg / s

We define important constants, and for comparison’s sake, we code the Planck distribution function

\[L_\nu (\nu)=\frac{8\pi R_\mathrm{phot}^2 h\nu^3}{c^2}\frac{1}{\exp\left(\frac{h\nu}{k_BT_\mathrm{phot}}\right)-1}\]

where \(L_\nu\) is the luminosity density with respect to frequency, \(\nu\) is frequency, \(h\) is Planck’s constant, \(c\) is the speed of light, and \(k_B\) is Boltzmann’s constant:

h = const.h.cgs
c2 = const.c.cgs**2
kB = const.k_B.cgs

def planck_function(nu):
    return 8 * np.pi**2 * r_phot**2 * h * nu**3 / (c2 * (np.exp(h * nu / (kB * temperature)) - 1))

We plot the Planck distribution and a histogram of the generated packet distribution:

# We set important quantites for making our histogram
bins = 200
nus_planck = np.linspace(min(nus), max(nus), bins)
bin_width = nus_planck[1] - nus_planck[0]

# In the histogram plot below, the weights argument is used
# to make sure our plotted spectrum has the correct y-axis scale
plt.hist(nus.value,
         bins=bins,
         weights=lumin_per_packet/bin_width)

# We plot the planck function for comparison
plt.plot(nus_planck, planck_function(nus_planck))

plt.xlabel('Frequency (Hz)')
plt.ylabel('Luminosity density w.r.t. frequency (erg/s/Hz)')
plt.show()

We finally plot the generated \(\mu\) density distribution, followed by the generated \(\theta=\arccos (\mu)\) density distribution, compared with the respective curves \(\rho = 2\mu\) and \(\rho = \sin(2\theta)\):

x = np.linspace(0, 1, 1000)

plt.hist(mus, bins=bins, density=True)
plt.plot(x, 2*x)
plt.xlabel('Propagation direction')
plt.ylabel('Probability density')
plt.show()

thetas = np.linspace(0, np.pi/2, 1000)

plt.hist(np.arccos(mus), bins=bins, density=True)
plt.plot(thetas, np.sin(2*thetas))
plt.xlabel('Angle with normal (rad)')
plt.ylabel('Probability density')
plt.show()

Custom Packet Source

TARDIS allows for the user to input a custom function that generates energy packets instead of the basic blackbody source described here. See Running TARDIS with a Custom Packet Source for more information.