MA Fynn Wawrzyniak

Initial Conditions

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

Problems

Mean Photon Density in Irradiation Update (belt/src/Makemake2.1/Irradiation.c ~ Line 570)

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. The 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);

Temperature-Based Initialization of Missing Ionization Fields

For restarts without stored ionization fields, the ionized and neutral fractions have to be initialized manually. If all cells are initialized as fully neutral, hot gas close to the inner boundary can become artificially optically thick to EUV radiation, because the ionizing optical depth depends directly on the neutral fraction,

\begin{align} \tau_\mathrm{ion} = n_\mathrm{H} y \sigma_\star \Delta x . \end{align}

To avoid this, the missing ionization fields are initialized once at startup using the gas temperature. This is only an initial condition for restarts without ionization data. A smooth temperature transition is used between \(T_0 = 1\cdot10^4 \,\mathrm{K}\) and \(T_1 = 2\cdot10^4 \,\mathrm{K}\):

\begin{align} s = \frac{T_\mathrm{gas} - T_0}{T_1 - T_0}, \qquad 0 \leq s \leq 1 , \end{align}

with \(s\) clipped to the interval \([0,1]\). The initial ionized fraction is then

\begin{align} x_\mathrm{init} = s^2(3-2s), \end{align}

and the neutral fraction is set consistently as

\begin{align} y_\mathrm{init} = 1 - x_\mathrm{init}. \end{align}

This prevents already hot gas from being initialized as fully neutral and therefore artificially opaque. After the initialization, \(x\) and \(y\) evolve normally through the ionization-recombination solver.