You can interact with this notebook online: Launch notebook

Energy Packet Initialization

Theory

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 the model section) known as the photosphere. Currently, the photosphere is modeled as a spherical blackbody with a radius \(r_\mathrm{boundary\_inner}\) and temperature \(T_\mathrm{inner}\). Both of these quantities are calculated as a part of the model, and \(T_\mathrm{inner}\) is additionally updated throughout the simulation as a part of the Updating Plasma and Convergence process.

In TARDIS, all packets are assigned identical energies in the lab frame (see Reference Frames), and the total (lab-frame) 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.

Note

Since all packets have the same total, and photon energy is proportional to frequency, higher-frequency packets will represent less real photons than lower-frequency packets.

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

\[L_\mathrm{inner}=\frac{N\varepsilon}{\Delta t}=4 \pi r_\mathrm{boundary\_inner}^2 \sigma_{\mathrm{R}} T_{\mathrm{inner}}^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_\mathrm{inner}}=\frac{1}{4 \pi r_\mathrm{boundary\_inner}^2 \sigma_{\mathrm{R}} T_{\mathrm{inner}}^4}.\]

Note

As will be shown in the code example, this will lead to unphysically small values for \(\Delta t\). It may be easier to think of the Monte Carlo packets not as packets of energy \(\epsilon\) going through a simulation of duration \(\Delta t\), but as packets of luminosity that carry an energy \(\epsilon\) over a time \(\Delta t\) (and thus truly being luminosity packets of luminosity \(\frac{\epsilon}{\Delta t}\)). Indeed, this view of the packets will be useful when deriving the Monte Carlo Estimators.

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 (see the image below). Using a pseudo-random number generator which generates numbers \(z\) uniformly distributed on the interval \([0,1]\), the propagation direction is determined (due to physical considerations beyond the scope of this documentation) according to

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

This sampling is shown in the “Code Example” section.

Finally, each packet is assigned an initial frequency (or more precisely, the initial frequency of its constituent 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 co-moving frequency of the packet (and thus the photons represented by the packet) may vary over the course of the simulation. Thus, a packet may represent several different numbers of real photons throughout their lifetimes.

Code Example

We now demonstrate the TARDIS packet initialization framework:

import numpy as np
from tardis.montecarlo.packet_source import BlackBodySimpleSource
from tardis.montecarlo.montecarlo_numba.packet_collections import (
    PacketCollection,
)
from astropy import units as u
from tardis import constants as const
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')

The following cell contains values that you can change to see how it affects the spectrum (in an actual simulation, the seed and number of packets are set in the Monte Carlo configuration, and the photospheric radius is calculated as a part of the model):

# Base seed for the packet source
base_seed = 1
# Seed offset resets the Pseudo Random Number Generator with a new seed at each iteration
seed_offset = 0

# Number of packets generated
n_packets = 40000

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

We set the temperature of the photosphere \(T_\mathrm{inner}\), which will determine the photospheric luminosity (in an actual simulation, \(T_\mathrm{inner}\) is initially calculated as a part of the model and updated as a part of the Updating Plasma and Convergence process):

# Temperature in K
temperature_inner = 10000 * u.K

luminosity_inner = (
    4
    * np.pi
    * (r_boundary_inner**2)
    * const.sigma_sb
    * (temperature_inner**4)
)

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

print("Luminosity of photosphere:", luminosity_inner)

Luminosity of photosphere: 7.125600863965567e+42 erg / s

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(base_seed=base_seed)

# Set radii and temperature from model
packet_source.radius = r_boundary_inner
packet_source.temperature = temperature_inner

# Create packets
packet_collection = packet_source.create_packets(n_packets)

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

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

print("Energies:", packet_collection.initial_energies)
print("Radii:", packet_collection.initial_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]

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_inner
print("Time of simulation:", t_simulation)

# Array of luminosity contribution by each packet
lumin_per_packet = packet_collection.initial_energies / t_simulation
print("Luminosity per packet:", lumin_per_packet)

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

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

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

where \(L_\nu\) is the luminosity density (see Basic Spectrum Generation) 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_boundary_inner**2
        * h
        * nu**3
        / (c2 * (np.exp(h * nu / (kB * temperature_inner)) - 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(packet_collection.initial_nus), max(packet_collection.initial_nus), bins).value
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(packet_collection.initial_nus.value, bins=bins, weights=lumin_per_packet / bin_width)

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

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(packet_collection.initial_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(packet_collection.initial_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 How to Run TARDIS with a Custom Packet Source for more information.