A parameter study for modeling MgII h and k emission during solar flares

TL;DR

This study investigates how various thermodynamic parameters affect the formation of MgII h and k emission lines during solar flares, using radiative transfer modeling to match observed spectral features.

Contribution

It demonstrates that adjusting densities, temperatures, and velocities in flare models can reproduce observed single-peaked MgII line profiles, advancing understanding of flare spectra formation.

Findings

Reproduced single-peaked MgII h and k spectral profiles.

Identified effects of increased densities, temperatures, and velocities on line formation.

Showed that unresolved flows can broaden line wings.

Abstract

Solar flares show highly unusual spectra, in which the thermodynamic conditions of the solar atmosphere are encoded. Current models are unable to fully reproduce the spectroscopic flare observations, especially the single-peaked spectral profiles of the MgII h and k lines. We aim at understanding the formation of the chromospheric and optically thick MgII h and k lines in flares through radiative transfer calculations. We take a flare atmosphere obtained from a simulation with the radiative hydrodynamic code RADYN as input for a radiative transfer modeling with the RH code. By iteratively changing this model atmosphere and varying thermodynamic parameters, such as temperature, electron density, and velocities, we study their effects on the emergent intensity spectra. We can reproduce the typical single-peaked MgII h and k flare spectral shape and their approximate intensity ratios to the subordinate MgII lines by either increasing densities, temperatures or velocities at the line core formation height range. Additionally, by combining unresolved up- and downflows up to ~250 km/s within one resolution element, we also reproduce the widely broadened line wings. While we cannot unambiguously determine which mechanism dominates in flares, future modeling efforts should investigate unresolved components, additional heat dissipation, larger velocities, and higher densities, and combine the analysis of multiple spectral lines.

Equations

$$I_{\nu }^0=\frac{1}{\mu }{\int }_z{S}_{\nu }{\chi }_{\nu }{e}^{-\frac{{\tau }_{\nu }}{\mu }}dz=\frac{1}{\mu }{\int }_z{C}_{i}dz,$$

Full text

1. Introduction

Although solar flares have been investigated for decades, there are still open questions, such as how the released energy is transported throughout the atmosphere and deposited in the lower atmospheric layers. By studying chromospheric spectral lines, we probe the underlying physical processes and the response of the lower atmosphere to flare heating. NASA’s most recent solar mission, the Interface Region Imaging Spectrograph (IRIS; De Pontieu et al., 2014) showed puzzling spectra of the near-UV (NUV) Mg II h and k resonance lines in flares, which could not be explained by any modeling efforts so far.

The Mg II h and k lines are an important contributor to the UV emission during flares (Lemaire et al., 1984). Their unexplained characteristics in flares are: 1) a lack of their central self-reversal, which is unusual for lines dominated by scattering; 2) very broad line wings (line center $\pm$ 1.5 Å); 3) the presence of the subordinate Mg II lines (3p-3d transition) in emission, which in the quiet Sun have been reported to be sensitive to heating in the low chromosphere (Pereira et al., 2015), but may form differently in flares; and 4) redshifts often occurring due to strong downward velocities (Kerr et al., 2015; Liu et al., 2015; Graham & Cauzzi, 2015).

The so far only direct comparison of observed flare Mg II spectra and hydrodynamic simulations (Rubio da Costa et al., 2016) showed that the peculiar shape of flare spectra cannot be reproduced yet and that simulations always obtain profiles with a central reversal, more similar to quiet Sun than to flares. There have been several attempts to understand the behavior of the Mg II spectra during solar flares (Machado et al., 1980; Lemaire & Gouttebroze, 1983; Kerr et al., 2016; Rubio da Costa et al., 2016; Liu et al., 2015; de la Cruz Rodríguez et al., 2016; Kowalski et al., 2017; Reid et al., 2017), but no simulation to date has reproduced any of the observed profiles and it is unclear, which physical mechanisms are responsible for the observed line profiles. Our goal for this paper is to carry out a parameter study to investigate the origin of the Mg II spectral shapes during flares and to understand the atmospheric conditions needed to reproduce the observed flaring UV spectra. Carlsson et al. (2015) performed a similar parameter study with the RH code, modifying the thermodynamic parameters in their model atmosphere, but enforcing hydrostatic equilibrium, to study the Mg II k profiles in plages with the objective of matching the observations. Here, we investigate even unlikely possibilities of variations, such as velocities or microturbulences significantly higher than those reported in observations, but we cannot provide an unambiguous solution to the conditions in a flaring atmosphere, which would require comparisons with multiple spectral lines and which will be our next step.

There are currently three codes commonly used for flaring hydrodynamic simulations: RADYN (Carlsson & Stein, 1997), FLARIX (Kašparová et al., 2009), and HYDRAD (Bradshaw & Cargill, 2013). The HYDRAD code does not account for optically-thick radiative losses, an important energy term in the chromosphere. Therefore, it is not a suitable tool to study the response of the lower atmosphere to flare heating and in particular the Mg II emission in the UV. The FLARIX code considers radiative transfer for hydrogen, calcium and more recently magnesium, but currently it does not consider the radiative losses of the helium atom, which may be up to 10% of the total energy radiative loss (Mariska et al., 1989; Liu et al., 2009), affecting the modeled chromospheric emission. In contrast, RADYN solves the hydrodynamic equations together with the detailed radiative transfer for the atoms dominating the radiative losses in the chromosphere; i.e. hydrogen, helium, calcium and magnesium. Here we therefore use a model atmosphere from a RADYN simulation to start with the most realistic flare atmosphere and vary its thermodynamic parameters to investigate their influence on the Mg II lines. The non-LTE radiative transfer code RH (Uitenbroek, 2001) is used to calculate the Mg II lines from the modified atmospheres. The advantage of RH is the proper treatment of the effects of angle-dependent partial frequency redistribution (PRD), improving the assumption of complete frequency redistribution (CRD) used by RADYN, which has previously been demonstrated to play an important role in the formation of the Mg II h and k line profiles (see e.g. Leenaarts et al., 2013a, b). Additionally, RH accounts for possible frequency overlap from bound-bound transitions (Rybicki & Hummer, 1991, 1992), whereas RADYN does not.

The paper is organized as follows: a description of the motivation of our Mg line spectra study during solar flares is shown in Section 2; a summary of the RADYN and RH codes is presented in Sections 3.1 and 3.2, followed by the study of the behavior of the Mg II h and k line profiles to different atmospheric conditions in Section 4; a more detailed analysis on the line formation is described in Section 5, and a summary and discussion of the results in Section 6.

2. Motivation

The observed shapes of the Mg II lines are very diverse and variations exist on pixel-to-pixel scales. Some examples from the X1 flare on March 29, 2014 (Kleint et al., 2015) are shown in Figure 1. Flare profiles range from single peaks with broad wings (e.g. middle panel red and purple spectra), to small central reversals (bottom panel purple), to strongly asymmetric red components (e.g. top panel red and blue, bottom panel blue). The locations of the profiles are indicated with color-coded horizontal lines in the image on top. Grey spectra show the quiet Sun for a comparison. RHESSI HXR contours from 32-100 keV are also drawn. Because the non-thermal counts were comparably low at this time, we needed to integrate for 30 s (17:47:20-17:47:50 UT), using detectors 2 to 7. The RHESSI contours are therefore not strictly simultaneous to the IRIS data, but the sources did not move significantly at this time (Battaglia et al., 2015). The blue profile of slit position 7 corresponds to a time just when accelerated electrons were detected at the same location. This spectral profile shows a reversal, though shifted to the blue, and additionally a downflow, and its intensity increases significantly. A few seconds later, such profiles turn into the broad single-peak profiles. For example, the HXR source passed the red and purple profiles of slit 4 about 30 to 90 seconds earlier.

The line ratios given in the images are calculated as integrals from 2795–2798 Å (k line), 2802–2805 Å (h line), and 2798–2800 Å (subordinate line at 2799 Å).

The magenta crosses are for reference only and indicate locations where the H$\alpha$ and Ca II 8542 Å were compared in Rubio da Costa et al. (2016). These spectral lines matched relatively well to the model atmosphere that is used for the parameter study in this paper.

Single-peaked profiles seem universal for all analyzed flares and are also commonly observed in sunspots, though in sunspot they appear without any broad wings and with only tiny, if any, emission from the subordinate lines. So our goal is to simulate single-peak broad Mg II h and k profiles (e.g. red or purple profile of slit position 4) and simultaneously produce emission in the subordinate 3p-3d Mg II lines at the correct intensity ratios, to understand the characteristics of the atmospheric parameters necessary to reproduce the observed typical line profiles during solar flares.

2.1. Formation of the Mg II h and k Spectral Lines

The Mg II h and k lines are resonance lines formed under conditions of non-LTE. Considering that the radiative damping is the largest contribution to the total damping in the chromosphere (at $z>500$ km), that Van der Waals broadening is the dominant contribution in the photosphere and that quadratic Stark broadening is only of minor importance, the radiative damping component seems to be the main contribution to the formation of the Mg II h and k line core, while the Van der Waals broadening seems to be the main contribution to the broad line wings (Leenaarts et al., 2013a). 3D effects have been reported to be very important, especially in the line core. Unfortunately, there are no simulations yet that are able to calculate Mg II realistically in 3D.

The 3p-3d level subordinate lines usually appear in absorption but show a significant enhancement during solar flares (e.g. Kerr et al., 2015). Their emission was analyzed by Pereira et al. (2015) who concluded them to be optically thick and forming in the lower chromosphere in the quiet Sun. They found them to turn into emission only rarely, especially during their quiet Sun simulations, and concluded that a large temperature gradient ($\geq~{}1500$ K) must be present in the lower chromosphere to cause the subordinate lines to turn into emission.

3. Modeling the line profiles

We start with a modeled flare atmosphere from a RADYN simulation and modify its parameters to investigate their influence on the Mg II lines, using the non-LTE radiative transfer code RH to properly treat the PRD effects. In the following sections 3.1 and 3.2 we will briefly introduce each code individually.

3.1. RADYN Code

We used the RADYN code of Carlsson & Stein (1997), including the modifications of Abbett & Hawley (1999) and Allred et al. (2005), to simulate the radiative-hydrodynamic response of the lower atmosphere to energy deposition by non-thermal electrons in a flare loop. We use the atmosphere from the run described in Rubio da Costa et al. (2016) at 17:45:31 UT and manually modify its temperature, density, non-thermal broadening or plasma velocity with the RH code (see Section 3.2). Additionally, we used RADYN simulations with electron beam fluxes from 109 to 1012 erg s*-1* cm*-2* to investigate the dependence of the spectral line shape on beam flux. We also compared the spectra resulting from the RADYN simulations described in Rubio da Costa et al. (2015b), including particle transport and stochastic acceleration with the ones assuming an ad-hoc single power-law, derived from fitting the non-thermal component of the hard X-ray spectra to a single power-law and applying the collisional thick-target modeling (Brown, 1971; Emslie, 1978).

Considering that the Mg II h and k lines are strongly affected by the effect of PRD, their simulations with RADYN, which only assume CRD, are not sufficiently realistic. We therefore recalculate the Mg II lines with the RH code (see Section 3.2) including PRD effects. As an input to RH we use snapshots of the atmosphere calculated with RADYN with manually modified parameters.

3.2. RH Code

We performed non-LTE radiative transfer computations with a modified version of the RH code (Uitenbroek, 2001; Pereira & Uitenbroek, 2015), which includes the heating due to the injection of electrons in the atmosphere. The modification (performed by J. Allred; private communication) includes the computation of the collisional ionization rate from non-thermal electrons, following Ricchiazzi & Canfield (1983), given the electron heating rate estimated from the RADYN code. RH can treat the effects of angle-dependent partial frequency redistribution using the fast approximation by Leenaarts et al. (2012). We use a 10 level plus continuum model atom for Mg II as described in Leenaarts et al. (2013a).

One of the differences between RH and RADYN is that RADYN uses the Uppsala opacity package (Gustafsson, 1973), while RH calculates the opacity of the transitions in the indicated atoms. Additionally, by using snapshots of the atmospheric model from RADYN as input into the RH code, we recalculate the ionization populations using statistical equilibrium, instead of using the non-equilibrium ionization already estimated by RADYN. This changes the plasma density and, as a consequence, shifts the height scale in comparison with the atmosphere produced by RADYN, leading to differences in the formation height in units of length. Although statistical equilibrium can be a poor assumption at the beginning of the impulsive phase of the simulation, the effect should be smaller at later times. Unfortunately there is currently no option to avoid statistical equilibrium in RH. We would like to emphasize that the goal is to study the behavior of the Mg II line for different atmospheric parameter changes. Even though we use the ionization population based on the statistical equilibrium assumption for excitation processes, we still keep the non-equilibrium electron density estimated from RADYN, rather than recalculating the equilibrium electron density, since it is the one resulting from the hydrodynamic equations. Therefore, we believe it to be more realistic and the electron density population is far from equilibrium in a flaring atmosphere. By not solving for hydrostatic equilibrium, we have not enforced the charge conservation in RH. We find that the velocity variations in the atmosphere do not displace the charge balance of positive vs. negative charges significantly ($<$ 4%). The temperature and density variations affect the charge conversation more with the exact numbers listed in the respective sections below.

These points may make a difference when comparing the line profiles with the observations; for instance, the continuum emission at 2992 Å resulting from RH is less than 20% lower than the one estimated from RADYN. Another possibility that can explain the lower continuum emission estimated with RH, in comparison with RADYN, can be the lower opacity in the wings of the Mg II h and k lines. Kowalski et al. (2017) have included the Mg II h and k wing opacity in the calculations of the excess continuum intensity finding that the opacity at these wavelengths is important for an accurate treatment of the Mg II wing emission as well as the continuum in the upper photosphere. An underestimation of the opacities can result in a lower continuum emission, by a factor of 15–-30%.

Each line profile is calculated with up to 100 frequency points across the transition using a Voigt profile and a line broadening due to the Stark effect. Throughout the paper we use $\mu$=0.77 to be consistent with the observations in Figure 1.

4. Synthetic line profiles

In this section we modify one thermodynamic parameter at a time. While this probably is not representative of real solar changes during flares, it allows us to investigate the dependence of the Mg II spectra on different parameters. We ran variations of several hundred atmospheres and show the main results in the following sections.

4.1. Varying temperature

We varied the temperature structure, including (1) temperature steps at different heights below the transition region, (2) higher chromospheric temperatures, and (3) localized heating at certain heights. The electron density and plasma velocities were kept constant for all runs. Figure 2 shows some representative runs in which the temperature increases just below the transition region, where the core of the Mg II k and h lines are formed. A strong temperature increase in this region leads to single-peaked profiles whose intensity is lower than the original atmosphere. However, the intensities of the subordinate lines are too low compared to the h and k lines. The dotted gray line corresponds to the red spectrum at slit 4 in Figure 1 with its intensity multiplied by 5. The observed line wings are clearly much broader.

We also increased the temperature of the chromosphere (red and blue lines in Figure 3) or introduced strong heating at limited heights (yellow and light blue lines in Figure 3). Even when unrealistically increasing the chromospheric temperature by 3000 K (blue), the line wings are only slightly broadened, far less than the observed broadening. Temperature spikes in the middle chromosphere (yellow) lead to sharp spikes in the line wings, which to our knowledge have not been found in observations. Temperature spikes in the upper chromosphere (light blue) lead to a similar behavior as shown in Figure 2.

As long as there is a steep temperature increase at the formation height of the line core, the profile turns into single peak and the line is insensitive to the temperature behavior above these heights.

Because the charge conservation has not been enforced, the charge balance is displaced. Considering the atoms of hydrogen, helium, and magnesium, the maximum ratio between positive and negative charges is 1.16 at 10000 K and 1.01 at 50000 K in the models of Figure 2.

4.2. Varying electron density

In this set of models, we kept the temperature constant and increased the electron density at different heights. Figure 4 shows that heights just below the transition region have no effect on the shape of the line core (yellow and red lines). The comparison of yellow and light blue lines shows that the density between log(column mass) between $-3.6$ and $-3.8$ g cm*-2*, which corresponds to 1.15-1.17 Mm, is responsible for the single peak of the h and k lines. By increasing the density at least by a factor of 10, the reversed profile turns into a single peak (black vs. blue vs. red profiles). The k/h ratio is similar to the observations, the k/2799 ratio differs by a factor of 2, however the intensity is significantly too high in the models. Therefore the observed spectrum (dotted gray) was multiplied by a factor of 10.

Figure 5 shows the effect of more smoothly increasing and decreasing densities. Similarly to the previous figure, only densities at the formation height of the line core matter for obtaining single-peaked profiles.

Kowalski et al. (2017) have recently studied the continuum emission together with the asymmetries observed by IRIS in the Fe II 2832.39 Å line profile during the impulsive phase of an X1.0 class flare. They suggest a model producing two flaring regions (a condensation and stationary layers) in the lower atmosphere. The condensation, due to the high densities in the chromosphere, compresses the chromosphere and therefore enhances the continuum emission; while the stationary flare layers in the chromosphere explain the bright redshifted Fe II emission component, reported to be observed during solar flares. This plausible scenario with two flaring regions in the chromosphere, could explain the sudden density increase within a narrow region.

Another possible physical mechanism that can explain the increase of the electron density without elevating the temperature is the non-thermal ionization resulting from a beam of non-thermal electrons, creating a compressed region in the chromosphere for short timescales. The lack of charge conservation leads to unrealistic values (increase of negative charges by an order of magnitude at T=10000 K) in the yellow and light blue examples in Figure 4, because we modified the density starting low in the atmosphere. For all other atmospheres in Figure 4 and 5, the ratios are maximally 1.19 and 1.01 at T=10000 K and 50000 K, respectively.

4.3. Varying Doppler velocity

The Doppler velocities may change strongly during flares, generally from strong upflows at coronal temperatures to downflows of few tens of km s*-1* in the chromosphere (e.g. Milligan & Dennis, 2009). For these tests, we vary the magnitude of the downflow at and near the line core formation height. Figure 6 shows the approximate Doppler velocities from Milligan & Dennis (2009) in light blue (negative velocities represent downflows), which lead to an asymmetric profile with a stronger blue peak. It is possible to obtain a single-peak profile with its center at the line core just by varying velocities (e.g. yellow). But such variations do not lead to a symmetric profile.

All single-peak profiles plotted in Figure 6 (yellow, blue, red) show an extra enhancement at the blue side of the line, at around 2795.9 Å, similarly to the blue wing enhancement observed by IRIS in the bottom panel of Figure 1 (blue profile). Such IRIS profiles seem to only occur during a few seconds and exactly at locations where HXR footpoints are observed. Figure 7 shows the “HXR-footpoint” IRIS profile in dotted gray. We noticed that by starting the velocity variation lower in the chromosphere the emission becomes red-shifted (dark blue vs. red in Figure 7), but the location of the dip at $\sim$2796 Å remains. We combined the red and blue profiles (35% blue + 65% red) into a new profile (light blue), simulating what would be observed if these two components existed within one pixel. The light blue profile qualitatively matches the “HXR-footpoint” IRIS profile very well, both in the shape and extent of the red wing and in the location of the dip. We re-ran RH with the velocities derived from this combination (yellow), but the resulting profile shows that a single component cannot explain the observations. While the non-perfect fit of the emission feature at $\sim$2796 Å and of the subordinate lines indicates that the middle-lower chromosphere is not fully reproduced yet, our results indicate that “HXR-footpoint” IRIS profiles are due to unresolved strong downflows.

To study how the emission of the Mg II k line responds to the combination of plasma motions in opposite directions, we combined a strong downflow (red line in Figure 8) with a moderate upflow (blue) into a new profile (yellow, consisting of 35% blue + 65% red). The existence of unresolved up- and downflows in the same pixel could explain the broad wings. We increased the turbulent velocity to 15 km s*-1*, otherwise the absorption feature (as seen in Figure 7) affected the shape of the line wings. Also, the downflow velocities have to be larger than the upflows, because the red wing is slightly more extended in observations. Since nearly all IRIS flare ribbon pixels have these typical profiles, this result may indicate that there are similar unresolved up- and downflows throughout the ribbons.

4.4. Varying the radiation field

By artificially including a UV radiation field in the upper corona in the RH code, including emission lines such as He I, He II and Mg X, we aim at compensating for energy losses of optically thin lines not considered in the RH code, which may influence the resulting synthetic emission during solar flares, especially in the UV.

We considered a wavelength range between 302 and 613 Å and estimated the total irradiance (Pietarila & Judge, 2004; Judge & Pietarila, 2004). Comparing the resulting Mg II h and k line profiles with the original ones without external irradiation field in the upper corona, we found that the differences are minimal, especially in the wings. The emission in the line core is 0.02% higher when including the external irradiation field in the upper corona. Therefore, the inclusion of coronal irradiance does not contribute enough to the total Mg II UV emission, to explain the differences between the synthetic line profiles and those observed by IRIS.

4.5. Varying the transition region

For this set of models we vary the location of the transition region by modifying the column mass at which the transition region is located. Figure 9 shows the modifications of temperature and density. The lower the transition region, the lower the intensity of Mg II h and k. These variations do not reproduce typical flare profiles. Even the yellow single-peaked profiles are not similar to observed profiles because the subordinate lines have higher intensities than the h and k lines.

4.6. Varying the non-thermal broadening

Since the resulting synthetic line profiles are narrower than the observed flaring Mg II profiles, we also tried to modify the non-thermal broadening parameter, by increasing it in the upper chromosphere, at a column mass higher than -3.1 g cm*-2*. Figure 10 shows the broadening of the line profile for different test runs, while increasing the non-thermal broadening and varying the region where this occurs. For these simulations, we also increased the electron density in a narrow layer at the upper chromosphere, where the line core is formed, to obtain a single core line emission to be able to better compare the effects of microturbulence to observed flare profiles.

Considering that typical values for microturbulence to compensate for the lack of small-scale random motions in the model are $\sim 3$ km s*-1* (de la Cruz Rodríguez et al., 2012, 2016), even unrealistically high velocities of 40 km s*-1* are not enough to explain the wide broad wings observed by IRIS. For instance de la Cruz Rodríguez et al. (2016) use 8 km s*-1*, describing this value as very high. Moreover, the line profile not only becomes broader, but also flatter in the core of the line, contrary to the observations.

4.7. Varying the flux of non-thermal electrons

We studied how the variation of the non-thermal electron flux influences the Mg II UV emission by using atmospheres from different RADYN simulations and calculating the profiles with RH. The resulting atmospheres correspond to four different simulations with a constant total non-thermal electron energy flux of $10^{9}$, $10^{10}$, $10^{11}$ and $10^{12}$ erg cm*-2* s*-1* and assuming that the electron spectra follow a single power law at non-thermal energies, with a cutoff energy of $E_{c}=25$ keV and a spectral index $\delta=4$ for the first three simulations and $\delta=3$ for the simulation with the highest electron energy flux.

Since the atmosphere evolves differently for each electron energy flux simulation, it is hard to make a direct comparison of the line profiles and assess a direct contribution of the electron flux. We tried to minimize this effect by selecting a timestep at the beginning of the simulation, when the atmosphere has not reached the “explosive” chromospheric evaporation described by Allred et al. (2005). We selected a timestep of $t=5$ seconds after heating the atmosphere for the first two runs (F09 and F10), at which point the atmosphere has not reached the “explosive” chromospheric evaporation; $t=6.1$ seconds in the case of F11, having very strong upflow velocities in the lower corona of almost 1200 km s*-1* as well as downflowing velocities close to 25 km s*-1*; and $t=2.8$ seconds for F12, with downflowing velocities of 70 km s*-1*. Figure 11 shows how much the atmospheres differ within the runs.

In general, increasing the beam strength leads to higher intensities in the Mg II UV wavelength range, especially in the continuum regime as well as the subordinate line profiles. The asymmetries (higher red peaks) of the F09 and F10 runs (red and blue lines of Figure 11) are mostly due to the downflowing velocities (negative values) at a column mass region of $\approx-4$ g cm*-2*, where the line core is formed. The line profile resulting from the F11 run (green line) is the strongest among the four cases, having also the highest electron density and temperature, more than one order of magnitude increase in the upper chromosphere, between -3.4 and -4 g cm*-2*. The strong slightly blueshifted line core is due to the upflowing (positive) velocities of almost 25 km s*-1* in the upper chromosphere. The F12 simulation (violet line) shows lower intensities in the line core and a lack of reversal at 2796.35 Å due to the sudden temperature increase below a column mass of -4 g cm*-2*, decoupled from the electron density, since the sudden temperature increase does not reflect such sudden changes in the electron density. The secondary emission peak at higher wavelengths (2796.8 Å) is due to the sudden changes in the velocity stratification, from upflowing to downflowing plasma motions, located at the column mass of -4 g cm*-2*. These changes are associated to the sudden temperature increase, located at that height.

Even though we only show the line profile from a specific time step of the RADYN simulations, we have performed a study of the temporal evolution of the line profiles. We found that none of the simulations showed Mg II h and k as single-peaked profiles for any beam strength at any time step.

4.8. Including particle transport and acceleration

By considering the coupling of particle transport and stochastic acceleration using the FLARE code (Petrosian et al., 2002) with the radiative transfer and hydrodynamics from the RADYN code, we obtain a spectrum formed by a quasi-thermal component and a power-law tail as explained by Rubio da Costa et al. (2015b). The atmospheric evolution of the run when coupling of the particle transport and acceleration equations with the hydrodynamics results in a stronger chromospheric evaporation as well as greater up- and downflowing velocities, when compared with the atmosphere resulting from the single power-law run.

By taking the atmosphere from the hydrodynamic simulations of Rubio da Costa et al. (2015b) after 60 seconds of heating, we studied how the Mg II h and k line profiles are affected when coupling particle transport and stochastic acceleration. We find that although the intensity decreases by almost a factor of 2, the line profiles still show narrow wings and a core reversal at the line center (see Figure 12).

The line profiles resulting from the power-law run (red line) have stronger intensities in the line core, but almost the same contribution in the continuum as the stochastic acceleration run. This is because, although the temperature and electron densities are about the same order of magnitude, the atmosphere resulting from the power-law run shows slightly higher densities and temperatures at a column mass between -3.4 and -4 g cm*-2*.

The inclusion of particle transport and acceleration results in higher upflowing plasmas close to the transition region (below to a column mass of -6 g cm*-2*) and stronger downflowing velocities in the upper chromosphere (at a column mass of -5 g cm*-2*). The line core is formed lower in the atmosphere, at a column mass close to -3.6 g cm*-2*; therefore the plasma motions that significantly differ in both runs do not affect the line formation. Instead, the line profiles are symmetric, as the synthetic spectra in Figure 12 show.

5. Formation of the Mg II k line profile: intensity contribution function

To investigate the behavior of the source function and where in the atmosphere the line is formed, we calculate the so-called contribution function. By writing the formal solution of the transfer equation for emergent intensity (Equation 1; Carlsson & Stein, 1997), we can investigate how the atmospheric stratification affects the formation of the line profile and where exactly the line is formed. The line emergent intensity $I_{\nu}^{0}$ can be written as:

[TABLE]

where $z$ is the atmospheric height; $\tau_{\nu}$ is the monochromatic optical depth; $\chi_{\nu}$ is the monochromatic opacity per unit volume; $S_{\nu}$ is the source function, defined as the ratio between the emissivity to the opacity of the atmosphere and $C_{i}$ is the intensity contribution function, indicating how much emergent intensity originates at a certain height $z$.

The intensity contribution function, $C_{i}$, can be expressed as the product of $S_{\nu}$, $\frac{\chi_{\nu}}{\tau_{\nu}}$ and $\tau_{\nu}e^{-\frac{\tau_{\nu}}{\mu}}$, as the panels in Figure 13 show (for a more detailed explanation see Leenaarts et al., 2013b). Since both Mg II h and k line profiles have a similar behavior, we will focus on the Mg II k line profile in the following sections.

5.1. Original RADYN atmosphere

The Mg II k line profile resulting from the simulation using the original RADYN atmosphere shows a reversal at the line center (top left panel of Figure 13). The asymmetry in the blue wing of the line is associated to the asymmetry in the $\tau_{\nu}=1$ layer at different frequencies, as the term $\tau_{\nu}e^{-\frac{\tau_{\nu}}{\mu}}$ in the upper right panel indicates, since its contribution to the total intensity at each frequency has an asymmetry towards blue wavelengths. The source function $S_{\nu}$ in the lower left panel, calculated for the line center wavelength, closely follows the Planck function (dotted line) in the lower atmosphere, up to 1.14 Mm in comparison with the sudden increase at higher layers, due to the temperature stratification. The departure of the source function from the Planck function means that the line core is formed under non-LTE conditions, because both profiles have a different temperature response.

The intensity contribution function $C_{i}$ is the product of the previous three panels: the main contribution comes from the upper chromosphere, between 1.06 and 1.14 Mm, since the ratio of the opacity to the optical depth ($\frac{\chi_{\nu}}{\tau_{\nu}}$) is very small at lower heights. The integration of the intensity contribution function at each frequency results in the black line profile displayed in the lower right panel. The core reversal is due to the sudden decrease of the source function due to the decoupling of the source and Planck functions at 1.14 Mm, where the photons of the line core are formed.

By studying the formation of the subordinate lines between the 3p${}^{2}P$ and 3d${}^{2}D$ states (figure not shown), we found that they are formed in the upper chromosphere, between 0.83 and 1.14 Mm, at higher layers than the values reported by Pereira et al. (2015) for the quiet Sun.

5.2. Atmosphere with increased temperature

By increasing the temperature in the upper chromosphere, as described by the yellow atmosphere of Figure 2, we find that the formation height of the line center is reduced by $\sim$200 km (see top right panel of Figure 13), explaining the steeper peak at the line core. The $\tau_{\nu}=1$ layer in the line wings has the same shape as in the top left panel of Figure 13 with the exception of the plateau at wavelengths close to the line core. The redshifted line core is due to the downflowing velocities of $\sim$ 14 km s*-1* present at the line core formation height (between 1.06 and 1.14 Mm). Although the ratio between the opacity and the optical depth is almost the same as in the previous atmospheric example, there is no clear emissivity above $\tau_{\nu}=1$ close to the line center. The temperature increase results in a coupling of the source function to the Planck function at higher heights, especially also where the line core is formed (lower left panel). This indicates LTE conditions for the whole line and explains the single-peaked profile.

The asymmetry of the line profile is reflected by the frequency distribution of the intensity contribution function. There is a secondary intensity contribution function component at a height of 1.13 Mm, and above $\tau_{\nu}=1$, which indicates that there is some contribution to the core emission, which is formed under optically thin conditions.

5.3. Atmosphere with increased density

The bottom left panel of Figure 13 shows the different terms of the total intensity contribution function for the yellow atmosphere of Figure 4. Here the electron density increase results in a slightly stronger opacity component $\chi_{\nu}$ at a height of 1.14 Mm, in comparison with the resulting opacity of Section 5.2, with an increased atmospheric temperature. The term $\chi_{\nu}/\tau_{\nu}$ is also stronger in the line center, in comparison with the original atmosphere of Section 5.1 because $\chi_{\nu}$ is proportional to the density of emitting particles. Therefore, $\frac{\chi_{\nu}}{\tau_{\nu}}$ is higher when there is a large number of emitters at low optical depth (i.e., the produced photons can escape), explaining the stronger emissivity component in comparison with the top right panel of Figure 13.

The main difference in comparison with the increase of temperature is that the coupling between the source function and the Planck function is located at slightly higher heights, as well as the formation height ($\tau_{\nu}=1$) of the line core. By increasing the electron density, the source function $S_{\nu}$ keeps increasing with the temperature, although it is decoupled from the Planck function $B_{\nu}$. As in Section 5.2, the source function is still coupled to the Planck function at the core formation height, indicating LTE conditions. The line core shows a flatter profile because the intensity contribution function has a more symmetric frequency distribution, with a fainter blue wing asymmetry.

5.4. Atmosphere with increased velocities

Previous simulations have shown that asymmetries in chromospheric line profiles during flares are strongly affected by the velocity field in the flaring atmosphere (Kuridze et al., 2015). The bottom right panel of Figure 13 shows that by increasing the downflowing plasma velocity to about 200 km s*-1*, as shown in the blue atmosphere of Figure 7, the line profile is formed in a broader region, in comparison with the original atmosphere of the top left panel of Figure 13. The intensity contribution function also covers a much broader frequency range, broadening the wings of the line profile (note the different scaling of the x-axis). The Planck and source functions have similar behavior, increasing with temperature in the lower atmosphere, departing from the Planck function above 1.12 Mm, similarly to the original atmosphere. At the core formation height, just above 1.12 Mm, the source function is starting to decouple from the Planck function, although both still increase with temperature, resulting in a single core in emission. Due to the plasma motions, the line shows an asymmetric profile, especially towards redder wavelengths, as a result of the increased downflowing plasma motions.

By increasing the plasma velocity, the formation height at which $\tau_{\nu}=1$ (green line) is not flat near the line core, as in the previous cases, but it linearly increases within the wavelength range $-0.2$ to $1.6$ Å.

For other variations of the velocity, we can either have still coupled Planck and source functions (e.g. red line in Figure 6), or an already decoupled and decreasing source function with temperature (e.g. red line in Figure 7). Yet we do not observe a reversal in any of the cases because even with a decreasing source function, the large velocities fill the intensity where usually a reversal is observed.

It is also worth noting the secondary contribution coming from heights close to 1.2 Mm (not shown), above the $\tau_{\nu}=1$ curve, indicating that there is a contribution to the intensity in the red wing (at wavelengths between 0.6 and 2.1 Å) formed under optically-thin conditions.

6. Summary and Discussion

6.1. Options to simulate flare-like Mg profiles

Our parameter study has shown several methods to obtain more flare-like Mg II h and k profiles in simulations with a single peak in emission. In Section 4.1 we have seen that by increasing the temperature at the line core formation height, we do not only obtain a single peak in emission, but also decrease the peak intensity. This is favorable, considering that synthetic line profiles tend to show higher intensities than IRIS observations (Rubio da Costa et al., 2016; Kowalski et al., 2017). This modification also results in symmetric line profiles, similar to the typical observed flare ribbon profiles. While the line ratio between the Mg II k and Mg II h lines looks reasonable for most of our simulations, the ratio between them and the subordinate triplet lines is too large for the atmosphere with a modified temperature (cf. Figure 2). This may indicate that the temperature or its gradient in the lower chromosphere may also need to be modified. However, this can be tricky as it might affect other lines, such as H\textalpha or Ca II 8542 Å, which already show a good agreement with observations (see Figure 12 in Rubio da Costa et al. (2015a), or Figure 11 in Rubio da Costa et al. (2016), for an easier comparison, we marked the locations in Figure 1).

The increase of the electron density results in a line core in emission as well, but as can be seen in Figure 4, the intensity of the line increases significantly, to values not reached by observations. The line ratio of Mg II k to subordinate lines is closer to the observations for the electron density variation than for the temperature variation.

In the case of increased downflowing velocities, the line core may show a single-peak in emission as well. The formation height at which the optical depth $\tau_{\nu}=1$ increases linearly near the center of the line, with its peak at 1.15 Mm and 1.6 Å. At this height, $S_{\nu}$ and $B_{\nu}$ are already decoupled. The source function at the line core formation heights for the velocity variations can either be coupled to or decoupled from the Planck function, depending on the velocity stratification.

The intensity contribution function $C_{i}$ shows a secondary component in emission, formed at heights above $\tau_{\nu}=1$ and therefore under optically-thin conditions. However, the resulting line profile is highly asymmetric, contrary to standard flare profiles. Only by combining profiles with different velocity structures, one can reproduce different IRIS profiles, thus indicating that the IRIS profiles may be unresolved. By combining downflows velocities, starting at different heights, the shape of the short-lived IRIS profile from HXR emission sites can be reproduced. The single peak and broad wings of typical IRIS flare profiles can be reproduced by combining up- and downflows. The unresolved components could also be explained by a combination of several threads overlaying in space (Warren, 2006; Rubio da Costa et al., 2016; Reep et al., 2016). The subordinate lines are comparably low in this case, and would require further modifications in the lower chromosphere as discussed above.

By increasing the temperature or density (Section 4.1 and 4.2 and purple and blue lines in Figure 14), the line core is formed under LTE conditions; while the velocity variations may result in a decoupling of the source function from the Planck function, therefore non-LTE conditions (see green line in Figure 14). In these cases we still have single peak profiles due to the large asymmetries of the line. It is possible that these conditions occur during flares, possibly through a combination of all parameters discussed above.

6.2. Options that fail to simulate flare-like Mg profiles

The increase of the chromospheric temperature (Figure 3) as well as the shift of the transition region towards lower heights (Figure 9) does not result in realistic line profiles. Either the line core is still reversed, or the subordinate lines have stronger intensities than those of the Mg II h and k lines, showing a discrepancy with respect to IRIS flare observations.

As discussed in Section 4.6, we also discarded a higher non-thermal broadening as a possibility to explain the broad observed line wings. Although the Van der Waals broadening seems to be the main contribution to the broad line wings, this term does not explain the widely broadened observed line wings. Although the RH code considers quadratic Stark broadening in the calculation, the adiabatic approximation used in RH for the quadratic Stark effect for the Mg atom underestimates the broadening (Mihalas, 1978). Therefore, we should keep in mind that the implementation of the quadratic Stark effect may not be sufficiently accurate. Very wide Mg II line profiles have also been observed in stellar flares (i.e. on YZ CMi, observed by Hawley et al., 2007). Hawley et al. (2007) discussed non-thermal beam excitations as the reason for the broadening of the Mg II wings. They argue that unresolved flows with high velocities ($\sim$ 200 km s*-1*), which would result in high kinetic energies that may exceed the radiated energy in the flare, are a possible explanation of the broad wings. Further studies with higher non-thermal electron fluxes will clarify this point.

A larger opacity in the photosphere would result in broader wings and this could be accomplished by increasing the temperature or the density, but would probably affect other spectral lines whose fit currently is reasonably good. 3D effects, as discussed by Leenaarts et al. (2013a), may have an important contribution on the formation of the Mg II h and k line emission, especially in the line core, but are currently not computationally feasible.

The importance of the advection term in the statistical equilibrium equation for Mg II it is still unclear. Cool flows flowing into hot plasma or hot flows into cool plasma could result in non-equilibrium of the ionization and excitation.

We would like to emphasize that this is purely a parameter study and we do not claim that the density or temperature must increase by an order of magnitude in flares. We do not have a physical explanation yet as to what is occurring in a flaring atmosphere. Nevertheless, we can conclude that a critical piece is missing in current hydrodynamic simulations, possibly more particles that may increase densities, more heat dissipation, or unresolved and stronger velocities. For future work, we plan to simultaneously compare many spectral lines, which will better constrain the atmospheric parameters at different heights and will potentially allow us to better constrain the conditions in flaring atmospheres.

We would like to thank P. Judge, P. Heinzel, W. Liu, R. Rutten, G. Kerr, M. Carlsson, T. Pereira, and the anonymous referee for their helpful discussions. Work performed by F.R.dC. is supported by NASA grants NNX13AF79G, NNX14AG03G, 8100003073 and NNX17AC99G. L.K. was partially supported by NASA grant NNX13AI63G. F.R.dC. and L.K. thank ISSI and ISSI-BJ for the support of the team “Diagnosing heating mechanisms in solar flares through spectroscopic observations of flare ribbons”. We gratefully acknowledge the use of supercomputer resources provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center. IRIS is a NASA small explorer mission developed and operated by LMSAL with mission operations executed at NASA Ames Research center and major contributions to downlink communications funded by ESA and the Norwegian Space Centre. CHIANTI is a collaborative project involving George Mason University, the University of Michigan (USA), and the University of Cambridge (UK).