we update the densities of each chemical species implicitly by using a backward difference formula (BDF). The equations to evolve the density of each species can be generally written as $$\frac{{{\rm d}}n_i}{{{\rm d}}t} = C_i (T, n_j) - D_i (T, n_j) n_i, 
\label{eq:chemistry}$$ where $n_i = \rho_i / (A_i m_{\rm H})$, $A_i$ is the atomic mass number of the $i$-th species, and $m_{\rm H}$ is the proton mass. This time $i$ is either e, ${\ensuremath{\mbox{{H$\,${\sc i}}}}}$, ${\mbox{{H$\,${\sc ii}}}}$, ${\mbox{{He$\,${\sc i}}}}$, ${\mbox{{He$\,${\sc ii}}}}$ or ${\mbox{{He$\,${\sc iii}}}}$. The first term of the right-hand side, $C_i$, is the collective source term responsible for the creation of the $i$-th species. The second term involving $D_i$ represents the destruction mechanisms for the $i$-th species and are thus proportional to $n_i$.

Since the timescales for the ionization and recombination differ by many orders of magnitude depending on chemical species, Eqn. (\[eq:chemistry\]) is a stiff set of differential equations. In numerically solving a stiff set of equations, implicit schemes are required unless an unreasonably small timestep is employed. As in @anninos97 we adopt a BDF. Discretisation of Eqn. (\[eq:chemistry\]) yields $$n_i^{t + \Delta t} = \frac{C_i^{t + \Delta t}\Delta t + n_i^t}{1 + D_i^{t + \Delta t} 