5.1.1 Thin case

The cooling for a given molecule/atom M depends on the transition between its excited levels. It is controlled by the transitions between upper to lower levels as \(\Lambda = \sum_j [n_j \sum_{i<j} (A_{ji}\Delta E_{ji})]\) in erg/cm3/s, where \(n_j\) is the number density in the jth excited level, \(A_{ji}\) the Einstein's coefficient for the transition \(j\to i\), and \(\Delta E_{ji}\) the difference in energy between the two levels.

The main problem is to retrieve the population of the different levels, which depends on the frequency of exciting collisions and hence on the gas temperature. In particular, a molecule can be excited by collisions with a rate coefficient \(C_{ij}\) where \(j>i\), but also de-excitated by the colliders (\(C_{ji},\,j>i\)). Note the value of the coefficients is a function of the temperature \(T\) and of the type of collider \(c\), hence it should be \({}^{c}C_{ji}(T)\). At the same time the levels are de-excitated by spontaneous transitions with a rate of \(A_{ji},\,j>i\).

Considering both the transitions between the levels we obtain a differential equation for each level as \(\frac{\mathrm dn_i}{\mathrm dt}=\sum_jn_j(A_{ji}+\sum_kn_k{}^kC_{ji}) - n_i\sum_j(A_{ij}+\sum_kn_k{}^kC_{ij}) = 0\) that should be equal to zero (i.e. we want to find the steady state of the system). The term \(n_k\) represents the abundance of the kth collider, which is assumed to be in the ground state.

For a N-level system we have N-1 independent equations, hence we replace one of the equations with a continuity (mass conservation) equation \(\sum_in_i=n_\mathrm{M}\), where \(n_\mathrm{M}\) is the total amount of the given coolant molecule.

This linear system can be algebraically solved for \(N<4\), while for larger systems a LAPACK-like library is required (e.g. see dgesv method).

5.1.2 Thin case with radiation-induced transitions

When an external radiation interacts with the gas, the population of the different levels is influenced by induced transitions. In particular, for a radiation field with intensity \(I(\nu)\) we have \(n_i=\sum_jn_j(A_{ji}+B_{ji}I(\nu_{ji})+C_{ji}) - n_i\sum_j(A_{ij}+B_{ij}I(\nu_{ij})+C_{ij}) = 0\), where \(B_{ij}=A_{ij}\frac{c^2}{2h\nu^3}\) and \(g_jB_{ji}=g_iB_{ij}\) where \(g_i\) is the multiplicity of the ith level. Note here that \(\nu_{ij}\) is the frequency that corresponds to the energy of the transition \(i\to j\), i.e. \(\Delta E_{ij}=h\nu_{ij}\). See e.g. eqn.(9) of http://arxiv.org/pdf/1307.1982v2.pdf paper.

The additional terms influence the population of the different levels and then the final amount of cooling.

5.1.3 Generalization to the thick case

When the gas reaches higher densities the opacity of the gas plays a key role. The number of emitted photons is influenced by the gas beta-escape and in particular by using the weighting term \(\beta_{ij}\) for \(A_{ij}\) defined as \(\beta_{ij}=\frac{1}{1+3\tau_{ij}}\). The opacity term is \(\tau_{ij}=\frac{hc}{4\pi\left|\mathrm dv/\mathrm dz\right|} (n_jB_{ji}-n_iB_{ij})\), where \(\left|\mathrm dv/\mathrm dz\right|\) is the velocity gradient (e.g. see bit.ly/1gYbBo5).

Note that the approximation employed here is the so-called plane-parallel large velocity gradient (LVG) flow approximation. For other geometries please refer to Tab.1 of bit.ly/1gYbBo5 paper.

The system is then \(\frac{\mathrm dn_i}{\mathrm dt}=\sum_jn_j[(A_{ji}+B_{ji}I(\nu_{ji}))\beta_{ij}+C_{ji}] - n_i\sum_j[(A_{ij}+B_{ij}I(\nu_{ij}))\beta_{ji}+C_{ij}] = 0\) and the total cooling is \(\Lambda = \sum_j[ n_j (\sum_{i<j} \beta_{ij}A_{ji}\Delta E_{ji})]\).

Note that the system is no longer linear and hence different integration methods are required, as the damped affine invariant Newton method (implemented in the routine NLEQ1, http://www.zib.de/en/numerik/software/ant/nleq1.html), or the MINPACK’s hybrd and hybrj algorithms (e.g. implemented in numpy as fsolve http://docs.scipy.org/doc/scipy-0.13.0/reference/generated/scipy.optimize.fsolve.html).