\section{Introduction}\label{sec:intro}
Solar flares, and other transient energy release events, can dramatically disturb the solar atmosphere, driving it out of equilibrium. The importance of non-equilibrium (NEQ) ionisation and recombination in determining the atomic level populations (and therefore the radiative response to the disturbance) will likely vary from species to species, and the charge states of those species. However, when modelling the radiative response of the atmosphere to flares it is often the case that the statistical equilibrium (SE) solution is sought, given the computational demands of performing time dependent simulations that include physical processes, such as partial frequency redistribution, that are required for certain transitions \citep[recent examples include][]{2017ApJ...842...82R,2019ApJ...879...19Z,kerr_2019}. This is partly mitigated by using the non-equilibrium electron density from dynamic simulations. We address the importance of non-equilibrium ionisation on the modelling of Mg~\textsc{ii} h \& k lines during solar flares, a routinely observed pair of strong spectral lines from the Sun.  

In the standard flare model, magnetic reconnection occurs in the corona, releasing vast amounts of magnetic energy, up to $10^{32}$~ergs, \citep{2011SSRv..159...19F} . This results in \textsl{in situ} heating and particle acceleration. Electrons (and likely ions) are accelerated at relativistic speeds along flare loops, loosing energy via Coulomb interactions as they travel into denser regions \citep[][]{1971SoPh...18..489B,2011SSRv..159..107H,2019ApJ...880..136J}. A broadband enhancement to the solar radiative output results from both thermal and non-thermal processes following plasma heating and ionisation. Chromospheric ablations (`evaporations'), with bulk upflows reaching several hundreds of km s$^{-1}$, and condensations, with bulk downflows of dense material reaching a few tens of km s$^{-1}$, are also driven during flares, revealed by Doppler shifts of spectral lines \citep[e.g.][]{1985ApJ...289..414F,1989ApJ...346.1019F,2009ApJ...699..968M,2015ApJ...807L..22G}.

The transition region and chromosphere are both the sites of energy deposition and the origin of the bulk of the radiative output during flares \citep{2011SSRv..159...19F,2014ApJ...793...70M}. They are crucial locations for both diagnosing the flaring plasma, and as test grounds for models of energy transport during flares. For the former, forward modelling radiation from different flare atmospheres allows us to understand how to extract physical properties of the flaring plasma from observations. For the latter, model-data comparisons allows us to understand how well simulated flare atmospheres driven by a particular energy transport mechanism compare with the actual flaring chromosphere, and therefore how consistent our models of energy transport are with the reality.

Since the launch of the Interface Region Imaging Spectrograph \citep[][IRIS]{2014SoPh..289.2733D}, observations of the Mg~\textsc{ii} h \& k resonance and subordinate lines have become commonly used in observational studies \citep[e.g.][]{2015A&A...582A..50K,2015SoPh..290.3525L,2015ApJ...807L..22G,2018PASJ...70..100T,2018ApJ...856...34T,2018ApJ...861...62P}, and are attractive as a means to diagnose or critically attack flare (and other) models \citep[e.g][]{2016ApJ...827...38R,2016ApJ...827..101K,2019ApJ...879...19Z}. These lines have been observed in hundreds of flare events since the launch of IRIS, show significant variation from the quiet Sun during flare heating, and form over a range of chromospheric layers meaning they have the potential to diagnose a large part of the lower solar atmosphere. Particularly if they are combined with other IRIS or ground based observables. 

These lines are optically thick, and as such are complex to interpret, requiring forward modelling to understand their formation \citep{2013ApJ...772...89L}. Modelling of the Mg~\textsc{ii} NUV spectrum has been a key complementary activity to the analysis of observations, and it is essential that we are modelling these lines in an accurate manner. This is both to provide an accurate physical interpretation of these complex optically thick lines, but also to facilitate high fidelity model-data comparison by which models of flare energy transport are critically interrogated. 

In Paper 1 \citep{kerr_2019} we used \texttt{RADYN} \citep{1992ApJ...397L..59C,1997ApJ...481..500C,2015ApJ...809..104A} radiation hydrodynamic flare atmospheres with the radiation transport code \texttt{RH}, \citep{2001ApJ...557..389U}, a common approach to modelling these lines in flares. This was to investigate how various setups affect the Mg~\textsc{ii} solution such as partial frequency redistribution (PRD) vs complete frequency redistribution (CRD), inclusion of other species, and effects of coronal irradiation.  In particular we noted that PRD is still required, resulting in substantial differences to the CRD solution. However, in that work the level populations were obtained under the assumption of statistical equilibrium, since RH is a stationary code. Each atmospheric snapshot was treated in isolation, and the statistical equilibrium solution obtained. This was necessary in order to include the more advanced radiation transfer from the RH code (in particular PRD), that would be computationally very demanding to include in a time-dependent flare simulation. The history of the atmosphere and any time dependent effects were neglected. In this paper we will investigate the impact of non-equilibrium processes on the formation of Mg~\textsc{ii} during solar flares, with the aim of determining if the usual practice of omitting these effects is safe. 

The chromosphere is dynamic, especially during solar flares. \cite{1999AIPC..471...23C,2002ApJ...572..626C} demonstrated, from \texttt{RADYN} simulations of propagating acoustic waves, that the ionisation and recombination timescale ($\tau\sim10^3-10^5$~s) for hydrogen is long compared to the dynamical timescale, and that if statistical equilibrium is assumed then the ionisation fraction is underestimated by several orders of magnitude in certain locations of the atmosphere. The electron density would consequently be very different.  Similarly, the requirements for considering non-equilibrium effects on other species have been noted, such as Helium \citep{2014ApJ...784...30G,2016ApJ...817..125G}, O~\textsc{iv} \citep{2013ApJ...767...43O}, and Si~\textsc{iv} \citep{2016ApJ...817...46M}. \cite{2013ApJ...772...89L} investigated if non-equilibrium processes were important for Mg~\textsc{ii} in the quiet Sun, concluding that whenever the temperature was large enough to produce significant amounts of Mg~\textsc{iii}, the relaxation time was shortened, and that using statistical equilibrium was sufficient. Of course flare chromospheres are very dynamic and the conclusions of \cite{2013ApJ...772...89L} might not apply in flaring conditions. We investigate the flaring scenario in this work. 

The non-equilibrium atomic level population equation is: 

\begin{equation}\label{eq:nequil}
	\frac{\partial n_{i}}{\partial t} + \frac{\partial n_{i}v}{\partial z} - \left(\sum^{N^\prime}_{j\ne i} n_{j}P_{j,i} - n_{i}\sum^{N^{\prime}}_{i\ne j} P_{i,j} \right) = 0,
\end{equation}

\noindent where $v$ is the atmospheric velocity, $N^{\prime}$ is the total number of states, $t$ is time, $z$ is the height in the atmosphere, $P_{i,j}$ describes the total rates (collisional plus radiative) from $i$ to $j$, and $P_{j,i}$ is the total rate from $j$ to $i$. The transition rates are functions of the local atmospheric conditions (including energy input) which vary with time in dynamical simulations. If the local thermodynamic state of the atmosphere or the radiation field vary faster than the timescale for ionisation and recombination then there is not enough time for the atmosphere to reach equilibrium - the `history' of the atmosphere becomes important. 

If the ionisation and recombination timescales are sufficiently fast then the populations can be approximated by statistical equilibrium (setting $\partial n_{i}/\partial t$ and $\partial n_{i}v/\partial z$ to zero in Equation~\ref{eq:nequil}):

\begin{equation}\label{eq:statequil}
	 \sum^{N^\prime}_{j\ne i} n_{j}P_{j,i}  =  n_{i}\sum^{N^{\prime}}_{i\ne j} P_{i,j} ,
\end{equation}

\noindent In situations where the history of the atmosphere is important but statistical equilibrium is used, there may consequently be errors in the population densities of atomic states, and in the synthetic spectra. 



\section{Flare Simulations}\label{sec:codes}
\begin{figure*}
	\centering 
	{\includegraphics[width = .95\textwidth, clip = true, trim = 0.cm 0.cm 0.cm 0.cm]{Flare_atmospheres.pdf}}	
	\caption{\textsl{Stratification of temperature (a,b,c), electron density (d,e,f), and macroscopic velocity (g,h,i; uplows are negative) in the three flare simulations. The first column shows the F9 simulation, second column shows F10 simulation, and third column shows the F11 simulation. Colour represents time. Recall that heating ceased at $t=10$~s}. Note that this is a reproduction of Figure~1 in \cite{kerr_2019}.}
	\label{fig:flare_atmos}
\end{figure*}
Three flare simulations were produced using the \texttt{RADYN} radiation hydrodynamic code \citep{1992ApJ...397L..59C,1997ApJ...481..500C,1999ApJ...521..906A,2005ApJ...630..573A,2015ApJ...809..104A,FP_inprep}, with the flares driven by non-thermal electron beams that were injected into a pre-flare atmosphere that spanned the sub-photosphere through to the corona. Both the code and details of these simulations are described in Paper 1 \citep{kerr_2019}. The electron beam fluxes were $F=[1\times10^{9}, 1\times10^{10}, 1\times10^{11}]$~erg~cm$^{-2}$~s$^{-1}$ (hereafter F9, F10, \& F11) which were modelled as power law distributions with spectra index $\delta = 5$ and a low-energy cutoff $E_{c} = 20$~keV. Energy was injected at a constant rate for $t=10$~s, and the simulations allowed to continue to evolve until $t=50$~s. Figure~\ref{fig:flare_atmos} shows the evolution of the atmospheres as a function of time in each simulations. The grey shaded portion is the pre-flare atmosphere. 

\texttt{RADYN} simulated the time-dependent response of the atmospheres to flare energy injection so that the atomic level populations, electron density, and temperature stratification are all NEQ, but with radiation computed using the simplifying assumption of CRD (necessary to make the dynamic flare problem computationally tractable). Species important for energy balance are considered in the main simulation (H, He, and Ca~\textsc{ii}), with other species included in an optically thin radiation loss function. While Mg~\textsc{ii} h \& k are very strong chromospheric lines, PRD is required to accurately model the radiative losses. Omitting the Mg~\textsc{ii} h \& k lines is likely safe to do from an energetic balance standpoint, given that the Ca~\textsc{ii}  H \& K lines \textsl{are} included, but in CRD. So, losses from Ca~\textsc{ii} H \& K are overestimated, but this is mitigated by ignoring Mg~\textsc{ii} h \& k.

In order to obtain the Mg~\textsc{ii} NEQ atomic level populations and emergent profiles we must therefore turn to \texttt{MS\_RADYN}, the minority species version of \texttt{RADYN}. This code uses each internal timestep of an existing \texttt{RADYN} RHD solution to solve the NEQ NLTE radiation transport problem for a desired minority species, such as Mg~\textsc{ii}. It includes the time dependent and advection terms when solving the atomic level populations. Since the variables required for this simulation are stored for each internal timestep (and not the output cadence of the main set of variables, which can be many times larger) \texttt{MS\_RADYN} can capture changes to the atomic level populations on very short timescales that result from the changing atmospheric state.  \texttt{MS\_RADYN} was used in this fashion to study C~\textsc{ii} emission \citep{2003ApJ...597.1158J} and more recently to investigate radiative transfer effects on Si~\textsc{iv} during flares \citep{2019ApJ...871...23K}.

Our three flare simulations were used as input in \texttt{MS\_RADYN} to simulate the atomic level populations and synthetic spectra of a ten-level-with-continuum model of Mg~\textsc{ii}, the same as used in \cite{2013ApJ...772...89L} and that we used in \cite{kerr_2019}. This model atom included the h \& k transitions, and the subordinate line transitions. To test the impact of NEQ effects we then repeated these simulations, but switched off the time-dependent and advection terms when computing the atomic level populations so that \texttt{MS\_RADYN} then used statistical equilibrium (Equation~\ref{eq:statequil}). In this latter series of simulations, the Mg~\textsc{ii} problem was solved using the non-equilibrium hydrogen and electron densities (similar to the more typical post-processing of RHD/HD snapshots through \texttt{RH}). 

Using \texttt{MS\_RADYN} in this fashion meant that the only differences were the terms included in the atomic level population equation (Equation~\ref{eq:nequil} vs \ref{eq:statequil}), with the same flare atmospheres and background opacities used throughout. Any differences in the level populations and emergent spectra are then due to non-equilibrium processes. 



\section{Line Profiles}\label{sec:emergent_intensity}

\begin{figure}
	\centering 
	{\includegraphics[width = 0.5\textwidth, clip = true, trim = 0.cm 0.cm 0.cm 0.cm]{MgII_kline_NEQ_stackplot_F11.pdf}}	
	\caption{\textsl{Mg~\textsc{ii} k line profiles in the F11 simulation. In panels (a,b) wavelength is shown on the x-axis, and time is stacked on the y-axis, so that the temporal evolution is shown. Panel (a) shows the NEQ profiles, and panel (b) the percentage change between the NEQ and SE profiles (positive means that SE is more intense, negative that NEQ is more intense, and the change is saturated on the scale indicated). Panels (c-g) show a comparison of the NEQ (black lines) to SE (red dashed lines) line profiles at selected times, where the grey lines are the percentage change.}}
	\label{fig:f11_kline_stack}
\end{figure}
\begin{figure}
	\centering 
	{\includegraphics[width = 0.5\textwidth, clip = true, trim = 0.cm 0.cm 0.cm 0.cm]{MgII_kline_NEQ_stackplot_F10.pdf}}	
	\caption{\textsl{Same as Figure~\ref{fig:f11_kline_stack}, but for the F10 simulation.}}
	\label{fig:f10_kline_stack}
\end{figure}
\begin{figure}
	\centering 
	{\includegraphics[width = 0.5\textwidth, clip = true, trim = 0.cm 0.cm 0.cm 0.cm]{MgII_kline_NEQ_stackplot_F09.pdf}}	
	\caption{\textsl{\textsl{Same as Figure~\ref{fig:f11_kline_stack}, but for the F9 simulation}}}
	\label{fig:f09_kline_stack}
\end{figure}
\begin{figure}
	\centering 
	{\includegraphics[width = 0.5\textwidth, clip = true, trim = 0.cm 0.cm 0.cm 0.cm]{MgII_sline_NEQ_stackplot_F09.pdf}}	
	\caption{\textsl{\textsl{Same as Figure~\ref{fig:f11_kline_stack}, but for the F9 simulation, and here we show the Mg~\textsc{ii} 2791\AA\ subordinate line.}}}
	\label{fig:f09_2791_stack}
\end{figure}

How the emergent profiles differ between each simulation is the most important consideration indicating if there is a negligible, tolerable, or significant difference between NEQ and SE when modelling Mg~\textsc{ii}, given our goal of comparing to observations. Here we discuss these differences in each scenario, but we do not address in detail the line features themselves and how they vary with flare strength or how well they compare to observations, which will be the focus of other investigations. 
 
To quantify the difference between the two results we compute the percentage change between the intensities across the resonance and subrodinate line profiles, $\frac{I_{\mathrm{SE}} - I_{\mathrm{NEQ}}}{I_{\mathrm{NEQ}}}\times100$, where $I_{\mathrm{SE}}$ is the intensity computed using SE, and $I_{\mathrm{NEQ}}$ is the intensity computed using NEQ. A positive percentage means that $I_{\mathrm{SE}}$ is more intense.

For the bulk of the duration of the flares the differences between the results is only minor or negligible, though there are exceptions. Figures~\ref{fig:f11_kline_stack},~\ref{fig:f10_kline_stack}, and \ref{fig:f09_kline_stack} show a comparison between the NEQ and SE Mg~\textsc{ii} k line profiles, and the percentage change between the two solutions, for the F11, F10, and F9 simulations respectively. The Mg~\textsc{ii} 2791~\AA\ line is shown for the F9 simulation in Figure~\ref{fig:f09_2791_stack}. In each figure panel (a) shows the line profile as a function of time (stacked on the y-axis), and panel (b) shows the percentage change. Panels (c-g) show profiles at selected times of interest. 

It is immediately clear from these figures that while assuming SE results in very little change during the main heating phase, it does result in under- or over-estimations of the line intensity across the whole k line during the initial energy injection and following cessation of the beam. The subordinate line differences are more limited to the line core and near wings. Line shapes and features are preserved, though in the decay phase the subordinate lines return from emission to absorption more rapidly in the NEQ simulations. 

Within the first second of the F11 and F10 simulations, there is an intensity change $-[5-50]\%$, with the strongest changes around the emission peaks. These differences rapidly subside in both F11 and F10, until the cessation of energy input after which a very narrow region around the line core shows approximately a $10-20\%$ change, again with NEQ more intense. In the F10 simulation there are several seconds following the cessation of energy input where the NEQ solution is significantly less intense, up to approximately $230\%$ in the blue emission peak. The F9 simulation evolves more gradually, and shows close agreement between the NEQ and SE solutions in the initial stages. Unlike the other two simulations, however, differences appear during the heating phase, with a larger NEQ result, with a difference on the order $5-20\%$. In the cooling phase the F9 behaves similar to the F10 simulation, with a stronger SE intensity. The h \& k lines are not affected by the same magnitude meaning the k:h line ratio will vary between the NEQ and SE solutions, as discussed in Section~\ref{sec:line_ratios}. 

The largest differences appear around the emission peaks in the decay phase, with the wings showing less significant changes. Integrating across the line profiles therefore reduces the discrepancy at these times substantially. The discrepancies in the initial heating phase for the F10 and F11 simulations affected the line wings also, so the percentage change is of a similar magnitude to the specific intensity case at those times. Figures~\ref{fig:lcurves}(a,b,c) show the lightcurves of the Mg~\textsc{ii} k line, integrated $\pm 0.5$~\AA\ around the line core. 

\begin{figure}
	\centering 
	\hbox{
		\subfloat{\includegraphics[width = 0.5\textwidth, clip = true, trim = 0.0cm 0.15cm 0cm 0cm]{MgII_NEQ_Lcurve_F11.pdf}}	
		 }
	\hbox{
	        \subfloat{\includegraphics[width = 0.5\textwidth, clip = true, trim = 0.0cm 0.15cm 0cm 0cm]{MgII_NEQ_Lcurve_F10.pdf}}	
                   }
        \hbox{
	        \subfloat{\includegraphics[width = 0.5\textwidth, clip = true, trim = 0.0cm 0.15cm 0cm 0cm]{MgII_NEQ_Lcurve_F09.pdf}}	
                   }
	\caption{\textsl{Lightcurves of the Mg~\textsc{ii} k line $\pm0.5$~{\AA} in each simulation (a) F11, (b) F10, and (c) F09. The h line and subordinate lines behave qualitatively similar. An inset in panel (a) shows $t < 1$~s in more detail. In each panel: the black line is the k line NEQ solution, the red is the k line SE solution, the yellow line is the $2791$\AA\ NEQ solution, the blue dashed line is the 2791\AA\ SE solution, the grey dot-dashed line is the k line percentage change, and the purple dot-dashed line is the 2791\AA\ percentage change.}}
	\label{fig:lcurves}
\end{figure}

In the moderate-to-strong flares the initial heating ($t<1$~s) shows a local maximum in the NEQ solution, but not in the SE solution. In the F10 flare, the NEQ solution shows a steeper decrease following the cessation of the beam. At other times the lightcurves generally show similar behaviour, with small intensity differences (note the logarithmic scales, so that while the lines do appear in close agreement they can actually differ by a few tens of percent). 


\section{Line Ratios \& Formation Heights}\label{sec:line_ratios}
\begin{figure}
	\centering 
	{\includegraphics[width = 0.5\textwidth, clip = true, trim = 0.cm 0.cm 0.cm 0.cm]{MgII_NEQ_ratio}}	
	\caption{\textsl{Ratio of the k to h line ($R_{k:h}$) in each simulation. The red lines are the SE solution, black are the NEQ solution and grey dot-dashed lines shows the percentage change.}}
	\label{fig:ratios}
\end{figure}

The k:h line ratio, $R_{k:h}$, is a useful metric that can indicate if the lines are optically thick, if there are relative changes between the lines during the flare, or if radiative excitation processes are significant \citep[e.g.][]{2014ApJ...792...93H}. \cite{2015A&A...582A..50K} noted that in some areas of a flare ribbon observed by IRIS, $R_{k:h}$ changed in response to the flare. The ratio decreased slightly, and the relative spread in values decreased (there was a tighter correlation of $R_{k:h}$ than in areas outside of the ribbon). They speculated that this could be caused by the h \& k lines forming closer together than in the quiet Sun, and/or sampling a more uniform chromosphere during the flare. In the optically thin limit this ratio is $R_{k:h} = 2$, the ratio of the statistical weights of the k \& h upper levels. It is typical that $R_{k:h} \approx 1.2$, both in the quiet Sun and in flares \citep[e.g.][]{2015A&A...582A..50K,2015SoPh..290.3525L}, indicating optically thick line formation. Note, though, that $R_{k:h}$ can theoretically have a value of two even in the optically thick case as the ratio of the source functions can in effect take \textsl{any} value \citep{2015ApJ...811...80R}. 

It is important to determine if this ratio is affected by non-equilibrium effects that will confuse the interpretation of $R_{k:h}$ variation in flares if forward modelled using statistical equilibrium. 

Figure~\ref{fig:ratios} shows $R_{k:h}$ in each simulation where again we show the percentage change between NEQ and SE. The magnitude of the variations between the NEQ and SE solutions differs for the h \& k lines, meaning that $R_{k:h}$ consequently shows differences over time. The magnitude of this difference is relatively small, on the order of $<10$~\%. Generally the temporal profile of $R_{k:h}$ is preserved, though in the weaker flares the rate of change at the start of the decay phase is smaller in the NEQ solution, and in the F10 \& F11 simulations the NEQ ratios decrease somewhat whereas the SE ratios increase within the first second. Both NEQ and SE solutions have $R_{k:h}< 2$, indicating optically thick line formation, which we confirmed from inspecting the detailed line formation properties. 

\begin{figure}
	\centering 
	\hbox{
        \subfloat{\includegraphics[width = 0.5\textwidth, clip = true, trim = 0.cm 0.cm 0cm 0cm]{MgII_NEQ_F11_zform_Comps.pdf}}
        }
        \hbox{
        \subfloat{\includegraphics[width = 0.5\textwidth, clip = true, trim = 0.cm 0.cm 0cm 0cm]{MgII_NEQ_F11_zdiff_Comps.pdf}}
        }
        \hbox{
        \subfloat{\includegraphics[width = 0.5\textwidth, clip = true, trim = 0.cm 0.cm 0cm 0cm]{MgII_NEQ_F11_Sfn_Comps.pdf}}	
        }
	\caption{\textsl{Panel (a) shows the formation height of the Mg~\textsc{ii}  h line (orange),  Mg~\textsc{ii} k line (black) and  Mg~\textsc{ii} 2791\AA\ line (blue) as a function of time in the F11 NEQ simulation. Panel (b) shows the formation height difference resulting from using the NEQ or SE solution, with the inset highlighting the first $t=0-1$~s of the simulation (note the change in scale to meters). Panel (c) shows the Mg~\textsc{ii} h line source function at $t = 0.1$~s, where black is the NEQ solution and green is the SE solution.}}
	\label{fig:f11lineform}
\end{figure}

\begin{figure}
	\centering 
	\hbox{
        \subfloat{\includegraphics[width = 0.5\textwidth, clip = true, trim = 0.cm 0.cm 0cm 0cm]{MgII_NEQ_F10_zform_Comps.pdf}}
        }
        \hbox{
        \subfloat{\includegraphics[width = 0.5\textwidth, clip = true, trim = 0.cm 0.cm 0cm 0cm]{MgII_NEQ_F10_zdiff_Comps.pdf}}
        }
        \hbox{
        \subfloat{\includegraphics[width = 0.5\textwidth, clip = true, trim = 0.cm 0.cm 0cm 0cm]{MgII_NEQ_F10_Sfn_Comps.pdf}}	
        }
	\caption{\textsl{Same as Figure~\ref{fig:f11lineform}, but for the F10 simulation. Note also the change in scale in panel (b) compared to Figure~\ref{fig:f11lineform}.}}
	\label{fig:f10lineform}
\end{figure}

The ratio is influenced by the relative separation of the h \& k lines in the flare atmosphere, and the formation height of the lines is also useful in relating model results to observations of the flaring plasma. Figures~\ref{fig:f11lineform} \& ~\ref{fig:f10lineform} show  the formation heights of the h, k and 2791\AA\ lines in the NEQ solution (panels a), the formation height differences of the lines ($\Delta z$) between the NEQ and SE solutions (panels b), and the h line source function in the NEQ and SE solutions (panel c), for the F11 and F10 simulations respectively. Note that the scales vary between the two flares. The formation height here is defined as the height at which the $\tau_{\nu} = 1$ surface is maximal (that is, we are defining the line core to the part of the line forming highest in the atmosphere, with the greatest opacity). 

In both cases during the main phase of the flare there is little difference in the formation heights between the NEQ or SE solutions. However, at the very start of the heating phase in both flares the formation heights can differ. In the F11 simulation this is on the order of kilometres, and in the F10 simulation this is on the order tens to a few hundred km. The $\Delta z$ of k and h lines are of different magnitudes, and can be of different directions in the F11 simulation (the k line forms deeper in the NEQ solution whereas the h line forms higher, compared to the SE case). In both cases the subordinate line at 2791\AA\ actually shows a much greater formation height variation than the h \& k lines, forming a few hundred km deeper in the NEQ solution in the F10 flare during the initial decay phase. The h \& k lines in that flare do reach $\Delta z\sim100$~km, but this rapidly decreases to $\Delta z \sim5$~km or smaller.

Given the generally small differences in the ratios and formation heights, and their short lived nature, we do not envisage using SE will result in any significant misinterpretations when relating line profile features to plasma properties. A possible exception is the subordinate line during the decay phase in weak or moderate flare simulations, or if extremely strong gradients are present near the formation heights of the Mg~\textsc{ii} lines. 

\section{Ion Fractions \& Relaxation Timescales}\label{sec:levelpops}
\begin{figure*}
	\centering 
	{\includegraphics[width = 0.9\textwidth, clip = true, trim = 0.cm 0.cm 0.cm 0.cm]{MgII_NEQ_ionfrac_comps}}	
	\caption{\textsl{Mg~\textsc{ii} ion fraction in each simulation at three different times. The top row is the initial heating phase ($t=0.3$~s), the middle row is the main heating phase ($t=6$~s), and the bottom row is the initial cooling phase ($t=10.3$~s). The black lines are the NEQ solutions, the red dashed line are the SE solutions, and the grey dot-dashed lines are the percentage change.}}
	\label{fig:ionfrac_comps}
\end{figure*}

The differences in the line profile results obtained from the NEQ and SE solutions discussed previously can be understood by studying the ion fraction stratification and the atomic level populations. 

Figures~\ref{fig:ionfrac_comps}(a,b,c) show the fraction Mg~\textsc{ii}/Mg, $\chi_{MgII}$, at the very start of the heating phase of the three simulations. In the two stronger flares, where the flare disturbs the atmosphere in a more impulsive and dramatic manner, $\chi_{MgII}$ is larger in the NEQ solution than the SE solution, and consequently there is a larger population in the resonance and subordinate line upper levels. This then results in more radiative decays, and more intense lines. In the F11 simulation there is a narrow region where there is up to an $80~\%$ change in $\chi_{MgII}$ between SE and NEQ. The discrepancy rapidly reduces, both in magnitude and spatial extent.  The F10 behaves in a similar manner. The weak flare, however, shows a difference only through the TR, slightly above the formation height of Mg~\textsc{ii}, so that the line intensity is not really affected.

During the main heating phase the NEQ and SE solutions give largely similar results with only small percentage changes in the F11 and F10 solutions.  Figures~\ref{fig:ionfrac_comps}(d.e,f) shows $t=6$~s, illustrating that in the main heating phase of the flares the NEQ and SE solutions have only marginal differences. In the F9 simulation, which evolved more slowly, differences have started to appear, explaining why there are line intensity differences at this stage in the flare. 

Finally, Figures~\ref{fig:ionfrac_comps}(g,h,i) shows the start of the cooling phase, shortly after the beam has ceased depositing energy, where the atmospheric temperature rapidly drops. Here the situation is reversed, with the SE solution predicting more $\chi_{MgII/Mg}$. The F11 simulation, in which the flare induced electron density and temperature enhancements were much larger, shows a smaller difference in comparison to the F10 and F9 simulations.

The ion fraction differences between the NEQ and SE solutions are due to the ionisation/recombination timescales in each simulation, which vary with atmospheric state (temperature and electron density).

We determined the timescales for the ionisation equilibrium of Mg~\textsc{ii}/Mg~\textsc{iii} following the methodology of \cite{2002ApJ...572..626C} and \cite{2013ApJ...772...89L}: For a selected time step of a flare simulation, the temperature was increased by $1\%$ throughout the atmosphere and the rate equations for the full magnesium model atom were solved as a function of time (keeping the hydrodynamic state constant at the perturbed value), following the relaxation of the population densities from the initial state towards the new equilibrium. The relaxation time scale was calculated from a fit of the time evolution of the population density of Mg~\textsc{iii} to the analytic solution for a two-level atom:

\begin{equation}\label{eq:relaxtime}
n(t)=n(\infty)+(n(0)-n(\infty))e^{-t/\tau_{\mathrm{relax}}},
\end{equation}

\noindent where $n(t)$ is the population density of Mg~\textsc{iii} at time t, $n(\infty)$ is the statistical equilibrium population density in the perturbed atmosphere, $n(0)$ is the equilibrium population density in the initial atmosphere, and $\tau_{\mathrm{relax}}$ is the relaxation timescale.

\begin{figure*}
	\centering 
	{\includegraphics[width = 0.75\textwidth, clip = true, trim = 0.cm 0.cm 0.cm 0.cm]{MgII_Relaxation_Timescale_F11.pdf}}	
	\caption{\textsl{\textsl{In each panel the relaxation timescale (a) and temperature (b) are shown at various snapshots in the F11 simulation. Several times are highlighted by coloured lines, as indicated in the legend.}}}
	\label{fig:f11_timescale}
\end{figure*}
\begin{figure*}
	\centering 
	{\includegraphics[width = 0.75\textwidth, clip = true, trim = 0.cm 0.cm 0.cm 0.cm]{MgII_Relaxation_Timescale_F10.pdf}}	
	\caption{\textsl{\textsl{Same as Figure~\ref{fig:f11_timescale}, but for the F10 simulation.}}}
	\label{fig:f10_timescale}
\end{figure*}
\begin{figure*}
	\centering 
	{\includegraphics[width = 0.75\textwidth, clip = true, trim = 0.cm 0.cm 0.cm 0.cm]{MgII_Relaxation_Timescale_F09.pdf}}	
	\caption{\textsl{\textsl{Same as Figure~\ref{fig:f11_timescale}, but for the F09 simulation.}}}
	\label{fig:f09_timescale}
\end{figure*}

Figures~\ref{fig:f11_timescale},~\ref{fig:f10_timescale}, \& ~\ref{fig:f09_timescale} show the relaxation timescales (top panels) and temperatures (bottom panels) at various times in the three simulations. The relaxation timescales are initially on the order of $\tau_{\mathrm{relax}} \sim 1-2$~s in the upper chromosphere. Through the first second of flare heating in the F11 and F10 simulations, the relaxation time decreases to the order of $\tau_{\mathrm{relax}} < 0.1$s in response to the substantial temperature and electron density increase in the chromosphere. NEQ effects are prominent within this initial second of heating, since even though the relaxation timescale is decreasing, the atmosphere evolves very impulsively. In the main heating phase the atmosphere evolves somewhat more slowly, so that Mg~\textsc{ii} ionisation equilibrium keeps pace given the small  $\tau_{\mathrm{relax}}$.

In the F9 simulation  $\tau_{\mathrm{relax}} = 1$ does not decrease as much as in the stronger flares, but the atmosphere also evolves more gradually so that NEQ effects only become apparent when the dynamics become faster. At these times  $0.1 < \tau_{\mathrm{relax}} < 1$~s, and Mg~\textsc{ii} is somewhat out of equilibrium.

During the decay phase of the F9 and F10 simulations  $\tau_{\mathrm{relax}}$ rises back to the order of 1s when the temperature drops, so that there is initially a large disagreement between NEQ and SE while the atmosphere is cools rapidly. After the initial sharp decrease in temperature (and electron density), the rate of change of the atmosphere where Mg II forms is smaller so that the differences between NEQ and SE get smaller (but are still present to some degree), and Mg~\textsc{ii} is only partially out of equilibrium.
 
In the F11 case the temperature is so much larger at the end of the heating phase that even through the decay phase the relaxation timescale is still sufficiently small that the discrepancies are reduced.
 


\section{Summary \& Conclusions}\label{sec:conclusions}
Non-equilibrium effects on the formation of Mg~\textsc{ii} spectra during solar flares of different magnitudes has been investigated using \texttt{RADYN} and \texttt{MS\_RADYN} simulations. The time-dependent NLTE NEQ atomic level populations and synthetic spectra were computed, and compared to the NLTE SE solution.

While line profile shapes are preserved in each solution, the intensities of lines can differ, sometimes substantially, if NEQ effects are taken into consideration. These largely appear in the initial heating phase ($t<1$~s), and in the initial seconds of the decay phase. Investigation of Mg~\textsc{ii} ionisation equilibrium showed that changes in the atmospheric state can lower the relaxation timescale, meaning that Mg~\textsc{ii} can be very close to equilibrium during the main heating phase, but that when the temperature rapidly falls Mg~\textsc{ii} is again driven out of equilibrium. 

In \cite{kerr_2019} we found that partial frequency redistribution is required to accurately forward model Mg~\textsc{ii} in flares, with intensity changes of $200-1000$~\% between the CRD and PRD solutions, and certain features only appearing due to redistribution. With \textsc{MS\_RADYN} we are able to include NEQ effects, but the line profiles are computed assuming CRD.  PRD is seemingly more important to include than NEQ effects. While the magnitudes of differences between NEQ and SE was generally smaller than the differences between CRD and PRD, and shorter lived (PRD was required throughout the duration of the flare), if one is interested in the initial flare heating, or in the decay phase, then NEQ effects should ideally be considered. 

A time dependent code capable of including both NEQ and PRD effects in flares should be developed, but this is a computationally demanding endeavour. For the moment, if we wish to study the initial heating and cooling, we can capture both NEQ and PRD by using a multi-step process. Radiative hydrodynamic flare atmospheres are produced by \texttt{RADYN}. The \texttt{RADYN} solutions are used in \texttt{MS$\_$RADYN} to obtain the NEQ CRD Mg~\textsc{ii} populations. Those NEQ populations are fed into \texttt{RH} which, using the modifications we described in \cite{kerr_2019} to fix the level populations, will solve the PRD Mg~\textsc{ii} radiation transfer.  
	
Another useful flare line observed by IRIS is Fe \textsc{xxi} 1354.1~\AA, which forms in 10MK plasma and is likely to experience non-equilibrium effects \citep{2011ApJS..194...26B}. It is commonplace to model this line using data from the \texttt{CHIANTI} atomic database \citep{1997A&AS..125..149D,2019ApJS..241...22D}, under the assumption of ionisation equilibrium and optically thin emission \citep[e.g.][]{2015ApJ...799..218Y,2019ApJ...879L..17P}. While the latter assumption is safe, ignoring non-equilibrium effects may not be correct. A similar comparison to that presented here for Mg~\textsc{ii} should be performed, comparing the predicted level Fe~\textsc{xxi} ion fraction from \texttt{MS\_RADYN} to those from \texttt{CHIANTI} given model atmospheres from \texttt{RADYN}. \\

\textsc{Acknowledgments:} \small{GSK was funded by an appointment to the NASA Postdoctoral Program at Goddard Space Flight Center, administered by USRA through a contract with NASA. This research was supported by the Research Council of Norway through its Centres of Excellence scheme, project number 262622, and through grants of computing time from the Programme for Supercomputing. JCA acknowledges funding support from the Heliophysics Supporting Research and Heliophysics Innovation Fund programs.}

\bibliographystyle{aasjournal}

