1D H II region and dust reprocessing test

This problem is a one-dimensional, planar H II region test — the plane-parallel analogue of the classical Strömgren sphere — run together with a dust reprocessing test on the same slab. A constant flux of hydrogen-ionizing photons is injected at the middle of a uniform, cold, neutral hydrogen slab and beamed outward in both directions. Each ionization front eats into the neutral gas, photoheats what it ionizes, and slows as it approaches the planar Strömgren column, where recombinations in the ionized column balance the arriving ionizing photons. A second, non-ionizing optical band is injected alongside it: dust absorbs that light, heats until emission balances absorption, and radiates the energy back out in the infrared, where the dust is nearly transparent, so the reprocessed light escapes.

The test therefore exercises photoionization and photoheating driven by a chemistry band, a mixed thermal + chemical multigroup configuration, the directed (beamed) radiation source, and the separate dust-temperature solver, all in one run. The 3D DTypeFront problem is the spherical counterpart.

Radiation groups

Three groups are used, with frequency boundaries \(\{10^{8},\, 10^{14},\, 3.29\times10^{15},\, 10^{19}\}\) Hz:

Chemistry bands must always be the last groups. The two outermost boundaries should be read as \(0\) and \(\infty\) rather than as physical band edges: ComputePlanckEnergyFractions accumulates the Planck integral from zero, so group 0 receives the whole blackbody below \(10^{14}\) Hz whatever the first entry says, and emission above the first chemistry-band boundary is dropped rather than assigned to group 2, so the last entry never enters the emission budget.

Initial conditions

The domain is \(x \in [0, L]\) with \(L = 2 \times 10^{19}\) cm, resolved by 128 uniform cells. The gas is uniform, neutral atomic hydrogen with number density \(n_0 = 100 \, \text{cm}^{-3}\) at \(T_0 = 100\) K. Every radiation group is initialized to the negligible energy-density floor \(E_{\text{rad,floor}}\). The reduced speed of light is \(\hat c = c/1000\) and the run integrates to \(t_{\text{end}} = 2 \times 10^{11}\) s, so each front travels \(\hat c\, t = 6 \times 10^{18}\) cm, about 60% of the way from the source to the wall.

Boundary conditions and source

Reflecting boundaries are applied at \(x = 0\) and \(x = L\), and the source sits at the middle of the domain, \(x_c = L/2\), occupying one cell on each side of it. Because the internal source array carries a luminosity volume density (\(\text{erg}\,\text{s}^{-1}\,\text{cm}^{-3}\)), the injected value for a photon flux \(F\) delivered to each side is

where \(E_\gamma\) is the mean photon energy, \(\Delta x = L/128\), and \(n_s\) is the number of source cells per side (photoionize.source_cells, shipped as 1). The companion reduced flux source is set to \(-1\) in the cells left of \(x_c\) and \(+1\) in those to the right, so each half of the slab injects radiation with reduced flux \(f = F_{\text{rad}}/(c E_{\text{rad}})\) of unit magnitude, free-streaming away from the centre. Setting photoionize.beamed = 0 drops the flux source and makes the injection isotropic instead, which is useful for the comparison in the answer check below.

Both choices — a central source and outward beaming — exist to keep the momentum budget clean. Both walls are reflecting, and a reflecting wall is a momentum source: it turns radiation around, and the reduced speed of light amplifies the bookkeeping value of what it turns around by \(c/\hat c\). Placing the source at the centre and stopping before either front arrives means no wall is ever reached. Beaming matters for a different reason: an isotropic source injects no net momentum at all, so the outward momentum would have to be generated by transport as the M1 closure beams the two wings, and the budget would then measure how fast that happens rather than what the solver did with the momentum.

A thermal group’s source is scaled internally by \(\hat c / c\) and a chemistry band’s is not, so the two shipped fluxes differ by exactly that factor and deliver equal energy: \(F_{\text{opt}} = 10^{13}\) and \(F_{\text{ion}} = 10^{10}\,\text{cm}^{-2}\,\text{s}^{-1}\) per side. Because the hook takes a reduced flux rather than a flux, the solver applies that same scaling when it converts to a real flux source, so unit reduced flux means “beamed” for either kind of band. Photochemistry is enabled, so the ionizing band is not merely transported: it is absorbed by photoionization, and each resulting front stalls at its Strömgren column well inside \(\hat c\, t\).

Ionization front and the Strömgren column

The ionizing band is transparent to the gray dust opacity, so photoionization is the only process that removes it. The network uses a photoionization cross-section of \(1.5 \times 10^{-18}\,\text{cm}^{2}\) and a case-B recombination coefficient

enabled by network.recombination_temperature_dependent = 1. Photoheating is on (network.energy_switch = 1), so the gas the front ionizes is heated toward \(\sim 10^4\) K rather than staying at its initial 100 K.

Each front advances outward until the recombinations behind it consume the arriving ionizing photons. In plane-parallel geometry that balance is \(F_{\text{ion}} = \alpha_B \, n_0^2 \, x_S\), so the equilibrium ionized column per side is

For \(F_{\text{ion}} = 10^{10}\,\text{cm}^{-2}\,\text{s}^{-1}\) and \(n_0 = 100\,\text{cm}^{-3}\) this is \(3.8 \times 10^{20}\,\text{cm}^{-2}\) at \(10^4\) K, falling to \(3.3 \times 10^{20}\) at 8000 K because \(\alpha_B\) is larger there. The measured ionized column at \(t_{\text{end}}\) is \(6.7 \times 10^{20}\,\text{cm}^{-2}\) over both wings, i.e. \(3.3 \times 10^{20}\) per side, which lands inside that range — so each front is at or close to its Strömgren column, the residual depending on how far the ionized gas has been photoheated.

That column sits well inside the light-travel distance: the neutral column out to each light front is \(n_0 \, \hat c \, t = 6 \times 10^{20}\,\text{cm}^{-2}\), nearly twice \(N_S\). The ionization front therefore stalls while a stable fraction of the injected ionizing photons is still in flight, which is what the surviving-fraction bound in the answer check relies on.

Dust model

Each thermal group carries a constant gray dust opacity, \(\kappa_{\text{IR}} = 10\) and \(\kappa_{\text{opt}} = 10^{3}\,\text{cm}^{2}\,\text{g}^{-1}\). With \(\rho = n_0 m_H = 1.66\times10^{-22}\,\text{g}\,\text{cm}^{-3}\) the optical depth from the source out to each front is \(1.0\) in the optical and \(0.01\) in the IR. Dust is therefore optically thick to the incoming light and thin to what it re-emits, which is what makes the reprocessing one-way. Opacity in Quokka is pure absorption, so an opaque group also emits its share of the local blackbody, and that is what supplies the re-emission.

A separate dust temperature is solved for, with the gas–dust collisional coupling switched off (radiation.dust_gas_interaction_coeff = 0). The dust is then fixed purely by radiative equilibrium with the local radiation field,

which gives \(T_d \approx 130\) K for these parameters. At that temperature \(h\nu/(k T_d) = 37\) at the IR/optical boundary, so the Planck function has nothing left above the boundary: essentially all re-emission lands in the IR group and none returns to the optical one. Note that this decoupling is thermal only — radiation momentum is a separate channel and is still deposited, so the radiation accelerates the gas outward and evacuates the cells nearest the source.

Analytic solution

Because the dust returns nothing to the optical band, behind each light front the optical group is in pure attenuation,

so the surviving unprocessed energy, summed over both wings, has the closed form

Everything the optical band loses reappears in the IR, so the two thermal bands together must retain the energy injected into the optical one. The ionizing band is budgeted separately, since photoionization drains it: photons injected equal photons still in the field, plus the ionized column, plus recombinations.

Momentum gives an independent statement. Each wing is injected beamed, so it arrives already carrying \(F E_\gamma / c\) per unit time and area, and the two together inject

Because the source is mirror-symmetric the signed total momentum is zero and says nothing on its own; the informative quantity is the outward momentum, each cell signed by \(\text{sgn}(x - x_c)\). Note the absence of any \(\hat c\) factor here, whereas the energy budget carries \(\hat c / c\), so the two are independent statements.

Running in 3D

The problem is planar, so it runs unchanged in 3D under the same test name and the same input file. The slab is uniform across y and z and the source is a plane at the middle of the x domain, so every transverse column repeats the one-dimensional solution; the point of running it in 3D is to check that the radiation, dust and photochemistry paths behave identically there. inputs/DTypeFront1D.toml serves both builds: only its x entries matter in 1D, and its eight transverse cells – the minimum Quokka allows – span an extent of 8 dx, so the cells stay cubic in 3D.

The domain sums the checks rely on are per unit area, so in more than one dimension they are divided by the number of transverse cells; otherwise each column would be counted once per transverse cell. With that normalization the 3D run reproduces the 1D answers to six or seven significant figures – momentum ratio 1.0038444 against 1.0038441, reprocessed fraction 0.3667441 against 0.3667438, and an identical front position – and takes a few seconds.

Answer check

The test passes if all of the following hold. Measured values at the reference resolution are given in the last column.

Check 2a is the one that exercises the radiation force, and it is a conservation law only because the source is beamed: it measures 1.0038 at 128 cells and 1.0041 at 1024, flat to \(4\times10^{-4}\) over an eightfold refinement. Run with photoionize.beamed = 0 and it falls to 0.78 and drifts with resolution, because the outward momentum then has to be generated by the M1 closure rather than injected. Deleting the gas momentum kick drops it to 0.644. Widening the source slab does not help either: photoionize.source_cells of 1, 2, 4, 8 gives 1.0038, 1.0083, 1.0105, 1.0118 beamed and 0.78, 0.74, 0.63, 0.47 isotropic, and a wide slab also pushes the measured light front outward, so one cell per side is shipped.

The IR band is excluded from check 2a for the same reason the ionizing band is excluded from check 1: nothing injects momentum into it. The dust creates the IR by re-emission, which is isotropic and carries no net momentum, and the outward flux it later develops is generated by transport down the radiation pressure gradient rather than by the source. Outward momentum is not conserved under transport — only the signed total is — so that term has no place in a budget against what was injected. It is not small either: including it would add 9.2%, because the reservoir \(F/(c\,\hat c)\) is inflated by \(c/\hat c\) relative to the physical \(E/c\). The residual \(+0.4\%\) in check 2a is the IR the dust reabsorbs, a fraction \(\tau_{\text{IR}} = \rho\,\kappa_{\text{IR}}\,\hat c\, t\) of the outward IR flux landing in the gas momentum.

Check 3 measures the front one cell beyond \(\hat c\, t\) at this resolution; one cell is 2.6% of the front distance, and the difference falls to 0.7% by 1024 cells. The front is located on the optical band using a threshold of 5% of the unattenuated plateau, rather than the more natural 50%, because dust absorption thins the light by \(e^{-\tau} \approx 0.37\) before it reaches the front; a 50% threshold would report the absorption depth instead.

Check 5 is the H II region check as well as a scaling check. It pins the chemistry-band source scaling, which is not multiplied by \(\hat c / c\) the way a thermal band’s is. The budget inequality 5b alone does not catch a mis-scaling, because fewer photons produce a proportionally smaller ionized column and the inequality still holds; the lower bound in 5a is what catches it, since the ionization front stalls at its Strömgren column well inside \(\hat c\, t\) and always leaves a stable fraction of the photons in flight. Mis-scaling the source by \(c/\hat c\) collapses the surviving fraction to \(8\times10^{-5}\).