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 (obsolete)

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.

\(\color{black}\)

Local ionization equilibrium solver (TODO)

Motivation

The ionization module Kuiper etal 2020 evolves the ionization fraction \(x\) and the neutral fraction \(y\) of hydrogen. In the current implementation, the local ionization and recombination rate equations are advanced in time using implicit Euler steps. Since the ionization and recombination timescales are assumed to be much shorter than the MHD timescale, this local time evolution is expected to relax the gas toward a quasi-stationary ionization state between two MHD updates.

The aim of this test is to investigate whether this quasi-stationary state can be found more directly. Instead of approaching the local equilibrium through a sequence of implicit time steps, the stationary rate equation is solved for its root in each cell.

This does not change the underlying local rate equation. The difference is only numerical: the current approach obtains the equilibrium state by integrating the rate equation in pseudo-time, while the tested approach tries to solve the stationary balance condition directly. During this local ionization update, the hydrodynamic transport terms are neglected and the gas density is treated as fixed.

Theory

Starting point

The rate equation for the ionization fraction \(x\) is

\begin{align} \partial_t(\rho_{\mathrm{gas}} x) + \vec{\nabla} \cdot (\rho_{\mathrm{gas}} x \vec{u}_{\mathrm{gas}}) = \rho_{\mathrm{gas}} \left[ (\sigma_{\mathrm{EUV}}u_{\mathrm{EUV}} + \sigma_{\mathrm{rec}}u_{\mathrm{rec}})c \, y + C(T_{\mathrm{gas}}) n_{\mathrm{H}} \, x y - \alpha^{(1)}(T_{\mathrm{gas}}) n_{\mathrm{H}} \, x^2 \right]. \end{align}

For the local ionization update, the advection term is neglected and the gas density is treated as constant:

\begin{align} \vec{\nabla} \cdot(\rho_{\mathrm{gas}} x \vec{u}_{\mathrm{gas}}) = 0, \qquad \partial_t(\rho_{\mathrm{gas}}x) = \rho_{\mathrm{gas}} \partial_t x . \end{align}

After canceling \(\rho_{\mathrm{gas}}\), this gives

\begin{align} \partial_t x = (\sigma_{\mathrm{EUV}}u_{\mathrm{EUV}} + \sigma_{\mathrm{rec}}u_{\mathrm{rec}})c \, y + C(T_{\mathrm{gas}}) n_{\mathrm{H}} \, x y - \alpha^{(1)}(T_{\mathrm{gas}}) n_{\mathrm{H}} \, x^2 . \end{align}

The terms correspond to direct and diffuse recombination photoionization, collisional ionization, and recombination.

Cell-averaged direct EUV photon density

Instead of using only the incoming photon density at the cell interface, a cell-averaged value is used for the direct EUV photon density:

\begin{align} \langle u_{\mathrm{EUV}}\rangle = u_L \frac{1-e^{-\Delta\tau}}{\Delta\tau}, \qquad \Delta\tau = n_{\mathrm{H}} \sigma_{\mathrm{EUV}} \Delta x \, y . \end{align}

Here, \(u_L\) is the incoming direct EUV photon number density at the left or inner cell interface. Inserting this into the direct photoionization term gives

\begin{align} \partial_t x = \sigma_{\mathrm{EUV}} c \, y \, u_L \frac{1-e^{-n_{\mathrm{H}}\sigma_{\mathrm{EUV}}\Delta x \, y}} {n_{\mathrm{H}}\sigma_{\mathrm{EUV}}\Delta x \, y} + \sigma_{\mathrm{rec}}u_{\mathrm{rec}}c \, y + C(T_{\mathrm{gas}}) n_{\mathrm{H}} \, x y - \alpha^{(1)}(T_{\mathrm{gas}}) n_{\mathrm{H}} \, x^2 . \end{align}

The factor \(y\sigma_{\mathrm{EUV}}\) cancels in the first term, resulting in

\begin{align} \partial_t x = \frac{c u_L}{n_{\mathrm{H}}\Delta x} \left( 1-e^{-n_{\mathrm{H}}\sigma_{\mathrm{EUV}}\Delta x \, y} \right) + \sigma_{\mathrm{rec}}u_{\mathrm{rec}}c \, y + C(T_{\mathrm{gas}}) n_{\mathrm{H}} \, x y - \alpha^{(1)}(T_{\mathrm{gas}}) n_{\mathrm{H}} \, x^2 . \end{align}

With \(y=1-x\), the local rate equation becomes

\begin{align} \partial_t x = \frac{c u_L}{n_{\mathrm{H}}\Delta x} \left[ 1-e^{-n_{\mathrm{H}}\sigma_{\mathrm{EUV}}\Delta x \, (1-x)} \right] + \sigma_{\mathrm{rec}}u_{\mathrm{rec}}c \, (1-x) + C(T_{\mathrm{gas}}) n_{\mathrm{H}} \, x(1-x) - \alpha^{(1)}(T_{\mathrm{gas}}) n_{\mathrm{H}} \, x^2 . \end{align}

Compact form

The following abbreviations are introduced:

\begin{align} A_{\mathrm{E}} = \frac{c u_L}{n_{\mathrm{H}}\Delta x}, \qquad q_{\mathrm{E}} = n_{\mathrm{H}}\sigma_{\mathrm{EUV}}\Delta x, \qquad A_{\mathrm{rec}} = \sigma_{\mathrm{rec}}u_{\mathrm{rec}}c, \qquad A_{\mathrm{C}} = C(T_{\mathrm{gas}}) n_{\mathrm{H}}, \qquad A_{\mathrm{R}} = \alpha^{(1)}(T_{\mathrm{gas}}) n_{\mathrm{H}} . \end{align}

This gives

\begin{align} \partial_t x = A_E \left[ 1-e^{-q_E(1-x)} \right] + A_{\mathrm{rec}} \, (1-x) + A_{\mathrm{C}} \, x(1-x) - A_{\mathrm{R}} \, x^2 . \end{align}

Stationary local equilibrium

The stationary local equilibrium corresponds to the ionization state for which

\begin{align} \partial_t x = 0 . \end{align}

Therefore, the stationary ionization fraction is obtained as a root of

\begin{align} f(x) = A_E \left[ 1-e^{-q_E(1-x)} \right] + A_{\mathrm{rec}} \, (1-x) + A_{\mathrm{C}} \, x(1-x) - A_{\mathrm{R}} \, x^2 . \end{align}

The solution must lie in the physical interval

\begin{align} 0 \le x \le 1 . \end{align}