%====================================================== \section{Numerical Solution Strategy}\label{chap:rt-numerical-strategy} %====================================================== %--------------------------------------------- \subsection{Discretization of Frequencies} %--------------------------------------------- The moments of the equation of radiative transfer, given in eqns.~\ref{eq:dEdt-freq}~-~\ref{eq:dFdt-freq}, are still frequency specific, and would need to be solved for each frequency between 0 Hz and infinity individually. This is obviously not a feasible approach. Instead, these equations need to be discretized in frequency first. I follow the approach of \cite{ramses-rt13} and for a rough approximation of multi-frequency, a relevant frequency range is split into a number $M$ of mutually exclusive groups of frequency ranges: \begin{align*} & [\nu_{00}, \nu_{01}\ : \ \nu_{10}, \nu_{11}\ : ... \ : \ \nu_{M0}, \nu_{M1}\ ] = [\nu_{0}, \infty [ \end{align*} The frequency ranges are chosen to be convenient for us. Specifically, since the interaction cross sections of ionizing species are zero below the ionizing frequency of their corresponding species, it makes sense to use the various ionizing frequencies given in eqs.~\ref{eq:nuIonHI}-\ref{eq:nuIonHeII} as the boundaries for the frequency groups. Rather than treating the photon energy densities and fluxes for individual frequencies, we now use their integrated averages over a frequency group. For any frequency interval, or group, $i$, we have \begin{align} E_i &= \int\limits_{\nu_{i0}}^{\nu_{i1}} E_\nu d\nu \\ \Fbf_i &= \int\limits_{\nu_{i0}}^{\nu_{i1}} \Fbf_\nu d\nu \end{align} giving us the following equations (discretized in frequency) to solve: \begin{align} \DELDT{E_i} + \nabla \cdot \Fbf_i &= - \sum\limits_{j}^{\absorbers} n_j \sigma_{i j}^N c E_i + \dot{E}_i^* + \dot{E}_i ^{rec} \label{eq:dEdt-group} \\ \DELDT{\Fbf_i} + c^2 \ \nabla \cdot \mathds{P}_i &= - \sum\limits_{j}^{\absorbers} n_j \sigma_{i j}^N c \Fbf_i \label{eq:dFdt-group} \end{align} The expression for the number weighted average cross section $\sigma_{ij}^N$ is given in eq. \ref{eq:sigma_N}. The radiation pressure tensor is discretized in the same manner: \begin{align} \mathds{P}_i &= \mathds{D}_i E_i \label{eq:pressure-tensor-group-start}\\ \mathds{D}_i &= \frac{1- \chi_i}{2} \mathds{I} + \frac{3 \chi_i - 1}{2} \mathbf{n}_i \otimes \mathbf{n}_i \label{eq:eddington-group} \\ \mathbf{n}_i &= \frac{\Fbf_i}{|\Fbf_i|} \\ \chi_i &= \frac{3 + 4 f_i ^2}{5 + 2 \sqrt{4 - 3 f_i^2}} \\ f_i &= \frac{|\Fbf_i|}{c E_i} \label{eq:pressure-tensor-group-end} \end{align} For the computation of the photo-heating and photo-ionization rates, which will be discussed subsequently, we need to introduce the mean photon energy $\overline{\epsilon}_i$ of the frequency bin $i$. The exact value cannot be known, as it would require the knowledge of the energy density at each frequency, which is not available as we are using quantities averaged over frequency intervals. Instead, we guess a reasonable spectrum. \cite{ramses-rt13} recommend taking an average value among all stellar radiation sources. GEAR-RT currently only works with a blackbody spectrum, which for a characteristic temperature $T_{bb}$ is given by: \begin{align} \mathcal{J}_\nu(T_{bb}) = \frac{2 \nu^2}{c^2} \frac{h \nu}{\exp\left(h\nu/k_B T_{bb}\right) - 1} \ . \label{eq:blackbody} \end{align} The temperature $T_{bb}$ in this case would be some effective temperature for stellar sources, not the local temperature of the gas. With an assumed spectral shape, the mean photon energy $\overline{\epsilon}_i$ in each frequency group $i$ can be estimated as \begin{align} \overline{\epsilon}_i \equiv \frac{E_i}{N_i} = \frac{ \int\limits_{\nu_{i0}}^{\nu_{i1}} E_\nu \ d\nu }{ \int\limits_{\nu_{i0}}^{\nu_{i1}} N_\nu \ d\nu } \approx \frac{ \int\limits_{\nu_{i0}}^{\nu_{i1}} \mathcal{J}_\nu \ d\nu }{ \int\limits_{\nu_{i0}}^{\nu_{i1}} \mathcal{J}_\nu / h\nu \ d\nu } \end{align} These integrals are performed using GSL integrator functions. Using the mean photon energy $\overline{\epsilon}_i$, the photo-heating rate $\mathcal{H}_{i,j}$ of the gas for the photon group $i$ and ionizing species $j$ then becomes: \begin{align} \mathcal{H}_{i, j} &= \left[ \int\limits_{\nu_{i0}}^{\nu_{i1}}\de \nu h \nu N_\nu \sigma_{j\nu} - h \nu_{ion,j}\ \int\limits_{\nu_{i0}}^{\nu_{i1}}\de \nu N_\nu \sigma_{j\nu} \ \right] \ c \ n_j \\ % &= \left[ \sigma_{ij}^E - h \nu_{ion,j}\ \sigma_{ij}^N /\ \overline{\epsilon}_i \right] E_i\ c \ n_j \label{eq:photoheating-group} \end{align} And the photo-ionization rate can be written as \begin{align} \Gamma_{i, j} &= c \ \int\limits_{\nu_{i0}}^{\nu_{i1}}\de \nu \ \sigma_{\nu j} N_\nu \\ &= c \ \sigma_{ij}^N N_i = c \ \sigma_{ij}^N E_i / \overline{\epsilon}_i \label{eq:photoionization-group} \end{align} Finally, the photon absorption (destruction) rates for a frequency bin $i$ and ionizing species $j$ are given by \begin{align} \DELDT{E_i} &= \deldt{(N_i \ \overline{\epsilon}_i)} = \overline{\epsilon}_i \DELDT{N_i} = -\overline{\epsilon}_i c \ \sigma_{i j} ^ N n_j \end{align} Here we have introduced the number- and energy-weighted average cross sections: \begin{align} \sigma_{ij}^N &= \frac{ \int\limits_{\nu_{i0}}^{\nu_{i1}}\de \nu \ N_\nu \ \sigma_{j\nu} } { \int\limits_{\nu_{i0}}^{\nu_{i1}}\de \nu \ N_\nu } \approx \frac{ \int\limits_{\nu_{i0}}^{\nu_{i1}}\de \nu \ \mathcal{J}(\nu) / (h \ \nu) \sigma_{j\nu} } { \int\limits_{\nu_{i0}}^{\nu_{i1}}\de \nu \ \mathcal{J}(\nu) / (h \ \nu) } \label{eq:sigma_N}\\ \sigma_{ij}^E &= \frac{ \int\limits_{\nu_{i0}}^{\nu_{i1}}\de \nu \ h \nu N_\nu \ \sigma_{j\nu} } { \int\limits_{\nu_{i0}}^{\nu_{i1}}\de \nu \ h \nu N_\nu } \approx \frac{ \int\limits_{\nu_{i0}}^{\nu_{i1}}\de \nu \ \mathcal{J}(\nu) \ \sigma_{j\nu} } { \int\limits_{\nu_{i0}}^{\nu_{i1}}\de \nu \ \mathcal{J}(\nu) } \label{eq:sigma_E} \end{align} %--------------------------------------------- \subsection{One RT Time Step} %--------------------------------------------- %--------------------------------------------- \subsubsection{Outline}\label{chap:rt-numerics-outline} %--------------------------------------------- For each photon frequency group, the equations \ref{eq:dEdt-group} and \ref{eq:dFdt-group} are solved with an operator-splitting strategy. Following the approach of \cite{ramses-rt13}, the equations are decomposed into three steps that are executed in sequence over the same time step $\Delta t$: \begin{enumerate} \item \textbf{Photon injection step}: the radiation from radiative sources is injected into the grid. \item \textbf{Photon transport step}: Photons are transported in space by solving the homogenized equations \ref{eq:dEdt-group} and \ref{eq:dFdt-group}, i.e. the right hand side of these equations is set to zero. \item \textbf{Thermochemistry step}: The rest of the source terms (the right hand side) of the equations \ref{eq:dEdt-group} and \ref{eq:dFdt-group} is solved. \end{enumerate} In what follows, the numerical and algorithmic aspects of these three steps are discussed individually. But first, let's list the fundamental approach upon which GEAR-RT is built: \begin{itemize} \item The underlying method to solve the hyperbolic conservation laws for hydrodynamics and for radiative transfer is a finite volume particle method. \item The particles are co-moving with the gas, and the motion of the particles is determined by the hydrodynamics. \item Radiation fields at particle positions are treated in a static frame of reference. (In contrast to the hydrodynamics, which are treated in a Lagrangian manner.) \item Particles are given individual time step sizes, for both the hydrodynamics and the radiative transfer individually. \item In a simulation step, the hydrodynamics are solved before the radiative transfer. Since RT is solved in a static frame of reference, this means that we can re-use all neighbour dependent data, in particular the smoothing lengths, partitions of unity (eq.~\ref{eq:psi}), and matrices $\mathcal{E}$ and $\mathcal{B}$ (eqs.~\ref{eq:matrix-E} and \ref{eq:matrix-B}). \end{itemize} %-------------------------------------------------------------------- \subsection{First step: Injection} \label{chap:injection-step} %-------------------------------------------------------------------- %------------------------------------------------ \subsubsection{Injecting the Energy Density} %------------------------------------------------ In the injection step, the radiation is gathered from radiating sources and injected into the gas. Radiation emitting sources are taken to be stars or entire stellar populations, which in SWIFT are represented by individual star particles. The equation to be solved in this step is \begin{equation*} \DELDT{E_i} = \dot{E}_i^* \label{eq:solve-dEdt} \end{equation*} where $\dot{E}_i^*$ is the total energy density in the frequency group $i$ emitted by all stars $k$ that are within compact support radius of the corresponding particle. If each star $k$ deposits some fraction $\xi_k$ of its respective total radiation energy density to be injected, then the total radiation energy density injected into a single gas particle $i$ is given by: \begin{align} \dot{E_i}^* = \sum_k E_{i,k}^* \xi_k \label{eq:injection_onto_particle} \ . \end{align} The exact choice how the fraction $\xi_k$ is determined will be discussed further below. Eq.~\ref{eq:injection_onto_particle} is solved using a simple finite difference discretization and first order forward Euler integration: \begin{equation} E_i(t + \Delta t) = E_i(t) + \Delta t \dot{E}_i^* \end{equation} for each particle. Star particles also have individual time step sizes, just like gas particles. Because both star and gas particles have individual time step sizes, the way the energy density $\dot{E}_i^*$ is deposited from stars onto particles is determined by the star particles' time steps. Each time a star particle is active, i.e. the simulation is at a step where the star particle finishes its time integration, the star particle ``pushes'' radiation onto both active and inactive gas particles. We now define the weights $\xi_s$ used in eq.~\ref{eq:injection_onto_particle} to distribute the total radiation energy ejected by a star $s$ over a time step $\Delta t$ onto the surrounding gas particles $p$. A natural way of depositing the energy density from the star on a neighboring gas particle is to make use of the already present partition of unity, $\psi_p(\x_p - \x_s)$ (see eqs. \ref{eq:psi} - \ref{eq:omega}), which is a fundamental building block for finite volume particle methods. % It guarantees to sum up to unity, and additionally is taking into account the particle % configuration due to its normalization. Then the total emitted energy of the star $s$, $E^*_i(\x_s)$, is fully and self-consistently distributed on the gas (here gas particles have the index $g$): \begin{align} \sum_g \psi_g(\x_s)\ E_i^*(\x_s) = \sum_g \psi(\x_s - \x_g, h_s) \ E_i^*(\x_s) = E_i^*(\x_s) \end{align} For the actual distribution of the radiation from stars to gas particles, the sum over neighboring stars from the gas particle's point of view looks like: \begin{align} E^*_i (\x = \x_{gas}) = \sum_{stars} E^*_i(\x_{star}) \ \psi(\x_{star} - \x_{gas}, h(\x_{star})) \end{align} suggesting \begin{align} \xi_s = \psi(\x_{s} - \x_{g}, h(\x_{s})) \ . \end{align} A potential problem here is that unless the particle distribution is sufficiently symmetric, the resulting injected energy density (and flux) will not be isotropic, which is what stellar radiation is typically assumed to be. So instead, we ``correct'' the weight $\psi_p$ in an attempt to improve the isotropy of the deposited net energy and flux. Suppose we divide space into 8 octants centered on a star. Each octant $i$ is assigned a weight $w_i$ defined as \begin{align*} w_i = \sum\limits_{\text{particles $j$ in $i$}} \xi (\x_j) \end{align*} We now apply the correction through a factor $\mu_i$ for each octant individually. Demanding that \begin{enumerate} \item The total weight must remain constant, i.e. $\sum_i \mu_i w_i = \sum_i w_i $ \item The modified weight of each quadrant should be equal, such that at least along these eight directions, the injected energy and flux are isotropic: $\mu_i w_i = \mu_j w_j $ for all $i, j$ \end{enumerate} gives us: \begin{equation} \mu_a = \frac{\sum_b w_b}{8 w_a} \end{equation} The correction term needs a small modification for cases where there is a octant which contains zero particles and hence zero weight: Let $q_{nz}$ be the number of octants with non-zero weight, i.e. with non-zero particles in them. Then \begin{equation} \mu_a = \frac{\sum_b w_b}{q_{nz} w_a} \label{eq:isotropy-correction-with-zero} \end{equation} this ultimately gives us \begin{align} \xi_s = \mu_a \ \psi(\x_{s} - \x_{g}, h(\x_{s})) \end{align} where $a$ is the octant in which the gas particle resides. In principle, we could also inject photon fluxes directly from stellar sources, just as we do with photon energy density. Test however have shown that best results are achieved by not injecting any net flux onto gas particles, so we don't do that in GEAR-RT. %----------------------------------------------------------------------- \subsection{Second Step: Transport}\label{chap:transport-step} %----------------------------------------------------------------------- %----------------------------------------------------------------------- \subsubsection{Solving the Transport Equations} %----------------------------------------------------------------------- In this step, we solve \begin{align} \DELDT{E_i} + \nabla \cdot \Fbf_i &= 0 \label{eq:homogenized-dEdt}\\ \DELDT{\Fbf_i} + c^2 \nabla \cdot \mathds{P}_i &= 0 \label{eq:homogenized-dFdt} \end{align} They are both in the form of a hyperbolic conservation law, so we can use the finite volume particle method given in eq.~\ref{eq:meshless-explicit}. We have the state vector $\U$ and flux tensor $\F$: \begin{align} \U = \begin{pmatrix} E_i \\ \Fbf_i \end{pmatrix} \quad , \quad \F = \begin{pmatrix} \Fbf_i \\ c^2 \mathds{P}_i \end{pmatrix} \end{align} The total sum $\sum_l \F_{kl} A_{kl}$ for each particle $k$ is collected during a neighbor interaction loop. Once the sum is accumulated, the final explicit update of the state is given by: \begin{align} \U_k(t + \Delta t_k) = \U_k (t) - \frac{1}{V_k} \sum\limits_{l} \min\{\Delta t_k, \Delta t_l\} \F_{\alpha,kl} \mathcal{A}_{kl}^\alpha \label{eq:transport-update-explicit} \end{align} $V_k$ are the associated particle volumes given in eq.~\ref{eq:volume}. \Aij are the effective surface terms given in eq.~\ref{eq:HopkinsAij}. It remains to determine the inter-particle fluxes $\F_{kl}$. We also follow the scheme outlined by the finite volume particle method. As it is done for the hydrodynamics, we approximate the problem by defining an ``interface'' at the position \begin{equation} \mathbf{x}_{kl} = \mathbf{x}_k + \frac{h_k}{h_k + h_l} ( \mathbf{x}_l - \mathbf{x}_k ) \end{equation} and extrapolate the value of the conserved variables at the position: \begin{equation} \U_{k, l}(\mathbf{x} = \mathbf{x}_{kl}) \approx \U_{k, l} + \DELDX{ \U_{k, l}}\ (\mathbf{x}_{kl} - \mathbf{x}_{k,l}) \label{eq:gradient-extrapolation} \end{equation} such that the fluxes $\F_{kl}$ in the update formula~\ref{eq:transport-update-explicit} between particles is estimated as the solution of the Riemann problem with ``left'' state $\U_k(\x_{kl})$ and ``right'' state $\U_l(\x_{kl})$ \begin{align} \F_{kl} = RP \left(\U_k + (\x_{k} - \x_{kl})_\alpha \DELDXALPHA{\U_k}, \ \U_l + (\x_{l} - \x_{kl})_\alpha \DELDXALPHA{\U_l} \right) \label{eq:rt-riemann-problem} \end{align} GEAR-RT offers two approximate Riemann solvers to find solutions to this Riemann problem. They are described in Section~\ref{chap:riemann-rt}. The gradients $\DELDX{\U_{k,l}}$ are estimated again using the least-squares second order accurate gradient expression given in eqn. \ref{eq:gradient}, for which a particle-particle interaction loop is required, which furthermore needs to be performed before the transport interaction loop, during which the sum of the exchanged (time integrated) fluxes of eq.~\ref{eq:transport-update-explicit} are accumulated. Since the RT is solved after the hydrodynamics and the particle positions haven't changed since the last hydrodynamics drift, the particles' smoothing lengths are still accurate, and we don't require an additional density loop for the radiative transfer. For that reason the gradient particle loop constitutes the first particle interaction loop in the RT scheme. %--------------------------------------------------- \subsubsection{Flux Limiters} %--------------------------------------------------- The gradient extrapolation (eq.~\ref{eq:gradient-extrapolation}) is equivalent from going from a piece-wise constant reconstruction of the radiation field to a piece-wise linear one, which makes the method second order accurate in space. For each quantity $Q_{k,l} \in \U = (E_{k,l}, \Fbf_{k,l})^T$, we extrapolate the value at the interface $\x_{kl}$ using the estimated gradient $dQ_{k,l}$. To prevent oscillations and instabilities from occurring, we need to employ a flux limiter at this point, which effectively reduces the gradient slopes $dQ_k$ and $dQ_l$ whenever necessary. The two limiters that worked well in tests are the \emph{minmod limiter} \begin{align} % \alpha_{\text{minmod}}(dQ_k, dQ_l) = % \begin{cases} % \text{if } dQ_k \times dQ_l > 0: % \begin{cases} % \text{if } |dQ_k| < |dQ_l| :\quad &\alpha = dQ_k / dQ_l \\ % \text{else: } &\alpha = dQ_l / dQ_k\\ % \end{cases}\\ % \text{else: } \alpha = 0 % \end{cases} dQ_{\text{minmod}}(dQ_k, dQ_l) = \begin{cases} dQ_k & \text{ if } |dQ_k| < |dQ_l| \text{ and } dQ_k dQ_l > 0 \\ dQ_l & \text{ if } |dQ_k| > |dQ_l| \text{ and } dQ_k dQ_l > 0 \\ 0 & \text{ if } dQ_k dQ_l \leq 0 \end{cases} \label{eq:rt-minmod} \end{align} and the \emph{monotonized central difference (MC) limiter}: \begin{align} r &= dQ_k / dQ_l \nonumber \\ \alpha_{\text{MC}}(r) &= \max \left\{ 0, \min\left[\frac{1 + r}{2}, 2, 2r \right] \right\} \label{eq:rt-MC}\\ dQ_{k, limited} &= \alpha_{MC}(r)\ dQ_k \nonumber\\ dQ_{l, limited} &= \alpha_{MC}(r)\ dQ_l \nonumber \end{align} The minmod limiter appears to be the most robust one with regards to ensuring stability for radiative transfer, and for this reason is the recommended choice. Other possible choices are the \emph{superbee} limiter: \begin{align} r &= dQ_k / dQ_l \nonumber \\ \alpha_{\text{superbee}}(r) &= \max \left\{ 0, \min (1, 2r), \min(2, r) \right\} \label{eq:rt-superbee}\\ dQ_{k,l,\ limited} &= \alpha_{superbee}(r)\ dQ_{k,l} \nonumber \end{align} And the \emph{van Leer} limiter: \begin{align} r &= dQ_k / dQ_l \nonumber\\ \alpha_{\text{vanLeer}}(r) &= \frac{r + |r|}{1 + |r|} \label{eq:rt-vanLeer}\\ dQ_{k,l,\ limited} &= \alpha_{vanLeer}(r)\ dQ_{k,l} \nonumber \end{align} %------------------------------------------------------------------------------------------- \subsubsection{Riemann Solvers for the Moments of the RT equation}\label{chap:riemann-rt} %------------------------------------------------------------------------------------------- It remains to find a solution for the Riemann problem~\ref{eq:rt-riemann-problem} which gives us the inter-particle (``inter-cell'') fluxes $\F_{kl}$. To find $\F_{kl}$, we use the extrapolated and flux limited states $E_{k,l}$ and $\Fbf_{k,l}$ at the interface position $\x_{kl}$ to compute the states and fluxes of the conservation law, $\U_{k,l}(\x_{kl})$ and $\F_{k,l}(\x_{kl})$. The components of the vector valued photon flux density $\Fbf_{k,l}$ are projected along the unit vector pointing from the particle towards the surface, allowing us to treat each component individually as a one dimensional problem. We then adapt the convention that particle $k$ is the left state of the Riemann problem, $\U_L$ and $\F_R$, while particle $l$ is the right state $\U_R$ and $\F_R$. The \emph{Global Lax Friedrich (GLF)} Riemann solver \citep{ramses-rt13} then gives the following solution for the flux $\F$ at an interface which separates a left state $\U_L$ and right state $\U_R$: \begin{equation} \F_{kl}(\U_L, \U_R) = \frac{\F_{L} + \F_{R}}{2} - \frac{c}{2} \left(\U_R - \U_L \right) \label{eq:riemann-GLF} \end{equation} which gives us an approximate solution for the flux $\F_{kl}$ we are looking for. An alternative is the \emph{Harten - Lax - van Leer (HLL)} Riemann solver (\cite{gonzalezHERACLESThreedimensionalRadiation2007}), given by: \begin{align} \F(\U_L, \U_R) &= \frac{ \lambda^{+} \F_{L} - \lambda^{-} \F_{R} + \lambda^{+} \lambda^{-} (\U_R - \U_L)}{ \lambda^{+} - \lambda^{-} } \label{eq:riemann-HLL} \\ \lambda^+ &= \max(0, \lambda_{max})\\ \lambda^- &= \min(0, \lambda_{min}) \end{align} Here, $\lambda^{max}$ and $\lambda^{min}$ are the minimum and the maximum of the Eigenvalues of the Jacobian matrix $\frac{\del \F(\U)}{\del \U}$. It turns out that the Eigenvalues can be described using only two parameters, the angle between the flux and the interaction interface, $\theta$, and the reduced flux $\mathbf{f} = \Fbf / (cE)$. Rather then continuously computing them on-the-fly for every interaction, they are tabulated\footnote{ A program to produce the tables of the Eigenvalues depending on $\theta$ and $\mathbf{f}$ is provided on \url{https://github.com/SWIFTSIM/swiftsim-rt-tools}. } and interpolated. The interest in the HLL solver is mainly because while it is more expensive than the GLF solver, it was shown \citep{ramses-rt13, gonzalezHERACLESThreedimensionalRadiation2007} to be less diffusive and is thus better suited to form shadows more correctly than the GLF solver, which is a known weakness of the moment based radiative transfer. %--------------------------------------------- \subsection{Third step: Thermochemistry} %--------------------------------------------- In this final step, we solve for the interaction between photons and the gas. The equations to be solved are: \begin{align*} \DELDT{E_i} &= - \sum\limits_{j}^{\absorbers} n_j \sigma_{i j} c E_i + \dot E_i ^{rec} \\ \DELDT{\Fbf_i} &= - \sum\limits_{j}^{\absorbers} n_j \sigma_{i j} c \Fbf_i \end{align*} The source term for recombination radiation, $\dot{E_i}^{rec}$, was added in for completeness. GEAR-RT currently (2023) neglects the emission of recombination. The actual thermochemistry is solved using the GRACKLE library\footnote{ GRACKLE is available on \url{https://github.com/grackle-project/grackle}.\\ Note that a frozen version of GRACKLE, tested and verified to work with GEAR-RT and other SWIFT projects, is being maintained on \url{https://github.com/mladenivkovic/grackle-swift}. It is strongly suggested to use that version of GRACKLE. } \citep{smithGrackleChemistryCooling2017}. This involves the evolution of the individual species number densities, as well as the heating of the gas. Currently (2023) GEAR-RT supports the ``6 species network'', composed of H$^0$, H$^+$, He$^0$, He$^+$, He$^{++}$, and $e^{-}$, provided by GRACKLE. GRACKLE requires us to provide it with the (total) photo-heating rate $\mathcal{H} = \sum_{i,j} \mathcal{H}_{ij}$ and the individual photo-ionization rates $\Gamma_{j} = \sum_i \Gamma_{i,j}$ for each species $j$ and photon group $i$. The expressions for the heating and photo-ionization rates for each photon group are given in eqns.~\ref{eq:photoheating-group} and \ref{eq:photoionization-group}. GRACKLE then solves the network of equations of thermochemistry for the 6 species and evolves the internal energy of the gas over some time step $\Delta t$. While GRACKLE evolves the state of the gas, we need to take care of the removal of the absorbed radiation energy over the thermochemistry time step manually. We solve this again using a simple first order forward Euler integration. A minor complication is that the number density of the ionizing species $n_j$ is not constant over the time step $\Delta t$, as we are currently in the process of ionizing the gas. This may lead to scenarios where the absorption rate is overestimated by not accounting for the ionization during the thermochemistry time step. To account for this effect, we take the average absorption rate over the time step instead: \begin{align*} \frac{E_i (t + \Delta t) - E_i}{\Delta t} = -\overline{\epsilon}_i \ N_i(t) c \ \sum_j \sigma_{i j} ^ N \ \frac{n_j(t + \Delta t) + n_j(t)}{2} \end{align*} which we can do since GRACKLE provides us with the updated value $n_j(t + \Delta t)$. The final equation used to remove absorbed photons from the radiation field is: \begin{align} E_i (t + \Delta t) = E_i(t) \left(1 - c \ \sum_j \sigma_{i j} ^ N \ \frac{n_j(t + \Delta t) + n_j(t)}{2} \ \Delta t \right) \end{align} The same is applied to the photon fluxes: \begin{align} \Fbf_i (t + \Delta t) &= \Fbf_i(t) \left(1 - c \ \sum_j \sigma_{i j} ^ N \ \frac{n_j(t + \Delta t) + n_j(t)}{2} \ \Delta t \right) \end{align} %--------------------------------------------- \subsection{Additional Topics} %--------------------------------------------- %----------------------------------------------------- \subsubsection{Computing the Time Step Size} %----------------------------------------------------- The time step size of the particles is computed in the same manner as we do it for hydrodynamics. A ``cell size'' $\Delta x$ of a particle is estimated using \begin{align} \Delta x \approx \left(\frac{V_i}{V_{S,\nu}} \right)^{1/\nu} \end{align} where $V_{S,\nu}$ is the volume of a $\nu$-dimensional unit sphere, i.e. $4 \pi / 3$ in 3D. The signal velocity in case of radiation is simply the speed of light $c$, which leads to the time step estimate \begin{align} \Delta t = C_{CFL} \frac{\Delta x}{c} \label{eq:rt-cfl} \end{align} %--------------------------------------------- \subsubsection{Dealing With Particle Drifts}\label{chap:rt-drift} %--------------------------------------------- While we assume particles are static w.r.t. the simulated volume in the context of RT, they are drifted for the purposes of hydrodynamics. As a consequence, we need to make corrections to the radiation fields when a particle is drifted for hydrodynamics purposes. We do that by simply extrapolating the particle value at the new point given its current state and the gradients $\nabla \U$ we obtained using the general gradient expression given in eq~\ref{eq:gradient}. Explicitly, for a particle $i$ and for each component of the radiation state vector $\U^k$, we first obtain the drift distance over the time $\Delta t_{i, drift}$ as \begin{align} \Delta \x_{i,\text{drift}} &= \V_{i} \Delta t_{i,\text{drift}} \end{align} where $\V_i$ is the particle velocity. We then obtain the value corrected at the drifted location \begin{align} \U_{i}^k (\x + \Delta \x_{\text{drift}}) \approx \U_{i}^k (\x) + \nabla \U_{i}^k (\x) \cdot \Delta x_{i,\nu,\text{drift}} \end{align} %--------------------------------------------- \subsubsection{Cosmology} %--------------------------------------------- TODO %--------------------------------------------- \subsubsection{Reducing The Speed of Light} %--------------------------------------------- In an attempt to reduce the computational expense associated with these tiny time step sizes, \cite{ramses-rt13} employ the Reduced Speed of Light Approximation (RSLA), which was introduced by \cite{gnedinMultidimensionalCosmologicalRadiative2001}. I also adapt the RSLA with GEAR-RT. The approach consists of globally reducing the speed of light everywhere by some user-set factor $f_c$, i.e. replace the speed of light $c$ in all equations involved with radiative transfer by \begin{align} \tilde{c} = f_c c \end{align} %--------------------------------------------- \subsubsection{Ion Mass Fluxes} %--------------------------------------------- Hydrodynamics with Finite Volume Particle Methods exchange mass fluxes $\Delta m$ between particles. This means that we need to pay attention to the individual ionizing and ionized species' mass fractions being exchanged. The mass of each individual species being exchanged needs to be traced individually, and depending on the direction: If a particle loses mass, then it loses mass and species mass fractions according to its own mass fractions. If it gains mass during an exchange, then it gains species mass fractions according to the interaction partner's mass fractions. The usual convention for a mass flux between particles ``$i$'' and ``$j$'' is to treat $i$ as the ``left'' particle and $j$ as the ``right'' particle in an analogy to cells. Let $X_k$ denote the mass fraction of each species $k$. During each \emph{hydrodynamics} flux exchange, for each species $k$ the total mass flux $\Delta x_{i,k}$ is being accumulated: \begin{align} \Delta x_{i,k} &= \sum_j \Delta m_{ij} X_{ij,k}\\ \text{where }\quad X_{ij,k} &= \begin{cases} X_{i,k} \quad \text{ if } \quad \frac{\Delta m_{ij}}{\Delta t} > 0 \\ X_{j,k} \quad \text{ if } \quad \frac{\Delta m_{ij}}{\Delta t} < 0 \end{cases} \end{align} The masses which particles carry are updated during the kick operation, during which the individual species' masses $x_{i,k}$ are also evolved. For a particle $i$ with mass $m_i$, the masses of individual species at the beginning of a step are given by: \begin{align} x_{i,k}^{t} &= m_i X_{i,k}^{t} \end{align} During the kick operation, they get updated as follows: \begin{align} x_{i,k}^{t + \Delta t} &= x_{i,k}^{t} + \Delta x_{i,k} \end{align} Finally, after the update, the new mass fractions of the species can be obtained using \begin{align} X_{i,k}^{t + \Delta t} &= \frac{x_{i,k}^{t + \Delta t}}{\sum_k x_{i,k}^{t + \Delta t}} \end{align}