Our goal is to investigate how well a simulation computed with the chromospheric extension of the MURaM code (MURaM-ChE) reproduces the spatially averaged  \(\lambda 854.2\) nm line profile. We aim to determine the contributions from the isotopic splitting versus the dynamics in the atmosphere to the resulting line width and asymmetry. For this purpose, we forward model the line based on a simulated enhanced network region. Our study builds on forward modeling of the  \(\lambda 854.2\) nm line in a series of MURaM-ChE simulation snapshots representing an enhanced network region. We solve the radiative transfer problem three times, once considering only the most abundant isotope of calcium in the atmosphere, once taking six calcium isotopes into account, and finally using a single “composite’’ atom model, which mimics the presence of all six isotopes. We find the forward modeled spatially and temporally averaged spectra to be in good agreement with the ftsatlas observation of the quiet Sun. In order to match the observed line width, the simulated atmosphere must be sufficiently dynamic. The typical red-asymmetry can only be reproduced by taking the isotopic splitting effect into account, as suggested in the literature. The closer match between the new model and the observations compared to earlier numerical models is a result of the higher rms-velocity in the MURaM-ChE chromosphere. The center of the spatially averaged line profile tends to be slightly red-shifted, which is a result of a net downflow velocity at the formation height of the line center intensity. This does however not imply average mass downflow. We find the composite atom model is a good approximation to the full isotope computation, but shows some differences in the line core and asymmetry. We show that forward modeling of the  \(\lambda 854.2\) nm line from an MURaM-ChE simulation can result in a close match to the line shape of an average quiet Sun observation. The atmosphere must be sufficiently dynamic to match observed line width. Our results confirm that it is important to include the isotopic splitting effect of calcium in modeling the  \(\lambda 854.2\) nm line.

1 Introduction↩︎

The chromosphere is a component of the solar atmosphere that connects the photosphere to the overlying hot corona. Spectral line formation in the chromosphere is complicated, and the interpretation of the observed fine structure is difficult [1], [2]. Observations of strong spectral lines that form in the chromosphere are the key to inferring thermodynamic quantities such as velocity and temperature, as well as the magnetic field. Interpreting the observations, however, requires a detailed understanding of line formation. The ne and nlte conditions in the chromosphere are reflected in the line formation. Among the available chromospheric diagnostics, the infrared triplet lines have advantageous properties. Their location in the infrared makes them observable from the ground. They are not subject to scattering as much as other chromospheric lines such as H&K and H\(\alpha\) [3] and to ne effects [4], which makes them easier to model and interpret. In addition, the infrared triplet lines, in particular  \(\lambda 854.2\) nm, are good candidates for magnetic field inference in the chromosphere through inversions [5]–[8].

Spatially averaged observations of the qs, for example with the ftsatlas [9], indicate a red-asymmetry in the core of the  \(\lambda 854.2\) nm line [10], which is in contrast to photospheric lines that show a blue-asymmetry [11], [12]. In the photosphere, a different area coverage, line strength, and composition of upflows in the granules and downflows in the intergranular lanes lead to the observed blue asymmetry. The asymmetry of spectral lines is typically inferred from the line bisector. In the case of the  \(\lambda 854.2\) nm line, the bisector has the shape of a mirrored “C’’ letter and is, therefore, often called an”inverse C-shape bisector’’. In contrast, other chromospheric lines such as H\(\alpha\) show a C-shaped bisector [13].

[10] suggested that the inverse C-shape is a result of acoustic waves that steepen into shocks at chromospheric heights, similar to the analysis of K2V bright grains [14]. The wave propagation results in an asymmetry in the time the plasma spends in up- and downflows. By using three-dimensional convective simulations from [11], [10] found a red asymmetry in the spatially averaged  \(\lambda 854.2\) nm profile, but of only roughly one-fifth of the observed amplitude. Numerical models of the solar chromosphere, in combination with three-dimensional rt computations, failed to reproduce the observed asymmetry [15]. The profiles instead showed symmetric line cores.

[16] included all isotopes of calcium that exist stably in the solar chromosphere in the rt calculation. By solving the rt problem simultaneously for six isotopes of calcium, the authors found significant differences compared to the computations where only the most abundant isotope was used. For their analysis, the authors used the FAL C model [17], a time series from a RADYN simulation [18], and a model computed with the Bifrost code [19]. The synthesized  \(\lambda 854.2\) nm profiles from all three model atmospheres resulted in a red asymmetry in the line core, that is an inverse C-shaped bisector.

While there has been progress in explaining the observed  \(\lambda 854.2\) nm line, the forward modeled spectra still show discrepancies. Although the spectra show a red asymmetry, the exact shapes do not match the observations, and the bisector amplitude is too low. For chromospheric lines, the spatially averaged line width of forward-modelled spectra is typically too low [7], [20], [21]. The mismatch of the line width is assumed to be due to a too-low magnitude of dynamical motions [22].

In recent work, [23] showed that new models of the chromosphere are able to match the line properties of the h&k lines relatively closely. The authors find that high maximum velocity differences in the chromosphere may explain the observed line width. The atmosphere model resembles an enhanced network region similar to the setup of the public Bifrost snapshot [22]. In this work, we study the spatially averaged  \(\lambda 854.2\) nm line in the same model of the chromosphere as in , which is simulated with the recently developed chromospheric extension of MURaM muramche¹, [24]. We aim to identify the requirements to match the observed spatially averaged line profile. While the effect of isotopic splitting on the line asymmetry was demonstrated by [16], it is not clear how important this effect is in the muramche model used here. We therefore study the role of the dynamic motions versus the isotopic splitting effect in the resulting spatially averaged line profile.

This paper is structured as follows. In Sect. 2 we present the atmosphere model and rt calculations. In Sect. 3 we present our results. In Sect. 4 we summarize and discuss the results, and in Sect. 5 conclusions are given.

2 Model atmosphere and forward modeling↩︎

We model the  \(\lambda 854.2\) nm spectral line in an atmosphere calculated with muramche, which is a rmhd code. The MURaM code originally included the physics required to simulate near-surface convection [25] but was limited to a lte approximation of the radiative and atomic physics. [26] extended the code to include optically thin losses and heat conduction along magnetic field lines to model the corona.

Recently, [24] developed the chromospheric extension of MURaM, which includes an ne treatment of hydrogen in and above the photosphere. The convection zone is modeled by a non-ideal eos generated with the free-EoS package [27]. Above the photosphere, a ne eos is used, which includes a time-dependent treatment of hydrogen ionisation following the prescriptions in [28]–[30]. In the eos, all non-hydrogen elements are treated in lte. The code includes chromospheric line losses and optically thin losses in the corona, as described in [31]. In addition, the code includes 3D euv back-heating of the chromosphere similar to [31]. Radiative losses in the photosphere and low chromosphere are calculated using a four-band multigroup short-characteristics rt scheme [32], [33], which was extended to include scattering effects as in [34] and [35]. The solar abundances are taken from [36]. The formation and dissociation of \(H_2\) molecules are treated in ne, while \(H_2^{+}\) and \(H^{-}\) are treated in chemical equilibrium. A slope-limited diffusion scheme is included as described in [37] and [26].

The simulation setup resembles an enhanced network region, similar to the public Bifrost snapshot, and is the same as the one used in . The simulation domain spans \(24 ~\mathrm{Mm}\times 24~\mathrm{Mm}\times 24 ~\mathrm{Mm}\) with a resolution of \(23.46 ~\mathrm{km}\) horizontally and \(20 ~\mathrm{km}\) vertically. The convection zone extends roughly \(7 ~\mathrm{Mm}\) below the  \(\tau_{500}=1\)  surface, and the atmosphere covers approximately \(17 ~\mathrm{Mm}\) in height. We use the convention that a positive vertical component of the velocity corresponds to an upflow in the atmosphere.

Details about the se up of the simulation can be found in [24] and . Here, we give a brief summary of how the system reached the state in which it was used for the forward modeling. The simulation was initially set up as an ssd simulation with a horizontal extent of \(12 ~\mathrm{Mm}\times 12 ~\mathrm{Mm}\) [38]. The horizontal domain was extended by tiling it \(2 \times 2\) and running it for another \(8\) hours of simulation time to break the periodicity. After that, a bipolar magnetic feature with roughly \(8 ~\mathrm{Mm}\) separation between the poles was added, similar to the enhanced network model of the public Bifrost simulation [22]. The simulation was then run for another \(1.5\) hours after which the ne treatment of hydrogen was turned on. Following this, the simulation was run for another \(10 \min\) to let the populations settle into a new equilibrium. We then computed a time series of \(10 \min\), which we use for our analysis. The simulation is overall in a steady state but shows signatures of the typical \(5 \min\) oscillations in the photosphere and \(3\)–\(5 \min\) oscillations in the chromosphere. For a detailed analysis of the oscillations, a longer time series is required.

We use the RH1.5D code [39], [40] to synthesize the  \(\lambda 854.2\) nm spectral line. In this code, each vertical column in the simulated atmosphere is treated independently as a plane-parallel atmosphere. This is also called the 1.5D rt approach. We perform the spectral line synthesis at a heliocentric angle \(\theta\) such that \(\mu=\cos(\theta)=1\) corresponding to a disk center observation. By doing so, the los components of vector quantities such as velocity or magnetic field correspond to the vertical component in the simulation grid.

We perform three sets of spectral synthesis. Firstly, following [16] we treat the all six stable isotopes of calcium, \(^{40}\)Ca,\(^{42}\)Ca,\(^{43}\)Ca,\(^{44}\)Ca,\(^{46}\)Ca,\(^{48}\)Ca as separate atoms that are all treated in nlte in the rt computation, hereafter im. For the rt computations five-level-plus-continuum model atoms of were constructed similar to [16] based on the experimental data from [41] and [42].

Secondly, we used a cm. This is an approximation to the first computation with the advantage that only one atom needs to be computed in nlte to solve the rte. The effect of multiple isotopes is mimicked by modifying the absorption coefficient of the model atom. This is achieved by summing up Voigt profiles centered at the rest wavelengths of the corresponding isotopes. The sum is weighted by the relative abundance of the isotopes.

Finally, we used a five-level-plus-continuum atom model where it is assumed that the most abundant isotope in the solar atmosphere (i.e., \(^{40}\)Ca) is the only present calcium isotope in the atmosphere, hereafter sim. By doing so, we can test the effects of multiple isotopes on the line width and asymmetry by comparing it to the computation where only the most abundant isotope was considered. The computations were conducted in the approximation of crd, which is sufficient for the  \(\lambda 854.2\) nm spectral line [3].

3 Results↩︎

We present our results as follows: first, we describe the intensity images at \(500 \, \mathrm{nm}\) and at the line profile minimum, together with the temperature, the vertical component of the velocity, and the vertical component of the magnetic field at the corresponding formation height of the line profile minimum. We then show the spatially averaged spectra with corresponding line bisectors based on the three sets of rt calculations. Following this, we analyze how the spatially averaged spectrum and its asymmetry depend on different regions in the atmosphere that are selected by the local vertical magnetic field and vertical velocity in each column.

3.1 Intensity and atmospheric properties at the formation height of the line core↩︎

In the here presented model, the  \(\lambda 854.2\) nm line profile minimum forms approximately between \(0.5 ~\mathrm{Mm}\) and \(2.5 ~\mathrm{Mm}\), where \(z=0 ~\mathrm{Mm}\) is the average height of the \(\tau_{500}=1\) surface. Within these heights, the individual profiles form under a wide range of physical conditions. We therefore present in Fig. 1 an overview of the continuum intensity at \(\lambda= 500 \, \mathrm{nm}\), the line profile minimum intensity \(I_{\mathrm{min}}\), and formation height at the wavelength position of the minimum intensity. We determine the formation height by first finding the wavelength of the intensity minimum \(\lambda_{\mathrm{min}}\) at each single spectrum, and then locating the height position \(h\) in the atmosphere where the optical depth at this wavelength reaches unity, that is \(h(\tau_{\lambda,\mathrm{min}})=1\). The wavelength position of \(\lambda_{\mathrm{min}}\) is obtained in a window of \(\pm 0.5 \, \AA\) around the rest wavelength of the most abundant calcium isotope \(^{40}\)Ca. We note that the synthetic profiles in the line core can be complex with additional emission peaks. A determination of the line profile minimum in the line core is thus not always unique. In addition, we show temperature, vertical velocity, and vertical magnetic field at the formation height of the minimum intensity.

In the continuum intensity map (Fig. 1, panel a), the granulation patterns at the bottom of the photosphere are visible. Small bright structures reveal the presence of the network magnetic field in the simulation domain. The intensity map of the line-core minimum (panel b) shows web-like structures throughout the whole simulation domain, except above the strong network fields, where larger bright structures are visible. The web-like structure is a result of shock waves in the lower atmosphere that expand horizontally with height and interfere with each other. By compressional (adiabatic) heating, the temperature (panel c) is locally increased, leading to the visible pattern [43].

The temperature map (panel c) underlines the correlation between the line core intensity and shock patterns in the quiet regions at roughly \(0\leq y / ~\mathrm{Mm}\leq 5\) and \(17\leq y / ~\mathrm{Mm}\leq 24\). Above the strong network fields at roughly \(5\leq y / ~\mathrm{Mm}\leq 17\), the temperature is higher and is often associated with lower formation heights (panel d). Small, bright pixels can be seen as the formation height can change sharply in 1.5D rt. In addition, some intensity profiles show core reversals with multiple peaks, such that a clear definition of the line profile minimum is not always straightforward.

The formation height (panel d) reveals the corrugated surface over which the  \(\lambda 854.2\) nm line core forms in the muramche model. In the horizontal center of the simulation domain, the formation height indicates loop-like regions forming in the mid-to-upper chromosphere and shows fine structure. The fibrilar structures are not visible in the intensity image (panel b), presumably because of missing horizontal rt. A similar effect was found for the H\(\alpha\) line [44] and the h&k lines [45] in the Bifrost public snapshot.

The vertical velocity (panel e) indicates where the line core forms in upflows (blue) or downflows (red). Qualitatively, it can be seen that structures associated with upflows are more concentrated than downflows. The area coverage of rays, where the line core forms in a downflow, is \(57\%\) vs. \(43\%\) in an upflow. The average upflow velocities are $3.42 , $ vs. $-3.84 , $ in downflows. Both upflows and downflows reach maximum values of approximately $20 , $ magnitude. The higher magnitude of the downflow velocities results from lower densities [10], which also agree with the larger area (and volume) coverage of downflows in the chromosphere.

Finally, in panel (e), we show the vertical component of the magnetic field at the formation height of the line core. The network fields expand in the chromosphere , forming a canopy. The bipolar magnetic field patterns roughly coincide with the increased intensity pattern in panel (b). Small irregular features such as at \((x,y)\approx (8\, ~\mathrm{Mm}, 12.5\, ~\mathrm{Mm})\) result from ambiguities in the determination of the line profile minimum due to complex line core profile shapes.

3.2 Spatially averaged line profiles↩︎

In the following, we present the spatially and temporally averaged  \(\lambda 854.2\) nm spectral line profile calculated from a series of four snapshots, where we included all six isotopes of calcium (im, see Sect. 2). For one instance in time, snapshot muram_en_518000_503s², we show additionally the sim and cm computations. We compare the spectra to an average qs observation from the ftsatlas [9]. We begin by discussing the line profiles in comparison to the observation. After this, we discuss the asymmetry of the line profiles.

In Fig. 2 (a), we show the qs disk center ftsatlas spectrum together with synthetic computations in a large wavelength window around the line center. The temporal average consists of four spatially averaged spectra from snapshots separated by approximately \(2 \min\) of simulation time. The spatially averaged line profiles of these snapshots are shown in App. 6. In the far wings (i.e., \(|\lambda - \lambda_{\mathrm{^{40}Ca,rest}}| > 2.5 \, \AA\)), which form at the lower photosphere, the synthetic intensity approximately matches the observed spectrum and the differences between the computations are small. The equivalent widths of the synthetic spatially averaged profiles are \(2.76\,\AA\) (im), \(2.74\, \AA\) (sim), and \(2.71\,\AA\) (cm), and thus are rather similar. The slightly higher synthetic intensity compared to the observation might be due to the absence of blend lines in the rt computation. For a more precise analysis of the line core, we present in panel (b) a zoom-in (as indicated by the dashed box in panel a). Here, the differences between the im, cm, and sim computations from snapshot muram_en_518000_503s become more visible. We determined the line widths of the line core (as described in App. 7) of the plotted profiles which are \(0.47 \, \AA\) (ftsatlas), \(0.48\, \AA\) (im), \(0.49 \,\AA\) (cm), and \(0.46\,\AA\) (sim). Thus, the line width of the im and cm computations are close to the observed profile, and even \(2 \,\%\) (im) and \(4\,\%\) (cm) larger. The computation with only the most abundant isotope (sim) results in an approximately \(2 \%\) smaller line width compared to the im computation. The line width of the temporally averaged snapshots is \(0.48\, \AA\), and hence there is no significant difference between the time-averaged spectrum (red color) and that from a single snapshot (orange color). We note, however, that more snapshots are needed with a higher time resolution for a better comparison, which comes with additional computational costs.

We find that the line profile minima of all computed spatially averaged line profiles are redshifted with respect to the rest wavelength of the most abundant calcium isotope \(\lambda_{\mathrm{^{40}Ca, rest}}\). In the forward modeled spectra, the vertical component of the velocity at the formation height of the line profile minima is strongly correlated with the resulting Doppler shift of the line profile minima. Due to the imbalance of upflows (\(43\%\)) and downflows (\(57\%\)) at the formation heights of the line profile minima (see Sect. 3.1 and Fig. 1 e), the minimum of the spatially averaged spectral line is redshifted. A similar result was found by [13] who analyzed observations that suggest that the  \(\lambda 854.2\) nm line forms slightly preferred in downflows. While the redshift of the im computation from muram_en_518000_503s seems to be in good agreement with the observation (see Fig. 2 d), we would like to note that the line shift varies by a small amount of \(\approx 16 \, \mathrm{m\AA}\) with time. This difference can be seen in Appendix 6.

We now compare the line width from the spatially averaged line profile from snapshot muram_en_518000_503s with the averaged line width of single profiles. To this end, we calculated the line width for the im computation from single profiles and found an average value of \(0.37 \, \AA\), which is approximately \(23 \, \%\) lower than what we found for the line width of the spatially averaged profile (\(0.48\, \AA\), im). We checked to what extent the line width of the spatially averaged line profile is due to the Doppler shifts of single profiles. We therefore calculated, similar to [15], the standard deviation of the Doppler shift of the line profile minima \(\sigma_{\Delta v, \mathrm{min}}\) from all spectra. We obtained a value of $_{v, }=4.56 , $, which might explain the larger line width of the spatially averaged spectrum³.

This value is approximately a factor of four times higher than what [15] found in their forward modeled spectra and even higher than what these authors found from their observational data taken by the CRISP instrument [46] at the Swedish \(1 \, \mathrm{m}\) Solar telescope. This suggests that in the muramche model, the line width of the spatially averaged line profile partly results from single lines that are Doppler-shifted by the dynamic velocities in the simulated chromosphere.

While the line width of the cm computation is larger than in the im computation, the line profile of the cm computation is slightly narrower near the minimum intensity. The location of this difference in the line profile coincides with the rest wavelengths of the less abundant calcium isotopes [16]. This hints that the fundamental assumption used to construct the composite model atom breaks down in the simulated chromosphere, which is indeed the case as we show in App. 8.

3.3 Asymmetry of spatially averaged line profiles↩︎

The asymmetry of a spectral line can be visualized by studying its bisector. In its simplest way, the bisector of a spectral line can be constructed via \[\begin{align} b(I) = \frac{1}{2}(\lambda_{\mathrm{red}}(I) + \lambda_{\mathrm{blue}}(I)). \label{ref:eq95bisector} \end{align}\tag{1}\] Here, \(I\) is the intensity, \(\lambda_{\mathrm{red}}\) and \(\lambda_{\mathrm{blue}}\) are the wavelengths on the red or blue side of the spectrum with respect to the line profile minimum \(\lambda_{\mathrm{min}}\). In the literature, the bisectors of the  \(\lambda 854.2\) nm line are shown at different reference wavelengths. They can be shown on an absolute wavelength scale [10], they can be centered on the bisector intensity in the far wings [13], or they can be shown centered on the wavelength of the line profile minimum \(\lambda_{\mathrm{min}}\) as in [16].

The bisector centered on an absolute wavelength scale, or centered on the far wing, allows for comparisons of line Doppler shifts in the line core. The bisector as defined in [16] is given in terms of the wavelength relative to \(\lambda_{\mathrm{min}}\) and not on the absolute wavelength scale. This allows for comparison of the shape and amplitude of different bisectors irrespective of the absolute wavelength calibration of the observation. For compatibility, we show our results in all three representations in Fig. 2. As an additional measure of line asymmetry, we determine for each bisector the amplitude \(a_{\mathrm{core}}\), which is the distance in wavelength between \(\lambda_{\mathrm{min}}\) and the red-most excursion. We call the intensity of the line profile at which the maximum red excursion is reached \(I_{\mathrm{bs}}\). The amplitude \(a_{\mathrm{wing}}\) is the distance between the red-most excursion and the bisector at \(I_{\mathrm{c}}/2\), where \(I_{\mathrm{c}}\) is the continuum intensity. For an overview of these quantities, see App. 7 and Fig. 6.

In Fig. 2 (panels c, d, and e), the line bisectors of the spatially averaged profiles shown in panel (b) are presented. The ftsatlas line profile shows the typical ‘’inverse-C shaped” line bisector with an amplitude of \(a_{\mathrm{core}} = 14.98 \, \mathrm{m\,\AA}\) at an intensity of \(I_{\mathrm{bs}} = 11.86 \, \mathrm{\,n\, J\, m^{-2} \, s^{-1}\, Hz^{-1}\, sr^{-1}}\), the corresponding line wing amplitude of the bisector is \(a_{\mathrm{wing}}=34.47 \, \mathrm{m\,\AA}\). The line profiles computed with the im or cm show amplitudes of \(a_{\mathrm{core}} = 19.44 \, \mathrm{m\,\AA}\) (im) and \(a_{\mathrm{core}} = 14.74 \, \mathrm{m\,\AA}\) (cm) at intensities of \(I_{\mathrm{bs}} = 10.98\) (im) and \(I_{\mathrm{bs}} = 11.83\) (cm). Thus, in this particular snapshot, the amplitude and intensity of the cm computation match the observed bisector better than the full isotope computation (im). We note, however, the minimum intensity of the line profile in the cm computation is higher than in the im computation, which explains the shift of the bisector intensity. The wing amplitudes of the im and cm computations are \(a_{\mathrm{wing}} = 56.52 \, \mathrm{m\,\AA}\) (im) and \(a_{\mathrm{wing}} = 45.16 \, \mathrm{m\,\AA}\) (cm) which are both larger than in the observation. By taking only the most abundant isotope of calcium into account, the bisector of the sim computation does not show a visible red asymmetry, and hence no clear’’inverse C-shape”. The bisector from the time-averaged im computations results in an \(a_{\mathrm{core}} = 16.16 \, \mathrm{m\,\AA}\), which is closer to the observed value and indicates that the amplitude is time-dependent. The wing amplitude of the time-averaged spectra is the largest with a value of \(a_{\mathrm{wing}}=61.19\, \mathrm{m\,\AA}\), displaying that even in the wings of the lines, the average spectrum is sensitive to time averaging.

In panel (d), the bisectors are shown on an absolute wavelength scale. The wavelength axis is centered on the rest wavelength of \(^{40}\mathrm{Ca}\) to avoid large wavelength numbers in the plot. It can be seen that above \(20 \mathrm{\,n\, J\, m^{-2} \, s^{-1}\, Hz^{-1}\, sr^{-1}}\) the computations accounting for isotopic splitting agree with each other. The sim computation is different, which suggests that the isotopic splitting effect is still visible outside the line core.

In Fig. 2 panel (e) the bisectors are shown centered on the wavelength where the bisector intensity is \(60 \%\) of the continuum intensity, which is similar to the representation in [13]. In this representation, the maximum redshift of the im bisector is \(61.8 \, \mathrm{m\,\AA}\), which is larger than the \(34.3 \, \mathrm{m\,\AA}\) for the sim computation, highlighting again the isotopic splitting effect on the asymmetry of the spatially averaged line profile. The im computation compares approximately to the value found by [13], which was \(\approx 51 \mathrm{\,m\AA}\) (\(5.1\,\mathrm{pm}\)) in their quiet region and \(\approx58 \mathrm{\,m\AA}\) (\(5.8\,\mathrm{pm}\)) in their internetwork region. We note that we have not degraded our spectra with an instrumental PSF. In addition, the observed region of [13] might contain a different amount of magnetic flux than our simulation.

3.4 Dependence of line profiles on atmospheric conditions↩︎

In Sects. 3.2 and 3.3, we discussed the shape and asymmetry of the synthetic line profiles that were averaged over the whole computational domain in comparison with the observation. There are two mechanisms that contribute to the width and asymmetry of the line, namely the atmospheric structure and the presence of multiple calcium isotopes. In this section, we study how the atmospheric structure affects the line profiles. We study the shape and asymmetry of the  \(\lambda 854.2\) nm line by preselecting columns based on their atmospheric properties before computing the averages. We distinguish between columns based on their velocity structure and magnetic field strength along the los. For each vertical column in the atmosphere, we computed the average vertical velocity \(v_{\mathrm{avg}}\) between the formation height of the inner wing intensity and the line core, such that: \[\begin{align} v_{\mathrm{avg}} = \frac{\int_{z1}^{z2} v_z (z') \,\mathrm{d} z'}{z_2-z_1} \label{eq:average95velocity95formation95region} \end{align}\tag{2}\] where \[\begin{align} z_1 = \frac{z(\tau_{\mathrm{wing,v}}=1) + z(\tau_{\mathrm{wing,r}}=1)}{2} \end{align}\] and \[\begin{align} z_2 = \max_{\{\lambda\; | \; ||\lambda - \lambda_0|| = 1 \AA\}}z(\tau_{\lambda}=1). \end{align}\]

Here \(v_z\) is the vertical velocity in the atmosphere, \(z_1\) is the average formation height of the wing intensity at \(\pm 1 \AA\) from the rest wavelength \(\lambda_0\) of \(^{40}\mathrm{Ca}\), and \(z_2\) is the maximum formation height in the wavelength window \(\pm 1 \AA\) around the rest wavelength. We compute the averaged unsigned vertical magnetic field \(|B_z|_{\mathrm{avg}}\) analogously to the average vertical velocity.

We categorize columns with \(|B_z|_{\mathrm{avg}}<60 ~\mathrm{G}\) as having weak los magnetic field and \(|B_z|_{\mathrm{avg}}>60 ~\mathrm{G}\) as having strong los magnetic field. This separation roughly distinguishes between columns outside and inside the large-scale network field (see Fig. 1, panel f). Similarly, we distinguish between columns having \(v_{z,\mathrm{avg}}>0\) as upflows and \(v_{z,\mathrm{avg}}<0\) as downflows. The different categorizations between vertical motions and magnetic fields are not meant to form stochastically independent sets.

In Fig. 3 we show how the average profiles and bisector shapes vary depending on magnetic field (panels a, b, and c) and vertical velocity (panels d, e, and f). The bisectors are presented once relative to the rest wavelength of \(^{40}\)Ca (panels b and e) and once relative to the line profile minimum (panels c and f). In Fig. 3 (a), it can be seen that the line profiles of the wfp all appear similar to their corresponding profiles in Fig. 2. This is expected because approximately \(93 \%\) of all profiles in this snapshot are formed in the weak field regions. The sfp, which are averaged over approximately \(7 \%\) of all spectra, have higher intensities across the whole wavelength range shown. The continuum at \(10 \, \AA\) away from the line core is at roughly \(41 \mathrm{\,n\, J\, m^{-2} \, s^{-1}\, Hz^{-1}\, sr^{-1}}\), which is comparable to the wfp. In App. 9 we show the line profiles of the im computation over a larger wavelength window. The corresponding bisectors of the sfp in Fig. 3 (b and c, solid lines) all show an inverse “C’’ shape and a blueward turn close to the wing intensity, irrespective of the different treatments of isotopic splitting in the rt computation.

In the bottom panels (d), (e), and (f), we show the resulting ufp in solid lines and dfp in dashed lines, see Eq. (2 ). In this snapshot, we find \(45.5\%\) of the profiles form in an upflow and \(54.5\%\) form in a downflow. These numbers compare approximately to the vertical velocity at the formation height of the line profile minimum (see Sect. 3.1 and Fig. 1 e).

In Fig. 3 (d), it can be seen that ufp, besides being blue-shifted, have slightly higher intensities in the minimum. The differences in the ufp between the im and sim rt computations are negligible bluewards of the line core, but there are significant differences redwards. The corresponding bisectors of the ufp (panels e and f, solid lines) are blue shifted (panel e) and show a red asymmetry close to the line core with respect to the line profile minimum (panel f). Even for the sim (green curve), a slight red asymmetry in the line core is visible. In the outer wings, that is for intensities above \(\approx 18 \mathrm{\,n\, J\, m^{-2} \, s^{-1}\, Hz^{-1}\, sr^{-1}}\), the averaged ufp become approximately symmetric around the rest wavelength of \(^{40}\)Ca (panel e).

Unsurprisingly, the dfp are red-shifted with respect to the rest wavelength of \(^{40}\)Ca. Similar to the ufp, in the dfp, the different rt computations are in good agreement bluewards to the line center and differ redwards. The corresponding bisectors (dashed lines) are redshifted (panel e), and indicate a much lower red asymmetry with respect to the line profile minimum (panel f). While the im computation shows a much lower bisector amplitude compared to the spatially averaged line profile computed over the whole simulation domain (see Fig. 2 c), the bisector from the sim remains blueshifted at all intensities with respect to the line profile minimum.

This study can be summarised as follows: The strengths and shapes of the line profiles and their corresponding bisectors can depend strongly on regions selected by the los velocity or the los magnetic field strength in the atmosphere. In strong los magnetic field regions, even the averaged profile of the computation with only the most abundant isotope (sim) shows an inverse “C’’ shape. But in regions of weaker los magnetic field or downflows, the sim computation does not result in an inverse C-shaped line bisector, that is, when isotopic splitting is not taken into account.

4 Summary and discussion↩︎

In this work, we presented spatially averaged synthetic spectra of the  \(\lambda 854.2\) nm line, computed from an muramche rmhd model, which we compared to an average qs observation from the ftsatlas. Including a large-scale bipolar magnetic feature, the atmosphere model resembles an enhanced network region. We computed three different sets of synthetic spectra in 1.5D rt. We calculated the spectra once taking all abundant isotopes of calcium in the solar atmosphere into account, with the transitions due to the various isotopes treated as blend lines (im), and once where only the most abundant isotope was considered (sim). In addition, we utilized an approximate model of isotopic splitting that combines the effects of isotopes in only one model atom (cm).

The two-dimensional image of the line-core intensity (obtained from the im computation) shows a typical shock-expansion pattern [43] in the quiet part of the simulation. Above the network magnetic field, the intensity is enhanced. In the intensity image, no fibrilar structures are seen, while they are visible in the formation height of the minimum intensity. This is likely attributed to the lack of horizontal rt similar to the H\(\alpha\) line [44]. The formation height of the line core shows large spatial variation with fine-structured details.

The vertical component of the velocity at the formation height of the line center shows thinner upflow structures that are surrounded by expanded downflow structures. This suggests upward mass flows in thin dense channels, while the mass falls back in blobs of lower density, which cover a larger volume in the atmosphere.

We then studied the spatially and temporally averaged profiles for the im computation and spatially averaged profiles from one snapshot for the im, sim, and cm computation. We find a good match between the spatially and temporally averaged im computation and the observation in terms of the line profile depth, width, and shape, including the asymmetry as expressed by the line bisector. When only the most abundant isotope of calcium is considered in the synthesis, the resulting average spectrum is slightly narrower and almost symmetric. These findings confirm the results of [16], who suggested that isotopes must be taken into account to reproduce the observed inverse C-shape bisector. The CM calculation similarly reproduces the overall line width and asymmetry of the IM calculation. In the innermost line core, the intensity in the cm computation is, however, slightly higher redwards of the line profile minimum. In the composite model approach, it is assumed that the ratios of population-level densities from different isotopes are constant throughout the whole atmosphere. In App. 8 we show that this approximation breaks down already in the lower to mid chromosphere where the  \(\lambda 854.2\) nm line forms.

We found that the line width of single profiles is, on average, approximately \(23\%\) smaller than the width of the spatially averaged line profile. Similar to [15], we measured the standard deviation of the Doppler shifts of the intensity minimum of the line profiles. We found a value of ${v, }=4.56 , \(, which might explain the difference between the line widths of single profiles and the width of the spatially averaged line profile. The obtained\){v, }$ is approximately four times higher than in [15]. This suggests the close match of the line width in the muramche simulation with the observation is due to the dynamic atmosphere. [15] compared their results with observations from a coronal hole, which showed a value of $_{v, }=2.2 , \(. The reason for the roughly two times higher value we find in the [muramche]{acronym-label="muramche" acronym-form="singular+abbrv"} simulation could be similar to the result of [@{2018ApJ...864...21K}] who found a larger line width of the k line in the [qs]{acronym-label="qs" acronym-form="singular+short"} compared to a coronal hole. This suggests the simulation presented here compares well to the [qs]{acronym-label="qs" acronym-form="singular+short"} but a different simulation might be needed for the comparison with coronal hole observations. The observations from the third flight of the Sunrise balloon-borne observatory [@{2025SoPh..300...75K}] will provide superior seeing-free observations of the  \)$ nm line arising from different regions on the Sun. A comparison between these observations and the computations presented here and future muramche simulations will help to better understand the dynamics of the chromosphere.

We also considered how the shape of the bisector not only depends on whether all isotopes of calcium are included but can also be strongly influenced by atmospheric conditions. To this end, we separated the atmosphere into four sets of columns with strong or weak vertical magnetic fields and up- or downflows measured over the formation height range of  \(\lambda 854.2\) nm. These sets are not disjoint but still show different features in the average spectra. We found that average spectra over columns where \(|B_{z}|_{\mathrm{avg}} < 60 ~\mathrm{G}\) (\(93\%\) of the spectra) match the observed ftsatlas spectra best in terms of the line shape and the amplitude of the bisector. When the average is taken over profiles from columns where \(|B_{z}|_{\mathrm{avg}} \geq 60 ~\mathrm{G}\), the line intensity is weakened, and even the sim computation shows an ’’inverse C shape” which is comparable in amplitude to the other synthetic computations. The higher intensity is a result of the higher temperature in the regions dominated by a stronger magnetic field. The asymmetry must be a result of the flow structure above the magnetic features.

The separation into up- and downflows results in multiple effects. The asymmetries of the line shapes are substantially different, as measured by the line bisectors. The averaged profiles from columns showing downflows are much less asymmetric at low intensities, that is, close to the line core. The profiles averaged over regions showing on average upflows show a red asymmetry close to the line core.

The comparison presented here between different atmospheric parameters and the resulting spectral line properties is intended to understand the effect of isotopic splitting at different locations in the simulation domain in an average sense. [47] demonstrated that line profiles from single resolution elements in a simulation can be more complicated than individual observed profiles. This trend even holds after the synthetic profiles were degraded to match the instrumental conditions. In a follow-up study, we plan to do a detailed comparison between the simulation and high-resolution observations.

5 Conclusions↩︎

We could show that the muramche code can produce an atmosphere model that results in a close match with observations. For this purpose, we forward modeled the chromospheric  \(\lambda 854.2\) nm spectral line and found a good agreement in comparison with observations from the ftsatlas. The reasonably good match of the line width is a result of the dynamic velocity field at chromospheric heights in the muramche model. The line profile minimum is shifted towards the red with respect to the rest wavelength and approximately matches the position of the observed ftsatlas line profile. In the muramche model, the net shift is a result of an imbalance between up and downflows in the chromosphere at the formation height of the  \(\lambda 854.2\) nm line core. Whether a similar trend is apparent in high-resolution observations will be investigated in a follow-up study with Sunrise III observations.

We confirmed the need for multiple calcium isotopes in the forward modeling to reproduce the observed inverse C-shape asymmetry of the line profile. The asymmetry in the spectral profile produced by isotopes can, however, be mimicked by other atmospheric properties such as the magnetic field and the velocity structure. As a consequence, if the effects of multiple isotopes are not taken into account during the inversion of the  \(\lambda 854.2\) nm line profile then the inversion code will likely compensate for the missing isotopes by tweaking velocity and/or magnetic field structures and thus introducing errors in the inferred atmosphere.

Overall, we conclude that the model atmosphere must be sufficiently dynamic and the isotopic splitting effect must be taken into account in the forward modeling to match the observed  \(\lambda 854.2\) nm line profile.

6 Time variation of spatially averaged line profile↩︎

The ftsatlas observation represents an average region of the quiet Sun that is averaged over space and time. In the main text in Sect. 3.2 we discussed the spatially averaged line profile from one snapshot of the simulation, and a profile that was averaged over four snapshots. To provide more context on the temporal variation of the spatially averaged line profiles from these four snapshots, we present them here explicitly in Fig. 4 on an absolute wavelength scale. The four snapshots are separated by approximately \(2 \min\) of simulation time. The snapshot at \(2.06 \min\) corresponds to snapshot muram_en_518000_503s, which is discussed in detail in the main text. The red curve in Fig. 2 corresponds to the temporal average of the four snapshots shown in Fig. 4. It can be seen that even outside the line core, which is roughly above \(20 \mathrm{\,n\, J\, m^{-2} \, s^{-1}\, Hz^{-1}\, sr^{-1}}\), the bisector shows time variations. This suggests that the simulated lower atmosphere displays oscillations. Using the bisector of these far wings as a reference might thus be time-dependent. For a reasonable estimate of such a reference wavelength, many more synthesized snapshots are required.

7 Determination of line width and bisector parameters↩︎

We determine the line width following the description in [48]. First, we determine the wavelength position of the profile minimum \(\lambda_\mathrm{min}\). From that wavelength position, we take the average intensity in the line wings at the positions \(\lambda_{\mathrm{wing}}\) that fulfill \(|\lambda - \lambda_{\mathrm{min}}| = 0.6 \AA\). The line width is then given by the wavelength difference of the intensity profile at half the intensity between the minimum and the wing intensity. This technique was used instead of just the fwhm to make sure that we are considering the width of the line core, which is the important part of the line for the chromosphere. In this way, we avoid the influence of the prominent line wings. The measured line widths are indicated in Fig. 5.

We describe the bisector by parameters similar to those used in [49] and [50]. The bisector amplitude in the line core \(a_{\mathrm{core}}\) is defined as the difference in wavelength between the profile minimum \(\lambda_\mathrm{min}\) and the red-most excursion of the bisector. The intensity at the red-most excursion is named \(I_{\mathrm{bs}}\). The wing amplitude \(a_{\mathrm{wing}}\) of the bisector is defined as the difference between the red-most excursion of the bisector and the bisector at an intensity of \(\approx I_{\mathrm{c}}/2\) where \(I_{\mathrm{c}}\) is the continuum intensity at \(10 \, \AA\) away from the profile minimum. An overview of the measured bisector amplitudes is shown in Fig. 6.

8 The composite model atom↩︎

Here, we briefly describe that the fundamental assumption behind the composite model is, in general, not fulfilled in the simulated chromosphere.

In the composite model, the absorption coefficient is constructed by summing up Voigt profiles centered at the corresponding rest wavelengths of the different isotopes. The sum is weighted by the relative abundance of the isotope with respect to the most abundant isotope \(^{40}\)Ca. This approach is an approximation because it is assumed that population-level ratios between different isotopes are constant throughout the line formation region [16], [51]. [16] discussed that this is only an approximation, but can be justified by insignificant differences from the full computation. As shown in Fig. 2, we find in the muramche model differences towards the red from the line profile minimum between the im and cm computation.

In Fig. 7, we compare the population ratios in a cut through the atmosphere. We compute this ratio for the two most abundant isotopes, which are \(^{40}\mathrm{Ca}\) and \(^{44}\)Ca, and show that the ratios of the populations in the atmosphere deviate with height. We computed the ratio of the population-level densities (here for the lower level “l” of the transition) via: \[\begin{align}\delta_{40,44}^{\mathrm{l}}(z) = \frac{n_{\mathrm{^{40}Ca,\,l}}(z)}{n_{\mathrm{\,^{44}Ca,\,l}}(z)}. \end{align}\] Here, \(n\) is the population-level density of the respective calcium isotopes. The population-level ratio for the upper level (not shown) is computed analogously and shows a similar behavior. We normalized the height-dependent ratio by the ratio at \(z=0\). In the lower atmosphere, at about \(z\leq 0.5 ~\mathrm{Mm}\), \(\delta_{40,44}^{\mathrm{l}}\) is roughly constant. Higher up in the atmosphere, but still below the formation height of the line profile minimum (as indicated by the orange curve), the \(\delta_{40,44}^{\mathrm{l}}\) starts to deviate. The assumption on which the composite model is based breaks down at these heights in the simulated chromosphere.

9  \(\lambda 854.2\) nm line with continuum↩︎

In the main text (Sect. 3.4), we presented the intensity profile averaged over regions depending on the vertical component of the magnetic field averaged over the height of formation between the inner wings and the line core. Here, we present an overview of how the line profiles compare in a larger wavelength window containing the continuum. Figure 8 shows the intensity of the ftsatlas observation and the averaged intensities over regions, which have \(|B_z|_{\mathrm{avg}}<60 ~\mathrm{G}\) (dashed lines) those averaged over regions which have \(|B_z|_{\mathrm{avg}}\geq60 ~\mathrm{G}\) (solid lines, see Sect. 3.4). While the line core intensity of the sfp profile is clearly enhanced, the difference between sfp and wfp profiles is small in the continuum.

References↩︎