You can interact with this notebook online: Launch interactive version

Volume-Based Estimators

When TARDIS runs, we enter into a loop with two main components: a Monte Carlo iteration occurs, and then the plasma state is updated based on the “estimators” described in this page. These estimators use the Monte Carlo packets to estimate how the light-matter interactions in the supernova affect the conditions in the ejecta. This concept was originally developed by [Lucy99a] and successively refined by [Lucy99b], [Lucy02] and [Lucy03].

Theory

J and nu_bar

Ordinarily, TARDIS is not concerned about the physical amount of time a packet spends traveling through the ejecta. Instead, we consider the “time of simulation” Δt which is chosen to be the amount of time in which the photosphere emits the ensemble of packets (see Energy Packet Initialization). When looking at the estimators, a slightly different interpretation of the packets is necessary. Here, we view the packets as not carrying a discrete amount of energy ε that is emitted in a time interval Δt, but as being a flow of energy that carries an energy ε over a time Δt – that is, each packet is carrying a luminosity (energy per unit time) of L=εΔt. Now, we can say that if a packet spends a time δt in the supernova’s ejecta, it contributes an energy of Lδt=εΔtδt into the radiation energy of the ejecta.

To account for the effects of the Monte Carlo packets on the ejecta, TARDIS uses the packets to first determine the average radiation energy density E throughout each shell, where the energy density is the total radiation energy in the shell divided by the volume of the shell V=43π(r3outerr3inner). Therefore, we add up the amount of energy each packet contributes to the radiation energy in that shell, and divide by the total volume of the shell:

E=1ViLiδti=1ViεiΔtδti=1VΔtiεiδti

where we sum over every Monte Carlo packet in the shell. Note that we are interested in the energy density in the co-moving frame (i.e. the energy density “according to the plasma,” see reference frames). Now, we note that the amount of time the Monte Carlo packet spends in a shell is δt=lic where l is the distance that the packet travels through the shell. Thus, our estimator is

E=1VΔtiεilic=1cVΔtiεili.

Using this energy density, we can then calculate the mean radiation intensity J in that shell using the relation J=c4πE, which gives us

J=14πVΔtiεili.

Since along any path the co-moving energy of the packet is continuously doppler shifted, we approximate this estimator using the co-moving energy at the beginning of the packet’s path (theoretically, the midpoint of the path would be a better option. However, we use the beginning of the path for computational ease at a very small cost to the estimator’s accuracy).

Next, we calculate the mean radiation frequency in each shell. For this, in each shell we add up the frequency of each packet weighted by the intensity they contribute to the shell. Remembering that intensity is c4π times the energy density, and as before each packet contributes an energy of εilicΔt and thus energy density of εilicVΔt to its shell, we have that each packet contributes an intensity of εili4πVΔt to its shell. So,

ˉν=iεili4πVΔtνi=14πVΔtiεiνili

where once again the co-moving energy and frequency of each packet are taken at the beginning of the packet’s path.

Note

Both estimators take on a different value in each shell.

These estimators allow us to calculate the radiative temperature Trad and dilution factor W in each shell using

Trad=hkBπ4360ζ(5)ˉνJ

and

W=πJσRT4rad

where h is Planck’s constant, kB is Boltzmann’s constant, σR is the Stefan–Boltzmann constant, and ζ is the Riemann zeta function. The equation for W comes from the fact that the dilution factor is the ratio of the actual mean intensity to that of a blackbody, which is Jblackbody=σRT4π.

The new Trad and W are then used as inputs to the updated plasma calculations which account for the effect of the Monte Carlo packets on the plasma state (precisely, these calculated Trad and W help determine the Trad and W used as inputs in the plasma calculations – see convergence for specifics).

J_blue

Another estimator, called the J_blue or Jblu estimator, is unlike the two previous estimators discussed. Instead of storing the mean intensity over the entire spectrum, it stores the intensity at a specific frequency. More precisely, since frequency is a continuum, it stores the intensity per unit frequency. In each shell, we record the intensity per unit frequency at the blue end (higher frequency end; this is where the “b” superscript in Jblu comes from) of each line transition – that is, if a line transition lu (from the lower energy level l to the upper energy level u, hence the lu in Jblu) has a frequency νlu, the mean intensity between νlu and νlu+dν is Jbludν. This means that the Jblu estimator has an entry for each atomic line in each shell. Now, using our previous J estimator, we have

Jbludν=14πVΔtiεidli

where dli is the infinitesimal distance that the packet travels while it has a co-moving frequency between νlu and νlu+dν.

Now, say the packet with lab frequency νlab has a co-moving frequency of νlu at a radius r1 and propagation direction μ1, and it has a co-moving frequency of νlu+dν at a radius r2 and propagation direction μ2. Then (see reference frames):

νlu=(1r1μ1ctexplosion)νlab

and

νlu+dν=(1r2μ2ctexplosion)νlab.

But then subtracting, we get

dν=(r2μ2r1μ1)νlabctexplosion=dlνlabctexplosion

(for the last equality, see propagation in a spherical domain).

But now inputting this into the equation for Jblu (using dlidν=ctexplosionνlab,i), we get

Jblu=ctexplosion4πVΔtiεiνlab,i.

Implementation

As previously discussed, a major component of each Monte Carlo iteration is the packet propagation process. During the packet propagation process this step, the J and ˉν estimators are updates every time a packet is moved to the next event location. Specifically, every time a packet is moved, εl is added to the “running total” J estimator in the shell where the packet is, and ενl is added to the “running total” ˉν estimator in the shell where the packet is (where l is the distance the packet is moved, and ε and ν are respectively the packet’s co-moving energy and frequency at the beginning of the packet’s path). The factor of 14πVΔt, for computational ease, is not attached to the estimators but is included during the calculation of Trad and W from the estimators. Specifically, we use

Trad=hkBπ4360ζ(5)iεiνiliiεili=hkBπ4360ζ(5)real nu_bar estimatorreal J estimator

and

W=iεili4σRT4radVΔt=real J estimator4σRT4radVΔt

After each Monte Carlo iteration, the advance_state() method is called on the Simulation object. The estimators are then used to update the radiative temperature and dilution factor according to the convergence strategy, and the plasma state is recalculated (see plasma calculations) using the updated radiative temperature and dilution factor as inputs. This process repeats until the final iteration or until convergence has been reached (see convergence).

Similarly, during the propagation process, every time a packet passes through a Sobolev point, meaning it comes in resonance with an atomic line (and thus reaches the frequency targeted by Jblu – not necessarily going through a line interaction), the Jblu for that atomic transition in the shell it is in is incremented by ενlab, where ε is the packet’s energy qnd νlab is the packet’s lab-frame frequency. As before, for computational ease, the factor ctexplosion4πVΔt is included in any calculations using the estimator.

Note

Since the J_blue estimator is updated every time a packet comes into resonance with an atomic line (not necessarily going through a line interaction), the estimator is only equal to zero in some shell for a specific line if no packet comes into resonance with that line within the shell.

If set to detailed mode (see plasma configuration), the J_blue plasma property will will be replaced with the value of the Jblu estimator (the raw estimator times the factor of ctexplosion4πVΔt). Otherwise, the J_blue in the plasma are calculated as they typically are in the plasma calculations, and the Jblu estimator is only used for the formal integral.

Code Example

We now show a detailed example of how the plasma is updated using the estimators after a Monte Carlo iteration. First, we import the needed packages and set up a simulation (see Setting Up the Simulation):

[1]:
from tardis.io.config_reader import Configuration
from tardis.simulation import Simulation
from tardis.io.atom_data.util import download_atom_data
import numpy as np

# We download the atomic data needed to run this notebook
download_atom_data('kurucz_cd23_chianti_H_He')
/usr/share/miniconda3/envs/tardis/lib/python3.8/site-packages/setuptools_scm/git.py:105: UserWarning: "/home/runner/work/tardis/tardis" is shallow and may cause errors
  warnings.warn(f'"{wd.path}" is shallow and may cause errors')
Iterations:
0/? [00:00<?, ?it/s]
Packets:
0/? [00:00<?, ?it/s]
[2]:
tardis_config = Configuration.from_yaml('tardis_example.yml')
sim = Simulation.from_config(tardis_config)

model = sim.model
plasma = sim.plasma
runner = sim.runner
/home/runner/work/tardis/tardis/tardis/plasma/properties/radiative_properties.py:93: RuntimeWarning: invalid value encountered in true_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 true_divide
  (g_lower * n_upper) / (g_upper * n_lower)

We show the initial radiative temperature, dilution factor, electron densities, and tau sobolevs in each shell:

[3]:
model.t_rad
[3]:
[9926.502, 9911.6354, 9896.8133, 9882.0354, 9867.3016, 9852.6117, 9837.9654, 9823.3627, 9808.8032, 9794.2868, 9779.8133, 9765.3825, 9750.9943, 9736.6484, 9722.3446, 9708.0828, 9693.8628, 9679.6844, 9665.5474, 9651.4516]K
[4]:
model.w
[4]:
array([0.40039164, 0.33245223, 0.28966798, 0.2577141 , 0.23224568,
       0.21120466, 0.19341188, 0.17811402, 0.16479446, 0.1530809 ,
       0.14269498, 0.13342262, 0.12509541, 0.11757849, 0.11076215,
       0.10455605, 0.09888494, 0.09368554, 0.08890421, 0.08449514])
[5]:
plasma.electron_densities
[5]:
0     2.865134e+09
1     2.182365e+09
2     1.678840e+09
3     1.303472e+09
4     1.020811e+09
5     8.059447e+08
6     6.411609e+08
7     5.137297e+08
8     4.144079e+08
9     3.364195e+08
10    2.747519e+08
11    2.256656e+08
12    1.863476e+08
13    1.546660e+08
14    1.289928e+08
15    1.080764e+08
16    9.094799e+07
17    7.685317e+07
18    6.520063e+07
19    5.552442e+07
dtype: float64
[6]:
plasma.tau_sobolevs
[6]:
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
atomic_number ion_number level_number_lower level_number_upper
14 5 0 36 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
35 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
34 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
1 35 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
0 33 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
20 0 176 191 6.110053e-12 2.765775e-12 1.290338e-12 6.191522e-13 3.049898e-13 1.539688e-13 7.953758e-14 4.198522e-14 2.261771e-14 1.241992e-14 6.944449e-15 3.949786e-15 2.283112e-15 1.340079e-15 7.980698e-16 4.818838e-16 2.948086e-16 1.826252e-16 1.144850e-16 7.258825e-17
177 191 4.581877e-12 2.074031e-12 9.676133e-13 4.642969e-13 2.287092e-13 1.154599e-13 5.964455e-14 3.148436e-14 1.696083e-14 9.313591e-15 5.207583e-15 2.961911e-15 1.712086e-15 1.004913e-15 5.984657e-16 3.613605e-16 2.210745e-16 1.369491e-16 8.585135e-17 5.443331e-17
178 191 3.048177e-13 1.379787e-13 6.437225e-14 3.088820e-14 1.521530e-14 7.681178e-15 3.967963e-15 2.094555e-15 1.128350e-15 6.196037e-16 3.464439e-16 1.970465e-16 1.138997e-16 6.685370e-17 3.981402e-17 2.404017e-17 1.470738e-17 9.110787e-18 5.711417e-18 3.621275e-18
138 144 2.205327e-10 9.983538e-11 4.658121e-11 2.235344e-11 1.101214e-11 5.559796e-12 2.872356e-12 1.516359e-12 8.169468e-13 4.486457e-13 2.508778e-13 1.427045e-13 8.249564e-14 4.842544e-14 2.884190e-14 1.741666e-14 1.065620e-14 6.601804e-15 4.138951e-15 2.624506e-15
177 190 8.129552e-14 3.679921e-14 1.716821e-14 8.237950e-15 4.057954e-15 2.048587e-15 1.058264e-15 5.586224e-16 3.009335e-16 1.652497e-16 9.239740e-17 5.255276e-17 3.037731e-17 1.783004e-17 1.061850e-17 6.411571e-18 3.922495e-18 2.429869e-18 1.523249e-18 9.658031e-19

29224 rows × 20 columns

We set the number of packets and we run one iteration of the Monte Carlo simulation:

[7]:
N_packets = 10000

# 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)
Iterations:
1/20 [00:00<00:04, 4.49it/s]
Packets:
6538/10000 [00:00<00:00, 65378.19it/s]

We now show the values of the three estimators previously mentioned (note that these are the raw estimators, and the factors of 14πVΔt etc are not included):

[8]:
runner.j_estimator
[8]:
array([1.13264880e+14, 9.80619708e+13, 8.77127720e+13, 7.93178023e+13,
       7.33452674e+13, 6.69288783e+13, 6.42825236e+13, 6.09589567e+13,
       5.76244911e+13, 5.49744149e+13, 5.27742170e+13, 5.11516006e+13,
       4.94278177e+13, 4.78690453e+13, 4.66549201e+13, 4.55199706e+13,
       4.46151244e+13, 4.38293494e+13, 4.31417403e+13, 4.22137375e+13])
[9]:
runner.nu_bar_estimator
[9]:
array([9.23669721e+28, 8.00440961e+28, 7.16701149e+28, 6.52233231e+28,
       6.04286574e+28, 5.48759864e+28, 5.29893354e+28, 5.02177670e+28,
       4.72808091e+28, 4.50654639e+28, 4.30918375e+28, 4.12731776e+28,
       3.97685768e+28, 3.81488727e+28, 3.68151763e+28, 3.55929590e+28,
       3.46700820e+28, 3.38472442e+28, 3.32955181e+28, 3.24241018e+28])
[10]:
# Because most rows of the j_blue estimatior are partially or mostly
# zero, we show just rows with all nonzero elements
runner.j_blue_estimator[np.all(runner.j_blue_estimator != 0, axis=1)]
[10]:
array([[2.54208693e-19, 3.98921452e-19, 3.51022121e-19, ...,
        4.78256257e-20, 1.83610133e-19, 4.62693098e-20],
       [5.25249079e-20, 3.47850697e-19, 3.53175242e-19, ...,
        4.68369084e-20, 4.76896351e-20, 1.42425885e-19],
       [5.25205919e-20, 3.47822114e-19, 3.03857848e-19, ...,
        4.68330598e-20, 4.76857164e-20, 1.42414182e-19],
       ...,
       [9.16402388e-19, 4.44096419e-19, 8.75723476e-19, ...,
        1.23926509e-18, 8.25398566e-19, 8.28841475e-19],
       [9.16135899e-19, 4.43967276e-19, 4.26276968e-19, ...,
        8.28490409e-19, 1.23557284e-18, 8.28600449e-19],
       [2.29780963e-18, 9.13630461e-19, 2.30523331e-18, ...,
        4.23437186e-19, 4.28206752e-19, 4.22703828e-19]])

We note that the shape of the j_blue estimator and the tau_sobolevs are the same: namely, each contain a value for each possible atomic line transition in each radial cell (as opposed to the other two estimators which just have one value for each cell):

[11]:
plasma.tau_sobolevs.shape, runner.j_blue_estimator.shape
[11]:
((29224, 20), (29224, 20))

We now advance the state of the simulation based on the estimators, and demonstrate this by showing the four quantities we showed before running the simulation. Compare them with their values above!

[12]:
# When advance_state is called, a brief summary of the updated t_rad's
# and w's is given
sim.advance_state();
Shell No. t_rad next_t_rad w next_w
0 9.93e+03 1.02e+04 0.4 0.482
5 9.85e+03 1.03e+04 0.211 0.193
10 9.78e+03 1.02e+04 0.143 0.114
15 9.71e+03 9.79e+03 0.105 0.0893
/home/runner/work/tardis/tardis/tardis/plasma/properties/radiative_properties.py:93: RuntimeWarning: invalid value encountered in true_divide
  (g_lower * n_upper) / (g_upper * n_lower)
[13]:
model.t_rad
[13]:
[10212.75, 10222.331, 10232.848, 10298.015, 10317.919, 10268.099, 10323.259, 10316.708, 10275.408, 10266.076, 10225.734, 10104.852, 10076.04, 9980.4059, 9882.1331, 9792.2701, 9731.819, 9671.1819, 9665.1674, 9619.1215]K
[14]:
model.w
[14]:
array([0.48154582, 0.3839523 , 0.31710814, 0.25991916, 0.22230488,
       0.19324106, 0.17011455, 0.15176297, 0.13706675, 0.12362115,
       0.11375669, 0.10928659, 0.10111279, 0.09644048, 0.09283222,
       0.0893001 , 0.08539126, 0.08195894, 0.0771517 , 0.07348597])
[15]:
plasma.electron_densities
[15]:
0     2.890015e+09
1     2.199909e+09
2     1.691446e+09
3     1.313845e+09
4     1.028790e+09
5     8.114482e+08
6     6.459942e+08
7     5.175646e+08
8     4.172918e+08
9     3.387728e+08
10    2.765659e+08
11    2.267949e+08
12    1.872454e+08
13    1.552250e+08
14    1.292984e+08
15    1.082114e+08
16    9.099949e+07
17    7.684337e+07
18    6.520026e+07
19    5.549692e+07
dtype: float64
[16]:
plasma.tau_sobolevs
[16]:
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
atomic_number ion_number level_number_lower level_number_upper
14 5 0 36 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
35 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
34 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
1 35 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
0 33 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
20 0 176 191 3.751288e-12 1.623835e-12 7.234954e-13 3.026676e-13 1.402851e-13 7.487558e-14 3.433806e-14 1.783878e-14 1.002376e-14 5.444603e-15 3.172899e-15 2.159921e-15 1.278110e-15 8.635643e-16 5.965582e-16 4.127740e-16 2.748029e-16 1.855430e-16 1.145663e-16 7.713120e-17
177 191 2.813059e-12 1.217700e-12 5.425431e-13 2.269679e-13 1.051986e-13 5.614856e-14 2.574982e-14 1.337715e-14 7.516734e-15 4.082862e-15 2.379330e-15 1.619706e-15 9.584440e-16 6.475795e-16 4.473539e-16 3.095357e-16 2.060723e-16 1.391371e-16 8.591229e-17 5.784002e-17
178 191 1.871440e-13 8.100975e-14 3.609368e-14 1.509946e-14 6.998534e-15 3.735387e-15 1.713055e-15 8.899396e-16 5.000647e-16 2.716200e-16 1.582893e-16 1.077539e-16 6.376225e-17 4.308142e-17 2.976104e-17 2.059243e-17 1.370934e-17 9.256347e-18 5.715472e-18 3.847913e-18
138 144 1.351661e-10 5.850663e-11 2.606589e-11 1.090035e-11 5.051694e-12 2.697050e-12 1.236483e-12 6.423822e-13 3.610450e-13 1.961192e-13 1.143170e-13 7.787523e-14 4.608972e-14 3.115875e-14 2.153771e-14 1.491086e-14 9.930657e-15 6.707647e-15 4.141899e-15 2.789347e-15
177 190 4.991160e-14 2.160543e-14 9.626242e-15 4.027047e-15 1.866519e-15 9.962334e-16 4.568742e-16 2.373482e-16 1.333680e-16 7.244146e-17 4.221601e-17 2.873817e-17 1.700551e-17 1.148990e-17 7.937332e-18 5.492048e-18 3.656313e-18 2.468691e-18 1.524331e-18 1.026248e-18

29224 rows × 20 columns