.. _photoionisation:

Photoionisation of neutral species
===============================================


Ionisation by cosmic rays
---------------------------------

To model the ionisation of the neutral species by cosmic rays, we use the approximation shown in [lehtinen2007]_

.. math::

    q = q_p \frac{N}{N_p} \exp (1 - N/N_p)

with :math:`N` the neutral density. :math:`N_p` is the neutral density at the peak of the production, at altitude :math:`h_p`. :math:`q_p` is the maximum value of :math:`q`. We take:

* :math:`Q_p = 10 cm^{-3} s^{-1}`
* :math:`h_p = 15 km`


Ionisation by solar flux
---------------------------------

Once the absorbed fluxes are computed (see :ref:`absorption`), the next step is to compute the number of ions/electron pairs produced. This is done through

.. math::

    Q = N \times \sigma_N^i \times \Phi \times k_\lambda

where :math:`N` is the number density of the neutral species being ionised, :math:`\sigma_i` is the ionisation cross-section and :math:`\Phi` is the absorbed flux. :math:`k_\lambda` is the ionisation efficiency, given by [nicolet1960]_ and [bourdeau1966]_ . It is 45 for soft X-ray flux and 165 for hard X-rays, and represents the number of electron-ion pairs per photon in the solar flux.

We consider here that:

- :math:`O_2` is ionised by soft and hard X-rays to give :math:`O_2^+`
- :math:`O_2(^1 \Delta_g)` is ionised by EUV and UV flux to give :math:`O_2^+` ([pavlov2014]_, [paulsen1972]_)
- :math:`N_2` is ionised by soft and hard X-rays to give :math:`N_2^+`, which is immediately converted into :math:`0_2^+`
- :math:`NO` is ionised by EUV fluxes to give :math:`NO^+`

.. important::

    The value of :math:`Q` also has to be converted into particles/cm-3, which means converting :math:`\Phi` by multiplying it by a factor :math:`\frac{\lambda}{2 \pi \hbar c}`. :math:`\lambda` is the average wavelength here

All the cross sections are given in [pavlov2014]_ (Table 3) and its Erratum. The ionisation, absorption cross-sections and :math:`k_\lambda` (see :ref:`photoionisation`) for hard X-rays are found in [siskind2022]_ . The only exception is the :math:`k_\lambda` factor (:math:`pe/pi + 1`, with the notations in [siskind2022]_) for the ionisation of :math:`O` is 217.2 for the lowest HXR bin (it was not precised in the article, so we took the averaged value for [solomon2005]_ ).

Ionisation of :math:`O^+` and :math:`N^+`
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

In addition to the species mentioned above, the solar X-ray flux also ionises :math:`N_2` and :math:`O_2` to give :math:`O^+` and :math:`N^+`. Those two species are not present in our ions scheme, but they contribute to the ionisation as they are converted into :math:`NO^+` and :math:`O2^+`.

From the reactions rate reported in Table 6 of [pavlov2014]_ , the ratio of :math:`NO^+` to :math:`O_2^+` created can be determined. This is done in the :py:func:`~lir_achem.compute_ionisation.conversion_OpNp` function. Those added contributions to the ion densities are recomputed from the solar X-ray flux at each time-step, similarly to any other ionisation source, for both :math:`N^+` and :math:`O^+`.

Ionisation of :math:`O_2(^1 \Delta_g)`
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

As mentioned above, this is done following [paulsen1972]_ . We use the formulation proposed by [swider1979]_ which explicitely involves :py:obj:`H` and :py:obj:`Ch`. Those two quantities are therefore computed once by the :py:func:`~lir_achem.compute_source.chapman` function and stored as an attribute of the :py:class:`~lir_achem.class_definition.radiation` class.

.. math::

    q(O_2^+) = [O_2 (^1 \Delta_g)] (0.549 \times 10^{-9} \exp(-2.406\times 10^{-20} [O_2] H Ch)

    + 2.614 \times 10^{-9} \exp(-8.508\times 10^{-20} [0_2] H Ch))

In practice, this is done by calling the different functions of :py:mod:`compute_ionisation` (\ ``ci``\ below):

.. code-block::

    # Compute photon fluxes, using the taus already computed
    Phi_SXR, Phi_HXR, Phi_EUV = ci.compute_photon_flux(rad_here)

    # Compute ionisation of different species

    # O+ from O2
    ionisation_Op = ci.ion_O2_to_Op(Phi_SXR, Phi_HXR, rad_here, n_here)

    # N2+ from N2
    ionisation_N2p = ci.ion_N2_to_N2p(Phi_SXR, Phi_HXR, rad_here, n_here)

    # N+ from N2
    ionisation_Np = ci.ion_N2_to_Np(Phi_SXR, Phi_HXR, rad_here, n_here)

    # NO+ from NO
    ionisation_NOp = ci.ion_NO_to_NOp(Phi_EUV, rad_here, n_here)

    # Cosmic-rays
    ionisation_cr = ionisation_cr = ci.ion_from_cosmicrays(n_here, rad_here)

For a quiet day in winter above Italy (solar zenith angle : 58°), this gives:

.. image:: ../Figures/Ionisation_layer.png

.. important::

    We consider that all :math:`N_2^+` are converted into :math:`O_2^+` (e.g. [mitra1975]_ ).