Numerical Simulations and Observations of Mg ii in the Solar Chromosphere

The solar chromosphere is the region where the energetics are controlled by a nonthermal "mechanical" heating that becomes dominant in setting the density, magnetic field, and temperature structure. In the lower chromosphere, especially in the quiet Sun, this mechanical heating, driven by granular convective motions, consists of primarily acoustic processes (Carlsson & Stein 1997), and 3D MHD simulations are capable of reproducing the most salient aspects of the spectral lines formed there (e.g., Leenaarts et al. 2009, and references cited within).

However, even while these advanced models reproduce many features of the chromosphere, they also show significant discrepancies for lines whose cores form above the β = 1 (β ≡ p_g /p_B , where p_g is the gas pressure and p_B the magnetic pressure B²/2μ₀) layer, some 750 km above the photosphere in the quiet Sun. Clearly, the magnetic field plays a vital role at greater heights, and this role is not fully understood. This issue has become particularly pressing when interpreting the vast amounts of data collected by the IRIS satellite (De Pontieu et al. 2014), especially in the Mg ii lines.

Mg ii lines are uniquely sensitive to the middle and upper chromospheric conditions and cannot be obtained from ground-based observatories. In radiative MHD Bifrost simulations, synthetic Mg ii lines often came out too faint, strongly asymmetric, or too narrow (Carlsson et al. 2019). These discrepancies could have many causes, such as indicating a lack of opacity, a lack of heating, a lack of small-scale motions in the chromosphere, and/or 3D effects on the radiative transfer (Judge et al. 2020).

The lack of spatial resolution in radiative MHD simulations will cause a reduction in the production of small-scale motions or turbulence and hence result in narrower synthetic line profiles than those observed. As an example, consider the Ca ii 854.2 nm profile discussed in Leenaarts et al. (2009), which was computed with a horizontal resolution of 64 km. This original calculation was redone with a resolution of 48 km, which also produced a profile significantly narrower than that observed. However, when the simulation was repeated with a spatial resolution of 31 km, the correspondence between the line core widths was found to be much better. This indicates both that vigorous large-scale dynamics and turbulent motions at smaller scales could be reproduced in higher-resolution models and that for this mid-chromospheric line it is the turbulent motions that determine the core width. This can be seen in Figure 1 (Mats Carlsson, private communication), which shows the Ca ii 854.2 nm profile at 48 and 31 km horizontal resolution, along with the observed profile taken from the FTS Solar Atlas (Brault & Neckel 1999). This also sets an upper limit on the amplitude of turbulent motions in the lower to mid-chromosphere and implies that a numerical resolution of, say, 20–30 km is good enough to capture relevant solar dynamics in this region.

Zoom In Zoom Out Reset image size

Effects that go beyond the standard MHD description of the chromospheric plasma, e.g., generalized Ohm's law (GOL; e.g., Khomenko et al. 2018), including ambipolar diffusion and nonequilibrium (NEQ) ionization of H and He will also alter the opacity of the Mg ii lines, by altering the heating and electron density profiles of the chromosphere and thus potentially the intensity and width of the emergent lines (Martínez-Sykora et al. 2017; Bose et al. 2021; Przybylski et al. 2022). There are other smaller-scale candidate physical processes for heating the chromosphere as well: the thermal Farley–Buneman instability (TFBI; Oppenheim et al. 2020). Preliminary multifluid simulations of the TFBI show a temperature increase due to heating (Evans et al. 2022). TFBI also provides a large increase of the turbulent motions during the nonlinear phase of the instability, which could broaden chromospheric spectral lines.

Furthermore, the lack of chromospheric heating and dynamics in 3D MHD models could be caused by not properly capturing the magnetic topology of the outer solar atmosphere correctly. The magnetic field has been an almost free parameter of 3D models, and at small spatial scales it is not well constrained by observations. This problem has been partially mitigated by recent models of photospheric, near-surface simulations (Rempel 2014, 2018) showing the importance and depth dependence of the local dynamo in generating quiet-Sun-like fields. The topology and evolution of the magnetic field impact the structure of the chromosphere both through the heating rate and through the Lorentz force, which can carry and support material in excess of hydrostatic equilibrium, changing the density, temperature structure, and opacity, especially in the β < 1 region of the mid- to upper chromosphere (as well as in the corona). Flux emergence is one agent that will change all of these parameters, as it can bring plasma up to chromospheric heights, where it can reside long after the emergence phase is over. However, it is still unclear how much small-scale flux emergence occurs regularly outside of newly forming active regions. This is true for both mature active regions and plage, where the rates are not known. Some progress has been made for typical quiet-Sun regions, where Gošić et al. (2022) measure 68 Mx cm⁻² day⁻¹, perhaps as a result of a local dynamo.

The problem of correctly reproducing Mg ii intensities and core widths is particularly pressing when chromospheric plage is considered. A series of "semiempirical" models (Carlsson et al. 2015) was constructed in order to look into what effect varying chromospheric parameters had on the Mg ii profile in plage-like conditions. The goal of this was to find how a chromospheric atmosphere can produce "single-peaked" profiles, where the k₂ peaks and the k₃ minimum share nearly the same intensity and where the line core width is as wide as observed. The paper gives a number of parameters that conspire to change the Mg ii profile: the electron density n _e , the temperature T _g , and the turbulent velocity v_turb. Furthermore, Carlsson et al. (2015) found that the total intensity of the Mg ii profile is dependent on the temperature, or equivalently pressure, of the overlying corona, which sets the column mass at which the transition region temperature rise occurs and hence the density of the upper chromosphere.

Inversions of observed IRIS spectra tell essentially the same story: in order to reproduce the observed Mg ii profiles, the (plage) chromosphere needs to be dense, extended, and hot, perhaps also with a large (5–6 km s⁻¹) turbulent velocity in the region of line core formation (da Silva Santos et al. 2020; Bose et al. 2022; Sainz Dalda & De Pontieu 2022).

In this paper we continue these studies by comparing synthetic profiles of Mg ii for a set of numerical simulations of varying spatial resolution, magnetic topology, and physics to IRIS observations of quiet Sun, active region, and plage.

We have chosen four large IRIS raster scans to serve as a basis of comparison between simulated and observed Mg ii spectra. Additionally, we present HMI magnetograms (Scherrer et al. 2012; Schou et al. 2012) of the same regions to present maps of the magnetic environment. The sizes of the regions chosen are roughly the same as the simulation box presented later in this paper; hence, they should in principle host magnetic field topologies of roughly the same scale. The raster scans cover two typical quiet-Sun regions, one long-lived apparently unipolar plage region, and an active region sunspot with surrounding plage. We find a range of average Mg ii profile intensities, shapes, and widths in subregions covering of order 10'' × 10''.

The top row of Figure 2 shows a "typical" quiet-Sun location. When considering quiet-Sun profiles in the following, we note the difference, but we largely refer to both network (NW) and internetwork (IN) profiles as quiet Sun (QS). When averaging, we implicitly assume that the spatial filling factor of NW versus IN is the same in the observations and the simulations. To the right of the co-temporal HMI magnetograms, we show the Mg ii k line spectra taken at Sun center on 2015 December 20. The magnetogram shows weak magnetic fields across the entire image, with small patches of stronger (>100 G) fields of both polarities and a diagonal band of strong negative-polarity field stretching from solar (x, y) ≈ (−40'', 30'') to (x, y) ≈ (20'', −30''). The average Mg ii profile is double peaked with a fairly deep k₃ core—half the strength of the k₂ peaks for the darker examples, somewhat less deep for the network ("IRIS NW") box in red. The network patch also shows an average profile that has a slight asymmetry, with the violet k₂ peak some 30% stronger than the red peak. The k₂/k₃ intensities are all of order 1 nW m⁻² Hz⁻¹ sr⁻¹, or less, for internetwork regions with FWHM widths measured to 0.051 nm (corresponding to roughly 50 km s⁻¹), while the network emission is stronger and wider, with intensities of 1.7 nW m⁻² Hz⁻¹ sr⁻¹ and width of 0.056 nm. The second quiet-Sun raster scan was taken on 2017 October 15 and previously analyzed by Martinez-Sykora et al. (2022). We refer to that paper for further details. As shown later, the peak separation and intensities of this observation resemble the quiet-Sun regions (though not the network) of the third panel on the top row of Figure 2.

Zoom In Zoom Out Reset image size

In the middle row the central portions of NOAA Active Region 12296, on 2015 March 8, are shown to consist mainly of negative-polarity plage, with some positive-polarity plage near (x, y) (−25'', −75'') stretching toward the northeast. The positive-polarity plage is actually mostly remnants of NOAA AR 12192, one of the largest sunspot groups of Cycle 24, which crossed the central meridian five rotations earlier, on 2014 October 23. In 2015 March, the negative-polarity plage has dispersed and stretches from the equator to −500'' S. The small bipole forming NOAA AR 12296 emerged into this band of negative polarity near the equator at least a week before the crossing of the central meridian and appears to be fading, but a new active region, NOAA AR 12298, forms in the same location 4 days later (2015 March 12), so there may be weak flux emergence occurring continuously near this location. The quiet-Sun (or canopy) box chosen lies between negative and positive polarities, and the average Mg ii profile in this box is very similar to the quiet-Sun profiles shown in the top row, although it is a factor of two brighter. We note that what we call QS/canopy may be significantly affected by the neighboring active region plage and associated canopy, which is why we refer to these profiles as "QS/canopy" rather than just "QS." The blue box, outlining the total extent of the active region, has a higher intensity with the k₂/k₃ peaks at intensities of 2 nW m⁻² Hz⁻¹ sr⁻¹ and with width 0.059 nm. The red box outlines an area of strong plage emission with maximum k₂/k₃ intensities of more than 4 nW m⁻² Hz⁻¹ sr⁻¹ and FWHM width of 0.059 nm, similar to the entire active region.

The bottom row shows data taken on 2016 January 14, just after 23:00 UT, focusing on NOAA AR 12480. This active region appeared on the east limb on 2016 January 6 and was fully emerged at that point. The IRIS raster covers both a positive-polarity plage region to the east of a large sunspot and, between these polarities, less bright quiet-Sun/chromospheric canopy structures. The average plage profile (red) is very similar to that seen in the NOAA AR 12296 raster, shown in the middle row, with intensities of more than 5 nW m⁻² Hz⁻¹ sr⁻¹ and width of 0.059 nm. Likewise, the quiet-Sun spectra are similar to the quiet-Sun cases discussed above—here, in blue, with intensities of 0.8 nW m⁻² Hz⁻¹ sr⁻¹ and FWHM width of 0.048 nm. The yellow box covering the penumbra shows intensities of 1.5 nW m⁻² Hz⁻¹ sr⁻¹ and width 0.041 nm. Finally, the umbral spectra in green are among the narrowest found, with FWHM width of 0.025 nm and relatively low intensity 0.85 nW m⁻² Hz⁻¹ sr⁻¹. Note that the umbra is on average single peaked with no clear k₃ minimum.

We have run several simulations using the Bifrost MHD code (Gudiksen et al. 2011) to model the photosphere and outer solar atmosphere. In this work we consider models with different resolutions, physics, and field topology. The extension of physics includes NEQ hydrogen and/or helium ionization (Leenaarts et al. 2007; Golding et al. 2016) and/or GOL (Nóbrega-Siverio et al. 2020).

In its base configuration Bifrost solves the equations of MHD using an energy equation that includes optically thick radiative transfer (Hayek et al. 2010), a tabulated form of effectively and optically thin radiative losses (Carlsson & Leenaarts 2012), and Spitzer thermal conductivity along magnetic field lines. The latter is solved either implicitly via operator splitting and utilizing a multigrid solver or explicitly along with the MHD equations using a hyperbolic formulation, which limits the speed of conduction fronts, allowing a reasonable time step, as described by Rempel (2017). The equation of state is, for the "simple" models, based on table lookups of the LTE ionization state of a plasma of solar abundance.

In short, three of the models presented here are run with relatively coarse (100 km) horizontal resolution and include only "simple" MHD physics extending over a 72 × 72 × 60 Mm box. The convection zone is modeled to a depth of 8.5 Mm, and the models reach more than 50 Mm above the photosphere. The size of this computational box allows the capture of granular to supergranular (or at least mesogranular) size scales.

Other horizontally extended domains are the three simulations already described in detail in Martínez-Sykora et al. (2017, 2020). These simulations are 2.5D and span a 90 × 42.8 Mm box. The model contains two plage-like regions connected with ∼40 Mm long loops and reveals features resembling type I and II spicules (the latter only with the presence of ambipolar diffusion), low-lying loops, and other physical processes. The simulation spans a vertical domain stretching from ∼3 Mm below the photosphere to 40 Mm above into the corona with a nonuniform vertical grid size of 12 km in the photosphere and chromosphere and 14 km grid size in the horizontal axis. The three models differ in the included physics, one with only the "simple" Bifrost configuration (called spicule nGOL in the following), one with GOL (spicule GOL), and the third with NEQ ionization of hydrogen and helium and GOL (spicule GOL NEQ(H,He)).

The other two models considered cover small domains (6 × 6 × 10 Mm) but with very high resolution (5 km horizontally) of an internetwork field as detailed in Martínez-Sykora et al. (2019, 2022). The two simulations differ in the physics included, one without GOL or NEQ ionization (called QS nGOL in the following), the other with NEQ ionization of hydrogen and GOL (QS GOL, NEQ(H)). The simulation spans a vertical extent stretching from ∼2.5 Mm below the photosphere to 8 Mm above. A nonuniform vertical grid is employed with a size of 4 km in the photosphere and chromosphere and somewhat larger outside these regions. Initially, the simulation box is seeded with a uniform weak vertical magnetic field of 2.5 G. From this starting point a local dynamo is active and generated a magnetic field that reaches a statistically steady state with B_rms = 57 G at photospheric heights (similar to that described by Vögler & Schüssler 2007; Rempel 2014; Cameron & Schüssler 2015). These models also generated an in situ magnetic field in the chromospheric portion of their domain.

In all models, the upper boundary is based on a characteristic extrapolation of the variables, which in principle allows waves to exit the computational box without reflection (see Gudiksen et al. 2011, for details). The temperature gradient, and hence conductive flux in the vertical direction, is set to zero so that no heat enters the box from above. At the lower boundary, the entropy of the material flowing into the computational box at the bottom boundary is set so that the effective temperature is close to solar.

In order to produce synthetic diagnostic profiles of the optically thick Mg ii lines (as well as the lines of Ca ii), we employ the RH1.5D code (Uitenbroek 2001; Pereira & Uitenbroek 2015). This code performs multilevel non-LTE calculations with partial frequency redistribution. RH calculates spectra from 3D atmospheric models on a column-by-column basis. Note that the 1.5D nature of the solution process will miss some of the effects of horizontal transfer (see Leenaarts et al. 2013; Sukhorukov & Leenaarts 2017, where 3D effects on Mg ii are discussed), which we do not believe are vital for the analysis performed in this paper, though see also Judge et al. (2020) for a slightly alternate view.

3.1. Quiet-Sun Model(s)

Let us begin our analysis by considering low-resolution, 3D MHD, "quiet-Sun" models (with the relationship between the quiet Sun and these models yet to be clarified). These large-scale models were initialized with a thermodynamically relaxed convection zone, photosphere, and lower chromosphere. Two different models were derived from this starting point: One, called nw072100 henceforth, has an average horizontal field of 100 G in the convection zone and photosphere, as well as a very weak horizontal field at coronal heights. In the second, called qs072100 in the following, a vertical field of 5 G is included in addition to the horizontal field used in nw072100. In nw072100, a strong horizontal flux sheet is imposed at the lower boundary but has not yet reached the photosphere for the profiles discussed in this section. Both models were allowed to evolve for several hours solar time, and a salt-and-pepper magnetic field is rapidly established in the photosphere with a mean unsigned B_z of 30 G (mean B of 60 G), which is smaller than the observed quiet-Sun field of 〈∣B_z ∣〉 = 60 G (e.g., Orozco Suárez & Bellot Rubio 2012) by roughly a factor of two.

Though the photospheric field strengths in these models are very similar, we find quite different coronal temperature structures as shown in Figure 3, presumably due to the differing magnetic field topologies: The nw072100 model, which initially had a nearly horizontal field, originally achieved very high temperatures as portions of the initial 100 G photospheric field expand into the corona. This hot corona took several hours to cool, since thermal conduction along the magnetic field could not form an effective cooling mechanism when the connection to the transition region and chromosphere is tenuous. However, the model cools eventually and finally reaches a minimum temperature state after some 3 hr solar time (a vertical slice of the temperature structure of the model is shown in Figure 3). Spatially highly variable, the average temperature lies between 200 and 300 kK from the top of the chromosphere, at 2 Mm, up to 15 Mm above the photosphere. The average temperature decreases to some 100 kK above this height. After this time, the average coronal temperature rises slowly, before rising rapidly when new flux emerges from the photosphere as described in Section 3.4.

Zoom In Zoom Out Reset image size

The qs072100 model was also initially hot, but with a significant amount of vertical field, it cools much more efficiently. The vertical field component also allows the spread of high temperatures more easily from localized braiding-caused heating events. Thus, the temperature structure appears much smoother than that found in nw072100, as shown in Figure 3. Of particular interest to this paper are the chromospheres of both models, both of which are of order 2 Mm in vertical extent, but with elements reaching up to 5 Mm above the photosphere, as seen in the figure.

Figure 4 shows the vertical photospheric magnetic field B_z , the Mg ii line core intensity, and the Mg ii line profiles. As mentioned above, the photospheric field forms a salt-and-pepper-like pattern, organized on a larger scale into what appear to be network cells with diameters of order 10–20 Mm. The fields are concentrated into small patches where we find field strengths up to ±2000 G, while, as mentioned above, the average vertical field in both models is of order 〈∣B_z ∣〉 = 30 G at the time the figure represents. The distribution of the magnetic field is similar in both models, though we do see one clear instance of flux emergence in the nw072100 model near (x, y) ∼ ([10–20], [20–30]). The Mg ii k₃ line core image reflects this network pattern to a certain extent: intensities are high in the same locations in which the field is strong.

Zoom In Zoom Out Reset image size

We also see some remnants of the horizontal field that initially filled the corona as weak horizontal stripes aligned with the y-axis where cool gas remains trapped in the corona. The amount of cool gas at great heights may be greater than what one finds in the typical quiet Sun outside of prominences or regions of strong flux emergence; while interesting, the effect this has on the Mg ii profile appears to be small, aside from a small broadening that disappears in the hours following this snapshot before flux emerges and the corona reheats.

We have chosen three boxes with dimensions of order 10 × 10 Mm in which we compute the average profiles from both quiet-Sun models. ⁵ The observed quiet-Sun k₂ and k₃ intensities have radiation temperatures of roughly T_rad = 5000 K, with the network brighter at T_rad = 5500 K. The model shown in the top row, nw072100, displays radiation temperatures as low as 4900 K and up to 5500 K for internetwork and network intensities, respectively. On the other hand, the model shown in the bottom row, qs072100, is somewhat hotter, with k₂ and k₃ intensities equivalent to radiation temperatures between 5300 and 5600 K. The qs072100 coronal temperature has stabilized, but we note that the coronal density at 10 Mm is (still) decreasing at the time of this snapshot, having fallen from 4 × 10⁹ cm⁻³ to 2.5 × 10⁹ cm⁻³ in the half hour preceding the displayed snapshot. The difference between these models' average chromospheres (below 2 Mm) is slight; both have an average Joule heating of 0.3 W m⁻³ at z = 1.0 Mm above the photosphere, but the qs072100 model, shown in the bottom row, has a slightly hotter corona with an average temperature of 600 kK at z = 10 Mm, while the nw072100 model, shown in the top row, has an average of only 250 kK at z = 10 Mm at this time. We will discuss how coronal temperatures can impact the Mg ii intensities further in Section 4, but note that the coronal temperature plays an important role in setting the Mg ii intensity level. This point is also made by Carlsson et al. (2015) and recently by Bose et al. (2022), who note the correspondence between bright Mg ii core emission and sites of high intensity Fe ix 17.1 nm in moss regions.

While the synthetic quiet-Sun Mg ii intensities lie well within the range of observed values, some interesting differences are also clear. Profiles in the cooling qs072100 model show stronger k_2r peaks than k_2v, whereas the nw072100 profiles, while still showing some asymmetry, are more balanced. More importantly, the Mg ii core widths are seen to be much narrower than what is observed. In the regions covered in both models of Figure 4 the line widths vary between 0.025 and 0.03 nm FWHM which is roughly half of the observed quiet Sun and network profile widths.

Let us now compare these models with a model where the resolution is greater (∼5 km vs. ∼100 km) and with increased physical complexity: the QS nGOL and QS GOL NEQ(H) models. These models were initially run for 51 minutes, at which point the magnetic field statistically reached a steady state with 〈∣B_z ∣〉 = 17 G and B_rms = 57 G, and with some flux concentrations reaching 2 kG at photospheric heights. These models have a simplified magnetic topology that mostly includes granular-scale magnetic fields (typical of internetwork regions) and do not include regions that resemble the quiet-Sun network. In both models magnetic field is generated in situ at chromospheric heights from conversion of kinetic energy to magnetic energy (Martínez-Sykora et al. 2019). As detailed in Martínez-Sykora et al. (2022), in the model that includes GOL and NEQ hydrogen ionization, the magnetic energy built in situ in the chromosphere produces high rates of heating that occur in regions where ions and neutral particles slip in relation to each other, i.e., the upper chromosphere, while NEQ hydrogen ionization has the consequence that in colder regions the electron density remains much higher than when treating hydrogen ionization in LTE. Figure 5 shows the chromospheric temperature structure of the QS GOL and QS GOL NEQ(H) models, including the locations of the height of τ_ν = 1 layers for the k₃ minimum and k₂ peaks.

Zoom In Zoom Out Reset image size

In Figure 6 we show the average Mg ii profiles obtained from these high-resolution quiet-Sun models and compare them with the quiet-Sun models discussed above and with IRIS quiet-Sun and network observations. The low-resolution profile is from a snapshot later in time in nw072100 where the line asymmetries are reversed and correspond better to what is usually observed. We find that the high-resolution models (QS nGOL and QS GOL NEQ(H)) have intensities somewhat higher but of the same order as those observed. However, the profiles are asymmetric, with k_2v substantially brighter than k_2r in both models, perhaps indicating rising coronal temperatures and densities with associated flows. The core profiles are broader than the previous low-resolution, non-GOL, LTE ionization nw072100 and qs072100 models, and though they are indeed a closer match to observed quiet-Sun profiles, we still find that the profiles are narrower than what is observed. The introduction of higher spatial resolution seems to have halved the discrepancy (with the observations) in line width, with QS GOL NEQ(H) giving a marginally broader profile. The synthetic intensity in this model gives a radiation temperature of T_rad ≈ 5000 K and a line width of 0.040 nm FWHM, as compared to 0.051 nm for IN and 0.056 for NW observations presented above. Finally, we note that while the results of Judge et al. (2020) and Sukhorukov & Leenaarts (2017) show a slight increase in the Mg ii h and k line widths owing to 3D radiative transfer, it is clear that the other effects considered produce more drastic variations.

Zoom In Zoom Out Reset image size

3.2. Spicule Model

In the upper chromosphere, plasma β <<1 and the energetics and dynamics of the chromosphere become increasingly dominated by the magnetic field with increasing height above the photosphere. This is particularly true for regions of stronger field such as those characteristic of spicules. To investigate the effect of stronger, more dynamic fields, capable of exciting spicules, on the Mg ii profile emission, the 2.5D atmospheres described by Martínez-Sykora et al. (2017, 2020) are utilized.

This simulated atmosphere spans from the upper convection zone 2.8 Mm below the photosphere to the corona 40 Mm above the photosphere. The model also has a large horizontal extent, spanning 96 Mm. A uniform spatial resolution of dx = 14 km along the horizontal axis and a nonuniform resolution of up to 12 km in the vertical direction ensure that the relevant physics is well resolved. Thus, though the model is 2.5D, convective motions in and below the photosphere are large enough to cause braiding of the magnetic field and a self-consistently heated chromosphere and corona through Joule heating, where we find coronal temperatures up to 2 MK. At the resolution used, ambipolar diffusion is larger than the artificial numerical diffusion by 3−5 orders of magnitude in extended regions of the chromosphere. The magnetic configuration allows the formation of structures that closely resemble spicules of type ii, especially in the simulations that include GOL, as discussed in Martínez-Sykora et al. (2017).

The initial magnetic field in this model contains two plage-like regions of high magnetic field strength and opposite polarity. These are connected by magnetic loops that are up to 50 Mm long, some footpoints of which are shown in Figure 7, which also shows the temperature structure of the chromosphere and transition region for all three models. The mean unsigned value of the vertical magnetic field 〈∣B_z ∣〉 = 86 G, B_rms = 271 G, with magnetic flux concentrations reaching up to 1400 G.

Zoom In Zoom Out Reset image size

Figure 8 shows the Mg ii profiles resulting from this model when run with "simple" MHD, with GOL, and with GOL and NEQ hydrogen and helium ionization (see also Bose et al. 2021). The profiles are computed in a ∼3000 km wide region containing the strongest magnetic fields, which is the site of strong spicule formation. The k₂ and k₃ intensities are stronger than that found in the quiet Sun and, though slightly weaker, are approaching the intensities found in active regions. This is especially true for the model where both GOL and NEQ ionization are considered. All three profiles are strongly asymmetrical, more so than the observed active region profiles, with the k_2v peak significantly brighter than k_2r. The synthetic line widths also approach or surpass those measured in the observed quiet Sun, with widths of 0.04 nm FWHM in the non-GOL and GOL models and 0.05 nm FWHM in the GOL and NEQ hydrogen and helium ionization model. However, these widths are still slightly narrower than those found in typical active regions and plage as shown with dashed lines in Figure 8 (see also Figure 2).

Zoom In Zoom Out Reset image size

3.3. Plage Model

When strengthening the magnetic field and expanding the areas covered by unipolar field, the dynamics of near surface convection, the structure of the chromosphere and corona, and the heating processes that determine the thermal and dynamic structure of the outer atmosphere are altered. A model typical of plage or small active region, named pl072100 in the following, illustrates these issues. This model is initiated with two opposite-polarity patches of vertical magnetic field covering a large percentage of the modeled photosphere. The model has the same dimensions and is run with the same resolution as the large-scale quiet-Sun simulations, nw072100 and qs072100 discussed above. In contrast to the quiet-Sun models, the field in pl072100 is quite strong, with an average unsigned vertical magnetic field strength of 〈∣B_z ∣〉 = 180 G. The field concentrations of the "plage" regions are of order ±1800 G, but slightly weaker in the positive-polarity plage. Though the field was originally confined, convective motions have caused the diffusion of quite strong fields also into the more quiet areas of the simulation at the time of the featured snapshot.

The left panel of Figure 9 shows the geometry of the photospheric magnetic field. The Mg ii k₃ intensity is featured in the second panel of this figure and shows bright speckled areas in locations overlying the strong magnetic plage that appear quite similar to those observed (see Figure 2). The third panel contains average Mg ii core intensities over four boxes covering regions of interest (ROIs) with different magnetic field configurations, placed both in regions of the strongest magnetic field and in areas between and away from the two main polarities. We find very bright emission in ROIs located in the strongest plage, with intensities >20 nW m⁻² Hz⁻¹ sr⁻¹, which is much greater than the typical measured plage intensities of 5 nW m⁻² Hz⁻¹ sr⁻¹.

Zoom In Zoom Out Reset image size

The typical observed plage intensity translates to a radiation brightness temperature of T_rad ≈ 6300 K, while ranging from average radiation temperatures of T_rad = 6100 K for the k₃ core in the "QS/Canopy" patch, through T_rad = 7000 K for the (blue box) region between the two polarities, to T_rad = 7700 K for the brightest, plage-like patch. The last profile is single peaked, and the small-scale structure of the emission is similar to what is observed in observed plage. However, when we consider the line widths, we again find that they are all much narrower than what is observed, of order 0.025 nm FWHM, which is less than half what is observed. Instead, the synthetic widths are very close to what is found in sunspot umbrae as can be seen by comparison with Figure 2 (see also Bose et al. 2022). Comparisons should also be made with the profiles found in the models seen in Section 3.2, which have strong fields, but at higher resolution and with other physical processes at work. Likewise, we find strong-field regions and single-peaked profiles in models that include flux emergence as described in the following section.

3.4. Emerging Field Models

The two low-resolution quiet-Sun models presented above had line core widths significantly narrower than what is observed, but also average unsigned magnetic field strengths of only 〈∣B_z ∣〉 ≈ 30 G. The line widths were greater in the quiet-Sun model that was run at high resolution and with GOL and NEQ hydrogen ionization and a stronger, local dynamo generated, field, but even at this resolution (5 km) they did not reproduce the observed widths. Let us now reconsider the low-resolution case, but in which the field strength gradually rises, e.g., in an emerging flux model.

The emerging flux model is a continuation of the nw072100 quiet-Sun model described above. Flux emergence is initiated by injecting into the original configuration a flux sheet, aligned with the y-axis, of strength B_y = 200 G at the entire bottom boundary for 95 minutes. After this initial period, the flux sheet strength was increased to first B_y = 1000 G for 70 minutes and thereafter 2000 G for another 150 minutes. Afterward, the strength of the injected field was reduced to B_y = 300 G, which is injected continuously thereafter. Figure 10 shows the evolution of the field in the middle of the computational domain (y = 35 Mm) over an 8 hr period. The initial 100 G field is kneaded and pulled down by convection below the photosphere in several locations, while some field rises through the chromosphere and into the corona in the first few hours of the simulation. In the corona the field is nearly horizontal and the coronal plasma becomes quite hot, only cooling slowly as the field diffuses through radiative cooling and losses through conduction. The latter is hampered by the weak thermal connection between the corona and chromosphere.

Zoom In Zoom Out Reset image size

This early coronal field has a strength of some 20 G, decreasing from the fourth to the fifth hour to a minimum of 10 G. From the fifth hour onward, portions of the strong injected flux sheet reach the photosphere, rising to the upper atmosphere layers above, in places where the emerging field is strong enough in the photosphere to overcome its lack of buoyancy (Acheson 1979; Archontis et al. 2004).

The expanding and increasing coronal field leads to a restructuring of the coronal field topology, the introduction of cool plasma carried along with the field, and rising coronal heating rates. The temperature structure in the model toward the end of the run (at 8 hr 21 minutes) is shown in Figure 11. As in the quiet-Sun models, we see that chromospheric temperatures are largely confined to the 2 Mm above the photosphere. However, there are several regions where cool, T_g < 30,000 K gas is present up to 5 or even 10 Mm above the photosphere. This cool gas has been carried up into the corona by the magnetic field as the field expands from the photosphere, and now forms cool fibrils. Furthermore, in regions where heating is strong, we find loop-shaped structures with temperatures T_g > 5 MK where the chromosphere is compressed and the transition region temperature rise occurs already at 1.5 Mm above the photosphere, very similar to the structure we find in the plage model pl072100.

Zoom In Zoom Out Reset image size

In Figure 12 the average field, coronal Joule heating, and Ca ii K, Mg ii k, and Fe ix 17.1 nm line intensities are shown as a function of time, starting roughly an hour before the emerging field reaches and pierces the photosphere and for the next 5 hr. Five hours into the simulation run, emerging flux that has collected just below the photosphere breaks through and interacts with the ambient coronal field; this causes a spike in (coronal) Joule heating, which increases by more than a factor of three for a short interval, to 2 × 10⁻² W m⁻³, with an accompanying spike in both the Mg ii and Fe ix intensities. There is no corresponding spike in the Ca ii intensity, and we note that the chromospheric Joule heating does not increase at this time either, but rather remains at Q_J = 0.3 W m⁻³. After the spike, the Fe ix emission remains high, while the Mg ii intensities fall back to their previous value as soon as the Joule heating subsides. The mean vertical magnetic field strength 〈∣B_z ∣〉 in the photosphere remains at some 30 G until 5.5 hr, at which point it rises with the increasing emergence of magnetic flux, passing 60 G at 6.75 hr, and reaching 100 G at the 8 hr mark, at which point the mean strength remains nearly constant.

Zoom In Zoom Out Reset image size

Figure 12 attests that as the average magnetic field strength rises, so does the Joule heating rate, both in the chromosphere at z = 1.0 Mm and in the corona at z = 10.0 Mm (and in points in between). The chromospheric heating remains at Q_J = 0.3 W m⁻³ until hour 6, thereafter rising rapidly to >3.0 W m⁻³ in the span of 2 hr. Except for the spike in heating around hour 5, mentioned above, the coronal heating rate mirrors the chromospheric rate, but at a level that is a factor of 100 lower in amplitude.

The increased heating rate impacts both the intensities of Ca ii K and Mg ii k in the chromosphere and the Fe ix 171.1 nm line in the corona. We compute the average line core intensities of Ca ii and Mg ii integrated over a 0.02 nm FWHM wide Gaussian filter and the total line intensity of the Fe ix line. The increase in the Mg ii line intensity becomes especially apparent when the average field becomes larger than 〈∣B_z ∣〉 ≈ 70 G, when we see a rapid rise in the core intensity.

Let us now consider the Mg ii profile in greater detail. In Figure 13 we find line plots of Mg ii k, as well as a map of the Mg ii k and the Mg ii 2799 triplet line along x = 32 Mm (dashed yellow line in the second panel from the left). These profiles are taken from periods just after the field has pierced the photosphere and at later stages as the field has risen into the outer atmosphere and pervades the steadily more magnetically active chromosphere and corona.

Zoom In Zoom Out Reset image size

The top row shows the Mg ii line at t = 5 hr 23 minutes, some 67 minutes after the previous nw072100 model "quiet-Sun" snapshot shown in Figure 4. The intensities at this time, 1 hr later, are roughly unchanged and of the same order as those obtained with the IRIS observations of quiet-Sun regions. The k₂ peaks are slightly asymmetric, but at this time with the k_2v peak brighter than the k_2r peak, as is often observed, and as opposed to the snapshot shown in Figure 4. However, the line profiles are still too narrow to resemble those observed. We note that the vertical unsigned magnetic field at this point in time is still only slightly larger than 30 G.

The middle row shows the Mg ii spectra at t = 7 hr 3 minutes when the unsigned magnetic field strength has risen to 〈∣B_z ∣〉 ≈ 60 G. At this time we find that the Mg ii core intensities (in two of the three ROIs shown) have increased by roughly a factor of 2—and equivalent to the radiation temperature T_rad increasing by 400 K from 5250 to 5650 K. The line widths have all increased significantly and are approaching the same widths as observed quiet-Sun line widths. The profiles in the "quiet" regions outlined by red and green boxes are twin peaked (or complex), while the profile formed in the region of strongest field is single peaked and with a high, and rising, intensity of 2 nW m⁻² Hz⁻¹ sr⁻¹. We note that, given the amount of flux emergence in the simulation, it is not fully clear which type of solar region (in terms of observations) is the best comparison point.

Finally, in the bottom row of Figure 13 the profiles at simulation time t = 8 hr 21 minutes are presented. At this time the magnetic field has reached 〈∣B_z ∣〉 ≈ 100 G with an accompanying increase in both chromospheric and coronal heating Q_J, as is visible in the middle panel of Figure 12. The two darker ROIs, outlined with red and green boxes, are centered on darker "quiet-Sun," or canopy, areas that both have line cores that are as wide as is observed with intensities that lie between those found in observed average quiet-Sun and active regions, i.e., between 0.5 and 2 nW m⁻² Hz⁻¹ sr⁻¹. The ROI that covers the region of most intense activity, outlined by the blue box, shows a very intense, of order 10 nW m⁻² Hz⁻¹ sr⁻¹, single-peaked average profile that is wide, though not as wide as the "quiet-Sun"/canopy profile derived from the area covered by the red box ROI. These darker profiles, along with the equivalent profiles from t = 7 hr 45 minutes, are shown in greater detail in Figure 14.

Zoom In Zoom Out Reset image size

The formation of Mg ii and, in particular, the details of how the intensity and width of the line core, k₁, k₂, and k₃, are formed are quite complex and dependent on several aspects of the chromosphere, corona, and magnetic field, both along the radiating ray and in the general vicinity of the emitting plasma. These include the chromospheric density and temperature structure, velocity flows, and turbulence, which all can and will alter the opacity and emissivity of Mg ii. The formation of Mg ii is discussed in detail in Leenaarts et al. (2013). Furthermore, the root cause of "opacity broadening" is discussed in depth in Rathore & Carlsson (2015). We here repeat some of the general formation properties most relevant for the context of this paper.

The emergent intensity I_ν of a spectral line can be approximated by the Eddington–Barbier relation

$$\begin{eqnarray}&&{I}_{\nu }(\mu =1)={\int }_{0}^{\infty }{{Se}}^{-{\tau }_{\nu }}{\kappa }_{\nu }{dz}\ \approx S({\tau }_{\nu }=1),\end{eqnarray} \tag{ 1 }$$

where S is the source function, κ_ν is the opacity, dz is a length element along the ray, and τ_ν is the optical depth at frequency ν. In principle, the source function is a function of frequency as well for lines in which partial frequency distribution is important, such as outside the k₂ peaks of the Mg ii line, but for the purposes and context of this discussion, we assume that this can be ignored. Also ignored are the effects of horizontal radiative transfer, which is important for the core of the Mg ii h and k lines (Leenaarts et al. 2013).

Instead, let us concentrate on the frequency-independent source function as a function of height in the outer solar atmosphere. Deep in the atmosphere, the mean free path for (all) photons in the line is very short, essentially no photons escape, and conditions are very close to local thermodynamic equilibrium, such that

$$\begin{eqnarray}&&S=J={B}_{\nu }({T}_{g}),\end{eqnarray} \tag{ 2 }$$

where J is the mean intensity and B_ν (T_g ) is the Planck function, which depends on the local gas temperature T_g . In general, the source function will depend on both the local temperature through B_ν and the radiation field J. As one moves outward in the atmosphere, B_ν will follow the chromospheric temperature structure, while the mean intensity J will decrease as photons can more readily escape the atmosphere. The source function will initially follow B_ν but will eventually fall toward J as the radiation field and local emissivity decouple.

Magnesium is an abundant element, with Mg ii being the dominant ionization state in the chromosphere at temperatures below some 15 kK. The line core (k₃) is therefore formed close to the transition region temperature rise (see Figures 5 and 7 or Leenaarts et al. 2013). At that point, the source function is often sufficiently decoupled from the Planck function such that the value of the source function is higher farther down in the atmosphere. The k₂ peaks are formed at a wavelength separation from the line core such that τ_ν = 1 falls at this local maximum of the source function. For a larger column mass difference between the locations of the ${\tau }_{{k}_{3}}=1$ and ${\tau }_{{k}_{2}}=1$ points, we need to go farther out in the absorption profile to have an opacity low enough, and we get a broader intensity profile.

A high density at the height where the k₃ is formed will ensure that S is close to B_ν (T_g ), and a high temperature in the upper chromosphere will give a high intensity. This will, for example, occur if the coronal temperature is high, >2–3 MK, and the thermal conduction forces the transition region to small geometric heights and greater densities as is the case in the plage model shown in Figure 9 and the later stages of the flux emergence model, where a small coronal bright point with associated high coronal temperatures is forming. On the other hand, a coronal temperature <1 MK leads to a transition region located at large heights, z > 3 Mm above the photosphere, and a source function S that is decoupled from the Planck function and approaching the steadily decreasing mean intensity J, leading to low k₃ intensities.

Note that it is not only the coronal temperature that can modify the transition region geometry and density; a strong horizontal magnetic field with associated Lorentz force will do the same, for example, in the form of spicules, or low-lying loops, some forming the fibrils seen in Hα, the Ca ii H and K lines, or Mg ii h and k. Likewise, emerging flux can raise cold photospheric material to great heights, modifying the chromospheric geometry, while at the same time being the source of enhanced reconnection activity and heating as the emerging field comes into contact with the preexisting ambient field.

The k₂ intensity is usually higher than k₃, as it is formed at the height where S is still strongly coupled to B_ν (T_g ) and thus the chromospheric temperature rise. There are two ways to obtain a single peak profile. If the upper chromosphere density is so high that the source function still remains coupled to the local temperature, there will be no local maximum of the source function and therefore no k₂ peaks and a "single-peaked" profile. An alternative could be for cases (e.g., in the umbra, where faint single-peaked profiles are common) where the density is low throughout the chromosphere so that the source function is never well coupled to the local temperature, so that there is no local peak in the source function with frequency.

The k₁ spectral feature, the point of minimum intensity just outside the core, is found at the frequency where photons can escape readily from the chromospheric temperature minimum, some few hundred kilometers above the photosphere. A schematic cartoon of the relationship between S, J, B_ν (T_g ), and the k₁, k₂, and k₃ frequencies is presented in Figure 15. We note that the temperature structure of the chromosphere can be highly time variable, and indeed, for the lower chromosphere the time variability can be its most salient feature (e.g., Carlsson et al. 2019, and references cited therein).

Zoom In Zoom Out Reset image size

Two processes set the width of the Mg ii k core: broadening of the atomic absorption/emission profile due to small-scale turbulent velocities in the chromosphere at the locations of core emission, and "opacity broadening" as also discussed by Carlsson et al. (2015) and Rathore & Carlsson (2015). The former could be the result of high-frequency waves, motions driven by episodic heating events due to magnetic reconnection, or motions resulting from instabilities such as the Kelvin–Helmholtz instability (Antolin et al. 2015) or the TFBI (Evans et al. 2022). However, the measured line width of the optically thin O i 135.4 nm line in plage regions (Carlsson et al. 2015; Lin & Carlsson 2015) is only of order 10–15 km s⁻¹ FWHM, ⁶ which is not enough to explain the Mg ii k core width. We note, though, that implicit in this argument is the assumption that O i is formed in the same region of the atmosphere as Mg ii. While numerical modeling of the quiet Sun suggests that this is true (Lin & Carlsson 2015), it is not clear whether that is the case for all regions on the Sun. Another way of obtaining a broad spectral core profile is by increasing the vertical extent of the dense, chromospheric temperature plateau. This will form a broad range of frequencies in which the source function is coupled to the local temperature and only far from the line center, ν₀, will S decrease to half of the peak value (where the FWHM intensity is formed).

Thus, in order to reproduce the observed intensities and line core widths of the Mg ii k line, a solar chromosphere that is hot and dense over an extended range in height, or more accurately in optical depth, is required. In quiet-Sun models, we do find an increase in core width in higher-resolution models, as shown in Figure 4, which could be caused by increased turbulent velocities and/or increased mass loading of the upper chromosphere because of more concentrated energy release and resulting stronger motions. It is also possible that the topological differences between the various QS models play a role in the different line widths. A further increase of the width is found when considering models where GOL, and especially GOL and NEQ hydrogen (and helium) ionization, are included. The GOL models feature higher average chromospheric temperatures, larger chromospheric scale heights, and hence greater mid- and upper chromospheric densities. However, while models including GOL and NEQ ionization do increase the width of the line core, reducing the discrepancy significantly, neither of these models sufficiently changes chromospheric structure enough to match observed quiet-Sun Mg ii core widths.

Models of the more active Sun, in which the magnetic field plays a more prominent role, are also capable of producing larger line core widths (Figure 8), especially when GOL and NEQ hydrogen and helium ionization are included. Spicule dynamics bring significant mass up into the upper chromosphere in these simulations. This increases the density, while currents associated with the spicule acceleration heat the upper chromosphere, transition region, and lower corona, leading to the atmospheric structure needed to produce large line widths. Yet, while the synthetic line core widths from these models are of the same order as what is observed in the quiet Sun, they are not quite wide enough to reproduce the line profiles seen in active regions or plage. We note, though, that these models are 2D models, and it remains unclear whether expansion into 3D would further reduce the discrepancies.

The importance of the magnetic field strength and structure is also apparent in the models featuring flux emergence. In the most active portion of the post-emergence nw072100 model (outlined in blue in Figures 13 and 16) hot coronal loops reach temperatures of >5 MK. This leads to high coronal pressures and a highly compressed chromosphere. This chromosphere, while geometrically foreshortened, is dense and hot in its upper portion (in similarity with the plage model pl072100). The heating driving the high chromospheric and coronal temperatures comes both from high angle reconnection of emerging and preexisting ambient field lines and from the small angle reconnection due to the braiding of already present field lines (e.g., as in Bose et al. 2022). Note that while Mg ii k line intensities are quite high, higher than those observed in plage, and while the profiles are single peaked, their average width is significantly smaller than that observed.

Zoom In Zoom Out Reset image size

We find that the column mass, m_c , can be a good proxy for optical depth and that ${\mathrm{log}}_{10}({m}_{c})=-2.2$ (kg m^–2) has a rough correspondence with the height of k₃ emission. An overview of the magnetic field at coronal heights, the column mass height, the Joule heating per unit mass in the chromosphere (1.5 Mm above the photosphere), and the emission in the Fe ix 17.1 nm spectral line showing sites of strong coronal heating and temperatures are displayed in Figure 16. The figure also contains plots of average and individual line profiles, as well as the chromospheric structure of the temperature, vertical velocity, and particle densities as a function of height above the photosphere. The average vertical fields in the larger green, red, and blue ROIs are 〈∣B_z ∣〉 = 9, 37, and 150 G, respectively, and the chromospheric heating is much stronger in the central blue region than in the relatively more quiet green and red areas. However, note the strong shock structures seen in the second (green) row and also the large density scale height in both the red and green rows. The Fe ix image implies that both of these locations are covered by cool canopy-like material. The emission from the central (blue) region is extremely bright and is formed over a compressed dense chromosphere, with the transition region to the corona placed only slightly above 1 Mm above the photosphere.

The role of the global magnetic field is illustrated by considering regions far away from the strongest fields. In Figure 14 profiles from late in the nw072100 model are shown. At this stage of the simulation, profiles that are nearly identical to those found in the quiet Sun and in the average spectra of a small active region are found. While the average vertical field strength 〈∣B_z ∣〉 at this stage of the simulation is high, higher than that measured in the quiet Sun, we find that the average photospheric field in the approximately 10 × 10 Mm regions (outlined in green and red in Figure 13), directly below the chromosphere forming the Mg ii lines, is only 10–20 G: it is the larger-scale field forming longer loops above the ROIs that plays the leading role in forming chromospheric structure by bringing cold material to great heights and thereafter holding it aloft. A key question, of course, is whether such a scenario represents the quiet-Sun environment on the Sun well, even if the spectral line properties are improved. After all, the sequence of events in nw072100 involves prior flux emergence on very large scales, which is key to bringing cold material to great heights, and subsequent weakening and dispersal of the magnetic field, while holding the previously injected cold material at great heights. Such large-scale flux emergence is not the typical cause for quiet-Sun magnetism, which is thought to be sustained through the continual emergence of much smaller scale ephemeral regions and, to a lesser extent, decaying active regions. Perhaps the relative success of nw072100 points instead to the key role of mass loading into the chromosphere, which likely is not properly captured by current simulations.

We have found that key contributors to the observed widths are the chromospheric heating, mass loading, and spatial extent, while the Mg ii intensities are strongly coupled to both chromospheric and coronal heating. While these conclusions were already clear from the semiempirical modeling of plage in Carlsson et al. (2015), the forward models presented here show possible physical conditions and processes required in order to form chromospheres that can reproduce the observations.

We present several quiet-Sun models, listed in Table 1; those with coarse resolution have widths of 0.025 nm, while both the non-GOL and GOL models with high 5 km resolution are substantially broader at 0.04 nm. This is a similar sensitivity to resolution as in the case for the Ca ii 854.2 nm line shown in Figure 1. However, we note that there still is a difference between observed quiet-Sun widths, which are of order 0.05 nm for internetwork regions and 0.056 nm for network regions.

It is first when we introduce more mass at greater heights into the chromosphere than what is typically found in semiempirical models, such as the 2 Mm high VAL3C semiempirical model (Vernazza et al. 1981) or the qs072100 model presented here, that we find widths approaching those observed.

The results of models featuring emerging flux regions or spicules, where mass is carried or thrown into the upper chromosphere and lower corona, show a much improved correspondence with IRIS observations. This is true in terms of both line intensities and line widths. The magnetic field is the key player in these phenomena, underscoring the importance of capturing both the magnetic field strength and topology in simulations. On the other hand, the plage model and the more active parts of the flux emergence simulation still produce Mg ii profiles with widths smaller than those observed, indicating that we are likely missing something important in our understanding of the more active Sun's and unipolar plage's geometry, dynamics, and heating processes.

It seems clear that turbulent motions alone are not likely to be the sole root cause of the observed Mg ii core line widths. This is because the observed discrepancy would require very large values of turbulence that appear to be incompatible with current observations from other lines, e.g., O i 135.5 nm (Carlsson et al. 2015). However, it should be noted that it is not fully clear whether the line formation region of this optically thin line is very similar to that of Mg ii (as suggested by Lin & Carlsson 2015 for quiet-Sun simulations) for all solar regions.

It is interesting to note that the good correspondence found in quieter regions, which have undergone significant mass loading through large-scale flux emergence, is possible even with the relatively coarse resolution of 100 km. This suggests that the heating of injected mass associated with the field, and hence opacity broadening, is a key agent of the large Mg ii widths. The question, then, is which process dominates this mass supply, with both spicules and flux emergence candidate processes considered in this paper. Neither of these processes is fully captured with current models. The spicule models we have shown include GOL and NEQ ionization but are limited to 2D. The current simulations provide a better match with observed spectral line properties but are not sufficient. It remains to be seen whether 3D models can further reduce the discrepancies with observations. Similarly, flux emergence models typically include emergence on very large spatial scales that are a good fraction of an active region size. Some of these models during late stages similarly provide a much better match with observed spectra, with discrepancies much reduced. While such large-scale emergence clearly plays a dominant role in the formation of active regions, it is not clear whether the scenario involved in the large-scale emergence applies to the formation of quiet Sun, for which smaller-scale emergence is thought to play a key role. Results from sufficiently high resolution models, such as the quiet-Sun model discussed in this paper (see also Martinez-Sykora et al. 2022), predict that a local dynamo will be active. When sufficiently deep, these models predict, even in the quiet Sun away from any contribution from the global dynamo, fields of order 〈∣B_z ∣〉 ≈ 60 G (e.g Rempel 2014, 2018). This may be sufficient to cause "continual" flux emergence, but whether the field will be strong enough to significantly perturb chromospheric structure remains to be seen. It is also not clear whether such models include sufficient numbers of medium-scale ephemeral regions, which are known to affect the chromosphere (Gošić et al. 2021) and which are thought to play, through continual emergence, a dominant role in supplying the quiet-Sun network.

It is also clear from our work that very strong magnetic field regions such as plage continue to represent a challenge in terms of reproducing the observed properties of Mg ii. It is unclear whether mass loading alone can resolve this issue. Our results indicate that NEQ ionization and non-MHD effects such as GOL play an important role in reducing discrepancies with observations. However, these effects are computationally expensive and have not been studied for all different types of simulations. In addition, it is still unclear whether these effects can further reduce remaining discrepancies, for example, at higher field strengths, which are known to occur in plage. Another unknown aspect is whether other multifluid or microphysics effects, such as the TFBI, which have not been studied under realistic chromospheric conditions but are thought to potentially play a role in the chromosphere (Oppenheim et al. 2020; Evans et al. 2022), can significantly reduce discrepancies.

Our results provide a path forward for further studies focusing on understanding the formation of Mg ii. The results of higher numerical resolution simulations indicate that a higher resolution appears to significantly reduce discrepancies with observations, suggesting that the numerical approach plays a key role. This likely goes beyond the direct effects of higher velocities on the broadening of the line, but it also includes the significant increase in heating and mass loading as a result of higher resolution. Furthermore, having a large spatial extent, allowing field topologies at many scales, is clearly important. The challenge, then, is to produce numerical simulations that combine many of these different effects in order to determine whether the combination of these various impacts leads to a full explanation of the average Mg ii profiles. Beyond that, it is also clear that much can be learned from studying the spatial distribution and temporal evolution of the simulated profiles, as well as from investigating where in the simulations the discrepancies are not present. That will all be part of future work.

We gratefully acknowledge support by NASA grants 80NSSC20K1272 and contract NNG09FA40C (IRIS). The simulations have been run on the Pleiades cluster through the computing project s1061 and s2278 from the High End Computing (HEC) division of NASA and on the Betzy cluster through computing project nn2834k from the Norwegian Sigma2 High-Performance Computing center. This research is also supported by the Research Council of Norway through its Centres of Excellence scheme, project No. 262622.