\Delta t}, 
  \label{eq:bdf}$$ where all source terms are evaluated at the advanced timestep. However, not all source terms can be evaluated at the advanced timestep due to the intrinsic nonlinearity of Eqn. (\[eq:chemistry\]). We hence sequentially update densities of all species in the order of increasing ionization states rather than updating them simultaneously; We evaluate the source terms contributed by the ionization from and recombination to the lower states at the advanced timesteps. This method has been found to be very efficient and accurate [e.g. @anninos97; @yoshikawa06].

Further improvements in accuracy and stability can be made by subcycling the rate solver over a single timestep with which the RT is solved. The subcycle timestep, which we call the ‘chemical timestep’, is determined so that the maximum fractional change in the electron density is limited to 10% per timestep: $$\Delta t_{\rm chem} = 0.1 \left|\frac{n_{\rm e}}{\dot{n}_{\rm e}}\right|. 
\label{eq:chemicaltimespte}$$

Photo-heating and radiative cooling {#sec:heatcool}
-----------------------------------

Similarly to the photoionization, photo-heating rate for each mesh due to the photoionization of the $i$-th species is given by $${\cal H}_i = \sum_\alpha {\cal H}_{i, \alpha},$$ where ${\cal H}_{i, \alpha}$ indicates the contribution from a radiation source (or a group of sources), $\alpha$, and $i =  