diff --git a/.gitignore b/.gitignore
index 2f2b42d95d656022649b3b77ef5c07e914cf21af..8feacfd8ae84bb9dccfd8a021ea8a5314715fc68 100644
--- a/.gitignore
+++ b/.gitignore
@@ -222,6 +222,7 @@ theory/Multipoles/mac_potential.pdf
theory/Cosmology/cosmology.pdf
theory/Cooling/eagle_cooling.pdf
theory/Gizmo/gizmo-implementation-details/gizmo-implementation-details.pdf
+theory/RadiativeTransfer/GEARRT/GEARRT.pdf
m4/libtool.m4
m4/ltoptions.m4
diff --git a/theory/RadiativeTransfer/GEARRT/1-meshless.tex b/theory/RadiativeTransfer/GEARRT/1-meshless.tex
new file mode 100644
index 0000000000000000000000000000000000000000..f1e3b95d2f1baef6ac0bef6980038c58cca40cee
--- /dev/null
+++ b/theory/RadiativeTransfer/GEARRT/1-meshless.tex
@@ -0,0 +1,115 @@
+%===================================================================================
+\section{Finite Volume Particle Method: Essentials} \label{chap:meshless}
+%===================================================================================
+
+The Finite Volume Particle methods (FVPM) are a way to solve hyperbolic systems of conservation
+laws of the form
+
+\begin{equation}
+ \DELDT{\mathcal{U}} + \nabla \mathcal{F}(\mathcal{U}) = 0 \label{eq:conservation-law}
+\end{equation}
+
+for some state vector $\mathcal{U}$ and flux tensor $\mathcal{F}(\mathcal{U})$
+and are used to solve the hydrodynamics equations. Since the M1 RT method can be written in the same
+form (see Section~\ref{chap:rt-equations}), the method will also be used to solve these equations.
+For completeness and quick reference, the most important formulae are listed below. For the actual
+implementation details, please refer to the `Gizmo' theory section in the SWIFT repository, and
+\citet{hopkinsGIZMONewClass2015}. For a step-by-step derivation and introduction, see
+\citet{ivkovicThesis}.
+
+
+
+The FVPM schemes rely on the ``partitions of unity''. The partition of unity assigned to a particle
+$i$ at some point $\x$ is defined as:
+
+\begin{align}
+ \psi_i(\x) &= \frac{1}{\omega(\x)} W(\x - \x_i, h(\x)) \label{eq:psi} \\
+ \omega(\x) &= \sum_j W(\x - \x_j, h(\x)) \label{eq:omega}
+\end{align}
+
+where $h(\x)$ is the smoothing length at $\x$ and $\omega(\x)$ is used to normalise the volume
+partition at any point $\x$.
+
+The associated particle volumes are given by
+
+\begin{align}
+ V_i &= \int_V \psi_i(\x) \de V \approx \frac{1}{\omega(\x_i)} \label{eq:volume}
+\end{align}
+
+
+The actual state updates are given by
+
+\begin{equation}
+ \frac{\D}{\D t} (V_i \mathcal{U}_{k,i}) + \sum_j \mathcal{F}_{k,ij} \cdot \Aijm = 0
+\end{equation}
+
+
+with
+
+\begin{equation}
+ \Aijm^\alpha = V_i \psitilde_j^\alpha (\x_i) - V_j\psitilde_i^\alpha (\x_j) \label{eq:HopkinsAij}
+\end{equation}
+
+for every component $k$ of the state vector $\mathcal{U}$ and for every coordinate dimension
+$\alpha$.
+
+The $\psitilde(\x)$ come from the $\mathcal{O}(h^2)$ accurate discrete gradient expression from
+\cite{lansonRenormalizedMeshfreeSchemes2008}:
+
+\begin{align}
+ \frac{\del}{\del x_{\alpha}} f(\x) \big{|}_{\x_i} &=
+ \sum_j \left( f(\x_j) - f(\x_i) \right) \psitilde_j^\alpha (\x_i) \label{eq:gradient} \\
+ \psitilde_j^\alpha (\x_i) &= \sum_{\beta = 1}^{\beta=\nu} \mathcal{B}_i^{\alpha \beta}
+ (\x_j - \x_i)^\beta \psi_j(\x_i) \\
+ \mathcal{B}_i &= \mathcal{E_i} ^ {-1} \label{eq:matrix-B}\\
+ \mathcal{E}_i^{\alpha \beta} &= \sum_j (\x_j - \x_i)^\alpha (\x_j - \x_i)^\beta \psi_j(\x_i)
+\label{eq:matrix-E}
+\end{align}
+
+
+where $\alpha$ and $\beta$ again represent the coordinate components for $\nu$ dimensions.
+
+
+The final explicit update formula is given by
+
+\begin{align}
+ \U_i^{n+1} = \U_i^{n} + \frac{\Delta t}{V_i} \sum_j \F_{\alpha, ij} \Aijm^\alpha
+\label{eq:meshless-explicit}
+\end{align}
+
+
+
+
+
+
+
+
+%===================================================================================
+\subsection{Individual Timestepping}
+%===================================================================================
+
+With particles being allowed to have individual time steps, the update formula for the finite volume
+particle methods (eq.~\ref{eq:meshless-explicit}) needs to be slightly updated in order to
+maintain its conservative properties. In particular, rather than exchanging fluxes, particles now
+need to exchange time integrated fluxes, which are integrated using a trivial explicit Euler scheme:
+
+\begin{align}
+\U_i^{n+1} =
+ \U_i^n + \frac{1}{V_i} \sum_j \left[
+ \min\{ \Delta t_i, \Delta t_j \} \F_{\alpha, ij} \Aijm^\alpha
+ \right]
+\label{eq:meshless-explicit-individual-timesteps}
+\end{align}
+
+The time integration needs to be the smaller of the two time step sizes $\{\Delta t_i,\ \Delta
+t_j\}$ since if one particle has a smaller time step size than the other, it also means that there
+will be several interactions between that particle pair, each with the time step size of the smaller
+of the two. Suppose for example $\Delta t_i = 4 \Delta t_j$. Then for a single update of particle
+$i$, there will be in total 4 updates of particle $j$. Each update of particle $j$ will entail a
+symmetric (time integrated) flux exchange with particle $i$. Therefore the time integration of the
+flux to be exchanged must be applied over the smaller time interval $\Delta t_j$ for each of
+the four (time integrated) flux exchanges.
+
+
+
+
diff --git a/theory/RadiativeTransfer/GEARRT/2-rt-equations.tex b/theory/RadiativeTransfer/GEARRT/2-rt-equations.tex
new file mode 100644
index 0000000000000000000000000000000000000000..ff8e1335ad4083c777185d4d853fedf733fac0b2
--- /dev/null
+++ b/theory/RadiativeTransfer/GEARRT/2-rt-equations.tex
@@ -0,0 +1,224 @@
+%=============================================================================
+\section{The Equations of Moment-Based Radiative Transfer}\label{chap:rt-equations}
+%=============================================================================
+
+
+%-----------------------------------------------------------------------
+\subsection{The Equations of Radiative Transfer and the M1 Closure}
+%-----------------------------------------------------------------------
+
+
+%-----------------------------------------------------------------------
+\subsubsection{The Equations of Radiative Transfer}
+%-----------------------------------------------------------------------
+
+
+Radiative transfer (RT) and radiation hydrodynamics (RHD) contain a plethora of variables and
+coefficients. For clarity, an overview of the relevant quantities and coefficients is given
+in Appendix~\ref{app:variables} along with their respective units.
+
+
+The equation of radiative transfer is given by:
+
+\begin{align}
+ \frac{1}{c} \DELDT{I_\nu} + \mathbf{n} \cdot \nabla I_\nu
+ &= \eta_\nu - \alpha_\nu I_\nu \nonumber \\
+ &= \eta_\nu - \sum_j^{\text{photo-absorbing\ species}} \sigma_{j,\nu} n_j I_\nu \ .
+\label{eq:RT-sigma}
+\end{align}
+
+\begin{itemize}
+\item $I_\nu$ is the specific intensity and has units of erg cm$^{-3}$ rad$^{-1}$ Hz$^{-1}$
+s$^{-1}$.
+\item $\eta_\nu$ is a source function of radiation, i.e. the term describing radiation being added
+along the (dimensionless) direction $\mathbf{n}$ due to yet unspecified processes, and has units of
+erg cm$^{-3}$ rad$^{-1}$ Hz$^{-1}$ s$^{-1}$, which is the same as the units of the specific
+intensity $I_\nu$ per cm.
+\item $\alpha_\nu$ is an absorption coefficient, describes how much radiation is being removed, and
+has units of cm$^{-1}$. Naturally only as much radiation as is currently present can be removed, and
+so the sink term must be proportional to the local specific intensity $I_\nu$.
+\end{itemize}
+
+The equation holds for any photon frequency $\nu$ individually. In eq.~\ref{eq:RT-sigma} we split
+the absorption coefficient $\alpha_\nu$ into the sum over the photo-absorbing species $j$, which
+for GEAR-RT will only be the main constituents of primordial gas, namely hydrogen, helium,
+and singly ionized helium. The photo-absorption process is expressed via interaction cross sections
+$\sigma_{j,\nu}$, which has units of cm$^2$, while $n_j$ represents the number density of
+photo-absorbing species $j$ in cm$^{-3}$.
+
+
+
+%-----------------------------------------------------------------------
+\subsubsection{Moments of The Equations of Radiative Transfer}
+%-----------------------------------------------------------------------
+
+
+GEAR-RT solves for the (angular) moments of the equation of radiative transfer, which are obtained
+through integrating eq.~\ref{eq:RT-sigma} over the entire solid angle to obtain the zeroth moment
+equation, and over the entire solid angle multiplied by the direction unit vector $\mathbf{n}$ for
+the first moment equation.
+Additionally, we make use of the following quantities:
+
+\begin{align*}
+ E_\nu (\x, t) &= \int_{4 \pi} \frac{I_\nu}{c} \de \Omega
+ && \text{total energy density }
+ && [E_\nu] = \frac{\text{erg}}{\text{cm}^3 \text{ Hz}}\\
+ \Fbf_\nu(\x, t) &= \int_{4 \pi} I_\nu \mathbf{n} \de \Omega
+ && \text{radiation flux }
+ && [\Fbf_\nu] = \frac{\text{erg}}{\text{cm}^2 \text{ s Hz}}\\
+ \mathds{P}_\nu (\x, t) &= \int_{4 \pi} \frac{I_\nu}{c} \mathbf{n} \otimes \mathbf{n} \de \Omega
+ && \text{radiation pressure tensor }
+ && [\mathds{P}_\nu ] = \frac{\text{erg}}{\text{cm}^3 \text{ Hz}}
+\end{align*}
+
+where $\mathbf{n} \otimes \mathbf{n}$ denotes the outer product, which in components $k$, $l$ gives
+%
+\begin{align*}
+ (\mathbf{n} \otimes \mathbf{n})_{kl} = \mathbf{n}_k \mathbf{n}_l \ .
+\end{align*}
+
+
+
+This gives us the following equations:
+
+
+\begin{align}
+ \DELDT{E_\nu} + \nabla \cdot \Fbf_\nu &=
+ - \sum\limits_{j}^{\absorbers} n_j \sigma_{\nu j} c E_\nu + \dot{E}_\nu
+ \label{eq:dEdt-freq} \\
+ \DELDT{\Fbf_\nu} + c^2 \ \nabla \cdot \mathds{P}_\nu &=
+ - \sum\limits_{j}^{\absorbers} n_j \sigma_{\nu j} c \Fbf_\nu
+ \label{eq:dFdt-freq}
+\end{align}
+
+
+
+Note that $E_\nu$ (and $\dot{E}_\nu$) is the radiation energy \emph{density} (and \emph{density}
+injection rate) in the frequency interval between frequency $\nu$ and $\nu + \de \nu$ and has
+units of $\text{erg / cm}^3 \text{ / Hz}$ (and $\text{erg / cm}^3 \text{ / Hz / s}$).
+$\Fbf$ is the radiation flux, and has units of $\text{erg / cm}^2 \text{ / Hz / s}$, i.e.
+dimensions of energy per area per frequency per time.
+
+Furthermore, it is assumed that the source term $\dot{E}_\nu$ stems from point sources which radiate
+isotropically. This assumption has the consequence that the vector net flux $\Fbf_\nu$ must sum up
+to zero, and hence the corresponding source terms in eq.~\ref{eq:dFdt-freq} are zero.
+
+
+
+%-----------------------------------------------------------------------
+\subsubsection{The M1 Closure}
+%-----------------------------------------------------------------------
+
+
+
+To close this set of equations, a model for the pressure tensor $\mathds{P}_\nu$ is necessary. We
+use the so-called ``M1 closure'' \citep{levermoreRelatingEddingtonFactors1984a} where we describe
+the pressure tensor via the Eddington tensor $\mathds{D}_\nu$:
+
+
+\begin{equation*}
+ \mathds{P}_\nu = \mathds{D}_\nu E_\nu
+\end{equation*}
+
+The Eddington tensor is a dimensionless quantity that encapsulates the local radiation field
+geometry and its effect in the radiation flux conservation equation. The M1 closure sets the
+Eddington tensor to have the form:
+
+\begin{align*}
+\mathds{D}_\nu &=
+ \frac{1- \chi_\nu}{2} \mathds{I} + \frac{3 \chi_\nu - 1}{2} \mathbf{n}_\nu \otimes
+ \mathbf{n}_\nu
+ \label{eq:eddington-freq} \\
+\mathbf{n}_\nu &=
+ \frac{\Fbf_\nu}{|\Fbf_\nu|} \\
+\chi_\nu &=
+ \frac{3 + 4 f_\nu ^2}{5 + 2 \sqrt{4 - 3 f_\nu^2}} \\
+f_\nu &=
+ \frac{|\Fbf_\nu|}{c E_\nu}
+\end{align*}
+
+
+
+
+
+
+
+
+
+
+
+
+%---------------------------------------------------------------
+\subsection{Photo-ionization and Photo-heating Rates}
+%---------------------------------------------------------------
+
+
+In the context of radiative transfer and photo-ionization, the photo-ionization rate $\Gamma_{\nu,
+j}$ in units of s$^{-1}$ for photons with frequency $\nu$ and a photo-ionizing particle species $j$
+is then given by
+
+\begin{align}
+ \DELDT{n_j} = -\Gamma_{\nu, j} \ n_j = - c \ \sigma_{\nu j} \ N_\nu \ n_j
+\end{align}
+
+where $N_\nu = E_\nu / (h \nu)$ is the photon number density, and $n_j$ is the number density of
+the photo-ionizing particle species. Note that both the interaction cross sections and the
+photo-ionizing species $j$ are specific to a frequency $\nu$.
+
+For the cross sections, we use the analytic fits for the photo-ionization cross sections from
+\cite{vernerAtomicDataAstrophysics1996} (via \cite{ramses-rt13}), which are given by
+
+\begin{align}
+\sigma(E) &= \sigma_0 F(y) \times 10^{-18} \text{ cm}^2 \label{eq:sigma-parametrizaiton}
+\\
+F(y) &= \left[(x - 1)^2 + y_w^2 \right] y ^{0.5 P - 5.5} \left( 1 + \sqrt{y / y_a} \right)^{-P}
+\\
+x &= \frac{E}{E_0} - y_0 \\
+y &= \sqrt{x^2 + y_1^2}
+\end{align}
+
+where $E$ is the photon energy $E = h \nu$ in eV, and $\sigma_0$, $E_0$, $y_w$, $y_a$, $P$, $y_0$,
+and $y_1$ are fitting parameters. The fitting parameter values for hydrogen, helium, and singly
+ionized helium are given in Table~\ref{tab:cross-sections}. The thresholds are given as frequencies
+in eqs.~\ref{eq:nuIonHI}-\ref{eq:nuIonHeII}. Below this threshold, no ionization can take place, and
+hence the cross sections are zero.
+
+
+\input{tables/fitting_parameters}
+
+
+
+
+
+
+Conversely, the rate at which photons are absorbed, i.e. ``destroyed'', must be equal to the
+photo-ionization rate, which means
+
+\begin{align*}
+\DELDT{E_\nu} &= h \ \nu \DELDT{N_\nu} = -h\ \nu \ c \ \sigma_{\nu j} n_j N_\nu
+\end{align*}
+
+
+
+Finally, the photo-heating rate is modeled as the rate of excess energy absorbed by the gas during
+photo-ionizing collisions. To ionize an atom, the photons must carry a minimal energy corresponding
+to the ionizing frequency $\nu_{ion,j}$ for a photo-ionizing species $j$. In the case of hydrogen
+and helium, their values are
+
+\begin{align}
+ \nu_{\text{ion,HI}} &= 2.179 \times 10^{-11} \text{ erg} = 13.60 \text{ eV} \label{eq:nuIonHI}\\
+ \nu_{\text{ion,HeI}} &= 3.940 \times 10^{-11} \text{ erg} = 24.59 \text{ eV}
+\label{eq:nuIonHeI}\\
+ \nu_{\text{ion,HeII}} &= 8.719 \times 10^{-11} \text{ erg} = 54.42 \text{ eV}
+\label{eq:nuIonHeII}
+\end{align}
+
+All excess energy is added to the gas in the form of internal energy, and the heating rate
+$\mathcal{H}$ (in units of erg cm$^{-3}$ s$^{-1}$ ) for photons of some frequency $\nu$ and some
+photo-ionizing species $j$ is described by
+
+\begin{align*}
+\mathcal{H}_{\nu, j} = (h \nu - h \nu_{ion,j}) \ c \ \sigma_{\nu j} \ n_j \ N_\nu
+\end{align*}
+
+
diff --git a/theory/RadiativeTransfer/GEARRT/3-rt-numerical-method.tex b/theory/RadiativeTransfer/GEARRT/3-rt-numerical-method.tex
new file mode 100644
index 0000000000000000000000000000000000000000..a3350bc7a62bb00ada5d9c0e7922304bfc094274
--- /dev/null
+++ b/theory/RadiativeTransfer/GEARRT/3-rt-numerical-method.tex
@@ -0,0 +1,857 @@
+%======================================================
+\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}
+
+
+
diff --git a/theory/RadiativeTransfer/GEARRT/GEARRT.tex b/theory/RadiativeTransfer/GEARRT/GEARRT.tex
new file mode 100644
index 0000000000000000000000000000000000000000..e90a2aad2a0524624763079e47f326395a4c66fe
--- /dev/null
+++ b/theory/RadiativeTransfer/GEARRT/GEARRT.tex
@@ -0,0 +1,57 @@
+\input{./header}
+
+
+
+%------------------------------------------
+%: Metadata
+%------------------------------------------
+
+\title{GEAR-RT}
+\subtitle{Theory Documentation}
+\author{Mladen Ivkovic}
+\date{September 2023}
+
+%------------------------------------------
+
+
+
+
+
+
+\begin{document}
+
+
+%Titlepage
+\maketitle
+\tableofcontents
+\clearpage
+
+
+
+
+
+\include{1-meshless}
+\include{2-rt-equations}
+\include{3-rt-numerical-method}
+
+
+
+
+
+
+
+
+\clearpage
+\appendix{}
+\include{appendix}
+
+
+%------------
+%Bibliography
+%------------
+\clearpage
+\bibliography{references}
+
+
+
+\end{document}
diff --git a/theory/RadiativeTransfer/GEARRT/appendix.tex b/theory/RadiativeTransfer/GEARRT/appendix.tex
new file mode 100644
index 0000000000000000000000000000000000000000..b739b91ac701c983549c20e6a54548671d4fb381
--- /dev/null
+++ b/theory/RadiativeTransfer/GEARRT/appendix.tex
@@ -0,0 +1,288 @@
+%------------------------------------------------
+\section{RT Related Quantities And Their Units}\label{app:variables}
+%------------------------------------------------
+
+It could be helpful to have the commonly used quantities of RT written down
+somewhere along with their units. So here you go.
+
+
+\include*{./tables/rt_variables}
+
+
+
+
+
+
+
+\clearpage
+%--------------------------------------------
+\section{Function Calls of an RT step}
+%--------------------------------------------
+
+
+
+\begin{landscape}
+{\footnotesize
+
+\begin{tabular}[l]{%
+ >{\raggedright\arraybackslash}p{2.6cm}%
+ >{\raggedright\arraybackslash}p{2.8cm}%
+ >{\raggedright\arraybackslash}p{7cm}%
+ >{\raggedright\arraybackslash}p{7cm}%
+}
+\textbf{task} & \textbf{task type} & \textbf{task purpose} & \textbf{function calls} \\[.5em]
+\hline
+\hline
+Injection Prep &
+ Interacting star and gas particles &
+ Gather gas particle neighbour data in preparation for the injection &
+ \texttt{runner\_iact\_rt\_inject\_prep} in \verb|src/rt/method/rt_iact.h| \\
+\hline
+Star Emission Rates &
+ Work on individual star particles &
+ Prepare everything necessary that needs to be done for radiation sources before the radiation sources and the gas interact. &
+ \texttt{rt\_compute\_stellar\_emission\_rate} in \verb|src/rt/method/rt.h| \\
+\hline
+\hline
+RT in &
+ Implicit&
+ Collect dependencies &
+ - \\
+\hline
+Injection &
+ Interacting star and gas particles &
+ Distribute the radiation from star particles onto gas particles &
+ \verb|runner_iact_rt_inject| in \verb|src/rt/method/rt_iact.h| \\
+\hline
+Ghost1 &
+ Work on individual gas particles &
+ Finish up any work that needs to be done for gas particles before the next gas $\leftrightarrow$ gas interaction begins &
+ \texttt{rt\_finalise\_injection} in \verb|src/rt/method/rt.h|\\
+\hline
+Gradient &
+ Interacting gas and gas particles &
+ Compute necessary gradients of the radiation quantities &
+ \verb|rt_gradients_collect| in \verb|src/rt/method/rt_gradients.h| and
+ \verb|rt_gradients_nonsym_collect| in \verb|src/rt/method/rt_gradients.h|\\
+\hline
+Ghost2 &
+ Work on individual gas particles &
+ Finish up any work that needs to be done for gas particles before the next gas $\leftrightarrow$ gas interaction begins &
+ \texttt{rt\_end\_gradient} in \verb|src/rt/method/rt.h|\\
+\hline
+Transport &
+ Interacting gas and gas particles &
+ Compute/exchange fluxes of homogenized equation of radiative transfer. &
+ \verb|runner_iact_rt_transport| in \verb|src/rt/method/rt_iact.h| and
+ \verb|runner_iact_nonsym_rt_transport| in \verb|src/rt/method/rt_iact.h|\\
+\hline
+Transport out &
+ Implicit&
+ Collect dependencies &
+ - \\
+\hline
+Thermochemistry &
+ Work on individual gas particles &
+ Finish up any work that needs to be done for gas particles before the thermochemistry part of the computation can begin. Then do the thermochemistry computation. &
+ \texttt{rt\_finalise\_transport} in \verb|src/rt/method/rt.h|,
+ \verb|rt_do_thermochemistry| in \verb|src/rt/method/rt_thermochemistry.h|\\
+\hline
+RT out &
+ Implicit&
+ Collect dependencies &
+ - \\
+\hline
+\end{tabular}
+
+
+}
+\end{landscape}
+
+
+
+
+
+
+
+
+
+\clearpage
+%-------------------------------------------------------------------------
+\section{Creating Collisional Ionization Equilibrium Initial Conditions}
+%-------------------------------------------------------------------------
+
+
+On startup, GEARRT offers the possibility to generate the ionization mass fractions of the gas
+particles assuming the gas is in collisional ionization equilibrium, composed of Hydrogen and
+Helium, and that there is no radiation present. In order to determine the ionization mass fractions
+of all species (H$^0$, H$^+$, He$^0$, He$^+$, He$^{++}$) for a given specific internal energy $u$,
+an iterative procedure needs to be applied because the gas variables are interconnected:
+
+
+
+Commonly the average particle mass $\overline{m}$ of a gas is expressed using the (unitless) mean
+molecular weight $\mu$ :
+
+\begin{align}
+ \overline{m} = \mu m_u
+\end{align}
+
+where $m_u$ is the atomic mass unit and $\mu$ is given by
+
+\begin{align}
+ \frac{1}{\mu} = \sum_j \frac{X_j}{A_j} (1 + E_j)
+\end{align}
+
+
+Specifically, ionization changes the mass fractions $X_j$ of the species, and therefore also the
+mean molecular weight $\mu$. In turn, the mean molecular weight determines the gas temperature at a
+given specific internal energy:
+
+\begin{align}
+ T = u (\gamma - 1) \mu \frac{m_u}{k}
+\end{align}`
+
+Lastly, the gas temperature determines the collisional ionization and recombination rates, which
+need to be balanced out by the correct number density of the individual species in order to be in
+ionization equilibrium, i.e. at a state where the ionization and recombination rates exactly cancel
+each other out.
+
+We take the ionization and recombination rates from \citet{katzCosmologicalSimulationsTreeSPH1996},
+which are given in Table \ref{tab:coll-ion-rates-katz}. For a gas with density $\rho$, Hydrogen
+mass fraction $X_H$ and Helium mass fraction $X_{He} = 1 - X_H$, the total number densities of all
+hydrogen and helium species are
+
+\begin{align}
+ n_H &= X_H \frac{\rho}{m_u} \\
+ n_{He} &= X_{He} \frac{\rho}{4 m_u}
+\end{align}
+
+and in equilibrium, the number densities of the individual species are given by
+
+\begin{align}
+ n_{H^0} &=
+ n_H \frac{A_{H^+}}{A_{H^+} + \Gamma_{H^0}} \\
+ n_{H^+} &=
+ n_H - n_{H^0} \\
+ n_{He^+} &=
+ n_{He} \frac{1}{1 + (A_{He^+} + A_d) / \Gamma_{He^0} + \Gamma_{He^+} / A_{He^{++}}} \\
+ n_{He^0} &=
+ n_{He^+} \frac{A_{He^+} + A_d}{\Gamma_{He^0}} \\
+ n_{He^{++}} &=
+ n_{He^+} \frac{\Gamma_{He^+}}{A_{He^+}}
+\end{align}
+
+
+To summarize, the tricky bit here is that the number densities determine the mean molecular weight,
+the mean molecular weight determines the temperature of the gas for a given density and specific
+internal energy, while the temperature determines the number densities of the species.
+
+To find the correct mass fractions, the iterative Newton-Raphson root finding method is used.
+Specifically, using some initial guesses for temperature and mean molecular weights, $T_{guess}$
+and $\mu_{guess}$, in each iteration step we determine the resulting specific internal energy
+$u_{guess} = k T_{guess} / (\gamma - 1) / (\mu_{guess} m_u)$. The function whose root we're looking
+for is
+
+\begin{align}
+ f(T) = u - u_{guess}(T)
+\end{align}
+
+with the derivative
+
+\begin{align}
+ \frac{\del f}{\del T} = - \frac{\del u}{\del T} (T = T_{guess}) = \frac{k}{(\gamma -
+1) / (\mu_{guess} m_u )}
+\end{align}
+
+where $u$ is the specific internal energy of the gas, which is fixed and provided by the initial
+conditions. We now look for the $T$ at which $f(T) = 0$. The Newton-Raphson method prescribes to
+find the $n+1$th $T_{guess}$ using
+
+\begin{align}
+ T_{n+1} = T_n + \frac{f(T_n)}{\frac{\del f}{\del T}(T_n)}
+\end{align}
+
+During each iteration, the new mass fractions and the resulting mean molecular weight given the
+latest guess for the temperature are computed. At the start, the first guess for the temperature
+$T_{guess}$ is computed assuming a fully neutral gas. Should that gas temperature be above $10^5$ K,
+the first guess is changed to a fully ionized gas. The iteration is concluded once $f(T) \leq
+\epsilon = 10^{-4}$.
+
+
+
+\include*{tables/thermochemistry_rates}
+
+
+
+
+
+
+
+
+
+\clearpage
+%-------------------------------------------------------------------------------
+\section{Converting Photon Number Emission Rates to Photon Energy Emission Rates}
+%-------------------------------------------------------------------------------
+
+Many other, in particular older, codes and papers use photon \emph{number} injection rates $\dot{N}_{\gamma}$ for their emission rates rather than the \emph{energy} injection rates $\dot{E}_\gamma$ or equivalently luminosities $L$.
+
+To convert between these two quantities, we need to assume that the emission follows some spectrum $J(\nu)$.
+
+In the case of a single photon group, the conversion is quite simple: We first need to compute the average photon energy $\overline{E}_\gamma$:
+
+\begin{align*}
+ \overline{E}_\gamma = \frac{\int J(\nu) \ \de \nu }{\int J(\nu) / (h \nu) \ \de \nu}
+\end{align*}
+
+then the emitted luminosity (energy per unit time) is
+
+\begin{align}
+ L = \overline{E}_\gamma \ \dot{N}_{\gamma}
+\end{align}
+
+Note that in many cases, the given emission photon number rate is the number rate of \emph{ionizing} photons. For us, this means that we need to start the integrals at the lowest ionizing frequency $\nu_{\text{ion, min}}$ in order to have the correct translation to the luminosity of the \emph{ionizing} energy:
+
+
+\begin{align}
+ \overline{E}_\gamma = \frac{\int\limits_{\nu_{\text{ion, min}}}^\infty J(\nu) \ \de \nu }{\int\limits_{\nu_{\text{ion, min}}}^\infty J(\nu) / (h \nu) \ \de \nu}
+\end{align}
+
+
+In the case of several photon groups being used, the conversion requires a little adaptation in order to preserve the correct number of photons emitted. For each photon group $i$, the average photon energy is given by
+
+\begin{align}
+ \overline{E}_i = \frac{\int\limits_{\nu_{i \text{, min}}}^{\nu_{i \text{, max}}} J(\nu) \ \de \nu }{\int\limits_{\nu_{i \text{, min}}}^{\nu_{i \text{, max}}} J(\nu) / (h \nu) \ \de \nu}
+\end{align}
+
+
+Secondly, we need to compute the fraction $f_i$ of ionizing photons in each bin, which is given by
+
+\begin{align}
+ f_i = \frac{\int\limits_{\nu_{i \text{, min}}}^{\nu_{i \text{, max}}} J(\nu) / (h \nu) \ \de \nu }{\int\limits_{\nu_{min}}^{\infty} J(\nu) / (h \nu) \ \de \nu}
+\end{align}
+
+Then the number of emitted photons in each bin is given by
+
+\begin{align}
+\dot{N}_i = f_i\ \dot{N}_\gamma
+\end{align}
+
+And the luminosities are given by
+
+\begin{align}
+ L_i &= \overline{E}_i \ \dot{N}_i \\
+ &= \frac{\int\limits_{\nu_{i \text{, min}}}^{\nu_{i \text{, max}}} J(\nu) \ \de \nu }{\int\limits_{\nu_{min}}^{\infty} J(\nu) / (h \nu) \ \de \nu} \ \dot{N}_\gamma
+\end{align}
+
+
+
+Python scripts to compute and convert photon number rates into luminosities are provided in
+\url{https://github.com/SWIFTSIM/swiftsim-rt-tools}.
+
+
+
+
+
+
+
diff --git a/theory/RadiativeTransfer/GEARRT/header.tex b/theory/RadiativeTransfer/GEARRT/header.tex
new file mode 100644
index 0000000000000000000000000000000000000000..4d70ae5486be35478305da801139eca81a3b59b1
--- /dev/null
+++ b/theory/RadiativeTransfer/GEARRT/header.tex
@@ -0,0 +1,204 @@
+\documentclass[12pt, a4paper, english, singlespacing, parskip]{scrartcl}
+
+%\documentclass[
+%11pt, % The default document font size, options: 10pt, 11pt, 12pt
+%oneside, % Two side (alternating margins) for binding by default, uncomment to switch to one side
+%chapterinoneline, % Have the chapter title next to the number in one single line
+%english, % ngerman for German
+%singlespacing, % Single line spacing, alternatives: onehalfspacing or doublespacing
+%draft, % Uncomment to enable draft mode (no pictures, no links, overfull hboxes indicated)
+%nolistspacing, % If the document is onehalfspacing or doublespacing, uncomment this to set spacing in lists to single
+%liststotoc, % Uncomment to add the list of figures/tables/etc to the table of contents
+%toctotoc, % Uncomment to add the main table of contents to the table of contents
+%parskip, % Uncomment to add space between paragraphs
+%nohyperref, % Uncomment to not load the hyperref package
+%headsepline, % Uncomment to get a line under the header
+%]{scrartcl or scrreprt or scrbook} % The class file specifying the document structure
+
+\usepackage{lmodern} % Diese beiden packages sorgen für echte
+\usepackage[T1]{fontenc} % Umlaute.
+
+\usepackage{amssymb, amsmath, color, graphicx, float, setspace, tipa}
+\usepackage[utf8]{inputenc}
+\usepackage[english]{babel}
+\usepackage[pdfpagelabels,
+ pdfstartview = FitH,
+ bookmarksopen = true,
+ bookmarksnumbered = true,
+ linkcolor = black,
+ plainpages = false,
+ hypertexnames = false,
+ citecolor = black,
+ breaklinks]{hyperref}
+\usepackage{url}
+\usepackage{array}
+
+\allowdisplaybreaks % allows page breaks in align/equation environment
+
+\usepackage{authblk} % titlepage stuff
+\usepackage[titletoc, title]{appendix}
+
+\usepackage{newclude} % use \include*{file} instead \include{} to omit pagebreak after include
+\usepackage{lscape}
+
+
+
+
+%===================
+% BIBLIOGRAPHY
+%===================
+
+\usepackage[]{natbib}
+\bibliographystyle{apalike}
+
+%DONT FORGET TO COMPILE THE BIBLIOGRAPHY WITH BIBTEX WHEN CHANGES ARE MADE.
+
+
+
+%--------------------------------------------
+% OPTIONAL
+%--------------------------------------------
+
+
+
+%% Change font
+%\newcommand{\changefont}[3]{
+%\fontfamily{#1} \fontseries{#2} \fontshape{#3} \selectfont}
+%\changefont{ppl}{m}{n} nach \begin{document} einsetzen
+
+% Fig. instead of Figure, Tab. instead of Table
+%\usepackage[footnotesize]{caption2}
+%\addto\captionsenglish{\renewcommand{\figurename}{Fig.}}
+%\addto\captionsngerman{\renewcommand{\figurename}{Fig.}}
+%\renewcommand{\tablename}{Tab.}%
+
+%\pagestyle{headings} % Write headings on each page
+
+%\usepackage{chngcntr} \counterwithout{figure}{section} % Integer only figure numbers, ignoring chapter numbers
+
+
+
+
+
+
+%-----------------------------------------------
+% FORMAT TITLE
+%-----------------------------------------------
+
+
+% Set fonts of document parts
+\setkomafont{title}{\rmfamily\bfseries\boldmath}
+\addtokomafont{section}{\rmfamily\bfseries\boldmath}
+\addtokomafont{subsection}{\rmfamily\bfseries\boldmath}
+\addtokomafont{subsubsection}{\rmfamily\bfseries\boldmath}
+\addtokomafont{disposition}{\rmfamily} % table of contents and stuff
+\setkomafont{descriptionlabel}{\rmfamily\bfseries\boldmath}
+
+
+\usepackage{dsfont} % for Pressure tensor P with \mathds{}
+
+
+
+
+%-----------------------------------------------
+% Document specific definitions
+%-----------------------------------------------
+
+\newcommand{\Aij}{$\mathcal{A}_{ij}$} % A_ij
+\newcommand{\Aijm}{\ensuremath{\mathcal{A}_{ij}}} % A_ij math
+\newcommand{\U}{\ensuremath{\mathcal{U}}} % State vector
+\newcommand{\W}{\ensuremath{\mathcal{W}}} % State vector
+\newcommand{\F}{\ensuremath{\mathcal{F}}} % Flux tensor
+\newcommand{\Ubf}{\ensuremath{\mathbf{U}}} % State vector
+\newcommand{\Fbf}{\ensuremath{\mathbf{F}}} % Flux tensor
+\newcommand{\psitilde}{\ensuremath{\tilde{\psi}}} % psi tilde
+\newcommand{\half}{1/2} % 1/2
+
+
+\newcommand{\absorbers}{\text{HI, HeI, HeII}}
+\newcommand{\Ndot}{\dot{N}}
+
+
+%-----------------------------------------------
+% General Lazyness
+%-----------------------------------------------
+
+\newcommand{\corresponds}{\mathrel{\widehat{=}}} % equals with hat
+
+\newcommand {\arctanh}{\mathrm{arctanh}} % Atanh
+\newcommand{\arccot}{\mathrm{arccot }} % Acotanh
+
+\newcommand{\limz}[1]{\lim\limits_{#1 \rightarrow 0}} % Limes of something towards zero
+
+\newcommand{\bm}{\boldmath} % Bold font in math
+\newcommand{\dps}{\displaystyle}
+
+\newcommand{\e}{\mbox{e}} % e noncursive in math mode
+
+\newcommand{\del}{\partial} % partial diff operator
+\newcommand{\de}{\mathrm{d}} % differential d
+\newcommand{\D}{\mathrm{d}} % differential d
+\newcommand{\GRAD}{\mathrm{grad}\ } % gradient
+\newcommand{\DIV}{\mathrm{div}\ } % divergence
+\newcommand{\ROT}{\mathrm{rot}\ } % rotation
+
+\newcommand{\CONST}{\mathrm{const.\ }} % constant
+\newcommand{\var}{\mathrm{var}} % variance
+
+\newcommand{\g}{^\circ} % degrees
+\newcommand{\degr}{^\circ} % degrees
+
+\newcommand{\msol}{M_\odot} % solar mass
+\newcommand{\Lsol}{L_\odot} % solar luminosity
+\newcommand{\Lsun}{L_\odot} % solar luminosity
+\newcommand{\order}{\mathcal{O}} % order, e.g. O(h^2)
+
+
+\newcommand{\x}{\mathbf{x}} % x vector
+\newcommand{\xdot}{\dot{\mathbf{x}}} % x dot vector
+\newcommand{\xddot}{\ddot{\mathbf{x}}} % x doubledot vector
+\newcommand{\R}{\mathbf{r}} % r vector
+\newcommand{\rdot}{\dot{\mathbf{r}}} % r dot vector
+\newcommand{\rddot}{\ddot{\mathbf{r}}} % r doubledot vector
+\newcommand{\vel}{\mathbf{v}} % v vector
+\newcommand{\V}{\mathbf{v}} % v vector
+\newcommand{\vdot}{\dot{\mathbf{v}}} % v dot vector
+\newcommand{\vddot}{\ddot{\mathbf{v}}} % v doubledot vector
+
+\newcommand{\dete}{\mathrm{d}t} % dt
+\newcommand{\delte}{\del t} % partial t
+\newcommand{\dex}{\mathrm{d}x} % dx
+\newcommand{\delx}{\del x} % partial x
+\newcommand{\der}{\mathrm{d}r} % dr
+\newcommand{\delr}{\del r} % partial r
+
+
+\newcommand{\deldx}{\frac{\del}{\del x}} % shortcut partial derivative, in line
+\newcommand{\ddx}{\frac{\de}{\de x}} % shortcut total derivative, in line
+\newcommand{\DELDX}[1]{\frac{\del #1}{\del x}} % shortcut partial derivative, on fraction
+\newcommand{\DDX}[1]{\frac{\de #1}{\de x}} % shortcut total derivative, on fraction
+
+\newcommand{\deldvecx}{\frac{\del}{\del \x}} % shortcut partial derivative, in line
+\newcommand{\ddvecx}{\frac{\de}{\de \x}} % shortcut total derivative, in line
+\newcommand{\DELDVECX}[1]{\frac{\del #1}{\del \x}} % shortcut partial derivative, on fraction
+\newcommand{\DDVECX}[1]{\frac{\de #1}{\de \x}} % shortcut total derivative, on fraction
+
+\newcommand{\deldr}{\frac{\del}{\del r}} % shortcut partial derivative, in line
+\newcommand{\ddr}{\frac{\de}{\de r}} % shortcut total derivative, in line
+\newcommand{\DELDR}[1]{\frac{\del #1}{\del r}} % shortcut partial derivative, on fraction
+\newcommand{\DDR}[1]{\frac{\de #1}{\de r}} % shortcut total derivative, on fraction
+
+\newcommand{\deldt}{\frac{\del}{\del t}} % shortcut partial derivative, in line
+\newcommand{\ddt}{\frac{\de}{\de t}} % shortcut total derivative, in line
+\newcommand{\DELDT}[1]{\frac{\del #1}{\del t}} % shortcut partial derivative, on fraction
+\newcommand{\DDT}[1]{\frac{\de #1}{\de t}} % shortcut total derivative, on fraction
+
+\newcommand{\deldxalpha}{\frac{\del}{\del x^\alpha}}
+\newcommand{\DELDXALPHA}[1]{\frac{\del #1}{\del x^\alpha}}
+\newcommand{\ddxalpha}{\frac{\de}{\de x^\alpha}}
+\newcommand{\DDXALPHA}[1]{\frac{\de #1}{\de x^\alpha}}
+
+
+
+
+
diff --git a/theory/RadiativeTransfer/GEARRT/my_defines.tex b/theory/RadiativeTransfer/GEARRT/my_defines.tex
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/theory/RadiativeTransfer/GEARRT/references.bib b/theory/RadiativeTransfer/GEARRT/references.bib
new file mode 100644
index 0000000000000000000000000000000000000000..957f446b944c2c6506e065aeeb60de4e615a9d1f
--- /dev/null
+++ b/theory/RadiativeTransfer/GEARRT/references.bib
@@ -0,0 +1,430 @@
+
+@article{hopkinsGIZMONewClass2015,
+ archivePrefix = {arXiv},
+ eprinttype = {arxiv},
+ eprint = {1409.7395},
+ title = {{{GIZMO}}: {{A New Class}} of {{Accurate}}, {{Mesh}}-{{Free Hydrodynamic Simulation Methods}}},
+ volume = {450},
+ issn = {0035-8711, 1365-2966},
+ shorttitle = {{{GIZMO}}},
+ abstract = {We present two new Lagrangian methods for hydrodynamics, in a systematic comparison with moving-mesh, SPH, and stationary (non-moving) grid methods. The new methods are designed to simultaneously capture advantages of both smoothed-particle hydrodynamics (SPH) and grid-based/adaptive mesh refinement (AMR) schemes. They are based on a kernel discretization of the volume coupled to a high-order matrix gradient estimator and a Riemann solver acting over the volume 'overlap.' We implement and test a parallel, second-order version of the method with self-gravity \& cosmological integration, in the code GIZMO: this maintains exact mass, energy and momentum conservation; exhibits superior angular momentum conservation compared to all other methods we study; does not require 'artificial diffusion' terms; and allows the fluid elements to move with the flow so resolution is automatically adaptive. We consider a large suite of test problems, and find that on all problems the new methods appear competitive with moving-mesh schemes, with some advantages (particularly in angular momentum conservation), at the cost of enhanced noise. The new methods have many advantages vs. SPH: proper convergence, good capturing of fluid-mixing instabilities, dramatically reduced 'particle noise' \& numerical viscosity, more accurate sub-sonic flow evolution, \& sharp shock-capturing. Advantages vs. non-moving meshes include: automatic adaptivity, dramatically reduced advection errors \& numerical overmixing, velocity-independent errors, accurate coupling to gravity, good angular momentum conservation and elimination of 'grid alignment' effects. We can, for example, follow hundreds of orbits of gaseous disks, while AMR and SPH methods break down in a few orbits. However, fixed meshes minimize 'grid noise.' These differences are important for a range of astrophysical problems.},
+ number = {1},
+ journal = {Monthly Notices of the Royal Astronomical Society},
+ doi = {10.1093/mnras/stv195},
+ author = {Hopkins, Philip F.},
+ month = jun,
+ year = {2015},
+ keywords = {Astrophysics - Cosmology and Nongalactic Astrophysics,Astrophysics - Instrumentation and Methods for Astrophysics,Physics - Computational Physics,Astrophysics - Astrophysics of Galaxies,Physics - Fluid Dynamics,_tablet},
+ pages = {53-110},
+ file = {/home/mivkov/Zotero/storage/23ALA3M9/Hopkins_2015_GIZMO.pdf;/home/mivkov/Zotero/storage/SHV6QDMB/1409.html}
+}
+
+@article{lansonRenormalizedMeshfreeSchemes2008,
+ title = {Renormalized {{Meshfree Schemes I}}: {{Consistency}}, {{Stability}}, and {{Hybrid Methods}} for {{Conservation Laws}}},
+ volume = {46},
+ issn = {0036-1429},
+ shorttitle = {Renormalized {{Meshfree Schemes I}}},
+ abstract = {This paper is devoted to the study of a new kind of meshfree scheme based on a new class of meshfree derivatives: the renormalized meshfree derivatives, which improve the consistency of the original weighted particle methods. The weak renormalized meshfree scheme, built from the weak formulation of general conservation laws, turns out to be \$L\^2\$ stable under some geometrical conditions on the distribution of particles and some regularity conditions of the transport field. A time discretization is then performed by analogy with finite volume methods, and the \$L\^1\$, \$L\^\textbackslash{}infty\$, and \$BV\$ stabilities of the obtained time discretized scheme are studied. From the same analogy with finite volume methods, a hybrid particle scheme is built using the Godunov method and is numerically compared to the weak renormalized scheme.},
+ number = {4},
+ journal = {SIAM J. Numer. Anal.},
+ doi = {10.1137/S0036142903427718},
+ author = {Lanson, Nathalie and Vila, Jean-Paul},
+ month = apr,
+ year = {2008},
+ keywords = {finite volume,hybrid particle schemes,meshfree methods,nonlinear conservation law,renormalization,SPH,weak scheme,_tablet},
+ pages = {1912--1934},
+ file = {/home/mivkov/Zotero/storage/QPBH9IJ6/Lanson_Vila_2008_Renormalized Meshfree Schemes I.pdf}
+}
+
+@article{ivanovaCommonEnvelopeEvolution2013,
+ archivePrefix = {arXiv},
+ eprinttype = {arxiv},
+ eprint = {1209.4302},
+ title = {Common {{Envelope Evolution}}: {{Where}} We Stand and How We Can Move Forward},
+ volume = {21},
+ issn = {0935-4956, 1432-0754},
+ shorttitle = {Common {{Envelope Evolution}}},
+ abstract = {This work aims to present our current best physical understanding of common-envelope evolution (CEE). We highlight areas of consensus and disagreement, and stress ideas which should point the way forward for progress in this important but long-standing and largely unconquered problem. Unusually for CEE-related work, we mostly try to avoid relying on results from population synthesis or observations, in order to avoid potentially being misled by previous misunderstandings. As far as possible we debate all the relevant issues starting from physics alone, all the way from the evolution of the binary system immediately before CEE begins to the processes which might occur just after the ejection of the envelope. In particular, we include extensive discussion about the energy sources and sinks operating in CEE, and hence examine the foundations of the standard energy formalism. Special attention is also given to comparing the results of hydrodynamic simulations from different groups and to discussing the potential effect of initial conditions on the differences in the outcomes. We compare current numerical techniques for the problem of CEE and also whether more appropriate tools could and should be produced (including new formulations of computational hydrodynamics, and attempts to include 3D processes within 1D codes). Finally we explore new ways to link CEE with observations. We compare previous simulations of CEE to the recent outburst from V1309 Sco, and discuss to what extent post-common-envelope binaries and nebulae can provide information, e.g. from binary eccentricities, which is not currently being fully exploited.},
+ number = {1},
+ journal = {The Astronomy and Astrophysics Review},
+ doi = {10.1007/s00159-013-0059-2},
+ author = {Ivanova, N. and Justham, S. and Chen, X. and De Marco, O. and Fryer, C. L. and Gaburov, E. and Ge, H. and Glebbeek, E. and Han, Z. and Li, X.-D. and Lu, G. and Marsh, T. and Podsiadlowski, Ph and Potter, A. and Soker, N. and Taam, R. and Tauris, T. M. and van den Heuvel, E. P. J. and Webbink, R. F.},
+ month = nov,
+ year = {2013},
+ keywords = {Astrophysics - High Energy Astrophysical Phenomena,Astrophysics - Solar and Stellar Astrophysics,_tablet},
+ file = {/home/mivkov/Zotero/storage/V6IL3NUY/Ivanova et al_2013_Common Envelope Evolution.pdf;/home/mivkov/Zotero/storage/MMFFSMV7/1209.html}
+}
+
+@phdthesis{vandenbrouckeAdvancedModelsSimulating,
+ title = {{Advanced models for simulating dwarf galaxy formation and evolution}},
+ language = {nl},
+ author = {Vandenbroucke, Bert},
+ year = {2016},
+ file = {/home/mivkov/Zotero/storage/T3UZJAWQ/Vandenbroucke - Advanced models for simulating dwarf galaxy format.pdf}
+}
+
+
+
+@article{ramses-rt13,
+ title = {{{RAMSES}}-{{RT}}: Radiation Hydrodynamics in the Cosmological Context},
+ shorttitle = {{{RAMSES}}-{{RT}}},
+ author = {Rosdahl, J. and Blaizot, J. and Aubert, D. and Stranex, T. and Teyssier, R.},
+ year = {2013},
+ month = dec,
+ volume = {436},
+ pages = {2188--2231},
+ issn = {0035-8711},
+ doi = {10.1093/mnras/stt1722},
+ abstract = {We present a new implementation of radiation hydrodynamics (RHD) in the adaptive mesh refinement (AMR) code RAMSES. The multigroup radiative transfer (RT) is performed on the AMR grid with a first-order Godunov method using the M1 closure for the Eddington tensor, and is coupled to the hydrodynamics via non-equilibrium thermochemistry of hydrogen and helium. This moment-based approach has the great advantage that the computational cost is independent of the number of radiative sources - it can even deal with continuous regions of emission such as bound-free emission from gas. As it is built directly into RAMSES, the RT takes natural advantage of the refinement and parallelization strategies already in place. Since we use an explicit advection solver for the radiative transport, the time-step is restricted by the speed of light - a severe limitation that can be alleviated using the so-called reduced speed of light approximation. We propose a rigorous framework to assess the validity of this approximation in various conditions encountered in cosmology and galaxy formation. We finally perform with our newly developed code a complete suite of RHD tests, comparing our results to other RHD codes. The tests demonstrate that our code performs very well and is ideally suited for exploring the effect of radiation on current scenarios of structure and galaxy formation.},
+ file = {/home/mivkov/Zotero/storage/M5QSAVLR/Rosdahl et al_2013_RAMSES-RT.pdf},
+ journal = {Monthly Notices of the Royal Astronomical Society},
+ keywords = {methods: numerical,radiative transfer}
+}
+
+@article{ramses-rt15,
+ title = {A Scheme for Radiation Pressure and Photon Diffusion with the {{M1}} Closure in {{RAMSES}}-{{RT}}},
+ author = {Rosdahl, J. and Teyssier, R.},
+ year = {2015},
+ month = jun,
+ volume = {449},
+ pages = {4380--4403},
+ issn = {0035-8711, 1365-2966},
+ doi = {10.1093/mnras/stv567},
+ abstract = {We describe and test an updated version of radiation-hydrodynamics (RHD) in the RAMSES code, that includes three new features: i) radiation pressure on gas, ii) accurate treatment of radiation diffusion in an unresolved optically thick medium, and iii) relativistic corrections that account for Doppler effects and work done by the radiation to first order in v/c. We validate the implementation in a series of tests, which include a morphological assessment of the M1 closure for the Eddington tensor in an astronomically relevant setting, dust absorption in a optically semi-thick medium, direct pressure on gas from ionising radiation, convergence of our radiation diffusion scheme towards resolved optical depths, correct diffusion of a radiation flash and a constant luminosity radiation, and finally, an experiment from Davis et al. of the competition between gravity and radiation pressure in a dusty atmosphere, and the formation of radiative Rayleigh-Taylor instabilities. With the new features, RAMSES-RT can be used for state-of-the-art simulations of radiation feedback from first principles, on galactic and cosmological scales, including not only direct radiation pressure from ionising photons, but also indirect pressure via dust from multi-scattered IR photons reprocessed from higher-energy radiation, both in the optically thin and thick limits.},
+ archivePrefix = {arXiv},
+ eprint = {1411.6440},
+ eprinttype = {arxiv},
+ file = {/home/mivkov/Zotero/storage/8WS8SH75/Rosdahl_Teyssier_2015_A scheme for radiation pressure and photon diffusion with the M1 closure in.pdf;/home/mivkov/Zotero/storage/U86AN9F6/1411.html},
+ journal = {Monthly Notices of the Royal Astronomical Society},
+ keywords = {_tablet_modified,Astrophysics - High Energy Astrophysical Phenomena,Astrophysics - Instrumentation and Methods for Astrophysics},
+ number = {4}
+}
+
+
+@article{ilievCosmologicalRadiativeTransfer2006,
+ title = {Cosmological {{Radiative Transfer Codes Comparison Project I}}: {{The Static Density Field Tests}}},
+ shorttitle = {Cosmological {{Radiative Transfer Codes Comparison Project I}}},
+ author = {Iliev, Ilian T. and Ciardi, Benedetta and Alvarez, Marcelo A. and Maselli, Antonella and Ferrara, Andrea and Gnedin, Nickolay Y. and Mellema, Garrelt and Nakamoto, Taishi and Norman, Michael L. and Razoumov, Alexei O. and Rijkhorst, Erik-Jan and Ritzerveld, Jelle and Shapiro, Paul R. and Susa, Hajime and Umemura, Masayuki and Whalen, Daniel J.},
+ year = {2006},
+ month = mar,
+ journal = {arXiv:astro-ph/0603199},
+ eprint = {astro-ph/0603199},
+ eprinttype = {arxiv},
+ doi = {10.1111/j.1365-2966.2006.10775.x},
+ abstract = {Radiative transfer simulations are now at the forefront of numerical astrophysics. They are becoming crucial for an increasing number of astrophysical and cosmological problems; at the same time their computational cost has come to the reach of currently available computational power. Further progress is retarded by the considerable number of different algorithms (including various flavours of ray-tracing and moment schemes) developed, which makes the selection of the most suitable technique for a given problem a non-trivial task. Assessing the validity ranges, accuracy and performances of these schemes is the main aim of this paper, for which we have compared 11 independent RT codes on 5 test problems: (0) basic physics, (1) isothermal H II region expansion and (2) H II region expansion with evolving temperature, (3) I-front trapping and shadowing by a dense clump, (4) multiple sources in a cosmological density field. The outputs of these tests have been compared and differences analyzed. The agreement between the various codes is satisfactory although not perfect. The main source of discrepancy appears to reside in the multi-frequency treatment approach, resulting in different thicknesses of the ionized-neutral transition regions and the temperature structure. The present results and tests represent the most complete benchmark available for the development of new codes and improvement of existing ones. To this aim all test inputs and outputs are made publicly available in digital form.},
+ archiveprefix = {arXiv},
+ langid = {english},
+ keywords = {Astrophysics},
+ file = {/home/mivkov/Zotero/storage/6L4AAYAN/Iliev et al. - 2006 - Cosmological Radiative Transfer Codes Comparison P.pdf}
+}
+
+
+
+@article{gonzalezHERACLESThreedimensionalRadiation2007,
+ title = {{{HERACLES}}: A Three-Dimensional Radiation Hydrodynamics Code},
+ shorttitle = {{{HERACLES}}},
+ author = {Gonz{\'a}lez, M. and Audit, E. and Huynh, P.},
+ year = {2007},
+ month = mar,
+ journal = {A\&A},
+ volume = {464},
+ number = {2},
+ pages = {429--435},
+ issn = {0004-6361, 1432-0746},
+ doi = {10.1051/0004-6361:20065486},
+ abstract = {Methods. The radiation transfer is modelled using a two-moment model and a closure
+relation that allows large angular anisotropies in the radiation field to be preserved and
+reproduced. The radiative equations thus obtained are solved by a second-order Godunov-type method
+and integrated implicitly by using iterative solvers. HERACLES has been parallelized with the MPI
+library and implemented in Cartesian, cylindrical, and spherical coordinates. To characterize the
+accuracy of HERACLES and to compare it with other codes, we performed a series of tests including
+purely radiative tests and radiation-hydrodynamics ones. Results. The results show that the physical
+model used in HERACLES for the transfer is fairly accurate in both the diffusion and transport
+limit, but also for semi-transparent regions. Conclusions. This makes HERACLES very well-suited to
+studying many astrophysical problems such as radiative shocks, molecular jets of young stars,
+fragmentation and formation of dense cores in the interstellar medium, and protoplanetary discs.},
+ langid = {english},
+ file = {/home/mivkov/Zotero/storage/NKTGBME2/González et al. - 2007 - HERACLES a three-dimensional
+radiation hydrodynam.pdf}
+}
+
+@article{ilievCosmologicalRadiativeTransfer2006,
+ title = {Cosmological {{Radiative Transfer Codes Comparison Project I}}: {{The Static Density
+Field Tests}}},
+ shorttitle = {Cosmological {{Radiative Transfer Codes Comparison Project I}}},
+ author = {Iliev, Ilian T. and Ciardi, Benedetta and Alvarez, Marcelo A. and Maselli, Antonella and
+Ferrara, Andrea and Gnedin, Nickolay Y. and Mellema, Garrelt and Nakamoto, Taishi and Norman,
+Michael L. and Razoumov, Alexei O. and Rijkhorst, Erik-Jan and Ritzerveld, Jelle and Shapiro, Paul
+R. and Susa, Hajime and Umemura, Masayuki and Whalen, Daniel J.},
+ year = {2006},
+ month = mar,
+ journal = {arXiv:astro-ph/0603199},
+ eprint = {astro-ph/0603199},
+ eprinttype = {arxiv},
+ doi = {10.1111/j.1365-2966.2006.10775.x},
+ abstract = {Radiative transfer simulations are now at the forefront of numerical astrophysics.
+They are becoming crucial for an increasing number of astrophysical and cosmological problems; at
+the same time their computational cost has come to the reach of currently available computational
+power. Further progress is retarded by the considerable number of different algorithms (including
+various flavours of ray-tracing and moment schemes) developed, which makes the selection of the most
+suitable technique for a given problem a non-trivial task. Assessing the validity ranges, accuracy
+and performances of these schemes is the main aim of this paper, for which we have compared 11
+independent RT codes on 5 test problems: (0) basic physics, (1) isothermal H II region expansion and
+(2) H II region expansion with evolving temperature, (3) I-front trapping and shadowing by a dense
+clump, (4) multiple sources in a cosmological density field. The outputs of these tests have been
+compared and differences analyzed. The agreement between the various codes is satisfactory although
+not perfect. The main source of discrepancy appears to reside in the multi-frequency treatment
+approach, resulting in different thicknesses of the ionized-neutral transition regions and the
+temperature structure. The present results and tests represent the most complete benchmark available
+for the development of new codes and improvement of existing ones. To this aim all test inputs and
+outputs are made publicly available in digital form.},
+ archiveprefix = {arXiv},
+ langid = {english},
+ keywords = {Astrophysics},
+ file = {/home/mivkov/Zotero/storage/6L4AAYAN/Iliev et al. - 2006 - Cosmological Radiative Transfer
+Codes Comparison P.pdf}
+}
+
+@article{ilievCosmologicalRadiativeTransfer2009,
+ title = {Cosmological {{Radiative Transfer Comparison Project II}}: {{The Radiation-Hydrodynamic
+Tests}}},
+ shorttitle = {Cosmological {{Radiative Transfer Comparison Project II}}},
+ author = {Iliev, Ilian T. and Whalen, Daniel and Mellema, Garrelt and Ahn, Kyungjin and Baek,
+Sunghye and Gnedin, Nickolay Y. and Kravtsov, Andrey V. and Norman, Michael and Raicevic, Milan and
+Reynolds, Daniel R. and Sato, Daisuke and Shapiro, Paul R. and Semelin, Benoit and Smidt, Joseph and
+Susa, Hajime and Theuns, Tom and Umemura, Masayuki},
+ year = {2009},
+ month = dec,
+ journal = {Monthly Notices of the Royal Astronomical Society},
+ volume = {400},
+ number = {3},
+ eprint = {0905.2920},
+ eprinttype = {arxiv},
+ pages = {1283--1316},
+ issn = {00358711, 13652966},
+ doi = {10.1111/j.1365-2966.2009.15558.x},
+ abstract = {The development of radiation hydrodynamical methods that are able to follow gas
+dynamics and radiative transfer self-consistently is key to the solution of many problems in
+numerical astrophysics. Such fluid flows are highly complex, rarely allowing even for approximate
+analytical solutions against which numerical codes can be tested. An alternative validation
+procedure is to compare different methods against each other on common problems, in order to assess
+the robustness of the results and establish a range of validity for the methods. Previously, we
+presented such a comparison for a set of pure radiative transfer tests (i.e. for fixed, non-evolving
+density fields). This is the second paper of the Cosmological Radiative Transfer (RT) Comparison
+Project, in which we compare 9 independent RT codes directly coupled to gasdynamics on 3 relatively
+simple astrophysical hydrodynamics problems: (5) the expansion of an H II region in a uniform
+medium; (6) an ionization front (I-front) in a 1/r2 density profile with a flat core, and (7), the
+photoevaporation of a uniform dense clump. Results show a broad agreement between the different
+methods and no big failures, indicating that the participating codes have reached a certain level of
+maturity and reliability. However, many details still do differ, and virtually every code has showed
+some shortcomings and has disagreed, in one respect or another, with the majority of the results.
+This underscores the fact that no method is universal and all require careful testing of the
+particular features which are most relevant to the specific problem at hand.},
+ archiveprefix = {arXiv},
+ langid = {english},
+ keywords = {Astrophysics - Cosmology and Nongalactic Astrophysics},
+ file = {/home/mivkov/Zotero/storage/2YNJDDST/Iliev et al. - 2009 - Cosmological Radiative Transfer
+Comparison Project.pdf}
+}
+
+@article{katzCosmologicalSimulationsTreeSPH1996,
+ title = {Cosmological {{Simulations}} with {{TreeSPH}}},
+ author = {Katz, Neal and Weinberg, David H. and Hernquist, Lars},
+ year = {1996},
+ month = jul,
+ journal = {The Astrophysical Journal Supplement Series},
+ volume = {105},
+ pages = {19},
+ issn = {0067-0049},
+ doi = {10.1086/192305},
+ abstract = {We describe numerical methods for incorporating gasdynamics into cosmological
+simulations and present illustrative applications to the cold dark matter (CDM) scenario. Our
+evolution code, a version of TreeSPH (Hernquist \& Katz 1989) generalized to handle comoving
+coordinates and periodic boundary conditions, combines smoothed-particle hydrodynamics (SPH) with
+the hierarchical tree method for computing gravitational forces. The Lagrangian hydrodynamics
+approach and individual time steps for gas particles give the algorithm a large dynamic range, which
+is essential for studies of galaxy formation in a cosmological context. The code incorporates
+radiative cooling for an optically thin, primordial composition gas in ionization equilibrium with a
+user-specified ultraviolet background. We adopt a phenomenological prescription for star formation
+that gradually turns cold, dense, Jeans-unstable gas into collisionless stars, returning supernova
+feedback energy to the surrounding medium. In CDM simulations, some of the baryons that fall into
+dark matter potential wells dissipate their acquired thermal energy and condense into clumps with
+roughly galactic masses. The resulting galaxy population is insensitive to assumptions about star
+formation; we obtain similar baryonic mass functions and galaxy correlation functions from
+simulations with star formation and from simulations without star formation in which we identify
+galaxies directly from the cold, dense gas.},
+ keywords = {_tablet,COSMOLOGY: DARK MATTER,COSMOLOGY: LARGE-SCALE STRUCTURE OF UNIVERSE,COSMOLOGY:
+THEORY,GALAXIES: FORMATION,HYDRODYNAMICS,METHODS: NUMERICAL},
+ file = {/home/mivkov/Zotero/storage/2YV7DC26/Katz et al_1996_Cosmological Simulations with
+TreeSPH.pdf}
+}
+
+@article{levermoreRelatingEddingtonFactors1984a,
+ title = {Relating {{Eddington}} Factors to Flux Limiters},
+ author = {Levermore, C. D.},
+ year = {1984},
+ month = feb,
+ journal = {Journal of Quantitative Spectroscopy and Radiative Transfer},
+ volume = {31},
+ number = {2},
+ pages = {149--160},
+ issn = {0022-4073},
+ doi = {10.1016/0022-4073(84)90112-2},
+ abstract = {Variable Eddington factors and flux-limiters have been introduced in the P-1 and
+diffusion equations, respectively, to handle situations when the underlying intensity is too
+anisotropic for the unmodified theories to remain valid. We present a derivation of a relation
+between the two for which a new approach to the diffusion approximation is used. Algebraic
+expressions for Eddington factors satisfying the moment conditions are not satisfactory for closing
+the P-1 equations but, by using the derived relation, yield acceptable flux-limited diffusion
+theories.},
+ langid = {english},
+ file = {/home/mivkov/Zotero/storage/WSX9USJ7/Levermore_1984_Relating Eddington factors to flux
+limiters.pdf;/home/mivkov/Zotero/storage/XAJD9TFY/0022407384901122.html}
+}
+
+@article{smithGrackleChemistryCooling2017a,
+ title = {Grackle: A {{Chemistry}} and {{Cooling Library}} for {{Astrophysics}}},
+ shorttitle = {Grackle},
+ author = {Smith, Britton D. and Bryan, Greg L. and Glover, Simon C. O. and Goldbaum, Nathan J. and
+Turk, Matthew J. and Regan, John and Wise, John H. and Schive, Hsi-Yu and Abel, Tom and Emerick,
+Andrew and O'Shea, Brian W. and Anninos, Peter and Hummels, Cameron B. and Khochfar, Sadegh},
+ year = {2017},
+ month = apr,
+ journal = {Mon. Not. R. Astron. Soc.},
+ volume = {466},
+ number = {2},
+ eprint = {1610.09591},
+ eprinttype = {arxiv},
+ pages = {2217--2234},
+ issn = {0035-8711, 1365-2966},
+ doi = {10.1093/mnras/stw3291},
+ abstract = {We present the Grackle chemistry and cooling library for astrophysical simulations and
+models. Grackle provides a treatment of non-equilibrium primordial chemistry and cooling for H, D,
+and He species, including H2 formation on dust grains; tabulated primordial and metal cooling;
+multiple UV background models; and support for radiation transfer and arbitrary heat sources. The
+library has an easily implementable interface for simulation codes written in C, C++, and Fortran as
+well as a Python interface with added convenience functions for semi-analytical models. As an
+open-source project, Grackle provides a community resource for accessing and disseminating
+astrochemical data and numerical methods. We present the full details of the core functionality, the
+simulation and Python interfaces, testing infrastructure, performance, and range of applicability.
+Grackle is a fully open-source project and new contributions are welcome.},
+ archiveprefix = {arXiv},
+ keywords = {Astrophysics - Astrophysics of Galaxies,Astrophysics - Cosmology and Nongalactic
+Astrophysics,Astrophysics - Instrumentation and Methods for Astrophysics},
+ file = {/home/mivkov/Zotero/storage/8YEYQAV2/Smith et
+al_2017_Grackle.pdf;/home/mivkov/Zotero/storage/QRXIKLGA/1610.html}
+}
+
+
+@article{ivkovicThesis,
+ author = {{Ivkovic}, Mladen},
+ title = "{GEAR-RT: Towards Exa-Scale Moment Based Radiative Transfer For Cosmological
+Simulations Using Task-Based Parallelism And Dynamic Sub-Cycling with SWIFT}",
+ journal = {arXiv e-prints},
+ keywords = {Astrophysics - Instrumentation and Methods for Astrophysics},
+ year = 2023,
+ month = feb,
+ eid = {arXiv:2302.12727},
+ pages = {arXiv:2302.12727},
+ doi = {10.48550/arXiv.2302.12727},
+archivePrefix = {arXiv},
+ eprint = {2302.12727},
+ primaryClass = {astro-ph.IM},
+ adsurl = {https://ui.adsabs.harvard.edu/abs/2023arXiv230212727I},
+ adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+
+@article{vernerAtomicDataAstrophysics1996,
+ title = {Atomic {{Data}} for {{Astrophysics}}. {{II}}. {{New Analytic Fits}} for {{Photoionization
+Cross Sections}} of {{Atoms}} and {{Ions}}},
+ author = {Verner, D. A. and Ferland, G. J. and Korista, K. T. and Yakovlev, D. G.},
+ year = {1996},
+ month = jul,
+ journal = {The Astrophysical Journal},
+ volume = {465},
+ eprint = {astro-ph/9601009},
+ pages = {487},
+ issn = {0004-637X, 1538-4357},
+ doi = {10.1086/177435},
+ urldate = {2022-11-11},
+ abstract = {We present a complete set of analytic fits to the non-relativistic photoionization
+cross sections for the ground states of atoms and ions of elements from H through Si, and S, Ar, Ca,
+and Fe. Near the ionization thresholds, the fits are based on the Opacity Project theoretical cross
+sections interpolated and smoothed over resonances. At higher energies, the fits reproduce
+calculated Hartree-Dirac-Slater photoionization cross sections.},
+ archiveprefix = {arxiv},
+ langid = {english},
+ keywords = {Astrophysics,Physics - Atomic Physics},
+ file = {/home/mivkov/Zotero/storage/GAQCFADD/Verner et al. - 1996 - Atomic Data for Astrophysics.
+II. New Analytic Fit.pdf}
+}
+
+
+
+@article{smithGRACKLEChemistryCooling2017,
+ ids = {smithGRACKLEChemistryCooling2017a,smithGrackleChemistryCooling2017a},
+ title = {{{GRACKLE}}: A Chemistry and Cooling Library for Astrophysics},
+ shorttitle = {{{GRACKLE}}},
+ author = {Smith, Britton D. and Bryan, Greg L. and Glover, Simon C. O. and Goldbaum, Nathan J. and
+Turk, Matthew J. and Regan, John and Wise, John H. and Schive, Hsi-Yu and Abel, Tom and Emerick,
+Andrew and O'Shea, Brian W. and Anninos, Peter and Hummels, Cameron B. and Khochfar, Sadegh},
+ year = {2017},
+ month = apr,
+ journal = {Monthly Notices of the Royal Astronomical Society},
+ volume = {466},
+ eprint = {1610.09591},
+ pages = {2217--2234},
+ issn = {0035-8711},
+ doi = {10.1093/mnras/stw3291},
+ urldate = {2021-10-27},
+ abstract = {We present the GRACKLE chemistry and cooling library for astrophysical simulations and
+models. GRACKLE provides a treatment of non-equilibrium primordial chemistry and cooling for H, D
+and He species, including H2 formation on dust grains; tabulated primordial and metal cooling;
+multiple ultraviolet background models; and support for radiation transfer and arbitrary heat
+sources. The library has an easily implementable interface for simulation codes written in C, C++
+and FORTRAN as well as a PYTHON interface with added convenience functions for semi-analytical
+models. As an open-source project, GRACKLE provides a community resource for accessing and
+disseminating astrochemical data and numerical methods. We present the full details of the core
+functionality, the simulation and PYTHON interfaces, testing infrastructure, performance and range
+of applicability. GRACKLE is a fully open-source project and new contributions are welcome.},
+ archiveprefix = {arxiv},
+ keywords = {astrochemistry,Astrophysics - Astrophysics of Galaxies,Astrophysics - Cosmology and
+Nongalactic Astrophysics,Astrophysics - Instrumentation and Methods for Astrophysics,galaxies:
+formation,methods: numerical},
+ annotation = {ADS Bibcode: 2017MNRAS.466.2217S},
+ file = {/home/mivkov/Zotero/storage/7JY7484L/Smith et
+al_2017_GRACKLE.pdf;/home/mivkov/Zotero/storage/5MWMNWRG/abstract.html;/home/mivkov/Zotero/storage/
+QRXIKLGA/1610.html}
+}
+
+
+@article{gnedinMultidimensionalCosmologicalRadiative2001,
+ title = {Multi-Dimensional Cosmological Radiative Transfer with a {{Variable Eddington
+Tensor}} Formalism},
+ author = {Gnedin, Nickolay Y. and Abel, Tom},
+ year = {2001},
+ month = oct,
+ journal = {New Astronomy},
+ volume = {6},
+ pages = {437--455},
+ issn = {1384-1076},
+ doi = {10.1016/S1384-1076(01)00068-9},
+ urldate = {2022-11-14},
+ abstract = {We present a new approach to numerically model continuum radiative transfer based on
+the Optically Thin Variable Eddington Tensor (OTVET) approximation. Our method insures the exact
+conservation of the photon number and flux (in the explicit formulation) and automatically switches
+from the optically thick to the optically thin regime. It scales as N log N with the number of
+hydrodynamic resolution elements and is independent of the number of sources of ionizing radiation
+(i.e. works equally fast for an arbitrary source function). We also describe an implementation of
+the algorithm in a Soften Lagrangian Hydrodynamic code (SLH) and a multi-frequency approach
+appropriate for hydrogen and helium continuum opacities. We present extensive tests of our method
+for single and multiple sources in homogeneous and inhomogeneous density distributions, as well as a
+realistic simulation of cosmological reionization.},
+ keywords = {Astrophysics},
+ annotation = {ADS Bibcode: 2001NewA....6..437G},
+ file = {/home/mivkov/Zotero/storage/FN64KVYB/Gnedin_Abel_2001_Multi-dimensional cosmological
+radiative transfer with a Variable Eddington.pdf}
+}
+
+
diff --git a/theory/RadiativeTransfer/GEARRT/run.sh b/theory/RadiativeTransfer/GEARRT/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..11be530f1e14df5686245a9fa5a07e8bf94d42eb
--- /dev/null
+++ b/theory/RadiativeTransfer/GEARRT/run.sh
@@ -0,0 +1,7 @@
+#!/bin/bash
+echo "Generating PDF..."
+
+pdflatex -jobname=GEARRT GEARRT.tex
+bibtex GEARRT
+pdflatex -jobname=GEARRT GEARRT.tex
+pdflatex -jobname=GEARRT GEARRT.tex
diff --git a/theory/RadiativeTransfer/GEARRT/tables/fitting_parameters.tex b/theory/RadiativeTransfer/GEARRT/tables/fitting_parameters.tex
new file mode 100644
index 0000000000000000000000000000000000000000..48aee571b4ba0c8bd3716444d86d35d281e2ee3e
--- /dev/null
+++ b/theory/RadiativeTransfer/GEARRT/tables/fitting_parameters.tex
@@ -0,0 +1,16 @@
+
+\begin{table}
+\centering
+\begin{tabular}{l | c c c c c c c }
+ & $E_0$ [eV] & $\sigma_0$ [cm$^2$] & $y_a$ & $P$ & $y_w$ & $y_0$ & $y_1$ \\
+\hline\\[-0.5em]
+H$^0$ & 0.4298 & 5.475$\times 10^{-14}$ & 32.88 & 2.963 & 0 & 0 & 0 \\
+He$^0$ & 13.61 & 9.492$\times 10^{-16}$ & 1.469 & 3.188 & 2.039 & 0.4434 & 2.136 \\
+He$^+$ & 1.720 & 1.369$\times 10^{-14}$ & 32.88 & 2.963 & 0 & 0 & 0 \\
+ \end{tabular}
+\caption{Fitting parameters used for the ionizing cross section parametrizations given in
+eq.~\ref{eq:sigma-parametrizaiton}. The parametrization and these parameters are taken from
+\cite{vernerAtomicDataAstrophysics1996}. Note that the table given in the appendix of
+\cite{ramses-rt13} contains a typo in $E_0$ for He$^0$, and is off by a factor of 100.}
+\label{tab:cross-sections}
+\end{table}
diff --git a/theory/RadiativeTransfer/GEARRT/tables/rt_variables.tex b/theory/RadiativeTransfer/GEARRT/tables/rt_variables.tex
new file mode 100644
index 0000000000000000000000000000000000000000..7284945af25f8ce4d406dd7126fd0e06652044ed
--- /dev/null
+++ b/theory/RadiativeTransfer/GEARRT/tables/rt_variables.tex
@@ -0,0 +1,125 @@
+%\begin{table}
+\begin{center}
+\begin{small}
+\begin{tabular}{p{0.32\textwidth} p{0.42\textwidth} p{0.26\textwidth}}
+
+
+variable & name & units (cgs) \\[.5em]
+%---------------------------------------------------------------------------------------------------
+\hline \\
+
+$I_\nu (\x, \mathbf{n}, t)$ &
+ specific intensity &
+ $[I_\nu] = \frac{\text{erg}}{\text{cm}^2 \text{ rad Hz s}}$
+\\[.5em]
+$u_\nu (\x, \mathbf{n}, t) = \frac{I_\nu}{c}$ &
+ radiation energy density &
+ $[u_\nu] = \frac{\text{erg}}{\text{cm}^3 \text{ rad Hz}}$
+\\[.5em]
+$E_\nu (\x, t) = \int_{4 \pi} \frac{I_\nu}{c} \de \Omega$ &
+ total energy density &
+ $[E_\nu] = \frac{\text{erg}}{\text{cm}^3 \text{ Hz}}$
+\\[.5em]
+$N_\nu (\x, t) = E_\nu / h$ &
+ photon number density &
+ $[N_\nu] = \frac{1}{\text{cm}^3 \text{ Hz}}$
+\\[.5em]
+$E_{rad} (\x, t) = \int_{0}^{\infty} E_\nu \de \nu$ &
+ total integrated energy density &
+ $[E_{rad}] = \frac{\text{erg}}{\text{cm}^3}$
+\\[.5em]
+$N_{rad} (\x, t) = E_{rad} / h$ &
+ integrated photon number density &
+ $[N_{rad}] = \frac{1}{\text{cm}^3}$
+\\[.5em]
+% $J_\nu (\x, t) = \int_{4 \pi} I_\nu \frac{\de \Omega}{4 \pi}
+% = \frac{c}{4 \pi} E_\nu$ &
+% mean radiation specific intensity &
+% $[J_\nu] = \frac{\text{erg}}{\text{cm}^2 \text{ Hz s}}$
+% \\[.5em]
+$\mathbf{F}_\nu(\x, t) = \int_{4 \pi} I_\nu \mathbf{n} \de \Omega$ &
+ radiation flux &
+ $[\mathbf{F}_\nu] = \frac{\text{erg}}{\text{cm}^2 \text{ s Hz}}$
+\\[.5em]
+$\mathbf{P}_\nu (\x, t) = \frac{\F_\nu}{c^2}$ &
+ radiation momentum density &
+ $[\mathbf{P}_\nu] = \frac{\text{erg}}{ \text{cm}^4 \text{ s}^{-1} \text{ Hz}}$
+\\[.5em]
+$\mathds{P}_\nu (\x, t) = \int_{4 \pi} \frac{I_\nu}{c} \mathbf{n} \otimes \mathbf{n} \de \Omega$ &
+ radiation pressure tensor &
+ $[\mathds{P}_\nu ] = \frac{\text{erg}}{\text{cm}^3 \text{ Hz}}$
+
+
+\\
+%---------------------------------------------------------------------------------------------------
+\hline\\
+
+
+$\eta_\nu(\x, t)$ &
+ source function (of $I_\nu$)&
+ $[ \eta_\nu ] = \frac{\text{erg}}{\text{cm}^3 \text{ rad Hz s}}$
+\\[.5em]
+$n(\x, t)$ &
+ number density &
+ $[ n ] = \text{cm}^{-3}$
+\\[.5em]
+$\mathcal{J}_\nu(T)$ &
+ (specific intensity of a) spectrum &
+ $[\mathcal{J}_\nu (T)] = \frac{\text{erg}}{\text{cm}^2 \text{ rad Hz s}}$
+\\[.5em]
+$\Gamma_{\nu, j}$ &
+ photoionization rate of species $j$ &
+ $[\Gamma] = \text{s}^{-1}$
+\\[.5em]
+$\mathcal{H}_{\nu,j}$ &
+ (photo)heating rate of species $j$ &
+ $[\mathcal{H}_{\nu,j}] = \frac{\text{erg}}{s\ cm}^{3}$
+\\[.5em]
+$\alpha_\nu = \frac{1}{\lambda_\nu} = \rho \kappa_\nu = n \sigma_\nu$ &
+ absorption coefficient &
+ $[\alpha_\nu] = \text{cm}^{-1}$
+\\[.5em]
+% $\alpha_E = \frac{\int_0^\infty \alpha_\nu E_\nu \de \nu}{\int_0^\infty E_\nu \de \nu}$ &
+% energy mean of abs. coeff. &
+% $[\alpha_E] = \text{cm}^{-1}$
+% \\[.5em]
+% $\alpha_F = \frac{\int_0^\infty \alpha_\nu F_\nu \de \nu}{\int_0^\infty F_\nu \de \nu}$ &
+% flux mean of abs. coeff. &
+% $[\alpha_F] = \text{cm}^{-1}$
+% \\[.5em]
+% $j_\nu $ &
+% emission coefficient &
+% $[j_\nu] = \frac{\text{erg}}{\text{cm}^3 \text{rad Hz s}}$
+% \\[.5em]
+$\lambda_\nu $ &
+ mean free path &
+ $[\lambda_\nu] = \text{cm}$
+\\[.5em]
+$\kappa_\nu $ &
+ opacity &
+ $[\kappa_\nu] = \frac{\text{cm}^2}{\text{g}}$
+\\[.5em]
+$\tau_\nu(s) = \int_0^s \alpha_\nu(x) \de x$ &
+ optical depth &
+ $[\tau_\nu] = 1$
+\\[.5em]
+$\sigma_{j\nu}$ &
+ interacton cross section (of species $j$)&
+ $[\sigma_{j\nu}] = \text{cm}^2$
+\\[.5em]
+$\sigma_{ij}^N = \int_{\nu_{i}} N_\nu \ \sigma_{j\nu} \de \nu / \int_{\nu_i} N_\nu \de \nu$ &
+ number weighted average cross section &
+ $[\sigma_{ij}^N] = \text{cm}^2$
+\\[.5em]
+$\sigma_{ij}^E = \int_{\nu i} E_\nu \ \sigma_{j\nu} \de \nu / \int_{\nu i} E_\nu \de \nu$ &
+ energy weighted average cross section &
+ $[\sigma_{ij}^E] = \text{cm}^2$
+%
+
+\end{tabular}
+\end{small}
+\end{center}
+%\caption{Common quantities appearing in the context of radiative transfer, and their units.}
+%\label{tab:rt-variables}
+%\end{table}
+
diff --git a/theory/RadiativeTransfer/GEARRT/tables/thermochemistry_rates.tex b/theory/RadiativeTransfer/GEARRT/tables/thermochemistry_rates.tex
new file mode 100644
index 0000000000000000000000000000000000000000..0b14e3a05e5e7f503ddc08fd72e78321fb4ae08f
--- /dev/null
+++ b/theory/RadiativeTransfer/GEARRT/tables/thermochemistry_rates.tex
@@ -0,0 +1,42 @@
+\begin{table}
+\begin{center}
+\begin{tabular}{lll}
+\hline\\
+$A_{H^+} $ &
+ $8.4 \times 10^{-11} T^{-1/2} T_3^{-0.2} (1 + T_6^{0.7})^{-1}$ &
+ RR for H$^{+}$
+\\[.5em]
+$A_{d} $ &
+ $1.5 \times 10^{-10} T^{-0.6353}$ &
+ dielectronic RR for He$^{+}$
+\\[.5em]
+$A_{He^+} $ &
+ $1.9 \times 10^{-3} T^{-1.5} \mathrm{e}^{-470000/T} (1 + 0.3 \mathrm{e}^{-94000/T})$ &
+ RR for He$^{+}$
+\\[.5em]
+$A_{He^{++}} $ &
+ $3.36 \times 10^{-10} T^{-1/2} T_3^{-0.2} (1 + T_6^{0.7})^{-1}$ &
+ RR for He$^{++}$
+\\[.5em]
+\hline\\
+$\Gamma_{H^{0}} $ &
+ $5.85 \times 10^{-11} T^{1/2} \mathrm{e}^{-157809.1/T} ( 1 + T_5^{1/2})^{-1}$&
+ CIR for H$^{0}$
+\\[.5em]
+$\Gamma_{He^{0}} $ &
+ $2.38 \times 10^{-11} T^{1/2} \mathrm{e}^{-285335.4.1/T} ( 1 + T_5^{1/2})^{-1}$&
+ CIR for He$^{0}$
+\\[.5em]
+$\Gamma_{He^{+}} $ &
+ $5.68 \times 10^{-12} T^{1/2} \mathrm{e}^{-631515/T} ( 1 + T_5^{1/2})^{-1}$&
+ CIR for He$^{+}$
+\\[.5em]
+\hline
+\end{tabular}
+\end{center}
+\caption{Temperature ($T$) dependent recombination rates (RR) and collisional ionization rates
+(CIR) for Hydrogen and Helium species, adapted from \citet{katzCosmologicalSimulationsTreeSPH1996}.
+All rates are in units of cm$^3$ s$^{-1}$. $T_n$ is shorthand for $T / 10^n$K.}
+\label{tab:coll-ion-rates-katz}
+\end{table}
+