{\ensuremath{\mbox{{H$\,${\sc i}}}}}$, ${\mbox{{He$\,${\sc i}}}}$, and ${\mbox{{He$\,${\sc ii}}}}$. The total photo-heating rate is defined by ${\cal H} = \sum_i {\cal H}_i n_i$. The contribution from a point-like source, $\alpha$, is written as $${\cal H}_{i, \alpha} = \frac{1}{4 \pi r_\alpha^2} \int_{\nu_i}^\infty 
\frac{{{\rm d}}\nu}{\nu} \sigma_i(\nu) L_\alpha(\nu) (\nu - \nu_i) 
  \exp\left[- \sum_j \sigma_j(\nu) {\cal N}_{j, \alpha} \right]. 
\label{eq:photoheat}$$ As for the photoionization, we generate a look-up table for each species when all sources have the identical spectral shape.

We solve the energy equation for each mesh implicitly as $$u^{t + \Delta t} = u^{t} + \frac{{\cal H}^{t + \Delta t} - \Lambda(n_i^{t + \Delta t}, T^{t+\Delta t})}{\rho^t} \Delta t, 
\label{eq:heatcool}$$ where $u$ and $T^t = T(n_i^t, u^t)$ are respectively the specific internal energy and temperature of the gas and $\Lambda$ is the cooling function. Although this implicit integration is always stable, we need to subcycle the energy solver with $\Delta t_{\rm chem}$ because both $C_i$ and $D_i$ in Eqn. (\[eq:bdf\]) are functions of the temperature. We thus perform the rate solver and the energy solver alternately. The chemical timestep $\Delta t_{\rm chem}$ is recalculated before every subcycle.

Chemical reaction and cooling rates
-----------------------------------

We try to use the chemical reaction and cooling rates as 