MA Fynn Wawrzyniak

Initial Conditions

\( \color{red} TODO \color{black} \)

Problems

Mean Photon Density in Irradiation Update (

During the iterative irradiation update, the local ionization update requires a representative photon density inside each cell for the photoionization term. The ray tracer knows the incoming photon density at the left cell face, \(u_\mathrm{L}\), and computes the outgoing photon density at the right cell face, \(u_\mathrm{R}\).

In the following, only the absorption inside one cell is considered. Attenuation factor in the code describes the geometrical dilution of the ray between the inner and outer cell face. Therefore, the derivation below focuses only on the local exponential absorption.

Assuming a constant absorption coefficient inside the cell, the photon density decreases exponentially with distance \(s\):

\begin{align} u(s) = u_\mathrm{L} e^{-\chi s} \qquad , \qquad 0 \leq s \leq \Delta x . \end{align}

Here \(\chi\) is the local ionizing absorption coefficient. For ionizing stellar photons,

\begin{align} \chi = n_\mathrm{H} \, y \, \sigma_\star , \end{align}

where \(n_\mathrm{H}\) is the hydrogen number density, \(y\) is the neutral fraction, and \(\sigma_\star\) is the ionization cross section. The optical depth across the full cell is

\begin{align} \tau = \chi \Delta x . \end{align}

The outgoing photon density at the right cell face is therefore

\begin{align} u_\mathrm{R} = u_\mathrm{L} e^{-\tau}. \end{align}

However, the local rate equations should not use only the right-face value. The photoionization term requires a representative photon density inside the cell. The physically motivated choice is the cell average,

\begin{align} \langle u \rangle = \frac{1}{\Delta x} \int_0^{\Delta x} u(s)\,\mathrm{d}s . \end{align}

Inserting the exponential profile gives

\begin{align} \langle u \rangle &= \frac{1}{\Delta x} \int_0^{\Delta x} u_\mathrm{L} e^{-\chi s} \,\mathrm{d}s \\ &= \frac{u_\mathrm{L}}{\Delta x} \frac{1 - e^{-\chi \Delta x}}{\chi} \\ &= u_\mathrm{L} \frac{1 - e^{-\tau}}{\tau}. \end{align}

Thus, the photon density used in the local ionization update should be

\begin{align} \boxed{ \langle u \rangle = u_\mathrm{L} \frac{1 - e^{-\tau}}{\tau} } \end{align}

instead of a simple arithmetic average between the left and right cell faces.

Previously, the mean photon density was approximated as

\begin{align} \langle u \rangle_\mathrm{old} &= \frac{1}{2} \left( u_\mathrm{L} + u_\mathrm{R} \right) \\ &= u_\mathrm{L} \frac{1 + e^{-\tau}}{2}. \end{align}

This approximation is reasonable only for optically thin cells. In this limit, both expressions agree to first order in \(\tau\). The cell-averaged expression gives

\begin{align} \frac{1 - e^{-\tau}}{\tau} = 1 - \frac{\tau}{2} + \mathcal{O}(\tau^2), \end{align}

while the old arithmetic average gives

\begin{align} \frac{1 + e^{-\tau}}{2} = 1 - \frac{\tau}{2} + \mathcal{O}(\tau^2). \end{align}

Therefore, for \(\tau \ll 1\), both estimates are very similar.

For optically thick cells, \(\tau \gg 1\), the outgoing photon density becomes nearly zero,

\begin{align} u_\mathrm{R} = u_\mathrm{L} e^{-\tau} \rightarrow 0 . \end{align}

The old approximation then approaches

\begin{align} \langle u \rangle_\mathrm{old} \rightarrow \frac{1}{2} u_\mathrm{L}. \end{align}

This overestimates the mean photon density in highly optically thick cells, because it assumes that roughly half of the incoming photon density is representative for the whole cell.

The cell-averaged expression has the correct optically thick limiting behavior. For

\begin{align} \tau \gg 1: \qquad \frac{1 - e^{-\tau}}{\tau} \approx \frac{1}{\tau}, \end{align}

one obtains

\begin{align} \langle u \rangle \approx \frac{u_\mathrm{L}}{\tau}. \end{align}

Therefore, the right cell face can be almost completely dark while the cell-averaged photon density remains finite. However, it decreases with increasing optical depth as \(1/\tau\), instead of approaching the constant value \(u_\mathrm{L}/2\) from the old arithmetic average.

In the crashing cell, the optical depth was extremely large,

\begin{align} \tau \simeq 3.7\times 10^4 . \end{align}

In this limit, the outgoing photon density is numerically zero,

\begin{align} u_\mathrm{R} = u_\mathrm{L} e^{-\tau} \approx 0 . \end{align}

The old arithmetic estimate therefore gives

\begin{align} \langle u\rangle_\mathrm{old} \approx \frac{1}{2}u_\mathrm{L}, \end{align}

whereas the cell-averaged expression gives

\begin{align} \langle u\rangle \approx \frac{u_\mathrm{L}}{\tau}. \end{align}

Thus, for the crashing cell, the old estimate overestimated the photon density entering the local photoionization update by approximately

\begin{align} \frac{\langle u\rangle_\mathrm{old}}{\langle u\rangle} \approx \frac{\tau}{2} \sim 2\times 10^4 . \end{align}

This was likely problematic for the ionization iteration: the ray tracer treated the cell as optically thick enough to absorb essentially all ionizing photons, while the local rate equations still used a comparatively large photon density inside the same cell.

The corresponding code change is

dtau_ion   = nH * y * sigmastar * dx;
Extinction = exp(-dtau_ion);
ustarxr    = ustarxl * Attenuation * Extinction;

/*
 * Old approximation:
 *     ustar = 0.5 * (ustarxl + ustarxr);
 */
if (dtau_ion > 1e-12) {
    ustar = ustarxl * (1.0 - Extinction) / dtau_ion;
} else {
    /*
     * Optically thin limit:
     *     (1 - exp(-dtau_ion)) / dtau_ion -> 1
     *
     * Avoid numerical division by zero.
     */
    ustar = ustarxl;
}

x = IonizationFractionRateEquation(ustar, urec, Tgas, nH, x);
y = NeutralFractionRateEquation(   ustar, urec, Tgas, nH, y);

This change makes the photon density used in the local update consistent with the optical depth used by the ray tracer. In optically thin cells, the result is essentially unchanged. In optically thick cells, however, the representative photon density is reduced from approximately \(u_\mathrm{L}/2\) to approximately \(u_\mathrm{L}/\tau\), which is physically more consistent with the exponential attenuation inside the cell.