up-to-date as possible. The sources of these rates are summarised in Table \[table:rates\]. Note that there are notable differences in the recombination cooling rates between literatures [see @iliev06].

Time stepping
-------------

Since the optical depth $\tau(\nu)$ at $t + \Delta t$ depends on densities of all species at $t + \Delta t$, we have to solve the static RT equation (Eqn. (\[eq:formal2\])), the chemical reactions (Eqn. (\[eq:bdf\])), and the energy equation (Eqn. (\[eq:heatcool\])) iteratively. We iterate these steps until the relative difference in the electron number density becomes sufficiently small: $|n_{\rm e}^{(n)} - n_{\rm e}^{(n - 1)}|/n_{\rm e}^{(n)} < \epsilon$, where superscripts indicate the number of iterations and we set $\epsilon$ to $10^{-4}$ throughout this paper. The timestep $\Delta t$, with which we solve the RT equation to obtain $\Gamma_i^{t + \Delta t}$ and ${\cal H}_i^{t + \Delta t}$, could be much larger than the chemical timestep $\Delta t_{\rm chem}$, with which we subcycle the rate and energy solvers.

We however choose to employ a timestep that is defined by the timescale of the chemical reactions: $$\Delta t_i = \epsilon_{\rm e} \left|\frac{n_{\rm e}}{\dot{n}_{\rm e}}\right|_i +  
\epsilon_{\rm HI} \left|\frac{n_{\rm HI}}{\dot{n}_{\rm HI}}\right|_i, 
  \label{eq:timestep_i}$$ where the second term in the right-hand side prevent 