Hydrogen Non-Equilibrium Ionisation Effects in Coronal Mass Ejections
P. Pagano A. Bemporad D. H. Mackay
astro-ph.SR
 


   

    

   

Pagano et al.
Ionisation CMEs

School of Mathematics and Statistics, University of St Andrews, North Haugh, St Andrews, Fife, Scotland KY16 9SS, UK 
            INAF – Turin Astrophysical Observatory, via Osservatorio 20, 10025 Pino Torinese (TO), Italy 
 

      pp25@st-andrews.ac.uk
      

   



 
  
   A new generation of coronagraphs to study the solar wind and Coronal Mass Ejections (CMEs) are being developed and launched. These coronagraphs will heavily rely on
   multi-channel observations where visible light (VL) and UV-EUV observations provide new plasma diagnostics.
   One of these instruments, Metis on board ESA-Solar Orbiter, will simultaneously observe VL and the UV Lyman-α line.
   The number of neutral Hydrogen atoms (a small fraction of coronal protons) is a key parameter for deriving plasma properties such as temperature from the observed Lyman-α line intensity.
   However, these measurements are significantly affected if non-equilibrium ionisation effects occur, which can be relevant during CMEs. 
  
   The aim of this work is to determine if non-equilibrium ionisation effects are relevant in CMEs and in particular
   when and in which regions of the CME plasma ionisation equilibrium can be assumed for data analysis.
  
   We use a magneto-hydrodynamic simulation of a magnetic flux rope ejection to generate a CME. From this we then reconstruct the ionisation state of Hydrogen atoms in the CME by evaluating both the advection of neutral and ionised Hydrogen atoms and the ionisation and recombination rates in the MHD simulation.
  
   We find that the equilibrium ionisation assumption holds mostly in the core of the CME, which is represented by a magnetic flux rope. In contrast
   non-equilibrium ionisation effects are significant at the CME front, where
   we find about 100 times more neutral Hydrogen atoms than prescribed by ionisation equilibrium conditions, even if this neutral Hydrogen excess might be difficult to identify due to projection effects.
  
   This work provides key information for the development of a new generation of diagnostic techniques that aim at combining visible light
   and Lyman-α line emissions. The results show that non-ionisation equilibrium effects need to be considered when we analyse CME fronts.
   To  incorrectly assume equilibrium ionisation in these regions would lead to a systematic underestimate of plasma temperatures.

   

   Hydrogen Non-Equilibrium Ionisation Effects in Coronal Mass Ejections
    P. Pagano<ref> A. Bemporad<ref> D.H. Mackay<ref>
          
     
=====================================================================







§ INTRODUCTION


Coronal Mass Ejections (CMEs) are very dynamic and violent phenomena occurring in the solar corona.
Despite continuous observations and monitoring of these events
from ground and space
over the last decades, their elusive nature still presents a major challenge for solar physics. One of the main reasons is that measuring properties of the embedded plasma is always an arduous task. It is complicated by the optical thinness of the coronal plasma, the large speeds and rapid evolution of the CME, and the rapidly changing temperatures.
Due to this, remote sensing observations acquired at different wavelengths are combined and analysed with a number of diagnostic techniques
to derive the thermodynamic properties of the coronal plasma ejected in CMEs.

New diagnostic techniques currently under development
are based on the forthcoming remote sensing observations of the ultraviolet (UV) 
Lyman-α line emitted from neutral Hydrogen atoms. 
These diagnostics will be applied to
future data that will be acquired for instance by Metis <cit.>, a new generation multi-channel coronagraph on board the ESA-Solar Orbiter mission
that will observe at the same time the solar corona in visible light (VL, polarised and non-polarised components) and in the Lyman-α line.
This will be done through the simultaneous acquisition of observations from a broad-band VL channel (580-640 nm)
and  a narrow-band one in the UV emission around 121.6 nm, where neutral Hydrogen atoms emit the Lyman-α line.

At the same time, other forthcoming missions will investigate the solar corona with UV observations centred on the Lyman-α line (the strongest line in the whole EUV coronal spectrum),
such as the LST instrument that will fly on-board the Chinese ASO-S mission in 2022 <cit.>.
The key strength of these new instruments is the switch to multi-channel observations that
will allow the investigation of plasma properties that so far have been impossible to access due to the limitation of single channel VL observations.

These new observations will pose unprecedented challenges, because it will be essential 
to have the proper interpretative tools to obtain accurate derivations of CME properties.
For example, the simultaneous observation of VL (mostly photospheric emission scattered by free electrons in the solar corona)
and UV Lyman-α (emitted by neutral Hydrogen atoms mostly excited by chromospheric radiation) can open up a number of diagnostic possibilities. In a plasma, the intensity of the emissions coming from different ionisation stages of the same atomic species are 
connected to the abundance of each ionisation level, and can thus be a useful proxy for the plasma temperature. On the other hand, the Lyman-α emission from neutral H atoms can be used to derive the plasma temperatures by estimating the fraction of neutrals, once the plasma densities are known from VL data. 

These techniques were extensively discussed for instance by <cit.> who demonstrated, by using SOHO/LASCO and UVCS observations, that CME plasma temperatures can be derived by a direct combination of VL and UV Lyman-α intensities
even when neglecting spectroscopic information. Similar data analysis methods were also tested by <cit.>, based on the analysis of synthetic data and under different assumptions mostly related with the integration along the line-of-sight (hereafter LOS) through the optically thin plasma. These techniques will likely provide the first ever 2D maps of plasma temperatures inside CMEs in the intermediate corona, and were also strengthened by the improved methods developed by <cit.>
and <cit.> to derive the column density of the CME plasma from 
polarised VL observations. Moreover, because of the Doppler dimming <cit.>, the UV Lyman-α emission by coronal atoms also depends on their radial outflow speed, and thus the derivation of temperatures in CMEs also requires the technique recently developed by <cit.> to reconstruct the 2D distribution of plasma speed in the CME bodies from a sequence of VL images.

For all of the techniques described above to be accurate, we need to understand 
how the abundance of various ionisation levels evolves in dynamic events such as CMEs.
In this context, the key question is whether the ionisation equilibrium conditions
hold for Hydrogen atoms throughout the propagation of CMEs in the solar corona 
or is there instead a significant departure from ionisation equilibrium. The ionisation equilibrium hypothesis corresponds to the assumption that a balance exists between all processes causing the stripping of electrons from atoms of a
given element with all the recombination processes leading to lower ionisation stages.
For a long period, this simplifying hypothesis has been made in the analysis of UVCS 
observations of CMEs and post-CME current sheets in the intermediate corona <cit.>, where ionisation equilibrium was assumed to derive plasma temperatures.

The issue of the ion distribution in plasmas that are  involved in dynamic events
has puzzled solar physicists for decades, as non-equilibrium ionisation effects can become important in a number of scenarios.
<cit.> and <cit.> addressed the ionisation of the chromosphere during solar flares.
<cit.> modelled the non-ionisation equilibrium effects following 
a heating event due to reconnection and found the locations where these effects can be relevant.
<cit.> found that non-equilibrium ionisation effects are relevant during nanoflares.
<cit.> studied the non-equilibrium ionisation in the solar wind.
<cit.> described how the non-equilibrium ionisation affects the formation of IRIS lines. Recently, the ionisation equilibrium assumption was (often tacitly) 
performed by many authors to derive the 2D distributions of plasma temperatures during CMEs 
<cit.> from the analysis of EUV images acquired with the multiple filters provided by the AIA telescopes on-board SDO. 
In fact, the applied emission measure techniques
require the construction of synthetic EUV spectra for all lines emitted by different ions
in each instrument band-pass at various temperatures, to derive the instrument 
temperature responses. These synthetic spectra can be constructed so far
<cit.> only by assuming ionisation equilibrium.

On the other hand, departures from ionisation equilibrium are expected to occur in many
dynamic phenomena in the solar atmosphere, as suggested 
by the significant spectral variability observed to occur on timescales shorter than the ion 
equilibration times. The principles of a numerical tool for the calculation of non-equilibrium 
ionisation states in the solar corona were illustrated by <cit.>, and
calculations of the ionisation equilibrium timescales have been performed by
several authors, such as <cit.> who derived the ionisation equilibrium timescales in the absence of flows. Recently <cit.> reviewed the state of the art for these methods.
The departure from ionisation equilibrium in CME shocks for Oxygen and Silicon was justified with a numerical model by <cit.>.
The violation of the 
ionisation equilibrium assumption in the post-CME EUV dimming region was 
reported by <cit.> in the higher temperature range. Later
<cit.> showed that, because of non-equilibrium ionisation effects, 
the analysis of EUV imaging observations of post-CME current sheets would lead 
to an under (over) estimate of plasma temperatures at low (large) 
heights by about a factor of two if equilibrium is assumed. On the other hand,
<cit.> reported that in a CME core, departures from ionisation 
equilibrium were within 10% for Hydrogen and Helium atoms in the erupting filament.
Very recently, possible observational consequences of non-equilibrium ionisation
in the EUV imaging of CME related phenomena (such as post-flare arcades and CME-driven shocks)
were investigated by <cit.>: the authors re-analysed previous dynamic events
studied with an emission measure technique applied to EUV and soft X-ray images and 
concluded that the temperature of hotter plasma not in ionisation equilibrium is 
lower than the temperature calculated by assuming ionisation equilibrium. These
effects can be important not only for the study of CMEs, but also for elemental
abundance determination in stationary solar wind flows <cit.>.

Similar analyses have often been performed for observations of erupting 
events in the transition region and in the inner corona (R < 1.5 R_⊙),
but not very often for CMEs observed higher up in the intermediate corona during the expansion phase (R > 1.5 R_⊙). At these altitudes the problem of ionisation equilibrium
was discussed for instance by <cit.>, dealing with UVCS observations of a CME 
and the reconstruction of the thermal evolution of the observed plasma. More recently <cit.> and <cit.> have reported similar observations of CMEs with UVCS 
spectra and other EUV imagers 
and spectrometers. In the interpretation of UVCS observations the problem was investigated by <cit.> also with numerical simulations, inferring the plasma heating rate in a CME by assuming ionisation equilibrium. Later, <cit.> simulated the UV spectroscopic emission from a CME (focusing on the evolution during a shock transit) as observed by UVCS 
by including also the effects of non-equilibrium ionisation. More recently, a non-equilibrium ionisation analysis
was performed by <cit.> to derive the physical parameters of an erupting prominence
observed spectroscopically with UVCS. While a significant amount of work has been carried out, none of these works focused on the 
evolution of neutral Hydrogen ionisation states for the analysis of the Lyman-α line
intensity evolution observed during CMEs.

Hence in summary, the interpretation of coronagraphic observations of CMEs in the Lyman-α line that will be provided by forthcoming instruments such as Metis on-board Solar Orbiter and LST on-board ASO-S missions will require a deep understanding of which regions in CME bodies can be considered in ionisation equilibrium.
In order to answer this question, in this paper we adopt a theoretical approach by using our realistic numerical modelling of the propagation of a CME to study the evolution of the ionisation state of
Hydrogen atoms during the CME. Through this we can relate our results with the ionisation state
that is prescribed by ionisation equilibrium conditions.

To reach our goals we start from the MHD simulation of <cit.> of a magnetic flux rope ejection.
Magnetic flux ropes are structures found in the solar corona
that usually lie on polarity inversion lines <cit.>
that are known to suddenly erupt and produce CMEs <cit.>. 
We first use a magnetofrictional model to describe the formation phase of magnetic flux rope using a quasi-static and magnetically dominated approach
and we switch to MHD simulations when the magnetic configuration is unstable.
In the MHD simulations, we solve the general form of the MHD equations
that are appropriate to describe a dynamic evolution where the dominant forces can be magnetic, pressure gradients, or gravity. 

Once we have solved the MHD equations for the eruption of the magnetic flux rope, 
we devise a technique to reconstruct the history of the ionisation state of Hydrogen atoms in the 3D domain of the MHD equations as a function of time.
Such an approach has already been successfully adopted by <cit.>
for the study of the emission of Oxygen and Silicon lines in shocks connected with CMEs.
In order to reconstruct the ionisation state of Hydrogen in our MHD domain,
we separately estimate the advection of neutral and ionised Hydrogen atoms
and the local rate of ionisation and recombination processes. Finally we combine
these two terms to obtain the abundances of neutral and ionised Hydrogen atoms 
as a function of space and time.

The paper is structured as follows.
In Sec.<ref> we summarise the key properties of the MHD simulation 
we use here and in Sec.<ref> we describe the technique we use to reconstruct the history of the ionisation state of Hydrogen atoms. Next
in Sec.<ref> we illustrate the results of our study and we draw our conclusions in Sec.<ref>.



§ MHD SIMULATION



In the present paper, we analyse 
the ionisation state of Hydrogen atoms during a CME
starting from the data
obtained in the MHD simulation of <cit.>.
In <cit.> the simulation of magnetic flux rope ejection is carried out by coupling the magnetofrictional model of <cit.> with a magnetohydrodynamic (MHD) simulation.
This technique, initiated by <cit.> and used in a number of studies <cit.>, has proven to be effective in modelling the full life span of magnetic flux ropes.

In that simulation, we used the MPI-AMRVAC software <cit.> to solve the MHD equations,
where solar gravity, anisotropic thermal conduction, and optically thin radiative losses are treated as source terms:

    ∂ρ/∂ t+∇⃗·(ρv⃗)=0,


    ∂ρv⃗/∂ t+∇⃗·(ρv⃗v⃗)
       +∇ p-(∇⃗×B⃗)×B⃗)/4π=+ρg⃗,


    ∂B⃗/∂ t-∇⃗×(v⃗×B⃗)=0,


    ∂ e/∂ t+∇⃗·[(e+p)v⃗]=ρg⃗·v⃗-n^2χ(T)-∇·F⃗_⃗c⃗,

where t is the time, ρ the density,
v⃗ velocity, p thermal pressure, B⃗ magnetic field,
e the total energy,
n number density, F⃗_⃗c⃗ the conductive flux according to <cit.>,
and χ(T) the radiative losses per unit emission measure <cit.>.
To close the set of Eqs. <ref>-<ref> we have a relation between internal, total, kinetic, and magnetic energy

    p/γ-1=e-1/2ρv⃗^2-B⃗^2/8π,

where γ=5/3 denotes the ratio of specific heat, and
the expression for solar gravitational acceleration is

    g⃗=-G M_⊙/r^2r̂,

where G is the gravitational constant, M_⊙ denotes the mass of the Sun,
r is the radial distance from the center of the Sun
and r̂ is the corresponding unit vector.
The spatial domain of this simulation extends over 3 R_⊙ in the radial direction starting from
r=R_⊙. The colatitude, θ, spans from θ=30^∘ to
θ=100^∘ and the longitude, ϕ, spans over 90^∘.
For more details,
we refer the reader to a series of works:
<cit.>
where the initial pre-eruptive configuration is derived using a magnetofrictional model,
<cit.> where for the first time we coupled
the eruptive initial conditions with the MHD simulations,
<cit.> where we study the effect 
of gravity on the propagation of CMEs in the corona, 
and finally <cit.> where
we included the effects of non-ideal MHD.

In the magnetofrictional model two neighbouring magnetic bipoles are allowed to 
evolve 
under the effect of differential rotation, meridional flows and surface diffusion.
As one of the two bipoles is initially more sheared than the other
a magnetic flux rope forms in the corona above its polarity inversion line (PIL) after 19 days of evolution <cit.>.
At the same time, due to the magnetic field configuration another PIL separates the two bipoles 
and above this one a null point is present.  With evolution and the formation of a flux rope the magnetic field configuration obtained is not stable and a significant Lorentz force is present underneath the magnetic flux rope.
In this MHD simulation, we also construct an atmosphere around the magnetic field configuration
where the magnetic flux rope is modelled initially to be colder and denser 
than its surroundings.
This is done by imposing a temperature distribution dependent of the θ-component of the magnetic field, so to have cold plasma where the we have a magnetic flux rope that develops with an θ-directed axial magnetic field in this configuration and, at the same time, maintaining an horizontal thermal pressure balance,
so that these cold regions are also denser than their surroundings. 
In the initial condition of this simulation the minimum temperature is set to T_min=10^4 in a corona at T=2 MK, leading to a magnetic flux rope structure that is up to 200 times denser than the background corona.
More details on how this atmosphere is constructed can be found in <cit.>.

The flux rope is then ejected out of the corona due to
the Lorentz force distribution that builds-up over 19 days of evolution.
Figure <ref> shows the simulation after 19.7 minutes of evolution for both the column density (Fig <ref>a) and the density-weighted temperature along the LOS, T_LOS (Fig <ref>b)

    T_LOS=∫_-∞^+∞ρ T   dz/∫_-∞^+∞ρ  dz
    .



The key features of this evolution are that we find plasma travelling outwards
because of the ejection and the plasma distribution closely resembles a 
three-component CME,
as a dense front is followed by a denser core with a void in between.
A similar pattern is found in the temperature map, 
where we find a hot front of the ejection due to the compression where
a colder region lags behind.

One important aspect of the simulation is that the direction along which the magnetic flux rope expands is not radial, and instead the magnetic flux rope is ejected in the direction of the null point
that lies between the two bipoles.
Therefore, knowing the direction of propagation and the location of the magnetic flux rope and the CME front on this direction as function of time, we can
follow the plasma density and temperature inside these structures.
Fig.<ref> shows the density and temperature at the centre of the flux rope and at the CME front as a function of time.


We find that the densities of both the centre of the flux rope and the CME front decrease in time
as all structures expand and move into regions of lower density higher in the solar corona.
The evolution of the plasma temperature is instead more complex.
As explained in <cit.>, the temperature of the magnetic flux rope immediately increases during
the very early stage of the MHD simulation due to the conversion of magnetic energy into thermal energy.
Later on, as the magnetic flux rope expands adiabatically, it also slowly cools down.
The CME front shows instead a constant temperature increase.
It basically coincides with the magnetic flux rope at the beginning of the simulation (t<6 min)
as the two structures are very close and it is not possible to identify a very distinct CME front.
Afterwards the compression occurring at the CME front heats the plasma to over 10 MK,
where it more or less plateaus till the end of the simulation.
While a temperature significantly above 10 MK may seem unrealistic in the solar corona, such values
are nevertheless in the correct order of magnitude 
for the tenuous plasma of the outer corona violently heated by the compression generated from a violent CME.



§ IONISATION STATE OF PLASMA



In this section we illustrate how we reconstruct the ionisation state of the Hydrogen atoms in the plasma 
as a function of space and time in our MHD simulation.
This is done by post-processing the results of the MHD simulation, 
thus assuming that the ionisation state of Hydrogen does not affect
the dynamic and thermal evolution of the CME in the MHD simulation.
This assumption is based on the fact that
the amount of neutral Hydrogen atoms is usually many orders of magnitude smaller 
than coronal protons and electrons, and in collisionless plasma even the very few 
neutral atoms are not interacting with other particles neither by collisions, nor by
magnetic fields. This strategy has been already been adopted in <cit.> 
for the derivation of the ionisation state of Oxygen and Silicon.

In order to reconstruct the ionisation state of Hydrogen we need to solve the following coupled differential equations

    ∂ n_HI/∂ t=-∇⃗·(n_HIv⃗)+n_HII(α_HIIn_HII-q_HIn_HI)


    ∂ n_HII/∂ t=-∇⃗·(n_HIIv⃗)+n_HII(q_HIn_HI-α_HIIn_HII)

where n_HI and n_HII are the numbers densities of neutral and ionised Hydrogen atoms respectively as function of space and time,
α_HII is the recombination rate by collisions of ionised Hydrogen,
q_HI is the ionisation rate by collisions of neutral Hydrogen,
t is the independent variable time and v⃗ is the velocity field that is input from the MHD simulation.
As we have the velocity field only at the output cadence of the MHD simulation, that is 34.8 second,
we linearly interpolate for intermediate times.
The values of α_HII and q_HI are taken from CHIANTI spectral code <cit.> as a function of the temperature T of the plasma
and number density n_e of the electrons.

The general solution of these equations provides the temporal and spacial dependent abundance of 
n_HI and n_HII, including solutions that describe the 
out of ionisation equilibrium configurations.
In order to solve these equations, we assume an initial condition at t=0
where the Hydrogen atoms are in ionisation equilibrium.

    n_HII=ρ/m_p


    [α_HIIn_HII-q_HIn_HI]=0.

Eq.<ref> expresses that the number of ionised Hydrogen atoms is the same as the number density of electrons as described by the MHD equations,
and in Eq.<ref> we derive the number of neutral Hydrogen atoms from the number of ionised Hydrogen atoms when the number of recombinations 
matches the number of ionisations.
As the temperature is not uniform at t=0 in our MHD simulation, neither is the ratio n_HI/n_HII.



 §.§ Solution of the ionisation state of Hydrogen


In order to solve the two differential equations (Eq.<ref>-<ref>) that govern the 
ionisation state of Hydrogen we need to evaluate 
the left hand side of both equations to obtain the rate of change
for n_HII and n_HI as a function of space at a given time t.
To do so, we assume that we can split the right hand side terms of 
Eq.<ref> and Eq.<ref>
into two terms, the advection term and the ionisation/recombination term,
and that accordingly we can integrate Eq.<ref> and Eq.<ref> in time
to obtain the finite variation of n_HII and n_HI over a finite integration time:

    Δ n_HI=[ADV_HI]+[IR_HI]


    Δ n_HII=[ADV_HII]+[IR_HII].

We then separately evaluate each term on the right hand side of Eq.<ref> and Eq.<ref>
and add their contributions to evaluate the left hand side.
In the next two sections we describe how we evaluate the advection terms and the ionisation/recombination terms.



  §.§.§ Advection term

The advection term describes the motion of the Hydrogen atoms (neutrals and ions) 
due to the velocity field derived from the solution of the MHD equations.
For neutral Hydrogen atoms this term can be written as:

    [ADV_HI]=-∇⃗·(n_HIv⃗) Δ t

where Δ t is the finite integration time step.
In order to estimate [ADV_HI] and [ADV_HII] we use an explicit Godunov scheme <cit.>
where we use boundary conditions consistent with the ones adopted in solving 
the MHD equations.
In the Godunov scheme we need to estimate the flux of ionised F⃗^⃗n⃗_⃗H⃗I⃗I⃗ and neutral F⃗^⃗n⃗_⃗H⃗I⃗ Hydrogen atoms  that cross each
face of the finite difference cells of the MHD simulation grid. 
In order for this explicit scheme to work, we need to limit
the integration time step, Δ t, so that the change of the number of neutral or ionised Hydrogen atoms
due to the fluxes is smaller than the number of atoms in each cell.
For instance, considering the r direction and the number density n_HII,
this constraint reads:

    Δ t< c_CFLΔ r/F_r^n_HII.

where c_CFL<1 is a parameter to satisfy the CFL (Courant-Friedrichs-Lewy) 
condition <cit.>. It is
chosen to be c_CFL=0.75 in this work.
We then take the minimum from all the cells, for all directions and for n_HII and n_HI to estimate the integration time Δ t.
However, because of the high heterogeneity in the initial n_HII and n_HI distributions
and the steep dependence of α_HII and q_HI on the temperature,
the explicit solution of the advection term can lead to values of Δ t
too small for reasonable computational times.

In order to circumvent this problem, we apply two measures.
First of all we use a smoothed distribution of ρ, T, and v⃗
as input from the MHD equations, and second we apply a flux limiter for the Godunov scheme.
The smoothing simply consists of replacing the values of each quantity at time t in each computational cell
with the average among a volume of 3×3×3 cells centred around the same cell.
As far as the flux limiter is concerned, we first estimate the Δ t resulting from
a given distribution of n_HII, n_HI, and v⃗.
If Δ t is smaller than 0.1 s in a cell, we redistribute a fraction of 10^-4 of the 
number of ionised and neutral Hydrogen atoms from that cell to the neighbouring cells.
This is iterated until c_CFLΔ t≥0.1 s.
When this condition is satisfied we pick Δ t=0.1 s that is certainly smaller than the Δ t required by Eq.<ref>.
This flux limiter is equivalent to a diffusion term that operates on a faster time scale of the advection of the bulk plasma motion,
such as turbulence.



  §.§.§ Ionisation/Recombination term

The ionisation/recombination term expresses the rate at which
neutral or ionised Hydrogen atoms change ionisation state
per unit of time.
As Hydrogen atoms have only two possible ionised state we can write that:

    [IR_HI]=-[IR_HII]=n_HII(α_HIIn_HII-q_HIn_HI)Δ t.

In order to estimate the value of [IR_HI] and [IR_HII] 
we use an implicit scheme.
In this implicit scheme we set Δ t=0.1, as we do for the advection term, to ensure
that the two terms are integrated over the same time interval and 
then we solve the implicit scheme associated with the finite difference expression for Eq.<ref> for [IR_HI] and the corresponding equation for  [IR_HII].

This corresponds to solving the system of quadratic equations

    n_HI^t+Δ t-n_HI^t  =n_HII^t+Δ t(n_HII^t+Δ tα_HII-n_HI^t+Δ tq_HI)Δ t


    n_HII^t+Δ t-n_HII^t=n_HII^t+Δ t(n_HI^t+Δ tq_HI-n_HII^t+Δ tα_HII)Δ t

for n_HI^t+Δ t and n_HII^t+Δ t that are the values for n_HI and n_HII at the time t+Δ t,
whereas n_HI^t and n_HII^t are the known values of n_HI and n_HII at the time t.


In a few cells near the lower boundary of the MHD simulation where boundary conditions 
lead to sharp gradients in the number density of ionised and neutral Hydrogen atoms
or rapid changes in the plasma temperature occur,
it is possible that n_HII or n_HI could become negative.
For this reason we always impose ionisation equilibrium over 
the first 10 layers of the computational domain in r,
till r<1.17 R_⊙,
and when this happens we 
impose equilibrium ionisation also on other cells to restore a physically meaningful solution.
It should be noted that as soon as the magnetic flux rope is ejected
from its initial position near the lower boundary at t=0,
these locations are no longer relevant for the analysis of the CME ionisation state as the magnetic flux rope is quickly located far from the lower boundary where these issues arise.



§ RESULTS



Having solved the MHD equations for a magnetic flux rope ejection
and then reconstructed the history of the ionisation state of Hydrogen atoms
we now discuss how the abundances of ionised and neutral Hydrogen atoms without
ionisation equilibrium compare to those derived in the usual hypothesis of 
ionisation equilibrium.

Fig.<ref> shows
the column density of n_HI and n_HII
at t=0 and t=19.7 min
from the same point of view as in Fig.<ref>.
In order to display these results in terms of a potential observations we use a 3D cartesian box, where x is the Sun East-West direction,
y is the North-South direction, and z is the direction along the LOS.

The evolution of both n_HI and n_HII follow a similar
pattern as in Fig.<ref>
where the distribution of n_HI and n_HII
change significantly as soon as the magnetic flux rope ejection
propagates through the corona.
The CME three-parts structure is hinted in the distributions of 
n_HI and n_HII.

The key aim of this work is to 
measure the out-of-equilibrium ionisation effects 
in different CME components.
Therefore, in Fig.<ref> we show 
the difference of the logarithms of the ionisation fraction for the column densities of neutral and ionised Hydrogen at t=19.7 min .

In these maps we find positive values when the non-equilibrium ionisation
effects lead to a cumulative higher number of neutral or ionised Hydrogen atoms 
compared to that of the ionisation equilibrium condition
along the LOS. Negative values occur when the number of neutral or ionised atoms is less than in the equilibrium condition.
We find extended regions where non-equilibrium ionisation effects become relevant.
This is particularly true for the abundance of neutral Hydrogen atoms where differences in logarithm between the 
number of atoms and the equilibrium condition along the LOS can be as large as 0.2, whereas
for the ionised Hydrogen atoms such differences are less than than 10^-8
making non-equilibrium ionisation effects rather marginal.
In particular for the abundance of neutral Hydrogen we find 
that overall there are more atoms of this species along the LOS
than prescribed by ionisation equilibrium conditions and only in some small localised regions do we find less atoms of this species with respect to ionisation (Fig.<ref>a).

As non-equilibrium ionisation effects are relevant for neutral Hydrogen in particular, 
we plot the number density of neutral Hydrogen in Fig.<ref>a
at the initial condition, at t=14.5 min, and at t=19.7 min
along the direction of propagation of the magnetic flux rope.
If we focus on the curve at t=19.7 min
we find that the number of neutral Hydrogen atoms
increases with respect to the initial conditions 
where the magnetic flux rope is located (s∼0.7 R_⊙),
as this is a structure colder than the ambient corona travelling upwards, 
thus carrying a higher abundance of 
neutral Hydrogen atoms through advection.
At the same time, the situation is more complex for the CME front (s∼1.6 R_⊙)
where the number of neutral Hydrogen atoms has not changed significantly with respect to the 
pre-CME corona. The number of neutral Hydrogen atoms changes more because of the local temperature variations, i.e. ionisation and recombination processes rather than neutral Hydrogen atoms advection.


Comparing the neutral Hydrogen atoms abundance with the ionisation equilibrium conditions,
the difference is much more evident at the CME front than at the magnetic flux rope location (Fig.<ref>b).
This is clearly visible in the cuts at t=14.5 min and t=19.7 min.
At the magnetic flux rope location
the abundance of neutral Hydrogen atoms can be more than 10 times larger than 
those of ionisation equilibrium, but
it is also highly variable and this occurs in extremely localised regions.
In contrast, at the CME front the abundance of neutral Hydrogen atoms
is almost 100 times larger than that of ionisation equilibrium and this occurs consistently through the CME front for an extension of about ∼0.3 R_⊙.
In Fig.<ref>c we show the change in the neutral Hydrogen relative abundance along the same direction, where we find an increase of neutral Hydrogen relative abundance with respect to the initial configuration at the flux rope location. As said, this is a natural consequence of the advection of the many neutral Hydrogen atoms from the initial position of the flux rope to higher altitude.
A small increase is found also at the CME front locations that is due to local
velocity fields near the front.

In Fig.<ref> we compare the neutral Hydrogen atoms abundance with the ionisation equilibrium conditions
for the centre of the magnetic flux rope location and the CME front as a function of time.
For the flux rope, the abundance of the neutral Hydrogen atoms remains very close to those of the one at ionisation equilibrium condition at all times, except for a very short lived peak at the beginning of the simulation when impulsive heating occurs. This shows that
the high density conditions in the flux rope allows for a quick settlement to ionisation equilibrium.
In contrast, the abundance of neutral Hydrogen atoms 
steadily departs from ionisation equilibrium conditions 
as soon as a CME front develops.
In our study we follow the ionisation state of Hydrogen atoms for approximately 20 minutes 
and at the end the abundance of the neutral Hydrogen atoms saturates
to a value that is about 100 times larger than the abundance expected from
ionisation equilibrium.


The physical explanation for this behaviour can be explained by comparing the
evolution of the number of neutral Hydrogen atoms in the magnetic flux rope and in the CME front 
as a function of time, which is shown in Fig.<ref>.
This is to be compared with Fig.<ref> where we show the evolution of density and temperature in the same structures.

For the centre of the magnetic flux rope the expected density of neutral Hydrogen atoms at ionisation equilibrium 
steadily decreases (blue dashed curve) as the plasma density there decreases faster than the temperature does. 
While a lower temperature would lead to more neutral Hydrogen atoms, the overall lower density 
prescribes less particles in general. 
However, when we reconstruct the ionisation state of the plasma we find less neutral Hydrogen atoms in the 
magnetic flux rope (blue solid line) than at the ionisation equilibrium
as the ionisation state needs to adapt while the temperature is decreasing.
Also for the CME front, the expected density of neutral Hydrogen atoms at ionisation equilibrium
steadily decreases because of the steady decrease of the plasma density (red dashed curve), whereas
the CME front temperature remain more or less constant.
However, the decrease of neutral Hydrogen atoms is not as rapid (red solid curve),
because the front is initially loaded with a higher number of neutral Hydrogen atoms
and as the temperature increases, the ionisation state needs time to adapt to the changing temperature.

Fig.<ref> shows the relative abundance of neutral Hydrogen atoms as a function of the temperature for the flux rope and the CME front.
This is compared with the abundance that would be measured at the ionisation equilibrium. Although this specific estimate is highly dependent on the CME model that is used here, 
it still bears some interesting results, as it happens that the flux rope is mostly near the ionisation equilibrium except in the range of 5.7<Log_10(T)<5.9.
This is the temperature range in which the flux rope shortly after the eruption when the resistive heating occurs and it then converge to equilibrium as it cools down.
In contrast, the front remain close to ionisation equilibrium for lower temperatures Log_10<6.9 and then it significantly departs from equilibrium at higher temperatures. 


It should be noted however that the unavoidable superposition of structures along the LOS can have a significant effect on these estimates.
Fig.<ref> shows the number density of neutral Hydrogen along the two lines of sight identified in Fig.<ref>, 
one cutting through the magnetic flux rope (Fig.<ref>a)
and the other through the CME front (Fig.<ref>b)
and it compares these densities with the values obtained assuming ionisation equilibrium.
The magnetic flux rope is a structure much denser than other structures along the LOS, therefore the measured ionisation state is not changing significantly whether we consider or not the integration along the LOS.
On the other hand, the CME front has a density comparable with the background density, hence the significantly larger amount of neutral Hydrogen atoms that we find at the front location is smoothed when the integration along the LOS is considered.
For this reason, the flux rope shows an abundance of neutral Hydrogen atoms that is on the order of 10 times (or less) the abundance in ionisation equilibrium and the same ratio is shown when we integrate along the LOS. The CME front, instead, shows an abundance of neutral Hydrogen atoms that is up to 100 times larger than the abundance in ionisation equilibrium, but when we integrate along the LOS the count of neutral Hydrogen atoms is only about 10 times larger than the ionisation equilibrium conditions.




§ DISCUSSION AND CONCLUSIONS



In this work we have used the MHD simulation of a magnetic flux rope ejection and CME propagation 
in the solar corona from <cit.> to study the abundances of ionised and neutral Hydrogen atoms 
during CMEs.
In particular, we focus on the non-ionisation equilibrium effects in order to verify to what extend the assumption of
ionisation equilibrium holds during these phenomena.
This work is especially important for the forthcoming measurements of Metis, the coronagraph on board the Solar Orbiter mission, and LST, the coronagraph on-board ASO-S,
that will observe the solar corona in the Lyman-α line that is produced by neutral Hydrogen atoms.
A number of measurements of plasma properties, such as the plasma temperature, crucially depend 
on the amount of neutral Hydrogen atoms and on the assumptions that are been made for the ionisation state of the plasma.
For example, one way to derive the plasma temperature is to compare the abundances of neutral and ionised Hydrogen atoms
and to find the temperature at which the measured abundances are in agreement with the theoretical ones estimated by the ionisation equilibrium hypothesis.
Such kinds of derivations are inherently prone to errors if the plasma is not in ionisation equilibrium.
To this end, this work sheds light on the regions and phenomena  that occur when the assumption of ionisation equilibrium holds and when it does not.

In order to reconstruct the ionisation state of Hydrogen atoms during a CME propagation in the solar corona,
we use the MHD simulation of <cit.> that provides the 3D plasma density and temperature distributions 
as a function of time, together with the plasma velocity fields.
Using these as inputs and assuming ionisation equilibrium at the beginning of the CME propagation,
we devised a technique to solve the differential equations that describe the change in the number of neutral and ionised Hydrogen atoms as a function of space and time, using the strategy already introduced in <cit.>.
More specifically, we separately estimate the advection and ionisation/recombination terms to 
express the rate of change of the number of ionised and neutral Hydrogen atoms in one cell of the MHD simulation.
By time-integrating the rate of change we can describe abundance of neutral and ionised Hydrogen atoms
as a function of space and time in our MHD simulation.
This technique bears no assumption on the ionisation state of the plasma 
and it can describe abundances of ionised and neutral Hydrogen atoms out of ionisation equilibrium.

Finally, to analyse the results of this study with regards to the temperature measurements in the solar corona,
we have compared the obtained distributions of neutral and ionised Hydrogen atoms
with respect to the distributions associated with ionisation equilibrium.
We do this through  the entire 3D domain of the MHD simulation and specifically for 
the centre of the magnetic flux rope and the CME front
that are two prominent features of CMEs.

We find that the departures from ionisation equilibrium are negligible for the abundance of ionised Hydrogen atoms, where any
departure is small enough to become irrelevant, as the vast majority of Hydrogen atoms are ionised in the solar corona.
However, this is not the case for  neutral Hydrogen since there are inherently few 
numbers of these atoms in the corona 
and therefore small variations can lead to a significant departure from the assumed abundance at ionisation equilibrium.
We crucially find that non-equilibrium ionisation effects for neutral Hydrogen atoms are marginal for the regions pertinent to the magnetic flux rope.
This means that the amount of neutral Hydrogen atoms in the magnetic flux rope is generally consistent with the number 
one would assume in ionisation equilibrium conditions.
Although we observe some small scale variations and rapid changes from one location to another within the magnetic flux rope,
overall we find that these effects are not detectable when the whole flux rope is considered.
This means that, for the future derivation of plasma temperatures based on observed Lyman-α line
intensities, the ionisation equilibrium hypothesis is still applicable.

In contrast, we find that non-equilibrium ionisation effects are relevant at the CME front,
where the number of neutral Hydrogen is much larger than the expected abundance from the assumption of ionisation equilibrium conditions.
We find this to be consistent over the whole extension of the CME front and that
the departure is as large as 2 orders of magnitude in excess with respect to the ionisation equilibrium conditions. In the region where we have an excess of neutral Hydrogen atoms, the derived temperature from the abundances would be much lower than 
the actual temperature value, and thus using the ionisation equilibrium assumption would lead to a significant plasma temperature underestimate at the CME front.
At the same time, we also find that the projection effect along the LOS can partially hide these effects as it happens here for the CME front where non-equilbrium ionisation effects are smoothed out after the integration along the LOS and the superposition of structures in ionisation equilibirium that have comparable densities.

In the future more studies need to be carried out, to investigate the parameter space of the CME velocities, densities and temperatures, along with considering more diverse configurations of the magnetic flux rope ejections. It is also important to apply this approach (the magnetofrictional simulation, the MHD simulation, and the reconstruction of the ionisation state of the Hydrogen atoms) 
to observed events to more accurately verify these results.
In particular, the analysis carried out in Fig.<ref> can potentially be extended to many more cases where MHD simulations can be used to derive a more accurate and reliable lookup table to invert the temperature of CME fronts and cores as a function of the relative abundance of neutral Hydrogen.
At the same time, this work already provides key information for a correct interpretation of future plasma properties derived from Lyman-α observations, which is
certainly crucial for the future observations of Metis on board Solar Orbiter and LST on board ASO-S telescopes.



This research has received funding from the Science and Technology Facilities Council (UK) through the consolidated grant ST/N000609/1 and the European Research Council (ERC) under the European Union Horizon 2020 research and innovation program (grant agreement No. 647214).
D.H.M. would like to thank both the UK STFC and the ERC (Synergy grant: WHOLE SUN, grant Agreement No. 810218) for financial support. D.H.M. and P.P. would like to thank STFC for IAA funding under grant number SMC1-XAS012.
This work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk. The equipment was funded by BEIS capital funding via STFC capital grants ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham University and STFC operations grant ST/R000832/1. DiRAC is part of the National e-Infrastructure.
We acknowledge the use of the open source (gitorious.org/amrvac) MPI-AMRVAC software, relying on coding efforts from C. Xia, O. Porth, R. Keppens.
CHIANTI is a collaborative project involving George Mason University, the University of Michigan (USA), University of Cambridge (UK) and NASA Goddard Space Flight Center (USA).



aa


