Deciphering Solar Coronal Heating: Energizing Small-scale Loops through Surface Convection

The initial configuration consists of a rotationally symmetric 3D nullpoint magnetic topology (see Figure A1 in Appendix A, where details concerning the initial condition and the numerical code are also given). The inner spine of the structure is rooted in a region with a parasitic positive polarity, while the fan structure is anchored in the surrounding field of the opposite sign. A general outlook of the subsequent evolution of the experiment can be gained through Figure 1 and its associated animation. The initial topology structure is rapidly disrupted by the stochastic granulation in the photosphere leading to a magnetic network pattern of the predominant polarity at the surface; the overall magnetic field topology above is stressed and magnetic reconnection is triggered at the nullpoint. Simultaneously, the magnetic loops underneath the fan surface experience intense heating, giving rise to a CBP with conspicuously enhanced emissivity in EUV, X-rays, and UV, as evidenced in Figure 1 through the synthetic observables of panels (f), (g), (h), and (i) (see Appendix C for the calculation details). The diameter of the CBP, estimated from the EUV response at its base, is approximately 20 Mm, which is within the typical range reported in observations (5–40 Mm; see Madjarska 2019, and references therein).

The CBP nullpoint undergoes continuous perturbations, as illustrated in panel (b) through the time-evolution of its coordinates (x_np, y_np, z_np), which can be shifted by up to a few megameters, similarly to the nullpoint displacements inferred from observations (Galsgaard et al. 2017). Around t = 111.67 minutes, the time of the snapshots in Figure 1, the current sheet around the nullpoint experiences a major disturbance, concerning both its location and shape. The reason for this perturbation is the approach of an elongated cool and dense structure to the reconnection site visible in panels (c) and (d). This structure is shown to be a Hα chromospheric fibril in Section 2.3. Throughout the evolution of the CBP, other chromospheric fibrils and spicule-like structures are found to develop in the lower atmosphere (see animation at t = 47, 107, 153, 221 minutes), reaching the reconnection site and affecting the CBP dynamics as well as its brightness. Such features could be the result of the interaction of chromospheric waves and shocks with the CBP magnetic structure. A detailed analysis will be covered in a follow-up paper.

Rapid upflows can be seen associated with the CBP all along the experiment, as visualized through the blue–purple volume rendering in panel (a), which shows velocities from 45 to 75 km s⁻¹, and through the vertical velocity u_z in panel (e). These upflows are mainly concentrated around the spines of the CBP; those along the outer spine are basically reconnection outflows emanating from the nullpoint, whereas the ones around the inner spine are probably due to the channeling of magnetosonic waves. The flows arising from the CBP nullpoint, like the one at the time of the figure, are of particular interest for the corona. These hot (1–2 MK) collimated ejections have low contrast in density compared to the surrounding plasma, and therefore they are not (or barely) distinguishable in the synthetic EUV images (see the discussion in Section 3).

To address the fundamental question about the sustained CBP heating and its location, we have first analyzed the heating per unit mass due to the Joule and viscous terms. Panels (a) and (b) of Figure 2 contain time-averaged results over approximately 4.5 hr of the simulation at three different heights, z = [8, 5, 2] Mm, thus encompassing the regions in close proximity to the CBP nullpoint downwards to the low atmosphere. Both the Joule heating (left column) and viscous heating (right column) primarily occur at lower heights, at a few megameters above the solar surface, concentrated around the inner spine, in the region dominated by the strong parasitic positive polarity. We found that in these regions the electric current vector subtends a small angle with respect to the concentrated nearly vertical magnetic field, indicating that the magnetic configuration is close to force-free there. Heating is also found, albeit to a lesser extent, in the magnetic patches outside of the fan surface in the surroundings of the CBP at z = 2 Mm (bottom panels), as well as in the region near the nullpoint at z = 8 Mm (top panels). Thus, our findings indicate that CBP loops are primarily heated in the low atmosphere, with a secondary contribution at coronal heights near the reconnection site. Furthermore, there is a remarkable spatial correlation between the Joule and viscous heating across various heights. This suggests that regions exhibiting the highest electric current intensity and velocity gradient typically coincide. As a consequence, our results concerning the location of the CBP heating could be general, regardless of which mechanism (Joule or viscous) is dominant, even if in our experiment the Joule heating term is approximately 1 order of magnitude greater than the viscous one (see the discussion in Section 3).

Zoom In Zoom Out Reset image size

Another important question for CBPs, and for the corona in general, concerns whether the heating occurs on short or long timescales (Reale 2014). To pursue this question, panels (c) and (d) of Figure 2 illustrate, as a function of time, the horizontal average within the CBP fan surface of the Joule (black) and viscous (orange) heating per unit mass at z = [2, 5] Mm, respectively. Both heating mechanisms seem to be continuously at work, showing fluctuations with different time scales. For instance, the pronounced peak observed mainly in the viscous heating profile at z = 5 Mm at t ∈ [110, 120] minutes corresponds to the Hα fibril eruption explained in the previous section. In general, there is a clear time correlation between the Joule and viscous dissipation mechanisms. In fact, the Pearson correlation coefficient r between them is greater than 0.8 for the different heights we have analyzed.

As a final point, we consider the energy injection through the Poynting flux. Panel (e) of Figure 2 contains the mean vertical Poynting flux computed within the CBP fan surface and averaged in time, $\overline{\langle {S}_{z}\rangle }$, as a function of height (blue curve). This plot strongly suggests that the heating of the plasma within the CBP domain is due to the deposition of the upward Poynting flux, which, itself, results from the stochastic photospheric motions stressing the CBP magnetic field. This implies that no major organized photospheric flows (leading, e.g., to spine dragging or convergence of magnetic polarities) are necessary to maintain the heating of the CBPs. Importantly, the average $\overline{\langle {S}_{z}\rangle }$ within the CBP is larger at all heights than the corresponding quantity in the surroundings (red curve), and the same applies to the corresponding gradient in height (in absolute value). The greater vertical Poynting flux obtained within the CBP is due to the strong magnetic concentration in the parasitic polarity. The Poynting flux values and the height gradient are in agreement with those indicated by other radiative-MHD numerical experiments as necessary to sustain heating in the quiet Sun corona (Rempel 2017; Chen et al. 2022). In addition, these values seem to fulfill the energy requirements for a hot corona inferred from observations (e.g., 0.8 × 10⁶ erg cm⁻² s⁻¹ for coronal holes; Withbroe & Noyes 1977). Moreover, computing the surface integral of the vertical Poynting flux within the fan surface domain over two horizontal planes in the corona (e.g., z = 2 Mm and z = 5 Mm), we obtain a magnetic energy deposition in the CBP of ≈7 × 10²⁴ erg s⁻¹, which is more than enough to explain the rough estimate of the CBP energy losses obtained through observations (10²³–10²⁴ erg s⁻¹; Habbal & Withbroe 1981). The loss of electromagnetic energy through the fan surface is most probably small by comparison, at any rate when averaging in time.

To reinforce the previous results, we present here SDO and SST observations (see Appendix B) of a CBP that occurred on 2022 July 1 located at heliocentric coordinates (x, y) = (177'', 147''). The million-K coronal response of the CBP, shown in the SDO/Atmospheric Imaging Assembly (AIA) images (panels (a) and (b) of Figure 3), has a compact and bright rowel-like shape that covers a horizontal domain of ∼40'' × 40''. The photospheric line-of-sight magnetic field image from the SST/CRisp Imaging Spectropolarimeter (CRISP) in panel (d) shows a parasitic negative polarity at the center of the CBP, surrounded by a dominant opposite polarity magnetic region, a configuration that can naturally lead to a fan-spine configuration with a nullpoint at coronal heights (Zhang et al. 2012; Mou et al. 2016; Galsgaard et al. 2017; Madjarska et al. 2021; Cheng et al. 2023). The chromospheric Hα core image from SST/CRISP in panel (c) indicates that the CBP contains a bunch of dark and long Hα fibrils that connect opposite polarities, analogous to fibrils reported in lower resolution observations (Madjarska et al. 2021), and there are strong brightenings mainly located at the base of the fibrils, possibly indicative of chromospheric heating. The spectral profiles at different locations of the fibrils and brightenings are shown in panel (e).

Zoom In Zoom Out Reset image size

In order to prove that our model satisfactorily reproduces the observations, Figure 4 and its associated animation contain the forward-modeling results from the simulations. Panel (a) shows that we indeed obtain a rowel-like structure composed of hot loops with enhanced EUV response similar to the SDO/AIA 193 Å observations, although our CBP is smaller than the observed one. Given its current interest, in panel (b) we also add the EUV response as would be observed by the Solar Orbiter's Extreme Ultraviolet Imager (EUI) instrument, with its higher spatial resolution. The comparison to observations is also successful for the chromospheric counterpart, as seen in panel (c). Our CBP model is realistic enough to reproduce a pattern of elongated dark fibrils, as observed in the Hα images, and the existence of strong brightenings at the base of these fibrils, surely a manifestation of the location of the heating in the lower atmosphere, as discussed in the previous section. The agreement is further illustrated by comparing panels (e) of Figures 3 and 4. The observational and synthetic spectral profiles show a similar behavior regarding the Hα line core intensity at different locations of the fibrils and brightenings, namely, the fibrils are characterized by a less bright Hα core than the average with a small Doppler shift, while the brightenings show an enhanced Hα core. The smaller width in the synthetic profiles may be due to the modest numerical resolution not capturing the highly dynamic complexity of the chromosphere (Hansteen et al. 2023). At the photosphere, the vertical magnetic field of the model (panel (d) of Figure 4) shows that the elongated Hα fibrils connect the central parasitic polarity with the strong, opposite-polarity magnetic patches surrounding the CBP, as shown in our observations (panel (d) of Figure 3) as well as in other CBP observations (Madjarska et al. 2021).

The experiment was carried out using the radiation-MHD Bifrost code (Gudiksen et al. 2011). This code includes, in a self-consistent manner, radiative transfer from the photosphere to the corona; approximate recipes for the main radiative losses in the chromosphere due to neutral hydrogen, singly ionized calcium, and singly ionized magnesium; field-aligned thermal conduction, which tends to an isotropic conductivity prescription for regions with a weak magnetic field, such as nullpoints; optically thin losses; and an equation of state considering the 16 most relevant atomic elements for the solar atmosphere either because of their high abundance or because of their importance as electron donors at low temperatures. Explicit resistivity and viscosity are included in Bifrost through hyper-diffusion terms (Gudiksen et al. 2011), following the idea originally described by Nordlund & Galsgaard (1995). A comprehensive description of the resistive and viscous hyper-diffusion operator and the associated free parameters can be found in the papers by Færder et al. (2023a, 2023b), including a comparative analysis of the hyper-diffusive resistivity versus other resistivity models.

The initial condition was built using a weakly magnetized (0.2 G at the corona) preexisting 3D numerical experiment that had reached a statistically stationary equilibrium. In the vertical direction, the numerical box covers a region from the uppermost layers of the solar interior to the corona, that is, −2.9 Mm ≤ z ≤ 32.3 Mm, with z = 0 corresponding to the solar surface (more precisely the horizontal surface, where τ_{500 nm} = 1). The horizontal extent is 32 × 32.0 Mm², with both x and y starting at the origin of coordinates (0, 0). The domain is resolved with 512 × 512 × 512 cells, employing a uniform grid in the horizontal directions with spacing Δx = Δy = 62.5 km, and a nonuniform mesh for z, using Δz = 50 km from the interior up to 12 Mm in the corona, and then progressively decreasing the resolution down to Δz = 150 km at the top of the corona. The boundary conditions are periodic in the horizontal direction. In the vertical direction, at the top, characteristic boundary conditions are applied; at the bottom, constant entropy is set for the incoming plasma to keep solar-like convection, while the rest of the variables are extrapolated. Panel (a) of Figure A1 shows the horizontal averages of the temperature T and mass density ρ of the background stratification.

On top of the relaxed snapshot mentioned above, we added a potential magnetic nullpoint distribution as illustrated in panel (b) of Figure A1. The potential configuration was calculated from a prescribed axisymmetric distribution at the bottom boundary containing a circular positive polarity (the parasitic polarity) embedded in a negative background. This imposed configuration is such that (a) there is a nullpoint at (x_np, y_np, z_np) = (16, 16, 8) Mm; (b) at large heights, the magnetic field becomes asymptotically vertical and uniform with B_z → − 10 G, thus mimicking a coronal hole structure; and (c) at z = 0, the total positive flux is Φ⁺ = 4.9 × 10¹⁹ Mx, the maximum positive vertical field strength is B_z = 141.8 G, the parasitic polarity covers a circular patch of radius 6.3 Mm, and the intersection of the fan surface with the horizontal plane has a radius of 12.4 Mm.

Zoom In Zoom Out Reset image size

We have used observations from the SST (Scharmer et al. 2003) obtained on 2022 July 1 from 08:04:57 to 08:13:03 UT centered at a CBP located at heliocentric coordinates (x, y) = (176'', 135''). To cover the whole CBP and surroundings within the field of view (FOV) of the SST instrumentation, we have created a mosaic constructed from 3 × 3 telescope pointings. The SST data set contains the following spectral scans from the CRisp Imaging Spectropolarimeter (CRISP; Scharmer et al. 2008): the Fe i 6173 Å line in spectropolarimetric mode sampled in 14 positions between −320 and +680 mÅ from the line core, and Hα 6563 Å in spectroscopic mode sampled in 31 positions ranging from −1500 to +1500 mÅ from the line core. Each of the CRISP scans of the mosaic has an FOV of ∼58'' × 58'' and a spatial sampling of 0057 per pixel. The CRISP data were processed using the SSTRED data reduction pipeline (de la Cruz Rodríguez et al. 2015; Löfdahl et al. 2021), which includes Multi-Object Multi-Frame Blind Deconvolution (van Noort et al. 2005) image restoration. Line-of-sight (LOS) magnetograms were obtained by Milne–Eddington inversions of the Fe i 6173 Å Stokes profiles (de la Cruz Rodríguez 2019).

The coronal response of the targeted CBP is obtained from the AIA (Lemen et al. 2012) onboard SDO (Pesnell et al. 2012) around 08:08:52 UT: the time at which the center of the SST mosaic was taken. The photospheric magnetic field is obtained from the Helioseismic and Magnetic Imager (HMI; Scherrer et al. 2012). Regarding the alignment, cross-alignment is carried out in all the AIA channels to the HMI continuum; then, the SDO images are manually aligned with the SST ones.

For the forward modeling in the EUV and UV, the expressions of optically thin radiative transfer are used, assuming, additionally, coronal abundances (Feldman 1992) and ionization equilibrium. The emissivity corresponding to any given transition between electron configuration levels i and j of a given ionic species of an atomic element (e.g., Fe ix, or Si iv) is given by

$$\begin{eqnarray}&&\epsilon ={n}_{{\rm{H}}}\ {n}_{e}\ G(T,{n}_{e},{\nu }_{{ij}}),\end{eqnarray} \tag{ C1 }$$

where n_e and n_H are the electron and hydrogen number densities, respectively; ν_ij is the radiation frequency of the emitted photon; and G(T, n_e , ν_ij ) is the corresponding gain function (or contribution function), which already contains the abundance of the atomic element. The emitted intensity is calculated by integrating Equation (C1) along the LOS in the numerical box.

For the SDO/AIA coronal synthetic images, we have used CHIANTI (Del Zanna et al. 2021) from the Solar Soft package, calculating the gain functions for the atomic transitions in the relevant spectral range of the different AIA filters and weighing with the corresponding effective Area Function. The result is then transformed into SDO/AIA count numbers. To account for the possible obscuration effects from cool and dense features in the EUV images, we have included absorption by neutral hydrogen, neutral helium, and singly ionized helium in the integration of the emissivity along the LOS adding an absorption factor e^−τ following the recipes and procedures in the literature (Anzer & Heinzel 2005; De Pontieu et al. 2009).

For the synthetic 174 Å images of the Extreme Ultraviolet Imager of the High Resolution Imager (EUI-HRI; Rochus et al. 2020) on Solar Orbiter (SoLO; Müller et al. 2020), we have employed a contribution function privately provided by Dr. Frédéric Auchère, a member of the Solar Orbiter team.

For the synthetic IRIS (De Pontieu et al. 2014) transition region counterpart, the contribution function of the Si iv 1393.755 Å line is obtained through a simple procedure (Chianti Solar Soft routine ch_synthetic.pro with the flag /goft), multiplying the output by the silicon abundance relative to hydrogen. This procedure and the conversion to IRIS count numbers are detailed in a recent paper (Nóbrega-Siverio et al. 2017).

After obtaining the intensity images, we degraded the results to match the spatial resolution of the instruments, namely, 15 for SDO/AIA (Lemen et al. 2012); 033 for IRIS (De Pontieu et al. 2014); and 200 km for EUI-HRI at perihelion (Rochus et al. 2020).

The X-ray synthetic images have been computed to mimic observations taken by the X-ray telescope (XRT; Golub et al. 2007) onboard Hinode (Kosugi et al. 2007). In particular, we have combined the make_xrt_wave_resp.pro and make_xrt_temp_resp.pro routines from Solar Soft to obtain the spectral and temperature response for the Al-poly channel of the telescope. The result is correspondingly multiplied by the squared electron number density and then degraded to the 20 instrument resolution (Golub et al. 2007).

To obtain Hα synthetic images, we have used the 3D non–local thermodynamic equilibrium radiative transfer code Multi3D (Leenaarts & Carlsson 2009). Multi3D evaluates the radiation field in full 3D taking the horizontal structure of 3D model atmospheres into account. To reduce the 3D radiative transfer computation cost, we have halved the horizontal resolution of the model atmosphere. We further optimized the depth scaling of the atmospheric quantities to improve the convergence time and to reduce the effect of an insufficiently sampled optical depth scale onto the synthesized Hα line profiles. Hα was synthesized from a five-level plus continuum hydrogen model atom assuming Doppler line profiles for the Lyα and Lyβ lines, which affect the level populations of the Hα transition (Leenaarts et al. 2012). We have reduced the spatial resolution of the synthetic images to the SST/CRISP instrumental values, namely, 013 at 6301 Å (Scharmer et al. 2008).