MA Fynn Wawrzyniak - Revision history

Fynn.W: K1 = 0 Special Case

2026-06-02T14:26:57Z

K1 = 0 Special Case

← Older revision		Revision as of 16:26, 2 June 2026
Line 700:		Line 700:

	[[File:MA FW 05 disk limit x Tgas.png]]		[[File:MA FW 05 disk limit x Tgas.png]]

			=== 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}

Fynn.W: Special Case x=0 and x=1

2026-06-01T12:15:23Z

Special Case x=0 and x=1

Show changes

Fynn.W: Ionization Shielded Gas Case

2026-05-31T22:43:47Z

Ionization Shielded Gas Case

← Older revision		Revision as of 00:43, 1 June 2026
Line 445:		Line 445:
	0 \le x \le 1 .		0 \le x \le 1 .
	\end{align}		\end{align}

			== Solver ==

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

			== Special Cases ==

			=== 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.

			[[File:MA FW 05 disk limit x Tgas.png]]

Fynn.W: M: Rename

2026-05-20T10:56:03Z

M: Rename

← Older revision		Revision as of 12:56, 20 May 2026
Line 387:		Line 387:

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

Line 403:		Line 403:
	\partial_t x		\partial_t x
	=		=
	~~A_E~~		K_1
	\left[		\left[
	1-e^{-~~q_E~~(1-x)}		1-e^{-Q_1(1-x)}
	\right]		\right]
	+		+
	~~A_{\mathrm{rec}}~~ \, (1-x)		K_2 \, (1-x)
	+		+
	~~A_{\mathrm{C}}~~ \, x(1-x)		K_3 \, x(1-x)
	-		-
	~~A_{\mathrm{R}}~~ \, x^2 .		K_4 \, x^2 .
	\end{align}		\end{align}

Line 428:		Line 428:
	f(x)		f(x)
	=		=
	~~A_E~~		K_1
	\left[		\left[
	1-e^{-~~q_E~~(1-x)}		1-e^{-Q_1(1-x)}
	\right]		\right]
	+		+
	~~A_{\mathrm{rec}}~~ \, (1-x)		K_2 \, (1-x)
	+		+
	~~A_{\mathrm{C}}~~ \, x(1-x)		K_3 \, x(1-x)
	-		-
	~~A_{\mathrm{R}}~~ \, x^2 .		K_4 \, x^2 .
	\end{align}		\end{align}

Fynn.W: M: WIP – Local ionization equilibrium solver

2026-05-19T19:27:48Z

M: WIP – Local ionization equilibrium solver

Show changes

Fynn.W: Plots + Initial Ion Fraction

2026-05-12T19:05:53Z

Plots + Initial Ion Fraction

@@ Line 148: / Line 148: @@
 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,
@@ Line 219: / Line 221: @@
 y = NeutralFractionRateEquation(   ustar, urec, Tgas, nH, y);
 </syntaxhighlight>

Fynn.W at 17:41, 12 May 2026

2026-05-12T17:41:11Z

← Older revision		Revision as of 19:41, 12 May 2026
Line 6:		Line 6:
	= Problems =		= Problems =

	== Mean Photon Density in Irradiation Update ( ==		== 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}\).		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}\).

Fynn.W at 17:36, 12 May 2026

2026-05-12T17:36:33Z

← Older revision		Revision as of 19:36, 12 May 2026
Line 219:		Line 219:
	y = NeutralFractionRateEquation( ustar, urec, Tgas, nH, y);		y = NeutralFractionRateEquation( ustar, urec, Tgas, nH, y);
	</syntaxhighlight>		</syntaxhighlight>

	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.

Fynn.W at 17:34, 12 May 2026

2026-05-12T17:34:09Z

← Older revision		Revision as of 19:34, 12 May 2026
Line 10:		Line 10:
	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}\).		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.		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\):		Assuming a constant absorption coefficient inside the cell, the photon density decreases exponentially with distance \(s\):

Fynn.W: Created page with "= 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..."

2026-05-12T17:33:32Z

Created page with "= 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..."

New page

= 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

<syntaxhighlight lang="c">
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);
</syntaxhighlight>

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.