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

\[T_{\mathrm{rad\ estimated}} = \frac{h}{k_{\mathrm{B}}} \frac{\pi^4}{360 \zeta(5)} \frac{\bar \nu}{J}\]

and

\[W_\mathrm{estimated} = \frac{\pi J}{\sigma_{\mathrm{R}} T_{\mathrm{rad\ estimated}}^4}\]

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

\[T_\mathrm{rad\ estimated}= \frac{h}{k_{\mathrm{B}}} \frac{\pi^4}{360 \zeta(5)} \frac{\sum_i \varepsilon_i \nu_i l_i}{\sum_i \varepsilon_i l_i} = \frac{h}{k_{\mathrm{B}}} \frac{\pi^4}{360 \zeta(5)} \frac{\mathrm{real\ nu\_ bar\ estimator}}{\mathrm{real\ j\ estimator}}\]

and

\[W_\mathrm{estimated} = \frac{\sum_i \varepsilon_i l_i}{4\sigma_{\mathrm{R}} T_{\mathrm{rad\ estimated}}^4V\Delta t} = \frac{\mathrm{real\ j\ estimator}}{4\sigma_{\mathrm{R}} T_{\mathrm{rad\ estimated}}^4V\Delta t}.\]

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:

\[T_\mathrm{rad\ updated} = T_\mathrm{rad\ current} + d_{T_\mathrm{rad}}(T_\mathrm{rad\ estimated}-T_\mathrm{rad\ current})\]

and

\[W_\mathrm{updated} = W_\mathrm{current} + d_W(W_\mathrm{estimated}-W_\mathrm{current}).\]

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

\[T_\mathrm{inner\ estimated} = T_\mathrm{inner\ current} * \left(\frac{L_\mathrm{output}}{L_\mathrm{requested}}\right)^{\mathrm{t\_inner\_update\_exponent}}\]

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):

\[T_\mathrm{inner\ updated} = T_\mathrm{inner\ current} + d_{T_\mathrm{inner}}(T_\mathrm{inner\ estimated}-T_\mathrm{inner\ current}).\]

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):

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')

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

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')

# 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

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:

simulation_state.t_radiative

$[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] \; \mathrm{K}$

simulation_state.dilution_factor

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

simulation_state.t_inner

$9933.952 \; \mathrm{K}$

plasma.electron_densities

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:

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):

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

$[1.1326488 \times 10^{14},~9.8061971 \times 10^{13},~8.7712772 \times 10^{13},~7.9317802 \times 10^{13},~7.3345267 \times 10^{13},~6.6928878 \times 10^{13},~6.4282524 \times 10^{13},~6.0958957 \times 10^{13},~5.7624491 \times 10^{13},~5.4974415 \times 10^{13},~5.2774217 \times 10^{13},~5.1151601 \times 10^{13},~4.9427818 \times 10^{13},~4.7869045 \times 10^{13},~4.665492 \times 10^{13},~4.5519971 \times 10^{13},~4.4615124 \times 10^{13},~4.3829349 \times 10^{13},~4.314174 \times 10^{13},~4.2213738 \times 10^{13}] \; \mathrm{cm\,erg}$

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

$[9.2366972 \times 10^{28},~8.0044096 \times 10^{28},~7.1670115 \times 10^{28},~6.5223323 \times 10^{28},~6.0428657 \times 10^{28},~5.4875986 \times 10^{28},~5.2989335 \times 10^{28},~5.0217767 \times 10^{28},~4.7280809 \times 10^{28},~4.5065464 \times 10^{28},~4.3091838 \times 10^{28},~4.1273178 \times 10^{28},~3.9768577 \times 10^{28},~3.8148873 \times 10^{28},~3.6815176 \times 10^{28},~3.5592959 \times 10^{28},~3.4670082 \times 10^{28},~3.3847244 \times 10^{28},~3.3295518 \times 10^{28},~3.2424102 \times 10^{28}] \; \mathrm{Hz\,cm\,erg}$

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

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

$[9.4551879 \times 10^{10},~7.567263 \times 10^{10},~6.2755999 \times 10^{10},~5.2761131 \times 10^{10},~4.5475679 \times 10^{10},~3.8772285 \times 10^{10},~3.4871487 \times 10^{10},~3.1030738 \times 10^{10},~2.7579736 \times 10^{10},~2.4784064 \times 10^{10},~2.2450014 \times 10^{10},~2.0565932 \times 10^{10},~1.8811671 \times 10^{10},~1.727086 \times 10^{10},~1.5979504 \times 10^{10},~1.4819967 \times 10^{10},~1.3824559 \times 10^{10},~1.2941252 \times 10^{10},~1.2151917 \times 10^{10},~1.1355541 \times 10^{10}] \; \mathrm{\frac{erg}{s\,cm^{2}}}$

nu_bar = prefactor * nu_bar_estimator
nu_bar

$[7.7106609 \times 10^{25},~6.1768566 \times 10^{25},~5.1277933 \times 10^{25},~4.3385674 \times 10^{25},~3.7467096 \times 10^{25},~3.1789975 \times 10^{25},~2.8745246 \times 10^{25},~2.5563009 \times 10^{25},~2.2629132 \times 10^{25},~2.0316821 \times 10^{25},~1.8331155 \times 10^{25},~1.6594229 \times 10^{25},~1.5135472 \times 10^{25},~1.3763881 \times 10^{25},~1.2609351 \times 10^{25},~1.1588024 \times 10^{25},~1.0742962 \times 10^{25},~9.9938903 \times 10^{24},~9.3784896 \times 10^{24},~8.7221185 \times 10^{24}] \; \mathrm{\frac{Hz\,erg}{s\,cm^{2}}}$

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.

nu_bar/J

$[8.1549525 \times 10^{14},~8.1626032 \times 10^{14},~8.171001 \times 10^{14},~8.2230371 \times 10^{14},~8.2389307 \times 10^{14},~8.1991493 \times 10^{14},~8.2431946 \times 10^{14},~8.2379637 \times 10^{14},~8.2049851 \times 10^{14},~8.1975341 \times 10^{14},~8.1653201 \times 10^{14},~8.0687949 \times 10^{14},~8.0457885 \times 10^{14},~7.9694242 \times 10^{14},~7.8909526 \times 10^{14},~7.8191964 \times 10^{14},~7.7709258 \times 10^{14},~7.7225066 \times 10^{14},~7.717704 \times 10^{14},~7.6809361 \times 10^{14}] \; \mathrm{Hz}$

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:

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

$[10212.753,~10222.334,~10232.851,~10298.018,~10317.922,~10268.102,~10323.262,~10316.711,~10275.41,~10266.079,~10225.736,~10104.854,~10076.043,~9980.4086,~9882.1358,~9792.2728,~9731.8217,~9671.1846,~9665.17,~9619.1242] \; \mathrm{K}$

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

$[10069.627,~10066.985,~10064.832,~10090.027,~10092.612,~10060.357,~10080.614,~10070.037,~10042.107,~10030.183,~10002.775,~9935.1185,~9913.5184,~9858.5285,~9802.2402,~9750.1778,~9712.8422,~9675.4345,~9665.3587,~9635.2879] \; \mathrm{K}$

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

$[0.4815458,~0.38395229,~0.31710812,~0.25991915,~0.22230487,~0.19324105,~0.17011454,~0.15176296,~0.13706675,~0.12362115,~0.11375668,~0.10928659,~0.10111279,~0.096440481,~0.092832218,~0.089300099,~0.085391259,~0.081958935,~0.077151699,~0.073485967] \; \mathrm{}$

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

$[0.42473789,~0.34790225,~0.29790002,~0.25837561,~0.22926344,~0.20581558,~0.18642268,~0.1702087,~0.15647615,~0.14424298,~0.13401349,~0.12618181,~0.11790062,~0.11123709,~0.10538317,~0.099979266,~0.094836833,~0.090167556,~0.085378458,~0.08119239] \; \mathrm{}$

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

# 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

$7.9806817 \times 10^{42} \; \mathrm{\frac{erg}{s}}$

L_requested = sim.luminosity_requested
L_requested

$1.0592764 \times 10^{43} \; \mathrm{\frac{erg}{s}}$

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

$11444.76 \; \mathrm{K}$

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

$10991.518 \; \mathrm{K}$

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.

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)

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!

simulation_state.t_radiative

$[10069.626,~10066.983,~10064.831,~10090.025,~10092.61,~10060.355,~10080.612,~10070.035,~10042.105,~10030.182,~10002.773,~9935.1171,~9913.5171,~9858.5271,~9802.2388,~9750.1765,~9712.8409,~9675.4331,~9665.3574,~9635.2866] \; \mathrm{K}$

simulation_state.dilution_factor

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

simulation_state.t_inner

$10991.518 \; \mathrm{K}$

plasma.electron_densities

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