numerical-radiation-raytracing.tex

\subsection{Radiation transport: ray tracing}
\label{sec.num.raytracing}

Stars and black holes strongly affect their surroundings through
radiation.  Radiation transport is a well-studied problem; however,
its treatment in multidimensional calculations is difficult because of
the dependence on seven variables -- three spatial, two angular,
frequency, and time.  The non-local nature of the thermal and
hydrodynamical response to radiation sources further adds to the
difficulty.  Here we briefly describe \enzo's ray tracing
implementation \moray, which is presented in full detail in
\citet{Wise11_Moray} with seven code tests and six applications.

We solve the radiative transfer equation in comoving coordinates
(given by Equation~\ref{eq:rteqn}).  We can make some appropriate
approximations to reduce the complexity of this equation in order to
include radiation transport in numerical calculations.  Typically
timesteps in dynamic calculations are small enough so that $\Delta a/a
\ll 1$, therefore $a_{em}/a \approx 1$ in any given timestep, reducing
the second term to $\hat{n} \partial I_\nu/\partial \mathbf{x}$.  To
determine the importance of the third term, we evaluate the ratio of
the third term to the second term.  This is $HL/c$, where $L$ is the
simulation box length.  If this ratio is $\ll 1$, we can ignore the
third term.  For example at $z=5$, this ratio is 0.1 when $L =
c/H(z=5)$ = 53 proper Mpc.  In large boxes where the light crossing
time is comparable to the Hubble time, then it becomes important to
consider cosmological redshifting and dilution of the radiation.  Thus
equation (\ref{eq:rteqn}) reduces to the non-cosmological form in this
local approximation,
%
\begin{equation}
  \frac{1}{c} \frac{\partial I_\nu}{\partial t} + 
  \hat{n} \frac{\partial I_\nu}{\partial \mathbf{x}} =
  -\kappa_\nu I_\nu + j_\nu .
\end{equation}
%
We choose to represent the source term $j_\nu$ as point sources of
radiation (e.g. stars, quasars) that emit radial rays that are
propagated along the direction $\hat{n}$.

Ray tracing is an accurate method to propagate radiation from point
sources on a computational grid as long as there are a sufficient
number of rays passing through each cell.  Along a ray, the radiation
transfer equation reduces to
%
\begin{equation}
\label{eqn:rtray}
\frac{1}{c} \frac{\partial P}{\partial t} + \frac{\partial P}{\partial
  r} = -\kappa P,
\end{equation}
where $P$ is the photon number flux along the ray.  To sample the
radiation field at large radii, ray tracing requires at least $N_{ray}
= 4\pi R^2 / (\Delta x)^2$ rays to sample each cell with one ray,
where $R$ is the radius from the source to the cell and $\Delta x$ is
the cell width.  If one were to trace $N_{ray}$ rays out to $R$, the
cells at a smaller radius $r$ would be sampled by, on average,
$(r/R)^2$ rays, which is computationally wasteful because only a few
rays per cell are required to provide an accurate calculation of the
radiation field (as we will show later).

We avoid this inefficiency by utilizing adaptive ray tracing
\citep{Abel02_RT}, which is based on Hierarchical Equal Area
isoLatitude Pixelation \citep[HEALPix;][]{HEALPix} and progressively
splits rays when the sampling becomes too coarse.  In this approach,
the rays are traced along normal directions of the centers of the
HEALPix pixels that evenly divide a sphere into equal areas.  The rays
are initialized at each point source with the photon luminosity
(photon s$^{-1}$) equally spread across $N_{\rm pix} = 12 \times 4^l$
rays, where $l$ is the initial HEALPix level.  We usually find $l$ = 0
or 1 is sufficient because these coarse rays will usually be split
before traversing the first cell.

The rays are traced through the grid in a typical fashion
\citep[e.g.][]{Abel99_RT}, in which we calculate the next cell
boundary crossing.  The ray segment length crossing the cell is
%
\begin{equation}
  \label{eqn:trace}
  dr = R_0 - \min_{i=1 \rightarrow 3} \left[(x_{\rm cell,i} - x_{\rm src,i}) /
    \hat{n}_{\rm ray,i} \right],
\end{equation}
%
where $R_0$, $\hat{n}_{\rm ray}$, $x_{\rm cell,i}$, and $x_{\rm
  src,i}$ are the initial distance traveled by the ray, normal
direction of the ray, the next cell boundary crossing in the $i$-th
dimension, and the position of the point source that emitted the ray,
respectively.  However before the ray travels across the cell, we
evaluate the ratio of the face area $A_{\rm cell}$ of the current cell
and the solid angle $\Omega_{\rm ray}$ of the ray,
%
\begin{equation}
  \label{eqn:split}
  \Phi_c = \frac{A_{\rm cell}} {\Omega_{\rm ray}} = 
  \frac{N_{\rm pix} (\Delta x)^2} {4\pi R_0^2}.
\end{equation}
%
If $\Phi_c$ is less than a pre-determined value (usually $>3$), the
ray is split into 4 child rays.  The pixel numbers of the child rays
$p^\prime$ are given by the ``nested'' scheme of HEALPix at the next
level, i.e. $p^\prime = 4 \times p + [0,1,2,3]$, where $p$ is the
original pixel number.  The child rays (1) acquire the new normal
vectors of the pixels, (2) retain the same radius of the parent ray,
and (3) get a quarter of the photon flux of the parent ray.
Afterward, the parent ray is discontinued.

A ray propagates and splits until at least one of the following
conditions is met: (1) the photon has traveled $c \times dt_P$, where
$dt_P$ is the radiative transfer timestep; (2) its photon flux is
almost fully absorbed ($>99.9\%$) in a single cell, which
significantly reduces the computational time if the radiation volume
filling fraction is small; (3) the photon leaves the computational
domain with isolated boundary conditions; or (4) the photon travels
$\sqrt{3}$ of the simulation box length with periodic boundary
conditions.  In the first case, the photon is halted at that position
and saved, where it will be considered in the solution of $I_\nu$ at
the next timestep.  In the next timestep, the photon will encounter a
different hydrodynamical and ionization state, hence $\kappa$, in its
path.  Furthermore any time variations of the luminosities will be
retained in the radiation field.  This is how this method retains the
time derivative of the radiative transfer equation.  The last
restriction prevents our method from considering sources external to
the computational domain. However, a uniform radiation background can
be used in conjunction with ray tracing that adds the background
intensity to the local radiation field.

The radiation field is calculated by integrating Equation
(\ref{eqn:rtray}) along each ray, which is done by considering the
discretization of the ray into segments.  In the following
description, we assume the rays are monochromatic for simplicity.  For
convenience, we express the integration in terms of optical depth
$\tau = \int \kappa(r,t) \; dr$, and for a ray segment
%
\begin{equation}
  \label{eqn:dtau}
  d\tau = \sigma_{\rm abs}(\nu) n_{\rm abs} dr.
\end{equation}
Here $\sigma_{\rm abs}$ and $n_{\rm abs}$ are the cross section and
number density of the absorbing medium, respectively.  In the static
case, equation (\ref{eqn:rtray}) has a simple exponential analytic
solution, and the photon flux of a ray is reduced by
%
\begin{equation}
  \label{eqn:flux}
  dP = P \times (1 - e^{-\tau})
\end{equation}
as it crosses a cell.  We equate the photo-ionization rate to the
absorption rate, resulting in photon conservation \citep{Abel99_RT,
Mellema06}.  Thus the photo-ionization and photo-heating rates
associated with a single ray ($k_{\rm ph}$ and $\Gamma_{\rm ph}$,
respectively) are
%
\begin{equation}
  \label{eqn:kph}
  k_{\rm ph} = \frac{P (1 - e^{-\tau})}{n_{\rm abs} \; V_{\rm cell} \; dt_P},
\end{equation}
\begin{equation}
  \label{eqn:gamma}
  \Gamma_{\rm ph} = k_{\rm ph} \; (E_{\rm ph} - E_i),
\end{equation}
where $V_{\rm cell}$ is the cell volume, $E_{\rm ph}$ is the photon
energy, and $E_i$ is the ionization energy of the absorbing material.
In each cell, the photo-ionization and photo-heating rates from each
ray in the calculation are summed. After the ray tracing is complete,
these rates are used as inputs to the solver described in
Section~\ref{sec.ov.chem} to update the ionization, chemical, and
energy states of the gas in each cell.