the timestep from becoming too short when the medium is almost neutral. Our choice for $\epsilon_{\rm e}$ and $\epsilon_{\rm HI}$ are 0.2 and 0.002, respectively. We follow the evolution of the system with the minimum of the timestep defined by Eqn. (\[eq:timestep\_i\]), i.e. $$\Delta t = \Delta t_{i, {\rm min}}. 
\label{eq:timestep}$$ The timestep $\Delta t$ is therefore only about twice as long as the shortest chemical timestep, $\Delta t_{{\rm chem}, {\rm min}}$. With this timestep, we find the solutions typically within 3 to 6 iteration steps. While we can of course use a longer timestep, a longer timestep requires more iterations and the total number of steps becomes similar or even larger than the case we employ the timestep defined by Eqn. (\[eq:timestep\]). With a longer timestep, the solutions sometimes never converge. This timestep is in general much shorter than the timestep defined by the Courant timestep criterion and therefore we have to subcycle the hydrodynamical timestep with this timestep when we couple the RT with the hydrodynamics.

When optically thick meshes exist, the solutions converge very slowly. We thus use smoothed photoionization rates, $\tilde{\Gamma}_i$, and heating rates, $\tilde{\cal H}_i$, instead of $\Gamma_i$ and ${\cal H}_i$. The smoothed rates for 