MA Fynn Wawrzyniak

Initial Conditions

The simulations are performed on a spherical grid in the coordinates \((r,\theta,\phi)\). However, for the description of the disk structure it is useful to additionally use the corresponding cylindrical coordinates \((R,z,\phi)\).

Disk

The disk is initialized as a smooth and axisymmetric model. The initial conditions do not represent a fully self-consistent equilibrium. Instead, they serve as a physically reasonable and numerically stable starting configuration. For this purpose, several approximations of a thin disk are used in the following.

Deviations from a complete equilibrium are adjusted during the first time steps. The aim is therefore to keep the transition from the prescribed initial condition to the physically evolved state as short and stable as possible.

Temperature distribution

Inside the disk, it is assumed that the initial temperature profile depends only on the cylindrical radius, \(T_{\mathrm{disk}} = T_{\mathrm{disk}}(R)\). This assumption is consistent with the thin-disk approximation. Close to the disk midplane, a geometrically thin disk satisfies \(|z| \ll R\). Therefore, from \(r = \sqrt{R^2+z^2}\), one obtains approximately \(r \approx R\). The temperature profile inside the disk can thus be described as a function of the cylindrical radius.

For the initial temperature profile, a simple power law is used.

\begin{align} T_{\mathrm{disk}}(R) = T_0 \left(\frac{R}{R_0}\right)^{-p_T} \end{align}

Here, \(T_0\) is the temperature at the reference radius \(R_0\). For the exponent, \(p_T = \frac{1}{2}\) is chosen. This choice is motivated by the typical temperature decrease of a disk irradiated from the inside. For a central radiation source, the radiative flux \(F\) decreases approximately as \(F \,\propto\, R^{-2}\) for geometrical reasons. Combining this with the rough estimate \(F \,\propto\, T^4\) gives \(T \,\propto\, R^{-1/2}\).

This estimate is strongly simplified, since the actual temperature structure additionally depends on the optical depth, the vertical structure of the disk, the irradiation geometry, and the local radiation transport. Nevertheless, it provides a useful qualitative approximation for a smooth initial temperature profile. In addition, the mean temperature profile inside the disk was analyzed from already evolved simulations. The comparison in Fig.~\ref{fig:initial_conditions_tgas_profile} shows that a power law with \(p_T = \frac{1}{2}\) approximately describes the radial temperature profile in the disk. Based on these simulation results, the normalization is chosen as \(T_0 = 1000\,\mathrm{K}\).

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

\( \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 photoionization by stellar EUV photons, diffuse photoionization by recombination photons, 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} K_1 = \frac{c u_L}{n_{\mathrm{H}}\Delta x}, \qquad Q_1 = n_{\mathrm{H}}\sigma_{\mathrm{EUV}}\Delta x, \qquad K_2 = \sigma_{\mathrm{rec}}u_{\mathrm{rec}}c, \qquad K_3 = C(T_{\mathrm{gas}}) n_{\mathrm{H}}, \qquad K_4 = \alpha^{(1)}(T_{\mathrm{gas}}) n_{\mathrm{H}} . \end{align}

This gives

\begin{align} \partial_t x = K_1 \left[ 1-e^{-Q_1(1-x)} \right] + K_2 \, (1-x) + K_3 \, x(1-x) - K_4 \, 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) = K_1 \left[ 1-e^{-Q_1(1-x)} \right] + K_2 \, (1-x) + K_3 \, x(1-x) - K_4 \, x^2 . \end{align}

The solution must lie in the physical interval

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

Stability of stationary solutions

A stationary solution of the local rate equation satisfies

\begin{align} f(x_\ast)=0 . \end{align}

However, not every mathematical root is necessarily the relevant equilibrium state. To check whether a stationary solution is stable, one can perturb it slightly,

\begin{align} x(t) = x_\ast + \delta x(t), \qquad |\delta x| \ll 1 . \end{align}

Inserting this into the local rate equation gives

\begin{align} \partial_t x = \partial_t (x_\ast+\delta x) = \partial_t \delta x = f(x_\ast+\delta x). \end{align}

The right-hand side can be expanded around \(x_\ast\):

\begin{align} f(x_\ast+\delta x) = f(x_\ast) + f'(x_\ast)\,\delta x + \mathcal{O}(\delta x^2). \end{align}

Since \(x_\ast\) is a stationary solution, \(f(x_\ast)=0\). Keeping only terms linear in the perturbation gives

\begin{align} \partial_t \delta x = f'(x_\ast)\,\delta x . \end{align}

Thus, a small perturbation evolves as

\begin{align} \delta x(t) = \delta x(0)\exp\left[f'(x_\ast)t\right]. \end{align}

Therefore, the stationary point is locally stable if

\begin{align} f'(x_\ast) < 0 , \end{align}

and unstable if

\begin{align} f'(x_\ast) > 0 . \end{align}

Solver

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

Special Cases

Boundary solutions \(x=0\) and \(x=1\)

Before considering special physical limits, it is useful to check when the boundary values \(x=0\) and \(x=1\) are stationary solutions of the local equilibrium equation,

\begin{align} f(x) = K_1 \left[ 1-e^{-Q_1(1-x)} \right] + K_2 \, (1-x) + K_3 \, x(1-x) - K_4 \, x^2 . \end{align}

For a fully neutral cell, \(x=0\), one obtains

\begin{align} f(0) = K_1 \left[ 1-e^{-Q_1} \right] + K_2 . \end{align}

Thus, \(x=0\) is only a stationary solution if the direct and diffuse photoionization terms vanish,

\begin{align} K_1 \left[ 1-e^{-Q_1} \right] + K_2 = 0 . \end{align}

Since these terms are non-negative, this requires

\begin{align} K_2 = 0 \end{align}

and additionally either

\begin{align} K_1 = 0 \qquad \mathrm{or} \qquad Q_1 = 0 . \end{align}

In other words, \(x=0\) can only be a stationary solution if there is effectively no photoionization, neither from the direct stellar EUV field nor from the diffuse recombination field. Whether this boundary solution is stable depends on the derivative

\begin{align} f'(x) = -K_1Q_1e^{-Q_1(1-x)} - K_2 + K_3(1-2x) - 2K_4x . \end{align}

At \(x=0\), this gives

\begin{align} f'(0) = -K_1Q_1e^{-Q_1} - K_2 + K_3 . \end{align}

Because we are interested in stationary solutions, we can combine this with the condition derived above for \(x=0\) to be a root. This requires \(K_2=0\) and either \(K_1=0\) or \(Q_1=0\). Under this condition, the derivative at the boundary reduces to

\begin{align} f'(0) = K_3 = C(T_{\mathrm{gas}}) n_{\mathrm{H}} . \end{align}

Therefore, if collisional ionization is present, \(K_3>0\), a small positive ionized fraction grows away from \(x=0\). In this case, the solution \(x=0\) is unstable. The case \(K_3=0\) is a special limit where the linear term vanishes and nonlinear terms determine the behavior. Physically, \(x=0\) can only be a stable endpoint if there is no photoionization and no collisional ionization.

For a fully ionized cell, \(x=1\), one obtains

\begin{align} f(1) = -K_4 . \end{align}

Thus, \(x=1\) is only a stationary solution if

\begin{align} K_4 = \alpha^{(1)}(T_{\mathrm{gas}})n_{\mathrm{H}} = 0 . \end{align}

This is generally not the case when recombination is active. Recombination therefore drives the solution away from the exact boundary value \(x=1\). If \(K_4=0\), one can also check the stability of this boundary solution. The derivative at \(x=1\) with \(K_4=0\) is

\begin{align} f'(1) = - \left( K_1Q_1 + K_2 + K_3 \right) \leq 0 . \end{align}

Thus, if recombination is absent and at least one ionization process is active, the boundary solution \(x=1\) is stable. However, this requires the special limit \(K_4=0\).

This means that the boundary values \(x=0\) and \(x=1\) are not generally valid equilibrium states. They can occur in special limits, but in the general case the physically relevant stationary solution is expected to lie inside the interval

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

For the numerical solver, this means that exact boundary roots should be treated with care. Unless the corresponding special conditions are fulfilled, the solver should not accept \(x=0\) or \(x=1\) as the final equilibrium state, but should instead search for a stable solution inside the physical interval.

Shielded gas: collisional ionization and recombination only

In very optically thick regions, for example inside the disk, the direct stellar EUV radiation can be completely shielded. The incoming photon density \(u_L\) is zero and therefore

\begin{align} K_1 = 0 . \end{align}

The diffuse recombination radiation term can also be neglected in this limit. If the on-the-spot approximation is used \(u_{\mathrm{rec}}=0\) and therefore

\begin{align} K_2 = 0 . \end{align}

Even without explicitly using the on-the-spot approximation, this limit is still physically plausible in the optically thick disk. The recombination photons have a very short mean free path and are absorbed close to where they are produced. Setting \(K_1=K_2=0\), the local rate equation reduces to

\begin{align} \partial_t x = K_3 \, x(1-x) - K_4 \, x^2 . \end{align}

The two stationary solutions are therefore

\begin{align} x_1 = 0, \qquad x_2 = \frac{K_3}{K_3+K_4}. \end{align}

The solution \(x_1=0\) is an exact mathematical equilibrium because the collisional ionization term is proportional to \(x\). However, if a small finite ionized fraction is present and \(K_3>0\), the collisional ionization term increases \(x\) and the solution moves toward \(x_2\). In this sense, \(x_2\) is the relevant stable equilibrium solution. Inserting \(K_3\) and \(K_4\), the density cancels out:

\begin{align} x_2(T_{\mathrm{gas}}) = \frac{C(T_{\mathrm{gas}})} {C(T_{\mathrm{gas}})+\alpha^{(1)}(T_{\mathrm{gas}})} . \end{align}

Therefore, in this shielded limit the stable ionization state depends only on the gas temperature.

No direct stellar EUV radiation: \(K_1=0\)

A slightly more general analytic limit is obtained if only the direct stellar EUV term vanishes,

\begin{align} K_1 = 0 . \end{align}

This corresponds to a cell that does not receive direct stellar EUV photons, while diffuse recombination photons may still contribute through \(K_2\). The local equilibrium equation then becomes

\begin{align} 0 = K_2 \, (1-x) + K_3 \, x(1-x) - K_4 \, x^2 . \end{align}

Expanding the terms gives

\begin{align} 0 = K_2 + (K_3-K_2)x - (K_3+K_4)x^2 . \end{align}

This is a quadratic equation in \(x\). The two roots are

\begin{align} x_{\pm} = \frac{ (K_3-K_2) \pm \sqrt{ (K_3-K_2)^2 + 4K_2(K_3+K_4) } }{ 2(K_3+K_4) }. \end{align}

The derivative in this limit is

\begin{align} f'(x) = (K_3-K_2) - 2(K_3+K_4)x . \end{align}

Evaluating the derivative at the two roots gives

\begin{align} f'(x_{\pm}) = \mp \sqrt{ (K_3-K_2)^2 + 4K_2(K_3+K_4) }. \end{align}

Thus, \(x_+\) is the stable root, while \(x_-\) is unstable. So \(x_+\) is the physically relevant solution in the interval \(0\le x\le 1\). In the additional limit \(K_2=0\), this result reduces to the shielded collisional-ionization case,

\begin{align} x_- = 0, \qquad x_+ = \frac{K_3}{K_3+K_4}. \end{align}