© ESO 2025

1. Introduction

It has long been debated if the propagation, damping, and dissipation of magnetohydrodynamic (MHD) waves is a significant contributor to the coronal heating problem, and if so, where these waves originate from. Magnetohydrodynamic theory (Ferraro & Plumpton 1958; Bel & Mein 1971; Nakagawa et al. 1973; Zhugzhda 1983; Zhugzhda & Dzhalilov 1984; Goedbloed & Poedts 2004; Goossens et al. 2011) has given us an extensive toolkit to study magnetoacoustic-gravity wave modes, including important insights into how linear and non-linear modes propagate, change, and interact with each other in a non-uniform medium. Specifically, it is important to note the evanescence (damping) of different waves due to criteria known as cut-off frequencies. For example, for acoustic modes this limit depends on the speed of sound, and greatly reduces the purely acoustic wave activity in the upper atmosphere. Furthermore, around the β = 1 transition, both magnetic and gas pressure can act as the dominant restoring force for magnetoacoustic waves. The distinction of modes in this region is less clear, and it is possible that the nature of the wave may change around this transition, and so mode conversion is usually assumed in this region. In the upper atmosphere, the magnetic field inclination also plays an important role in the trapping and reflection of waves (Nakagawa et al. 1973; Newington & Cally 2010). As the magnetic field becomes more vertical with height, strict constraints develop on the propagation of magnetoacoustic waves into the corona, and most wave energy is instead reflected back into the chromosphere.

Whether or not magnetosonic waves have enough energy to heat the corona, it is still unclear how or if they can transport it from the photosphere upwards to the corona. Idealised numerical simulations (Bogdan et al. 2003; Chatterjee 2020; Kumar & Kumar 2020; Riedl et al. 2021a,b; Cally 2022; Yelles Chaouche et al. 2023) have strengthened our understanding of these modes of oscillation, by tracking and studying driven waves through inhomogeneous media. Previous studies have shown how energy may be transferred across boundaries via mode conversion, and estimates of wave flux energy have been used to support or disband the significance of specific wave modes in heating the corona (Carlsson & Bogdan 2006; Fossum & Carlsson 2005; Newington & Cally 2010; Cally 2016; Tarr et al. 2017; Liu et al. 2023; McMurdo et al. 2023). From this, it is concluded that mainly acoustic modes dominate in the photosphere, whilst magnetoacoustic and Alfvén modes are present in the upper chromosphere and corona. These idealised simulations, however, are still far from the physical reality, lacking crucial physics and self-consistency that may influence wave behaviour. Ultimately, the next steps are to confirm such mode conversions and damping in realistic simulations, and eventually observations. Furthermore, the estimations of wave energy deposition have important implications for the boundary conditions of solar wind predictions (Van der Holst et al. 2014; Sishtla et al. 2022; Perri et al. 2022). It is possible to estimate the Alfvén wave energy transferred into the solar wind in MHD simulations using the Elssässer variables, and more recently, modified Q-variables that also encompass fast and slow magnetoacoustic wave energy (Van Doorsselaere et al. 2024). Alternatively, these wave fluxes can be observationally constrained (Morton et al. 2012; Sharma & Morton 2023); however, an accurate determination is the subject of ongoing work.

In realistic simulations, wave modes interact with each other, form shocks, experience turbulence, and change nature, making it challenging to follow wave activity through the computational domain. The altitude of the β = 1 transition is inhomogeneous, and contains both chromospheric and coronal plasma, such that waves travelling in any direction may cross this transition more than once. Furthermore, the shocks themselves can form transient structures, such as jets, which further excite or transmit waves through the transition region (Heggland et al. 2011; Kuniyoshi et al. 2024; Skirvin et al. 2024). When wave modes are self-driven by convective motions, it is difficult to differentiate individual waves from other time-dependent dynamics. Yadav et al. (2022) and Enerhaug et al. (2024) suggest possible methods of detecting wave activity across field lines and flux surfaces, respectively, by capturing characteristics specific to each magnetoacoustic wave mode. However, these typically do not extract the actual oscillations of wave activity. Raboonik et al. (2024) propose a method of decomposition using the eigenenergies of magnetoacoustic waves. This is able to capture both non-linear effects and phase-mixing in realistic 3D simulations, but still has limitations, such as the assumption of adiabaticity (no dissipation mechanisms), and the exclusion of gravity in the derivation of the method. Here, we elaborate on the work of Finley et al. (2022), using a combination of discrete Fourier transforms (DFTs) and component decomposition to infer wave activity in a physically realistic 3D simulation.

2. Numerical simulation

2.1. The Bifrost code

The simulation we have analysed was produced using the 3D radiative magnetohydrodynamic (rMHD) code, Bifrost, described in detail in Gudiksen et al. (2011). Bifrost incorporates a shallow but self-consistent convection zone that drives the motions above, up to the mid-corona. Within each atmospheric level, additional physics modules aim to capture the physical processes and create a realistic simulation. Thermal conductivity is applied in the upper atmosphere, and follows Spitzer conductivity, calculated through an explicit hyperbolic approximation (Spitzer & Seeger 1963; Rempel 2017; Cherry et al. 2024). In optically thin regions (upper chromosphere and corona), radiative transfer is assumed to only depend on density and temperature, with a transfer function calculated with ionisation, recombination, and collisional excitation rates from pre-computed data given by Shull & van Steenberg (1982), Arnaud & Rothenflug (1985), and Judge & Meisner (1994). Full radiative transfer in four opacity bins (Hayek et al. 2010) is then used for optically thick regions, predominantly in the photosphere and lower chromosphere. These calculations assume a static medium, and local thermal equilibrium (LTE), except for hydrogen ionisation in the chromosphere, which is treated as non-LTE. The radiative transfer equation is solved iteratively in order to include important scattering effects. The main MHD equations are solved using a third-order Hyman predictor-corrector scheme, from a staggered Cartesian grid with fifth- and sixth-order spatial operatives. The simulation was run with variable time-stepping, and localised hyper-diffusive terms to increase stability.

2.2. Set-up

We analysed a quiet Sun simulation that consists of a 12 × 12 Mm² horizontal extent, with a resolution of dx = dy = 11.7 km. The computational domains extends from the convection zone, approximately 2.5 Mm below the photosphere, to the mid-corona, ∼8 Mm above with a varying vertical resolution of dz = 6 − 35 km. In total, the domain consists of 1024³ grid points; a higher resolution than most published simulations from Bifrost. There are periodic boundaries in x and y, with open upper and lower boundaries. The upper boundary uses characteristic boundary conditions (see Gudiksen et al. 2011) that aim to transfer material through the boundary with minimal reflection. Even so, we avoided analysis within 2 Mm of this boundary in order to remove boundary effects. On the lower boundary, the entropy was maintained such that the solar effective temperature reaches 5773 K, although in reality it varies between 5717 K and 5761 K.

This study analyses the simulation starting at 1.5 hours (referred to as t = 0 seconds) of solar time, with snapshots available at a cadence of 10 seconds, continuing until 2.1 hours (t = 1890 seconds) of solar time. This pertains to 189 snapshots of data, spanning 36 minutes.

We have chosen three horizontal slices in which to proceed with analysis in this paper, shown in Figure 1. These correspond roughly to the upper chromosphere (2 Mm), the end of the transition region into the corona (4 Mm), and the low corona (6 Mm).

Fig. 1. Average vertical temperature profile (black) over the whole domain. The dashed lines display the horizontal planes that were used and analysed in this study, at z = 2, 4, and 6 Mm. The dotted line represents the plane used in Figure 2, where τ = 1. The grey lines show the evolution of the temperature profile across the analysed time span.

Fig. 2. Vertical magnetic configuration where τ = 1, at time = 136.8 minutes.

2.3. Overview

The simulation evolves a highly non-idealised structure, visible in the vertical magnetic field structure in Figure 2. The self-consistent convection zone continually drives the motions further up in the simulation box such that the magnetic field is constantly evolving. The average structure of the atmosphere consists of photospheric granulation at approximately 0 Mm, a temperature minimum at 500 km, and the plasma-β,

$$\begin{array}{c}\hfill \beta =\frac{P}{B^2/2{\mu }_0}\propto \frac{c_s^2}{v_{\mathrm{A}}^2},\end{array}$$(1)

passes unity (β = 1) in the lower chromosphere on average (see Figure 3). In time and space, the β = 1 level varies considerably between the photosphere, z = 0 Mm, and the corona, z >  6 Mm, and forms a highly corrugated surface. Thus, it is not so easy to reach general conclusions about the nature of a wave travelling along or through a horizontal plane. This is especially true for two prominent magnetic features in the domain (see Figure 3), which exhibit torsional ‘swirling’ velocities and twisted magnetic fields, expanding from the upper photosphere into the low corona. The motions of such features are described in depth in Tziotziou et al. (2023). These structures, henceforth referred to as swirls, elevate the β = 1 layer up to 6 Mm, on the boundaries of the swirl. These thin boundary areas could be important as a possible route for waves to be transferred to the low-β plasma from the lower atmosphere. It is also important to note that adiabaticity does not hold everywhere. This is seen by the differences in p and ρ throughout the upper atmosphere, showing that the ideal gas assumption generally does not hold. When subject to non-adiabatic conditions, both linear and non-linear waves can lose energy through processes such as damping (Carbonell et al. 2004).

Fig. 3. Ratio of the sound and Alfvén speeds at heights z = 2, 4, and 6 Mm at t = 1000 seconds. Regions around equipartition, outlined by the solid line, are most susceptible to containing phase mixing and damping processes. Two prominent swirling features are highlighted by the boxed regions.

3. Methodology

Inspired by the analysis in Finley et al. (2022), we used the DFT on the temporal and spatial evolution of a quantity, g(x, y, t), in a horizontal plane. The quantity was mapped to Fourier space, in terms of the horizontal wavenumbers, $\stackrel{\sim }{{\nu }_i}=k_i/2\pi$, and the frequencies, f = ω/2π,

$$\begin{array}{cc}& \mathcal{F}({\stackrel{\sim }{\nu }}_x,{\stackrel{\sim }{\nu }}_y,f)=\hfill \\ \hfill & {\displaystyle \sum _{n_x=1}^{N_x}}\left({\displaystyle \sum _{n_y=1}^{N_y}}\left({\displaystyle \sum _{n_t=1}^{N_t}}g(n_x,n_y,n_t)e^{-i\frac{2\pi f{n}_t}{N_t}}\right)e^{-i\frac{2\pi {\stackrel{\sim }{\nu }}_x{n}_x}{N_x}}\right)e^{-i\frac{2\pi {\stackrel{\sim }{\nu }}_y{n}_y}{N_y}},\hfill \end{array}$$(2)

where n_x, n_y, n_t represent the indices of the spatial and temporal directions, and N_x, N_y, N_t denote the total number of points in each direction. The range of the spatial and temporal frequencies are given by the spatial and temporal step-size, respectively, for the maximum frequency, and the size of the spatial or temporal domain for the minimum frequency. Therefore, the minimum frequency sampled is 1/L = 0.083 Mm⁻¹ for the spatial directions, and 1/T = 0.52 mHz for the temporal direction. These are also the respective sampling resolutions.

For a given set of filtered horizontal wavenumbers or temporal frequencies (frequency bin), the localised contributions to the quantity, g, in the original spatial-temporal domain may be reconstructed using an inverse DFT. We propose that localised wave activity may be detected in regions of the simulation domain by mapping the frequencies that reconstruct the most accurate strength at each point. Per index, the significance of a frequency bin, b, is determined through the absolute difference of the reconstructed signal, S​_b, to the original signal,

$$\begin{array}{c}\hfill {\mathrm{Diff}}_b[n_x,n_y,n_t]=|S_b[n_x,n_y,n_t]-g[n_x,n_y,n_t]|.\end{array}$$(3)

Therefore, the most significant signal, has the smallest absolute error. We note that if the signal is not centred on zero, it is important to combine the zeroth frequency bin with the other individual bins, in order to accurately compare the signals.

Figure 4 displays the somewhat trivial result for a signal with a single frequency, where the most significant frequency (MSF) calculation is comparing individual frequencies. When there is only one dominating frequency, the MSF will be exactly the input frequency, since all other frequencies will have zero amplitudes. When two or more waves travel through a domain, it is important to repeat the MSF calculation to capture the second-most significant frequency (SSF), as seen in Figure 5. Here, the signal comprises two frequencies: the main wave has a large amplitude and low frequency, whilst a smaller-amplitude higher-frequency wave runs across it, creating a more complex signal. The MSF calculation predominantly captures the low frequency wave. To repeat this process for the SSFs, we removed the signal produced by the composite MSFs from the original signal, leaving the point-wise residual, Diff_MSF. The MSF algorithm may then be applied to the residual, resulting in the SSF. We can see from the third panel of Figure 5 that the output of the MSF and SSF signals combined returns a signal very similar to the original one, except for in a few places around turning points. At these turning points, the distribution of MSFs is found to be either the zero frequency (relating to the equilibrium value), or the frequency bin that contains all frequencies higher than 1.50 Mm⁻¹, which could relate to numerical noise in the signal, for example. In more complex instances, the amplitudes of these higher frequencies can build up in the highest bin and give a false dominant signal. It is therefore important to acknowledge this in the results and do further analysis before drawing conclusions on the activity of higher frequencies. We note that by increasing the range of the individual frequencies considered in the MSF decomposition, the effect of false dominant signals in the highest bin will be decreased but the decomposition demands a greater computational expense. We note that the spread of frequencies detected in each MSF calculation may vary dependent on the difference of frequency and amplitude between the two waves, but in any case, after two calculations, we find that the two dominant frequencies are clear.

Fig. 4. ‘Proof of concept’ with an isolated signal in the xy plane, with wavenumber kx = 0.58 Mm−1, and no frequency. The first panel shows a 1D slice in the y axis, taken at an arbitrary time. The middle panel shows the amplitude spectrum from the DFT, averaged across the y and time axes, where there is a sharp peak at the 0.58 Mm−1, whilst the third panel shows the most significant signal per pixel in the domain. In grey, the MSF, ω, is shown, and is zero as expected.

Fig. 5. Composite signal of two frequencies, kx = 0.166 Mm−1 and kx = 0.751 Mm−1. The first panel shows the original signal, the second panel shows the results of the most and second most significant calculations, and the final panel shows the final output of the significant frequencies. The original frequency is shown by the dotted lines for reference.

This method allows us to go beyond the standard DFT power spectrum in two ways. Firstly, it can show where a frequency becomes dominant in a specific region of the simulation domain where the dynamics are continuously changing. Secondly, it allows us to create inverse ω-frequency outputs in a spatial plane, such that we can locate regions of the domain space that contain fluctuations at certain frequencies and wavenumbers.

Finally, it is a general characteristic of DFTs that they only decompose a signal into frequencies with constant amplitudes. Therefore, it is difficult to capture decaying signals, such as the ones imperative to energy transfer in the upper solar atmosphere, through damping and dissipation. In Figure 6, we show how damped signals can be detected by calculating the MSFs. When the amplitude of the wave is large, the MSF is confidently set to the correct frequency, but as the amplitude diminishes, the spread of ‘significant’ frequencies increases, with the power shifting away from the actual frequency, and towards the lowest (0 Mm⁻¹) and highest (> 1.5 Mm⁻¹) frequencies. This is a distinct characteristic for a damped wave; therefore, we can infer the damping of signals in regions where a strong signal dissolves into pixellated high and low frequencies.

Fig. 6. Most significant frequencies calculated for a signal with a damped amplitude.

4. Results

We performed a preliminary analysis using the MSF method on a realistic simulation to demonstrate its potential to detect localised wave activity and go beyond some of the limitations of the DFT decomposition. We exploited several variables in the simulation data to draw preliminary conclusions on the type of wave as purely acoustic, magnetoacoustic or purely magnetic (Alfvénic). We did not search for specific modes, which would have required a more extensive analysis and so was left as a topic of future work. Instead, we propose that oscillations found in pressure, p, and velocity parallel to the magnetic field,

$$\begin{array}{c}\hfill u_{\Vert }=\mathbf{u}·\mathbf{B},\end{array}$$(4)

should relate to longitudinal (acoustic) wave characteristics, whilst oscillations in the magnetic field, B, and perpendicular velocity,

$$\begin{array}{c}\hfill {\mathbf{u}}_{\perp }=\mathbf{v}-{\mathbf{v}}_{\Vert },\end{array}$$(5)

will, in general, relate to magnetic waves. We chose the velocities to be parallel and perpendicular to the magnetic field in light of the complex magnetic structuring found in the simulation. In this way, the velocities have some relation to the results from analysis on idealised magnetic structures with fully orthogonal decompositions such as Riedl et al. (2019) and Skirvin et al. (2024), which use simplified magnetic geometries. However, both transverse variables remain multi-dimensional and must be mapped to scalar variables for the analysis below. Therefore, we used the vertical magnetic field, B_z, since this component will oscillate for the majority of oscillations of the full magnetic field, and we used the magnitude of the perpendicular velocity to capture perpendicular oscillations to the magnetic field. However, this does not allow us to distinguish between rotational and lateral oscillations. Therefore, we also discuss fluctuations in the magnitude of the magnetic field, which may hint at sausage or non-linear modes, in which the magnetic field is compressed or expanded due to oscillations. We note that fast and slow modes may display oscillations in all variables, or just in the dominant characteristics (i.e. acoustic and magnetic, respectively, in the low-β regime). Future work will explore the physical mechanisms at play for each result in a more in-depth analysis that is outside the scope of this paper.

4.1. Power spectral density and pixel plots

Usually, the results of a DFT are displayed as a power spectral density (PSD), where

$$\begin{array}{c}\hfill \mathrm{PSD}({\stackrel{\sim }{\nu }}_x,{\stackrel{\sim }{\nu }}_y,f)={|\mathcal{F}({\stackrel{\sim }{\nu }}_x,{\stackrel{\sim }{\nu }}_y,f)|}^2\end{array}$$(6)

shows which frequencies contain the most power over the entire domain. High values of the PSD either correspond to the sum of many fluctuations covering a large area of the domain, or more localised fluctuations with large amplitudes in the domain. Figure 7 concerns the velocity parallel to the magnetic field at z = 4 Mm, and shows the frequency against the horizontal wavenumber,

Fig. 7. Power spectral density plot for u∥ (top), and corresponding MSF counts (bottom). Taken for the horizontal plane z = 4 Mm. On the right plot, the outermost values represent the total counts for all other frequencies.

$$\begin{array}{c}\hfill {\stackrel{\sim }{\nu }}_h={({\stackrel{\sim }{\nu }}_x^2+{\stackrel{\sim }{\nu }}_y^2)}^{1/2}.\end{array}$$(7)

The DFT PSD (top panel) is compared to an equivalent plot for the MSF method (bottom panel). For this, we counted the number of pixels with each combination of MSF and the most significant wavenumber (MSW) for a specific snapshot. In all of the examples in this paper, unless explicitly stated, we have taken a snapshot at time t = 1000 seconds of the analysed data and considered each discrete frequency sampled by the DFT individually (each frequency bin contains only one frequency, with width 0.525 mHz). In the case of the perpendicular wavenumber, the bins were created such that each wavenumber bin was centred on the individual wavenumbers in x, ${\stackrel{\sim }{\nu }}_x$, sampled by the DFT. Each bin therefore contains 2−3 wavenumbers and has a width of 0.08 Mm⁻¹.

The MSF method captures the most powerful areas of the PSD in Figure 7, confirming that the method works on the most important fluctuations in the domain. In the mid-power range, the MSF method highlights a group of larger wave numbers with a frequency around 4 mHz whose importance is obscured by the PSD. The right-most column (${\stackrel{\sim }{\nu }}_h>2.4$ Mm⁻¹) and similarly the upper-most row (f >  19.5 mHz) also show strong signals. However, this is because they contain all the higher frequencies (up to 50 mHz) and wavenumbers (up to 42.6 Mm⁻¹), respectively, which are not considered individually, and instead the associated power stacks up in a single bin. The advantage of the MSF calculations is that a pixel plot, such as in Figure 8, may be used to show the location of the MSFs in the domain, and by extension, the nature by which the power of the PSD is created. Figure 8 highlights the areas where fluctuations of a specific frequency begin to dominate in the domain, either locally or globally. For example, we see an overwhelming development of MSFs in the range 3−6 mHz as height increases. By z = 4 Mm, there are expansive and continuous patches where these frequencies dominate, covering 45% of the domain. This compliments the results from the PSD in Figure 7 that the most power is centred on these frequencies. By z = 6 Mm, the coverage has reached 50%, whilst in the photosphere (z = 0, not shown) it is only ∼25%. The same activity is not observed in the vertical magnetic field in Figure 9, and so it is likely that these fluctuations relate to longitudinal pressure modes. Although outside the scope of this paper, we note that since we now have location information, it would be possible to estimate the wave flux energy more accurately, by only integrating over the regions covered by the specific frequencies, instead of just taking an average over the whole domain. Thus the wave energy for longitudinal waves would be calculated with

$$\begin{array}{c}\hfill F_{\mathrm{long}}=\frac{1}{2}{\displaystyle \int {\int }_{\mathbf{S}}}\rho v_p^2{c}_s\mathrm{d}\mathbf{S},\end{array}$$(8)

Fig. 8. MSF in the velocity parallel to the magnetic field, u∥, for the chromosphere and upper atmosphere.

Fig. 9. MSF for the vertical magnetic field, Bz, in the chromosphere and upper atmosphere. The rectangles represent the locations shown in later figures: the dashed rectangle relates to Figure 12, and the solid line represents the slice in Figure 14.

where v_p is the velocity in the direction of propagation, and S is the area covered by the specific frequencies. Since the direction of propagation is unknown, an upper limit could be estimated with the magnitude of the velocity vector, |v|.

4.2. Matched frequencies

We can speculate on the nature of other fluctuations using the correlations between MSFs in the magnitude of the magnetic field, parallel and perpendicular velocity, and pressure. We use the example of magnetoacoustic waves in the upper atmosphere. Magnetoacoustic waves have both pressure and magnetic variations, as well as transverse and longitudinal oscillations. The top panel of Figure 10 counts the number of pixels where both parallel velocity and pressure have the same MSF at that location, ±one frequency step (0.53 mHz), representing longitudinal fluctuations. The bottom panel counts the number of pixels where the same is true between the perpendicular velocity and the magnitude of the magnetic field, |B|, representing transverse oscillations.

Fig. 10. Number of matched MSFs between pressure and parallel velocity (upper), and the magnitude of the magnetic field and perpendicular velocity (lower), at different heights. The frequencies are considered matched if they are within ± ∼ 0.53 mHz (1 frequency step) of each other. 3 and 5 mHz frequencies are highlighted by the vertical dotted lines.

There are strong enhancements – and therefore strong correlations – in the longitudinal and transverse variables at all heights in the 3 − 6 mHz range. At z = 2 Mm there are two distinct peaks in both the longitudinal and transverse components at around 3.5 and 5.26 mHz. This may be relevant in the context of observational results of global oscillations at 3 and 5 minute periods. The longitudinal peaks broaden and become one large enhancement centred on 3.16 mHz by 6 Mm. In contrast, the two peaks remain more distinct in the transverse component at all heights, although the enhancement seems to narrow such that the two peaks lie closer together, centred around 4 mHz. By z = 4 Mm, most of the domain lies above the equipartition level, and hence in the low-β regime (see Figure 3). It is therefore likely that the longitudinal oscillations are connected to magnetoacoustic waves above this height. The continuity of enhancement of frequencies in the range 3 − 6 mHz through all the heights suggest that the global 3 and 5 mHz pressure modes from the lower atmosphere are the drivers of the 3 − 6 mHz enhanced magnetoacoustic waves seen in the upper atmosphere.

From Figures 8 and 9, we have previously shown that the locations of the fluctuations in the acoustic and magnetic components in general do not match exactly in location in the upper atmosphere and this is confirmed in Figure 11. It is an interesting feature that we nevertheless see a similar number of pixels showing transverse oscillations as longitudinal at this height if they do not appear in similar locations. Using the correlated pixel plots in Figure 11, we can distinguish areas where the magnetic and acoustic MSFs do overlap. Matching the results of all four variables is an ambitious task, and we would expect that the less dominant oscillations of fast and slow magnetoacoustic modes would more likely be captured in the secondary or even tertiary iterations of the MSF calculations. Considering this, it is not unexpected that there are only a few overlapping regions to be found in Figure 11, mostly in the 3−6 mHz range. Nevertheless, there are large areas that have dense – although not necessarily overlapping – populations of both components (some examples are highlighted in Figure 11), suggesting the presence of magnetoacoustic waves in these areas. There is a notable coupling between both longitudinal and transverse components on the boundaries of both swirls. However, the matched MSF here corresponds to the highest frequency bin. As is discussed in Section 3, this bin is capable of giving false MSFs due to the combined result of the many frequencies contained in this bin. Therefore it is important to consider the physical characteristics on the swirl boundaries. We note that these regions do not have weak magnitudes in the velocities or in pressure, relative to the average of the domain, as is shown for parallel velocity in the right-most panel of Figure 11. Therefore, these high-frequency MSFs are unlikely to be due to noise or other artefacts from the method. On the other hand, the magnitude of the magnetic field can be low in these areas. Even so, the matching of three strong-valued variables inspires further analysis.

Fig. 11. Regions of matched MSFs between the gas pressure and parallel velocity (left) and the magnitude of the magnetic field and perpendicular velocity (middle) at z = 4 Mm. The frequencies are considered matched if they are within ± ∼ 0.53 mHz (1 frequency step) of each other. The right panel shows the results from the left two panels overlaid to show overlapping regions. ‘Magnetic’ fluctuations are shown in red, whilst ‘acoustic’ fluctuations are in blue. Here, the magnitude of the frequencies is omitted. The magnitude of parallel velocity is shown for reference in the background. For all plots, MSF’s of 0 mHz have been removed as a trivial result.

4.3. Small-scale processes

In contrast to the MSFs of the parallel velocity, the vertical magnetic field, shown in Figure 9, shows global fluctuations at the lowest frequencies in the upper atmosphere. We see concentrations of fluctuations with a broad range of frequencies localised along the magnetic field boundaries in the domain. In fact, for the 3 − 6 mHz frequencies in B_z, the coverage steadily decreases from ∼30% at z = 2 Mm, to ∼20% by z = 6 Mm, with strong concentrations around the edges of both swirls at z = 4 Mm. The behaviour of these concentrations is readily identified using the MSF method as it separates structure in the spatial and frequency domains.

Furthermore, lines of ‘speckled’ high-frequency MSFs are present in both the parallel velocity and vertical magnetic field, seen in the orange-yellow range in Figure 9. These networks of higher frequencies are observed in many of the variables, but most structured in the vertical magnetic field, B_z, along the edges of the swirls. In these regions, the MSF changes over very short distances (∼per pixel) in the spatial domain, creating the speckled pattern. It is unlikely that this is a physical phenomenon; instead it is a ‘shortcoming’ of the method that we can exploit to detect other phenomena. There are a number of reasons these patterns could arise. Firstly, in regions where the original signal is close to zero relative to the average of the domain, it is important to acknowledge the presence of noise in the simulation. The noise is most likely to be captured by the top ‘highest-frequencies’ bin (i.e. displayed in yellow), although it is also possible in the high individual frequency bins. However, around the edge of the swirls, for example, the magnetic field is only near zero in very small regions where the polarity changes, and the velocities remain strong, so the presence of these speckled regions is likely to be the result of a different phenomenon. Instead, we refer to the behaviour discussed in Section 3. Aside from noise, these regions could either indicate a smaller, secondary fluctuation with a different frequency, as is described in Figure 5, where the MSF method sometimes also captures the SSF, with some loss of precision, or it indicates the dissipation of a wave, as in Figure 6. There is a slight difference to the idealised example of the latter used in Section 3: in this case we are taking the temporal frequency, f, and relating it to spatial location, and so we rely on the assumption that the MSFs will show damping across the orthogonal direction in a similar fashion to how they react in the direction parallel to the given frequency.

In the case of the thicker areas of speckled MSFs for B_z on the boundary of the large swirl, at z = 4 Mm, the SSFs (not shown) provide no evidence of a secondary fluctuation in these regions. Further analysis favours the dissipation theory: the B_z signal appears to be damped in a quasi-oscillatory motion from the centre of the swirl, at ∼(8 Mm, 5 Mm), to the outside boundary located vertically downwards at Y ≈ 2 Mm. This is displayed in the top right panel of Figure 12 as a 1D vertical slice through the lower part of the swirl. The bottom right panel shows that this observed oscillation coincides with chaotic switching of MSFs, in a similar fashion to the distinct pattern for damping described in Section 3. In this region we also see several peaks in the viscous heating term, Q_visc, which measures the transfer of energy from kinetic to heat in the plasma through viscous forces. There is an increased spatial average across the slice compared to other areas in the domain, as is seen in the lower left panel. In the region of the observed dissipation, the peaks align exactly with each dip in the B_z oscillation, suggesting that energy is transferred from the wave into heat in this region. In contrast, other regions close to the swirl boundary display a thinner region of much stronger viscous heating, coinciding with the thin speckled regions. In these regions, we do not see oscillations, but only a sharp decrease in the magnitude of B_z to near zero. This observation, coupled with the strong viscous heating suggests that the damping rate here is much greater than the frequency of the wave, such that no oscillations are observed.

Fig. 12. 1D slice at X ≈ 8.20 Mm through the vertical magnetic field, Bz, at z = 4 Mm shows that the speckled high frequencies capture wave dissipation across the swirl boundary between Y ≈ 1.76 Mm and Y ≈ 4.69 Mm. The viscous heating, Qvisc, also correlates with this conclusion. The sliced region is emphasised by the yellow box in the left panels, and the damped oscillation region is highlighted in red on the right panels.

The above example demonstrates the kind of further analysis that is needed in order to distinguish the nature of any speckled region observed. It is quite likely that there are more instances of damping to be found in this simulation. Similarly, it would be possible to track these damping mechanisms through the atmosphere. If the above speculations are validated through further analysis, this could provide an invaluable tool for measuring the contribution to coronal heating from the dissipation of waves.

5. Discussion

5.1. Viscous heating in swirl boundaries

Despite this paper primarily serving as an introduction to the MSF method, the results from Sections 4.2 and 4.3 can be used to understand how energy is transported via waves and deposited in the low corona. In the swirl boundaries, the plasma-β remains greater than 1, as is seen in Figure 3. Therefore, the plasma displays characteristics of chromospheric plasma, but is pulled up to coronal heights by the motions within the swirls. Although greater than 1, the plasma-β remains close to unity, providing ideal conditions in the plasma for mode conversion, via phase-mixing or otherwise. This allows waves with both strong magnetic and acoustic characteristics to exist. As the oscillations diffuse outwards from the swirl boundaries, the plasma approaches low-β once more, and so the conditions for magnetosonic waves to propagate change, causing damping of the wave amplitude from external forces such as viscosity. The excess energy is transferred into heat, and found in the viscous heating term, Q_visc, which ultimately raises the temperature of the plasma. This conclusion requires further analysis to validate; namely, by following plasma along the spiral boundary and ‘tracking’ specific waves. This method provides the tools to do this in future work.

5.2. Contributions of the MSFs

Figure 13 shows the contribution,

$$\begin{array}{c}\hfill C=(1-{\mathrm{Diff}}_{\mathrm{MSF}})\times 100\%,\end{array}$$(9)

of the MSF signal to the original signal, where Diff_MSF is the difference between the two signals, calculated in the same way as in Equation (3), but including the sign of the difference. We do not expect one signal to account for the entire strength of the original signal, due to the extent to which neighbouring processes are continuously interacting. In theory, both a contribution greater than 100% and a negative contribution are possible, since the MSF signal may overestimate the original signal, or be counteracting the effect of an overestimated signal through destructive interference (with an opposite sign to the original signal). In our results, there are no negative contributions. This is because in most cases the overestimated signal will be closer to the original than the counteracting signal. The likelihood that the method provides an accurate result is lower for areas where the contribution is further from 100% (i.e. close to 0% or 200%). Nevertheless, Figure 13 shows that there are only a few areas where this is the case (shown by the black patches), and that these areas do not coincide with the locations of the results addressed in this study.

Fig. 13. Contribution of the MSF signal to the original signal of parallel velocity (upper) and magnitude of the magnetic field (lower). The insets provide examples of the small patches of negligible contributions by the MSF signal.

For structures where the original signal at a specific time is much larger than the average deviation, such as in the strong velocities and magnetic field of the swirls (that move in space and time), the MSF method can accurately depict the spatial and timescales, but not necessarily the correct contribution of the MSF or MSW signal to the original. Figure 14 is an example of this for the large swirl. This is because the MSF relies on the constant amplitudes taken from the DFT, which will be an average of the fluctuations in the domain for any given wavenumber. Therefore it is not the purpose of the MSF method to accurately detect the contribution of large, time-evolving structures in a global domain. If this is required, one could apply the DFT to a smaller region around the structure in question.

Fig. 14. Same 1D slice as in Figure 12, at X ≈ 8.20 Mm, for the vertical magnetic field, Bz, across the largest swirl. The signal created from the MSFs is shown by the dashed line, whilst the signal created when the SSFs are added is represented by the dotted line. The oscillations along the peak from X ≈ 6 Mm to X ≈ 7.5 Mm are not captured by either signal.

5.3. Future applications

In this paper, we have established an algorithm for detecting localised wave activity based on the results of DFTs. This method provides a new way to probe wave activity in the simulation domain, and has already produced results that show potential for improving our understanding of heating in the solar atmosphere. It surpasses the global results given by a DFTs power spectrum, by providing information on dominant frequencies in localised areas. This gives us valuable insights into the small-scale activity surrounding areas of interest such as swirl boundaries and the β = 1 transition. Furthermore, the behaviour of the MSFs for damped oscillations may be exploited for detection of damping mechanisms, which is not possible in the DFT alone. The algorithm is an alternative to the wavelet method (Jess et al. 2023) that does not require as many assumptions about the wave modes. It gives more precise results for the frequency and location than decomposition methods such as empirical mode decomposition (Huang et al. 1998) and single value decomposition (Santolík et al. 2003).

We have focused our attention on the MSFs and MSWs in the horizontal spatial domain. In the upper atmosphere, the majority of oscillations propagate vertically upwards, driven by the convective motions below, and more analysis in the vertical direction would be an important addition to future work. This could help to refine the results found on the boundaries of the swirl. A 3D MSF analysis could be used to find the direction of propagation in the boundary regions, and track fluctuation origin, propagation, and dissipation in time and as a function of height. The results of this could be compared to the ones of traditional wave guides, such as the ones described by Enerhaug et al. (2024). A similar application would be detailing the wave activity along coronal loops in a 3D simulation. This could extend the work of, for example, Riedl et al. (2021a) by including external physics in the plasma around the loop. Advanced 3+-dimensional calculations require significant memory and computational power due to the highly resolved spatio-temporal coverage in this simulation. By adapting the DFTs or domain-space used, we could reduce this requirement, which could also improve the method in areas of large structures, as was discussed above. Further analysis of the MSFs and their distribution in the domain will give us valuable insights into the origin and heating of the solar wind.

Acknowledgments

This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie [grant agreement N° 945371]. It is also supported by the Research Council of Norway through its Centres of Excellence scheme, project number 262622, and through grants of computing time from the Programme for Supercomputing. A.J.F. received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 810218 WHOLESUN).

References

All Figures

	Fig. 1. Average vertical temperature profile (black) over the whole domain. The dashed lines display the horizontal planes that were used and analysed in this study, at z = 2, 4, and 6 Mm. The dotted line represents the plane used in Figure 2, where τ = 1. The grey lines show the evolution of the temperature profile across the analysed time span.
In the text

	Fig. 2. Vertical magnetic configuration where τ = 1, at time = 136.8 minutes.
In the text

	Fig. 3. Ratio of the sound and Alfvén speeds at heights z = 2, 4, and 6 Mm at t = 1000 seconds. Regions around equipartition, outlined by the solid line, are most susceptible to containing phase mixing and damping processes. Two prominent swirling features are highlighted by the boxed regions.
In the text

Fig. 4. ‘Proof of concept’ with an isolated signal in the xy plane, with wavenumber kx = 0.58 Mm−1, and no frequency. The first panel shows a 1D slice in the y axis, taken at an arbitrary time. The middle panel shows the amplitude spectrum from the DFT, averaged across the y and time axes, where there is a sharp peak at the 0.58 Mm−1, whilst the third panel shows the most significant signal per pixel in the domain. In grey, the MSF, ω, is shown, and is zero as expected.

In the text

	Fig. 5. Composite signal of two frequencies, kx = 0.166 Mm−1 and kx = 0.751 Mm−1. The first panel shows the original signal, the second panel shows the results of the most and second most significant calculations, and the final panel shows the final output of the significant frequencies. The original frequency is shown by the dotted lines for reference.
In the text

	Fig. 6. Most significant frequencies calculated for a signal with a damped amplitude.
In the text

	Fig. 7. Power spectral density plot for u∥ (top), and corresponding MSF counts (bottom). Taken for the horizontal plane z = 4 Mm. On the right plot, the outermost values represent the total counts for all other frequencies.
In the text

	Fig. 8. MSF in the velocity parallel to the magnetic field, u∥, for the chromosphere and upper atmosphere.
In the text

	Fig. 9. MSF for the vertical magnetic field, Bz, in the chromosphere and upper atmosphere. The rectangles represent the locations shown in later figures: the dashed rectangle relates to Figure 12, and the solid line represents the slice in Figure 14.
In the text

	Fig. 10. Number of matched MSFs between pressure and parallel velocity (upper), and the magnitude of the magnetic field and perpendicular velocity (lower), at different heights. The frequencies are considered matched if they are within ± ∼ 0.53 mHz (1 frequency step) of each other. 3 and 5 mHz frequencies are highlighted by the vertical dotted lines.
In the text

Fig. 11. Regions of matched MSFs between the gas pressure and parallel velocity (left) and the magnitude of the magnetic field and perpendicular velocity (middle) at z = 4 Mm. The frequencies are considered matched if they are within ± ∼ 0.53 mHz (1 frequency step) of each other. The right panel shows the results from the left two panels overlaid to show overlapping regions. ‘Magnetic’ fluctuations are shown in red, whilst ‘acoustic’ fluctuations are in blue. Here, the magnitude of the frequencies is omitted. The magnitude of parallel velocity is shown for reference in the background. For all plots, MSF’s of 0 mHz have been removed as a trivial result.

In the text

Fig. 12. 1D slice at X ≈ 8.20 Mm through the vertical magnetic field, Bz, at z = 4 Mm shows that the speckled high frequencies capture wave dissipation across the swirl boundary between Y ≈ 1.76 Mm and Y ≈ 4.69 Mm. The viscous heating, Qvisc, also correlates with this conclusion. The sliced region is emphasised by the yellow box in the left panels, and the damped oscillation region is highlighted in red on the right panels.

In the text

	Fig. 13. Contribution of the MSF signal to the original signal of parallel velocity (upper) and magnitude of the magnetic field (lower). The insets provide examples of the small patches of negligible contributions by the MSF signal.
In the text

	Fig. 14. Same 1D slice as in Figure 12, at X ≈ 8.20 Mm, for the vertical magnetic field, Bz, across the largest swirl. The signal created from the MSFs is shown by the dashed line, whilst the signal created when the SSFs are added is represented by the dotted line. The oscillations along the peak from X ≈ 6 Mm to X ≈ 7.5 Mm are not captured by either signal.
In the text