© The Authors 2025

1. Introduction

The solar chromosphere is a highly dynamic and structured layer of the Sun’s atmosphere, bridging the gas pressure-dominated photosphere and the magnetically dominated corona. Its intricate interplay of radiation, magnetic fields, and plasma makes it a challenging region to study (Carlsson et al. 2019). Among its most striking features are spicules, rapid blueshifted excursions (RBEs), and rapid redshifted excursions (RREs), which potentially play a crucial role in transporting mass and energy to the upper solar atmosphere.

A key tool for studying chromospheric dynamics is the Hα line, one of the most widely used spectral lines for solar observations (Rutten 2007, 2008). Its strong absorption features on the solar disc and emission features near the limb make it well suited for imaging and spectroscopic techniques. The morphology and evolution of several large features such as filaments and prominences are best observed in Hα. In its wings, Hα works well to detect fine structures such as spicules on much smaller length scales – down to several hundred kilometres (Secchi 1871; Beckers 1968, 1972; Pikel’Ner 1971).

Spicules appear as elongated, fine multi-threaded strands that abound in the solar chromosphere at the limb. They are easily observed in the wings of the Hα and Ca II lines (De Pontieu et al. 2007; Martínez-Sykora et al. 2009). These dynamic features invite interest because their ubiquity and dynamic motion make them good candidates for mass and heat transport to the upper solar atmosphere. Spicules are categorised into types I and II. Type I spicules are observed to move up and fall back down with apparent velocities of 15–40 km/s in height-time plots, tracing a parabolic trajectory (Rouppe van der Voort et al. 2009; Pereira et al. 2012). They exhibit lifetimes of 3–10 minutes and are driven by magneto-acoustic shocks (De Pontieu et al. 2004). The on-disc counterparts of type I spicules are the dynamic fibrils. Type II spicules are more dynamic, showing rapid evolution with lifetimes of a few seconds up to 4 minutes. Some of the type II spicules are seen to disappear with high velocities (80–300 km/s), with nothing falling back to the surface after the rising phase. The on-disc counterparts of type II spicules are identified as RBEs ((Langangen et al. 2008; Rouppe van der Voort et al. 2009; Sekse et al. 2013a; Shetye et al. 2016, and others)). The RREs share similar morphologies and statistical properties to the RBEs, appearing in the red wing of Hα and Ca II lines. The RREs could be return flows of type II spicules from the transition region (Chaurasiya et al. 2024).

The rapid excursions (RBEs and RREs) are concentrated around the magnetic network and plage, with their footpoints rooted near strong magnetic field concentrations (Bose et al. 2021; Danilovic et al. 2023). First observed in the Ca II IR line (Langangen et al. 2008), they are characterised by lifetimes of up to a few minutes, linear trajectories, upflow speeds for RBEs (20–50 km/s), and distinct lengths (up to ∼6 Mm) and widths (a few hundred kilometres) in observations. Although we have a grasp on the statistical properties of type II spicules, RBEs and RREs, our understanding of their driving mechanisms and energetics remains incomplete (Pereira et al. 2012; Sekse et al. 2013b).

Spicules are particularly intriguing due to the velocities that govern their motion. Apparent velocities of type II spicules can be exaggerated by line-of-sight density effects, while actual plasma motions may be much slower. The use of Doppler velocities further complicates the interpretation, as rapid flow direction changes can cause features to disappear from one spectral passband and reappear in another (Pereira et al. 2016). The complex transverse and torsional motions observed in rapid excursions (Sekse et al. 2013b) also mean that RBEs can switch to RREs and vice versa, depending on the line of sight and the direction of plasma flows.

Simulations are an important tool to help us interpret the observations. Synthetic Hα images can be generated from the simulations to look for similar events to those which are observed. In simulations we have the full 3D structure which is important for understanding the dynamics. One approach to synthesising Hα, given a background atmosphere, is to assume that the hydrogen population levels are in statistical equilibrium with the non-local thermodynamic equilibrium (NLTE) radiation field (Leenaarts et al. 2012). This implicitly assumes that the NLTE radiation field and populations reach equilibrium on timescales shorter than those of the dynamics. The resulting problem then involves solving simultaneously for the radiation field and populations: this is the approach adopted by codes such as RH and MULTI3D (Uitenbroek 2001; Leenaarts & Carlsson 2009).

An alternative approach is to use the populations calculated by MURaM-ChE using the rate equations (Przybylski et al. 2022). The rate equations implemented in MURaM-ChE include the dynamical effects but include only a very approximate treatment of the NLTE radiation field (Sollum 1999). These MURaM-ChE populations can be used in the Hα synthesis in several ways with different levels of sophistication. One approach investigated by Krikova et al. (2023) is to synthesise Hε using a modified version of the RH code that includes the non-equilibrium hydrogen ionisation from a Bifrost simulation. In this paper we introduce the simplest way. We show that our approach with MURaM-ChE and MULTI3D synthesis (which assumes instantaneous equilibrium between the populations and the NLTE radiation field) leads to the identification of the same features, which is sufficient for our study.

Numerical modelling of solar spicules dates back to the late 1960s (see the review by Sterling 2000). Multiple models describing the driving mechanism of spicules have been proposed such as Alfvén waves, pressure pulses in the high chromosphere, shocks, etc. More recent simulations performed by Martínez-Sykora et al. (2017) outline the role of ambipolar diffusion in generating a solar spicule via the build-up of magnetic tension. There have been conjectures regarding propagating heating events (PHEs) associated with the formation of type-II spicules, RBEs and RREs (Rutten 2017). Using Hα observations from the Goode Solar telescope at the Big Bear Solar Observatory Samanta et al. (2019) show evidence of magnetic reconnection driving spicules and subsequent coronal heating. Using radiative-magnetohydrodynamic (radiative-MHD) simulations Druett et al. (2022) demonstrate how the destabilisation of low-lying twisted magnetic field lines and Lorentz forces contribute to the formation, mass-loading, and drainage of dynamic fibrils. Danilovic et al. (2023) compute the Hα line in 3D radiative-MHD simulations to reproduce RBEs and RREs and assert that Alfvénic waves are the primary driver. Studies also suggest that the apparent motions of these jet-like features correspond to propagating heating fronts rather than actual mass motions (De Pontieu et al. 2007). In a study by Bose et al. (2021), the heating in downflowing RREs is ascribed to ambipolar diffusion (Martínez-Sykora et al. 2017).

Synthetic observables are crucial for interpreting spicules, RBEs, and RREs in simulations. There have only been three attempts aimed at modelling RBEs and RREs: Kuridze et al. (2015), Srivastava et al. (2017), and Danilovic et al. (2023). These studies also show that it is not straightforward to model these on-disc fine structures such that they reproduce properties (e.g. lengths, lifetimes, and motion) as seen in observations. For our purpose, we develop a synthetic observable based on an escape probability (Hummer & Rybicki 1982) for the Hα spectral line.

In this work we introduce a Doppler-shifted Hα proxy (Sect. 3), which was computed using the populations of the different levels of hydrogen calculated in a 3D MURaM-ChE simulation. These populations were calculated using the rate equations and include information from the history of the plasma. This is a different approach from assuming instantaneous statistical equilibrium which is usually assumed in radiative transfer models such as RH and MULTI3D. Because we do not solve for the population levels, under the assumption of statistical equilibrium with the radiation field, the proxy synthesis is computationally cheap (a few CPU seconds for each snapshot) and can be done on the fly with MURaM-ChE with almost no additional costs. Despite the differences in the two approaches, we get very similar results in the morphology of on-disc Hα wing features for both the proxy and the MULTI3D synthesis for the two snapshots we study. The proxy acts as a tool for identifying fine-structures that can then be studied in various ways. Here we concentrate on a proof of concept based on a detailed study of a single feature. A statistical investigation of many such features in a MURaM-ChE simulation will be presented in an upcoming paper.

In this study we focus on the driving mechanism of an example RBE feature and associated heating. The properties of the feature are evaluated and compared to those of a typical RBE, as seen in observations (see Sect. 4). Section 4.1 focuses on the forces and flows driving the feature during its lifetime. The identified feature exhibits jet-like behaviour. The origin of this jet is discussed. We also observe an oscillatory behaviour during its evolution. Finally, we provide further discussions and the conclusions in Sect. 5.

2. Model atmosphere and synthesis

2.1. Enhanced network model with MURaM-ChE

The MURaM radiative-MHD code (Vögler et al. 2005; Rempel 2017) models the near-surface convection, photosphere, chromosphere, and corona in 3D Cartesian geometry. With the chromospheric extension to MURaM – hereby MURaM-ChE (Przybylski et al. 2022), the code now includes an equation of state that treats hydrogen level populations in non-equilibrium (NE). This treatment involves solving the NE rate equations for hydrogen and is essential to forward model the formation of Hα.

A simulation of an enhanced network region, performed using MURaM-ChE is used for this study. The grid dimensions are 24 Mm  ×  24 Mm  × 24 Mm with a resolution of 23.46 km  × 23.46 km  × 20 km in the x, y and z (vertical) directions, respectively. The simulation spans from approximately 7 Mm below the photosphere (at z ≈ 0 Mm), to 17 Mm above. It has a bipolar magnetic field configuration aimed to model a network field patch on the Sun (the same model is described in greater detail in Ondratschek et al. 2024). We used 101 snapshots spanning 606 s, with a cadence of ∼6 s, in order to perform the analysis. An overview of the magnetic field at different heights is shown in panels (a) and (b) of Fig. 1. An enhanced magnetic network with the bipolar configuration is clearly visible at the photosphere, and it extends into the chromosphere. In panels (c) and (d) of Fig. 1 we show synthetic Hα intensity maps generated using MULTI3D (see Sect. 2.2).

Fig. 1. Overview of the magnetic field and synthetic Hα intensity in a snapshot of the enhanced network simulation. In panels (a) and (b) we show the vertical component of the magnetic field (Bz) at the surface (τ500 = 1) and in the chromosphere (z = 1 Mm), respectively. Panels (c) and (d) are synthetic Hα intensities in the blue (37 km/s) wing and red (−37 km/s) wing, respectively. vD in the title of panels (c) and (d) denotes the Doppler shift.

2.2. Synthesis with MULTI3D

We used MULTI3D (Leenaarts & Carlsson 2009) with the modelled atmosphere to forward model the Hα spectral line. The MULTI3D code solves the radiation transfer problem in NLTE on a 3D Cartesian grid. We treated the Hα line under the approximation of complete-frequency-redistribution (CRD). The viewing angle was set to μ = 1, to study on-disc features. Every second column was sampled in the two horizontal directions to make the computation faster (see for example Danilovic et al. 2023). The vertical height limit for this computation was set from −1 Mm to 11 Mm. Most of the cool chromospheric material lies below 10 Mm. This determined our upper limit for the synthesis. The computation for Hα was done using a five-level plus continuum hydrogen model atom (Leenaarts et al. 2012). MULTI3D solves the RT equation over multiple frequency points (order of 10² points), with the assumption of statistical equilibrium.

3. Identifying fine structures in Hα

We focus on identifying features in the wings of the average Hα spectral line within the domain where Hα is synthesised. We accounted for the NLTE physics essential for Hα formation using MURaM-ChE. We then used MULTI3D to forward model the Hα line from the MURaM-ChE atmosphere. This process is computationally expensive (∼10⁵ CPU hours per snapshot). To study temporal evolution efficiently, we employed a faster approximation that generates proxy Hα maps in seconds using NE hydrogen populations, enabling feature identification in the line wings.

3.1. Description of the proxy: Two-level hydrogen model atom

The Hα proxy is designed to be modelled based on the probability of photon escape from the photosphere and chromosphere. On the solar disc Hα is in absorption (n = 2 → n = 3 transition). MURaM-ChE computes the hydrogen populations for a five-level plus continuum hydrogen model atom. In our fast method we used a two-level approximation (involving energy levels: n = 2, n = 3) of the hydrogen atom. We also assumed uniform, time-independent radiation coming from the photosphere. The simulated atmospheric parameters – temperature and velocity – were utilised to estimate the line profile function (Φ). Since the Hα line is broadened predominantly by the small-scale thermal and turbulent motions (Doppler broadening) on the Sun, Φ takes the form

$$\begin{array}{c}\hfill \mathrm{\Phi }({\mathsf{v}}_{\mathrm{D}})=\frac{c}{{\nu }_0}\sqrt{\frac{m}{2\pi k_{\mathrm{B}}T}}\exp (-\frac{m{(\mathsf{v}-{\mathsf{v}}_{\mathrm{D}})}^2}{2k_{\mathrm{B}}T}),\end{array}$$(1)

where ‘v’ is the line-of-sight velocity (v_z in our case). v_D is the Doppler velocity. T is the temperature, m is the atomic mass of hydrogen, and k_B is the Boltzmann constant.

For Hα we took the wavelength (λ₀) of 656.28 nm at the line core (in air), which corresponds to a frequency, ν₀. Following Rybicki & Lightman (1979) we computed a measure for the opacity in the Hα line from the absorption coefficient (κ), defined as

$$\begin{array}{c}\hfill \kappa =\frac{h{\nu }_0}{4\pi }(n_2{B}_{23}-n_3{B}_{32})\mathrm{\Phi }({\mathsf{v}}_{\mathrm{D}}),\end{array}$$(2)

where n₂, n₃ correspond to hydrogen populations in the (n = 2) and (n = 3) energy levels, respectively. B₂₃, B₃₂ are the Einstein coefficients. The first component (n₂B₂₃) corresponds to absorption, while the second term (n₃B₃₂) corresponds to stimulated emission. We have the following Einstein’s relations:

$$\begin{array}{c}\hfill \frac{A_{32}}{B_{32}}=\frac{2h{\nu }_0^3}{c^2},\&\frac{B_{23}}{B_{32}}=\frac{g_3}{g_2}.\end{array}$$(3)

Fig. 6. 2D histogram plots showing the classification of pixels. Panels (a), (b): 2D histograms of the proxy vs. the synthetic Hα intensity at Doppler velocity = 37 km/s (a) and −37 km/s (b). Panels (c), (d): The row normalised distribution of the same pixels corresponding to panels (a) and (b), respectively. We mark the dynamic threshold for the synthetic intensity with vertical dashed black lines. The static threshold for the proxy for feature identification is marked by the horizontal dashed black lines. The four pixel categories: true positive (TP), true negative (TN), false positive (FP), and false negative (FN) are also indicated on each of the panels.

For Hα, the Einstein’s coefficient A₃₂ = 4.41 × 10⁷ s⁻¹, g₃ = 18 and g₂ = 8 (see NIST¹). Plugging the computed Einstein’s coefficients in Eq. (2), we obtained the absorption coefficient. We then calculated the optical depth for Hα. Under the assumption that there is no scattering, the absorption coefficient equals the extinction coefficient. Note that while this assumption holds for the line wings, the line core is dominated by scattering. The optical depth for a distance, z, is given by

$$\begin{array}{c}\hfill \tau (z)={\displaystyle {\int }_{z=z_1}^{z=z_2}}\kappa dz,\end{array}$$(4)

where dz is the length element in the line of sight. z₁, z₂ denote the vertical extent (top to bottom) over which the integration is performed. The expression κdz gives the probability that a photon will be absorbed when travelling through distance dz. Let P(z) be the probability of a photon traversing a distance z without getting absorbed. Then it can be shown that

$$\begin{array}{c}\hfill P(z)=\exp (-\tau ).\end{array}$$(5)

This treats the absorption in Hα as a photon escape probability. For our simulation box, dz = 20 km. We chose z₁ = 6 Mm and z₂ = 0.7 Mm in Eq. (4) to cover heights where most of the fine-structures observed in Hα are formed (White & Wilson 1966). We used 6 Mm as an upper limit because integrating the atmosphere to higher heights does not affect our results. The average height of the temperature minimum lies close to 0.7 Mm. The average Doppler profile for this proxy is shown in Fig. 2. A shift of the proxy from Δv = 0 is visible in this figure. This shift varies in time with an approximate 3-minute period (based on a 10-minute time sequence).

Fig. 2. Spatially averaged profile of the Hα proxy, modelled as an escape probability (orange). The same quantity from the MULTI3D MURaM (-ChE) synthesis for comparison (purple).

In Fig. 2 we also show a comparison using the computed opacities (τ) from MULTI3D for the synthetic Hα spectral line to construct the photon escape probability, exp(-τ). The differences between the two escape probability profiles stem from several factors. MURaM-ChE accounts for the non-equilibrium hydrogen populations but assumes a very simple treatment of the radiation field (Sollum 1999). MULTI3D assumes statistical equilibrium for computing the hydrogen populations but has a more detailed treatment of the RT problem.

The proxy works as a physically motivated tool. It helps isolate Hα absorption features in the simulation domain based on the hydrogen populations at the n = 2 level (for absorption), temperature, and the line-of-sight velocities.

3.2. The synthetic Hα intensity maps

To test the performance of the proxy against 3D NLTE synthesis, we used MULTI3D to generate synthetic Hα intensity maps over a range of wavelengths – 6543.5 Å to 6591.3 Å. Two snapshots were synthesised with the computation described in Sect. 2.2; they are 12 seconds apart. Throughout Sect. 3 we show the analysis using only one snapshot. The results remain unchanged for the other snapshot.

We find numerous small-scale structures in the generated synthetic Hα maps (see panels c and d in Fig. 1) that appear to be rooted close to the network patches. The SST (Swedish Solar telescope) dataset observed on 25 May 2017 at 09:16:58 UTC shows such an enhanced network region (see for example Fig. 1 in Bose et al. 2021) with numerous RBEs and RREs. Our synthesised Hα maps appear to be similar to their observational counterparts.

We further compare the observed average Hα profile from the FTS Atlas of Neckel & Labs (1984) to the synthetic intensity average profile from our MURaM-ChE simulation (see Fig. 3). This helps to benchmark the synthetic Hα maps created using MULTI3D-MURaM against observations. The FTS Atlas observations correspond to quiet Sun regions spatially (and temporally) averaged over the disc centre whereas the simulation is that of an enhanced network. Nevertheless, the comparison shows a good match in the line profile, especially in the inner wings. The shift ( = 0.035 Å or 1.59 km/s) that we introduced to align the line core arises from two major sources. We calculated the gravitational redshift to be 0.63 km/s. Secondly, the average velocity at the formation height is a downflow of 0.88 km/s. These factors introduce a slight shift of the line towards the red.

Fig. 3. Spatially averaged Hα line profile from the FTS Atlas of Neckel & Labs (1984) (green) compared with the Hα profile synthesised by applying MULTI3D to the simulated MURaM-ChE atmosphere (purple). The dashed black line marks the centre of the line core. We have shifted the synthetic profile (MULTI3D MURaM) by 0.035 Å(≈1.59 km/s) towards the blue wind to align the line cores for comparison.

In the far wings of Hα, collisional or pressure broadening becomes significant. Consequently, a Voigt profile is generally more suitable for modelling the line than a simple Doppler profile. We fitted the average Hα line profile from the MULTI3D synthesis with the best-fit Voigt profile and Gaussian profile (see Fig. 4). The Gaussian fit aligns more closely with the MULTI3D average profile, while the Voigt profile appears broader in the far wings, as expected. This supports the approximation of using the Doppler profile to generate our proxy. The maximum discrepancy between the Voigt and Gaussian fits is on the order of 10⁻¹.

Fig. 4. Average escape probability profile from the MULTI3D synthesis on MURaM-ChE (black circles) fitted with a Voigt profile (maroon) and a Gaussian profile (teal).

3.3. Results: Goodness of the proxy

The Hα proxy effectively identifies the location and morphology of chromospheric absorption features at μ = 1 in the line wings (Fig. 5) based on the MULTI3D synthesis. However, it becomes less reliable closer to the line core, where numerous feature pixels saturate the maps and the structuring of individual features is no longer visible (see movie associated with Fig. 5). This arises because of the simplifying assumptions that went into the construction of the proxy (see Section 3.1). Despite the simplifications, our assumptions are sufficient for the scope of this study to capture absorption features in the Hα wings.

Fig. 5. Features in the blue wing at a Doppler velocity of 37 km/s. The leftmost panel shows the average line profile for the synthesised Hα intensity (using MULTI3D) with wavelength. The vertical blue line marks the wavelength at which the Doppler velocity is 37 km/s with respect to the line centre (Δλ = 0). The middle panel shows the synthetic Hα intensity at this wavelength. The third panel shows the proxy at the same wavelength. A scan through the line is available online.

To quantify the goodness of our proxy, we performed a pixel-by-pixel comparison. This was done by generating 2D intensity distributions for the synthetic MULTI3D image and the proxy (see Fig. 6). We find a positive correlation in the pixel distribution shown in Fig. 6. This means lower synthetic intensities (with MULTI3D) correspond to lower photon escape probabilities in the proxy and vice versa. In the wings there are fewer dark absorption features and these occupy a smaller area compared to the background intensity coming from the photosphere. Therefore, the major concentration of the pixels (bright pink patch on the top right of the 2D intensity distributions in panels a and b in Fig. 6) indicate the background or photospheric contribution in the wings. We chose thresholds in the 2D intensity maps for the classification of pixels. The threshold for identifying a feature pixel in the synthetic intensity was determined by the background contribution, such that the pixels corresponding to the background (bright pink patch), were identified as negative pixels. This threshold was determined automatically based on the 2D intensity maps and moves with different positions in the line wings. On the contrary our proxy, being an escape probability, is limited to values from 0 to 1. So, we fix a threshold at 0.8. Any pixel having a value below the thresholds belongs to an Hα feature. Based on this we obtain four regions in the 2D intensity maps and classify the pixels as true positives (TPs), true negatives (TNs), false positives (FPs), and false negatives (FNs). The classification of pixels into the four categories is further shown in Table 1 and in Fig. 6.

After categorising the pixels for a given map, we quantified the goodness of our proxy. In Fig. 7 we show the quantification of the goodness at a particular Doppler velocity of v_D = 37 km/s. As we move from the wings towards the line core, the absorption features cover larger areas and the background becomes less prominent. So we have a higher percentage of TPs and fewer TNs closer to the line core and vice versa. From the scan in the blue wing (see movie associated with Fig. 7), it is clear that the proxy reasonably identifies blueshifted Hα on-disc features. We obtain consistent results in the red wing.

Fig. 7. Goodness of the proxy at Doppler velocity = 37 km/s. Panel (a) gives a spatial distribution of (in)correctly identified pixels. The pixels are colour-coded to show true positives (TPs: red), true negatives (TNs: white), false positives (FPs: green) and false negatives (FNs: blue). Panel (b) represents the percentage of (in)correctly identified pixels as a function of the Doppler velocity (vD) in the blue wing of the line. The dashed black line marks the position in the line wing (vD = 37 km/s) at which panel (a) is obtained. A scan through the blue wing of the line exhibiting the goodness of the proxy is available online.

We applied this method to both synthesised snapshots with the calculated synthetic Hα intensities (with MULTI3D) and obtain consistent results. Our proxy correctly identifies the Hα wing features observed in the synthetic images. It is also computationally efficient as it takes a few CPU seconds to obtain the proxy for each snapshot. Thus, we have a quick and robust proxy which can be used to study features that appear to be RBEs and RREs in Hα. With this proxy we can locate and identify synthetic observables directly from the MURaM-ChE simulations.

Over and beyond using the line-of-sight integrated proxy maps, we can also apply the proxy to vertical cuts. This is then the normalised absorption coefficient (κ). We find that this works well to identify at which heights the feature forms in the atmosphere. In a way this serves as a proxy to the formation height of the identified features. This enables a 3D perception of the feature. In the following section, we show the analysis of a particular blueshifted feature using the devised method.

4. Analysing a blueshifted feature

We now present the analysis of an example feature identified in the blue wing of our proxy in the same snapshot as used for the previous analysis. The feature is visible at a Doppler shift of 37 km/s. It appears as a straight feature at μ = 1 with a maximum length of 3.4 Mm and displays a very similar appearance in the image generated by MULTI3D at the same Doppler shift (see panels a and b of Fig. 8). We take a vertical slice (panel d in Fig. 8) through the feature to better gauge how the structure looks in 3D. Although it appears as a straight feature in the on-disc line-of-sight integrated maps, it has a loop-like morphology when viewed in the vertical cuts. Furthermore, the formation height obtained from the Hα MULTI3D synthesis (panel c in Fig. 8) shows that the feature forms at about 2–4 Mm above the average solar surface. This is consistent with what we see in the vertical cut of the proxy feature (panel d in Fig. 8). The proxy is used as the normalised absorption coefficient (κ) for the vertical slices, which shows the density structure with height. In panel (e) we find that the feature (overplotted in black) lies in an upflow region, which explains why it appears in the blue wing. In panel (f) we see that it appears to be rooted at a negative magnetic field patch from the vertical component of the magnetic field.

Fig. 8. Overview of one blueshifted feature (vD = 37 km/s) detected by the proxy. Panel (a) shows the proxy map with the feature of interest highlighted by the red box. The green line passing through the feature marks the base for the vertical cut through the feature (hereby x′ axis). The 0 of x′ corresponds to (x, y) = (8.6, 0) Mm in the top view (in panel a). Panels (b) and (c) are a zoomed-in view of the feature showing the synthetic intensity in Hα and formation height, respectively (obtained with the MULTI3D synthesis). In panel (d) we take a vertical cut through the feature in the proxy (normalised κ). Panels (e) and (f) show the same vertical cut for the line-of-sight (i.e. vertical) component of the velocity and magnetic field, respectively, with the feature over-plotted in black. A time-lapse of the proxy map in the top view and vertical cuts of various MHD quantities is available online.

The movie associated with Fig. 8 shows that the feature forms in a region outside the chosen vertical slice. The feature forms close to time = 0 s. We use t = 0 as the time for feature formation henceforth. As the feature evolves, it moves into this vertical plane and exits it again before disappearing. The lifetime of the feature is approximately 4 mins. While the feature is still in the plane, it is observed to expand outwards. It can be seen that any feature detected by the proxy at this Doppler velocity always lies in upflow regions and has cool, chromsospheric temperatures. Furthermore, the hydrogen (n = 2) population densities are locally enhanced in the region of these features.

A close inspection of the movie (in Fig. 8) reveals that as the feature enters the vertical slice, a positive polarity intrusion is visible in the same region, in a predominantly negative polarity field. This hints at a possibility of magnetic reconnection. We explore this further in the following Section 4.1.

4.1. Driving mechanism

We traced magnetic field lines close to the site of feature formation to better understand possible magnetic reconnection through the rearrangement of fields. The 3D rendering of magnetic field lines reveals an overarching bundle of negative polarity field lines that overlie the site of flux cancellation (see Fig. 9). Another bundle of field lines with opposite polarity is seen to bend and cross these overarching field lines from below. This is seen as the opposite polarity intrusion in the vertical cuts shown in Fig. 8. This secondary bundle of field lines connects closer to the main footpoint of the feature and are seen to undergo rearrangement as flux emergence and cancellation occurs there. This is highlighted by the red circle. The rearrangement of magnetic field lines at this site is visible in the movie associated with Fig. 9.

Fig. 9. 3D rendering of magnetic field lines rooted at the footpoint of the feature at t = 0 s, the first appearance of the feature. The line-of-sight component of the magnetic field is shown as a surface map at the τ500 = 1 surface. The 3D field lines are coloured by the vertical (Bz) component. The site of magnetic reconnection lies across the main footpoint and is marked by the red circle. The dashed green line shows the xy projection of the 3D field line carrying the feature in this snapshot. The magnetic field lines are clipped for better visualisation. A time-series showing the rearrangement of magnetic field lines at the time of feature formation is available online.

We further looked at the evolution of the line-of-sight (vertical) component of the magnetic field at different optical depths. In Fig. 10 we show the evolution at two τ₅₀₀ surfaces zooming into the feature. At τ₅₀₀ = 1, the evolution of photospheric fields show signatures of flux emergence and consequent cancellation. The dashed red box in the leftmost panel in Fig. 10 highlights the site where magnetic reconnection is likely to occur. This is the same region as marked by the red circle in Fig. 9. It is clear that the magnetic patch outlined by the dashed red box begins with a predominant positive polarity. A small negative polarity patch emerges close to it and cancels the pre-existing positive polarity magnetic field. Eventually, the positive polarity is completely cancelled and negative polarity fields remain in that region. This is also visible in the middle photosphere (τ₅₀₀ = 0.01). This indicates that flux cancellation and reconnection also extend into the low chromosphere. Additionally, regions of enhanced current density |J|, integrated along the line of sight, are highlighted through cyan shading. The threshold is set to 1.75 times the mean value. This analysis further bolsters the hypothesis that flux cancellation and consequent magnetic reconnection, in the photosphere and chromosphere, drives the formation small-scale features (Shibata et al. 2007; Samanta et al. 2019; Duan et al. 2023).

Fig. 10. Zoom-in on the evolution of the line-of-sight magnetic field at various τ levels in the region of the studied feature. Left: Evolution of the vertical component of the magnetic field, Bz at the photosphere (τ500 = 1). The feature from the proxy is overplotted in purple. The dashed red box indicates the region of flux cancellation. Centre: Evolution of Bz in the middle photopshere (τ500 = 0.01). The grey scale for the magnetic field (Bz) ranges from −500 G to 500 G at the τ500 = 1 surface and from −100 G to 100 G at the τ500 = 0.01 surface. The regions of enhanced current density (75% above the mean) are plotted with the cyan contours in both the left and central panels. Right: Zoom-in on the evolution of the feature in the proxy as viewed on disc. A movie spanning the lifetime of the feature is available online.

The movie associated with Fig. 10 also shows some physical characteristics of the feature. As mentioned before it has a maximum length of 3.4 Mm and lifetime of 246 s. It also exhibits a lateral displacement (leftward as viewed on-disc) of up to 400 km during its lifetime. This is consistent with observations of RBEs and RREs (Sekse et al. 2013b; Kuridze et al. 2015; Bose et al. 2021).

4.2. Tracing field lines containing the feature

To better follow the evolution of the feature, we trace one field line threading the feature in each snapshot. We emphasise that the field line is not being traced through time and it is not necessarily the same field line in different snapshots. Rather, we follow the 3D density feature in the proxy for each snapshot and trace the field line based on its location. The seed point in each snapshot for tracing the field line is chosen based on the proxy feature, such that it lies approximately at the centroid of said feature. The integration is then done bi-directionally to trace the field line. A trilinear interpolation scheme is implemented with the RK45 integration. Fig. 11 shows the feature in the proxy and various MHD quantities in the plane containing the traced field line. This plane has a curved base (dashed green line in Fig. 9) based on the grid points where the field line is traced. In our analysis, we treat the negative polarity footpoint as the starting point (L = 0). The other footpoint shifts and changes length from one snapshot to the next.

Fig. 11. Vertical plane passing through the field line containing the feature. The field line is shown by the dotted black line. A vertical plane with a curved base, containing the field line is shown for various quantities. The base of this plane on the τ500 = 1 surface is indicated by the dashed green curve in Fig. 9. Panel (a) shows the proxy (normalised κ) in the vertical plane. Panel (b) is the hydrogen number density at the n = 2 atomic level. Panels (c) and (d) are the vertical velocity and magnetic field components. Panel (e) is the temperature. Panel (f) shows the viscous and resistive heating rate. A time-series for the evolution of this feature and its MHD properties and heating in this plane is available online.

4.3. Heating and MHD properties

We aim to track the evolution of the feature over time, focusing on its MHD properties. As is outlined in Sect. 4.2, the feature is identified using the proxy in vertical cuts, and its movement is traced from one frame to the next. To study the feature in detail, we followed the field line it resides on and analysed itsproperties through a vertical plane intersecting this line, as is shown in Fig. 11.

The feature, as seen in the proxy, appears elongated and aligned with the magnetic field line. In this region, we observe an increase in the hydrogen populations at the n = 2 energy level, which is expected for enhanced absorption in Hα. The temperature in this area is characteristic of chromospheric temperatures (8 kK–12 kK). Additionally, the feature is located in an upflow region with a vertical velocity of 30–40 km/s.

In panel (d) of Fig. 11, the footpoint of the field line carrying the feature is rooted at a negative polarity patch. The footpoints of the overarching field lines discussed in Sect. 4.1 and Fig. 9 lie in the same magnetic patch. However, the positive flux emerges much closer – at a distance of about 2 Mm from the negative polarity footpoint. This is the same region highlighted by the red circle and red box in Fig. 9 and Fig. 10, respectively.

In the bottom right panel of Fig. 11 we plot the viscous and resistive heating rate. Heating is visible at the footpoints and partly in the region of the feature for this snapshot. When we track the heating through time in the linked movie, we can see strong heating at the onset of the feature, associated with flux emergence. There is also intermittent heating at the footpoints as the feature evolves. While the heating appears to show correlation with the feature’s development, this is not as apparent in the temperature distribution, which remains relatively homogeneous.

4.4. Lorentz forces associated with flux emergence

We explore the role of the Lorentz forces associated with the flux emergence in the region of our interest. A detailed description of the Lorentz force treatment in MURaM can be found in (Rempel 2017). Figure 12 and the associated animation show the evolution of the vertical component of the magnetic field in the plane of the traced field line. We also show the evolution of the Lorentz forces and the net force from the momentum equation in the vertical (z) direction. The emergence of positive polarity fields is visible close to the negative polarity footpoint of the feature. As the field emerges, it drives reconnection and the traced field line (dashed green line) is pushed upwards (see the movie associated with Fig. 12). Strong positive Lorentz forces act close to the L = 0 footpoint of the feature. As the emerged flux meets the traced field line, negative Lorentz forces temporarily hold down the field line from expanding. At t ∼ 90 s, the field line is suddenly seen to expand significantly and the Lorentz force then becomes positive. The net force is dominated by the Lorentz force component (F_L) above the solar surface.

Fig. 12. Snapshot showing the vertical field and forces after flux emergence with a highly kinked field. The vertical component of the magnetic field (Bz) and the Lorentz (FL, z), and net forces (Fnet, z) is shown. The traced field line containing the feature is shown by the dashed green line. The evolution of these forces and fields is available online.

4.5. Frame of reference of the field line

To analyse the temporal behaviour of the feature, we move to the reference frame of the traced 3D magnetic field line. This allows us to examine how the feature evolves both along the field line’s direction and perpendicular to it.

4.5.1. Geometry for our analysis

We define three mutually orthogonal directions: tangents (T), normals (N) and binormals (B). Such a decomposition has been done in previous works to study flows and forces along and perpendicular to the field lines (e.g. Leenaarts et al. 2015). The tangents follow the magnetic field. The direction was chosen from the negative footpoint (L = 0) directed towards the positive footpoint. This was done to ensure that the starting point for our analysis is the negative polarity patch where the feature is rooted. The orthogonal set of basis vectors for this frame of reference was constructed along each grid point, s = (x, y, z) of the traced field line. The set of tangent vectors is defined as

$$\begin{array}{c}\hfill \left\}\right\}\widehat{t}(s)=-\left\}\right\}\frac{B}{|B|}(\mathrm{s}).\end{array}$$(6)

We obtained the set of binormal vectors, $\widehat{\mathbf{b}}(\mathrm{s})$, with the following relation, keeping in mind the periodic boundary condition in the horizontal directions:

$$\begin{array}{c}\hfill \widehat{}\}b(\mathrm{s})=\pm \widehat{}\}t(\mathrm{s})\times \frac{\mathrm{d}\mathbf{t}}{\mathrm{ds}}.\end{array}$$(7)

The normals, $\widehat{\mathbf{n}}(\mathrm{s})$, are obtained with

$$\begin{array}{c}\hfill \widehat{}\}n(\mathrm{s})=\widehat{}\}b(\mathrm{s})\times \widehat{}\}t(\mathrm{s}).\end{array}$$(8)

We have a closed field line at all times, where consecutive tangent vectors are not parallel. This allows this basis construction to be used effectively to move to the frame of reference of the traced field line.

We chose the sign in Eq. (7) such that $\widehat{\mathbf{b}}.\widehat{\mathbf{x}}>0$. This results in $\widehat{\mathbf{n}}.\widehat{\mathbf{z}}>0$. It is to be noted that the field line is always roughly oriented along $\widehat{\mathbf{y}}$. Figure 13 shows a cartoon of the directions T-N-B along a traced field line in panel (a), with the feature in the proxy on the same field line shown in panel (b).

Fig. 13. The velocity components driving the feature and its visibility. Panel (a) shows the direction of the tangent (T), normal (N) and binormal (B) along the field line for this snapshot. Panel (b) shows the feature on the field line viewed in the proxy. Panel (c) shows the line-of-sight (vertical) velocity and outlines the feature with the magenta contour. Panels (d)-(f) show the flow components along a bunch of field lines traced through the feature.

4.5.2. Time-distance analysis

Using the set-up to analyse quantities within the field line’s reference frame, we constructed time-distance plots spanning the feature’s entire lifetime. Figure 14 illustrates the feature’s evolution in the proxy, along with various MHD properties along the field line. It is important to note that the length of the field line changes over time. This is also indicated in Figs. 14 and 15.

Fig. 14. Time-distance plots of various MHD quantities along the selected field line (dashed green line in Fig. 12). Panel (a) shows the evolution of the feature in the proxy. Panels (b), (c), and (d) show the gas pressure, density, and temperature, respectively. Panel (e) shows the divergence of flows, highlighting regions of compression (red)/expansion (blue). Panel (f) shows the viscous and resistive heating rate where the dashed green line shows the Alfvén speed in the region of the feature at ∼100 km/s. Contour lines highlight the location of the feature in panels (b) – (f). The length of the traced field lines change with time within the domain, so the data is padded with zeros on the right. This is shown with the light cyan colour.

The feature in our proxy exhibits an oscillatory behaviour over time, as is illustrated in panel (a) of Fig. 14 which shows a clear S-shaped evolution. This indicates that the feature moves back and forth along the field line during its lifetime, a motion that is also visible in the close-up movie linked to Fig. 10. The period of these oscillations is about 3 minutes. While turbulent motions near the surface displace the footpoint of this feature, the feature’s movement along the field line is notably larger and thus significant. Panels (b) and (c) in Fig. 14 reveal increased gas pressure and density within the region of the feature. However, panel (d) shows no corresponding increase in temperature. In panel (e), we see the flow field’s divergence, with the black-contoured feature region distinctly dominated by compression.

We also follow the heating of the feature in panel (f) of Fig. 14. This is the viscous and resistive heating rate along the field line. We find clear signatures of co-spatial and co-temporal heating at the birth of the feature. There is a propagating heating front which is marked by the dashed green line in panel (f) at a speed of 100 km/s. The Alfvén speed in that region also has a similar value.

We analysed the net forces and flows projected in the tangent (T), normal (N) and binormal (B) directions as described in Sect. 4.5.1. Figure 15 presents the component-wise analysis using vector fields – the flow field (panels a-c), pressure gradient force field (panels d-f), and Lorentz force field (panels g-i).

Fig. 15. Time-distance plots for the projection of flows and forces in the tangential (T), normal (N) and binormal (B) directions, respectively. Panels (a)-(c): flows projected in the T, N, B directions in the frame of reference of the selected field line. Panels (d)-(f) are the same for the net pressure gradient forces. Panels (g)-(i) show the Lorentz forces projected in the T, N, B directions. The dashed purple line in panel (b) marks the apparent velocity of propagation along the field line which is 50 km/s. The apparent velocity of the feature itself marked by the dashed magenta line is 38 km/s. Black contours highlight the location of the feature in all the panels. The padding on the right is indicated by the light orange colour.

The tangential velocity (V_T; panel a, Fig. 15) is predominantly positive within the feature (outlined by black contours). This means that the flow consistently moves away from the footpoint. The pressure gradient force along the tangential direction (panel d) is also mostly positive, indicating that pressure differences help drive the flow along the field line.

The normal velocity component (V_N; panel b) is strongly positive in the region of the feature, suggesting an upward and outward flow. This again is driven by the pressure gradient force (panel e), which acts to push the field line. The propagation front (marked by the dashed purple line in panel b) is about 50 km/s, which is faster than the apparent velocity of the feature at 38 km/s (dashed magenta line).

The transverse or binormal component of the flow field V_B exhibits a mix of positive and negative values, indicating a more complex lateral motion. The pressure gradient force (panel f) shows negative values at the feature’s tip, suggesting that this component could contribute to leftward lateral movement (as seen in the movie associated with Fig. 10).

The treatment of the Lorentz force is similar to that of the curved vertical slices described in Sect. 4.4. The forces acting normal to the field line (panel h, Fig. 15) show that strong positive forces precede the prominent appearance of the feature in the proxy. As the feature forms prominently, we find some negative forces associated with the emerged flux interacting with the field line. At about t ∼ 90 s, the force becomes positive again, aiding in the field line expansion (also see Sect. 4.4). The binormal component (panel i, Fig. 15) of the Lorentz force shows some domination of negative forces.

5. Discussion and conclusions

We have developed a computationally efficient proxy for identifying Doppler-shifted Hα fine structure in MURaM-ChE simulations. The features identified by the proxy were validated with MULTI3D synthesis. By approximating key aspects of the NLTE RT problem, the proxy bypasses the need for full 3D NLTE synthesis. This reduces the computational costs from several hundred thousand CPU hours per snapshot to a few CPU seconds. The proxy allows us to capture the morphology and position of chromospheric features appearing in the wings of the spectral line, effectively covering Doppler shifts from approximately 77 km/s down to 17 km/s in either wing. The proxy is physically motivated in terms of capturing the interaction of the plasma with the radiation to identify features in the solar atmosphere visible in the Hα wings. Although the proxy focuses on features in the wings of Hα and does not model the line core, it has shown reliable accuracy in identifying TPs and TNs within its scope. Because the proxy is very fast, it allows us to identify many such features, so that the statistical properties of these features can be studied, which will be the topic of a follow-up paper.

The proxy assumes a Gaussian profile for the spectral dependence of the escape probability. To accurately model the features in the far wings where pressure broadening becomes important, a Voigt profile is potentially more fitting for the construction of the proxy. The average escape probability profile from the MULTI3D synthesis matches the Gaussian profile closely. Thus the approximation of using a Doppler profile to model the proxy works well. We also found that the maximum difference of the fits to the average exp(−τ) profile of the synthetic Hα line (Fig. 4) with a Voigt profile and a Gaussian profile is less than 10% of the maximum.

The RBEs and RREs identified using the proxy are similar in morphology and evolution to those observed. These are fine needle-like projections seen on the solar disc and appear to be jets. We also identify some fibrillar structures which are broader and more smudged. All these features appear to be rooted close to the network patch in our simulations, similar to observations in Hα and Ca II H & K.

The example feature that was investigated is an RBE at a Doppler shift of 37 km/s. We infer the morphology of the feature from two different vantage points – (1) the on-disc or top view and (2) the vertical cuts. In essence, we use the absorption coefficient directly (in vertical cuts) or to construct the escape probability (in the top view) to understand the feature in 3D. The on-disc map shows a nearly straight feature evolving on the solar disc. This gives an impression that it is a straight, jet-like feature. The proxy in the vertical cuts shows that we are actually looking near the top of a closed loop structure. Given that the feature is pushed along the field line, we observe a jet-like behaviour. The loop itself also moves upwards and expands carrying the feature along with it. However, it is different from a traditional type-II spicule, where we would expect a straight feature to shoot up in a linear trajectory. Whether the feature does appear as a type II spicule from a different viewing angle is a matter of future investigation. These aspects further highlight that studying the 3D morphology of these features is important for a comprehensive understanding of their behaviour. It also highlights the need to determine if there are fundamental differences between spicules and RBEs/RREs.

The feature has a maximum length of 3.4 Mm during its lifetime spanning 4 minutes. It also exhibits lateral displacement (∼400 km). The height of formation for this feature is 2 − 4 Mm, corresponding to heights within the chromosphere. The proxy identifies the correct height of formation, consistent with the result from MULTI3D. Our time-distance analysis shows that in the region of the feature the temperature ranges from 8 kK to 10 kK on average, corresponding to the chromosphere and low transition region. The density ranges from 10⁻¹³ g/cm³ to 10⁻¹² g/cm³. The pressure in the region of the feature is approximately 0.1 dyne/cm² to 1 dyne/cm². While the density, pressure and hydrogen (n = 2) densities are enhanced in the region of the feature, the temperature does not show this local enhancement. This concurs to the fact that the Hα opacity is only weakly sensitive to the temperature (Leenaarts et al. 2012).

The formation of the RBE is driven by flux emergence and magnetic reconnection. The flux emergence process pushes the plasma upwards. Reconnection follows and a convergence of flows in the region of the feature is clearly visible (panel f in Fig. 14). The signature of a bi-directional flow is visible at the site of reconnection (boundary of the feature in panel a, Fig. 15). The primary force behind the field-aligned flows acting on the feature is the pressure gradient. The reconnection outflow (panel a, Fig. 15) hits previously slow-moving gas and leads to compression, resulting in the local enhancement of the plasma density and pressure. Pressure gradient forces build up and this drives a flow along the field line carrying the feature. This makes our example feature behave like a jet. The pressure gradient forces also lead to strong flows acting on the field line itself pushing it upwards, causing it to expand. We also find strong Lorentz forces associated with the flux emergence process. The transverse component of the forces may be responsible for the lateral motion (∼400 km) seen in the feature.

The time-distance analysis of the flow component perpendicular (outwards) to the field line reveals actual mass flows at speeds of 37 − 39 km/s associated with the feature. In addition, we find a propagation front at 50 km/s. This front may correspond to faster features appearing in the same region as our target feature. It could also indicate flows which are not associated with mass flows, but rather wave modes propagating at sound speeds.

It is to be noted that we have identified features based on the line-of-sight, that is, vertical component of the velocity. Hence we detect the feature only in upflow regions of ∼37 km/s (see Fig. 13). However, the flow speeds in Fig. 13 along the field line (v_T) are much slower (∼7 km/s). Although there are strong upflows, the region occupied by the upflows and the plasma does not move significantly. Hence the jet appears relatively stationary.

We can see that bidirectional flows occur throughout the feature’s lifetime (panel a, Fig. 15), and the field remains kinked (movie associated with Fig. 12). Thus, the reconnection appears to be powering the feature for the most part of its lifetime.

An intriguing propagation front that we observe in our analysis comes from the propagation of the viscous and resistive heating that travels at Alfvénic speeds of ∼100 km/s. The strong viscous and resistive heating close to the first appearance of the feature exhibits this heating front. The footpoints of the field line carrying the feature are also heated intermittently. While this heating could possibly be unrelated to the feature itself, it shows a good correlation with the flux emergence event driving the feature formation (see the movie linked to Fig. 11). The heating at the footpoints show no associated feature formation. This could partly be because there are no associated features present, or because we do not identify them at the Doppler velocity used for this study. There have been conjectures on rapidly propagating heating fronts associated with spicules that extend into the corona (De Pontieu et al. 2011; Bose et al. 2021). What we observe could be one of these propagating heating fronts responsible for transporting the energy from the chromosphere into the corona, thus being one of the viable mechanisms of coronal heating. However, we do not find a multi-thermal nature associated with the feature unlike in Bose et al. (2021), although they study a downflowing feature.

We do not find the heating to manifest itself in the temperature. Radiative losses and ionisation processes may be responsible for this. The energy could also be used to increase the kinetic energy of the plasma, driving mass motions rather than a temperature increase. We surmise the latter to be the case.

Rutten (2017) argues that the high opacity of Hα fibrils is caused by hot precursors called propagating heating events or PHEs. The fibrils themselves are contrails of such fast-propagating heating events. RBEs also fall in the family of features that PHEs give rise to. Although we do not trace back to the conditions before the formation of the feature, we find the strong heating front only after the flux has emerged and interacts with the field line. This happens co-temporally with the appearance of the studied feature. We argue that magnetic reconnection events lead to viscous dissipation generating the heating in this feature. The formation of the feature happens simultaneously with – not following the heating event.

We also find an oscillatory behaviour: the feature moves back and forth along the magnetic field line carrying it. The period of these oscillations is about 3 minutes. This is visible in the S-shaped evolution of the feature in panel a, Fig. 14 and the movie associated with Fig. 10. We find a lateral displacement of a few hundred kilometres of the feature in a leftward direction in the on-disc maps. This is consistent with observations of RBEs (Sekse et al. 2013a,b). Transverse displacements of RBEs throughout the entire length of the feature have also been observed (Rouppe van der Voort et al. 2009; Kuridze et al. 2015). In our case, it appears to be the displacement of the entire flux tube at all heights as speculated by Kuridze et al. (2015). In a numerical study of fibrils by Druett et al. (2022), they argue that the gas pressure gradient plays a dominant role in such lateral displacements. Our findings are in sync with this. We find that it is a combination of both the gas pressure gradient and the Lorentz forces acting on the flux tube that moves it in a general leftward (-x) direction in the on-disc view.

We successfully observed many features similar to RBEs and RREs without the use of ambipolar diffusion in the simulations. Bose et al. (2021) found ambipolar diffusion (Martínez-Sykora et al. 2017) to play a major role in heating with their Bifrost (Gudiksen et al. 2011) simulations. It could be that the exclusion of ambipolar diffusion led to an insignificant impact on our temperatures. However, the formation of jet-like features in simulations done with the MURaM-ChE code, with properties similar to observed RBEs, does not necessitate ambipolar diffusion to be incorporated.

The RBEs and RREs analysed by Danilovic et al. (2023) exhibit several observed properties – such as spatial distribution, lateral motions, lengths, and lifetimes – that closely align with our findings. Additionally, their reported line-of-sight velocities (30 − 40 km/s) and formation heights (2 − 4 Mm) are consistent with the values we observe.

The proxy can be developed further to work at the solar limb to study spicules. This will be presented in a future work. Devising a proxy for the line core of Hα to study dynamic fibrils is not straightforward. As was elucidated by Leenaarts et al. (2012), the line core has nuances which are difficult to model. Modelling the entire Hα spectral line at a low computational cost remains a challenge.

Data availability

Movies associated with Figs. 5, 7–12 are available at https://www.aanda.org

Acknowledgments

We thank the anonymous referee for valuable comments that helped improve the paper. The authors also thank Yajie Chen for useful discussions. This work was carried out in the framework of the International Max Planck Research School (IMPRS) at the Technical University of Braunschweig. We are grateful for the computational resources provided by the Cobra and Raven supercomputer systems of the Max Planck Computing and Data Facility (MPCDF) in Garching, Germany. This research has received financial support from the European Union’s Horizon 2020 research and innovation program under grant agreement No. 824135 (SOLARNET) and through the European Research Council (ERC) No. 101097844 (WINSUN). This work was supported by the Deutsches Zentrum für Luft und Raumfahrt (DLR; German Aerospace Center) by grant DLR-FKZ 50OU2201. This project has received funding from the Swedish Research Council (2021-05613) and the Swedish National Space Agency (2021-00116). We acknowledge resources provided by the National Academic Infrastructure for Supercomputing in Sweden (projects NAISS 2023/1-15 and NAISS 2024/1-14) at the PDC Centre for High Performance Computing (PDC-HPC) at the Royal Institute of Technology.

References

All Tables

All Figures

Fig. 1. Overview of the magnetic field and synthetic Hα intensity in a snapshot of the enhanced network simulation. In panels (a) and (b) we show the vertical component of the magnetic field (Bz) at the surface (τ500 = 1) and in the chromosphere (z = 1 Mm), respectively. Panels (c) and (d) are synthetic Hα intensities in the blue (37 km/s) wing and red (−37 km/s) wing, respectively. vD in the title of panels (c) and (d) denotes the Doppler shift.

In the text

Fig. 6. 2D histogram plots showing the classification of pixels. Panels (a), (b): 2D histograms of the proxy vs. the synthetic Hα intensity at Doppler velocity = 37 km/s (a) and −37 km/s (b). Panels (c), (d): The row normalised distribution of the same pixels corresponding to panels (a) and (b), respectively. We mark the dynamic threshold for the synthetic intensity with vertical dashed black lines. The static threshold for the proxy for feature identification is marked by the horizontal dashed black lines. The four pixel categories: true positive (TP), true negative (TN), false positive (FP), and false negative (FN) are also indicated on each of the panels.

In the text

	Fig. 2. Spatially averaged profile of the Hα proxy, modelled as an escape probability (orange). The same quantity from the MULTI3D MURaM (-ChE) synthesis for comparison (purple).
In the text

	Fig. 3. Spatially averaged Hα line profile from the FTS Atlas of Neckel & Labs (1984) (green) compared with the Hα profile synthesised by applying MULTI3D to the simulated MURaM-ChE atmosphere (purple). The dashed black line marks the centre of the line core. We have shifted the synthetic profile (MULTI3D MURaM) by 0.035 Å(≈1.59 km/s) towards the blue wind to align the line cores for comparison.
In the text

	Fig. 4. Average escape probability profile from the MULTI3D synthesis on MURaM-ChE (black circles) fitted with a Voigt profile (maroon) and a Gaussian profile (teal).
In the text

Fig. 5. Features in the blue wing at a Doppler velocity of 37 km/s. The leftmost panel shows the average line profile for the synthesised Hα intensity (using MULTI3D) with wavelength. The vertical blue line marks the wavelength at which the Doppler velocity is 37 km/s with respect to the line centre (Δλ = 0). The middle panel shows the synthetic Hα intensity at this wavelength. The third panel shows the proxy at the same wavelength. A scan through the line is available online.

In the text

Fig. 7. Goodness of the proxy at Doppler velocity = 37 km/s. Panel (a) gives a spatial distribution of (in)correctly identified pixels. The pixels are colour-coded to show true positives (TPs: red), true negatives (TNs: white), false positives (FPs: green) and false negatives (FNs: blue). Panel (b) represents the percentage of (in)correctly identified pixels as a function of the Doppler velocity (vD) in the blue wing of the line. The dashed black line marks the position in the line wing (vD = 37 km/s) at which panel (a) is obtained. A scan through the blue wing of the line exhibiting the goodness of the proxy is available online.

In the text

Fig. 8. Overview of one blueshifted feature (vD = 37 km/s) detected by the proxy. Panel (a) shows the proxy map with the feature of interest highlighted by the red box. The green line passing through the feature marks the base for the vertical cut through the feature (hereby x′ axis). The 0 of x′ corresponds to (x, y) = (8.6, 0) Mm in the top view (in panel a). Panels (b) and (c) are a zoomed-in view of the feature showing the synthetic intensity in Hα and formation height, respectively (obtained with the MULTI3D synthesis). In panel (d) we take a vertical cut through the feature in the proxy (normalised κ). Panels (e) and (f) show the same vertical cut for the line-of-sight (i.e. vertical) component of the velocity and magnetic field, respectively, with the feature over-plotted in black. A time-lapse of the proxy map in the top view and vertical cuts of various MHD quantities is available online.

In the text

Fig. 9. 3D rendering of magnetic field lines rooted at the footpoint of the feature at t = 0 s, the first appearance of the feature. The line-of-sight component of the magnetic field is shown as a surface map at the τ500 = 1 surface. The 3D field lines are coloured by the vertical (Bz) component. The site of magnetic reconnection lies across the main footpoint and is marked by the red circle. The dashed green line shows the xy projection of the 3D field line carrying the feature in this snapshot. The magnetic field lines are clipped for better visualisation. A time-series showing the rearrangement of magnetic field lines at the time of feature formation is available online.

In the text

Fig. 10. Zoom-in on the evolution of the line-of-sight magnetic field at various τ levels in the region of the studied feature. Left: Evolution of the vertical component of the magnetic field, Bz at the photosphere (τ500 = 1). The feature from the proxy is overplotted in purple. The dashed red box indicates the region of flux cancellation. Centre: Evolution of Bz in the middle photopshere (τ500 = 0.01). The grey scale for the magnetic field (Bz) ranges from −500 G to 500 G at the τ500 = 1 surface and from −100 G to 100 G at the τ500 = 0.01 surface. The regions of enhanced current density (75% above the mean) are plotted with the cyan contours in both the left and central panels. Right: Zoom-in on the evolution of the feature in the proxy as viewed on disc. A movie spanning the lifetime of the feature is available online.

In the text

Fig. 11. Vertical plane passing through the field line containing the feature. The field line is shown by the dotted black line. A vertical plane with a curved base, containing the field line is shown for various quantities. The base of this plane on the τ500 = 1 surface is indicated by the dashed green curve in Fig. 9. Panel (a) shows the proxy (normalised κ) in the vertical plane. Panel (b) is the hydrogen number density at the n = 2 atomic level. Panels (c) and (d) are the vertical velocity and magnetic field components. Panel (e) is the temperature. Panel (f) shows the viscous and resistive heating rate. A time-series for the evolution of this feature and its MHD properties and heating in this plane is available online.

In the text

	Fig. 12. Snapshot showing the vertical field and forces after flux emergence with a highly kinked field. The vertical component of the magnetic field (Bz) and the Lorentz (FL, z), and net forces (Fnet, z) is shown. The traced field line containing the feature is shown by the dashed green line. The evolution of these forces and fields is available online.
In the text

Fig. 13. The velocity components driving the feature and its visibility. Panel (a) shows the direction of the tangent (T), normal (N) and binormal (B) along the field line for this snapshot. Panel (b) shows the feature on the field line viewed in the proxy. Panel (c) shows the line-of-sight (vertical) velocity and outlines the feature with the magenta contour. Panels (d)-(f) show the flow components along a bunch of field lines traced through the feature.

In the text

Fig. 14. Time-distance plots of various MHD quantities along the selected field line (dashed green line in Fig. 12). Panel (a) shows the evolution of the feature in the proxy. Panels (b), (c), and (d) show the gas pressure, density, and temperature, respectively. Panel (e) shows the divergence of flows, highlighting regions of compression (red)/expansion (blue). Panel (f) shows the viscous and resistive heating rate where the dashed green line shows the Alfvén speed in the region of the feature at ∼100 km/s. Contour lines highlight the location of the feature in panels (b) – (f). The length of the traced field lines change with time within the domain, so the data is padded with zeros on the right. This is shown with the light cyan colour.

In the text

Fig. 15. Time-distance plots for the projection of flows and forces in the tangential (T), normal (N) and binormal (B) directions, respectively. Panels (a)-(c): flows projected in the T, N, B directions in the frame of reference of the selected field line. Panels (d)-(f) are the same for the net pressure gradient forces. Panels (g)-(i) show the Lorentz forces projected in the T, N, B directions. The dashed purple line in panel (b) marks the apparent velocity of propagation along the field line which is 50 km/s. The apparent velocity of the feature itself marked by the dashed magenta line is 38 km/s. Black contours highlight the location of the feature in all the panels. The padding on the right is indicated by the light orange colour.

In the text