Three-dimensional modeling of chromospheric spectral lines in a simulated active region

TL;DR

This paper uses 3D MHD simulations and radiative transfer modeling to study chromospheric features like fibrils and flare ribbons, revealing their formation and structure in active regions.

Contribution

It presents a self-consistent 3D simulation reproducing chromospheric features and analyzes their formation, providing new insights into the role of thermal conduction in flare energy transport.

Findings

Fibrils span up to 4 Mm and are visible in multiple lines.

Flare ribbons are formed by dense chromosphere and transition region structures.

Simulated ribbons show profiles similar to observations, suggesting thermal conduction's role.

Abstract

Because of the complex physics that governs the formation of chromospheric lines, interpretation of solar chromospheric observations is difficult. The origin and characteristics of many chromospheric features are, because of this, unresolved. We focus here on studying two prominent features: long fibrils and flare ribbons. To model them, we use a 3D MHD simulation of an active region which self-consistently reproduces both of them. We model the H$\alpha$, Mg II k, Ca II K, and Ca II 8542 {\AA} lines using the 3D non-LTE radiative transfer code Multi3D. This simulation reproduces long fibrils that span between the opposite-polarity sunspots and go up to 4 Mm in height. They can be traced in all lines due to density corrugation. Opposite to previous studies, H$\alpha$, Mg II h&k, and Ca II H&K, are formed at similar height in this model. Magnetic field lines are aligned with the H$\alpha$ fibrils, but the latter holds to a lesser extent for the Ca II 8542 {\AA} line. The simulation shows structures in the H$\alpha$ line core that look like flare ribbons. The emission in the ribbons is caused by a dense chromosphere and a transition region at high column mass. The ribbons are visible in all chromospheric lines, but least prominent in Ca II 8542 {\AA} line. In some pixels, broad asymmetric profiles with a single emission peak are produced, similar to the profiles observed in flare ribbons. They are caused by a deep onset of the chromospheric temperature rise and large velocity gradients. The simulation produces long fibrils similar to what is seen in observations. It also produces structures similar to flare ribbons despite the lack of non-thermal electrons in the simulation. The latter suggests that thermal conduction might be a significant agent in transporting flare energy to the chromosphere in addition to non-thermal electrons.

Equations

$$ϵ=\frac{C_{ji}}{C_{ji}+A_{ji}+B_{ji}{B}_{\nu }},$$

$$F_{\!i}=\dfrac{\mathrm{d}n_{i}}{\mathrm{d}t}\equiv n_{\text{p}}\bigl{[}R_{\text{p},i}(n_{\text{e}})+C_{\text{p},i}(n_{\text{e}})\bigr{]}+\sum_{\makebox[0.0pt]{$\scriptstyle{j=1,j\neq i}$}}^{3}n_{\!j}\bigl{[}R_{\!j,i}+C_{\!j,i}(n_{\text{e}})\bigr{]}\\[-4.0pt] \qquad-n_{i}\Bigl{[}R_{i,\text{p}}(n_{\text{e}})+C_{i,\text{p}}(n_{\text{e}})+\sum_{\makebox[0.0pt]{$\scriptstyle{j=1,j\neq i}$}}^{3}\bigl{[}R_{i,j}+C_{i,j}(n_{\text{e}})\bigr{]}\Bigr{]}=0\\[-4.0pt] \quad\mbox{for }i=\{1,2,3\},$$

$$F_{\text{\!p}}=\dfrac{\mathrm{d}n_{\text{p}}}{\mathrm{d}t}\equiv\sum_{j=1}^{3}n_{\!j}\bigl{[}R_{\!j,\text{p}}(n_{\text{e}})+C_{\!j,\text{p}}(n_{\text{e}})\bigr{]}\\[-4.0pt] -n_{\text{p}}\Bigl{[}\sum_{j=1}^{3}\bigl{[}R_{\text{p},j}(n_{\text{e}})+C_{\text{p},j}(n_{\text{e}})\bigr{]}\Bigr{]}=0,$$

$$n_{\text{H}}^{\text{tot}}=n_1+n_2+n_3+n_{\text{p}}.$$

$$n_{\text{e}}=n_{\text{p}}+(\text{mbox}chargeofallpositiveionsbuthydrogen),$$

$$F_{\text{\!e}}=\dfrac{\mathrm{d}n_{\text{e}}}{\mathrm{d}t}\equiv n_{\text{e}}-\!\!\!\underbrace{n_{\text{p}}}_{\text{non-LTE}}\!\!\!-\,n_{\text{H}}^{\text{tot}}\underbrace{\biggl{[}\sum_{s=2}^{\text{species}}\!\!\alpha_{\!s}\!\sum_{i=2}^{\text{ions}}f_{\!s,i}(n_{\text{e}},T)\biggr{]}}_{\text{LTE}}=0,$$

$$\pmatrix{}{\partial F_{\!1}}{\partial n_{1}}&\frac{\partial F_{\!1}}{\partial n_{2}}\frac{\partial F_{\!1}}{\partial n_{3}}\frac{\partial F_{\!1}}{\partial n_{\text{p}}}\frac{\partial F_{\!1}}{\partial n_{\text{e}}}\\ \frac{\partial F_{\!2}}{\partial n_{1}}\frac{\partial F_{\!2}}{\partial n_{2}}\frac{\partial F_{\!2}}{\partial n_{3}}\frac{\partial F_{\!2}}{\partial n_{\text{p}}}\frac{\partial F_{\!2}}{\partial n_{\text{e}}}\\ \frac{\partial F_{\!3}}{\partial n_{1}}\frac{\partial F_{\!3}}{\partial n_{2}}\frac{\partial F_{\!3}}{\partial n_{3}}\frac{\partial F_{\!3}}{\partial n_{\text{p}}}\frac{\partial F_{\!3}}{\partial n_{\text{e}}}\\ \frac{\partial F_{\text{\!p}}}{\partial n_{1}}\frac{\partial F_{\text{\!p}}}{\partial n_{2}}\frac{\partial F_{\text{\!p}}}{\partial n_{3}}\frac{\partial F_{\text{\!p}}}{\partial n_{\text{p}}}\frac{\partial F_{\text{\!p}}}{\partial n_{\text{e}}}\\ \frac{\partial F_{\text{\!e}}}{\partial n_{1}}\frac{\partial F_{\text{\!e}}}{\partial n_{2}}\frac{\partial F_{\text{\!e}}}{\partial n_{3}}\frac{\partial F_{\text{\!e}}}{\partial n_{\text{p}}}\frac{\partial F_{\text{\!e}}}{\partial n_{\text{e}}}\pmatrix{\delta}n_{1}\vphantom{\frac{\partial F_{1}}{\partial_{p}}}\\ \delta n_{2}\vphantom{\frac{\partial F_{1}}{\partial_{p}}}\\ \delta n_{3}\vphantom{\frac{\partial F_{1}}{\partial_{p}}}\\ \delta n_{\text{p}}\vphantom{\frac{\partial F_{1}}{\partial_{p}}}\\ \delta n_{\text{e}}\vphantom{\frac{\partial F_{1}}{\partial_{p}}}=-\pmatrix{F}_{\!1}\vphantom{\frac{\partial F_{1}}{\partial_{p}}}\\ F_{\!2}\vphantom{\frac{\partial F_{1}}{\partial_{p}}}\\ F_{\!3}\vphantom{\frac{\partial F_{1}}{\partial_{p}}}\\ F_{\text{\!p}}\vphantom{\frac{\partial F_{1}}{\partial_{p}}}\\ F_{\text{\!e}}\vphantom{\frac{\partial F_{1}}{\partial_{p}}},$$

Full text

*Aims. * We focus here on studying two prominent features: long fibrils and flare ribbons. To model them, we use a 3D MHD simulation of an active region which self-consistently reproduces both of them.

*Methods. * We model the H$\alpha$, Mg II k, Ca II K, and Ca ii 8542 Å  lines using the 3D non-LTE radiative transfer code Multi3D. To obtain non-LTE electron densities, we solve the statistical equilibrium equations for hydrogen simultaneously with the charge conservation equation. We treat the Ca II K and Mg II k lines with partially coherent scattering.

*Results. *This simulation reproduces long fibrils that span between the opposite-polarity sunspots and go up to 4 Mm in height. They can be traced in all lines due to density corrugation. Opposite to previous studies, H$\alpha$, Mg II h&k, and Ca II H&K, are formed at similar height in this model. Although, some of the high fibrils are also visible in the Ca ii 8542 Å line, this line tend to sample loops and shocks lower in the chromosphere. Magnetic field lines are aligned with the H$\alpha$  fibrils, but the latter holds to a lesser extent for the Ca ii 8542 Å line.

The simulation shows structures in the H$\alpha$ line core that look like flare ribbons. The emission in the ribbons is caused by a dense chromosphere and a transition region at high column mass. The ribbons are visible in all chromospheric lines, but least prominent in Ca ii 8542 Å line. In some pixels, broad asymmetric profiles with a single emission peak are produced, similar to the profiles observed in flare ribbons. They are caused by a deep onset of the chromospheric temperature rise and large velocity gradients.

*Conclusions. *The simulation produces long fibrils similar to what is seen in observations. It also produces structures similar to flare ribbons despite the lack of non-thermal electrons in the simulation. The latter suggests that thermal conduction might be a significant agent in transporting flare energy to the chromosphere in addition to non-thermal electrons.

††offprints: J. P. Bjørgen,

1 Introduction

The chromosphere above active regions exhibits various dynamic phenomena. The most prominent features in quiescent periods of active region evolution are long fibrils. Their existence was revealed already with the early solar observations that sampled H$\alpha$ line center (Beckers 1964; Marsh 1976). Later, with improved spectral and spatial resolution, a similar scene was found in other chromospheric lines such as Ca ii 8542 Å  and Ca II K (Pietarila et al. 2009; Reardon et al. 2011; Robustini et al. 2019).

It is still, however not clear how they come about. Are they visible in chromospheric lines because of the density corrugations (Leenaarts et al. 2012) or do they outline magnetic tangential discontinuities where current sheets are formed (Judge et al. 2011)? The most recently proposed scenario is that the long fibrils are actually contrails where hot material shoots up as a consequence of energetic events happening lower down (Rutten & Rouppe van der Voort 2017). It is also unclear if and how well fibrils follow the magnetic field lines (de la Cruz Rodríguez & Socas-Navarro 2011; Asensio Ramos et al. 2017; Schad et al. 2013)? Answering this question is crucial for chromospheric seismology (e.g. Kuridze et al. 2012) as well as techniques to constrain extrapolation of the magnetic field from photospheric magnetograms (e.g. Aschwanden et al. 2016).

Another prominent feature in the chromosphere during the flaring phase in active regions, are flare ribbons. According to the standard flare model (Benz 2008; Fletcher et al. 2011), these are locations where energy is deposited in the chromosphere predominantly by beams of energetic electrons precipitating down from the coronal reconnection site. Strong chromospheric lines typically show broad emission peaks at those locations (Kerr et al. 2015; Rubio da Costa et al. 2015; Tei et al. 2018; Panos et al. 2018) which models have difficulty to reproduce.

The main problem in modeling both long fibrils and flare ribbons is that the radiation originating from the solar chromosphere forms under complex non-Local Thermodynamic Equilibrium (non-LTE) conditions. On the one hand, this makes observations of the chromosphere difficult to interpret; on the other hand, the forward synthesis of chromospheric lines becomes quite complex and computationally expensive (Leenaarts et al. 2012; Štěpán et al. 2015; Schmit et al. 2017; Leenaarts et al. 2013b; Bjørgen et al. 2018; Sukhorukov & Leenaarts 2017) if one wants to reproduce the spectral lines in detail.

The backbone of theoretical studies of 3D line formation in the solar chromosphere was so far the ”enhanced network” run (Carlsson et al. 2016) computed with the 3D radiation-MHD code Bifrost (Gudiksen et al. 2011). A simulation including flux emergence performed with the same code has also been used (Hansteen et al. 2017).

In case of flares, modeling of chromospheric lines has so far been limited to studies using 1D hydrodynamic models with heating by a prescribed electron beam (e.g. Allred et al. 2005; Rubio da Costa & Kleint 2017). Line profiles from these simulations fail to reproduce the line widths. However they can sometimes be achieved by adding artificial large velocity gradients to the simulation results. Some observed profiles, typically in Mg II h&k, have a single emission peak with a very wide symmetric base. Such profiles cannot be simulated by velocity gradients or microturbulence only (Rubio da Costa & Kleint 2017).

In this paper, we use a 3D MHD model that self-consistently shows both long fibrils and flare ribbons. It was computed using the radiation-MHD code MURaM (Vögler et al. 2005; Rempel 2017). The computational domain spans over an entire active region, and contains a bipolar sunspot pair. The run produces a magnetic reconnection event in the simulated corona resembling observational characteristics of a flare and an eruption of cool chromospheric material (Cheung et al. 2018).

The radiative transfer computations were performed in full 3D including partially coherent scattering (PRD) as well as charge conservation. (which is particularly important for hydrogen and has a clear effect in the Ca ii and Mg ii intensities). We investigate in detail the formation of the Ca II K/ $8542$ Å, Mg II k, and H$\alpha$ lines, all important diagnostics of the chromosphere. We are motivated to do this because we want to test to what extent the MURaM simulation produces a chromosphere that resemble real observations, specifically, we are interested in whether the large spatial scales lead to long chromospheric fibrils, and how the coronal flare affects the chromosphere, despite the lack of non-thermal electrons in the model.

We compare the synthetic observations with observations taken with the Swedish 1-meter Solar Telescope (Scharmer et al. 2003, SST;) and the Interface Region Imaging Spectrograph (De Pontieu et al. 2014).

The paper is structured as follows: the modeling and methods are presented in Section 2, an overview of the observations in Section 3, the results in Section 4, and finally the discussion and conclusions are given in Section 5.

2 Modeling

2.1 Model atmosphere

We ran a simulation with the same setup as in (Cheung et al. 2018), but with double the spatial resolution. We use a snapshot from this simulation taken at $71$ seconds before the flare peak, which is defined as the peak of the synthetic soft X-ray flux as it would be measured by the GOES satellite.

The simulation spans in the vertical direction from $-$7.5 Mm below the photosphere to 41.6 Mm above it and includes the top of the convection zone, the photosphere, chromosphere and a part of the corona. The simulation box has $1024\times 512\times 1536$ grid points, which correspond to $98.304\leavevmode\nobreak\ \text{Mm}\times 49.152\leavevmode\nobreak\ \text{Mm}\times 49.152\leavevmode\nobreak\ \text{Mm}$ in physical size, with a uniform horizontal grid of 96 km and a uniform vertical grid of 32 km. The simulation setup contains a bipolar active region that was present in the initial state. Near one of the sunspots a strongly twisted small bipolar magnetic structure was emerged by feeding in horizontal flux at the bottom boundary located about 7.5 Mm beneath the photosphere.

The simulation is run with gray LTE radiative transfer in the photosphere and chromosphere, optically thin radiative losses from the transition region and corona and included heat conduction along the magnetic field. The EOS uses a realistic solar composition and assumes local thermodynamic equilibrium.

The primary goal of this numerical experiment was to simulate the evolution of coronal magnetic field self-consistently in 3D in a flare. For this reason, the processes relevant for chromospheric layers were not accounted for in detail. The gray approximation is chosen because it resulted in a chromospheric temperature stratification which agreed with 1D standard semi-empirical models of solar atmosphere much better than the standard LTE multi-group opacity scheme with four frequency bins (Nordlund 1982). Furthermore, the opacity and source function are set to zero at the base of the transition region, or the first grid cell with $T>20$ kK, in order to avoid contributions from the corona in the downward directed rays.

A slope-limited diffusion scheme was chosen so that overshoots are minimized in locations where the gradients can be significant, like near the flare ribbons. The problem is however not completely removed. For example, it results in cold pockets at the base of the transition region in some areas of the computational domain. Finally, to relax the time-step constraints, the Boris correction and hyperbolic heat conduction are used. For more details about the code, see Rempel (2017).

Figure 1 displays the model atmosphere. Panel (a) shows the intensity in the continuum at 5 000 Å, with granulation and a pair of sunspot-like structures (but note the lack of penumbra around the spots), separated by 50 Mm. Panel (b) displays the vertical field strength, revealing that the upper sunspot is in fact a $\delta$-spot. In panel (c) we show the lowest height where the temperature is 50 kK, and can be thought of as a map of the height of the onset of the transition region. It shows elongated structures that reach up to 8 Mm and connect the sunspot pair.

In the model atmosphere, there is a long arm of material at chromospheric temperatures that has erupted up to a height of 18 Mm into the corona, reminiscent of a coronal mass ejection. For brevity we refer to this arm as CME.

Since the CME has chromospheric temperatures, it should ideally be included in our radiative transfer calculations. We choose not to do so because this would increase the computation time by a factor of 3, which we could not afford. Fortunately, the CME covers only 2.3% of the surface area of the simulation domain. The surface area of the removed CME arm is shown as red patches in Fig. 1c.

We removed the upper part of the CME, so that the model atmosphere size is reduced to $1024\times 512\times 290$ grid points, corresponding to heights from $-$1.3 Mm beneath the photosphere up to 8 Mm above it. The zero point of our height scale is defined as the height where the average optical depth at 5 000 Å is unity in a patch of the model atmosphere without strong vertical field strength. Due to the Wilson depression, there are columns where the height of optical depth unity at 5 000 Å is below the zero point of our height scale. The largest offset is $-$0.9 Mm located in the upper sunspot in Fig. 1.

For our non-LTE calculations with hydrogen we resampled the vertical resolution of the model from 32 km to 16 km around the transition region to better resolve large gradients. This resulted in a grid size of $1024\times 512\times 499$ points. For magnesium and calcium, we kept the original grid spacing because the spectra of these atoms are formed below the transition region.

In the MURaM simulation there exist locations where the temperature is as low as 1 100 K. Our background opacities are not correct at such low temperatures, and therefore we set the minimum temperature to 3 250 K in the model atmosphere. This has a minimal effect on the resulting spectra.

The simulation used a LTE equation of state. The electron densities in the real chromosphere are however far from LTE, and should ideally be computed using non-equilibrium ionization of at least hydrogen, and ideally also helium (e.g., Leenaarts et al. 2007). Computing the non-equilibrium ionization balance in post processing is not possible. Therefore, we chose to compute the electron density from non-LTE statistical equilibrium for the contribution from hydrogen, and in LTE from all other elements since they add a minimal number of electrons to the total in the upper chromosphere. The details for this procedure are given in Section 2.4.

2.2 Radiative transfer computations

We numerically solve the non-LTE radiative transfer problem with the Multi3D code (Leenaarts & Carlsson 2009) in a model atmosphere discretized on a Cartesian three-dimensional (3D) grid. The code treats one atom in non-LTE by simultaneously integrating the transfer equation and solving the statistical equilibrium equations for the atomic level populations. This is done iteratively using the multi-level accelerated $\Lambda$-iteration (M-ALI) method with preconditioned radiative rates following Rybicki & Hummer (1991, 1992). The transfer equation is integrated at spectral points covered by radiative transitions of the atom using a domain-decomposed and parallelized method of short characteristics (e.g., Olson & Kunasz 1987). For angle integrations we use a quadrature with 24 rays (set A4 from Carlson 1963). The opacities and the source function are approximated along the characteristics using linear interpolation, which is the most stable but the least accurate method. Our criterion for convergence is when the relative change in the populations is ${\leq}10^{-3}$.

In order to have a more accurate electron density, we first solve the non-LTE problem for hydrogen, where the statistical equilibrium equations are additionally constrained by the charge conservation equation (see Section 2.4). The obtained non-LTE electron density is re-used for solving the non-LTE problem for Ca II and Mg II.

We treat the Ca II K and the Mg II k lines in partial redistribution following the method by Sukhorukov & Leenaarts (2017).

2.3 Model atoms

We used three different model atoms to compute the various spectral lines from the model atmosphere.

To model H$\alpha$ we used a three-level plus continuum model atom of H I, that was constructed by removing the $n=4$ and $n=5$ levels from the atom used by Leenaarts et al. (2012). Following them, we treated Ly$\alpha$ in complete redistribution (CRD) with a Doppler absorption profile, to avoid the large computational cost of treating this line in partial redistribution (PRD).

First, we experimented with the original five-level plus continuum model atom, which produced population inversions in the Ly$\gamma$ and Ly$\delta$ transitions at several places of the model atmosphere. Such inversions cause negative opacities as well as negative source functions. Although the analytical formal solution of the transfer equation allows this, such inversions often produce numerical instabilities, especially, in the method of short characteristics. Algorithms used for interpolating the source function and the opacity either to facets of grid voxels or along the short characteristic, have inherent numerical errors, which might produce at the same time opposite signs of the source function and optical depth. Such an unphysical situation destroys the convergence of the M-ALI iteration scheme.

We tried the method of Koesterke et al. (2002), who proposed to branch the integration of the transfer equation in terms of emissivities and opacities instead of source functions and optical depths, but found this to be unstable too. Finally, we simply removed the upper levels of the transitions that had population inversions from the model atom.

For Mg II we used a four-level plus continuum model atom of Leenaarts et al. (2013a), where we treat only the Mg II k line in PRD and Mg II h in CRD to lower the computational costs.

Finally, for Ca II we used the five-level plus continuum model atom of Bjørgen et al. (2018) where we similarly treat only Ca II K in PRD and Ca II H in CRD.

We tested using the 1D FAL-C model atmosphere (Fontenla et al. 1993) to see whether the neglection of PRD in one of the doublet’s transitions changes the emergent intensity in the other line. In Ca II K only the inner core was marginally affected. In Mg II k the core was changed by 1.7% and the emission peaks were changed by less than 8.5%. In each doublet, Ca II K and H as well as Mg II k and h, both lines share similar formation properties, so in the following we focus on the Ca II K and the Mg II k lines only.

The emergent intensities as well as the formation heights slightly change owing to all these simplifications in our model atoms. These changes are relatively small. Since we are not concerned about detailed properties, they do not change the overall results and conclusions of the paper.

2.4 Charge conservation

A correct electron density is important for non-LTE modeling of chromospheric lines because the collisions with electrons affect the coupling of the source function to the local conditions of the gas. The model atmosphere does not explicitly contain the electron density, but MURaM implicitly assumes LTE electron densities through its equation-of-state. LTE is however a bad approximation of the electron density in the chromosphere and transition region, and the electron density should ideally be computed including non-equilibrium ionization of both hydrogen and helium (Carlsson & Stein 2002; Leenaarts et al. 2007; Golding et al. 2016).

Computing non-equilibrium ionization in post-processing is not possible, so instead we derive the non-LTE electron density through solving the non-LTE problem for hydrogen together with the charge conservation equation.

First, we tried an approach similar to what the RH code uses (Uitenbroek 2001). The electron density is updated after every M-ALI iteration to be consistent with the proton density, while the ionization of all other elements is computed in LTE. We found this approach to unstable in our calculations and the NLTE populations for hydrogen did not converge for our problem.

Instead, we decided to solve simultaneously in each M-ALI iteration the equation of charge conservation and the system of statistical equilibrium equations for preconditioned rates. This new system of equations is non-linear with respect to the unknown electron density and hydrogen level populations as some of the rates depend on the electron density. We solve it iteratively using the Newton-Raphson method. A similar technique has been used in other radiative transfer codes (e.g., Carlsson & Stein 1992; Paletou 1995; Heinzel 1995). Our implementation is explained in Appendix A.

We note that this is a substantial improvement over previous radiative transfer calculations with the Multi3D code for hydrogen, where the proton density could be larger than the electron density (Leenaarts et al. 2012; Hansteen et al. 2017).

2.5 Broadening of H$\alpha$

The default Stark broadening mechanism for H$\alpha$ in Multi3D is based on the approximation of Sutton (1978). When the electron density exceeds $10^{13}$ cm*-3* in the chromosphere, this approximation becomes inaccurate. We tested whether we should include the more accurate unified Stark broadening theory as recently proposed by Kowalski et al. (2017) to model the line profile. We performed 1D tests with several columns from the model atmosphere. For columns with chromospheric electron densities above $10^{13}$ cm*-3* we indeed found that the line profile computed with the unified theory broadens more. However, most profiles do not show a substantial difference since the typical electron densities are not as high in the line formation region. Since the unified theory is computationally expensive owing to a convolution operation, we decided not to include the unified Stark broadening theory for H$\alpha$.

2.6 H$\alpha$ in 1D and 3D

Previous H$\alpha$ studies using a 3D evaluation of the radiation field were performed using quiet-Sun conditions (Leenaarts et al. 2012). Under those circumstances, the intensities computed in 1.5D, where the effect of horizontal scattering is neglected, showed a strong imprint of photospheric granulation. A 3D evaluation of the radiation field was needed to obtain chromospheric features in the emerging intensities. In the present study, we are analyzing an active region, and we wanted to assess the importance of horizontal scattering compared to the quiet-Sun case.

Figure 2 shows the comparison between the vertically emergent intensity in 1D and 3D using our active region snapshot. The 1D computation shows chromospheric structures and traces of the large-scale structures, but much less enhanced than in the 3D case. We note that an imprint of the photosphere is still visible through the fibril-like structures in the 1D case, but not as prominent as in the quiet-Sun case.

The brightest features have a similar morphology in the 1D and 3D computations. Analysis shows that this is caused by the high mass density, and thus strong collisional coupling of the source function to the temperature, in the chromosphere at these locations (see Sec. 4.6.1). Therefore, this can potentially open the possibility to include H$\alpha$ in inversions that assume 1D geometry (Milić & van Noort 2018; de la Cruz Rodríguez et al. 2018) in active regions.

3 Observations

We use two different observational datasets to compare with the synthetic data. The first dataset is an observation of active region NOAA 12593, located near to the disk center at $\theta_{x}=-68$″and $\theta_{y}=-6$″(helioprojective-Cartesian coordinates), taken at the Swedish 1-meter Solar Telescope (Scharmer et al. 2003, SST;) on September 19, 2016, at 09:31–09:57 UT.

The data was recorded simultaneously by the CHROMIS and CRISP (Scharmer et al. 2008) instruments, which observed Fe i 6301 Å, Fe i 6032 Å, H$\alpha$, Ca ii 8542 Å, and Ca ii K. We only use the latter three spectral lines. The data were processed with the CRISPRED and CHROMISRED pipelines (de la Cruz Rodríguez et al. 2015; Löfdahl et al. 2018). For a more detailed description of the SST observation we refer to Leenaarts et al. (2018), here it suffices to mention that we obtained near-simultaneous spectral scans in all three lines.

Unfortunately, there is no co-observation with the Interface Region Imaging Spectrograph (IRIS De Pontieu et al. 2014). Therefore, we use IRIS data observed on February 5, 2016, at 08:21–09:12 UT in order to be able to compare our computations of the Mg II k line. The IRIS target is an emerging active region NOAA 12494, observed near the disk center at $\theta_{x}=-96$″and $\theta_{y}=-98$″.

Figure 3 shows the line core images for all four spectral lines. Panels a-c show an active region covered with fibrils and a sunspot in the upper right corner. Panel d shows a slit jaw-image at the nominal line center of Mg II k, of an active region with one sunspot located on the left side and a sunspot outside the FOV on the right side. Long elongated fibrils are connected with the sunspots and pores located in the middle of the FOV.

4 Results

4.1 The effect of non-LTE electron density on the emergent intensity

We compared how different are electron densities computed either using the LTE equation-of-state for all chemical species or by treating hydrogen in non-LTE with imposed conservation of charge. We found severe differences between the LTE and non-LTE values in the middle and upper chromosphere.

On average the inclusion of charge conservation in the calculations, leads to higher electron densities compared to LTE in cold areas of the chromosphere, owing to an increased hydrogen ionization degree through the Balmer continuum. This difference was largest (by a factor $\sim 10^{3}$) in cold pockets (3250 K) right below the transition region. In hot areas in the mid chromosphere we typically find a somewhat lower electron density than in LTE.

Figure 4 illustrates how the non-LTE electron densities change the emergent intensity in Mg II k. For the spatially-averaged intensity, the maximum relative difference reaches 30% near the line core (Fig. 4a). Individual profiles show even larger differences depending on which chromospheric structure they are emerging from. In dark fibrils the non-LTE intensity drops below the LTE one (Fig. 4b), but in flare ribbons it increases significantly near the line core (Fig. 4c).

Therefore, in the rest of the paper we only present results using the non-LTE electron density.

4.2 Formation heights

Calcium and magnesium are both predominantly singly ionized in the chromosphere. Mg II ionizes to Mg III at slightly higher temperatures owing to the higher ionization potential of Mg II (15.0 eV) compared to Ca II (11.87 eV). Both the Mg II h&k and Ca II H&K lines have the ground state as their lower level, so their opacity scales approximately with the density as long as Mg II and Ca II are the dominant ionization stage. As Mg II is 18 times more abundant than Ca II, the line cores of Mg II h&k are formed a few scale heights higher up in the chromosphere than the line cores of Ca II H&K (Leenaarts et al. 2013b; Bjørgen et al. 2018).

In contrast, the lower level of H$\alpha$ is an excited state at 10.2 eV above the ground state. In LTE this makes the line opacity very sensitive to the temperature. In the solar chromosphere, the H$\alpha$ line opacity depends in a complicated way on the temperature, density, and the time history of the atmosphere owing to the non-equilibrium ionization of hydrogen. Leenaarts et al. (2012) investigated the formation of H$\alpha$ in a radiation-MHD simulation of the quiet Sun and found that in those circumstances the line can be reasonably accurately modeled assuming statistical equilibrium. In that simulation they also found that the line opacity weakly depended on temperature and mainly depended on mass density. Carlsson & Leenaarts (2012) investigated the ionization of H I, Ca II, and Mg II, and found that at temperatures above 15 kK, the fraction of H I is much larger than the fractions of Ca II and Mg II. Rutten & Rouppe van der Voort (2017) and Rutten (2017) argued similarly, albeit based on LTE, that H$\alpha$ can have an appreciable opacity at temperatures well above 20 kK.

Studies of the formation heights of the line cores based on the quiet-Sun-like simulation of Carlsson et al. (2016) indicate that H$\alpha$ forms 300–1000 km below Ca II K and Mg II k (Leenaarts et al. 2013b). The similarity of the line-core images in Fig. 3a,b casts doubt on the validity of this result in active regions. We therefore investigated the formation heights in the MURaM simulation.

We measured the formation heights at optical depth unity in the corresponding line cores of H$\alpha$, Mg II k, Ca II K, and Ca II 8 542 Å. Figure 5 shows differences between the formation height of H$\alpha$ and the formation heights of the rest of the lines. On average, the core of H$\alpha$, k3, and K3 features are formed at similar heights and roughly 150 km above the core of 8 542 Å throughout the field-of-view.

We note that the distributions have long tails towards positive values meaning that H$\alpha$ can be formed much higher than the other lines. These points, where the height difference exceeds 1 000 km, constitute for the field-of-view ${\sim}$2% for both Mg II k and Ca II K, and $\sim$5% for 8 542 Å. Figure 6 illustrates this case. At 3.22 Mm where the H$\alpha$ core reaches optical depth unity, H I is only 97.5% ionized and at the same location the metals (Ca II and Mg II) are 99.9% ionized. H$\alpha$ is formed highest in the atmosphere. It retains opacity where the lines from the metals have none because H I does not ionize as quickly as Ca II and Mg II when the temperature increases.

4.3 Synthetic images

In Fig. 7 we show the vertically emergent intensity in the line core of Ca II 8542 Å, Ca II K, Mg II k, and H$\alpha$. We first describe what we see at nominal line center (Fig. 7, top panels). In all four spectral lines we see large elongated structures emanating from the two sunspots (see Fig. 1b), similar to large-scale fibrils seen in observations. From now on we call the structures on our synthetic images fibrils too. The fibrils are much more pronounced in H$\alpha$ than in the other lines. The longest fibrils extend up to 35 Mm.

Careful comparison of the panels reveals that the fibrils appear somewhat different in Ca II 8 542 Å, looking translucent with the background chromosphere shining through. This is especially visible in the upper-right corner. This is consistent with the lower formation height of this line, as discussed in Section 4.2.

In H$\alpha$, we observe several very bright patches close to the sunspots, which connect with a thinner bright arc. They look similar to flare ribbons. Along these ribbons the H$\alpha$ profiles are in emission and very broad. In Ca II K and Mg II k we see only a weak signature of the ribbon, while in Ca II 8 542 Å it is only visible as a small brightness enhancement to the right of the upper sunspot.

Now we describe what we see at wavelengths where the intensity in each pixel individually is taken at the wavelength where the formation height is largest. (Fig. 7, bottom panels). These images thus compensate for Doppler shifts of the line core (Leenaarts et al. 2013b; Bjørgen et al. 2018).

The H$\alpha$ image is relatively unchanged, while the other images are more affected. This difference is caused by the large thermal broadening of H$\alpha$ compared to the other lines, which makes it less sensitive to Doppler motions. In the Ca II and Mg II lines the carpet of fibrils now appears thicker and sharper than in the nominal line-core images. The flare ribbons are now visible in these lines, even in Ca II 8 542 Å.

The simulated images exhibit features that quite resemble observational features in our Fig. 3 or Fig. 1 by Reardon et al. (2011). The observations show that Ca II K and H$\alpha$ trace roughly the same structures. The Ca II 8 542 Å often shows different fibrils, with more visible patterns of waves and shocks from the lower chromosphere.

The SST/CRISP, SST/CROMIS, and, to a lesser extent, IRIS observations show finer structures than the synthetic images. This is most likely caused by the finite horizontal grid spacing (96 km) of the model atmosphere compared to the resolution of the observations ($\sim$100 km for SST, $\sim$250 km for IRIS). Bjørgen et al. (2018) showed for radiation-MHD models of the enhanced network that synthetic fibrils get thinner and closer in appearance to the observed ones if the horizontal grid spacing of the simulation decreases from $48$ km to $31$ km.

Leenaarts et al. (2012, 2013b) and Bjørgen et al. (2018) computed synthetic images of the chromosphere in H$\alpha$, Ca II K, and Mg II k from the quiet Sun simulation by Carlsson et al. (2016). Their line-core images showed ordered fibrils stretching between two opposite magnetic polarities separated by about 8 Mm, which also limited the fibril length to 8 Mm. Hansteen et al. (2017) computed H$\alpha$ and Mg II h for a flux emergence simulation, also computed with the Bifrost code. Their H$\alpha$ line-core image showed some thick unordered fibrils with a maximum length of 12 Mm, while the corresponding Mg II h image does not show a clear fibril structure.

The current MURaM simulation is the first simulation including a chromosphere that contains the spatial scales of an entire active region, and succeeds in reproducing long fibrils connecting the opposite polarities of an active region as seen in observations.

In the following subsection we analyse some of the structures seen in the synthetic intensity images in more detail.

4.4 Synthetic fibrils

The horizontal domain size of the current MURaM simulation is 98 Mm${}\times{}$49 Mm. This large size allows for $\sim$40 Mm long magnetic loops in the chromosphere that connect the two sunspots. The synthetic images in Fig. 7 show long fibrils that follow the large-scale structure of the field. We investigated whether the fibrils trace the magnetic field, both the vertical and horizontal components.

Figure 8 shows the line-core intensity of H$\alpha$, with the horizontal component of the magnetic field at the formation height overplotted. The fibrils follow the horizontal magnetic field direction in most places, but at some locations (such as the arrow at $(X,Y)=(54,34)$ Mm) they do not.

Figure 9 illustrates the magnetic field component in the plane of a model atmosphere slice that follows the blue line in Fig. 8. As expected for active regions, the magnetic pressure exceeds the gas pressure at low heights in the atmosphere, as illustrated by the plasma $\beta=1=P_{\text{gas}}/P_{\text{B}}$ curve, which goes from $-$1 Mm to $0.5$ Mm. The lowest height of $\beta=1$ curve is located at the sunspot ($X>69$ Mm), where the magnetic field strength reaches 6 kG in the photosphere. Horizontally the fibril spans from $X=51$ Mm to $X=68$ Mm and for $X<51$ Mm the cut intersects an area of the atmosphere that shows irregular short structures in the H$\alpha$ image of Fig. 8.

In the sunspot, the transition region is very low, 1 Mm above the quiet photosphere and roughly 2 Mm above the local height of continuum optical depth unity. so that all four spectral lines are formed mostly at the same height. The angle between the field vector and the $\tau=1$ surface is large. In the area where $x<51$ Mm we see intrusions of high mass density sticking into the low-density corona, and the $\tau=1$ surface follows these mass intrusions. We speculate that these could be Type I spicules as simulated before in 2D (e.g., Hansteen et al. 2006), but this is impossible to confirm without analyzing a time series.

Now we turn our attention to the long loop-like structure between $X=51$ Mm and $X=68$ Mm in Figure 9. The magnetic field vector and the $\tau=1$ height curve of H$\alpha$ are parallel to each other for a large fraction of the loop. It is not shown in the figure, but this is true to a lesser extent for the other lines, which all exhibit jumps in their formation height. We analyze a number of other fibrils not illustrated here and found that for clearly defined fibrils, H$\alpha$ indeed tends to follow the magnetic field lines. This is in contrast with the results from Leenaarts et al. (2015), who found in a simulation of quiet sun fibrils that H$\alpha$ fibrils typically follow the horizontal component, but not the vertical component, of the magnetic field.

To investigate the magnetic field we traced two fibrils to get a more detailed overview. They are marked by blue and red curves in Fig. 8. Figure 10 shows the offset between the inclination of the $\tau=1$ curve and magnetic field vector, and the offset between the azimuths of the fibril and magnetic field. We only show results for H$\alpha$ and Ca II 8 542 Å, as they are the lines with the highest and lowest formation heights.

Figure 10a shows the same fibril as in Fig. 9 inbetween 51 Mm and 68 Mm. The inclination offset is around $10^{\circ}$ for both lines. However, the formation height curve of H$\alpha$ is more smooth, while the Ca II 8 542 Å curve has a number of sudden jumps, which cause spikes in the offset curve. Physically, these jumps mean that depending on where you look along the fibril, you see different field lines, but locally the formation height curve follows the magnetic field inclination. The H$\alpha$ fibril follows the azimuth very well, while the Ca II 8 542 Å shows offsets up to 10°.

Figure 10b shows the offsets of a curved fibril indicated with a red line in Fig. 8. We see almost the same behavior as the previous fibril shows: in H$\alpha$ it follows both the inclination and the azimuth of the magnetic field vector very well, while in Ca II 8542 Å it does it less so owing to large changes in the formation height.

In addition to the fibrils analyzed here we manually investigated a number of other fibrils, and they show similar behavior. We thus conclude that most fibrils in this active region simulation do trace the horizontal magnetic field, especially for spectral lines formed in the upper chromosphere (Ca II H&K, Mg II h&k, and H$\alpha$). Fibrils seen in 8 542 Å mostly, but not always, trace the horizontal component of the magnetic field. This agrees with the observational studies (de la Cruz Rodríguez & Socas-Navarro 2011; Asensio Ramos et al. 2017). We note that this simulation does not include ambipolar diffusion. Inclusion of this process allows the magnetic field to be less aligned with density and temperature structures in the chromosphere because it relaxes the frozen-in condition (Martínez-Sykora et al. 2016). However, simulations with ambipolar diffusion have only been performed in 2.5D simulations (Cheung & Cameron 2012) or simplified 3D simulations (Arber et al. 2007) or radiative 3D simulation with only photosphere (Danilovic 2017). Therefore, only a full 3D MHD simulation that includes chromosphere with ambipolar diffusion should be used to investigate the possible misalignment properly (see also Sec. 3.8.3 of Cheung & Isobe 2014).

4.5 Spatially averaged spectral profiles

We compared the spatially-averaged synthetic spectra with observations from SST and IRIS to see how far the model is from the real Sun. Although the solar features represented in the observed and simulated samples are quite different, the comparison allows us also to place this model in the context of previous studies (e.g Leenaarts et al. 2013b; Bjørgen et al. 2018). Figure 11 shows the spatially-averaged profiles of Ca II 8542 Å, Ca II K, Mg II k, and H$\alpha$ from the model atmosphere and observations (see Section 3). Two different averages are shown for the synthetic spectra: the entire model atmosphere (panels a–d), and in a patch with relatively weak magnetic field (panels e–h). We note that the observations are not in any way chosen to represent all active regions, they just illustrate how average line profiles in active regions look like.

The observations of H$\alpha$ and Ca II 8542 Å show absorption profiles, while Mg II k, and Ca II K show a double-peaked central reversal (marginally visible in Ca II K owing to the spectral resolution of SST/CHROMIS). The synthetic spectra for the atmosphere outside the flare ribbons show qualitatively the same behavior, but substantial differences in a quantitative sense.

The averaged spectra from the model atmosphere gives a peak separation of 16 km s*-1* (0.15 Å) for Mg II k and 10.5 km s*-1* (0.14 Å) for Ca II K. These numbers are $\sim$30% larger than those reported for Bifrost simulations of the quiet Sun (Leenaarts et al. 2013b; Bjørgen et al. 2018), but still $\sim$1.5–2 times smaller than observed for the quiet Sun.

The averaged profiles of Ca II K and 8 542 Å show that the model atmosphere is colder in the lower chromosphere and photosphere than the observations. The synthetic line wings of Mg II k have a higher intensity than the observations. We also note that the inner wings of the Ca II 8542 Å, Ca II K, and H$\alpha$ line lie below the observed intensity, while the Mg II k wings lie above it.

Inspection of Fig. 11d shows that the averaged H$\alpha$ profile has a central emission and a weak and blue-shifted central absorption core, in stark contrast to the observations. Figure 11h shows that the core behaves as usual in the less magnetized patch of the simulation. We show in Section 4.6.1 that the filling of the H$\alpha$ core in the averaged profile is caused by a strong H$\alpha$ emission in the flare ribbons.

4.6 Spatially resolved synthetic profiles

In Figs. 12–15 we explain the formation of line intensity profiles in three different structures in the simulation: flare ribbons (with and without strong velocity gradients), a fibril, and a sunspot umbra. The top row of each figure shows the profiles for each line, while the bottom row shows the temperature and line source functions for H$\alpha$ and Ca II 8542 Å in the left panel, the two-level photon destruction probability in the middle panel, the height where $\tau_{\nu}=1$ for each line together with the vertical velocity in the right panel. The photon destruction probability is defined as

[TABLE]

with the downward collisional rate $C_{\!ji}$, the Einstein coefficients $A_{\!ji}$ for spontaneous deexcitation, $B_{\!ji}$ for induced deexitation, and the Planck function $B_{\nu}$. We do not display the line source functions for Ca II K and Mg II k because they are frequency-dependent owing to PRD effects. We show $\epsilon$ because it gives an indication of the sensitivity of the line source function to the local temperature using the two-level source function approximation: $S_{\nu}=(1-\epsilon)\bar{J}_{\nu 0}+\epsilon B_{\nu}$.

4.6.1 Flare ribbons

We investigated our model atmosphere and line formation in the flare ribbons, and show two columns in Figs. 12–13. The upper panels of Fig. 12 show an example of broad, slightly redshifted emission peaks in all four lines. The chromosphere has a high mass density as evidenced by the large values of $\epsilon$ (which is roughly proportional to the electron density). The transition region lies at a high column mass of $10^{-4.2}$ g cm*-2*, and the chromospheric temperature rise starts deep in the atmosphere at 0.2 g cm*-2*. This partly explains why the ribbons look so similar when the intensities are computed in 1D and 3D (see Fig. 2). The strong chromospheric temperature rise leads to a source function increasing with height for all lines, leading to the emission profiles, while the deep location of the onset of the rise causes the wide symmetric base of the emission peaks, especially for Mg II k. The downflow just below the transition region causes the redshift of the emission peak.

These relatively smooth flare ribbon profiles are actually rather rare in our FOV. We often find more complex spectral profiles, such as shown in Fig. 13. They typically show multiple-peaked profiles in Ca II K and Mg II k, with strongly Doppler-shifted emission peaks. The H$\alpha$ profile is very broad and has its strongest emission peak blueshifted by 130 km s*-1*. The temperature in the chromosphere shows multiple peaks, with low temperatures between the peaks. The mass density, and thus the photon destruction probability are again high and the transition region is again located at a high column mass. The atmosphere differs from the one shown in Fig. 12 in that the transition region harbors a strong velocity gradient. The cores of Mg II k and Ca II K form completely in the chromosphere in this gradient. The strongest emission peak in H$\alpha$ reaches optical depth unity between the last chromospheric grid point with a temperature of 5.5 kK, and the first transition-region/coronal grid point which has a temperature of 450 kK (i.e., in the temperature gradient indicated by the nearly vertical black line in the lower-left panel of Fig. 13.

The finite vertical resolution of the model atmosphere produces the peculiar, multiple-peaked line shapes of Mg II k, and Ca II K. The difference in vertical velocity between adjacent grid points is larger than the Doppler width, and this causes spurious emission peaks, as explained by Ibgui et al. (2013). This is reflected in the peaks in the $z(\tau_{\nu}=1)$ curves in the lower right panel. The formation height is large at the velocities of the grid points (black filled circles), but drops to low values when more than a Doppler width away from the velocities at the grid points.

The multiple emission peaks are thus artefact produced by the limited spatial resolution in the radiative transfer computation. Had the velocity gradient been properly resolved, then these line profiles would show a broad asymmetric emission peak.

The H$\alpha$ profile does not show such strong emission peaks. The emission peak at $\Delta\varv=-120$ km s*-1* is, however, badly resolved. The frequency grid is not fine enough to resolve the right flank, and as stated before, the peak forms within a single grid interval, which constitutes the interface between the chromosphere and transition region (compare the brightness temperature of the emission peak of 13 kK to the source function in the lower right panel). In this interval the optical depth increases from $10^{-5}$ to 1.2. The contribution function is thus not properly resolved, and the resulting intensity depends on the chosen solution scheme. We investigated the source function, and exactly in this grid interval its changes from two-level behavior in the chromosphere to recombination-radiation-dominated in the transition region.

4.6.2 Sunspot

Figure 14 shows the line formation in a sunspot in the model atmosphere. The H$\alpha$ core is in emission, and the three other spectral lines show double-peaked profiles. This is common in the model atmosphere, we find only a few locations where H$\alpha$ is in absorption and the other lines have single emission peaks.

The temperature stratification as shown in Fig. 14 is common across the sunspots, from the temperature minimum the gas temperature monotonically increases towards the transition region. The H$\alpha$ line is in emission because the line core is formed in the steep temperature gradient of the transition region, with a high source function. The line cores of the other lines tend to be formed just below the transition region. The chromosphere has a temperature rise that is deep enough to cause emission peaks.

Comparing to many available observations (e.g. Panos et al. 2018; Maurya et al. 2013; de la Cruz Rodríguez et al. 2013; Felipe et al. 2011; Engvold 1967), we see that the simulated sunspot profiles are somewhat unusual: observed H$\alpha$ profiles are typically in absorption; the Ca II K and Mg II k profiles are typically, but not always, single peaked; Ca II 8542 Å is in absorption and only shows emission when an umbral flash passes through.

4.6.3 Fibril

Figure 15 shows the typical line formation in a fibril. Here H$\alpha$ is in absorption, Ca II K, Mg II k, and Ca II 8542 Å each have an asymmetric double-peaked emission core. We find that this situation is common for H$\alpha$, Ca II K, and Mg II k. The density in the chromosphere is relatively low, and $\epsilon$ is an order of magnitude smaller than in the flare-ribbon columns. The H$\alpha$ line core is formed below the transition region and therefore is not influenced by the increase in the source function there.

The Ca II 8542 Å line shows much more variation, ranging from pure absorption to wide and asymmetric profiles with three emission peaks. The line is mainly in absorption in fibrils in our data, while we also find some instances of profiles with a double-peaked central reversal as in Fig. 15.

5 Discussion and conclusions

Cheung et al. (2018) performed a radiation-MHD simulation of an active region including emergence of a new flux that produced a flare-like reconnection event in the corona. We took a snapshot from a rerun of their simulation with twice the original numerical resolution.

This numerical experiment, aimed only at investigating the corona, produced a unique data set suitable for investigating the chromosphere as well. The simulation domain is geometrically $4\times 2\times 3$ times bigger and models a much more active chromosphere than the commonly-used publicly-available Bifrost simulation (Carlsson et al. 2016),

The trade-off for getting the large domain size is that the grid spacing, 96 km horizontally and 32 km vertically, is relatively coarse. Radiative energy losses in the photosphere and corona are approximated using the gray LTE transfer. The absorption of coronal radiation in the chromosphere is not included. The equation of state assumes LTE, whereas hydrogen and helium both require non-equilibrium ionization (Leenaarts et al. 2007; Golding et al. 2014). There are no interactions between ions and neutral species, which significantly structure and heat the chromosphere (Martínez-Sykora et al. 2012, 2017). These simplifications have an impact on the structure of the chromosphere, both for the temperature and density in the model atmosphere. However, these simplifications do not affect the large-scale magnetic structure in the chromosphere. Neither do they distort the formation of flare ribbons, chromospheric evaporations through heat conduction, even though the exact location of the transition region changes if more accurate radiative losses are included.

In this paper we have taken one snapshot of this MURaM simulation and have solved numerically the problem of 3D radiative transfer for H I, Ca II, and Mg II in non-LTE including partially-coherent line scattering (PRD) in the Mg II k and Ca II K lines. We have corrected the electron density by treating hydrogen in statistical equilibrium together with an equation for charge conservations. On the radiative transfer side, this is a large improvement over the electron density evaluated in LTE. For the sake of numerical stability and computational speed, we have set the minimum temperature in the model atmosphere to 3 250 K, we have used a simple approximation for the pressure broadening of H$\alpha$, and we have limited the size of our H I model atom.

The finite grid spacing resulted in a poor sampling in both frequency and optical depth that led to artefacts in the line profiles when the line-forming region contains large gradients in temperature, velocity and/or density. While these minor weaknesses should be addressed in future work to allow detailed quantitative comparison with observations, they have no consequences for the qualitative results shown this paper.

We have shown that 3D non-LTE radiative transfer computations including charge conservation and PRD can handle a large variation in all physical quantities, in a different parameter regime (active region) than studied before (more quiet-sun-like, Carlsson et al. 2016).

The synthetic line-core images have shown long ($\sim$35 Mm) strands connecting the opposite-polarity sunspots, that resemble fibrils seen in observations. We have established that the strong fibrils seen in H$\alpha$ are mostly aligned with the magnetic field direction. Fibrils seen in the other lines are less aligned, especially those in Ca II 8 542 Å, which suffer from sudden changes of the height of optical depth unity so. This contrasts with the findings in Leenaarts et al. (2015), who found that in more quiet-Sun circumstances only a fraction of the fibrils seen in H$\alpha$ follow the same field line.

In 3D MHD simulations of quiet Sun, Leenaarts et al. (2012, 2013b) found that the cores of Mg II k and Ca II are both formed in the upper chromosphere in a wide range of heights following the transition region while the cores of H$\alpha$ and Ca II 8 542 Å are formed at lower heights, in the middle chromosphere.

We have found two important differences with this in the MURaM simulation. Firstly, the H$\alpha$ line is formed in a much larger range of heights. Secondly, the cores of H$\alpha$, Mg II k, and Ca II K are formed at very similar heights while the core of Ca II 8 542 Å is formed only some 150 km below them. This qualitatively agrees with observations of active regions, which show a very similar appearance in the cores of H$\alpha$ and Ca II K Leenaarts et al. (2018). Given that the model used in this paper manages to reproduce overal apparence in all chromospheric observables we synhesized, indicates that the model gives a decent representation of long solar fibrils.

The flare ribbons in this simulation are caused by thermal conduction only. They appear bright in the line cores because the thermal conduction and dense coronal loops above the ribbon cause a hot and dense upper chromosphere, with sufficient electron density to couple the source function to the local temperature. This mechanism is similar to what Carlsson et al. (2015) proposed for Mg II k in plages. We find synthetic line profiles that are broad, asymmetric, and with a single emission peak, similar to the profiles observed in flare ribbons.

The flare ribbons are poorly visible in Ca II 8542 Å, in contrast to observations (e.g., Kleint 2012). This agrees with the standard flare model that requires beams of non-thermal electrons or Alfvénic waves to deposit thermal energy in the middle chromosphere. Nevertheless, the visibility of the flare ribbons hints that the fraction of the flare energy transported down by thermal conduction is enough to significantly alter the structure of the chromosphere.

This study is just a first step to understand the chromosphere in active regions using 3D radiation-MHD simulations. In our study, we have only touched on the formation of the strongest chromospheric lines. Future studies should look into the formation of the lines in more detail.

Acknowledgements.

The Swedish 1-m Solar Telescope is operated on the island of La Palma by the Institute for Solar Physics of Stockholm University in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias. The computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at the High Performance Computing Center North at Umeå University and the PDC Centre for High Performance Computing (PDC-HPC) at the Royal Institute of Technology in Stockholm. JL was supported by a grant (2016.0019) of the Knut och Alice Wallenberg foundation. SD and JdlCR were supported by a grant from the Swedish Civil Contingencies Agency (MSB). This study has been discussed within the activities of team 399 ‘Studying magnetic-field-regulated heating in the solar chromosphere’ at the International Space Science Institute (ISSI) in Switzerland. JdlCR is supported by grants from the Swedish Research Council (2015-03994), the Swedish National Space Board (128/15). This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (SUNMAG, grant agreement 759548). AVS acknowledges financial support from the Spanish Ministry of Economy and Competitiveness (MINECO) under the 2015 Severo Ochoa Program MINECO SEV-2015-0548. The National Center for Atmospheric Research is sponsored by the National Science Foundation. We would like to thank Adam F. Kowalski and Joel C. Allred for providing us their modified RH code, which includes the unified theory of electric pressure broadening. JPB would thank Gregal J. M. Vissers for the help in calibrating the IRIS data, and Carolina Robustini for discussions.

Appendix A Charge conservation in non-LTE

In the Multi3D code, the solution of the system of statistical equilibrium equations and the integration of the transfer equation are iterated self-consistently using the multi-level approximate $\Lambda$-operator scheme (M-ALI). In the particular case described in this paper the electron density is neither known a priori nor provided by the MURaM model atmosphere and must be inferred along the solution of the radiative transfer problem. For this purpose, in each M-ALI iteration we sub-iterate a system of non-linear equations, which join statistical equilibrium and charge conservation equations together. So far the code can treat in kinetic equilibrium (non-LTE) only one chemical element of interest while the others are considered in LTE. Hydrogen is treated in non-LTE using a model atom that has three levels of H I with the principal quantum numbers $n=\{1,2,3\}$ and H II continuum for protons. Here is the numerical implementation of this scheme.

The hydrogen population densities $n_{i}$ with $i=\{1,2,3\}$ and the proton density $n_{\text{p}}$ satisfy the system of 3+1 equations of statistical equilibrium for H I,

[TABLE]

and for H II,

[TABLE]

where $R_{a,b}$ and $C_{\!a,b}$ correspondingly denote preconditioned radiative and collisional rates for a transition $a{\rightarrow}b$ (for details see Rybicki & Hummer 1991, 1992). Except for the bound-bound radiative rates, $R_{i,j}$, all other rates explicitly depend on the electron density. We indicated this by the parentheses, $(n_{\text{e}})$, and emphasized by separating the bound-free rates having one proton subscript p from the remaining bound-bound rates having both numerical subscripts $i,j=\{1,2,3\}$.

As this homogeneous system has a trivial (zero) solution, one of the equations must be replaced with the condition of particle conservation for the total number density of hydrogen,

[TABLE]

We substitute the equation for the level having the largest population. This makes the matrix of rates well conditioned.

If the electron density $n_{\text{e}}$ is known and fixed then Eqs. (2–4) are linear with respect to the unknown populations and can be solved with any appropriate numerical method for algebraic linear systems.

If the electron density is not known then all rates explicitly dependent on $n_{\text{e}}$ are unknown and the system (2–4) is non-linear with respect to the unknowns $n_{i}$, $n_{\text{p}}$, and $n_{\text{e}}$. We add one more equation to constrain $n_{\text{e}}$ by stipulating charge neutrality, that is, the net negative and positive charges must be equal,

[TABLE]

which we express as

[TABLE]

where for the atomic species $s$, $\alpha_{\!s}$ is the element-to-hydrogen relative abundance, and $f_{\!s,i}$ is the corresponding ion fraction in LTE at the ionization stage $i>1$ relative to the neutral stage $i=1$, which is obtained from the Saha-Boltzmann law depending on the electron density $n_{\text{e}}$ and temperature $T$.

We use the multi-dimensional Newton-Raphson method (see Press 2007, Sect 9.6) to solve this non-linear system of joined statistical equilibrium, particle, and charge conservation equations by linearizing it to $\@vec{\hat{J}}\delta\@vec{n}=-\@vec{F}$ or, in expanded form,

[TABLE]

where $\delta\@vec{n}=(\delta n_{1},\delta n_{2},\delta n_{3},\delta n_{\text{p}},\delta n_{\text{e}})^{\top}$ is the column vector of corrections to unknowns, $\@vec{F}=(F_{\text{\!1}},F_{\text{\!2}},F_{\text{\!3}},F_{\text{\!p}},F_{\text{\!e}})^{\top}$ is the column vector of free terms given by Eqs. (2), (3), and (6), and $J_{i,j}=\partial F_{\!i}/\partial n_{\!j}$ is the Jacobian matrix with corresponding partial derivatives. As we already explained, one of the $F_{\!1}$, $F_{\!2}$, $F_{\!3}$, or $F_{\text{\!p}}$ equations is replaced by Eq. (4). We analytically computed the partial derivatives in a fashion similar to Leenaarts et al. (2007, see Appendix A).

This linear algebraic system of five equations with five unknowns is solved for the corrections $\delta\@vec{n}$ using a standard Gaussian elimination scheme and the unknowns are updated to $\@vec{n}+\delta\@vec{n}$. The initial values of $\@vec{n}$ are taken from the previous M-ALI iteration and are used to initialize $\@vec{F}$ and $\@vec{\hat{J}}$. The process is sub-iterated until the relative change in $\max|\delta\@vec{n}|$ reaches the desired tolerance of $10^{-4}$. At each sub-iteration, only the radiative bound-bound rates $R_{i,j}$ are kept fixed, the remaining rates and the ionization fractions, which depend on $n_{\text{e}}$, are updated accordingly. Algorithm 1 outlines the entire scheme.