# Simulating ionization feedback from young massive stars:
impact of numerical resolution

Yunwei Deng 1,2,3, Hui Li 1,4, Rahul Kannan 5, Aaron Smith 6, Mark
Vogelsberger 2,7 and Greg L. Bryan 4
1Department of Astronomy, Tsinghua University, Beijing 100084, People’s
Republic of China
2Department of Physics, Kavli Institute for Astrophysics and Space Research,
Massachusetts Institute of Technology, Cambridge, MA 02139, USA
3School of Astronomy and Space Science, Nanjing University, Nanjing 210093,
People’s Republic of China
4Department of Astronomy, Columbia University, New York, NY 10027, USA
5Department of Physics and Astronomy, York University, 4700 Keele Street,
Toronto, ON M3J 1P3, Canada
6Department of Physics, The University of Texas at Dallas, Richardson, Texas
75080, USA
7The NSF AI Institute for Artificial Intelligence and Fundamental
Interactions, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Modelling galaxy formation in hydrodynamic simulations has increasingly
adopted various radiative transfer methods to account for photoionization
feedback from young massive stars. However, the evolution of H ii regions
around stars begins in dense star-forming clouds and spans large dynamical
ranges in both space and time, posing severe challenges for numerical
simulations in terms of both spatial and temporal resolution that depends
strongly on gas density ($\propto n^{-1}$). In this work, we perform a series
of idealized H ii region simulations using the moving-mesh radiation-
hydrodynamic code arepo-rt to study the effects of numerical resolution. The
simulated results match the analytical solutions and the ionization feedback
converges only if the Strömgren sphere is resolved by at least $10$–$100$
resolution elements and the size of each time integration step is smaller than
$0.1$ times the recombination timescale. Insufficient spatial resolution leads
to reduced ionization fraction but enhanced ionized gas mass and momentum
feedback from the H ii regions, as well as degrading the multi-phase
interstellar medium into a diffuse, partially ionized, warm ($\sim 8000$ K)
gas. On the other hand, insufficient temporal resolution strongly suppresses
the effects of ionizing feedback. This is because longer timesteps are not
able to resolve the rapid variation of the thermochemistry properties of the
gas cells around massive stars, especially when the photon injection and
thermochemistry are performed with different cadences. Finally, we provide
novel numerical implementations to overcome the above issues when strict
resolution requirements are not achievable in practice.

###### keywords:

H ii regions – methods: numerical – radiative transfer – hydrodynamics –
galaxies: evolution

††pubyear: 2023††pagerange: Simulating ionization feedback from young massive
stars:  impact of numerical resolution–D

##  1 Introduction

Feedback from young massive stars plays a key role in the evolution of giant
molecular clouds (GMCs), galaxies, and the intergalactic medium throughout the
Universe. Born in their natal clouds, massive stars provide feedback to the
ambient interstellar medium (ISM) with ultra-violet (UV) radiation and stellar
winds, and die as core-collapse supernovae (SNe). Ionizing radiation, together
with stellar winds, from these massive stars alters the ionization state of
the ambient gas, and injects comparable thermal energy and momentum as SNe
(e.g. Agertz et al., 2013; Geen et al., 2015; Jeffreson et al., 2021). Such an
energetic early feedback mechanism can disperse GMCs on a short time scale of
$\sim 1.5$ Myr (Kruijssen et al., 2019) prior to the SN explosion (e.g. Fall
et al., 2010; Dale et al., 2012; Dale et al., 2014; Raskutti et al., 2016; Li
et al., 2019) and thus is crucial to regulate the star formation in galaxies
(e.g. Hopkins et al., 2012; Emerick et al., 2018; Hopkins et al., 2020;
Chevance et al., 2022). It is also a key ingredient of galaxy formation
simulations and significantly affects predicted galactic properties (e.g.
Shapiro, 1986; Whalen et al., 2004; Kannan et al., 2014; Kannan et al., 2020a;
Kannan et al., 2022a, b; Garaldi et al., 2022; Smith et al., 2022a).

Recently, with the surging interest in the role of radiative feedback in the
formation and evolution of galaxies, several radiation-hydrodynamics
implementations have been developed (e.g. Aubert & Teyssier, 2008; Petkova &
Springel, 2009; Wise & Abel, 2011; Rosdahl et al., 2013; Jaura et al., 2018;
Kannan et al., 2019; Smith et al., 2020; Chan et al., 2021; Peter et al.,
2023). Sophisticated numerical algorithms have significantly reduced the cost
of solving the radiative transfer (RT) equations, making it more affordable
for cosmological simulations. However, radiative feedback from individual
massive stars originates from very small $\sim 10^{4}$ K fully ionized
regions, namely initial Strömgren spheres (Strömgren, 1939) . Typically, the
mass of the initial Strömgren sphere (initial Strömgren mass) is only on the
order of several solar masses and its formation occurs within a short timespan
of hundreds of years. Driven by these ionized cores, compact H ii regions
expand outward and deposit momentum and kinetic energy to the turbulent ISM
until the massive stars die in a few million years. Capturing such small-scale
physics requires high spatial and temporal resolution and thus tremendous
computational resources, which is often impractical to treat directly in
galaxy formation simulations. Therefore, modelling the impact of individual H
ii regions (photo-heating and expansion) from massive stars is numerically
challenging.

The tiny mass and rapid formation of the initial Strömgren spheres make them
extremely challenging to model in galaxy formation simulations. The largest-
volume cosmological simulations of galaxy formation can only afford mass
resolutions larger than $10^{6}\,\text{M}_{\astrosun}$ for baryonic cells
(e.g. Schaye et al. 2015; Pillepich et al. 2018; Kannan et al. 2022a, 2023;
Hernández-Aguayo et al. 2022, and see Vogelsberger et al. 2020 for a recent
review). Even zoom-in simulations of single Milky Way-like galaxies (e.g.
Rosdahl et al., 2015; Hopkins et al., 2018; Marinacci et al., 2019) usually
have mass resolutions spanning from $10^{3}$–$10^{5}\,\text{M}_{\astrosun}$
per gas resolution element, which is far larger than the initial Strömgren
mass. Although recent dwarf galaxy simulations (e.g. Emerick et al., 2019;
Agertz et al., 2020; Lahén et al., 2020; Gutcke et al., 2021) have begun to
treat feedback from individual stars in dwarf galaxies with a much higher mass
resolution of several solar masses, directly resolving the individual H ii
regions still remains a numerical issue. Scaling inversely with density, the
Strömgren spheres are normally unresolved in gas with $n_{\text{H}}\gtrsim
10\,\text{cm}^{-3}$ for high-resolution galaxy radiation-hydrodynamic
simulations (e.g. Rosdahl et al., 2015; Rosdahl et al., 2022). Even for GMC
simulations, resolving H ii regions can still be challenging in high-density
regimes as both the mass and formation timescale of Strömgren spheres decrease
significantly at high density ($\propto n^{-1}$).

Inadequate resolution often leads to problematic results in galaxy formation
simulations. Historically, the inability to resolve the Sedov-Taylor phase of
SN explosion leads to artificially rapid radiative cooling and forming too
many stars too early in cosmological simulations (Katz, 1992). This so-called
over-cooling problem has been addressed by employing a variety of sub-grid
prescriptions (see Naab & Ostriker 2017 for a recent review). In terms of
ionization feedback, several works have explored the consequences of
insufficient resolution when simulating idealized H ii regions, including the
impact on radial profiles and momentum feedback (e.g. Wise & Abel 2011; Bisbas
et al. 2015; Petkova et al. 2021; Pittard et al. 2022; see also Ivkovic 2023).
However, the physics behind these numerical issues needs to be further
clarified, especially in the context of specific RT methods and applications.
Furthermore, most studies focus solely on spatial resolution and disregard the
influence of temporal resolution. Therefore, there is an urgent need for well-
understood solutions to address these resolution problems to enhance the
convergence of current RT implementations in a more physical manner.

In this article, we perform a suite of radiation-hydrodynamic simulations of
idealized H ii regions with the moment-based M1 closure RT implementation
(Kannan et al., 2019) in the moving-mesh code arepo (Springel, 2010). We study
both the effects of spatial (mass) and temporal (time-stepping) resolution on
the behaviour of simulated H ii regions. For simplicity, our tests are all
performed in uniform pure hydrogen media, but the discussions and conclusions
based on the hydrogen-only case are also valid for the other species by
replacing the coefficients correspondingly. We address three issues arising
from the insufficient spatial (mass) and temporal resolution on stellar
ionization feedback:

  * •

Insufficient spatial (mass) resolution

    1. 1.

Over-ionization: overestimation of the total ionized gas mass due to gas cells
that can not be fully ionized, shifting the ionization-recombination balance
to the ionization side;

    2. 2.

Over-heating: overestimation of energy and momentum deposition due to the
artificial heating of a substantial amount of partially ionized gas;

  * •

Insufficient temporal resolution

    1. 1.

Missing photons: underestimation of the overall feedback due to excess photons
being absorbed without any effect on the ionization state of the gas.

We also emphasize that missing photons can be a very serious problem when
photon injection and thermochemistry have different cadences. We explain the
physical reasons for these numerical issues and describe novel numerical
implementations to solve these issues at an acceptable computational cost.

This paper is organized as follows: We review the physics and analytic results
of H ii region expansion in different phases in Section 2. In Section 3, we
describe our simulation setup and RT implementations. In Sections 4 and 5, we
present the results of our tests for spatial and temporal resolution
dependence, respectively. We then introduce several solutions in Section 6.
Lastly, we discuss and summarize our results in Section 7.

##  2 Anlytical theory of H ii regions

In this section, we review the physics of idealized H ii regions in a uniform
medium. Because of the differences in the dominant physical processes and
timescales, the evolution of an idealized H ii region around a massive star
can be divided into two stages, namely the formation and expansion phase (e.g.
Yorke, 1986). Realistic H ii regions expanding in the turbulent environment
with a roughly power-law density profile (e.g. Shu, 1977; McKee & Tan, 2003;
Lee & Hennebelle, 2018) can expand preferentially toward the rarefied regions
and launch a “champagne flow" that rapidly ionizes and heats the surrounding
cloud (Franco et al., 1990; Zamora-Avilés et al., 2019; Geen et al., 2020,
2021). Detailed discussions of H ii regions in realistic environments are
outside the scope of this work on the impact of numerical resolution, we refer
the reader to the aforementioned literature.

###  2.1 Formation phase

Once a hot massive star ignites in neutral gas, a so-called R(rarified)-type
ionization front (I-front), characterized by a rapid expansion and a
negligible density change across the I-front, is driven through the gas. This
I-front leaves the gas behind the front hot ($\sim 10^{4}$ K) and ionized but
otherwise almost undisturbed. The speed of this I-front is initially
supersonic. It slows down gradually until the radius of the H ii region
reaches the “initial” Strömgren radius.

This idealized problem of a fully ionized, spherical region of uniform density
with the ionization maintained by constantly emitted ionizing photons from a
central massive star is known as a Strömgren sphere (Strömgren, 1939).
Assuming the gas is pure hydrogen, a steady solution for the radius of the
Strömgren sphere (Strömgren radius) can be found by requiring ionization
equilibrium,

$R_{\text{S}}=\left(\frac{3Q}{4\pi
n_{\text{H}}^{2}\alpha_{\text{B}}}\right)^{1/3}\approx
0.315\,Q_{48}^{1/3}n_{3}^{-2/3}{\alpha_{0}}^{-1/3}\,\text{pc}\,,$ (1)

where $\alpha_{\rm B}$ is the effective hydrogen radiative recombination rate
for Case B approximation (Baker & Menzel, 1938), $\alpha_{0}\equiv\alpha_{\rm
B}/(2.59\times 10^{-13}\,\text{cm}^{3}\text{s}^{-1})$, $Q$ is the rate of
emission of ionizing photons with $Q_{48}\equiv Q/(10^{48}\,\text{ s}^{-1})$,
and $n_{\text{H}}$ is the hydrogen density with $n_{3}\equiv
n_{\text{H}}/(10^{3}\,\text{cm}^{-3})$. Since the initial Strömgren sphere is
almost fully ionized, the initial mass of ionized gas, initial Strömgren mass,
is

$M_{\text{S}}=\frac{4}{3}\pi
R_{\text{S}}^{3}n_{\text{H}}m_{\text{H}}=\frac{m_{\text{H}}}{\alpha_{\text{B}}}\frac{Q}{n_{\text{H}}}\approx
3.25\,Q_{48}n_{3}^{-1}{\alpha_{0}}^{-1}\,\text{M}_{\astrosun}\,.$ (2)

Assuming the I-front is infinitely thin ($l_{\text{mfp}}\lesssim
10^{-4}n_{3}^{-1}$ pc), in the non-relativistic limit, the evolution of the
initial Strömgren sphere in the formation phase follows (Spitzer, 1978)

$r_{i}=R_{\text{S}}\left(1-e^{-t/t_{\text{rec}}}\right)^{1/3}\,,$ (3)

where $r_{i}$ is the radius of the I-front and $t_{\rm rec}$ is the
recombination timescale given by

$t_{\text{rec}}=\frac{1}{\alpha_{\text{B}}n_{\text{H}}}\approx
122.3\,n_{3}^{-1}{\alpha_{0}}^{-1}\,\text{yr}\,.$ (4)

Equations (1), (2), and (4) show that the typical initial Strömgren spheres
have sizes under a few parsecs, masses of several solar masses, and form
within several hundreds of years. Such spatial and time scales are much
smaller and shorter than the scales of the subsequent expansion phase in most
astrophysical environments.

###  2.2 Expansion phase

Once the I-front reaches the initial Strömgren radius ($R_{\text{S}}$), the
high pressure of the ionized gas drives the expansion of the H ii region. It
undergoes a rapid transition from R-type to D(dense)-type front: The ionized
gas in photoionization equilibrium expands, driven by the pressure gradient
relative to the neutral gas and causes a shock to separate from the I-front
and proceed into the surrounding neutral gas. The neutral gas thus begins to
accumulate between the shock front and the I-front. This expansion phase
itself can also be divided into two phases: the early phase, during which the
pressure between the ionized gas and surrounding neutral gas is large; and a
later phase when the internal pressure asymptotically balances with the
thermal pressure of the surrounding gas.

Assuming a thin swept-up shell ($v_{i}\sim v_{\text{shock}}$) and applying the
conservation of momentum flux across the I-front and shock front, Spitzer
(1978) gave the classic solution of the time evolution of the expanding
I-front. However, this solution deviates from the expansion at very early time
because it ignores the inertia of shocked gas. Incorporated the inertia,
Hosokawa & Inutsuka (2006) obtained a solution which shows better agreement
with the numerical results (Bisbas et al., 2015):

$r_{i}=R_{\text{S}}\left(1+\frac{7}{4}\sqrt{\frac{4}{3}}\frac{c_{i}t}{R_{\text{S}}}\right)^{4/7}\,.$
(5)

As the Strömgren sphere expands into the uniform medium, its ionized mass
continues to increase following a simple scaling relation (Bisbas et al.,
2015):

$M_{i}(t)=M_{\text{S}}\left(\frac{r_{i}(t)}{R_{\text{S}}}\right)^{3/2}\,.$ (6)

H ii regions embedded in very dense ($\sim 10^{6}$ cm-3) or warm ($>100$ K)
gas can enter their later phase of expansion (Raga et al., 2012). In the later
phase, the expansion of the H ii region is stagnant while the pressure is
equilibrated with the ambient gas so the classic analytic formulas are no
longer valid. In hydrodynamic simulations, there is no clear boundary between
the early and later phases because we follow the gas dynamics self-
consistently. Since the momentum feedback from H ii regions is injected mostly
in the early phase, we will therefore focus on the early phase expansion in
this work.

##  3 Methods

In this section, we briefly describe the setup of the initial conditions and
the terminology used to infer the definition of the mass and temporal
resolution. We also recap the arepo-rt implementation that is relevant for
later discussions.

###  3.1 Initial conditions

To test the dependence of simulation results on the spatial and temporal
resolution, we simulate idealized H ii regions formed by an individual massive
star in uniform pure neutral hydrogen gas. This is a purely radiation-
hydrodynamical test, i.e. gravity and magnetic fields are not included. For
each set of tests, a stellar particle with a steady ionizing photon emission
rate $Q$ is placed at the centre of a box with size $L$ at least 3.6 times
larger than the radius of the H ii region at the end of the simulations. The
gas cells are initially arranged as a regular staggered mesh, which is
constructed by overlapping two Cartesian meshes with a displacement of
$0.45\Delta x$ along each axis from each other. Such a mesh configuration is
adopted to stabilize the construction of the Voronoi tessellation and minimize
density fluctuations in arepo.

We initialize the gas cells as pure hydrogen to directly compare with
analytical results introduced in Section 2. The initial temperature and
density are different in our three tests, we will introduce each of them at
the beginning of the corresponding section.

####  3.1.1 Definitions of mass and temporal resolution

Since the main goal of this paper is to explore the effects of numerical
resolution on the ionization feedback of H ii regions, here we define
explicitly the meaning of the mass and temporal resolution in our numerical
experiments.

In quasi-Lagrangian codes like arepo where the mass of gas cells are close to
the target gas mass $M_{\text{cell}}$, we define the mass resolution
$\mathcal{R}_{i}$ of the initial Strömgren sphere as the ratio between the
initial Strömgren mass and the target gas mass:

$\mathcal{R}_{i}=\frac{M_{\text{S}}}{M_{\text{cell}}}=\frac{m_{\text{H}}}{\alpha_{\text{B}}}\frac{Q}{n_{\text{H}}M_{\text{cell}}}\,.$
(7)

Since the expansion of H ii regions is driven by the ionized gas inside the
initial Strömgren sphere, $\mathcal{R}_{i}$ is critical to determine whether
the Strömgren sphere and its feedback are resolved.

The mass resolution has a direct connection to the physics of the evolution of
H ii regions so it is much more useful and convenient than the spatial
resolution in this work. Still, we define the conversion between mass and
spatial resolution via

$\Delta x=\left(\frac{Q}{\alpha_{\rm
B}n^{2}_{\text{H}}\mathcal{R}_{i}}\right)^{1/3}\,.$ (8)

For convenience, we use the terms ‘spatial resolution’ and ‘mass resolution’
interchangeably.

The temporal resolution is determined by the size of the time integration
steps $\Delta t$. To examine how the results depend on temporal resolution, we
fix the time steps to a given value in each run for all the cells, and both RT
and hydro solvers in the main tests. In realistic simulations, gas cells and
operations (e.g. hydrodynamics and RT) are allowed to have different time
steps, enabling subcycling of processes with short characteristic timescales.
For simplicity, we assume a uniform and fixed time step $\Delta t$ for all
cells and operations throughout the simulation, except in Section 5.1.2 and
related sections. In these sections, we will illustrate how the temporal
resolution issues can be notably amplified due to the discrepancy between the
sizes of photon injection steps ($\Delta t_{\star}$) and RT/thermochemistry
steps ($\Delta t_{\text{RT}}$).

###  3.2 Radiative transfer implementation

We use arepo-rt (Kannan et al., 2019), the radiation-hydrodynamic extension of
the moving-mesh hydrodynamic code arepo (Springel, 2010) to handle the
propagation of radiation and calculation of thermochemistry. arepo-rt adopts a
moment-based RT scheme with the M1 closure relation (Levermore, 1984) to solve
the set of hyperbolic conservation equations for photon density and flux. The
UV radiation continuum is usually divided into three separate bins between
energy intervals of [13.6, 24.6, 54.4, $\infty$) eV mirroring the hydrogen and
helium ionizing photon groups. With these three bins, $4\times 3$ additional
variables are stored in each cell to describe the volume density of photons
$n_{\gamma}^{i}$ and the flux vector ${\bm{F}}_{i}$. The index $i$ here is the
energy band corresponding to ionic species H i, He i, and He ii. For
simplicity, these frequency-dependent coefficients are reduced to their
effective value as only one single $>13.6$ eV bin in our discussion if not
specified. Moreover, the ionization state of each species in each cell is also
stored as $x_{\text{\text{H\,{ii}}}}$, $x_{\text{\text{He\,{ii}}}}$, and
$x_{\text{\text{He\,{iii}}}}$, with $x_{\text{\text{H\,{i}}}}$ and
$x_{\text{\text{He\,{i}}}}$ set by detailed balancing. We adopt the on-the-
spot approximation (OTSA) which assumes all recombination photons are re-
absorbed immediately in the surrounding medium processes (case B, Osterbrock &
Ferland, 2006). We note that although the OTSA may not capture the ionization
profiles accurately in the optically thin regime near the source (Raičević et
al., 2014), it is a reasonable approximation in the high-density mostly
neutral regions where the resolution issues occur (see Section 7.1).

The RT equations are solved with an operator-splitting strategy over a pre-
determined time-step $\Delta t$. Each RT step is decomposed into three sub-
operations executed successively. Firstly, inject photons from the star into
the neighbouring cells. Secondly, solve the RT equations and propagate photons
in space. Lastly, solve the thermochemistry equations to determine the
ionization states and conduct the heating and cooling processes in each cell.
Since arepo uses a two-stage second-order Runge–Kutta integration scheme (i.e.
Heun’s method, Pakmor et al., 2016), the propagation and thermochemistry are
in fact split into two halves within an RT step.

####  3.2.1 Photon injection

In this work, we focus on the ionization feedback from an individual star.
Each test or simulation contains a single stellar particle with a fixed
ionizing photon emission rate $Q$ and ionization spectrum, which corresponds
to an idealized massive star. For a steady stellar source, its ionizing photon
emission rate $Q$ is divided into three frequency bins weighted by a given
spectrum and $\sum_{i}Q_{i}=Q$. At time-step $\Delta t$, the number of
injected photons from a stellar particle in the UV bin $i$ is $\Delta
N_{\gamma}^{i}=Q_{i}\Delta t$. These $\Delta N_{\gamma}^{i}$ photons from the
stellar particle are then weighted among a given number of its nearest
neighbouring gas cells ($N_{\rm nb}=32$ by default). The weight factor $w_{k}$
for cell $k$ is calculated by the solid angle of the cell opened to the
stellar particle. Therefore, the number of ionizing photons dumped into a cell
$k$ close to a star at time-step $\Delta t$ is $\Delta
N_{\gamma,k}^{i}=w_{k}\Delta N_{\gamma}^{i}$ where $\sum_{k}w_{k}=1$.

####  3.2.2 Radiation transport

arepo-rt simulates the propagation of radiation by solving the set of
hyperbolic conservation equations consisting of the zeroth- and first-order
moments, which takes the form of

$\displaystyle\frac{\partial n_{\gamma}^{i}}{\partial
t}+\nabla\cdot{\bm{F}}_{i}=0\,,$
$\displaystyle\frac{\partial{\bm{F}}_{i}}{\partial
t}+\tilde{c}^{2}\nabla\cdot\mathbb{P}_{i}=0\,,$ (9)

where $\mathbb{P}_{i}$ is the pressure tensor related to $n_{\gamma}^{i}$ by
the Eddington tensor. The true speed of light is reduced to $\tilde{c}$ to
prevent too small time-steps, which is known as the reduced speed of light
approximation (Gnedin & Abel, 2001). These equations are closed with the M1
closure relation. There are two frequently-used schemes to evaluate the
interface fluxes, the Harten-Lax-van Leer (HLL, Harten et al., 1983) flux
function and the global Lax-Friedrich (GFL, Rusanov, 1961) function. The GFL
function is more straightforward but makes the RT more diffusive, while the
HLL function has its inherent directionality but makes isotropic radiation
from stars asymmetric (Rosdahl et al., 2013). For the sake of simplicity and
symmetry, we adopt the GFL implementation in this work.

####  3.2.3 Thermochemistry, cooling, and heating

In arepo-rt, to calculate the number density of different ionic species in
each cell, a series of equations are solved by a semi-implicit time
integration approach based on the method outlined in Petkova & Springel (2009)
after the photon injection and transport are finished in each step. Only if
the internal energy changes by more than 10 per cent during a time-step, the
equations will be solved implicitly by calling the functions from the SUNDIALS
CVODE package (Hindmarsh et al., 2005).

For the pure hydrogen case, as the simplest example, the thermochemistry
equation is

$\frac{\mathrm{d}n_{\text{\text{H\,{ii}}}}}{\mathrm{d}t}=-\alpha_{\text{\text{H\,{ii}}}}n_{\text{\text{H\,{ii}}}}n_{e}+\sigma_{e\text{\text{H\,{i}}}}n_{e}n_{\text{\text{H\,{i}}}}+\tilde{c}n_{\text{\text{H\,{i}}}}\sum_{i}\sigma_{i\text{\text{H\,{i}}}}n_{\gamma}^{i}$
(10)

where $n_{\text{\text{H\,{i}}}}$, $n_{\text{\text{H\,{ii}}}}$, and
$n_{\gamma}^{i}$ are the volume densities of neutral, ionized hydrogen, and
photons in bin $i$, respectively, and
$n_{\text{\text{H\,{i}}}}=(1-x_{\text{\text{H\,{ii}}}})n_{\text{H}}$,
$n_{\text{\text{H\,{ii}}}}=x_{\text{\text{H\,{ii}}}}n_{\text{H}}$,
$n_{\gamma}^{i}=N_{\gamma}^{i}/V$. $\sigma_{i\text{\text{H\,{i}}}}$ is the
mean photoionization cross-sections for bin $i$, while
$\sigma_{\text{e\text{H\,{i}}}}$ and $\alpha_{\text{\text{H\,{ii}}}}$ are the
collisional ionization cross-section and recombination rate taken from Katz et
al. (1996). The three terms on the RHS correspond to the number density rates
for H ii recombinations, H i collisional ionizations, and H i
photoionizations, respectively.

When conducting the heating and cooling, the internal energy per unit mass $u$
is stored in each cell. arepo-rt calculates the heating ($\Gamma$) and cooling
($\Lambda$) rate per unit volume with parameters stored in each cell. The
variation of thermal energy $\Delta u$ for a cell in a time-step $\Delta t$ is

$\Delta u=\frac{1}{\rho}(\Gamma-\Lambda)\Delta t\,.$ (11)

The main heating process in the context of idealized H ii regions is
photoionization heating by the deposited energy of the ionizing photons. The
photoheating rates in a finite-frequency bins form can be given as

$\Gamma_{j}=n_{j}\sum_{i}\int_{\nu_{i1}}^{\nu_{i2}}\frac{4\pi
J_{\nu}}{h\nu}\sigma_{j\nu}(h\nu-h\nu_{j})\,\mathrm{d}\nu\,,$ (12)

where $h$ is the Planck constant and $h\nu_{j}$ is the ionization potential of
the ionic species $j$. The total photoheating rate can therefore be given as
$\Gamma=\sum_{j}\Gamma_{j}$.

We mainly consider the following cooling processes in the pure hydrogen H ii
regions: recombination ($\Lambda_{\text{rec}}$), collisional excitation
(bound-bound, $\Lambda_{\text{bb}}$), collisional ionization (bound-free,
$\Lambda_{\text{bf}}$), and Bremsstrahlung (free-free, $\Lambda_{\text{ff}}$).
The cooling equations we use follow Cen (1992) and are listed in Appendix A.
The total cooling rate is therefore calculated by

$\Lambda=\Lambda_{\text{rec}}+\Lambda_{\text{bb}}+\Lambda_{\text{bf}}+\Lambda_{\text{ff}}\,.$
(13)

##  4 Effects of different mass (spatial) resolution

In this section, we study the resolution effects on the evolution of idealized
H ii regions in their formation and expansion phases, respectively.

###  4.1 Test 1 – Formation phase and over-ionization

In the first test, we investigate the impact of resolution on the formation of
the initial Strömgren sphere. To do this, we perform static Strömgren sphere
tests by placing a steady source of hydrogen ionizing photons with a rate of
$Q=10^{63}$ Gyr${}^{-1}=3.17\times 10^{46}$ s-1 (unity rate in arepo-rt) in a
uniform pure hydrogen neutral medium of number density $n_{\text{H}}=1$ cm-3.
The effective H i ionization cross-section
$\sigma_{\text{\text{H\,{i}}}}=3.0\times 10^{-18}$ cm2. For this test, we turn
off the cooling and heating of gas and the temperature is fixed at $T=10^{4}$
K. The speed of light is reduced to $\tilde{c}=0.001c$. The chosen initial
condition prioritizes simplicity, while the results can be scaled to any
arbitrary initial conditions using the dimensionless resolution
$\mathcal{R}_{i}$. To test the effect of mass resolution $\mathcal{R}_{i}$, we
set a series of initial conditions with cell mass from about
$10^{-3}\,\text{M}_{\astrosun}$ to $10^{3}\,\text{M}_{\astrosun}$, equivalent
to $\mathcal{R}_{i}$ from $10^{5}$ to 0.1 or $\Delta x$ from $0.34$ pc to $34$
pc. All the cells are initially at rest and arranged in a regular staggered
mesh.

Figure 1: Structure of the Strömgren spheres (projected H ii fractions) in
different mass resolution from the highest (${\cal R}_{i}=10^{5}$) to lowest
(${\cal R}_{i}=0.1$) resolution simulations (Test 1, Section 4.1). The high-
resolution Strömgren spheres are nearly spherical, demonstrating the ability
of our code to obtain well-defined H ii regions with correct Strömgren radius
that is consistent with the analytical result as shown in the blue dotted
circles. Notice that the spatial axes are rescaled for each map, while the H
ii fractions are shown at the same scale.

In Figure 1, we present the projected maps of the H ii fractions for
simulations with different resolutions from ${\cal R}_{i}=10^{5}$ to ${\cal
R}_{i}=0.1$ resolution with a spacing of $10^{2}$. The ${\cal R}_{i}=10^{5}$
and ${\cal R}_{i}=1000$ Strömgren spheres exhibit an ideal spherical
morphology, indicating the ability of our code to capture the characteristic
features of H ii regions. Compared with the analytic Strömgren radius
($R_{\text{S}}$, blue dotted circles), the high-resolution Strömgren spheres
are almost fully ionized and concentrated inside $r\lesssim R_{\text{S}}$ with
a sharp boundary. On the contrary, the low-resolution Strömgren spheres are
much larger than $R_{\text{S}}$, filled with partly ionized gas with much
lower H ii fractions.

Figure 2: Results of the formation test (Test 1, Section 4.1). Top (middle)
panels: profiles of neutral (ionized) fraction as a function of radius for
different mass resolutions at $t\approx 15\,t_{\text{rec}}$; Bottom panel:
evolution of ionized mass as a function of time for different mass
resolutions.

In the top (middle) panels of Figure 2 we show the steady neutral (ionized)
fraction profiles at $t\approx 15\,t_{\text{rec}}$. The profile of the
ionization fraction can be described by (Osterbrock & Ferland, 2006)

$\frac{x_{\text{\text{H\,{i}}}}(r)}{4\pi
r^{2}}Qe^{-\tau(r)}\sigma_{\text{\text{H\,{i}}}}=x_{\text{\text{H\,{ii}}}}^{2}(r)n_{\text{H}}\alpha_{\text{\text{H\,{ii}}}}\,,$
(14)

where $x_{\text{\text{H\,{i}}}}(r)=n_{\text{\text{H\,{i}}}}(r)/n_{\text{H}}$,
$x_{\text{\text{H\,{ii}}}}(r)=1-x_{\text{\text{H\,{i}}}}(r)$, the optical
depth
$\tau(r)=n_{\text{H}}\sigma_{\text{\text{H\,{i}}}}\int_{0}^{r}x_{\text{H\,{i}}}(r^{\prime}){\rm
d}r^{\prime}$, and the cross-section $\sigma_{\text{H\,{i}}}$ here is $3\times
10^{-18}$ cm2. This analytic solution is plotted in Figure 2 with the grey
dashed curves as a benchmark. The analytic profile is well-reproduced only in
the $\mathcal{R}_{i}\geq 10^{5}$ simulation. With decreasing resolution, the
profile becomes more and more smooth, the fractional ionized region becomes
more extended while the fully ionized core region shrinks. In high-resolution
cases, the ionization fraction $x_{\text{\text{H\,{ii}}}}$ reaches unity in
the fully ionized core and there is no significant difference until about
$\mathcal{R}_{i}=10$. However, once the H ii region transitions to being
marginally resolved and then unresolved, the situation is very different. The
effective size of one cell is larger than the initial Strömgren sphere so that
the total ionized mass is only a small portion of the mass of one cell.
Consequently, the ionization fraction is reduced to $\sim
M_{\text{S}}/M_{\text{cell}}={\cal R}_{i}$. As we presented in Figure 2,
$x_{\text{\text{H\,{ii}}}}$ falls down quickly to $\lesssim 0.1$ in the
$\mathcal{R}_{i}=0.1$ case.

Although the ionization fractions are averaged down and ionization profiles
are smoothed away, the total mass of ionized gas $M_{i}$ is not conserved
during the resolution smoothing. Assuming the Strömgren sphere is fully
ionized, the time evolution of the I-front follows equation (3) and the total
ionized mass can be described as

$M_{i}(t)=M_{\rm S}\left(1-e^{-t/t_{\text{rec}}}\right)\,.$ (15)

In the bottom panel of Figure 2, we present the evolution of $M_{i}$ for
different resolutions $\mathcal{R}_{i}$, comparing with the analytic result of
equation (15). There is no significant difference among the
$\mathcal{R}_{i}=10^{5}$, $10^{4}$, and $10^{3}$ simulations, and the initial
Strömgren sphere can be regarded as well-revolved when the mass resolution is
larger than $100$, at least from the perspective of the ionized mass. However,
the total mass of ionized gas is enhanced in lower-resolution cases. In
marginally resolved cases like $\mathcal{R}_{i}=10$ or $1$, the equilibrium
mass of ionized gas is slightly larger than the actual mass, while in the
unresolved case of $\mathcal{R}_{i}=0.1$ it is more than ten times larger. We
refer to such overestimation of total ionized gas mass in low resolution cases
as over-ionization. The time for the Strömgren sphere to reach a steady state
also becomes longer with decreasing resolution, especially for the
$\mathcal{R}_{i}\leq 1$ simulations.

####  4.1.1 Over-ionization

The time lag to reach a steady state is not very important to the subsequent
evolution of H ii regions, since most processes we are interested in occur on
time-scales far longer than the formation of the initial Strömgren sphere.
However, the over-ionization issue we mentioned before can significantly
increase the total mass of ionized gas in low-resolution galactic and
cosmological simulations if we adopt a direct RT feedback implementation.

Equation (10) regulates the thermochemistry in the pure hydrogen medium, and
in the context of UV ionization, the second term on the RHS, corresponding to
the collisional process, can be neglected, while the third term can be
simplified to a single H i ionization photon band. The balance equation of the
ionization state of the gas can thus be written as

$\alpha_{\text{B}}n_{\text{\text{H\,{ii}}}}n_{e}\approx\tilde{c}n_{\text{\text{H\,{i}}}}\sigma_{\text{\text{H\,{i}}}}n_{\gamma}\,,$
(16)

where
$n_{\text{\text{H\,{ii}}}}=1-n_{\text{\text{H\,{i}}}}=n_{e}=x_{\text{\text{H\,{ii}}}}n_{\text{H}}$.
Substituting the number densities with ionization fractions, we obtain

$\alpha_{\text{B}}n^{2}_{\text{H}}x_{\text{\text{H\,{ii}}}}^{2}\approx\tilde{c}(1-x_{\text{\text{H\,{ii}}}})n_{\text{H}}\sigma_{\text{\text{H\,{i}}}}n_{\gamma}\,,$
(17)

where the LHS corresponds to the recombination rate $R_{\text{rec}}$ and the
RHS corresponds to the photoionization rate $R_{\text{ion}}$.

In high-resolution cases, we suppose that the initial Strömgren sphere can be
resolved by a series of thin shells. For a cell in the $k$th shell with radius
$r_{k}$, the equilibrium solution of photon density is
$n_{\gamma}=Qe^{-\tau(r_{k})}/4\pi r_{k}c$, where
$\tau(r_{k})=n_{\text{H}}\sigma_{\text{\text{H\,{i}}}}\sum_{1}^{k}x_{\text{\text{H\,{i}}}}(r_{k})\Delta
r$ and $\Delta r\sim\Delta x$. With sufficiently high resolution, equation
(17) will approach the analytical solution given by equation (14). However, in
the low-resolution cases with mass resolution $\mathcal{R}_{i}<1$, the cell
size $\Delta x$ will be larger than $R_{\text{S}}$. Assuming the entire
initial Strömgren sphere is embedded in a single cell, the equilibrium balance
equation becomes

$\alpha_{\text{B}}n^{2}_{\text{H}}x^{2}_{\text{\text{H\,{ii}}}}\approx
Q/V_{\text{cell}}$ (18)

substituting $V_{\text{cell}}=M_{\text{cell}}/n_{\text{H}}m_{\text{H}}$ with
the definition of mass resolution ${\cal R}_{i}$ (equation 7), we have

$x_{\text{\text{H\,{ii}}}}\approx{\cal R}_{i}^{1/2}.$ (19)

Thus, in fully unresolved cases (${\cal R}_{i}\ll 1$) the equilibrium
ionization fraction is proportional to ${\cal R}_{i}^{1/2}$ and the mass of
ionized gas is $\sim M_{\text{cell}}{\cal R}_{i}^{1/2}=M_{\text{S}}{\cal
R}_{i}^{-1/2}$, overestimated by a factor of ${\cal R}_{i}^{-1/2}$. Therefore,
over-ionization is a direct consequence of recombinations and ionizations
attempting to become balanced in a Strömgren sphere diluted by the impact of
insufficient spatial resolution.

In other words, in low resolution cases, the ionization fraction is supposed
to saturate at $x_{\text{\text{H\,{ii}}}}\sim\mathcal{R}_{i}$ to obtain the
correct ionized mass $M_{S}$. However, with this ionization fraction, the
ionization rate $R_{\text{ion}}$ is still much larger than the recombination
rate $R_{\text{rec}}$ at this time. More gas will be ionized before the
$R_{\text{ion}}$ and $R_{\text{rec}}$ are finally balanced at
$x_{\text{\text{H\,{ii}}}}\sim\mathcal{R}_{i}^{1/2}$, resulting in the over-
ionization. In section 6.1.2, we will introduce a method to prevent this issue
by enforcing the correct balance between recombination and ionization rates.
In section 7.1.1, we will discuss the conditions for which this correction is
needed.

Similar resolution issues make a significant difference in the expansion phase
of the H ii region as well – we will turn to the early time expansion and
discuss these numerical issues in the next section.

###  4.2 Test 2 – Early time expansion and over-heating

In the next series of tests, we simulate the expansion of an H ii region in a
uniform density medium. The simulation box is initialized with pure hydrogen
gas of density $n_{\text{H}}=100$ cm-3 and temperature $T=100$ K. A steady
source is placed at the centre of the box that emits a blackbody spectrum with
$T_{\text{eff}}=35\,000$ K at a rate of $Q=10^{48}$ s-1, which is similar to a
$\sim 30\,\text{M}_{\astrosun}$ massive star and the effective H i ionization
cross-section is $3.8\times 10^{-18}$ cm2. We use a reduced speed of light
with $\tilde{c}=0.01\,c$ and a multifrequency RT scheme as described in
Section 3.2. We vary the mass resolution $\mathcal{R}_{i}$ across a range of
$0.01$ to $1000$ corresponding to $M_{\text{cell}}(\Delta x)$ from
$4800\,\text{M}_{\astrosun}(12\,\text{pc})$ to
$0.048\,\text{M}_{\astrosun}(0.27\,\text{pc})$. To temporally resolve the
formation of the initial Strömgren sphere, we force the hydrodynamic time-step
to be shorter than $0.1\,t_{\text{rec}}$ (0.12 kyr).

Figure 3: Results of the expansion test (Test 2, Section 4.2). Shock
parameters at 1 Myr: the ionization (top left panel), temperature (top middle
panel), density (top right panel), pressure (bottom left panel), and velocity
(bottom middle panel) profiles at 1 Myr with different resolutions. Figure 4:
Results of the expansion test (Test 2, Section 4.2). Evolution of expanding H
ii regions: Top left panel: the Strömgren radius as a function of time for
simulations with different resolutions; Top right, bottom left, and bottom
right panels: the time evolution of ionized mass, momentum injection, and
energy injection with different resolutions.

In Figure 3, we show the profiles of ionization fraction, temperature,
density, pressure, and velocity at $1$ Myr for all the simulations in this
test. At this time, the Stömgren sphere has expanded to $r_{i}\approx 5R_{\rm
S}$. In high-resolution simulations like the $\mathcal{R}_{i}=1000$ run, the
profiles show characteristics of a typical D-type front. After the D-type
expansion begins, the I-front is always preceded by a hydrodynamical shock
front sweeping outward. The interior of the H ii region is evacuated by the
expansion and becomes much more rarefied than the surrounding neutral gas. The
$\mathcal{R}_{i}=10$ simulation retains these features with visually
acceptable flattening, while lower resolutions gradually flatten these away.
In the $\mathcal{R}_{i}=0.1$ simulation, the gas density is almost uniform,
implying that the shock is not captured.

Some delicate features can be captured in the highest-resolution simulation,
like a secondary pressure peak behind the shock front and a knee of
temperature in the transient region. However, only the $\mathcal{R}_{i}\geq
100$ simulations can reproduce these secondary structures, as the
$\mathcal{R}_{i}=10$ simulation has already smoothed them out. Nonetheless,
the profiles can be roughly reproduced in marginally-resolved cases
($\mathcal{R}_{i}=10$ and $1$). In simulations conducted at lower resolutions,
the initial Strömgren sphere is completely unresolved. The insufficient
resolution causes all the profiles to appear flatter and results in incomplete
ionization mirroring the scenario discussed in Section 4.1. Moreover, a
substantial fraction of ionized gas also fails to be heated to $10^{4}$ K, and
the density inside the I-front stays within the same order of magnitude as the
background density, and the pressure becomes much higher.

In Figure 4, we present the evolution of the ionization radius, ionized mass,
injected momentum, and injected energy. The ionization radius $r_{i}$ (top
left panel) is defined as the radius of the shell where
$x_{\text{\text{H\,{i}}}}=x_{\text{\text{H\,{ii}}}}$ since the profile of the
I-front is flattened. The $\mathcal{R}_{i}\leq 1$ runs have ionization
fractions that never exceed 0.5 so they are unable to be plotted in this
panel. The grey dashed curves in the top-left and top-right panels indicate
the Hosokawa-Inutsuka solution (equation 5). The ionization radius $r_{i}$
matches well with the analytic solution even for the marginally-resolved case
($\mathcal{R}_{i}=1$), and the ionized mass tends to converge to this solution
as the H ii region expands. However, the over-ionization issue makes a
significant difference in the ionized mass at very early times, as discussed
in Section 4.1. A surging of the ionized mass can be found in the
$\mathcal{R}_{i}=10$ to $0.01$ results, and it takes an increasing amount of
time to relax to the analytic solution with decreasing $\mathcal{R}_{i}$. This
suggests that the expansion of H ii regions can counteract the overestimation
of ionized mass caused by the over-ionization to some extent. Even so, the
expansion only increases the ionized mass slowly and it can take a
considerable time, even longer than the lifetime of the massive star, to
counteract the initial surging if the over-ionization is too severe.

In the bottom panels of Figure 4 we present the cumulative radial momentum and
injected energy summed over the entire simulation box. In the low-resolution
cases, the momentum and energy injections are also enhanced. Unlike the
overestimated ionized mass which tends to asymptote to the analytic solution,
the momentum and energy remain consistently larger than the high-resolution
solution. This is because once momentum and energy are over-injected, no
mechanism can eliminate the excess momentum or energy, unlike the ionization
state of the gas, which can regulate itself by balancing the photoionization
and recombination process.

####  4.2.1 Over-heating

In this section, we illustrate that the enhancement of momentum and energy
injected in low-resolution cases is a result of the artificial heating of a
substantial amount of partially ionized gas, raising its temperature to
several thousand Kelvin. We refer to this phenomenon as over-heating.

The hot gas in the initial Strömgren sphere is heated by photoionization and
reaches equilibrium with cooling processes at $T\gtrsim 10^{4}$ K. Then the
expansion of the H ii region is driven by the large pressure gradient between
the $10^{4}$ K hot ionized gas and the surrounding background neutral gas. As
the expansion is a hydrodynamic response to the initial ionization, the over-
heating issue is a secondary disaster of the over-ionization issue we have
discussed in Section 4.1.1.

When the star ignites, the neighbouring gas cells will begin to be ionized and
the deposited energy of ionization photons is converted to thermal energy in
each step, following equation 12. On the other hand, the hot gas will lose
internal energy efficiently through various cooling processes. The difference
between photoheating and radiative cooling results in a net increase of the
internal energy (equation 11).

In low-resolution cases, the over-ionization issue leads to an overestimation
of the mass of ionized gas but this gas can only be partially ionized because
ionization balance cannot be established properly (Section 4.1). Similarly, we
now discuss the balance between cooling and heating in the initial Strömgren
sphere. In equilibrium, the heating and cooling rates satisfy

$\Gamma=\Lambda=\Lambda_{\text{rec}}+\Lambda_{\text{bb}}+\Lambda_{\text{bf}}+\Lambda_{\text{ff}}\,,$
(20)

where each rate has the following dependence on the ionization fraction
(equations 49, 50, 51, and 52):

$\Gamma\propto(1-x_{\text{\text{H\,{ii}}}})\,,\;\Lambda_{\rm rec/ff}\propto
x^{2}_{\text{\text{H\,{ii}}}}\,,\;\text{and}\;\Lambda_{\text{bb/bf}}\propto
x_{\text{\text{H\,{ii}}}}(1-x_{\text{\text{H\,{ii}}}})\,.$ (21)

In low-resolution cases ($\mathcal{R}_{i}\ll 1$),
$x_{\text{\text{H\,{ii}}}}\sim\mathcal{R}_{i}$ and the heating and cooling
rates can be approximated as follows:

$\Gamma\propto 1\,,\;\Lambda_{\rm
rec/ff}\propto\mathcal{R}_{i}^{2}\,,\;\text{and}\;\Lambda_{\text{bb/bf}}\propto\mathcal{R}_{i}\,.$
(22)

The heating rate stays high in the gas cells even though it should drop down
in high-resolution cases where these cells are fully ionized. On the other
hand, the cooling rates are all reduced by insufficient ionization, especially
those for recombination and Bremsstrahlung.

This scenario bears resemblance to the over-ionization issue, yet an
additional numerical concern exacerbates matters. Besides the ionization
fraction, all cooling functions involved are sensitive to the temperature
(equations 49, 50, 51, and 52), while the temperature of a gas cell is
calculated based on the internal energy $u$ and mean molecular weight $\mu$ by

$T=\frac{2u}{3k_{\text{B}}\mu m_{\text{H}}}.$ (23)

In our pure hydrogen case, the mean molecular weight equals to

$\mu=1/(1+x_{\text{\text{H\,{ii}}}})\sim
1-\mathcal{R}_{i}\,,\;\text{when}\,\mathcal{R}_{i}\ll 1\,.$ (24)

If the initial ionization cannot be resolved, the over-ionization will lead to
a larger ionized mass but lower ionization fraction, thus higher molecular
weight $\mu$ and lower temperature $T$. In the $\mathcal{R}_{i}\ll 1$
unresolved cases the ionization fraction $x_{\text{\text{H\,{ii}}}}$ falls
from unity to $\sim 0$ so that $T$ can be reduced by half for the same
internal energy. For example, in the $\mathcal{R}_{i}=0.1$ case, the molecular
weight is 0.91 and the temperature of the ionized gas only reaches 0.55 of
that found in the high-resolution case, assuming the same net energy
injection. As we presented in Figure 3, in the $\mathcal{R}_{i}=10$ to $1000$
simulations, the central temperature reaches $\gtrsim 1.5\times 10^{4}$ K,
while in unresolved cases it falls to less than $8000$ K.

The most efficient cooling mechanism in the $10^{4}$ K primordial ISM is the
bound-bound collision (equations 51), which has the largest rate coefficients
($7.5\times 10^{-19}$ erg cm-3 s-1). However, it is also the most sensitive to
the temperature

$\Lambda_{\text{bb}}\propto\frac{e^{-118348/T}}{1+T^{1/2}_{5}}$ (25)

Therefore, $\Lambda_{\text{bb}}$ will be drastically reduced if the
temperature cannot reach its proper value. For instance, if the temperature
inside the H ii region is $T=8000$ K, the bound-bound collision cooling rate
will be reduced with a factor of $5.5\times 10^{-4}$ compared to the
$T=16\,000$ K case.

Consequently, all four primordial cooling mechanisms fail to cool the gas
under the combined effect of the above two factors. The total cooling rate is
thus drastically underestimated while the heating rate is enhanced. Therefore,
the injected thermal energy will be significantly overestimated and leads to
an over-heating problem.

It is crucial to avoid over-ionization and over-heating problems in
simulations because they both lead to divergent feedback when H ii regions are
unresolved, i.e. unresolved H ii regions will provide more feedback rather
than less feedback, which hinders the convergence of simulations with
different resolutions. In Section 7.1.1, we will discuss the specific physical
conditions in which over-ionization and over-heating problems become
particularly relevant.

####  4.2.2 Multi-phase gas

Figure 5: Results of the expansion test (Test 2, Section 4.2). Mass-weighted
phase diagrams for the $\mathcal{R}_{i}=1000$, $10$, and $0.1$ simulations in
Test 2 at $t=1$ Myr.

In Figure 5, we present the mass-weighted phase diagram for the
$\mathcal{R}_{i}=1000$, $10$, and $0.1$ simulations at $t=1$ Myr. The gas
shows a continuous distribution over the phase space in the high- and medium-
resolution simulations, which separates into several “islands” at low
resolution. The gas density in the high-resolution case spans a range of two
orders of magnitude. The inflection at $\sim 10^{2.5}$ cm-3 corresponds to the
shock-compressed gas that has accumulated in front of the D-type I-front.
Lowering the resolution flattens this inflection, and the gas in the ${\cal
R}_{i}=0.1$ run is distributed close to the vertical line $n_{\text{H}}=100$
cm-3 (initial density) with $T<10^{4}$ K. This suggests that the shock
transition zone and the compressed gas layer between the shock and I-front are
not properly resolved.

Figure 6: Results of the expansion test (Test 2, Section 4.2). Top panel:
mass-weighted histograms of gas mass as a function of temperature for the
$\mathcal{R}_{i}=1000$, $10$, and $0.1$ simulations in Test 2 at $t=1$ Myr.
Middle and bottom panels: the H ii fraction and cooling rate as functions of
temperature obtained from these simulations.

In Figure 6 top panel, we present mass-weighted histograms of the gas mass as
a function of temperature for the high ($\mathcal{R}_{i}=1000$), medium (10),
and low (0.1) resolution runs at 1 Myr. As a comparison, we also present the H
ii fraction (middle panel) and cooling rate (bottom panel) as functions of
temperature obtained from these runs.

In the high- and medium-resolution runs, the temperature has a continuous
distribution on the $100$ K to $16\,000$ K interval of temperature and there
are four populations of the different gas phases. The largest population on
the left is the undisturbed background gas. The hottest gas populated from
about $\sim 10\,000$–$15\,000$ K is the fully ionized gas inside the Strömgren
sphere. The population at $\sim 8000$–$10\,000$ K is partially ionized and
there is a ramp of heated warm neutral gas ranging from $\sim 100$–$8000$ K.
These two populations are the shock-heated gas between the I-front and shock
front during the D-type expansion, and they are accumulating as a bump around
$9000$ K because of the inefficient cooling. The cooling rate reaches its
maximum value at $10\,500$ K, which results in a valley between the two bumps
of populations.

At lower resolutions, the temperature distribution becomes discrete and the
valleys are unpopulated. In the $\mathcal{R}_{i}=0.1$ run, only the $<2500$ K
neutral population and $\sim 8000$ K low ionization population still exist and
the $T>10^{4}$ K fully ionized population disappears. This gives an intuition
for the impact of the over-ionization and over-heating issues, as discussed in
Section 4.2.1. The neutral gas fails to become fully ionized and heat to
$>10^{4}$ K but stops at $\sim 8000$ K with a low ionization fraction, which
is inefficient for cooling. Notice the y-axis is in log-scale so the total
mass of ionized gas at low resolution is actually larger than that at high
resolution though the H ii fraction is low at that temperature.

In summary, the insufficient spatial resolution flattens out the shock
structures and the over-heating issue changes the multi-phase gas structure of
the ISM by turning the highly-ionized, hot ($>10^{4}$ K) gas within H ii
regions into partially-ionized, warm ($\sim 8000$ K) gas. Such errors in the
gas phase structure can be problematic in line luminosity estimates (Smith et
al., 2022b; Tacchella et al., 2022). Failure to resolve the low-density bubble
evacuated by the ionization feedback and the swept-up shell leads to an
inability to reproduce the enhanced and weakened effects of SN explosions
(e.g. Walch & Naab, 2015; Haid et al., 2016). This leads to uncertainties in
the total energy/momentum injection from stellar feedback and the amount of
gas blown out of a galaxy.

##  5 Effects of different Temporal Resolution

Choosing a proper time step is crucial to obtain correct simulation results
and minimise computational cost. We now examine the dependence of the
simulation results on the choice of the time step. arepo-rt is a moment-based
RT implementation treating photons with finite speed of light. The most well-
known time-stepping criterion for convergent time integration is the Courant
criterion (Courant et al., 1928), and it is modified to involve the reduced
speed of light (Kannan et al., 2019)

$\Delta t_{\rm C}<\eta\,\frac{\Delta
x}{\tilde{c}+|{\bm{v}}_{\text{cell}}|}\,,$ (26)

where $\Delta x$ and ${\bm{v}}_{\text{cell}}$ are the width and velocity (in
the lab frame) of a cell and $\eta\sim 0.3$ is the Courant factor. In
expanding H ii regions, the velocity of a D-type shock is $v_{s}\lesssim
2c_{i}\ll\tilde{c}$ so

$\displaystyle\Delta t_{\rm C}$ $\displaystyle\approx\eta\frac{\Delta
x}{\tilde{c}}=\frac{\eta}{\tilde{c}}\left(\frac{Q}{\alpha_{\rm
B}n_{\text{H}}^{2}\mathcal{R}_{i}}\right)^{1/3}$ (27) $\displaystyle\approx
230\,\text{yr}\,\left(\frac{\eta}{0.3}\right)\left(\frac{\tilde{c}}{10^{-3}c}\right)^{-1}\left(\frac{{\cal
R}_{i}}{10}\right)^{-1/3}Q_{48}^{1/3}n_{3}^{-2/3}\,.$

It turns out in practice that such a criterion ensures a stable convergence
for hydrodynamic simulations. However, since the mean free path (MFP) of
photon ionization is very hard to spatially resolve ($\Delta
x<l_{\text{mfp}}=1/\sigma_{\text{\text{H\,{i}}}}n_{\text{H}}$), the Courant
criterion almost always fails to catch the ionization timescale

$\displaystyle t_{\rm
ion}=\frac{1}{\tilde{c}n_{\text{\text{H\,{i}}}}\sigma_{\text{\text{H\,{i}}}}}=\frac{\alpha_{\rm
B}}{\tilde{c}\sigma_{\text{\text{H\,{i}}}}}t_{\rm rec}\,$ (28)
$\displaystyle\approx 1.4\times
10^{-3}\left(\frac{\tilde{c}}{10^{-3}c}\right)^{-1}t_{\rm
rec}=0.17\,\text{yr}\left(\frac{\tilde{c}}{10^{-3}c}\right)^{-1}n_{3}^{-1}\,.$

In a forming H ii region, the ionization state of gas is changing drastically,
and in a microscopic view, the gas ionization timescale is essentially the
$t_{\text{ion}}$ presented above. Such a short timescale is difficult to
temporally resolve, which can lead to errors in the solutions. We emphasize
that this is a general issue across numerical RT methods and explicit time-
stepping schemes though we introduce it from a moment-based scheme with a
finite speed of light.

###  5.1 Missing photons issue

In most thermochemistry solvers (including arepo-rt), to reduce the
computational cost and ensure numerical stability, the photon transport and
the thermochemistry are solved with an operator-splitting strategy (e.g.
Petkova & Springel, 2009; Wise & Abel, 2011; Rosdahl et al., 2013; Kannan et
al., 2019; Chan et al., 2021). In the thermochemistry networks, the photon
number density $n_{\gamma}$ is not coupled in as a variable but is considered
to be a known quantity as the solution of the last photon transport step. The
absorption of photons is simply calculated independently of solving the
thermochemistry equation with the optical depth in a cell by

$\frac{{\rm d}n_{\gamma}}{{\rm
d}t}=-\tilde{c}n_{\text{\text{H\,{i}}}}\sigma_{\text{\text{H\,{i}}}}n_{\gamma}=\frac{n_{\gamma}}{t_{\text{ion}}}\,,$
(29)

where $n_{\text{\text{H\,{i}}}}$ takes the value at the beginning of the
current step. Once a package of photons is dumped or propagated into a neutral
gas cell, the solution for the above equation is simply

$n_{\gamma}(\Delta t)=n_{\gamma,0}\,e^{-\Delta t/t_{\text{ion}}}\,,$ (30)

where $n_{\gamma,0}$ is the photon density at the beginning of the timestep.

In this paper, we refer to these simplified thermochemistry solvers where
photon density is not coupled in the chemistry network as “uncoupled solvers”.
This simplification has the advantage of efficiency especially when multiple
RT bins and transfer of momentum from radiation to gas are implemented and it
has been commonly used in various RT implementations. However, we will show it
can lead to a severe numerical error by overestimating the number of absorbed
photons, turning them to “missing photons”.

####  5.1.1 Missing photons when radiation encounters neutral gas

Figure 7: Number of photons absorbed ($\Delta N_{\gamma}$) in a fully neutral
neighbouring cell during the first step after the star ignites, as a function
of the time step size $\Delta t$ (assuming the sizes of all the relative time
steps match each other). The red and green curves are the accurate solution
obtained with the coupled solver (see Section 6.2.1) and the problematic
solution obtained with the uncoupled solver. The pink-filled region presents
the missing photons described by equation (33).

When the star ignites, ionizing photons are dumped into the predominantly
neutral gas and propagate outwards. We use a one-zone case to demonstrate the
missing photons issue when photons are absorbed in a fully neutral
neighbouring cell of $n_{\text{H}}=10^{9}$ cm-3 during the first step after a
$Q=10^{48}$ s-1 star ignites. In Figure 7, we plot the number of photons
absorbed ($\Delta N_{\gamma}$) in a fully neutral neighbouring cell during the
first step after the star ignites, as a function of the time step size $\Delta
t$ (assuming the sizes of all the relative time steps match each other). The
red curve is obtained by an accurate thermochemical solver solving the photon
and ion numbers simultaneously (see Section 6.2.1), we regard this solution as
the accurate one for our question. On the other hand, the green curve is
obtained by the widely used uncoupled solver.

In each step when radiation encounters neutral gas, the number of photons
absorbed in this step is (equation 30)

$\Delta N_{\gamma}=N_{\gamma}(1-e^{-\Delta t/t_{\rm ion}})\approx
N_{\gamma}\,,$ (31)

where $N_{\gamma}$ is the number of photons prior to absorption, $\Delta
t/t_{\rm ion}\gg 1$ in most cases so the second equality holds. This equation
demonstrates that almost all photons dumped or propagated into the gas will be
consumed when the $t_{\text{ion}}$ (or MFP) is unresolved. As shown by the
green curve in Figure 7, the number of absorbed photons converges to the
number of injected photons $Q\Delta t/N_{\text{nb}}$ when the time step size
is larger than several $t_{\text{ion}}$.

On the other hand, the total number of neutral hydrogen atoms in a cell is

$N_{\text{\text{H\,{i}}}}=\frac{M_{\rm cell}}{m_{\rm H}}=\frac{Q}{\alpha_{\rm
B}n_{\rm H}{\cal R}_{i}}\,.$ (32)

Assuming one photon can ionize one neutral hydrogen atom (MFP unresolved), in
this case, if $\Delta N_{\gamma}>(1+f_{\text{rec}})N_{\text{\text{H\,{i}}}}$,
then $\Delta N_{\gamma}$ photons will be absorbed but only
$N_{\text{\text{H\,{i}}}}$ atoms can be ionized, the photon number is no
longer conserved. Here $f_{\text{rec}}$ is a factor accounting for the
recombination of ions. In other words, $\delta N_{\gamma}$ photons simply
disappear without any effect on the gas state, they become the missing photons
(pink-filled region in Figure 7). Specifically, the number of missing photons
is

$\delta N_{\gamma}=\Delta
N_{\gamma}-(1+f_{\text{rec}})N_{\text{\text{H\,{i}}}}\,,$ (33)

where $f_{\text{rec}}$ can be estimated as $\Delta t/t_{\text{rec}}$ because
$\Delta t\sim t_{\text{rec}}\gg t_{\text{ion}}$. If the size of the time step
is too large, $\Delta N_{\gamma}\propto Q\Delta
t\gg(1+f_{\text{rec}})N_{\text{\text{H\,{i}}}}$, almost all the photons are
lost, having no effect on the ionization state when they enter a neutral gas
cell.

In Appendix B, we track the step-by-step events for the missing photons issue
by considering an idealized scenario (Figure 15). We find that the following
criterion should be satisfied to avoid the occurrence of missing photons if
the photon density is not coupled in the thermochemistry solver:

$\Delta t\lesssim t_{\text{rec}}\frac{N_{\rm nb}}{{\cal
R}_{i}}\qquad\left(\text{or}\quad{\cal R}_{i}\lesssim N_{\rm
nb}\quad\text{if}\quad\Delta t\gtrsim t_{\text{rec}}\right)\,.$ (34)

This criterion is stringent because it suggests the mass resolution should not
be higher than $N_{\text{nb}}$ or the size of the time step should be far
smaller than $t_{\text{rec}}$. Nonetheless, the missing photons issue may not
be as troublesome as long as the photon injection and thermochemistry have the
same cadence. This only occurs when a large number of photons are being dumped
or propagated into a substantially neutral cell. However, realistic
implementations usually conduct photon injection, radiative transport, and
thermochemistry with different cadences. In the next section, we will show
that the missing photons issue can lead to disastrous consequences in these
cases.

####  5.1.2 Mismatching cadences make things worse

Figure 8: Similar to Figure 7, but assume $\Delta t_{\star}=2\Delta
t_{\text{RT}}$ and evolve it to the end of the second injection step after the
star ignites. The red and green curves are the cumulative number of photons
absorbed in these two inject steps (four RT steps) obtained by the accurate
and uncoupled solver. The slightly thinner curves are the number of photons
absorbed in the four RT steps with the uncoupled solver, the summation of
these four curves is the green line.

In realistic simulations, it is common to encounter mismatching cadences
between photon injection, radiative transport, and thermochemistry. Typically,
the frequency of photon injection is determined by the star particle time
step, while the RT and thermochemistry time steps of the neighbouring cells
can be further refined. Consequently, photon dumping is typically performed
with a relatively low cadence compared to RT, i.e., each photon injection step
can be associated with several RT steps. Further subcycling for thermochemical
integration is also common when the gas properties change drastically in a
single RT step, leading to a more severe mismatch between photon injection and
thermochemistry cadences.

To illustrate the problem, we consider the case of multiple subcycles of RT
and thermochemistry conducted within one injection time step (as described in
Kannan et al., 2019). In our tests, we use an equal time-stepping scheme, so
the injection time step of the star $\Delta t_{\star}$ is equal to the
hydrodynamic time step $\Delta t$, and the subcycled RT step $\Delta
t_{\text{RT}}=\Delta t/N_{\text{sub}}$. Even if one does not intentionally
subcycle the RT, the two-stage Heun’s integration method introduces at least
two “substeps” for RT, making this a potential issue for the similar variant
of second-order implementations, regardless of whether they adopt explicit RT
subcycling or not.

In Figure 8, we present a similar one-zone test as the one in the last section
but now assuming $\Delta t_{\star}=2\Delta t_{\text{RT}}$. We evolve the
system to the end of the second injection step after the star ignites. It is
important to note that the missing photons issue now occurs not only in the
first RT step (‘RT step 1’) immediately after the first injection of photons
but also in subsequent injection steps. The difference is that fresh photons
are not injected immediately after this RT step. However, in the remaining
second RT step (‘RT step 2’), the photon density $n_{\gamma}$ is almost $0$ if
the missing photons issue has occurred. As a result, ionization processes
become negligible and the gas begins to recombine. The H ii fraction before
the next injection step $x_{\text{inj}}$ can be estimated as

$x_{\text{inj}}=\frac{1}{1+\Delta t_{\text{RT}}/t_{\text{rec}}}\,.$ (35)

Consequently, the gas will not be fully ionized even after the occurrence of
the missing photons issue. The partially neutral gas will keep eliminating
photons that enter this cell in subsequent time steps. Similar to equation
(33), the number of missing photons in the first substep of the next sourcing
step (denoted as step 2) is

$\displaystyle\delta
N_{\gamma,2}=N_{\gamma,2}\left[1-e^{-(1-x_{\text{inj}})\Delta
t_{\text{RT}}/t_{\text{ion}}}\right]$
$\displaystyle-\left[1-x_{\text{inj}}+f_{\text{rec}}(\Delta
t_{\text{RT}})\right]N_{\text{\text{H\,{i}}}}\,,$ (36)

where $N_{\gamma,2}$ is the number of photons in the second step before
absorption.

As we present in Figure 8 with the ‘RT step 3’ curve, the second injection
step will meet a severe missing photons issue if $\Delta t_{\star}\gtrsim
0.1t_{\text{rec}}$. This leads to all $2Q\Delta t_{\star}$ photons injected in
these two steps being absorbed after the third RT step (‘RT step 3’), leaving
almost no photons in the gas. If the number of survived photons is
insufficient to suppress the recombination in the fourth RT step (‘RT step
4’), the H ii fraction will decrease again as equation (35). The missing
photons issue has now led to an instability: the H ii fraction oscillates
between 1 and $x_{\text{inj}}$ (equation 35) with a period of two RT steps and
the photon density does not increase after the first occurrence of missing
photons but remains near zero. In the most severe case (e.g. the
$t_{\star}/t_{\text{rec}}\gtrsim 0.2$ part of Figure 8), all photons injected
subsequently suffer the same severe missing photons issues caused by the
unstable H ii fraction as that in the first two injection steps. Even though a
fraction of photons may escape, they potentially also meet the same issue in
subsequent neutral cells. Consequently, a substantial amount of photons
emitted by the star can be dissipated (see Figure 16) and lead to an incorrect
H ii region and significantly reduced feedback.

Recall that $\Delta t_{\text{C}}$ is typically $10^{3}$ times larger than
$t_{\text{ion}}$ (equations 27 and 28), therefore even a relatively high but
not fully ionized H ii fraction can consume a large fraction of photons in the
first RT step. In Figure 16, we present the number of surviving photons as a
function of $\Delta t_{\text{RT}}$, where the survival ratio decreases rapidly
as $\Delta t_{\text{RT}}$ increases.

Generalizing equation (35), if we have $N_{\text{sub}}>1$ but only inject
photons once at the beginning, the H ii fraction after subcycling is
$x_{\text{inj}}=1/[1+(N_{\text{sub}}-1)\Delta t_{\text{RT}}/t_{\text{rec}}]$.
Thus, we obtain a time-stepping criterion $\Delta t_{\star}\lesssim
0.1t_{\text{rec}}N_{\text{sub}}/(N_{\text{sub}}-1)$, which can ensure
$x_{\text{inj}}>0.9$ to avoid oscillations in the H ii fraction and repeated
missing photons. We validate this criterion quantitatively in Appendix C, and
we also show this leads to converged results in Section 5.2. This also
suggests that having a different cadence between photon injection and
thermochemistry may lead to serious consequences, but increasing the number of
subcycles should not alter the behaviour significantly.

###  5.2 Test 3 – Consequences of repeating missing photons issue

In this experiment, we simulate the formation and expansion of an H ii region
in the uniform medium with a density of $10^{9}\,\text{cm}^{-3}$. Such a high
density is set to break the Courant criterion and it has been observed in
hyper compact H ii regions (e.g. Moscadelli et al., 2021). The simulation box
is initialized with pure hydrogen gas of temperature $T=100$ K and the
radiation source is the same as that used in Section 4.2 with $T=35\,000$ K
and $Q=10^{48}$ s-1. The reduced speed of light is set as $0.001c$ for a
larger $\Delta t_{\rm ion}$. The mass resolution $\mathcal{R}_{i}$ is set as
$10^{4}$ to exclude the effects of low spatial resolution. We minimize the
number of RT subcycles as $N_{\text{sub}}=2$, which is an intrinsic
requirement of the second-order scheme of arepo. We use a forced equal time-
stepping scheme so that the RT time steps for all the cells are the same
throughout the simulation and equal to the system integration time step.

With this simulation setup, we will present the reduced feedback issue when
the repeating missing photons issue happens and test the maximum available
time-step to obtain correct momentum feedback with an uncoupled solver. The
recombination time is $t_{\text{rec}}=1.22\times 10^{-4}$ yr such that the
ionization and Courant times are
$\\{t_{\text{ion}},t_{\text{C}}\\}=\\{2.23\times
10^{-3},18.9\eta\\}\,t_{\text{rec}}=\\{2.77\times 10^{-7},2.31\times
10^{-3}\eta\\}$ yr.

Figure 9: Results of the temporal resolution test (Test 3, Section 5.2). The
left panel shows the radial profiles of the H ii fraction for the $\Delta
t_{\star}=0.5\,t_{\rm ion}$ (fiducial) run and $\Delta
t_{\star}=\\{0.1,0.2,0.4,0.8,1.6,3.2\\}\,t_{\rm rec}$ runs, and the right
panel shows the evolution of injected momentum over $1000\,t_{\rm rec}$.

To determine the maximum available time step for simulating the feedback from
the H ii region, we run simulations with different time steps $\Delta t$ to
investigate the effect of temporal resolution. We focus on the radial momentum
injection as a quantitative way to evaluate the impact on feedback.

In the left panels of Figure 9, we present radial profiles of the H ii
fraction and temperature from the $\Delta t_{\star}=0.5\,t_{\rm ion}$
(fiducial) and $\Delta t_{\star}=\\{0.1,0.2,0.4,0.8,1.6,3.2\\}\,t_{\rm rec}$
runs at $t=1000\,t_{\rm rec}$. In the right panel, we show the evolution of
the cumulative injected momentum over $1000\,t_{\rm rec}$. We find that only
runs with a time-step of $\Delta t_{\star}\leq 0.2\,t_{\rm rec}$ can correctly
reproduce the fiducial result. When increasing the time-step to $0.4\,t_{\rm
rec}$, the Strömgren sphere becomes smaller and the momentum injection is
reduced by a factor of $1.5$. With larger $\Delta t_{\star}$, the size of the
H ii region and momentum injection suddenly plunge to below half of the
fiducial solution. This is because the $0.8\,t_{\rm rec}$ and larger $\Delta
t$ runs fall due to the severe repetition of the missing photons issue, as
presented in Section 5.1.2.

In summary, we find that $\Delta t_{\star}\lesssim
0.1\,t_{\text{rec}}N_{\text{sub}}/(N_{\text{sub}}-1)$ could be a useful
criterion for the RT time-step to obtain the correct H ii region feedback
under the current implementation of arepo-rt. A larger RT time step will
result in too little momentum injection, ionization, and heating.

We note that while our ideal experiment can present convergent results with a
time step that satisfies this criterion, in real simulations, gas cells are
also allowed to have different time steps. This creates the possibility of
photons flowing from a cell with a longer time step into one with a shorter
time step, resulting in the occurrence of recurring missing photon issues.
Therefore, it is crucial to devise a method to correct this issue. In section
6.2, we will discuss how to deal with this difficulty with the photon density
coupled thermochemistry solver and a corrected uncoupled solver.

##  6 Correcting for inadequate resolution

In this section, we will introduce several methods to correct these spatial
and temporal resolution problems.

###  6.1 Spatial correction

####  6.1.1 Lowering number of neighbours during photon injection

In Section 4, we have found that ${\cal R}_{i}\gtrsim 10$ is needed to ensure
the correct ionized mass and feedback while lower resolution leads to
increased ionized mass and enhanced feedback. The simplest way to correct the
resolution issues is to reduce the number of neighbouring cells that photons
from the star are injected into.

In the photon injection routine, $\Delta N_{\gamma}$ photons from the star
particle are dumped into a given number of its nearest neighbouring gas cells
in a weighted fashion. In our default setup, it is set to 32, in order to
enshroud the stellar particle. If one does not care about the morphology of
individual H ii regions but only their cumulative feedback to the ISM and
galaxy, the spatial resolution problem can be partly solved by lowering the
value of $N_{\rm nb}$ to $N_{\rm nb}<{\cal R}_{i}$. In practice, $N_{\rm nb}$
can be simply set as $2$ to maximize the probability of resolving the
Strömgren spheres (e.g. Kannan et al., 2020b).

However, this correction cannot eliminate the possibility of over-ionization
so that highly unresolved H ii regions will still lead to divergent feedback.
Also, this correction is a global method which is hard to adjust dynamically
based on the local resolution. Therefore, we will introduce a novel correction
method in the next section, which corrects the resolution effects and provides
convergent results by dynamically reducing the ionization rate in the
neighbouring cells.

####  6.1.2 Accurate ionization and heating balances

In this section, we describe an alternative way to correct spatial resolution
issues based on enforcing the correct balance between recombination and
ionization rates (equation 17).

Recall in Section 4, when we discussed the resolution effects in low-
resolution cases, we assumed that the entire initial Strömgren sphere is
embedded in a single cell. In this case, when the ionized mass reaches the
Strömgren mass $M_{\text{S}}$, the ionization fraction in this cell only
reaches ${\cal R}_{i}$ rather than $\sim 1$, which leads to an ionization rate
$R_{\text{ion}}$ larger than the recombination rate $R_{\text{rec}}$. Our
modification aims to reduce the ionization rate by a factor $f_{\text{cor}}$
so that the ionization and recombination are balanced once the ionized mass
reaches $M_{\text{S}}$, i.e.

$\alpha_{\text{B}}n^{2}_{\text{H}}x_{\text{\text{H\,{ii}}}}^{2}=f_{\text{cor}}\tilde{c}(1-x_{\text{\text{H\,{ii}}}})n_{\text{H}}\sigma_{\text{\text{H\,{i}}}}n_{\gamma}\,,$
(37)

when $M_{i}=M_{\text{S}}$.

In practice, photons from the star are dumped into $N_{\text{nb}}$ cells with
a weight factor $w_{k}$ for each neighbour, e.g. based on solid angle
weighting. Here we assume that photons from the star are always trapped in the
nearest neighbours of the star when the resolution is lower than ${\cal
R}_{i}=N_{\text{nb}}$. In this case, each neighbouring cell $k$ can thus be
regarded as an individual unresolved H ii region with a source rate of
$Q_{k}=w_{k}Q$. By definition of the mass resolution (equation 7), we have
${\cal R}_{i,k}$ for cell $k$ of

${\cal
R}_{i,k}=\frac{m_{\text{H}}}{\alpha_{\text{B}}}\frac{w_{k}Q}{n_{\text{H}}M_{\text{cell}}}=w_{k}{\cal
R}_{i}\,.$ (38)

When the initial Strömgren sphere is formed, the ionized gas of mass $M_{i}$
is hosted by these $N_{\text{nb}}$ cells, and
$M_{i}=\sum_{k}x_{\text{\text{H\,{ii}}},k}M_{\text{cell}}$ where
$x_{\text{\text{H\,{ii}}},k}$ is the ionization fraction of cell $k$. Because
$\sum_{k}w_{k}=1$, we have $x_{\text{\text{H\,{ii}}},k}\approx{\cal R}_{i,k}$
at the time when $M_{i}=M_{\text{S}}$. This is similar to the idealized one-
zone case we discussed in Section 4, except that we use ${\cal R}_{i,k}$ to
account for every neighbour. Replacing $x_{\text{\text{H\,{ii}}}}$ with ${\cal
R}_{i,k}$ in equation (37), we have the target balance equation for each cell

$\alpha_{\text{B}}n_{\text{H}}{\cal R}_{i,k}^{2}\approx
f_{\text{cor}}\tilde{c}\sigma_{\text{\text{H\,{i}}}}n_{\gamma}\,.$ (39)

Therefore, the correction factor $f_{\text{cor}}$ that is required to
establish an accurate balance once $M_{i}$ reaches $M_{\text{S}}$ is

$f_{\text{cor}}=\frac{\alpha_{\text{B}}n_{\text{H}}{\cal
R}_{i,k}^{2}}{\tilde{c}\sigma_{\text{\text{H\,{i}}}}n_{\gamma}}\,.$ (40)

The steps to implement this correction are listed below:

  1. 1.

We use an attribute ${\cal R}_{i,l}$ to store the local mass resolution of
each gas cell $l$ and initialize it as ${\cal R}_{i,l}=-1$ at the beginning of
each step.

  2. 2.

When looping through neighbouring gas cells around the star in the photon
injection routine, evaluate the local mass resolution ${\cal R}_{i,k}$ for
each cell $k$ based on its weight factor $w_{k}$ and gas attributes (equation
38).

  3. 3.

In the thermochemistry routine, for each cell $l$, we check whether $0<{\cal
R}_{i,l}<1$. If so, we go to the next step to conduct the resolution
correction. Otherwise the cell is resolved (${\cal R}_{i,l}>1$) or not a
neighbouring cell (${\cal R}_{i,l}=-1$).

  4. 4.

If cell $l$ requires a resolution correction, active the correction by
checking if $x_{\text{\text{H\,{ii}}},l}>{\cal R}_{i,l}$, calculate the
correction factor $f_{\text{cor}}$ (equation 40), and reduce the ionization
and heating rates by multiplying them with $f_{\text{cor}}$.

In this way, the spatial resolution correction is implemented dynamically only
for cells in the neighbour lists of the stars and only based on their own gas
attributes. This correction algorithm can thus be easily applied to realistic
simulations with non-uniform gas density. The analogous correction factors for
more species can also be obtained by modifying the detailed balance equations.
We derive and describe the correction factors for He and H2 in Appendix D.

Figure 10: Evolution of the ionized mass from low mass resolution static
Stömgren sphere simulations (Test 1, Section 4.1) with (dash curves) and
without (solid curves) our correction

In Figure 10, we present the low resolution runs in the static Strömegren
sphere test (Test 1, Section 4.1) with (dash curves) and without (solid
curves) our correction. With the correction, both the marginally resolved
(${\cal R}_{i}=1$ and $100$) and unresolved (${\cal R}_{i}=0.1$ and $0.01$)
Strömegren spheres can obtain equilibrium ionized masses within a factor of
$2$ of the expected analytic result. We only conduct resolution corrections
for cells in the neighbour list of the star, the ionized mass is increasing
over time slowly because a portion of photons leaks out of these cells.
However, we will soon see that this only leads to a minor effect on the
overall momentum feedback.

Figure 11: The momentum (top) and energy (bottom) feedback at $1$ Myr obtained
with (black triangles) and without (blue circles) the spatial resolution
correction as a function of spatial resolution $\mathcal{R}_{i}$. The magenta
triangle shows the results of the ${\cal R}_{i}=0.01$ run combining reducing
the number of neighbours to 2 and our spatial resolution correction. The
y-axis quantities are normalized to the fiducial results (${\cal
R}_{i}=1000$).

In Figure 11, we present the momentum and energy feedback with and without the
spatial resolution correction at $1$ Myr. For the marginally-resolved cases
($0.1\lesssim{\cal R}_{i}\lesssim 10$), the momentum and energy feedback after
correction perfectly match the high resolution results. In the case of
significant unresolved simulations (${\cal R}_{i}=0.01$), the inadequate mass
resolution results in an artificial enhancement of momentum and energy
feedback by over one order of magnitude arising from the effects of over-
ionization and over-ionization. With the resolution correction, the momentum
and energy feedback are much under control and are even reduced by a factor of
a few compared to the true solution. This is because our correction not only
ensures the correct initial Strömgren mass but also partly fixes the over-
heating issue by limiting the photoheating rate and correcting the balance
between heating and cooling so the internal energy is approximately conserved
(this is not guaranteed because the reduced cooling rate is not corrected). In
return, the unresolved Strömgren spheres reach even lower temperatures,
several hundred Kelvin for the ${\cal R}_{i}=0.01$ run, leading to
insufficient pressure gradients to drive the expansion by converting thermal
energy to kinetic energy. Therefore, the highly unresolved H ii regions simply
dissolve into the background ISM and provide less amount of feedback.

To maximize the possibility of resolving ionization feedback, we reduce
$N_{\text{nb}}$ to 2 (Section 6.1.1) along with our resolution correction
(Section 6.1.2). The magenta triangle in Figure 11 shows the corrected results
of the ${\cal R}_{i}=0.01$ run combining these two methods. By concentrating
the deposition of ionizing photons in fewer cells, we achieve more accurate
momentum feedback and energy injection compared to only using the correction.

####  6.1.3 Caveat: Trade-off between correcting ionization feedback and
achieving accurate recombination line emission

Although enforcing the correct balance between recombination and ionization
rate can help us obtain the correct ionized mass and feedback from unresolved
H ii regions, this can lead to additional issues in the post-processing
prediction of recombination line emission, as it preserves the conversion of
ionizing photons into line photons (e.g. Smith et al., 2022b).

Without our correction, the balance between the total ionizing flux and
integral recombination rate

$Q=\int R_{\text{rec}}{\rm
d}V=4\pi\int\alpha_{\text{\text{H\,{ii}}}}n_{\text{H}}^{2}x_{\text{\text{H\,{ii}}}}^{2}r^{2}{\rm
d}r\,,$ (41)

which holds for all resolutions because the photon number is conserved.
However our correction equivalently reduces $Q$ with a factor of ${\cal
R}_{i,k}$ to ensure

$M_{i}=4\pi m_{\text{H}}\int n_{\text{H}}x_{\text{\text{H\,{ii}}}}r^{2}{\rm
d}r=M_{\text{S}}\,.$ (42)

Thus, there is a trade-off between obtaining the correct ionized
momentum/energy on the fly based on the correct $M_{i}$ and obtaining precise
recombination line emission with post-processing based on the recombination
integral (equation 41).

Figure 12: Radial integration of the recombination rate (RHS of equation 41,
normalized to the ionizing photon rate $Q$) as a function of radius obtained
from the runs presented in Figure 10. The grey dashed lines indicate the value
of ${\cal R}_{i,k}$, to which the corrected runs asymptote.

In Figure 12, we present the radial integration of recombination rate as the
function of radius obtained from the runs presented in Figure 10. The
recombination integrals in the uncorrected runs converge to unity at all the
resolutions while those in the corrected runs decrease with a scaling
$\propto{\cal R}_{i,k}$. Consequently, employing our correction will lead to
an underestimate in line emissivity with a factor of ${\cal R}_{i,k}$.
Correcting ionization feedback while simultaneously achieving accurate
recombination line emission presents a challenging trade-off. A possible
solution is to use the ${\cal R}_{i,k}$ values stored in the gas cells to
further correct the recombination rate in the cells where the spatial
resolution correction is activated. Ongoing investigations aim to improve
solutions to this dilemma.

###  6.2 Temporal correction

####  6.2.1 Coupling photon density into thermochemistry calculations

The most direct (but expensive) solution to solve the temporal resolution
problems is coupling the photon density into the ODE system of thermochemistry
calculations. For hydrogen-only cases with a single UV bin, there are three
variables: the H ii fraction $x_{\text{\text{H\,{ii}}}}$, the photon density
normalized to its initial value $r_{\gamma}=n_{\gamma}/n_{\gamma,0}$, and
specific internal energy $u$. Thus, we need to solve the following three
equations simultaneously

$\displaystyle\frac{\mathrm{d}x_{\text{\text{H\,{ii}}}}}{\mathrm{d}t}$
$\displaystyle=-\alpha_{\text{\text{H\,{ii}}}}n_{\text{H}}x^{2}_{\text{\text{H\,{ii}}}}+\sigma_{\text{e\text{H\,{i}}}}n_{\text{H}}x_{\text{\text{H\,{i}}}}x_{\text{\text{H\,{ii}}}}+\tilde{c}\sigma_{\text{\text{H\,{i}}}}x_{\text{\text{H\,{i}}}}n_{\gamma,0}r_{\gamma}\,,$
$\displaystyle\frac{\mathrm{d}r_{\gamma}}{\mathrm{d}t}$
$\displaystyle=-\tilde{c}n_{\text{H}}\sigma_{\text{\text{H\,{i}}}}x_{\text{\text{H\,{i}}}}\,,\quad\text{and}$
(43) $\displaystyle\frac{\mathrm{d}u}{\mathrm{d}t}$
$\displaystyle=\frac{1}{\rho}\Bigr{[}\bar{\Gamma}x_{\text{\text{H\,{i}}}}r_{\gamma}-(\bar{\Lambda}_{\text{ff}}+\bar{\Lambda}_{\text{rec}})x^{2}_{\text{\text{H\,{ii}}}}-(\bar{\Lambda}_{\text{bb}}+\bar{\Lambda}_{\text{bf}})x_{\text{\text{H\,{i}}}}x_{\text{\text{H\,{ii}}}}\Bigr{]}\,,$

where the $\bar{\Gamma}$ is the constant part of the photoionization heating
rate and $\bar{\Lambda}$ is the temperature-dependent part of the heating and
cooling rates listed in Appendix A, respectively.

In principle, the photon density coupled thermochemistry solver should give an
accurate solution no matter how long the integration step is. In Figure 13, we
present the result of Test 3 (Appendix 5.2) solved by a coupled solver (red
curve). The integration time step is decided by the code as $\Delta
t=7.8\,t_{\text{rec}}$. This solution matches the $\Delta
t=0.5\,t_{\text{ion}}$ fiducial run (grey dash curve) perfectly except at the
very beginning, where the too-long time step leads to a delay in the formation
of the initial Strömgren sphere. As a comparison, the green curve shows the
solution of the uncoupled solver with the same time step, whose momentum
feedback is erroneously reduced by over two orders of magnitude.

Figure 13: Evolution of momentum injection for Test 3 (Section 5.2) runs
obtained with the coupled solver (red), the uncoupled solver without
correction (green), and the uncoupled solver with limited absorption
correction (blue), using a Courant time step decided by the code. The grey
dashed curve is the fiducial result obtained with a time step $\Delta
t=0.5\,t_{\text{ion}}$.

####  6.2.2 Limited photon absorption in an uncoupled solver

Figure 14: Left panel: Similar to Figure 7, but including results from the
uncoupled solver with the limited absorption correction (blue curve), which is
in agreement with the accurate coupled solver (red curve) avoiding the missing
photons issue seen in the uncoupled solver without correction (green curve).
Right panel: Similar to the left panel but $\Delta t=\Delta t_{\star}=\Delta
t_{\text{RT}}$ and the y-axis is normalized to the number of photons dumped in
an injection step with size $\Delta t_{\star}$. The red curve shows the
accurate solution of the coupled solver and the blue (green) dash, dot, and
solid curves are the solution of the corrected (uncorrected) uncoupled solver
for the first and second RT steps, and their summation, respectively.

Although the coupled solver provides an accurate solution for the temporal
resolution problem, it can be unnecessarily expensive if He and H2 are
included and adding frequency bins accounting for their photo-chemistry. Here
we introduce an approximate solution to correct the missing photons issue with
an uncoupled solver, based on the suggestion in Jaura et al. (2020).

In Section 5.1.1, we have seen that for a neutral cell, the maximum ability to
absorb photon is $(1+f_{\text{rec}})N_{\text{\text{H\,{i}}}}$, where
$f_{\text{rec}}=\Delta t/t_{\text{rec}}$. Jaura et al. (2020) proposed a
simple method to solve this problem, where they force the number of absorbed
photons to be smaller than or equal to the number of atoms that can be
ionized. For hydrogen, the number of photons absorbed in a thermochemistry
step is

$\Delta
N_{\gamma}=\min{}\left[(1+f_{\text{rec}})N_{\text{\text{H\,{i}}}}\,,N_{\gamma,0}e^{-\Delta
t/t_{\text{ion}}}\right]\,,$ (44)

where $N_{\text{\text{H\,{i}}}}$ and $N_{\gamma,0}$ are the number of neutral
hydrogen and the number of hydrogen ionizing photons in a gas cell,
respectively.

In the left panel of Figure 14, we present the number of photons absorbed in a
fully neutral neighbouring cell during the first step after the star ignites,
as a function of the size of the time step (with a setup of Test 3, assuming
the sizes of all the relative time steps match each other). The solution of
the coupled solver (red curve) is regarded as the accurate solution. Compared
to this, the uncoupled solver (green curve) deviates from the solution when
the time step is larger than $t_{\text{ion}}$ exhibiting a severe missing
photons issue, while the limited photon absorption correction (equation 44,
blue curve) matches the solution in both large and small $\Delta t$ intervals.

In Section 5.1.2, we discussed how the difference in cadence between the
photon injection and thermochemistry routines can exacerbate the missing
photons issue by repeatedly consuming photons in subsequent RT substeps. As
the green lines in Figure 14 show, if $\Delta t_{\star}$ is larger than a few
$0.01t_{\text{rec}}$, almost all photons injected will be absorbed in the
first RT step, leaving the gas photon-free to recombine in the second RT step.
This issue is also solved by the limited photon absorption correction. As
shown in the right panel of Figure 14, our correction method ensures that only
a small fraction of the total number of photons dumped in an injection step
$\Delta t_{\star}$ are absorbed in the first RT substep. This allows a
substantial number of photons to remain in the second substep and suppress gas
recombinations. Therefore, the total number of photons absorbed during an
injection step shows a perfect agreement with the accurate solution,
demonstrating the effectiveness of our correction.

In Figure 13, the blue curve shows the momentum feedback obtained with a
corrected uncoupled solver, which presents a perfect match to the result of
the coupled solver (red curve). Such a perfect correction is based on the fact
that the recombination rate, as a function of temperature, is not changing
dramatically in our idealized situation. For H ii regions, this holds in most
cases because they usually have an equilibrium temperature of about $10^{4}$
K, while one should be cautious about other environments with more drastic
temperature variations or fluctuations.

##  7 Discussion and conclusions

###  7.1 When are resolution corrections needed?

We have demonstrated that resolving the formation phase of H ii regions both
spatially and temporally is crucial to couple their momentum and energy
feedback to the ISM in a convergent manner. We have also proposed several
methods to correct these resolution problems. Here we briefly discuss in what
situations we need to include these resolution corrections.

####  7.1.1 Spatial resolution issues

First of all, we emphasise that eliminating the over-ionization and over-
heating problems is crucial to obtaining convergent simulation results. Here
we will discuss in what situations these two problems become severe numerical
issues.

Most RT-coupled Lagrangian codes with similar stellar and gas particle masses
$M_{\rm cell}\sim M_{*}\gtrsim 10^{3}\,\text{M}_{\astrosun}$ treat stars born
in a GMC as a single stellar particle (e.g. Hopkins et al., 2018; Marinacci et
al., 2019). In this case, the ionizing flux from stellar particles depends
only on the IMF-averaged mass-to-light ratio, which is $\langle q\rangle\sim
5\times 10^{46}\,\text{photons\,s}^{-1}\text{M}_{\astrosun}^{-1}$ (Rosdahl et
al., 2015). Thus, the resolution effects become important when

$n_{\text{H}}\gtrsim 16\,\text{cm}^{-3}\left(\frac{{\cal
R}_{i}}{10}\right)^{-1}\left(\frac{\langle q\rangle}{5\times
10^{46}\,\text{s}^{-1}\text{M}_{\astrosun}^{-1}}\right)\,.$ (45)

Regions with higher density in such simulations must be treated with
resolution corrections to avoid enhanced ionization feedback.

Higher-resolution simulations with $M_{\rm cell}\lesssim
10\,\text{M}_{\astrosun}$ must treat each massive star as a single star
particle (e.g. Emerick et al., 2019). In this case, the density threshold
requiring correction becomes

$n_{\text{H}}\gtrsim 32\,\text{cm}^{-3}\left(\frac{{\cal
R}_{i}}{10}\right)^{-1}\left(\frac{M_{\text{cell}}}{10\,\text{M}_{\astrosun}}\right)^{-1}Q_{48}\,.$
(46)

This result looks similar to that of low-resolution simulations. However, the
steep mass-luminosity function of individual massive stars suggests that the H
ii regions of low-mass B-type stars with low ionizing photon flux are
susceptible to the low-resolution effects. For example, if a
$7\,\text{M}_{\astrosun}$ star is still treated as an individual stellar
particle ($Q=8\times 10^{43}$ s-1, adopting the fitting function of Schaerer,
2002), the density threshold of its Strömgren sphere obtained from equation
(46) is only $0.0025$ cm-3. Therefore, choosing a reasonable mass floor to
conduct explicit RT feedback is a nontrivial problem. Fainter, less-massive
stars are also ineffective in providing momentum feedback. Estimated by the
early phase analytic results (equation 5 and 6), the total momentum injection
of a star follows $p\propto Q^{6/7}t^{3/7}$. For instance, the 7–15
$\text{M}_{\astrosun}$ stars only provide 1/10 of the momentum feedback
compared to the entire 7–120 $\text{M}_{\astrosun}$ interval (taking the mass-
luminosity and mass-lifetime functions of Schaerer (2002) weighted by the
Kroupa (2001) IMF), while the ionizing rate of a $15\,\text{M}_{\astrosun}$
can already boost this density threshold to $3$ cm-3.

Simulations focusing on the formation, accretion, and evolution of individual
stars in GMCs have much higher mass resolution and treat stars as sink
particles (e.g. Krumholz et al., 2011; Grudić et al., 2021). The typical mass
of gas cells in these simulations is $M_{\text{cell}}\sim
10^{-3}\,\text{M}_{\astrosun}$, allowing a density threshold $>10^{5}$ cm-3
following equation (46), so the spatial resolution issues may not be as
problematic for these simulations. However, these simulations usually aim to
resolve the physics in ultra-dense regions of $n>10^{7}$ cm-3, where we still
need to carefully deal with similar issues for the low-mass ionizing stars.

In summary, the initial Strömgren sphere can only be spatially resolved in
regions with relatively low density or around a particularly strong ionizing
source. Therefore, in simulations with large dynamic ranges, such as galaxy
formation, multiphase ISM, and GMC, it becomes essential to consider
implementing a correction to prevent divergent feedback from unresolved H ii
regions.

####  7.1.2 Temporal resolution issues

The requirement to ignore the temporal resolution correction is much more
stringent. As presented in Appendix B, if $\Delta t_{\text{RT}}\gg
t_{\text{rec}}$ the missing photons issue will appear as long as the spatial
resolution is high (i.e. ${\cal R}_{i}\gg 100$, equation 34) so that the
number of H atoms in a cell is much less than the number of photons. Failure
to meet this requirement is not problematic because the missing photons issue
only happens when the radiation meets neutral gas unless there are differences
between the sizes of the photon injection step and the thermochemistry step.

When the injection cadence is not the same as thermochemistry, the simplified
time-stepping criterion of photon injection step $\Delta t_{\star}\lesssim
0.1t_{\text{rec}}N_{\text{sub}}/(N_{\text{sub}}-1)$ implies a density
threshold of (combining equations 4 and 27)

$n_{\text{H}}\gtrsim
0.004\,\text{cm}^{-3}\,\eta^{-3}\left(\frac{N_{\text{sub}}}{N_{\text{sub}}-1}\right)^{3}\left(\frac{\tilde{c}}{10^{-3}c}\right)^{3}\left(\frac{{\cal
R}_{i}}{10}\right)Q_{48}^{-1}\,.$ (47)

This value is too small for most classes of simulations unless the spatial
resolution is extremely high. Thus, we strongly suggest that all RT
implementations with different sourcing and absorption cadences should
consider adopting a special time-stepping scheme or a temporal resolution
correction to avoid the missing photons issue. If the correction is not
implemented, choosing a smaller Courant factor $\eta$ or faster speed of light
$\tilde{c}$ is the easiest way to significantly increase the density threshold
($\sim\eta^{-3}\tilde{c}^{3}$) and ensure the correct simulation behaviour.

###  7.2 H2 and helium

In this work, we mainly focused on pure hydrogen gas, while the
thermochemistry of molecular hydrogen and helium can also play a significant
role in realistic astrophysical environments.

H2 in the H ii regions can be converted to H atoms through photodissociation
(PD) by 11.2–13.6 eV Lyman-Werner (LW) band photons and photoionization by
$>15.2$ eV UV photons. Theoretically, for typical massive stars, the PD-front
is merged with the I-front in the early phase, i.e. the H2 molecules are
dissociated to H atoms and ionized to H+ ions almost simultaneously (Bertoldi
& Draine, 1996). On the other hand, the processes forming H2 (dust catalysis
and gas phase formation) have very low coefficients compared to the
dissociation rates (Nickerson et al., 2018). For example, comparing the ratio
between the rate of the fastest H2 formation channel, the dust catalyst rate
$\alpha^{\text{D}}_{\text{\ce{H2}}}$, and LW band dissociation rate
$\sigma^{\text{LW}}_{\text{\ce{H2}}}$ with that of the ionization and
recombination rates of hydrogen, we have
$(\alpha^{\text{D}}_{\text{\ce{H2}}}/\sigma^{\text{LW}}_{\text{\ce{H2}}})/(\alpha_{\text{\text{H\,{ii}}}}\sigma_{\text{\text{H\,{i}}}})\approx
0.0001\,$ (assume a Milky-way dust-to-gas ratio of $\sim 0.01$, Draine et al.,
2007). Thus, both spatially and temporally resolving the dissociation of H2
when the H ii region is forming is much easier than doing this for H atoms.

For helium, the case B radiative recombination rate for He ii is similar to H
i ($2.72\times 10^{-13}T_{4}^{-0.789}$ cm3 s-1, Benjamin et al., 1999).
Adopting a primordial abundance of $X=0.76$, the number density ratio of He to
H is 0.08, so the initial Strömgren mass of helium is about
$3.5Q_{\text{He}}/Q_{\text{H}}$ times as much as that of hydrogen, where
$Q_{\text{H}}$ and $Q_{\text{He}}$ are the ionization photon rates of H and
He, respectively (equation 64). Therefore, for a given star, the mass
resolution of its helium Strömgren sphere is $15Q_{\text{He}}/Q_{\text{H}}$
times that of hydrogen. Resolving the initial Strömgren mass of helium is even
easier for most O-type stars with $Q_{\rm He}>Q_{\text{H}}/15$ so long as
their $T_{\text{eff}}\gtrsim 3.5\times 10^{4}$ K (Sternberg et al., 2003;
Martins et al., 2005), so the mass resolution for helium can meet this
requirement in most situations as long as hydrogen ionization is spatially
resolved.

The ionization of He ii requires photons with $h\nu>54.4$ eV.
$Q_{\text{\text{He\,{ii}}}}$ is usually very small for ordinary massive stars
and the recombination rate for He iii is about 10 times larger ($2.19\times
10^{-12}\,\text{cm}^{-3}$ s-1 Verner & Ferland, 1996). This makes the He iii
Strömgren sphere with very little mass difficult to be resolved. Additionally,
resolving or correcting the He ii ionization needs to obtain the correct He ii
fraction at first, which is not always guaranteed. Thus, the mass and physical
states of He iii content might be highly uncertain in most simulations but
overall He iii has a negligible contribution to the feedback.

For the temporal resolution of helium, the characteristic time-scale should be
similar to that of hydrogen, the recombination time-scale
$\\{t_{\text{rec,\text{He\,{ii}}}},t_{\text{rec,\text{He\,{iii}}}}\\}=\\{1/\alpha_{\text{\text{He\,{ii}}}}(n_{\text{H}}+n_{\text{He}}),1/\alpha_{\text{\text{He\,{iii}}}}(n_{\text{H}}+2n_{\text{He}})\\}$.
As the recombination coefficient ($\alpha_{\text{\text{He\,{ii}}}}$) of He ii
is comparable to that of H, $t_{\text{rec,\text{He\,{ii}}}}$ is also
comparable to that of H, making it similarly difficult to be temporally
resolved. However, for He iii, $\alpha_{\text{\text{He\,{iii}}}}$ is one order
of magnitude larger, making $t_{\text{rec,\text{He\,{iii}}}}$ one order of
magnitude smaller. This, in turn, makes the He ii–He iii thermochemistry more
challenging to be temporally resolved.

In summary, resolving the hydrogen Strömgren sphere presents the greatest
numerical challenge and dynamical significance compared to H2 and He i. As
such, the resolution requirements and correction methods proposed in this
paper remain relevant even in more complex chemical networks. For primordial
chemical networks, we outline the implementation of our spatial resolution
correction scheme, which includes He and H2 chemistry, in Appendix D.
Additionally, readers can refer to Jaura et al. (2020) for a limited
absorption correction that also considers He and H2.

###  7.3 H ii regions in inhomogeneous density structures

Realistic MCs exhibit highly filamentary structures, where protostars form in
overdensities characterized by power-law density profiles ($\rho\propto
r^{-w}$, Larson 1981; Kauffmann et al. 2010, see André et al. 2014 for a
review). Franco et al. (1990) demonstrated that H ii regions in media with
$w>3/2$ can expand in an accelerating manner, driving a “champagne flow". In
complex structures, H ii regions tend to expand toward the rarefied regions
and stall in the dense regions as blister-type regions (Mellema et al., 2006;
Gendelev & Krumholz, 2012; Zamora-Avilés et al., 2019; Jin et al., 2022).

In numerical simulations, these effects can be self-consistently modelled on a
scale larger than the numerical resolution, but the small-scale density
gradients are smoothed due to the inadequate resolution. In the absence of a
sub-grid turbulence model, gas represented by a cell is assumed to be uniform.
Consequently, numerically unresolved H ii regions suffer from over-ionization
and over-heating issues similar to those in the uniform medium, while our
corrections described in Section 6.1.2 ensure the correct ionized mass
predicted by the uniform solutions. However, the Strömgren mass evaluated
based on the mean density may overestimate the ionized mass since stars are
embedded in unresolved overdense cores where the H ii regions are trapped by
the high density. This potentially leads to an overestimation of feedback from
these trapped H ii regions.

Nonetheless, once the H ii regions break through these dense cores due to
expansion or the assistance of stellar winds and radiation pressure (e.g. Geen
et al., 2020, 2021), they are easier to resolve (${\cal R}_{i}\propto n^{-1}$)
and the local density fluctuations will have a less significant impact on the
global properties of the H ii regions as they expand toward the low-density
gas. High-resolution radiation-hydrodynamic simulations by Mellema et al.
(2006) showed that though the expansion of H ii region in turbulent MCs is
slower than that in uniform medium with the same mean density due to the
overdensity around the star, over time, the mean radius can eventually catch
up with the uniform case. Hence, our corrections are considered reasonable
approximations to capture the feedback from unresolved H ii regions in
realistic environments.

###  7.4 Other early feedback channels and their combination

In addition to ionization feedback, luminous massive stars also provide
energetic early (pre-SN) feedback through radiation pressure and high-velocity
stellar winds. However, resolving the effects of radiation pressure in
numerical simulations presents its own challenges due to the high spatial and
temporal resolution requirements, especially in dense dusty regions (e.g.
Krumholz, 2018; Hopkins & Grudić, 2019). Similarly, modelling the effects of
stellar winds also requires high resolution simulations (e.g. Pittard et al.,
2021). The rapid energy loss due to the turbulent mixing at the wind
bubble/shell interface poses additional challenges that need to be addressed
to accurately simulate the wind feedback (Lancaster et al., 2021a, b). The
combination of all these early feedback channels in simulations is nonlinear
and can provide a complicated collective impact rather than their simple
superposition (e.g. Haid et al., 2018; Geen et al., 2021). Additional
challenges on the numerical resolution can be introduced by such physics and
requires further investigation to better understand their collective impact on
galaxy formation simulations.

###  7.5 Summary

We have performed a suite of radiation hydrodynamic simulations of idealized H
ii regions with the moment-based M1 closure RT solver in the moving-mesh code
arepo. Below, we list our main findings.

  1. 1.

Sufficient mass (spatial) resolution (${\cal
R}_{i}=M_{\text{S}}/M_{\text{cell}}$) is critical to accurately capture the
evolution and feedback of simulated H ii regions. With a spatial resolution
higher than $\mathcal{R}_{i}=10$ (more than 10 cells inside the initial
Strömgren sphere), the momentum and energy feedback from an expanding H ii
region can be reproduced with an acceptable numerical error. If we want to
reproduce the profiles of shock parameters, the resolution should be higher
than $\mathcal{R}_{i}=100$.

  2. 2.

In the formation phase of H ii regions, insufficient spatial resolution lowers
the ionization fraction of gas cells but can lead to an overestimation of the
ionized mass when the H ii region is forming, i.e. the over-ionization
problem. In the expansion phase, the insufficient spatial resolution will fail
to heat the ionized gas cells to a correct temperature and result in
insufficient cooling, the total thermal energy of gas will thus also be
overestimated and lead to enhanced momentum feedback, i.e. the over-heating
problem. It is crucial to avoid over-ionization and over-heating problems in
simulations because they both lead to divergent feedback when the H ii regions
are unresolved. To correct the spatial resolution problems, we can consider
both lowering the number of neighbours $N_{\text{nb}}$, as well as introducing
a correction to reestablish the accurate balance between ionization and
recombination. However, although these corrections can correct the momentum
feedback from H ii regions, they may not help enough for dynamic multi-phase
gas structures. The over-heating problem will turn the highly-ionized, hot
($>10^{4}$ K) gas concentrated in the H ii regions into diffuse, partially
ionized, warm ($\sim 8000$ K) gas, changing the multi-phase gas structure of
the ISM in star-forming regions.

  3. 3.

If the size of the time step is too large, almost all photons dumped in the
first step will be lost and there is no mechanism to compensate for these
missing photons, leading to the missing photons issue. Most simulations of
real astrophysical systems (e.g., galaxy formation, multi-phase ISM, GMC)
cannot temporally resolve the ionization timescale of the gas. Thus, the
missing photons issue should be treated carefully if the cadences of photon
injection and thermochemistry are different in these simulations. A rough
time-stepping criterion for the injection time step is $\Delta
t_{\star}\lesssim 0.1t_{\text{rec}}N_{\text{sub}}/(N_{\text{sub}}-1)$ when
each injection time step is associated with $N_{\text{sub}}$ RT steps.

  4. 4.

If the ratio of injection step size to RT step size
($t_{\star}/t_{\text{RT}}$) is too large, the missing photons issue will
repeat after each injection. Ionizing fluxes reduced in this way will be
insufficient to ionize and heat the gas and reduce the feedback from the H ii
regions. Using a thermochemistry solver coupled with the photon densities can
obtain an accurate solution and solve the missing photons problem
fundamentally, but this is an expensive choice when including many species and
frequency bins. Alternatively, the uncoupled solver with the limited photon
absorption approximation proposed by Jaura et al. (2020) is a much cheaper but
very powerful correction for this issue.

  5. 5.

Spatially and temporally resolving the thermochemistry of hydrogen is
generally the most challenging part compared to that of H2 and He i. While
resolving and correcting the He ii–He iii thermochemistry also presents its
own difficulties, it is of relatively lesser importance for stellar feedback.
Once the ionization feedback from hydrogen is spatial and temporally resolved,
those for H2, He i will be already solved, alleviating the complexity when
designing resolved ionization feedback models. We also outline a method to
implement our spatial resolution correction scheme including He and H2
chemistry in Appendix D, and we refer the readers to Jaura et al. (2020) for
the limited absorption correction including He and H2.

## Acknowledgements

We thank Volker Springel for giving us access to arepo. YD is grateful to Josh
Borrow and Yang Ni for useful discussions. YD was a visiting student at the
Massachusetts Institute of Technology sponsored by the ZhengGang Fund of NJU
through the ZhengGang Scholarship for Overseas Study. AS acknowledges support
under an Institute for Theory and Computation Fellowships at the Center for
Astrophysics $|$ Harvard & Smithsonian. MV acknowledges support through NASA
ATP 19-ATP19-0019, 19-ATP19-0020, 19-ATP19-0167, and NSF grants AST-1814053,
AST-1814259, AST-1909831, AST-2007355 and AST-2107724. GLB acknowledges
support from the NSF (AST-2108470, XSEDE grant MCA06N030), NASA TCAN award
80NSSC21K1053, and the Simons Foundation through their support of the Learning
the Universe collaboration. The simulations of this work were run on the MIT
Engaging cluster, the Anvil cluster at Purdue University as part of XSEDE
through TG-PHY220025, and the Stampede2 HPC resource at Texas Advanced
Computing Center as part of XSEDE through TG-AST200007 and TG-MCA06N030. We
use python packages NumPy (Harris et al., 2020), SciPy (Virtanen et al.,
2020), astropy (Astropy Collaboration et al., 2013, 2018), and matplotlib
(Hunter, 2007) to analyze and visualize the simulation data.

## data availability

The data that support the findings of this study are available from the
corresponding author, upon reasonable request.

## References

  * Agertz et al. (2013) Agertz O., Kravtsov A. V., Leitner S. N., Gnedin N. Y., 2013, ApJ, 770, 25
  * Agertz et al. (2020) Agertz O., et al., 2020, MNRAS, 491, 1656
  * André et al. (2014) André P., Di Francesco J., Ward-Thompson D., Inutsuka S. I., Pudritz R. E., Pineda J. E., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. pp 27–51 (arXiv:1312.6232), doi:10.2458/azu_uapress_9780816531240-ch002
  * Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, A&A, 558, A33
  * Astropy Collaboration et al. (2018) Astropy Collaboration et al., 2018, AJ, 156, 123
  * Aubert & Teyssier (2008) Aubert D., Teyssier R., 2008, MNRAS, 387, 295
  * Baker & Menzel (1938) Baker J. G., Menzel D. H., 1938, ApJ, 88, 52
  * Benjamin et al. (1999) Benjamin R. A., Skillman E. D., Smits D. P., 1999, ApJ, 514, 307
  * Bertoldi & Draine (1996) Bertoldi F., Draine B. T., 1996, ApJ, 458, 222
  * Bisbas et al. (2015) Bisbas T. G., et al., 2015, MNRAS, 453, 1324
  * Cen (1992) Cen R., 1992, ApJS, 78, 341
  * Chan et al. (2021) Chan T. K., Theuns T., Bower R., Frenk C., 2021, MNRAS, 505, 5784
  * Chevance et al. (2022) Chevance M., et al., 2022, MNRAS, 509, 272
  * Courant et al. (1928) Courant R., Friedrichs K., Lewy H., 1928, Mathematische Annalen, 100, 32
  * Dale et al. (2012) Dale J. E., Ercolano B., Bonnell I. A., 2012, MNRAS, 424, 377
  * Dale et al. (2014) Dale J. E., Ngoumou J., Ercolano B., Bonnell I. A., 2014, MNRAS, 442, 694
  * Draine et al. (2007) Draine B. T., et al., 2007, ApJ, 663, 866
  * Emerick et al. (2018) Emerick A., Bryan G. L., Mac Low M.-M., 2018, ApJ, 865, L22
  * Emerick et al. (2019) Emerick A., Bryan G. L., Mac Low M.-M., 2019, MNRAS, 482, 1304
  * Fall et al. (2010) Fall S. M., Krumholz M. R., Matzner C. D., 2010, ApJ, 710, L142
  * Franco et al. (1990) Franco J., Tenorio-Tagle G., Bodenheimer P., 1990, ApJ, 349, 126
  * Garaldi et al. (2022) Garaldi E., Kannan R., Smith A., Springel V., Pakmor R., Vogelsberger M., Hernquist L., 2022, MNRAS, 512, 4909
  * Geen et al. (2015) Geen S., Rosdahl J., Blaizot J., Devriendt J., Slyz A., 2015, MNRAS, 448, 3248
  * Geen et al. (2020) Geen S., Pellegrini E., Bieri R., Klessen R., 2020, MNRAS, 492, 915
  * Geen et al. (2021) Geen S., Bieri R., Rosdahl J., de Koter A., 2021, MNRAS, 501, 1352
  * Gendelev & Krumholz (2012) Gendelev L., Krumholz M. R., 2012, ApJ, 745, 158
  * Gnedin & Abel (2001) Gnedin N. Y., Abel T., 2001, New Astron., 6, 437
  * Grudić et al. (2021) Grudić M. Y., Guszejnov D., Hopkins P. F., Offner S. S. R., Faucher-Giguère C.-A., 2021, MNRAS, 506, 2199
  * Gutcke et al. (2021) Gutcke T. A., Pakmor R., Naab T., Springel V., 2021, MNRAS, 501, 5597
  * Haid et al. (2016) Haid S., Walch S., Naab T., Seifried D., Mackey J., Gatto A., 2016, MNRAS, 460, 2962
  * Haid et al. (2018) Haid S., Walch S., Seifried D., Wünsch R., Dinnbier F., Naab T., 2018, MNRAS, 478, 4799
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Harten et al. (1983) Harten A., Lax P. D., Leer B. v., 1983, SIAM review, 25, 35
  * Hernández-Aguayo et al. (2022) Hernández-Aguayo C., et al., 2022, arXiv e-prints, p. arXiv:2210.10059
  * Hindmarsh et al. (2005) Hindmarsh A. C., Brown P. N., Grant K. E., Lee S. L., Serban R., Shumaker D. E., Woodward C. S., 2005, ACM Transactions on Mathematical Software (TOMS), 31, 363
  * Hopkins & Grudić (2019) Hopkins P. F., Grudić M. Y., 2019, MNRAS, 483, 4187
  * Hopkins et al. (2012) Hopkins P. F., Quataert E., Murray N., 2012, MNRAS, 421, 3488
  * Hopkins et al. (2018) Hopkins P. F., et al., 2018, MNRAS, 480, 800
  * Hopkins et al. (2020) Hopkins P. F., Grudić M. Y., Wetzel A., Kereš D., Faucher-Giguère C.-A., Ma X., Murray N., Butcher N., 2020, MNRAS, 491, 3702
  * Hosokawa & Inutsuka (2006) Hosokawa T., Inutsuka S.-i., 2006, ApJ, 646, 240
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Ivkovic (2023) Ivkovic M., 2023, arXiv e-prints, p. arXiv:2302.12727
  * Jaura et al. (2018) Jaura O., Glover S. C. O., Klessen R. S., Paardekooper J. P., 2018, MNRAS, 475, 2822
  * Jaura et al. (2020) Jaura O., Magg M., Glover S. C. O., Klessen R. S., 2020, MNRAS, 499, 3594
  * Jeffreson et al. (2021) Jeffreson S. M. R., Krumholz M. R., Fujimoto Y., Armillotta L., Keller B. W., Chevance M., Kruijssen J. M. D., 2021, MNRAS, 505, 3470
  * Jin et al. (2022) Jin Y., Kewley L. J., Sutherland R. S., 2022, ApJ, 934, L8
  * Kannan et al. (2014) Kannan R., et al., 2014, MNRAS, 437, 2882
  * Kannan et al. (2019) Kannan R., Vogelsberger M., Marinacci F., McKinnon R., Pakmor R., Springel V., 2019, MNRAS, 485, 117
  * Kannan et al. (2020a) Kannan R., Marinacci F., Simpson C. M., Glover S. C. O., Hernquist L., 2020a, MNRAS, 491, 2088
  * Kannan et al. (2020b) Kannan R., Marinacci F., Vogelsberger M., Sales L. V., Torrey P., Springel V., Hernquist L., 2020b, MNRAS, 499, 5732
  * Kannan et al. (2022a) Kannan R., Garaldi E., Smith A., Pakmor R., Springel V., Vogelsberger M., Hernquist L., 2022a, MNRAS, 511, 4005
  * Kannan et al. (2022b) Kannan R., Smith A., Garaldi E., Shen X., Vogelsberger M., Pakmor R., Springel V., Hernquist L., 2022b, MNRAS, 514, 3857
  * Kannan et al. (2023) Kannan R., et al., 2023, MNRAS, 524, 2594
  * Katz (1992) Katz N., 1992, ApJ, 391, 502
  * Katz et al. (1996) Katz N., Weinberg D. H., Hernquist L., 1996, ApJS, 105, 19
  * Kauffmann et al. (2010) Kauffmann J., Pillai T., Shetty R., Myers P. C., Goodman A. A., 2010, ApJ, 716, 433
  * Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
  * Kruijssen et al. (2019) Kruijssen J. M. D., et al., 2019, Nature, 569, 519
  * Krumholz (2018) Krumholz M. R., 2018, MNRAS, 480, 3468
  * Krumholz et al. (2011) Krumholz M. R., Klein R. I., McKee C. F., 2011, ApJ, 740, 74
  * Lahén et al. (2020) Lahén N., Naab T., Johansson P. H., Elmegreen B., Hu C.-Y., Walch S., Steinwandel U. P., Moster B. P., 2020, ApJ, 891, 2
  * Lancaster et al. (2021a) Lancaster L., Ostriker E. C., Kim J.-G., Kim C.-G., 2021a, ApJ, 914, 89
  * Lancaster et al. (2021b) Lancaster L., Ostriker E. C., Kim J.-G., Kim C.-G., 2021b, ApJ, 914, 90
  * Larson (1981) Larson R. B., 1981, MNRAS, 194, 809
  * Lee & Hennebelle (2018) Lee Y.-N., Hennebelle P., 2018, A&A, 611, A88
  * Levermore (1984) Levermore C. D., 1984, J. Quant. Spectrosc. Radiative Transfer, 31, 149
  * Li et al. (2019) Li H., Vogelsberger M., Marinacci F., Gnedin O. Y., 2019, MNRAS, 487, 364
  * Marinacci et al. (2019) Marinacci F., Sales L. V., Vogelsberger M., Torrey P., Springel V., 2019, MNRAS, 489, 4233
  * Martins et al. (2005) Martins F., Schaerer D., Hillier D. J., 2005, A&A, 436, 1049
  * McKee & Tan (2003) McKee C. F., Tan J. C., 2003, ApJ, 585, 850
  * Mellema et al. (2006) Mellema G., Arthur S. J., Henney W. J., Iliev I. T., Shapiro P. R., 2006, ApJ, 647, 397
  * Moscadelli et al. (2021) Moscadelli L., Cesaroni R., Beltrán M. T., Rivilla V. M., 2021, A&A, 650, A142
  * Naab & Ostriker (2017) Naab T., Ostriker J. P., 2017, ARA&A, 55, 59
  * Nickerson et al. (2018) Nickerson S., Teyssier R., Rosdahl J., 2018, MNRAS, 479, 3206
  * Osterbrock & Ferland (2006) Osterbrock D. E., Ferland G. J., 2006, Astrophysics of gaseous nebulae and active galactic nuclei
  * Pakmor et al. (2016) Pakmor R., Springel V., Bauer A., Mocz P., Munoz D. J., Ohlmann S. T., Schaal K., Zhu C., 2016, MNRAS, 455, 1134
  * Peter et al. (2023) Peter T., Klessen R. S., Kanschat G., Glover S. C. O., Bastian P., 2023, MNRAS, 519, 4263
  * Petkova & Springel (2009) Petkova M., Springel V., 2009, MNRAS, 396, 1383
  * Petkova et al. (2021) Petkova M. A., Vandenbroucke B., Bonnell I. A., Kruijssen J. M. D., 2021, MNRAS, 507, 858
  * Pillepich et al. (2018) Pillepich A., et al., 2018, MNRAS, 473, 4077
  * Pittard et al. (2021) Pittard J. M., Wareing C. J., Kupilas M. M., 2021, MNRAS, 508, 1768
  * Pittard et al. (2022) Pittard J. M., Kupilas M. M., Wareing C. J., 2022, MNRAS, 510, 2797
  * Raga et al. (2012) Raga A. C., Cantó J., Rodríguez L. F., 2012, MNRAS, 419, L39
  * Raičević et al. (2014) Raičević M., Pawlik A. H., Schaye J., Rahmati A., 2014, MNRAS, 437, 2816
  * Raskutti et al. (2016) Raskutti S., Ostriker E. C., Skinner M. A., 2016, ApJ, 829, 130
  * Rosdahl et al. (2013) Rosdahl J., Blaizot J., Aubert D., Stranex T., Teyssier R., 2013, MNRAS, 436, 2188
  * Rosdahl et al. (2015) Rosdahl J., Schaye J., Teyssier R., Agertz O., 2015, MNRAS, 451, 34
  * Rosdahl et al. (2022) Rosdahl J., et al., 2022, MNRAS, 515, 2386
  * Rusanov (1961) Rusanov V., 1961, J. Comput. Math. Phys. USSR, 1, 267
  * Schaerer (2002) Schaerer D., 2002, A&A, 382, 28
  * Schaye et al. (2015) Schaye J., et al., 2015, MNRAS, 446, 521
  * Shapiro (1986) Shapiro P. R., 1986, PASP, 98, 1014
  * Shu (1977) Shu F. H., 1977, ApJ, 214, 488
  * Smith et al. (2020) Smith A., Kannan R., Tsang B. T. H., Vogelsberger M., Pakmor R., 2020, ApJ, 905, 27
  * Smith et al. (2022a) Smith A., Kannan R., Garaldi E., Vogelsberger M., Pakmor R., Springel V., Hernquist L., 2022a, MNRAS, 512, 3243
  * Smith et al. (2022b) Smith A., et al., 2022b, MNRAS, 517, 1
  * Spitzer (1978) Spitzer L., 1978, Physical processes in the interstellar medium, doi:10.1002/9783527617722.
  * Springel (2010) Springel V., 2010, MNRAS, 401, 791
  * Sternberg et al. (2003) Sternberg A., Hoffmann T. L., Pauldrach A. W. A., 2003, ApJ, 599, 1333
  * Strömgren (1939) Strömgren B., 1939, ApJ, 89, 526
  * Tacchella et al. (2022) Tacchella S., et al., 2022, MNRAS, 513, 2904
  * Verner & Ferland (1996) Verner D. A., Ferland G. J., 1996, ApJS, 103, 467
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Vogelsberger et al. (2020) Vogelsberger M., Marinacci F., Torrey P., Puchwein E., 2020, Nature Reviews Physics, 2, 42
  * Walch & Naab (2015) Walch S., Naab T., 2015, MNRAS, 451, 2757
  * Whalen et al. (2004) Whalen D., Abel T., Norman M. L., 2004, ApJ, 610, 14
  * Wise & Abel (2011) Wise J. H., Abel T., 2011, MNRAS, 414, 3458
  * Yorke (1986) Yorke H. W., 1986, ARA&A, 24, 49
  * Zamora-Avilés et al. (2019) Zamora-Avilés M., et al., 2019, MNRAS, 487, 2200

##  Appendix A Cooling Rates

In this section, we list the fitting formulas accounting for the cooling
processes in pure hydrogen H ii regions used in this work (Cen, 1992).

The total cooling rate can be split as follows

$\Lambda=\Lambda_{\text{rec}}+\Lambda_{\text{bb}}+\Lambda_{\text{bf}}+\Lambda_{\text{ff}}\,,$
(48)

where $\Lambda_{\text{rec}}$, $\Lambda_{\text{bb}}$, $\Lambda_{\text{bf}}$,
and $\Lambda_{\text{ff}}$ are the cooling rates accounting for recombination,
bound-bound collision, bound-free, collision, and free-free collision
(bremsstrahlung), respectively We adopt the following approximation results

$\Lambda_{\text{rec}}=8.7\times
10^{-27}n_{e}n(\text{\text{H\,{ii}}})T^{1/2}\frac{T_{3}^{-0.2}}{1+T_{6}^{0.7}},$
(49)

$\Lambda_{\mathrm{ff}}=1.42\times
10^{-27}T^{1/2}g_{\mathrm{ff}}n_{e}n_{\text{\text{H\,{ii}}}}.$ (50)

$\Lambda_{\text{bb}}=7.5\times
10^{-19}n_{e}n(\text{\text{H\,{i}}})\frac{e^{-\frac{\epsilon_{21,1}}{k_{\text{B}}T}}}{1+T_{5}^{1/2}},$
(51) $\Lambda_{\text{bf}}=1.27\times
10^{-21}n_{e}n(\text{\text{H\,{i}}})T^{1/2}\frac{e^{-\frac{I_{1}}{k_{\text{B}}T}}}{1+T_{5}^{1/2}}.$
(52)

where $k_{B}$ is the Boltzmann constant, $I_{1}=13.6$ eV is the ionization
potential of H, $\epsilon_{21}=10.2$ eV is the energy level difference between
the ground state and first excited state of the hydrogen atom.

##  Appendix B Step-by-step tracking of the missing photons issue

Figure 15: Cartoon illustrating the step-by-step photon propagation and
thermochemistry after a massive star ignites. Far left: The 1st step
($t=0\rightarrow 1\Delta t$), when the star ignites and dumps $Q\Delta t$
photons into its $N_{\text{nb}}$ neighbouring cells. The neighbouring cells
are all fully ionized and the missing photons issue happens for the first
time. Middle left: The 2nd step ($t=1\Delta t\rightarrow 2\Delta t$), when the
photons stay in these neighbouring cells and the star dumps $Q\Delta t$
photons again. The photon density in these cells begins to increase linearly.
Middle right: The 3rd step ($t=2\Delta t\rightarrow 3\Delta t$), when the
photons dumped in the first step will propagate outwards into the cells
enclosing the $N_{\text{nb}}$ neighbouring cells but do not lead to ionization
because most of them are absorbed as missing photons. Far right: The 4th step
($t=3\Delta t\rightarrow 4\Delta t$), when the photons dumped in the 2nd step
now propagate into the cells enclosing the nearest neighbouring cells and the
second missing photons issue happens.

In this section, we demonstrate the missing photons issue by step-by-step
tracking the evolution of photon density immediately after the star ignition
with the aid of Figure 15.

Like previous tests, a stellar particle with ionizing photon rate $Q$ sits at
the centre of a pure hydrogen uniform simulation box which we take to be a
regular Cartesian mesh here for convenience. Once the star ignites, $Q\Delta
t$ photons are dumped into its $N_{\text{nb}}$ neighbouring cells
isotropically at the beginning of each step and then propagate outward. After
photon dumping and propagating, the code will perform the thermochemistry and
photon absorption by solving equation (10) and equation (29). We simply assume
the photon dumping, propagating, and thermochemistry is conducted successively
during one integration step, meaning that the performance of these steps
shares the same cadence.

The 1st step ($t=0\rightarrow 1\Delta t$): Once the star ignites,
$N_{\gamma}=Q\Delta t$ photons are dumped into the neighbouring cells during
the first $\Delta t$. As described in Section 3.2.1, these photons are
allocated to $N_{\rm nb}$ cells in the neighbourhood list found by the
neighbour search. The other cells are physically impossible to receive any
photons because $\Delta t$ must be smaller than the Courant criterion $\Delta
t_{\rm C}$ so that no photons can escape from these neighbouring cells
(assuming the neighbouring cells can completely wrap the centre star). Thus,
the number of photons in each neighbouring cell (at $t=\Delta t$) is

$N_{\gamma,1}=\frac{Q}{N_{\rm nb}}\Delta t\,,$ (53)

where the subscript of $N_{\gamma}$ accounts for the step.

The number of photons consumed in this step is (equation 30)

$\Delta N_{\gamma,1}=\frac{Q}{N_{\rm nb}}\Delta t(1-e^{-\Delta t/t_{\rm
ion}})\approx\frac{Q}{N_{\rm nb}}\Delta t\,,$ (54)

which exhibits that almost all photons dumped into the gas will be consumed
when the $t_{\text{ion}}$ (or MFP) is unresolved. Again, we emphasize that the
neutral density $n_{\rm HI}$ involved here is taken as the value at the
beginning of the step (i.e. $1\cdot n_{\rm H}$) because the photon absorption
equation is not coupled with the thermochemistry equations.

On the other hand, the total number of neutral hydrogen atoms in a cell is

$N_{\text{\text{H\,{i}}}}=\frac{M_{\rm cell}}{m_{\rm H}}=\frac{Q}{\alpha_{\rm
B}n_{\rm H}{\cal R}_{i}}\,.$ (55)

Assuming one photon can ionize one neutral hydrogen atom (MFP unresolved), in
this case, if $\Delta N_{\gamma}>(1+f_{\text{rec}})N_{\text{\text{H\,{i}}}}$,
then $\Delta N_{\gamma}$ photons will be absorbed but only
$N_{\text{\text{H\,{i}}}}$ atoms can be ionized. Here $f_{\text{rec}}$ is a
factor accounting for the recombination of ions. In other words, $\delta
N_{\gamma}$ photons simply disappear without any effect on the gas state. The
photon number is no longer conserved. For each neighbouring cell, the number
of missing photons is

$\delta N_{\gamma,1}=\frac{Q}{N_{\rm nb}}\Delta
t-(1+f_{\text{rec}})\frac{Q}{\alpha_{\rm B}n_{\rm H}{\cal R}_{i}}\,,$ (56)

where $f_{\text{rec}}$ can be estimated as $\Delta t/t_{\text{rec}}$ because
$\Delta t\sim t_{\text{rec}}\gg t_{\text{ion}}$.

On the contrary, if the MFP is still unresolved but $\Delta
N_{\gamma}<N_{\text{\text{H\,{i}}}}$, the ionization will be overestimated in
this step. Nonetheless, this error can be relieved in the subsequent steps by
the shifting of ionization and recombination rates. Thus, it has minor effects
on the results compared to the missing photons issue.

The 2nd step ($t=1\Delta t\rightarrow 2\Delta t$): If the missing photons
issue happened in the first step, the neighbouring cells must be fully ionized
in this step ($\Delta N_{\gamma}>N_{\text{\text{H\,{i}}}}$). With $\Delta
t<\Delta t_{\rm C}$, the photons will stay in these neighbouring cells during
the next few steps as long as $t<\Delta x/\tilde{c}=\Delta t_{\rm C}/\eta$,
where $\eta$ is the Courant factor. Again, the change of photon density due to
thermochemistry is calculated by equation (29) with a zero $n_{\rm HI}$
because the cell is fully ionized. Thus, the photon number density will
increase linearly without any missing photons issues in these steps.

Here we define $M=\lceil\Delta t_{\rm C}/\eta\Delta t\rceil$, the period of
this linear increasing stage is thus $M\Delta t$. After the $M$th step, the
photons dumped in the first step will propagate outwards into the cells
enclosing the $N_{\rm nb}$ neighbouring cells. For convenience, we simply
assume $M=3$ (i.e. $2<\Delta t_{\rm C}/\eta\Delta t<3$).

The 3rd step ($t=2\Delta t\rightarrow 3\Delta t$): In the third step, the
photons dumped in the first step will propagate outwards into the cells
enclosing the $N_{\rm nb}$ neighbouring cells. However, the photons dumped in
the first step have all been absorbed and become missing photons. Thus, what