Skip to content
Snippets Groups Projects
Commit cb625845 authored by Mladen Ivkovic's avatar Mladen Ivkovic Committed by Matthieu Schaller
Browse files

GEAR-RT: theory

parent 6b6eea0d
No related branches found
No related tags found
2 merge requests!1766added read Vz factor from the yml files.,!1744GEAR-RT: theory
Showing
with 2366 additions and 0 deletions
......@@ -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
......
%===================================================================================
\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.
%=============================================================================
\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*}
This diff is collapsed.
\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}
%------------------------------------------------
\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}.
\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}}
This diff is collapsed.
#!/bin/bash
echo "Generating PDF..."
pdflatex -jobname=GEARRT GEARRT.tex
bibtex GEARRT
pdflatex -jobname=GEARRT GEARRT.tex
pdflatex -jobname=GEARRT GEARRT.tex
\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}
%\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}
\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}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment