You can interact with this notebook online: Launch notebook

# Updating Plasma and Convergence¶

As light travels through a real plasma, it has effects on the properties of the plasma due to light-matter interactions as well as the presence of extra energy from the light. Additionally, as previously discussed, properties of the plasma affect how light travels through it. This is a typical equilibrium problem. We solve for the plasma properties by finding a steady-state solution; that is, the actual plasma will be in a state such that the plasma state will not change as light propagates through it, because the effects of the light on the plasma and the effects of the plasma on the light are in equilibrium.

One of TARDIS’s main goals is to determine this plasma state (as we need the actual plasma properties in order to get an accurate spectrum). This is done in an iterative process. After each Monte Carlo iteration (which sends light through the supernova ejecta), TARDIS uses the Monte Carlo estimators how the propagating light affects the plasma state, after which the plasma state is updated (as will be explained below and demonstrated in the code example). We do this many times, and attempt to have the plasma state converge to the steady-state we are looking for. In fact, all but the last Monte Carlo iteration is used for this purpose (after which TARDIS will have the necessary plasma state for its last iteration which calculates the spectrum).

Note

For all but the last iteration, TARDIS uses the number of Monte Carlo packets specified in the Monte Carlo configuration under the `no_of_packets`

argument. This is because a different number of packets may be necessary to calculate the spectrum as opposed to calculate the plasma state.

After each iteration the values for radiative temperature and dilution factor are updated by calling the `advance_state`

method on a `Simulation`

object as will be shown below in the code example. The goal of this is to eventually have the radiative temperature and dilution factor converge to a single value so that the steady-state plasma state can be determined. To ensure that the simulation converges, TARDIS employs a convergence strategy. Currently, only one convergence strategy is
available: damped convergence. This will be described in the following sections.

Note

Unless otherwise noted, all user-supplied quantities mentioned on this page are supplied in the convergence section of the Monte Carlo configuration, which will be referenced as the convergence configuration.

## \(T_\mathrm{rad}\) and \(W\)¶

The main way in which the plasma state is updated is by using the `j`

and `nu_bar`

estimators (see here) to update two of the key plasma inputs (see Plasma): the radiative temperature \(T_\mathrm{rad}\) and dilution factor \(W\). Recall that each of these quantities takes on a different value in each shell.

Using the estimators \(J\) and \(\bar \nu\), we estimate these quantities using

and

where \(h\) is Planck’s constant, \(k_{\mathrm{B}}\) is Boltzmann’s constant, \(\sigma_{\mathrm{R}}\) is the Stefan–Boltzmann constant, and \(\zeta\) 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 \(J_{\mathrm{blackbody}}=\frac{\sigma_{\mathrm{R}} T^4}{\pi}\).

Recall (see Implementation), however, that when TARDIS stores these estimators, it leaves out the prefactor of \(\frac{1}{4\pi V\Delta t}\) for computational ease. That is, \(J=\frac{1}{4\pi V\Delta t}\sum_i \varepsilon_i l_i=\frac{1}{4\pi V\Delta t}*\mathrm{real\ j\ estimator}\), and \(\bar \nu=\frac{1}{4\pi V\Delta t}\sum_i \varepsilon_i \nu_i l_i=\frac{1}{4\pi V\Delta t}*\mathrm{real\ nu\_ bar\ estimator}\). These factors are then included in our calculations; specifically, using the previous relations we have

and

While TARDIS can then update the plasma state using the estimated values, there is a good chance that these estimated values would “overshoot” the true value we want to converge to (for example, if the current value of the dilution factor in some cell the dilution factor is 0.4, and the true steady-state value TARDIS wants to find is 0.45, there is a good chance that the estimated value will be greater than 0.45). This could make the simulation take longer to converge or, at worst, make it so
the simulation does not converge at all. To account for this, users can set (in the convergence configuration) a “damping constant” for both the radiative temperature (\(d_{T_\mathrm{rad}}\)) and the dilution factor (\(d_W\)). When `advance_state`

is called, these quantities update as follows:

and

This means, for example, if the damping constant is 0.5, the updated value is halfway between the current value and the estimated value. If the damping constant is 0.7, the updated value is 70% of the way between the current value and the estimated value, and so on. **If the damping constant is 1, then the updated value is exactly the estimated value, and if the damping constant is zero, the value stays the same throughout the simulation and is not updated.**

The updated \(T_\mathrm{rad}\) 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.

## \(T_\mathrm{inner}\)¶

The temperature of the inner boundary, \(T_\mathrm{inner}\), plays a unique role in the simulation, as it is the primary determiner of the output luminosity. This is because the the luminosity of the inner boundary is proportional to \(T_\mathrm{inner}^4\) (see Energy Packet Initialization). Thus, \(T_\mathrm{inner}\) is updated throughout the simulation in order to match the output luminosity to the requested luminosity specified in the supernova configuration between the bounds specified in the supernova configuration. However, there is not necessarily a quartic relationship between \(T_\mathrm{inner}\) and the output luminosity, as changing \(T_\mathrm{inner}\) also changes the frequency distribution of the initialized packets (once again see Energy Packet Initialization). This then affects the light-matter interactions, affecting which packets make it to the outer boundary, which also affects the output luminosity. Because of this, there is not an exact way to estimate \(T_\mathrm{inner}\). To do this estimation, we use

where \(L_\mathrm{output}\) is the output luminosity calculated by adding up the luminosity of each packet (see Basic Spectrum Generation) between the bounds specified in the supernova configuration, \(L_\mathrm{requested}\) is the luminosity requested also in the supernova configuration (requested between those bounds previously mentioned), and `t_inner_update_exponent`

is provided by the user in the
convergence configuration. Note that what we are doing is “correcting” the previous value of the inner temperature by a factor of \(\left(\frac{L_\mathrm{output}}{L_\mathrm{requested}}\right)^{\mathrm{t\_inner\_update\_exponent}}\). Note that if \(\frac{L_\mathrm{output}}{L_\mathrm{requested}}\) is greater than 1, we want to lower \(T_\mathrm{inner}\) as the output luminosity is too high, and vice versa if the ratio is less than 1. Thus `t_inner_update_exponent`

should be negative.
Naively one might set `t_inner_update_exponent=-0.25`

, however as a default TARDIS uses `t_inner_update_exponent=-0.5`

as -0.25 may undershoot the correct \(T_\mathrm{inner}\) because of its previously mentioned effects on the initial frequency distribution.

After calculating the estimated \(T_\mathrm{inner}\), the quantity is updated using damped convergence with its own damping constant (once again set in the convergence configuration):

Once again, If the damping constant is 1, then the updated value is exactly the estimated value, and if the damping constant is zero, the value stays the same throughout the simulation and is not updated.

Additionally, because of the vast impact of \(T_\mathrm{inner}\) on the simulation, one may want to update it less frequently – i.e. allow \(W\) and \(T_\mathrm{rad}\) to reach a steady-state value for a particular \(T_\mathrm{inner}\) before updating \(T_\mathrm{inner}\). To do this, in the convergence configuration we set `lock_t_inner_cycles`

, which is the number of iterations to wait before updating \(T_\mathrm{inner}\). It is set to 1 by default, meaning
\(T_\mathrm{inner}\) would be updated every iteration.

## Convergence Information¶

During the simulation, information about the how \(T_\mathrm{rad}\), \(W\), and \(T_\mathrm{inner}\) are updated as well as a comparison of the total output luminosity and the requested luminosity are logged at the INFO level (see Configuring the Logging Output for TARDIS) as shown in the code below, to give users a better idea of how the convergence process is working.

In addition, TARDIS allows for the displaying of convergence plots, which allows users to visualize the convergence process for \(T_\mathrm{rad}\), \(W\), \(T_\mathrm{inner}\), and the total luminosity of the supernova being modeled. For more information, see Convergence Plots.

## Convergence Criteria¶

TARDIS also allows users to stop the simulation if the simulation reaches a certain level of convergence, which is checked upon the call of the `advance_state`

method. To enable this, users must set `stop_if_converged=True`

in the convergence configuration. Also in the configuration, the quantities `hold_iterations`

, `threshold`

, and `fraction`

are be specified to determine convergence as follows:

For the simulation to be considered to have converged, for `hold_iterations`

successive iterations, the estimated values of \(T_\mathrm{rad}\), \(W\), and \(T_\mathrm{inner}\) may differ from the previous value by a fraction of at most `threshold`

in at least `fraction`

fraction of the shells (for \(T_\mathrm{inner}\), since there is only one value, the `fraction`

part does not apply). For example, if `hold_iterations=3`

, `threshold=0.05`

for all three quantities, and
`fraction=0.8`

, the simulation will be considered to have converged if for 3 successive iterations the estimated values of \(T_\mathrm{rad}\) and \(W\) differ from the current respective values by at most 5% in at least 80% of the shells, *and* the estimated \(T_\mathrm{inner}\) differs by at most 5%. See the convergence configuration schema for default values of these quantities.

Note

To determine convergence, we compare the estimated value, **not** the updated value (which is related to the estimated value via the damping constant), with the previous value. If \(T_\mathrm{inner}\) is locked (see the previous section), the estimated value will still be calculated so convergence can be checked as usual.

Note

`hold_iterations`

and `fraction`

are universal quantities, i.e. they are each a single value that applies to \(T_\mathrm{rad}\) and \(W\), and for `hold_iterations`

also \(T_\mathrm{inner}\). `threshold`

, on the other hand, is supplied for each quantity separately, so for instance you could require \(T_\mathrm{rad}\) to differ by less than 1%, \(W\) to differ by less than 3%, and \(T_\mathrm{inner}\) to differ by less than 5% for convergence to be reached.

## Updating Other Quantities¶

Other quantities in the plasma state can depend on other estimators. Currently, this is only implemented for the `j_blue`

estimator: If `radiative_rates_type`

in the plasma configuration is set to `detailed`

, the `j_blues`

plasma property will will be replaced with the value of the \(J^b_{lu}\) estimator (the raw estimator times the factor of
\(\frac{ct_\mathrm{explosion}}{4\pi V \Delta t}\), once again see Implementation), which would then affect other plasma properties that depend on the `j_blues`

values (see Plasma). Otherwise, the `j_blues`

values in the plasma are calculated as they typically are in the plasma calculations, and the \(J^b_{lu}\) estimator is only used for the formal integral. Even in
the former case, while the estimator does contribute to the updating of the plasma when the `advance_state`

method is called, it does **not** contribute to the determination of if the simulation has converged, and it does **not** show up in the convergence logging or convergence plots.

## 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 necessary packages and set up a simulation (see Setting Up the Simulation):

```
[1]:
```

```
from tardis.io.configuration.config_reader import Configuration
from tardis.simulation import Simulation
from tardis.io.atom_data.util import download_atom_data
import numpy as np
from scipy.special import zeta
from astropy import units as u, constants as const
# 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.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')
```

```
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.
```

```
[2]:
```

```
tardis_config = Configuration.from_yaml('tardis_example.yml')
# We manually put in the damping constants and t_inner_update_exponent for
# illustrative purposes:
damping_t_radiative = 0.5
damping_dilution_factor = 0.3
damping_t_inner = 0.7
t_inner_update_exponent = -0.5
# We set the above values into the configuration:
tardis_config.montecarlo.convergence_strategy.t_rad.damping_constant = damping_t_radiative
tardis_config.montecarlo.convergence_strategy.w.damping_constant = damping_dilution_factor
tardis_config.montecarlo.convergence_strategy.t_inner.damping_constant = damping_t_inner
tardis_config.montecarlo.convergence_strategy.t_inner_update_exponent = t_inner_update_exponent
```

```
[3]:
```

```
sim = Simulation.from_config(tardis_config)
simulation_state = sim.simulation_state
plasma = sim.plasma
transport = sim.transport
```

```
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 show the initial radiative temperature and dilution factor in each shell, the initial inner boundary temperature, and the initial electron densities in each shell:

```
[4]:
```

```
simulation_state.t_radiative
```

```
[4]:
```

```
[5]:
```

```
simulation_state.dilution_factor
```

```
[5]:
```

```
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])
```

```
[6]:
```

```
simulation_state.t_inner
```

```
[6]:
```

```
[7]:
```

```
plasma.electron_densities
```

```
[7]:
```

```
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
```

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

```
[8]:
```

```
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.montecarlo.no_of_packets
sim.iterate(N_packets)
```

We now show the values of the \(J\) and \(\bar \nu\) estimators, attaching their proper units (note that these are the raw estimators, and the factors of \(\frac{1}{4\pi V \Delta t}\) etc are not included):

```
[9]:
```

```
j_estimator = transport.transport_state.j_estimator * (u.erg * u.cm)
j_estimator
```

```
[9]:
```

```
[10]:
```

```
nu_bar_estimator = transport.transport_state.nu_bar_estimator * (u.erg * u.cm * u.Hz)
nu_bar_estimator
```

```
[10]:
```

We show the values of \(J\) and \(\bar \nu\) including the prefactor:

```
[11]:
```

```
V = simulation_state.volume
Delta_t = transport.transport_state.time_of_simulation
prefactor = 1 / (4 * np.pi * V * Delta_t)
J = prefactor * j_estimator
J
```

```
[11]:
```

```
[12]:
```

```
nu_bar = prefactor * nu_bar_estimator
nu_bar
```

```
[12]:
```

As mentioned here, \(\bar \nu\) is not truly the mean frequency (as you can see by units). We show the true mean frequency, which will, as expected, be in units of Hz.

```
[13]:
```

```
nu_bar/J
```

```
[13]:
```

We show the calculations of the estimated and updated \(T_\mathrm{rad}\) and \(W\). Note that the `decompose()`

method is used to simplify the units:

```
[14]:
```

```
t_rad_estimated = ( (const.h.cgs / const.k_B.cgs)
* (np.pi**4 / (360 * zeta(5, 1)))
* (nu_bar_estimator / j_estimator) ).decompose()
t_rad_estimated
```

```
[14]:
```

```
[15]:
```

```
t_rad_updated = simulation_state.t_radiative + damping_t_radiative * (t_rad_estimated - simulation_state.t_radiative)
t_rad_updated
```

```
[15]:
```

```
[16]:
```

```
dilution_factor_estimated = ( j_estimator / (4 * const.sigma_sb.cgs * t_rad_estimated**4 * V * Delta_t) ).decompose()
dilution_factor_estimated
```

```
[16]:
```

```
[17]:
```

```
dilution_factor_updated = simulation_state.dilution_factor + damping_dilution_factor * (dilution_factor_estimated - simulation_state.dilution_factor)
dilution_factor_updated
```

```
[17]:
```

We show the output and requested luminosities, and use these to calculate the estimated and updated \(T_\mathrm{inner}\):

```
[18]:
```

```
# The output luminosity is calculated between the bounds specified below
nu_lower = 0
nu_upper = np.inf
# Using the commented out code below, we can also get the frequency bounds
# from the configuration -- try it out (note we convert between wavelength
# and frequency, and thus the order switches):
#nu_lower = tardis_config.supernova.luminosity_wavelength_end.to(u.Hz, u.spectral)
#nu_upper = tardis_config.supernova.luminosity_wavelength_start.to(u.Hz, u.spectral)
L_output = transport.transport_state.calculate_emitted_luminosity(0,np.inf)
L_output
```

```
[18]:
```

```
[19]:
```

```
L_requested = sim.luminosity_requested
L_requested
```

```
[19]:
```

```
[20]:
```

```
t_inner_estimated = simulation_state.t_inner * (L_output / L_requested)**t_inner_update_exponent
t_inner_estimated
```

```
[20]:
```

```
[21]:
```

```
t_inner_updated = simulation_state.t_inner + damping_t_inner * (t_inner_estimated - simulation_state.t_inner)
t_inner_updated
```

```
[21]:
```

We now advance the state of the simulation based on the estimators. This will also display a summary of the updated values of \(T_\mathrm{rad}\) and \(W\). Additionally, the `advance_state()`

method checks for convergence and returns the convergence status as shown below.

```
[22]:
```

```
sim.advance_state()
```

Shell No. | t_rad | next_t_rad | w | next_w |
---|---|---|---|---|

0 | 9.93e+03 K | 1.01e+04 K | 0.4 | 0.425 |

5 | 9.85e+03 K | 1.01e+04 K | 0.211 | 0.206 |

10 | 9.78e+03 K | 1e+04 K | 0.143 | 0.134 |

15 | 9.71e+03 K | 9.75e+03 K | 0.105 | 0.1 |

```
/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)
```

```
[22]:
```

```
False
```

Finally, we show the full updated \(T_\mathrm{rad}\), \(W\), and \(T_\mathrm{inner}\), as well as the updated electron densities which were updated along with the rest of the plasma based on the new values for \(T_\mathrm{rad}\), \(W\), and \(T_\mathrm{inner}\). Compare these with our calculations above and with the initial values at the beginning of the code example!

```
[23]:
```

```
simulation_state.t_radiative
```

```
[23]:
```

```
[24]:
```

```
simulation_state.dilution_factor
```

```
[24]:
```

```
array([0.42473789, 0.34790225, 0.29790003, 0.25837562, 0.22926344,
0.20581558, 0.18642268, 0.1702087 , 0.15647615, 0.14424298,
0.1340135 , 0.12618181, 0.11790063, 0.11123709, 0.10538317,
0.09997927, 0.09483683, 0.09016756, 0.08537846, 0.08119239])
```

```
[25]:
```

```
simulation_state.t_inner
```

```
[25]:
```

```
[26]:
```

```
plasma.electron_densities
```

```
[26]:
```

```
0 2.879194e+09
1 2.192349e+09
2 1.686013e+09
3 1.309369e+09
4 1.025252e+09
5 8.089293e+08
6 6.437049e+08
7 5.156985e+08
8 4.158632e+08
9 3.375895e+08
10 2.756466e+08
11 2.262224e+08
12 1.867888e+08
13 1.549413e+08
14 1.291439e+08
15 1.081434e+08
16 9.097366e+07
17 7.684827e+07
18 6.520044e+07
19 5.551064e+07
dtype: float64
```

```
[ ]:
```

```
```