Coupling to a Radiative Transfer Solver¶
By default, CHIMES will evolve the thermo-chemistry with fixed UV fluxes (i.e. in the presence of a constant radiation field). However, CHIMES also has the option to evolve the radiation fluxes to account for the depletion of photons due to photoionisation, and the injection of additional photons from the recombination of HII and HeII. This option can be used to couple CHIMES to a Radiative Transfer (RT) code.
This option is controlled by the globalVariables.rt_update_flux parameter. This is an integer flag, if it is set to 0 we use a constant radiation field, if it is set to 1 then we evolve the photon fluxes as described below.
RT Coupling Implementation¶
When the globalVariables.rt_update_flux is set to 1, we increase the system of Ordinary Differential Equations (ODEs) that we are integrating to include ODEs for the evolution of the photon densities and fluxes. RT codes typically follow photons in multiple energy bins (we will discuss how we define the energy bins in CHIMES further below). For each energy bin, we add two ODEs, one to evolve the photon density and one to follow the relative reduction in photon fluxes. Following Equation 38 of Chan et al. (2021) (which we generalise to include photoionisation of all ions in the CHIMES network), the rate of change of photon density is given by:
where the \(\alpha\) ‘s are the case A and case B recombination coefficients of HII and HeII, as indicated by the subscripts, and \(\Theta_{i, \, j}\) is equal to unity if recombinations of species \(j\) lie in the energy band \(i\) or zero otherwise. The term \(S_{\rm{\gamma, i}}\) captures the external sources of photons. We include this term in the thermo-chemistry solver equations here so that the ionisation state will tend towards the correct equilibrium solution.
In practice, we actually integrate \(\ln(n_{\rm{\gamma}, \, i})\) in CVode, as this makes the integration more stable.
The inclusion of photons from the recombination of HII and HeII is optional, and is controlled by the globalVariables.rt_use_on_the_spot_approx parameter. This is an integer flag, if it is set to 0 then recombination photons are included as shown above, and the recombinations of HII and HeII in the network use the Case A rates. However, if this parameter is set to 1 we exclude the terms for photons from recombinations in the above equation (i.e. the \(\Theta_{i, \, j}\) are all set to zero), and we use the Case B recombination rates in the network, as we use the On-The-Spot approximation.
In addition to the photon densities, we also need to evolve the photon flux vector (Equation 39 in Chan et al. 2021). In principle, this would involve three ODEs (for a three-dimensional problem), one for each component of the flux vector. However, if we exclude the external source term in Equation 39 (which we can treat outside the thermo-chemistry solver in an operator-split fashion), then we see we can re-write the flux vector as follows:
where \(g_{\rm{\gamma}}(t)\) gives the relative reduction in each component of the flux vector, which is the same for each component. We therefore only need to add one ODE per energy bin for \(g_{\rm{\gamma}}(t)\), which we track in gasVariables.flux_reduction_factor. Note that \(g_{\rm{\gamma}}(t = 0) = 1\) by definition. After running CHIMES over the hydrodynamic time-step, we can then use the final value of \(g_{\rm{\gamma}}(t)\) to update each component of the photon flux vector \(\mathbf{f}_{\rm{\gamma}}(t)\).
Finally, we use CVode to integrate the extended system of ODEs to evolve the species abundances, temperature, photon densities and flux reduction factors. For the photon densities and fluxes we use the same relative tolerances as for the thermo-chemistry equations. However, as the typical absolute values of \(n_{\rm{\gamma}}\) and \(g_{\rm{\gamma}}\) can be dramatically different from the species abundances and thermal energies, we specify separate absolute tolerances for these variables via globalVars.rt_density_absoluteTolerance and globalVars.rt_flux_absoluteTolerance.
Be aware: we do not currently account for photons that are used up in dissociating molecules, we only include ionisation reactions in the depletion of photons. Molecules such as H2 absorb dissociating photons in discrete lines. Therefore, if the energy bins do not resolve the individual lines, accurately capturing how these molecules shield one another is not as straight-forward as for the ionising radiation. We intend to add options for treating the dissociating radiation in the future.
Reduced Speed of Light¶
We see in the above equations that, to calculate the photoionisation rates from the photon densities, we require the speed of light, \(c\). In CHIMES the thermo-chemistry equations always use the true speed of light. However, some RT codes use a reduced speed of light approximation, to improve the performance of the code. Nevertheless, it is important to still use the true speed of light in the thermo-chemistry equations, as this will affect the time-scale on which the photons are depleted. Therefore, if the RT code uses a reduced speed of light, the photon densities that are passed to CHIMES will need to be re-scaled to give the same photoionisation rates in the CHIMES solver (as these depend on the product of the photon density and the speed of light).
Defining Multiple Energy Bins in CHIMES¶
In its default mode, without coupling to RT, CHIMES supports photoionisation from multiple UV spectra. We first specify the shape of each spectrum, and then we compute the spectrum-averaged photoionisation cross sections for each species, following Equation 2.5 in Richings et al. 2014. These cross sections are stored in HDF5 data files (one for each spectrum) that are then passed to CHIMES, so they are only pre-computed once and don’t need to be computed by CHIMES on the fly.
We then specify the normalisation of the spectrum in the ionising (EUV; \(13.6 \, \rm{eV} - \infty\)) and non-ionising (FUV; \(6 - 13.6 \, \rm{eV}\)) bands separately. For the EUV band, this is specified as the number density of photons in the band (\(n_{\rm{\gamma}}\)). For the FUV band, the normalisation is specified as the energy density of photons in this band relative to that in the Habing (1968) interstellar radiation field (sometimes known as the \(G_{0}\) parameter).
When we couple CHIMES to the RT, we can use this same infrastructure for multiple UV spectra but apply it to the multiple photon energy bins from the RT. So previously one spectrum was defined over all frequencies, but now when we couple to the RT one spectrum will only correspond to one particular energy bin. However, we then need to change the definition of the photon density that is passed to CHIMES for each spectrum (i.e. energy bin). Previously, it was integrated over \(13.6 \, \rm{eV} - \infty\), but when coupling to RT it is now only integrated over the energy range of that particular bin. Similarly when we calculate the cross section \(\sigma_{i, \, j}\), previously we normalised the frequency-averaging by \(n_{\rm{\gamma}}\) over the \(13.6 \, \rm{eV} - \infty\) band, but with RT coupling it now uses that particular energy band.
The cross sections data files for each energy bin can be created using the Generate Cross Sections python script. This script has a parameter, rt_coupling_flag, that determines whether RT coupling is being used (as this will affect how the photon densities and average cross sections are defined). This flag is also recorded in the Header of the resulting HDF5 cross sections file that the script produces. It must match the rt_update_flux parameter that is being used in CHIMES, otherwise CHIMES will exit with an error message. The python script will also record whether this particular energy bin includes photons from recombinations of HII and/or HeII (via the integer flags rt_HII_recombination_flag and rt_HeII_recombination_flag). When CHIMES is initialised, if RT coupling is being used and the On-The-Spot Approximation is disabled (i.e. recombination radiation is being included), then it will check that HII and HeII recombinations are included in exactly one energy bin each.