© M. V. Sieyra et al. 2022

1. Introduction

The study of magnetohydrodynamic (MHD) waves has two major applications: coronal seismology and coronal heating. In regards to the former, the properties of MHD waves depend on the properties of the plasma environment where they propagate, and therefore it is possible to use waves to probe macroscopic parameters of the plasma in the vicinity of the wave-propagation region (e.g., De Moortel & Nakariakov 2012). In relation to coronal heating, waves could carry the energy through different layers of the solar atmosphere and deposit it in the corona via different dissipative processes (see Van Doorsselaere et al. 2020), contributing to the maintenance of its high temperature.

In this study, we focus on propagating slow magnetoacoustic waves. These compressive waves are in general related to periodic intensity perturbations mostly observed in extreme-ultraviolet (EUV) or X-ray emission moving along the field-aligned plasma. They have been detected in several coronal structures such as polar plumes (e.g., Ofman & Davila 1997; DeForest & Gurman 1998) and active region fan loops (e.g., Berghmans & Clette 1999; De Moortel et al. 2002; Marsh & Walsh 2006). However, there is an alternative interpretation of the observed intensity disturbances based on data from the EUV Imaging Spectrometer (EIS) on board Hinode and the Interface Region Imaging Spectrograph (IRIS) that associates them with intermittent outflows and spicules produced at loop footpoint regions (for more details see De Pontieu & McIntosh 2010; Wang 2016). Numerical efforts have been made to tackle this issue by for example De Moortel et al. (2015). There, the authors performed forward modelling of the coronal emission from simulations using wave and periodic outflow drivers, and did not find any significant difference in the synthetic observational characteristics between the different drivers. In the current article, we deal with propagating intensity perturbations observed in active region fan loops, especially the three-minute oscillations. Such perturbations are considered to be relatively less ambiguous in their interpretation because of their regular periodic pattern lasting for hours to days, and because of their possible connection with the upward-propagating three-minute waves in the lower atmosphere. Therefore, we assume that they are the observational signatures of slow waves and hence we revisit some of the properties related to them.

Slow magnetoacoustic waves exhibit periods ranging from 2 to 10 min. In particular, in loops rooted in a sunspot umbra, their periodicity is about 3 min, and, for those rooted in the plage regions it is about 5 min (De Moortel et al. 2002). There is also a distinction between the periodicity found in active regions and that found in polar regions (for a detailed discussion see Banerjee et al. 2021, and references therein). The amplitude of the perturbations corresponds to a few percent of the background and their phase speeds range from 40 to 200 km s⁻¹, values that are comparable to the local sound speed. The perturbations also decay rapidly with height, with typical decay length scales of the order of tens of megameters (Nakariakov & Kolotkov 2020).

Slow waves are thought to originate in the lower atmosphere (Sych et al. 2009; Botha et al. 2011; Jess et al. 2012a; Krishna Prasad et al. 2015). For the three-minute waves some of the hypotheses on the origin include: externally driven by the global (broad-band) photospheric p-modes (Marsh & Walsh 2006; Jess et al. 2012a; Krishna Prasad et al. 2015; Zhao et al. 2016), magneto-convection (e.g., Jess et al. 2012b; Chae et al. 2017; Cho et al. 2019) and excitation above the cut-off frequency (Fleck & Schmitz 1991; Bogdan 2000). The generation and propagation of slow waves has also been the subject of many numerical studies. In particular, the mode conversion and energy transport have received considerable attention in 2D (Khomenko et al. 2009; Fedun et al. 2011; Santamaria et al. 2015) and 3D (Felipe et al. 2010; Mumford et al. 2015; Riedl et al. 2021) simulations. In these latter works, the authors obtained that slow waves that propagate along the magnetic field reach the upper atmosphere and the fast waves are reflected back, resulting in more acoustic energy propagating into the corona. Nevertheless, according to Riedl et al. (2021), only about 2% of the initial energy from the driver reaches the corona because of the strong damping due to the cut-off and to the geometric effects.

All these studies bring us to an important question about the slow magnetoacoustic waves specifically regarding the role played by the magnetic field in helping them reach coronal heights. The magnetic field plays an important role in the propagation of slow waves, not only acting as a wave guide but also modifying the cut-off frequency (Zhukov 2002). There is a strong consensus that the magnetic field inclination plays a crucial role in the frequency determination of for example penumbral waves because of the consequent reduction in the cut-off frequency (Reznikova et al. 2012; Reznikova & Shibasaki 2012; Kobanov et al. 2013; Jess et al. 2013; Yuan et al. 2014). Recently, Zurbriggen et al. (2020) showed, through an analytical calculation, that the cut-off frequency of slow magnetoacoustic-gravity waves not only depends on the inclination of the magnetic field but also on its intensity. The analytical results showed good agreement with the dominant periods found in an active region. Therefore, the cut-off frequency can be used as a seismological tool to determine the magnetic field and the temperature stratification.

Another seismological tool is the phase speed of the propagating disturbances in the plane of sky (POS). A typical method for studying them uses the so-called time–distance plot. Applying this technique to the stereoscopic observations of a coronal loop from the Extreme-Ultraviolet Imager instruments (EUVI) from STEREO A and B spacecraft, Marsh et al. (2009) estimated the true propagation phase speed from each spacecraft. The obtained values were $132.0_{-8.5}^{+9.9}\mathrm{km}{\mathrm{s}}^{-1}$ and $132.2_{-8.6}^{+13.8}\mathrm{km}{\mathrm{s}}^{-1}$. Assuming that these waves correspond to slow modes and under the condition of low plasma β, the corresponding plasma temperatures were estimated as $0.{84}_{-0.11}^{+0.13}\mathrm{MK}$ and $0.{84}_{-0.11}^{+0.18}\mathrm{MK}$, respectively. Shortly afterwards, Marsh & Walsh (2009) used spectroscopic diagnostic methods with data from Hinode/EIS and confirmed the temperature that was determined seismologically in the previous work, strengthening the slow mode interpretation. Almost simultaneously, Wang et al. (2009) determined the projected and Doppler speeds of propagating disturbances along a fan-like coronal structure observed in the 195 Å line using the same instrument. With these measurements, and applying linear wave theory, these latter authors obtained the inclination of the magnetic field, namely 59 ± 8°, a true propagation speed of 128 ± 25 km s⁻¹ , and a corresponding temperature of the loop of 0.7 ± 0.3 MK near the footpoint. Yuan & Nakariakov (2012) designed several techniques to measure the apparent propagation speeds from the time–distance maps in coronal imaging in a systematic way. They applied the techniques to quasi-periodic EUV disturbances propagating at a coronal fan-structure of an active region observed in the 171 Å bandpass from the Atmospheric Imager Assembly (AIA) on board Solar Dynamics Observatory (SDO), resulting in values that range from 48 to 52 km s⁻¹.

The observed apparent speed was constant in all of the works mentioned above. On the other hand, Sheeley et al. (2014), employing a running difference technique to track proper motions in XUV disturbances, detected accelerated patterns, from [35 − 45] km s⁻¹ to [60 − 100] km s⁻¹, in 2.6-min sunspot waves at coronal heights. These authors attributed this accelerated behaviour to the fading length of the height-time track as well as variations in the inclination of the plume above the sky plane. Also, Krishna Prasad et al. (2017) found accelerating slow magnetoacoustic waves along a coronal loop that show differential propagation speeds in two distinct temperature channels, namely the 131 Å and 171 Å from SDO/AIA. Using the inclination of the loop with respect to the line of sight (LOS) obtained from non-linear force-free magnetic field extrapolations these authors deprojected the measured phase speeds and estimated the corresponding local plasma temperatures. They found that the deprojected speeds still show an accelerated profile and that they are different for the two channels. Furthermore, the temperature in both channels increases with the distance along the loop and they found an appreciable difference between the two channels. Krishna Prasad et al. (2017) suggested that their results imply a multi-thermal, and consequently multi-stranded, structure of the loop.

With these findings in mind, in this work we analyse an event consisting of semi-circular, recurrent, arc-shaped intensity fronts that emerged from the active region (AR hereafter) NOAA 11243 on 2011 July 6. The intensity disturbances were captured by SDO/AIA in several EUV channels. The event was observed to persist for at least the first 8 h of the day, with neither flare activity nor coronal mass ejections reported nearby prior to or during the occurrence of the phenomenon. A periodicity analysis reveals the presence of a 2.5-min period at different heights of the AR near the footpoint of the loop rooted in the sunspot umbra. A kinematic analysis performed on EUV channels of SDO/AIA indicates an accelerated profile along the loop. This event was first analysed by Stekel et al. (2014), suggesting that the observed acceleration could be due to projection effects and the coronal oscillations may be originated in an umbral dot. With the aim of studying the projection effect caused by the inclination of the magnetic field, we performed numerical simulations. We modelled the coronal oscillations in a sunspot magnetic field configuration embedded in a gravitationally stratified atmosphere by inducing perturbations with an oscillatory driver near the photosphere. We calculated the projected propagation speeds from the forward-modelled intensity emission and compared the results with those from observations.

The paper is organised as follows. In Sect. 2, we describe the observations and give a detailed account of the observed event along with the periodicity analysis (Sect. 2.2) and the characterisation of wave propagation (Sect. 2.3). In Sect. 3, we introduce the numerical model and in Sect. 4 we display the results obtained from the synthetic intensity emission. We present a discussion on both the observational and numerical results in Sect. 5. Finally, we summarise our results and conclusions in Sect. 6.

2. Observations

For this study we consider space-based data from the AIA (Lemen et al. 2012) and the Helioseismic and Magnetic Imager (HMI, Schou et al. 2012) both on board the SDO (Pesnell et al. 2012). The AIA observations consist of cropped images of 124 × 80 arcsec² in five EUV bands (304 Å, 131 Å, 171 Å, 193 Å and 211 Å) and two UV bands (1700 Å and 1600 Å) from NOAA AR 11243 on 2011 July 6. The time sequence is from 00:00 UT to 02:00 UT with a temporal cadence of 12 s for the images taken in the EUV channels and 24 s for those in the UV spectrum. The leading spot of the active region is located at 14° N and 34° W at the beginning of the day (μ ≈ 0.83). The plate scale is 0.6 arcsec for both EUV and UV images. The HMI data correspond to the continuum intensity (near Fe I 6173 Å) and the LOS magnetogram, with a cadence of 45 s and a plate scale of 0.5 arcsec. All the images have been co-aligned, corrected for differential rotation, and brought to a common plate scale using the standard AIA software package available in IDL Solarsoft.

2.1. Data analysis

While inspecting a time-lapse sequence of 171 Å AIA images taken during the first 8 hours of 2011 July 6, we found a train of semi-circular and recurrent intensity fronts apparently emerging from NOAA AR 11243 towards the west (see movie in the online version). The arc-shaped, recurrent fronts appear to propagate through a bundle of loops apparently rooted in the umbra of the sunspot. Figure 1 shows snapshots of the AR as recorded by HMI and AIA on 2011 July 6, just at the beginning of the day. In the HMI images (top row, left and middle panels) and in the UV images of AIA (1700 Å, top-right panel, and 1600 Å, middle-left panel) of Fig. 1 we can see the sunspot associated to NOAA AR 11243. The HMI LOS magnetogram shows the value of the magnetic field intensity. The coronal channels (131 Å, in the right-middle panel and the 171 Å, 193 Å, and 211 Å in the bottom row) show the fan-like loop structure extending from the centre to the right of each snapshot. No flare activity or coronal mass ejections were reported in association with this phenomenon.

Fig. 1. Images of NOAA AR 11243 observed in HMI Continuum, HMI LOS Magnetogram, and AIA channels. The HMI and AIA UV channels show the sunspot associated with the active region while the coronal EUV channels exhibit the fan-like loop structure extending towards the right side. The primary ions of individual channels are listed in the title of each panel following Lemen et al. (2012). The temporal evolution of AIA 171 Å is available in the online movie.

2.2. Periodicity analysis

In order to identify the source and to characterise the observed periodicity of the phenomenon, we tracked the oscillations throughout the atmosphere, from the photosphere up to the corona. To that aim, we performed a periodicity analysis comprising different heights by computing the Lomb-Scargle periodograms (based on the fasper routine described in Press et al. 2002) of the intensity time series at each pixel of the region shown in Fig. 1 (ROI hereafter) considering 120 minutes of observations. From the periodograms obtained, we determined the maximum power value between 1.5 and 15 min and normalised the rest of the power spectrum to this value. Then we searched for all the local maxima within intervals of 1 min and chose the largest value among them. In Figs. 2–3 (from the second to the last column) we display this value of maximum power for every pixel in each range of periods between 2 and 5 min for the different wavelengths (from top to bottom row). The boundaries of the umbra and the penumbra are also shown in these plots. In the first column we show the corresponding observed intensities within the ROI.

Fig. 2. Snapshots of the ROI (first column) taken by SDO/AIA at approximately 00:00 UT for 1700 Å, 1600 Å, 304 Å, and 171 Å (from top to bottom). Maps (from second to fourth column) showing the power distribution obtained from the periodograms calculated for each pixel of the ROI over three one-minute intervals between 2 and 5 min. The second column corresponds to the main powers obtained in the interval [2–3] min, the third to [3–4] min, and the last one to [4–5] min. The colour scales in the different columns are normalised to the same value for comparison.

Fig. 3. Same as Fig. 2 but for wavelengths 193 Å and 211 Å.

We note that from the second column of Figs. 2–3, inside the umbra, there is a localised and intense power in the [2–3] min interval, whose area is growing as the atmospheric height increases (from top to bottom). In the lower atmosphere (1700 Å and 1600 Å), this area is coincident with the footpoints of the fan loops. As the height and the temperature increase through the 304 Å, 131 Å, and 171 Å channels, the area of the power enhancement expands, but a reduction in it is observed in the 193 Å and 211 Å channels. The growth of this area in the [2–3] min power could be attributed to the expansion of the flux tube through which the perturbation travels (see e.g., Jess et al. 2012a). In contrast, the shrinking of it in the 193 Å and 211 Å channels is likely a temperature effect as these channels observe relatively hotter plasma. It may be noted that the visibility of the associated loop structures is also limited in these channels (see Fig. 3). Surrounding the umbra, the power within the [3–4] min period band (third column) is enhanced, almost as much as that in the [2–3] min band in the photosphere and chromosphere. The area of this enhancement is also expanding outwards with increasing height but in the corona the power is much less compared to that in the [2–3] min period band. In the lower atmosphere (1700 Å and 1600 Å), the wave power within the [4–5] min band (fourth column) is more intense and ubiquitous compared to the other bands in both the penumbra and quiet region. On the other hand, in the higher atmosphere, especially in the coronal channels, the power within this band does not seem to display any clear enhancements. These results indicate that the oscillations we observe in the corona are localised and have dominant periods of between 2 and 3 min. Moreover, they also suggest that these oscillations could have originated in the lower atmosphere from where they propagate up to the corona through an expanding flux tube, which seems to be rooted in the sunspot umbra.

2.3. Kinematic characterisation

To determine the propagation speed of the observed propagating intensity disturbances, we tracked the intensity perturbations along a slightly curved slit of finite width located along one of the loops that seem to emanate from the active region. The chosen loop matches the region where the oscillation power is the largest (see Fig. 2). The extent and location of the slit are marked with the white lines in Fig. 4. The longitudinal extension of the loop is ≈10 150 km. The boundaries of both the umbra and the penumbra, as determined from the HMI continuum image, are also overlaid in this figure with red contours.

Fig. 4. Location of the analysed loop marked on the 171 Å image. The solid white lines outline the boundaries of the selected loop. The red contours represent the external boundaries of the umbra and penumbra of the sunspot.

By averaging the intensities across the loop, we construct time–distance (t − d) maps as displayed in Fig. 5 for the AIA 131 Å, 171 Å, 193 Å, and 211 Å channels. We refer to Krishna Prasad et al. (2012) for a detailed description of the method followed. As can be noted from the figure, a regular and recurrent ridge pattern is evident in all the channels. The inclination of the ridges in each t − d map reveals the projected propagation speed of the brightness disturbances in the POS as they propagate along the region defined by the slit. The apparent curvature in the inclination, prominent in 171 Å where the signal is strong up to the top end of the slit, is a signature of their apparent acceleration.

Fig. 5. Time–distance (t − d) maps obtained for the loop shown in Fig. 4 for the 131 Å, 171 Å, 193 Å, and 211 Å channels of AIA (from top to bottom).

To give an account of the kinematics of the phenomenon, a quantitative analysis of the t − d map patterns is required. We therefore estimated the POS speeds of the fronts using the method described in Tomczyk & McIntosh (2009). First, for each wavelength, we determined the time lags as a function of distance along the loop. For that purpose, we built a reference light curve from the average intensities across three consecutive rows where the oscillatory pattern is clearly visible. We then calculated the cross-correlation at all spatial locations with respect to that reference. Additionally, we estimated the errors on the cross-correlation (ΔC′) following Eq. (13) in Misra et al. (2018) which states:

$$\begin{array}{c}\hfill \mathrm{\Delta }{C}_{\mathit{XY}}^{\prime }=\frac{1-C_{\mathit{XY}}^2}{N^2\sqrt{{\sigma }_X^2{\sigma }_Y^2}}\sqrt{\sum _{k=-N/2}^{N/2}|{\stackrel{\sim }{X}}_k{|}^2{|{\stackrel{\sim }{Y}}_k|}^2},\end{array}$$(1)

where X and Y represent two light curves of N elements, σ_X and σ_Y are the standard variations of the respective light curves, $C_{\mathit{XY}}=c_{\mathit{XY}}/\sqrt{{\sigma }_X^2{\sigma }_Y^2}$ where c_XY is the cross-correlation between the two time series and $\stackrel{\sim }{X_k}$ and $\stackrel{\sim }{Y_k}$ are the discrete Fourier transforms of X and Y, respectively, in the frequency domain k. Therefore, the uncertainty in the correlation function depends on the value of the correlation and the amplitude of the correlated light curves. For each pair of correlated light curves, we took the location of maximum correlation along with the two preceding and the two following points. The correlation values at these locations were then fitted with a parabola function using the non-linear least squares method¹ considering the correlation errors. The vertex of the parabola corresponds to the time lag between that pair. Afterwards we fit a cubic B-spline function² to the obtained time lags to estimate the POS speeds. The time lags in all the channels and the corresponding fitted curves are shown in Fig. 6 (top panel) as a function of the distance along the loop. The vertical bars on each point represent the respective errors calculated from the uncertainties on the parameters of the parabola fitting. The spline fits shown by the solid lines are weighted by the inverse of these errors. Note that the error in the time lags is minimum near the reference row, at ≈3600 km, because the correlation is maximum there.

Fig. 6. Time lags and POS speeds from t − d maps. Top panel: time lags as a function of distance along the loop for the 131 Å, 171 Å, 193 Å, and 211 Å channels. The individual values are represented by coloured symbols, the associated errors are shown by the vertical bars, and the solid lines denote the corresponding spline fits to the data. Bottom panel: projected POS speeds derived from the spline fits to the time lag values.

As follows, we estimate the propagation speeds by computing the derivative of the distance along the loop with respect to the fitted time lags. The obtained POS speed values are exhibited in Fig. 6 (bottom panel) as a function of the distance. The different colours represent the values for the different wavelengths and the vertical bars correspond to the propagated errors estimated as v σ_tlag/|t_lag|, where v is the estimated speed, t_lag is the respective time lag, and σ_tlag is the corresponding error. Note that the uncertainties near the reference row are relatively large because t_lag ≈ 0. The speeds in all the wavelengths exhibit an accelerated profile in accordance with our visual prediction. Similar results were previously reported by Sheeley et al. (2014) and Krishna Prasad et al. (2017). Except for a few locations near the top of the loop segment (where the signal is low), the speed values mainly vary from 40 to 100 km s⁻¹. It is worth noting that the exact value of the speed (and in particular of the acceleration) is strongly sensitive to the fitting function used for the time lags (a small change in the fit could result in a larger change in the acceleration profile and hence of the local speed). This may indicate that it is difficult to distinguish the POS speeds in different channels. Nevertheless, while emphasising the accelerated behaviour, we must also bear in mind that in 193 Å and 211 Å the signal is low and the loop is not clearly discernible in the upper half of the slit. Therefore, the speeds calculated from these channels above ≈5000 km should be considered carefully.

3. Numerical model

In order to interpret the accelerated propagation observed in multiple channels, especially to decouple the effect of projection on the propagation speeds, we performed numerical simulations. To this end, we used the code Multi Advanced Non-ideal Code for High-resolution simulation of the solar Atmosphere (MANCHA3D, see Felipe et al. 2010; González-Morales et al. 2018; Khomenko et al. 2018), which solves non-linear non-ideal equations of the MHD for perturbations. The thermodynamic variables (density ρ, gas pressure p) and the magnetic field B were split into an equilibrium state and their departures from it, while the flow speed, v, was always treated as a (non-linear) perturbation. Consequently, the magnetohydrostatic (MHS) equilibrium is explicitly removed from the equations. This strategy conveys certain advantages for the simulations of waves in magnetic structures, because tiny numerical deviations from the MHS equilibrium do not grow in time, unlike in setups that rely on the full variables, guaranteeing a stable initial equilibrium configuration. We can summarise the equations that we used in this work as follows:

$$\begin{array}{cc}\hfill \frac{\partial {\rho }_1}{\partial t}& +\mathrm{\nabla }·\left[({\rho }_0+{\rho }_1){\mathit{v}}_1\right]={\left(\frac{\partial {\rho }_1}{\partial t}\right)}_{\mathrm{diff}}\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill \frac{\partial \left[({\rho }_0+{\rho }_1){\mathit{v}}_1\right]}{\partial t}& +\mathrm{\nabla }·[({\rho }_0+{\rho }_1){\mathit{v}}_1{\mathit{v}}_1^{\mathrm{T}}+(p_1+\frac{|{\mathit{B}}_1{|}^2}{2{\mu }_0}\hfill \end{array}$$(3)

$$\begin{array}{cc}& +\frac{{\mathit{B}}_0·{\mathit{B}}_1}{{\mu }_0})\mathbb{I}-\frac{1}{{\mu }_0}({\mathit{B}}_1{\mathit{B}}_1^{\mathrm{T}}+{\mathit{B}}_1{\mathit{B}}_0^{\mathrm{T}}+{\mathit{B}}_0{\mathit{B}}_1^{\mathrm{T}})]\hfill \\ \hfill & ={\rho }_1\mathit{g}+{\left(\frac{\partial [({\rho }_0+{\rho }_1){\mathit{v}}_1]}{\partial t}\right)}_{\mathrm{diff}}+\mathit{S}(t)\hfill \\ \hfill \frac{\partial e_1}{\partial t}& +\mathrm{\nabla }·[(e_0+e_1+p_0+p_1+\frac{|{\mathit{B}}_0+{\mathit{B}}_1{|}^2}{2{\mu }_0}){\mathit{v}}_1\hfill \end{array}$$(4)

$$\begin{array}{cc}& -\frac{1}{{\mu }_0}({\mathit{B}}_0+{\mathit{B}}_1)({\mathit{v}}_1·({\mathit{B}}_0+{\mathit{B}}_1))]\hfill \\ \hfill & =({\rho }_0+{\rho }_1){\mathit{v}}_1·\mathit{g}+{\left(\frac{\partial {\mathit{e}}_1}{\partial t}\right)}_{\mathrm{diff}}\hfill \\ \hfill \frac{\partial {\mathit{B}}_1}{\partial t}& =\mathrm{\nabla }\times [{\mathit{v}}_1\times ({\mathit{B}}_0+{\mathit{B}}_1)]+{\left(\frac{\partial {\mathit{B}}_1}{\partial t}\right)}_{\mathrm{diff}},\hfill \end{array}$$(5)

where g = −274 k̂ m s⁻² is the gravitational acceleration in the vertical direction z, S(t) is a time-dependent external force, and 𝕀 is the identity matrix. The subscripts “0” and “1” indicate the background variables and their perturbations, respectively. The background state should satisfy the MHS equilibrium. The total energy per unit volume is defined as $e=\frac{1}{2}\rho v^2+\frac{p}{\gamma -1}+\frac{B^2}{2{\mu }_0}$ and satisfies that e = e₀ + e₁. The terms subscripted with “diff” are numerical diffusion terms required for numerical stability and are computed according to Vögler et al. (2005) and Felipe et al. (2010).

In the numerical study presented here, we performed 2D simulations on a Cartesian grid of 1344 × 672 cells distributed over a domain of [4 × 10⁴]×[2 × 10⁴] km in the horizontal (x) and vertical (z) directions, respectively. This results in a resolution of ≈30 km. The horizontal coordinate varies from −2 × 10⁴ to 2 × 10⁴ km and the vertical coordinate from 0 to 2 × 10⁴ km, where 0 km coincides with the base of the photosphere.

3.1. Magnetohydrostatic equilibrium

As mentioned above, the MANCHA3D code ensures the MHS equilibrium of both stratified atmosphere and magnetic configuration when the initial condition is in such a state. In order to achieve background equilibrium, we need to satisfy −∇p₀ + ρ₀g + j₀ × B₀ = 0, where j₀ is the current density (other variables have already been defined). To model the initial atmosphere stratification we chose a smooth temperature profile given by:

$$\begin{array}{c}\hfill T_0(z)=a_{0}tanh\left(\frac{z-a_1}{a_2}\right)+a_3,\end{array}$$(6)

where a₀ = 5 × 10⁵ K, a₁ = 3 × 10³ km, a₂ = 5 × 10² km, and a₃ = 5.1 × 10⁵ K. While the parameters a_0, 3 constrain the maximum and minimum values of the temperature, a₁ controls the height at which the transition region begins and a₂ regulates the ‘sharpness’ of it. As a result of the combination of the parameters, we obtain a minimum temperature of 10⁴ K in the lower atmosphere and 10⁶ K in the corona. This profile can resemble the thermal inversion in the corona and was previously used, for example, in the analytical and observational work by Zurbriggen et al. (2020) where the authors studied how the different stratifications and magnetic field can influence the cut-off periods of magnetoacoustic-gravity waves. Once we had established the temperature stratification, we were able to obtain the pressure and the density profiles from the hydrostatic equilibrium equation:

$$\begin{array}{c}\hfill \mathrm{\nabla }{p}_0={\rho }_0\mathit{g},\end{array}$$(7)

accompanied by the ideal gas law with constant mean molecular weight of μ = 0.5. A value of 1.5 × 10¹ dyn cm⁻² was taken for the gas pressure at the base of the atmosphere. As the stratification is along the vertical direction, Eq. (7) becomes a unidimensional equation along z and in the horizontal direction, x, the thermodynamic variables are homogeneous. We show the initial temperature, density, and pressure in Fig. 7 (left panel) as a function of z in the middle of the domain.

Fig. 7. Initial conditions for the simulation. Left panel: temperature (red line), density (green line), gas (blue line), and magnetic pressure (black line) stratification along the z direction considering a constant value of 274 m s−2 for gravity. Right panel: symmetric potential magnetic field configuration. The background colour denotes the intensity of the magnetic field while the solid black lines show some representative magnetic field lines.

Furthermore, to ensure the MHS equilibrium, that is, to satisfy the force equation, we consider a force-free (j₀ × B₀ = 0) magnetic field as follows:

$$\begin{array}{cc}& B_{x0}(x,z)=-B_{\mathrm{base}}exp(-k_{z}z)sin(k_{z}x)\hfill \end{array}$$(8)

$$\begin{array}{cc}& B_{z0}(x,z)=-B_{\mathrm{base}}exp(-k_{z}z)cos(k_{z}x),\hfill \end{array}$$(9)

where k_z = π/L is the wave number, L = 2 × 10⁴ km, and B_base = 1 × 10³ G is the value of the magnetic field in the lower atmosphere. The magnetic field, displayed in Fig. 7 (right panel), is symmetric with respect to the central vertical axis and it decreases exponentially with height, resulting in a vanishing field at a large distance. This simple model represents a first-order approximation of a sunspot (Aschwanden 2004) and is similar to that used in Santamaria et al. (2015; without the null point) and to the extrapolated magnetic field in Jess et al. (2013) derived from the observations.

3.2. Wave generation

Based on the oscillation properties described in Sect. 2.2, a localised region with dominant periods in the [2–3] min interval and the oscillations in the corona appear to be generated in the lower atmosphere. Therefore, we employed a source S(t) in the vertical direction of the momentum Eq. (3) to drive the waves. This source acts as an external force in a localised region close to the bottom boundary of the simulation domain, namely the lower atmosphere. The source term depends on several free parameters that are free to vary during the simulations. In the case presented here, we used the following expression:

$$\begin{array}{c}\hfill \mathit{S}(t)=A{\rho }_0\mathit{g}exp(-\frac{{(x-x_0)}^2}{{\sigma }_x^2}-\frac{{(z-z_0)}^2}{{\sigma }_z^2})sin(\omega t),\end{array}$$(10)

where A = 5 × 10⁻⁵, x₀ = 0, z₀ = 20 × dz ≈ 595 km, σ_x = 33 × dx ≈ 982 km, σ_z = 3 × dz ≈ 89 km, and ω = 2π/180 s⁻¹, with dz = dx ≈ 30 km. The source self-consistently produces the perturbations in the velocity and other thermodynamic variables besides the magnetic field. In Fig. 8, we show the variation of the density (top panel) and the velocity (bottom panel) at t = 198 s.

Fig. 8. Snapshots of the simulation at t = 198 s. Upper panel: shows the perturbation of the density over the background density (ρ1/ρ0) and the bottom panel the vertical velocity (vz1). The solid black lines represent some of the background magnetic field lines.

As the source is located at the axis of the magnetic structure where the magnetic field is almost vertical, it produces slow magneto-acoustic waves as one would expect in a low-plasma-β atmosphere. These waves propagate along the field (the magnetic field lines are indicated as black lines in Fig. 8) through the chromosphere to the transition region and to the corona, albeit with a significant reflection from the transition region due to the sharp gradient in the sound speed.

It is important to mention that we kept the source amplitude very low in the simulations (though fully non-linear equations were solved), because our aim is to reproduce a linear wave behaviour. Observations do not show significant evidence for non-linear effects in the propagation of slow waves along coronal loops. Non-linearities affect the wave propagation speed, possibly making it larger than the local sound speed. By keeping the driver amplitude low, we avoid this non-linear effect in our simulations.

The boundary conditions on the horizontal direction are periodic. For the vertical boundaries, both up and down, we use a sponge layer to absorb the coming waves. While MANCHA3D has a perfectly matched layer (PML) boundary condition (Berenger 1994, 1996; Hu 1996, 2001; Parchevsky & Kosovichev 2007, 2009), we find that it does not work well for the extremely long-wavelength perturbations in the corona. Instead, we find a simple sponge layer to be a better choice (see González-Morales et al. 2019). For the sponge, we used 10 grid points at the bottom and 60 at the top.

4. Numerical results

In order to facilitate a direct comparison with the observations, we generated synthetic images in the EUV spectrum from the simulation results. As the oscillation amplitudes in the simulation are relatively low, we first re-scaled the perturbations in all the quantities by a factor of ten. This scaling is allowed for perturbations in a linear regime, and brings the oscillation amplitudes close to the observed values. After that we extracted a 2D flux tube from the whole numerical domain to emulate one of the fan loops observed in the EUV channels (see Fig. 1). The flux tube shown in Fig. 9 is defined by all the plasma contained between two magnetic field lines defined by Eqs. (8)–(9).

Fig. 9. Location of the flux tube extracted from within the simulation domain. The external boundaries are defined by two field lines from the magnetic field Eqs. (8)–(9). The LOS is aligned with the z-axis of the numerical domain and the POS with the x-axis. The observer location is displayed at the top of the domain.

Once we had extracted the flux tube, we computed the specific intensity in the EUV spectral lines using the Leuven Forward Modelling code for coronal emission (FoMo, Van Doorsselaere et al. 2016). Briefly, given the numerical density, temperature, and velocity, this code calculates the specific intensity for a certain monochromatic spectral line in the POS integrating the emissivity along the LOS. The POS together with the LOS are not necessarily aligned with the axes of the numerical simulations. Therefore, for generating the synthetic images, it is possible to choose the spectral line from a specific telescope and to give different viewing angles. In this study, we chose to forward model the intensity for the four EUV wavelength channels, the 131 Å, 171 Å, 193 Å, and 211 Å of AIA. We assumed that the LOS is aligned with the z-axis of the numerical domain and one of the axes of the POS is parallel to the x-axis. In terms of FoMo parameters, this implies that the viewing angles selected in the code are (0,0), and therefore the observer is located perpendicularly above the simulation domain (see Fig. 9). We then obtained the intensity t − d maps in the four wavelengths, as presented in Fig. 10. As can be seen, all these maps exhibit an accelerated recurrent pattern with no significant differences between the wavelengths. This is due to the narrow temperature distribution of the isolated fluxtube (around 1 × 10⁶ K) that we considered for calculating the synthetic intensity.

Fig. 10. Synthetic time–distance maps (t − d) in the 131 Å, 171 Å, 193 Å, and 211 Å channels (from top to bottom) obtained from the forward modelled intensities along the selected flux tube considering a projection direction parallel to the x-axis as shown in Fig. 9.

To study the propagation speeds in each spectral line, we first determined the time lags as a function of distance in the same way as for the observations. The only difference here is that the reference light curve is selected from a single row rather than averaging three rows, as these data are not noisy. The top panel of Fig. 11 shows the time lags computed for each passband as a function of the projected distance along the x-axis, together with the fitted spline curves. The different coloured symbols represent the values obtained for each wavelength and the vertical lines on them show the associated errors calculated in the same way as for observations. Note that the error bars here are smaller (comparable to the size of the symbol) than in observations mainly because the correlation values are higher. The solid lines mark the respective spline fits to the data but they are not discernible due to the high density of points (because of the high resolution of the simulations). Subsequently, taking the time derivative on the distance using the fitted time lag values, we determine the propagation speeds for each spectral line. We display the obtained values along with their respective errors in Fig. 11 (bottom panel). The errors are again large near the reference row (≈3600 km) as the time lags are close to zero and the propagated errors are inversely proportional to them. The speed values range from 40 to 160 km s⁻¹ and represent the apparent propagation speeds seen by an observer located above the simulated active region. They reach a maximum value relatively close to the sound speed at a distance of about 7300 km from the base of the flux tube. However, note that the distance is given in the projected coordinate, that is, in the x-direction, and, as can be seen from Fig. 9, the flux tube already reaches the top of the domain by this location. Therefore, the apparent speed beyond this location does not correspond to the upward propagation of waves but rather represents the speed at which the waves arrive at the adjacent locations at the top of the domain. This explains why the apparent propagation speed decreases at distances greater than 7300 km. In addition, the time lags obtained for different wavelengths are almost the same and therefore there is no distinction in speeds between them. However, the acceleration dependence on the distance is evident.

Fig. 11. Time lags and POS speeds from synthetic t − d maps. Top panel: time lags (coloured symbols) obtained from the synthetic t − d maps for each wavelength as a function of the projected distance along the loop. Vertical bars (comparable to the size of the symbols) on these values indicate the associated errors from the quadratic fit performed to obtain the time lags. The solid lines denote spline fits to the data, which are almost indistinguishable. Bottom panel: propagation speeds obtained from the spline fits to time lag values along with the associated errors.

5. Discussion

The observed apparent speeds calculated in Sect. 2.3 for the 131 Å, 171 Å, 193 Å, and 211 Å passbands vary from 40 to 100 km s⁻¹ and exhibit accelerated profiles. These values are smaller than the slow magnetosonic wave speeds calculated from the characteristic temperature of each channel. Interpreting the observed oscillations as slow waves that are generated in the lower atmosphere and propagate upwards following the magnetic field (see analysis from Sect. 2.2), their accelerated profiles may imply an increase in the local plasma temperature with distance. However, a changing inclination of the magnetic field with respect to the LOS would result in a similar behaviour because of the projection effect. Moreover, this would also explain why the speed values are lower compared to the sound speed. Therefore, a combination of both the physical and the non-physical effects is possible in observations. The observed speed values are similar to the 171 Å passband measurements described by Krishna Prasad et al. (2017) and to the sunspot waves events analysed by Sheeley et al. (2014). However, in contrast to Krishna Prasad et al. (2017), there is no clear distinction between the speeds determined in 131 Å and 171 Å channels. This suggests that the differences we find in the observed speeds across different wavelengths perhaps do not indicate a multi-thermal structure of the coronal loop. Nevertheless, this is difficult to determine conclusively, as the speeds are subject to large uncertainties that depend on the method and the differences were marginal. In any case, the absence of a clear temperature dependence in the propagation speeds and the diminished visibility of the loop itself (see Fig. 3) in the hotter channels (193 Å, and 211 Å) imply a narrow temperature distribution within the loop.

In our simulations, we excite slow waves that propagate parallel to the field lines at the sound speed c_s. To demonstrate this, we display in Fig. 12 the background sound speed c_s (red symbols) for all the locations along the middle field line of the flux tube together with its projection (in black symbols) considering the inclination of the magnetic field with respect to the LOS. The speed calculated from the synthetic t − d map for 171 Å is also shown for reference.

Fig. 12. Propagation speed for 171 Å (orange symbols), background sound speed cs (red symbols), and its projection (black symbols) considering the inclination of the magnetic field as a function of the projected distance along the loop.

We note that for the coronal conditions employed in the model, c_S = 167 km s⁻¹ (assuming a mean molecular weight μ = 0.5, for a temperature of 1 × 10⁶ K). As there is hardly any variation in temperature in the coronal part of the model (see from the figure that c_s is constant for the locations that correspond to the flux tube), the lower (than c_s) propagation speeds along with the accelerated pattern obtained from the synthetic t − d map highlight the projection effect caused by the inclination of the waveguide with respect to the LOS. Indeed, the speed calculated in 171 Å is comparable to the projected sound speed at most of the locations. The discrepancy near the bottom of the flux tube is because the spline fit departs from the data in these locations. In addition, the absence of any difference in speeds between different passbands (see Fig. 11) is also a consequence of having a uniform (coronal) temperature in our setup, which again supports our interpretation of a near-isothermal cross-section for the observed loop. Furthermore, the values that we calculate from the synthetic maps are comparable to the values found in observations. Therefore, we can attribute the accelerated profile from observations to the projection effect.

6. Summary and conclusions

We performed a comprehensive analysis of an oscillatory event observed in an active region fan loop system in the EUV channels of SDO/AIA. We calculated the distribution of periods in the ROI at different heights of the solar atmosphere, from near the photosphere to the corona. We obtain that the oscillations observed in the corona, especially in 171 Å, have a dominant period of between 2 and 3 min and are confined to a small region along the fan. By tracking these oscillations down to the transition region to the chromosphere and near the photosphere, we find that this period is also present and dominant in the sunspot umbra where apparently the footpoints of the coronal fan are rooted. This was previously suggested by Stekel et al. (2014). In addition, the region where the oscillations are localised is expanding towards the corona, indicating an expansion of the waveguide. We also performed a kinematic analysis along one of the fan loops observed in the EUV channels 131 Å, 171 Å, 193 Å, and 211 Å of AIA. We calculated the propagation speed of the perturbation in the POS as a function of the distance along the loop, with the origin of the loop in the footpoints located in the umbra of the associated sunspot. We find that the apparent propagation speed is slower than the local sound speed but exhibits an accelerated pattern in all four channels. Similar acceleration profiles in EUV wavelengths were previously identified and measured by Sheeley et al. (2014) and Krishna Prasad et al. (2017). However, in contrast to Krishna Prasad et al. (2017), there is no clear distinction between the speeds determined in the different wavelengths, suggesting a narrow temperature distribution across the loop.

With these observational features as constraints, we performed 2D numerical simulations using the MANCHA3D code. In our numerical model, we considered an initial configuration in MHS equilibrium. To achieve this, we considered a gravitationally stratified atmosphere using a smooth temperature profile from the photosphere up to the corona. For the magnetic field, we used a force-free configuration with a strong magnetic field in the photosphere resembling an active region. We drove waves in the vertical direction of the momentum equation using a source term with a Gaussian shape along the horizontal x and the vertical z axes near the photosphere. The period of the driver was selected to be 3 min. Using the forward modelling code FoMo, we calculated the intensity in the same four EUV channels, 131 Å, 171 Å, 193 Å, and 211 Å, from SDO/AIA for a 2D flux tube extracted from within the simulation domain. We then computed the propagation speed of the waves along the magnetic flux tube as seen on the projected x-axis, that is, the apparent speed seen by an observer located above the active region looking along z-axis. We obtain accelerated patterns in the wave propagation comparable to the ones detected in the observations with the speed values also in a similar range.

In observations, both thermal and projection effects could be present, producing an accelerated pattern, but in simulations the acceleration is purely due to the projection effect caused by the inclination of the magnetic field. Therefore, given the similarities of our numerical and observational results, we attribute the observed accelerated propagation speed to projection effects. Nevertheless, we emphasise that care should be taken when assigning a physical meaning to the acceleration in the observed propagation speeds of slow waves.

Movies

Movie 1 associated with Fig. 1 (AIA_171) Access here

Acknowledgments

MVS and TVD acknowledge support from the European Research Council (ERC), grant agreement No. 724326 under the European Union’s Horizon 2020 research and innovation programme. SKP is grateful to FWO Vlaanderen for a senior postdoctoral fellowship (No. 12ZF420N). TVD was also supported by C1 grant TRACESpace of Internal Funds KU Leuven. TVD has benefited from the funding of the FWO Vlaanderen through a Senior Research Project (G088021N). GS was supported by NASA Grant 80NSSC19K1261. EK acknowledges support from the Spanish Ministry of Science through the project PGC2018-095832-B-I00 and from the European Research Council through the project ERC-2017-CoG771310-PI2FA. The data used here are courtesy of NASA/SDO and the AIA and HMI science teams. MANCHA3D is a public collaborative code that can be requested from http://research.iac.es/proyecto/PI2FA/pages/codes.php

References

All Figures

Fig. 1. Images of NOAA AR 11243 observed in HMI Continuum, HMI LOS Magnetogram, and AIA channels. The HMI and AIA UV channels show the sunspot associated with the active region while the coronal EUV channels exhibit the fan-like loop structure extending towards the right side. The primary ions of individual channels are listed in the title of each panel following Lemen et al. (2012). The temporal evolution of AIA 171 Å is available in the online movie.

In the text

Fig. 2. Snapshots of the ROI (first column) taken by SDO/AIA at approximately 00:00 UT for 1700 Å, 1600 Å, 304 Å, and 171 Å (from top to bottom). Maps (from second to fourth column) showing the power distribution obtained from the periodograms calculated for each pixel of the ROI over three one-minute intervals between 2 and 5 min. The second column corresponds to the main powers obtained in the interval [2–3] min, the third to [3–4] min, and the last one to [4–5] min. The colour scales in the different columns are normalised to the same value for comparison.

In the text

	Fig. 3. Same as Fig. 2 but for wavelengths 193 Å and 211 Å.
In the text

	Fig. 4. Location of the analysed loop marked on the 171 Å image. The solid white lines outline the boundaries of the selected loop. The red contours represent the external boundaries of the umbra and penumbra of the sunspot.
In the text

	Fig. 5. Time–distance (t − d) maps obtained for the loop shown in Fig. 4 for the 131 Å, 171 Å, 193 Å, and 211 Å channels of AIA (from top to bottom).
In the text

Fig. 6. Time lags and POS speeds from t − d maps. Top panel: time lags as a function of distance along the loop for the 131 Å, 171 Å, 193 Å, and 211 Å channels. The individual values are represented by coloured symbols, the associated errors are shown by the vertical bars, and the solid lines denote the corresponding spline fits to the data. Bottom panel: projected POS speeds derived from the spline fits to the time lag values.

In the text

Fig. 7. Initial conditions for the simulation. Left panel: temperature (red line), density (green line), gas (blue line), and magnetic pressure (black line) stratification along the z direction considering a constant value of 274 m s−2 for gravity. Right panel: symmetric potential magnetic field configuration. The background colour denotes the intensity of the magnetic field while the solid black lines show some representative magnetic field lines.

In the text

	Fig. 8. Snapshots of the simulation at t = 198 s. Upper panel: shows the perturbation of the density over the background density (ρ1/ρ0) and the bottom panel the vertical velocity (vz1). The solid black lines represent some of the background magnetic field lines.
In the text

	Fig. 9. Location of the flux tube extracted from within the simulation domain. The external boundaries are defined by two field lines from the magnetic field Eqs. (8)–(9). The LOS is aligned with the z-axis of the numerical domain and the POS with the x-axis. The observer location is displayed at the top of the domain.
In the text

	Fig. 10. Synthetic time–distance maps (t − d) in the 131 Å, 171 Å, 193 Å, and 211 Å channels (from top to bottom) obtained from the forward modelled intensities along the selected flux tube considering a projection direction parallel to the x-axis as shown in Fig. 9.
In the text

Fig. 11. Time lags and POS speeds from synthetic t − d maps. Top panel: time lags (coloured symbols) obtained from the synthetic t − d maps for each wavelength as a function of the projected distance along the loop. Vertical bars (comparable to the size of the symbols) on these values indicate the associated errors from the quadratic fit performed to obtain the time lags. The solid lines denote spline fits to the data, which are almost indistinguishable. Bottom panel: propagation speeds obtained from the spline fits to time lag values along with the associated errors.

In the text

	Fig. 12. Propagation speed for 171 Å (orange symbols), background sound speed cs (red symbols), and its projection (black symbols) considering the inclination of the magnetic field as a function of the projected distance along the loop.
In the text