\documentclass{report}

\begin{document}

\newcommand{\Ri}{\Sigma_{jl}\zeta_{i\rightarrow jl}}
\newcommand{\FE}{\frac{FE}}
\newcommand{\Parn}{\partial n_k}
\title{Jacobians for X-ray reactions}

% ****************************************************************************************************
% ****************************************************************************************************
% ****************************************************************************************************
% ****************************************************************************************************
% ****************************************************************************************************
\chapter{X-ray reactions}
\subsection*{Primary Reactions}
Reactions for primary ionization leading to one fast electron have the form:
\begin{equation}
A + \gamma_X \rightarrow A^+ + e^- \nonumber 
\end{equation}
This reaction is set for H, He, C$^+$, N$^+$, O$^+$, Ne$^+$, Na$^+$, Mg$^+$, Si$^+$, S$^+$, Ar$^+$ and Fe$^+$.\\
While reactions for primary ionization leading to two fast electrons have the form:
\begin{equation}
A + \gamma_X \rightarrow A^{2+} + e^- + e^- \nonumber 
\end{equation}
This reaction is set for Ne, Na, Mg, Si, S, Ar and Fe.

For C, N and O the treatment is more detailed, if L (K) shell is hit, 1 (2) electron is released. Carbon, Oxygen and Nitrogen are then ionized either way.


%\begin{eqnarray}
%\rm H + \gamma_X &\rightarrow& \rm H^+ + e^- \nonumber \\ 
%\rm He + \gamma_X &\rightarrow& \rm He^+ + e^- \nonumber \\
%\rm C + \gamma_X &\rightarrow& \rm C^{2+} + 2e^- \nonumber \\
%\rm N + \gamma_X &\rightarrow& \rm N^{2+} + 2e^- \nonumber \\
%\rm O + \gamma_X &\rightarrow& \rm O^{2+} + 2e^- \nonumber \\
%\rm Si + \gamma_X &\rightarrow& \rm Si^+ + e^- \nonumber \\
%\rm S + \gamma_X &\rightarrow& \rm S^{2+} + 2e^- \nonumber \\
%\rm Cl + \gamma_X &\rightarrow& \rm Cl^+ + e^- \nonumber \\
%\rm Fe + \gamma_X &\rightarrow& \rm Fe^{2+} + 2e^- \nonumber \\
%\rm C^+ + \gamma_X &\rightarrow& \rm C^{2+} + 2e^- \nonumber \\
%\rm N^+ + \gamma_X &\rightarrow& \rm N^{2+} + 2e^- \nonumber \\
%\rm O^+ + \gamma_X &\rightarrow& \rm O^{2+} + 2e^- \nonumber \\
%\rm S^+ + \gamma_X &\rightarrow& \rm S^{2+} + 2e^- \nonumber \\
%\rm Fe^+ + \gamma_X &\rightarrow& \rm Fe^{2+} + 2e^- \nonumber \\
%\rm H_2 + \gamma_X &\rightarrow& \rm H^+ + \rm H^+ + 2e^- \nonumber \\
%\rm CH + \gamma_X &\rightarrow& \rm C^{2+} + \rm H + 2e^-\nonumber \\
%\rm NH + \gamma_X &\rightarrow& \rm N^{2+} + \rm H + 2e^-\nonumber \\
%\rm OH + \gamma_X &\rightarrow& \rm O^{2+} + \rm H + 2e^-\nonumber \\
%\rm CN + \gamma_X &\rightarrow& \rm C^{2+} + \rm N + 2e^-\nonumber \\
%\rm CN + \gamma_X &\rightarrow& \rm C^+ + \rm N^+ + 2e^-\nonumber \\
%\rm CN + \gamma_X &\rightarrow& \rm C + \rm N^{2+}+ 2e^-\nonumber \\
%\rm CO + \gamma_X &\rightarrow& \rm C^{2+} + \rm O + 2e^-\nonumber \\
%\rm CO + \gamma_X &\rightarrow& \rm C^+ + \rm O^+ + 2e^-\nonumber \\
%\rm CO + \gamma_X &\rightarrow& \rm C + \rm O^{2+} + 2e^-\nonumber \\
%\rm N_2 + \gamma_X &\rightarrow& \rm N^{2+} + \rm N + 2e^-\nonumber \\
%\rm N_2 + \gamma_X &\rightarrow& \rm N^+ + \rm N^+ + 2e^-\nonumber \\
%\rm SiH + \gamma_X &\rightarrow& \rm Si^+ + \rm H + 2e^-\nonumber \\
%\rm NO + \gamma_X &\rightarrow& \rm N^{2+} + \rm O + 2e^-\nonumber \\
%\rm NO + \gamma_X &\rightarrow& \rm N^+ + \rm O^+ + 2e^-\nonumber \\
%\rm NO + \gamma_X &\rightarrow& \rm N + \rm O^{2+} + 2e^-\nonumber \\
%\rm O_2 + \gamma_X &\rightarrow& \rm O + \rm O^{2+} + 2e^-\nonumber \\
%\rm O_2 + \gamma_X &\rightarrow& \rm O^+ + \rm O^+ + 2e^-\nonumber \\
%\rm SiO + \gamma_X &\rightarrow& \rm Si + \rm O^{2+} + 2e^-\nonumber \\
%\rm SiO + \gamma_X &\rightarrow& \rm Si^+ + \rm O^+ + 2e^-\nonumber \\
%\rm CH_2 + \gamma_X &\rightarrow& \rm C^{2+} + \rm H_2 + 2e^-\nonumber \\
%\rm NH_2 + \gamma_X &\rightarrow& \rm N^{2+} + \rm H_2 + 2e^-\nonumber \\
%\rm H_2O + \gamma_X &\rightarrow& \rm O^{2+} + \rm H+2 + 2e^-\nonumber \\
%\rm HCN + \gamma_X &\rightarrow& \rm C^{2+} + \rm NH + 2e^-\nonumber \\
%\rm HCN + \gamma_X &\rightarrow& \rm N^{2+} + \rm CH + 2e^-\nonumber \\
%\rm CO_2 + \gamma_X &\rightarrow& \rm C^{2+} + \rm O_2 + 2e^-\nonumber \\
%\rm CO_2 + \gamma_X &\rightarrow& \rm O^{2+} + \rm CO + 2e^-\nonumber 
%\end{eqnarray}
%\subsection*{Secondary Reactions (16)}
%\begin{eqnarray}
%\rm H + e^- &\rightarrow&  \rm H^+ + e^-+ e^-\nonumber \\
%\rm He + e^- &\rightarrow&  \rm He^+ + e^-+ e^-\nonumber \\
%\rm H_2 + e^- & \rightarrow& \rm  H_2^+ + e^-+ e^-\nonumber \\
%\rm H_2 + e^- & \rightarrow& \rm  H^+ + H + e^-+ e^-\nonumber \\
%\rm C + e^- &\rightarrow& \rm C^+ +e^-+ e^-\nonumber \\
%\rm N + e^- &\rightarrow&  \rm N^+ +e^-+ e^-\nonumber \\
%\rm O + e^- &\rightarrow&  \rm O^+ +e^-+ e^-\nonumber \\
%\rm Si + e^- &\rightarrow& \rm   Si^+ +e^-+ e^-\nonumber \\
%\rm S + e^- &\rightarrow& \rm S^+ +e^-+ e^-\nonumber \\
%\rm Cl + e^- &\rightarrow& \rm Cl^+ +e^-+ e^-\nonumber \\
%\rm Fe + e^- &\rightarrow&  \rm Fe^+ +e^-+ e^-\nonumber \\
%\rm C^+ + e^- &\rightarrow&  \rm C^{2+}+ e^-+ e^-\nonumber \\
%\rm N^+ + e^- &\rightarrow& \rm N^{2+}+ e^-+ e^-\nonumber \\
%\rm O^+ + e^- &\rightarrow&  \rm O^{2+}+ e^-+ e^-\nonumber \\
%\rm S^+ + e^- &\rightarrow&  \rm S^{2+} +e^-+ e^-\nonumber \\
%\rm Fe^+ + e^- &\rightarrow& \rm Fe^{2+}+ e^-+ e^-\nonumber
%\end{eqnarray}
\chapter{Primary Ionization}
For reactions of the kind:
\begin{eqnarray}
 A + \gamma_X &=& A^{2+} + 2e^- \\
 AB + \gamma_X &=& A^{2+} + B + 2e^- \\
               &=& A +B^{2+} +  2e^- \\
               &=& A^{+} +B^{+} +  2e^- 
\end{eqnarray}
The dimensions of all the following quantities are also listed in the appendix.
\section{Primary rates}
\subsection{Destroying species i}
The rate of ionization for the generic reaction written above is\footnote{The quantities $\sigma$, $\tau$ and F depend on the energy. The explicit dependence is omitted for notation purposes.}:
\begin{equation}\label{primaria}
 \zeta_{i\rightarrow jl} = \int \sigma_{i\rightarrow jl} \frac{f(r)}{E} dE\qquad [s^{-1}]
\end{equation}
Where $f(r)$ is the X-ray flux (expressed in erg s$^{-1}$ cm$^{-2}$) at the radial distance $r$. Assuming exponential dumping $f(r)$ is:
\begin{equation}
f(r) = f_*\,\rm{e}^{-\tau_{\rm X}(r)}
\end{equation}
Where $f_*$ is the stellar X-ray flux in erg s$^{-1}$ cm$^{-2}$ and $\tau_{\rm X}(r)$ is the X-ray optical depth at the radius $r$. Hence the rate is:
\begin{equation}\label{primaria}
 \zeta_{i\rightarrow jl} = \int \sigma_{i\rightarrow jl} \frac{f_*}{E} e^{-\tau_{\rm X}(r)} dE\qquad [s^{-1}]
\end{equation}
Where $\sigma_i$ is the cross section of the $i$-$th$ particle which absorbes the X-ray photon in the reaction. 

We now consider the temporal behavior of the species $i$ and defining $F_i$ as $\frac{dn_i}{dt}$ we consider only the contribution of primary X-ray ionization (destruction) to the abundance of the species $i$:
\begin{equation}
 F_i = - \Ri n_i\qquad [\rm{cm^{-3}s^{-1}}]
\end{equation}
The sum in $j$ and $k$ is actually performed when the primary ionization of species $i$ has different output channels (i.e. CO). The rate does not depend on $n_{\rm{el}}$, the variation of this rate with respect to the generic species $k$ is then calculated as follows:
\begin{eqnarray}\label{steps}
 \frac{dF_i}{dn_k} &=& - \Ri \frac{\partial n_i}{\partial n_k} - \frac{\partial \left(\Ri\right)}{\partial n_k} n_i \nonumber\\
                   &=& - \Ri \delta_{ik} 
\end{eqnarray}
Since $\partial_k (\Ri)=0$.
\subsubsection{Interpretation}
\begin{itemize}
%  \item If $i \neq k$ and $n_k$ is not a species included in the total cross section we find $\frac{\partial \sigma_{tot}}{\partial n_k}=0$, then:
% \begin{equation}
%   \frac{dF_i}{dn_k} =  0
% \end{equation}
% Trivially, changing the density of a species different than i and that is not involved in the calculation of the total cross section does not change the ionization rate for species i.
\item If $i\neq k$:
\begin{equation}
  \frac{dF_i}{dn_k} =  0
\end{equation}
A small change in density for the test species $k$ does not affect the primary ionization of the $i$-$th$ species, as the process is pure photon absorption. 
\item If $i=k$:
\begin{equation}
  \frac{dF_k}{dn_k} = - \Sigma_{jl} \zeta_{k \rightarrow jl}
\end{equation}
The rate for the particle $i$ only changes with a small density change of the particle $i$ itself.
\end{itemize}
\subsection{Creating species $i$}
The rate for the ionization of the species $j$ that produces the species $i$ is:
\begin{equation}\label{pop}
 \zeta_{j\rightarrow il} = \int \sigma_{j\rightarrow il} \frac{f_*}{E} e^{-\tau_X(r)} dE
\end{equation}
Then:
\begin{equation}
 F_i = \Sigma_j \zeta_{j\rightarrow il}n_j
\end{equation}
Where the sum on the index $j$ is performed when there are multiple reactions which create the species $i$ (i.e. H$^+$ is created from primary ionization of H and H$_2$). The rate does not depend from $n_{\rm{el}}$, hence the Jacobian can be calculated as follows:
\begin{eqnarray}
 \frac{dF_i}{dn_k} &=& \Sigma_j \zeta_{j\rightarrow il} \delta_{jk} + \frac{\partial \left(\zeta_{j\rightarrow il}\right)}{\partial n_k} n_j \nonumber\\
                   &=& \zeta_{k\rightarrow il}  
\end{eqnarray}
Where we repeated the same steps as in equation (\ref{steps})
\subsubsection{Interpretation}
It must be $i\neq k$. The only way to change $F_i$ is changing the density of the species $k$.

% ****************************************************************************************************
% ****************************************************************************************************
% ****************************************************************************************************
% ****************************************************************************************************
% ****************************************************************************************************

\chapter{Secondary Ionization}
For reactions of the kind:
\begin{equation}
 A + e^{-} = A^{+} + 2e^- 
\end{equation}
\section{Species i with $i \neq \rm He, H_2$}
\subsection{Depopulating $n_i$}
The rate of collisional ionization for species $i$ is given by:
\begin{equation}\label{second}
 \zeta_{i\rightarrow j} = r_i \frac{H_x}{W_Hn_H/n_{tot}} [s^{-1}]
\end{equation}
Where $r_i$ is a parameter that takes into account the geometrical cross section of species $i$ compared to the hydrogen cross section, $n_{tot}$ is the total hydrogen nuclei density. $H_x$ is the energy deposition and $W_H$ is the mean energy per ion (Dalgarno et al., 1999). 

$W_H$ for hydrogen is given by:
\begin{equation}\label{W}
W_H = W_{0,H} \left(1+c_1\left(\frac{n_{el}}{n_{tot}}\right)^{\alpha}\right)\left(1+c_2\frac{n_{H_2}}{n_H}\right)
\end{equation}
$W_H$ is the energy needed to collisionally ionize hydrogen. The values for other elements are scaled up considering the geometrical factor $r_i = \sigma_i ^{coll}/\sigma_H ^{coll}$. W has been calculated for $He$ and $H_2$ as well, these cases will be treated separately further on.
The X-ray rate of energy  deposition is defined as:
\begin{equation}\label{hx}
 H_x = \int \sigma_{tot}Fe^{-\tau_X}dE
\end{equation}
Where $\sigma_{tot}$ is given by:
\begin{equation}
 \sigma_{tot} = \sum_{j=1}^{41} \sigma_{j} n_j/n_{tot}
\end{equation}

The rate implicitly depends on $n_{el}$. This must be taken into account when the Jacobian matrix is filled:
\begin{equation}\label{tutte}
 \frac{dF_i}{dn_k} = \frac{\partial F_i}{\partial n_k} + \frac{\partial F_i}{\partial n_{el}}\frac{\partial n_{el}}{\partial n_k}
\end{equation}

Again we consider only the contribution of the secondary ionization to the statistical equilibrium of the species $i$. $F_i$ is defined as:
\begin{equation}\label{eia}
 F_i^{ides}= - R_{i\rightarrow j} n_i = - r_i \frac{H_x}{W_Hn_H}n_{tot} n_i
\end{equation}
Where the apex ''ides'' stands for ''destruction of the species $i$ due to secondary ionization''.
Deriving $F_i^{ides}$ with respect to $n_k$ and considering that $n_{el} = \Sigma_j n_j q_j$, where $q_j$ is the charge of the $Jth$ species, we obtain:
\begin{equation}\label{big}
\frac{dF_i^{ides}}{dn_k} = -r_in_{tot}\left[\frac{n_i}{W_Hn_H}\frac{\partial H_x}{\partial n_k}-\frac{H_xn_i}{W_H^2n_H}\frac{\partial W_H}{\partial n_k}-\frac{H_xn_i}{W_Hn_H^2}\delta_{H,k}+\frac{H_x}{W_Hn_H}\delta_{ik}\right]
\end{equation}
Here we develop separately the partial derivative for all the involved terms:
\begin{eqnarray}\label{hxder}
 \frac{\partial H_x}{\partial n_k} &=& \frac{1}{n_{tot}}\frac{\partial}{\partial n_k} \int \Sigma_j n_j \sigma_j Fe^{-\tau_X}dE \nonumber \\
                                   &=& \frac{1}{n_{tot}} \int \sigma_k F e^{-\tau_X}dE
\end{eqnarray}
\begin{eqnarray}
 \frac{\partial W_H}{\partial n_k} &=& \frac{\partial}{\partial n_k} \left[W_{0,H}\left(1+c_1\left(\frac{n_{el}}{n_{tot}} \right) ^{\alpha}\right)\left(1+c_2\frac{n_{H_2}}{n_H}\right) \right]  \nonumber\\
                                   &=&  W_{0,H}\alpha c_1 \frac{n_{el}^{\alpha-1}}{n_{tot}^\alpha} q_k \left(1+c_2\frac{n_{H_2}}{n_H}\right)\nonumber\\
                                   &+& W_{0,H}\left(1+c_1\frac{n_{el}^{\alpha}}{n_{tot}^{\alpha}}\right)\frac{c_2}{n_H} \delta_{H_2,k}\nonumber\\
                                   &-& W_{0,H}\left(1+c_1\frac{n_{el}^{\alpha}}{n_{tot}^{\alpha}}\right)c_2\frac{n_{H_2}}{n_H^2} \delta_{H,k}\nonumber
\end{eqnarray}
Where $ \frac{\partial n_{el}}{\partial n_{k}} = q_k$.

Substituting these results in equation (\ref{big}) we obtain:
\begin{eqnarray}
\frac{dF_i^{ides}}{dn_k} &=& -r_i\frac{H_x}{W_Hn_H}n_{tot}\delta_{ik} + \nonumber\\
                  &-&  r_i\frac{H_x}{W_H^2n_H^2}\frac{n_{tot}n_{H_2}n_i}{n_{H}}c_2W_{0,H} \left(1+c_1\frac{n_{el}^{\alpha}}{n_{tot}^{\alpha}}\right)\delta_{H,k}\nonumber\\
                  &+& r_i\frac{H_x}{W_H^2n_H^2}n_{tot}\,n_i\,c_2\,W_{0,H}\left(1+c_1\frac{n_{el}^{\alpha}}{n_{tot}^{\alpha}}\right)
                  \delta_{H_2,k}\nonumber\\
                  &+&  r_i \frac{H_x}{W_H^2}\frac{n_i}{n_H} W_{0,H}\,\alpha\, c_1  \left(\frac{n_{el}}{n_{tot}}\right)^{\alpha-1} q_k\left(1+c_2\frac{n_{H_2}}{n_H}\right) \nonumber\\
                  &+& r_i\frac{H_xn_i}{W_Hn_H^2}n_{tot}\delta_{H,k} \nonumber\\
                  &-&  r_i\frac{n_i}{W_Hn_H}\int \sigma_k F e^{-\tau_X}dE
\end{eqnarray}
\subsubsection{Interpretation}
The previous equation is composed of 6 terms, which are described as follows:
\begin{enumerate}
 \item This is the derivative of reaction rate. It is only non--zero for the species $i$ itself.
 \item The same as the previous point, only considering $H_2$ instead of $H$.
 \item This term only appears if the density of hydrogen changes, since it explicitly contains the dependence of $W_H$ on $n_H$.
 \item This term only appears if $k=H$, but if $i=k=H$ this term compensates with the first and the sum is zero. This is due to the fact that $F_i$ (\ref{eia}) does not depend explicitly (but only via $W_H$) on $n_H$ when $i=H$.

 \item Charge effect: this term contains the dependence of $n_{el}$ on $n_k$. If species $k$ is not charged this contribute is zero, since a change of density of a neutral species does not affect the electron density. 

If we increase the density of a positively charged species the electron density increases as well. This causes $W_H$ to increase (eq. \ref{W}), and then the total rate to decrease. Increasing the number of electrons in the gas phase causes an increase in the mean energy that is needed to form an ion pair, since thermal electrons also interact with fast electrons. This results in a positive contribution to $F_i$.

If the density of a negative species is increased the inverse effect takes place. It means that $n_{el}$, and then $W_H$, decrease. The total rate for the ionization of species $i$ increases and the contribution to $F_i$ is negative. Less electrons in the gas phase decreases Coulomb interaction in favor of collisional interactions with atoms and molecules, thus causing more secondary ionization.
 \item If the species $k$ contributes in building the cross section for X-ray absorption, then $H_x$ changes as well.
\end{enumerate}
\subsection{Populating $n_i$, with $i \neq He, H_2$}
The rate of collisional ionization for species $j$ that leads to the formation of species $i$ is given by:
\begin{equation}\label{second}
 \zeta_{j\rightarrow i} = r_j \frac{H_x}{W_Hn_H/n_{tot}} [s^{-1}]
\end{equation}
The analytical expression for the element $ik$ of the Jacobian matrix is exactly the same as in the previous case, but with a sign inversion since we are now considering the formation process for the products. Taking into account the correct indexes we obtain:
\begin{eqnarray}
\frac{dF_i^{ifor}}{dn_k} &=& r_j\frac{H_x}{W_Hn_H}n_{tot}\delta_{jk} + \nonumber\\
                  &+&  r_j\frac{H_x}{W_H^2n_H^2}\frac{n_{tot}n_{H_2}n_j}{n_{H}}c_2W_{0,H} \left(1+c_1\frac{n_{el}^{\alpha}}{n_{tot}^{\alpha}}\right)\delta_{H,k}\nonumber\\
                  &-& r_j\frac{H_x}{W_H^2n_H^2}n_{tot}\,n_j\,c_2\,W_{0,H}\left(1+c_1\frac{n_{el}^{\alpha}}{n_{tot}^{\alpha}}\right)
                  \delta_{H_2,k}\nonumber\\
                  &-&  r_j \frac{H_x}{W_H^2}\frac{n_j}{n_H} W_{0,H}\,\alpha\, c_1  \left(\frac{n_{el}}{n_{tot}}\right)^{\alpha-1} q_k\left(1+c_2\frac{n_{H_2}}{n_H}\right) \nonumber\\
                  &-& r_j\frac{H_xn_j}{W_Hn_H^2}n_{tot}\delta_{H,k} \nonumber\\
                  &+&  r_j\frac{n_j}{W_Hn_H}\int \sigma_k F e^{-\tau_X}dE
\end{eqnarray}
Where the apex ''ifor'' stands for ''formation of the species $i$ via ionization of species $j$''.
\section{$H_2$}
\subsection{Depopulating $n_{H_2}$}
The rate for $H_2$ collisional ionization is given by:
\begin{equation}\label{rateh2}
 \zeta_{H_2}^{sec} = \frac{H_x}{W_{H_2}n_{H_2}}n_{tot}
\end{equation}
Where $W_{H_2}$ is given by:
\begin{equation}
 W_{H_2} = W_{0,H_2}\left(1 + c_1 \left(\frac{n_{el}}{n_{tot}}\right)^{\alpha}\right)\left(1 + c_2\frac{n_H}{n_{H_2}}\right)
\end{equation}
Then $F_{H_2}^{des}$ for the depopulating case is:
\begin{equation}\label{fh2}
 F_{H_2}^{des} = - \zeta_{H_2}^{sec}n_{H_2}
\end{equation}
Substituting (\ref{rateh2}) in (\ref{fh2}) we obtain:
\begin{equation}
 F_{H_2}^{des} = - \frac{H_x}{W_{H_2}}n_{tot}
\end{equation}
Deriving $F_{H_2}^{des}$ with respect to $n_k$ we get:
\begin{equation}
 \frac{dF_{H_2}^{des}}{dn_k} = - n_{tot}\left[\frac{1}{W_{H_2}}\frac{\partial Hx}{\partial n_k} - \frac{H_x}{W_{H_2}^2}\frac{\partial W_{H_2}}{\partial n_k}\right]
\end{equation}
Here we compute the terms listed above:
\begin{eqnarray}
\frac{\partial W_{H_2}}{\partial n_k} &=& W_{0,H_2}\left(1+c_2 \frac{n_H}{n_{H_2}}\right)\alpha c_1 \frac{n_{el}^{\alpha-1}}{n_{tot}^{\alpha}}q_k  \nonumber\\
                                      &+& W_{0,H_2}\left(1+c_1\left(\frac{n_{el}}{n_{tot}}\right)^{\alpha}\right) c_2\frac{1}{n_{H_2}}\delta_{H,k}\nonumber\\
                                      &-& W_{0,H_2}\left(1+c_1\left(\frac{n_{el}}{n_{tot}}\right)^{\alpha}\right)c_2 \frac{n_H}{n_{H_2}^2}\delta_{H_2,k}
\end{eqnarray}
Using (\ref{hxder}) we obtain:
\begin{eqnarray}
 \frac{dF_{H_2}^{des}}{dn_k} &=& \frac{H_x}{W_{H_2}^2}\frac{n_{tot}}{n_{H_2}}c_2\,W_{0,H_2}\left(1+c_1\left(\frac{n_{el}}{n_{tot}}\right)^{\alpha}\right) \delta_{H,k}\nonumber\\
                       &-& \frac{H_x}{W_{H_2}^2}\frac{n_{tot}n_H}{n_{H_2}^2}c_2\,W_{0,H_2}\left(1+c_1 \left(\frac{n_{el}}{n_{tot}}\right)^{\alpha}\right)\delta_{H_2,k}\nonumber\\
                       &+& \frac{H_x}{W_{H_2}^2}W_{0,H_2}\alpha c_1\left(\frac{n_{el}}{n_{tot}}\right)^{\alpha-1}q_k\left(1+c_2 \frac{n_H}{n_{H_2}}\right)  \nonumber\\
                       &-& \frac{1}{W_{H_2}} \int \sigma_k(E)F(E,r) dE\nonumber
\end{eqnarray}
\subsection{Populating $H_2$}
It can be easily shown that:
\begin{equation}
 F_{H_2}^{form} = -F_{H_2}^{des}
\end{equation}
Then:
\begin{eqnarray}
 \frac{dF_{H_2}^{form}}{dn_k} &=& - \frac{H_x}{W_{H_2}^2}\frac{n_{tot}}{n_{H_2}}c_2\,W_{0,H_2}\left(1+c_1\left(\frac{n_{el}}{n_{tot}}\right)^{\alpha}\right) \delta_{H,k}\nonumber\\
                       &+& \frac{H_x}{W_{H_2}^2}\frac{n_{tot}n_H}{n_{H_2}^2}c_2\,W_{0,H_2}\left(1+c_1 \left(\frac{n_{el}}{n_{tot}}\right)^{\alpha}\right)\delta_{H_2,k}\nonumber\\
                       &-& \frac{H_x}{W_{H_2}^2}W_{0,H_2}\alpha c_1\left(\frac{n_{el}}{n_{tot}}\right)^{\alpha-1}q_k\left(1+c_2 \frac{n_H}{n_{H_2}}\right)  \nonumber\\
                       &+& \frac{1}{W_{H_2}} \int \sigma_k(E)F(E,r) dE\nonumber
\end{eqnarray}
\section{$He$}
\subsection{Depopulating $He$}
The only reaction involved is:
\begin{equation}
 He + e^- \rightarrow He^+ + 2e^-
\end{equation}
The rate is:
\begin{equation}
\zeta_{He}^{sec} = \frac{H_x}{W_{He}n_{He}}n_{tot}
\end{equation}
Then:
\begin{equation}
 F_{He}^{des} = - \frac{H_x}{W_{He}}n_{tot}
\end{equation}

$W_{He}$ is given by:
\begin{equation}
 W_{He} = W_{0,He}\left(1 + c_1\left(\frac{n_{el}}{n_{tot}}\right)^{\alpha}\right)
\end{equation}
The Jacobian term can be written as:
\begin{equation}
 \frac{dF_{He}^{des}}{dn_k} =  -n_{tot}\left[\frac{1}{W_{He}}\frac{\partial H_x}{\partial n_k}-\frac{H_x}{W_{He}^2}\frac{\partial W_{He}}{\partial n_k}\right]
\end{equation}
Where:
\begin{eqnarray}
\frac{\partial W_{He}}{\partial n_k} &=& W_{0,He}\,\alpha\,c_1 \frac{n_{el}^{\alpha-1}}{n_{tot}^{\alpha}}q_k  \nonumber\\
\end{eqnarray}
Using again (\ref{hxder}) we obtain:
\begin{eqnarray}
\frac{dF_{He}^{des}}{dn_k} &=& \frac{H_x}{W_{He}^2} W_{0,He}\,\alpha\,c_1 \left(\frac{n_{el}}{n_{tot}}\right)^{\alpha-1}q_k \nonumber\\
                           &-& \frac{1}{W_{He}} \int \sigma_k(E)F(E,r) dE\nonumber
\end{eqnarray}
\subsection{Populating $He$}
Given that:
\begin{equation}
 F_{He}^{form} = -F_{He}^{des}
\end{equation}
We obtain:
\begin{eqnarray}
\frac{dF_{He}^{form}}{dn_k} &=& -\frac{H_x}{W_{He}^2} W_{0,He}\,\alpha\,c_1 \left(\frac{n_{el}}{n_{tot}}\right)^{\alpha-1}q_k \nonumber\\
                           &+& \frac{1}{W_{He}} \int \sigma_k(E)F(E,r) dE\nonumber
\end{eqnarray}
\chapter*{Appendix}
Dimensions:
\begin{eqnarray}
 \sigma &=& [cm^2]  \nonumber\\
 F      &=& [erg\,s^{-1} cm^{-2} erg^{-1}]\nonumber\\
 N_H      &=& [cm^{-2}]\nonumber\\
 H_x      &=& [erg s^{-1} H^{-1}]\nonumber\\
 F_i      &=& [cm^{-3}s^{-1}]\nonumber\\
 W_H      &=& [erg]\nonumber\\
 n_i      &=& [cm^{-3}]\nonumber\\
 W_{0,H}      &=& [erg]\nonumber
\end{eqnarray}
\tableofcontents
\end{document}