Waves in the lower solar atmosphere: the dawn of next-generation solar telescopes

Abstract

Waves and oscillations have been observed in the Sun’s atmosphere for over half a century. While such phenomena have readily been observed across the entire electromagnetic spectrum, spanning radio to gamma-ray sources, the underlying role of waves in the supply of energy to the outermost extremities of the Sun’s corona has yet to be uncovered. Of particular interest is the lower solar atmosphere, including the photosphere and chromosphere, since these regions harbor the footpoints of powerful magnetic flux bundles that are able to guide oscillatory motion upwards from the solar surface. As a result, many of the current- and next-generation ground-based and space-borne observing facilities are focusing their attention on these tenuous layers of the lower solar atmosphere in an attempt to study, at the highest spatial and temporal scales possible, the mechanisms responsible for the generation, propagation, and ultimate dissipation of energetic wave phenomena. Here, we present a two-fold review that is designed to overview both the wave analyses techniques the solar physics community currently have at their disposal, as well as highlight scientific advancements made over the last decade. Importantly, while many ground-breaking studies will address and answer key problems in solar physics, the cutting-edge nature of their investigations will naturally pose yet more outstanding observational and/or theoretical questions that require subsequent follow-up work. This is not only to be expected, but should be embraced as a reminder of the era of rapid discovery we currently find ourselves in. We will highlight these open questions and suggest ways in which the solar physics community can address these in the years and decades to come.

1 Introduction

Understanding the energy flow through the Sun’s dynamic and tenuous atmosphere has long been a scientific interest for the global astrophysical community. The challenge of identifying the source(s) responsible for the elevated multi-million Kelvin temperatures in the solar corona has produced two main theoretical mechanisms. The first is via magnetic reconnection—the so-called ‘DC’ heating mechanism. Here, the continual re-configuration of the omnipresent magnetic fields that populate the Sun’s atmosphere allow the production of intense thermal heating as the magnetic energy is converted through the process of reconnection, producing dramatic flares that often release energies in excess of $10^{31}$ ergs during a single event (Priest 1986; Priest and Schrijver 1999; Shibata and Magara 2011; Benz 2017). However, such large-scale solar flares are relatively rare, and hence cannot supply the global background heating required to continuously maintain the corona’s elevated temperatures. Instead, there is evidence to suggest that the frequency of flaring events, as a function of their energy, is governed by a power-law relationship (Shimizu and Tsuneta 1997; Krucker and Benz 1998; Aschwanden et al. 2000; Parnell and Jupp 2000), whereby smaller-scale micro- and nano-flares (with energies $\sim 10^{27}$ ergs and $\sim 10^{24}$ ergs, respectively) may occur with such regularity that they can sustain the thermal inputs required to maintain the hot corona. Many modern numerical and observational studies have been undertaken to try and quantify the ubiquity of these faint reconnection events, which often lie at (or below) the noise level of current-generation facilities (Terzo et al. 2011). Due to the difficulties surrounding the extraction of nanoflare characteristics embedded within the noise limitations of the data, only tentative evidence exists to support their global heating abilities of the outer solar atmosphere (Viall and Klimchuk 2013, 2015, 2016, 2017; Jess et al. 2014, 2019; Bradshaw and Klimchuk 2015; Tajfirouze et al. 2016a, b; Ishikawa et al. 2017, to name but a few recent examples).

The second energy-supplying mechanism for the Sun’s outer atmosphere involves the creation, propagation, and ultimately dissipation of wave-related phenomena—often referred to as the ‘AC’ heating mechanism (Schwarzschild 1948). The specific oscillatory processes responsible for supplying non-thermal energy to the solar atmosphere have come under scrutiny since wave motions were first discovered more than 60 years ago (Leighton 1960; Leighton et al. 1962; Noyes and Leighton 1963a). Of course, such early observations were without the modern technological improvements that enhance image quality, such as adaptive optics (AO; Rimmele and Marino 2011) and image reconstruction techniques, including speckle (Wöger et al. 2008) and multi-object multi-frame blind deconvolution (MOMFBD; van Noort et al. 2005). As a result, many pioneering papers documenting the characteristics of wave phenomena in the lower solar atmosphere relied upon the study of large-scale features that would be less effected by seeing-induced fluctuations, including sunspots and super-granular cells, captured using premiere telescope facilities of the time such as the McMath-Pierce Solar Telescope (Pierce 1964) at the Kitt Peak Solar Observatory, USA, and the National Science Foundation’s Dunn Solar Telescope (DST; Dunn 1969), situated in the Sacramento Peak mountains of New Mexico, USA (see Fig. 1).

Fig. 1

Images depicting the construction of National Science Foundation facilities some 50 years apart. Panels b, d, e display construction stages of the Dunn Solar Telescope, which was first commissioned in 1969 in the Sacramento Peak mountains of New Mexico, USA. Panels a, c, f depict similar stages of construction for the Daniel K. Inouye Solar Telescope, which acquired first-light observations in 2019 at the Haleakal$\bar{\text{a}}$ Observatory on the Hawaiian island of Maui, USA. Images courtesy of Doug Gilliam (NSO) and Brett Simison (NSO)

Even at large spatial scales, Doppler velocity and intensity time series from optical spectral lines, including Fe i (Deubner 1967), H$\alpha$ (Deubner 1969), Ca ii (Musman and Rust 1970), C i (Deubner 1971), and Na i (Slaughter and Wilson 1972) demonstrated the ubiquitous nature of oscillations throughout the photosphere and chromosphere. Through segregation of slowly-varying flows and periodic velocity fluctuations, Sheeley and Bhatnagar (1971) were able to map the spatial structuring of wave power in the vicinity of a sunspot (see Fig. 2), and found clear evidence for ubiquitous photospheric oscillatory motion with periods $\sim 300$ s and velocity amplitudes $\sim 0.6$ km s$^{-1}$. Such periodicities and amplitudes were deemed observational manifestations of the pressure-modulated global p-mode spectrum of the Sun (Ulrich 1970; Leibacher and Stein 1971; Deubner 1975; Rhodes et al. 1977), where internal acoustic waves are allowed to leak upwards from the solar surface, hence producing the intensity and velocity oscillations synonymous with the compressions and rarefactions of acoustic waves.

Fig. 2

Image reproduced with permission from Sheeley and Bhatnagar (1971), copyright by Springer

Observations of the photospheric Fe i absorption line, showing the sum of blue- and red-wing intensities (displayed in a negative color scale; top), the total measured Doppler velocities across the field-of-view (middle-top), the slowly varying component of the plasma flows (middle-bottom), and the Doppler velocity map arising purely from oscillatory motion (bottom). The region of interest includes a large sunspot structure (left-hand side), and shows ubiquitous oscillatory signatures with periods $\sim 300$ s and velocity amplitudes $\sim 0.6$ km s$^{-1}$.

Difficulties arose in subsequent work, when the measured phase velocities of the waves between two atmospheric heights were too large to remain consistent with a purely acoustic wave interpretation (Osterbrock 1961; Mein and Mein 1976). It was not yet realized that the 5-min oscillations are not propagating acoustic waves, but instead are evanescent in character since their frequency was lower than the associated acoustic cut-off value (see Sect. 3.1 for further details). Researchers further hypothesized that the magnetic fields, which were often synonymous with the observed oscillations, needed to be considered in order to accurately understand and model the wave dynamics (Michalitsanos 1973; Nakagawa 1973; Nakagawa et al. 1973; Stein and Leibacher 1974; Mein 1977, 1978, to name but a few examples). The field of magnetohydrodynamics (MHD) was introduced to effectively link the observed wave signatures to the underlying magnetic configurations, where the strong field strengths experienced in certain locations (e.g., field strengths that can approach approximately 6000 G in sunspot umbrae; Livingston et al. 2006; Okamoto and Sakurai 2018) produce wave modes that are highly modified from their purely acoustic counterparts.

The importance of the magnetic field in the studies of wave phenomena cannot be overestimated, since both the alignment of the embedded magnetic field, $B_0$, with the wavevector, k, and the ratio of the kinetic pressure, $p_0$, to the magnetic pressure, $B_{0}^{2}/2\mu _{0}$, play influential roles in the characteristics of any waves present (see the reviews by, e.g., Stein and Leibacher 1974; Bogdan 2000; Mathioudakis et al. 2013; Jess et al. 2015; Jess and Verth 2016). Commonly, the ratio of kinetic to magnetic pressures is referred to as the plasma-$\beta$, defined as,

$$\begin{aligned} \beta = \frac{2\mu _{0}p_{0}}{B_{0}^{2}} , \end{aligned}$$

(1)

where $\mu _{0}$ is the magnetic permeability of free space (Wentzel 1979; Edwin and Roberts 1983; Spruit and Roberts 1983). Crucially, by introducing the local hydrogen number density, $n_{\text{H}}$, the plasma-$\beta$ can be rewritten (in cgs units) in terms of the Boltzmann constant, $k_{B}$, and the temperature of the plasma, T, giving the relation,

$$\begin{aligned} \beta = \frac{8{\pi }n_{\text{H}}Tk_{B}}{B_{0}^{2}} . \end{aligned}$$

(2)

In the lower regions of the solar atmosphere, including the photosphere and chromosphere, temperatures are relatively low ($T \lesssim 15,000\,\text{K}$) when compared to the corona. This, combined with structures synonymous with the solar surface, including sunspots, pores, and magnetic bright points (MBPs; Berger et al. 1995; Sánchez Almeida et al. 2004; Ishikawa et al. 2007; Utz et al. 2009, 2010, 2013a, b; Keys et al. 2011, 2013, 2014), all of which possess strong magnetic field concentrations ($B_{0} \gtrsim 1000\,\text{G}$), presents wave conduits that are inherently ‘low-$\beta$’ (i.e., are dominated by magnetic pressure; $\beta \ll 1$). Gary (2001) has indicated how such structures (particularly for highly-magnetic sunspots) can maintain their low-$\beta$ status throughout the entire solar atmosphere, even as the magnetic fields begin to expand into the more volume-filling chromosphere (Gudiksen 2006; Beck et al. 2013b). Using non-linear force-free field (NLFFF; Wiegelmann 2008; Aschwanden 2016; Wiegelmann and Sakurai 2021) extrapolations, Aschwanden et al. (2016) and Grant et al. (2018) provided further evidence that sunspots can be best categorized as low-$\beta$ wave guides, spanning from the photosphere through to the outermost extremities of the corona. As can be seen from Eq. (2), the hydrogen number density ($n_{\text{H}}$) also plays a pivotal role in the precise local value of the plasma-$\beta$. As one moves higher in the solar atmosphere, a significant drop in the hydrogen number density is experienced (see, e.g., the sunspot model proposed by Avrett 1981), often with an associated scale-height on the order of 150–200 km (Alissandrakis 2020). As a result, the interplay between the number density and the expanding magnetic fields plays an important role in whether the environment is dominated by magnetic or plasma pressures.

Of course, not all regions of the Sun’s lower atmosphere are quite so straightforward. Weaker magnetic elements, including small-scale MBPs (Keys et al. 2020), are not able to sustain dominant magnetic pressures as their fields expand with atmospheric height. This results in the transition to a ‘high-$\beta$’ environment, where the plasma pressure dominates over the magnetic pressure (i.e., $\beta > 1$), which has been observed and modeled under a variety of highly magnetic conditions (e.g., Borrero and Ichimoto 2011; Jess et al. 2013; Bourdin 2017; Grant et al. 2018). This transition has important implications for the embedded waves, since the allowable modes become effected as the wave guide passes through the $\beta \sim 1$ equipartition layer. Here, waves are able to undergo mode conversion/transmission (Schunker and Cally 2006; Cally 2007; Hansen et al. 2016), which has the ability to change the properties and observable signatures of the oscillations. However, we note that under purely quiescent conditions (i.e., related to quiet Sun modeling and observations), the associated intergranular lanes (Lin and Rimmele 1999) and granules themselves (Lites et al. 2008) will already be within the high plasma-$\beta$ regime at photospheric heights.

Since the turn of the century, there has been a number of reviews published in the field of MHD waves manifesting in the outer solar atmosphere, including those linked to standing (van Doorsselaere et al. 2009; Wang 2011), quasi-periodic (Nakariakov et al. 2005), and propagating (de Moortel 2009; Zaqarashvili and Erdélyi 2009; Lin 2011) oscillations. Many of these review articles focus on the outermost regions of the solar atmosphere (i.e., the corona), or only address waves and oscillations isolated within a specific layer of the Sun’s atmosphere, e.g., the photosphere (Jess and Verth 2016) or the chromosphere (Jess et al. 2015; Verth and Jess 2016). As such, previous reviews have not focused on the coupling of MHD wave activity between the photosphere and chromosphere, which has only recently become possible due to the advancements made in multi-wavelength observations and data-driven MHD simulations. Here, in this review, we examine the current state-of-the-art in wave propagation, coupling, and damping/dissipation within the lower solar atmosphere, which comprises of both the photosphere and chromosphere, which are the focal points of next-generation ground-based telescopes, such as DKIST.

In addition, we would also like this review to be useful for early career researchers (PhD students and post-doctoral staff) who may not necessarily be familiar with all of the wave-based analysis techniques the solar physics community currently have at their disposal, let alone the wave-related literature currently in the published domain. As a result, we wish this review to deviate from traditional texts that focus on summarizing and potential follow-up interpretations of research findings. Instead, we will present traditional and state-of-the-art methods for detecting, isolating, and quantifying wave activity in the solar atmosphere. This is particularly important since modern data sequences acquired at cutting-edge observatories are providing us with incredible spatial, spectral, and temporal resolutions that require efficient and robust analyses tools in order to maximize the scientific return. Furthermore, we will highlight how the specific analysis methods employed often strongly influence the scientific results obtained, hence it is important to ensure that the techniques applied are fit for purpose. To demonstrate the observational improvements made over the last $\sim 50$ years we draw the readers attention to Figs. 2 and 3. Both Figs. 2 and 3 show sunspot structures captured using the best techniques available at that time. However, with advancements made in imaging (adaptive) optics, camera architectures, and post-processing algorithms, the drastic improvements are clear to see, with the high-quality data sequences shown in Fig. 3 highlighting the incredible observations of the Sun’s lower atmosphere we currently have at our disposal.

Fig. 3

Observations of a sunspot (top row) and a quiet-Sun region (middle row) in the lower solar atmosphere, sampled at three wavelength positions in the Ca ii 8542 Å spectral line from the 1 m Swedish Solar Telescope (SST). The wavelength positions, from left to right, correspond to $-900$ mÅ, $-300$ mÅ, and 0 mÅ from the line core, marked with vertical dashed lines in the bottom-right panel, where the average spectral line and all sampled positions are also depicted. The bottom-left panel illustrates a photospheric image sampled with a broadband filter (centered at 3950 Å; filter width $\approx 13.2$ Å). For better visibility, a small portion of the observed images are presented. All images are squared. Images courtesy of the Rosseland Centre for Solar Physics, University of Oslo

After the wave detection and analysis techniques have been identified, with their strengths/weaknesses defined, we will then take the opportunity to summarize recent theoretical and observational research focused on the generation, propagation, coupling, and dissipation of wave activity spanning the base of the photosphere, through to the upper echelons of the chromosphere that couples into the transition region and corona above. Naturally, addressing a key question in the research domain may subsequently pose two or three more, or pushing the boundaries of observational techniques and/or theoretical modeling tools may lead to ambiguities or caveats in the subsequent interpretations. This is not only to be expected, but should be embraced as a reminder of the era of rapid discovery we currently find ourselves in. The open questions we will pose not only highlight the challenges currently seeking solution with the dawn of next-generation ground-based and space-borne telescopes, but will also set the scene for research projects spanning decades to come.

2 Wave analysis tools

Identifying, extracting, quantifying, and understanding wave-related phenomena in astrophysical time series is a challenging endeavor. Signals that are captured by even the most modern charge-coupled devices (CCDs) and scientific complementary metal-oxide-semiconductor (sCMOS) detectors are accompanied by an assortment of instrumental and noise signals that act to mask the underlying periodic signatures. For example, the particle nature of the incident photons leads to Poisson-based shot noise, resulting in randomized intensity fluctuations about the time series mean (Terrell 1977; Delouille et al. 2008), which can reduce the clarity of wave-based signatures. Furthermore, instrumental and telescope effects, including temperature sensitivity and pointing stability, can lead to mixed signals either swamping the signatures of wave motion, or artificially creating false periodicities in the resulting data products. Hence, without large wave amplitudes it becomes a challenge to accurately constrain weak wave signals in even the most modern observational time series, especially once the wave fluctuations become comparable to the noise limitations of the data sequence. In the following sub-sections we will document an assortment of commonly available tools available to the solar physics community that can help quantify wave motion embedded in observational data.

2.1 Observations

In order for meaningful comparisons to be made from the techniques presented in Sect. 2, we will benchmark their suitability using two observed time series. We would like to highlight that the algorithms described and demonstrated below can be applied to any form of observational data product, including intensities, Doppler velocities, and spectral line-widths. As such, it is important to ensure that the input time series are scientifically calibrated before these wave analysis techniques are applied.

2.1.1 HARDcam: 2011 December 10

The Hydrogen-alpha Rapid Dynamics camera (HARDcam; Jess et al. 2012a) is an sCMOS instrument designed to acquire high-cadence H$\alpha$ images at the DST facility. The data captured by HARDcam on 2011 December 10 consists of 75 min (16:10–17:25 UT) of H$\alpha$ images, acquired through a narrowband 0.25 Å Zeiss filter, obtained at 20 frames per second. Active region NOAA 11366 was chosen as the target, which was located at heliocentric coordinates ($356''$, $305''$), or N17.9W22.5 in the more conventional heliographic coordinate system. A non-diffraction-limited imaging platescale of $0{\,}.{\!\!}{''}138$ per pixel was chosen to provide a field-of-view size equal to $71''\times 71''$. During the observing sequence, high-order adaptive optics (Rimmele 2004; Rimmele and Marino 2011) and speckle reconstruction algorithms (Wöger et al. 2008) were employed, providing a final cadence for the reconstructed images of 1.78 s. The dataset has previously been utilized in a host of scientific studies (Jess et al. 2013, 2016, 2017; Krishna Prasad et al. 2015; Albidah et al. 2021) due to the excellent seeing conditions experienced and the fact that the sunspot observed was highly circularly symmetric in its shape. A sample image from this observing campaign is shown in the right panel of Fig. 4, alongside a simultaneous continuum image captured by the Helioseismic and Magnetic Imager (HMI; Schou et al. 2012), onboard the Solar Dynamics Observatory (SDO; Pesnell et al. 2012).

Fig. 4

An SDO/HMI full-disk continuum image (left), with a red box highlighting the HARDcam field-of-view captured by the DST facility on 2011 December 10. An H$\alpha$ line core image of active region NOAA 11366, acquired by HARDcam at 16:10 UT, is displayed in the right panel. Axes represent heliocentric coordinates in arcseconds

In addition to the HARDcam data of this active region, we also accessed data from the Atmospheric Imaging Assembly (AIA; Lemen et al. 2012) onboard the SDO. Here, we obtained 1700 Å continuum (photospheric) images with a cadence of 24 s and spanning a 2.5 h duration. The imaging platescale is $0{\,}.{\!\!}{''}6$ per pixel, with a $350\times 350$ pixel$^{2}$ cut-out providing a $210''\times 210''$ field-of-view centered on the NOAA 11366 sunspot. The SDO/AIA images are used purely for the purposes of comparison to HARDcam information in Sect. 2.3.1.

2.1.2 SuFI: 2009 June 9

The Sunrise Filter Imager (SuFI; Gandorfer et al. 2011) onboard the Sunrise balloon-borne solar observatory (Solanki et al. 2010; Barthol et al. 2011; Berkefeld et al. 2011) sampled multiple photospheric and chromospheric heights, with a 1 m telescope, in distinct wavelength bands during its first and second flights in 2009 and 2013, respectively (Solanki et al. 2017). High quality, seeing-free time-series of images at 300 nm and 397 nm (Ca ii H) bands (approximately corresponding to the low photosphere and low chromosphere, respectively) were acquired by SuFI/Sunrise on 2009 June 9, between 01:32 UTC and 02:00 UTC, at a cadence of 12 sec after phase-diversity reconstructions (Hirzberger et al. 2010, 2011). The observations sampled a quiet region located at solar disk center with a field of view of $14''\times 40''$ and a spatial sampling of $0{\,}.{\!\!}{''}02$ per pixel. Figure 5 illustrates sub-field-of-view sample images in both bands (with an average height difference of $\approx 450$ km; Jafarzadeh et al. 2017d), along with magnetic-field strength map obtained from Stokes inversions of the Fe i 525.02 nm spectral line from the Sunrise Imaging Magnetograph eXperiment (IMaX; Martínez Pillet et al. 2011). A small magnetic bright point is also marked on all panels of Fig. 5 with a circle. Wave propagation between these two atmospheric layers in the small magnetic element is discussed in Sect. 2.4.1.

Fig. 5

A small region of an image acquired in 300 nm (left) and in Ca ii H spectral lines (middle) from SuFI/Sunrise, along with their corresponding line-of-sight magnetic fields from IMaX/Sunrise (right). The latter ranges between $-1654$ and 2194 G. The circle includes a small-scale magnetic feature whose oscillatory behavior is shown in Fig. 25

2.2 One-dimensional Fourier analysis

Traditionally, Fourier analysis (Fourier 1824) is used to decompose time series into a set of cosines and sines of varying amplitudes and phases in order to recreate the input lightcurve. Importantly, for Fourier analysis to accurately benchmark embedded wave motion, the input time series must be comprised of both linear and stationary signals. Here, a purely linear signal can be characterized by Gaussian behavior (i.e., fluctuations that obey a Gaussian distribution in the limit of large number statistics), while a stationary signal has a constant mean value and a variance that is independent of time (Tunnicliffe-Wilson 1989; Cheng et al. 2015). If non-linear signals are present, then the time series displays non-Gaussian behavior (Jess et al. 2019), i.e., it contains features that cannot be modeled by linear processes, including time-changing variances, asymmetric cycles, higher-moment structures, etc. In terms of wave studies, these features often manifest in solar observations in the form of sawtooth-shaped structures in time series synonymous with developing shock waves (Fleck and Schmitz 1993; Rouppe van der Voort et al. 2003; Vecchio et al. 2009; de la Cruz Rodríguez et al. 2013; Houston et al. 2018). Of course, it is possible to completely decompose non-linear signals using Fourier analysis, but the subsequent interpretation of the resulting amplitudes and phases is far from straightforward and needs to be treated with extreme caution (Lawton 1989).

On the other hand, non-stationary time series are notoriously difficult to predict and model (Tunnicliffe-Wilson 1989). A major challenge when applying Fourier techniques to non-stationary data is that the corresponding Fourier spectrum incorporates numerous additional harmonic components to replicate the inherent non-stationary behavior, which artificially spreads the true time series energy over an uncharacteristically wide frequency range (Terradas et al. 2004). Ideally, non-stationary data needs to be transformed into stationary data with a constant mean and variance that is independent of time. However, understanding the underlying systematic (acting according to a fixed plan or system; methodical) and stochastic (randomly determined; having a random probability distribution or pattern that may be analyzed statistically but may not be predicted precisely) processes is often very difficult (Adhikari and Agrawal 2013). In particular, differencing can mitigate stochastic (i.e., non-systematic) processes to produce a difference-stationary time series, while detrending can help remove deterministic trends (e.g., time-dependent changes), but may struggle to alleviate stochastic processes (Pwasong and Sathasivam 2015). Hence, it is often very difficult to ensure observational time series are truly linear and stationary.

The upper-left panel of Fig. 6 displays an intensity time series (lightcurve) that has been extracted from a penumbral pixel in the chromospheric HARDcam H$\alpha$ data. Here, the intensities have been normalized by the time-averaged quiescent H$\alpha$ intensity. It can be seen in the upper-left panel of Fig. 6 that in addition to sinusoidal wave-like signatures, there also appears to be a background trend (i.e., moving average) associated with the intensities. Through visual inspection, this background trend does not appear linear, thus requiring higher order polynomials to accurately model and remove. It must be remembered that very high order polynomials will likely begin to show fluctuations on timescales characteristic of the wave signatures wishing to be studied. Hence, it is important that the lowest order polynomial that best fits the data trends is chosen to avoid contaminating the embedded wave-like signatures with additional fluctuations arising from high-order polynomials. Importantly, the precise method applied to detrend the data can vary depending upon the signal being analyzed (e.g., Edmonds and Webb 1972; Edmonds 1972; Krijger et al. 2001; Rutten and Krijger 2003; de Wijn et al. 2005a, b). For example, some researchers choose to subtract the mean trend, while others prefer to divide by the fitted trend then subtract ‘1’ from the subsequent time series. Both approaches result in a more stationary time series with a mean value of ‘0’. However, subtracting the mean preserves the original unit of measurement and hence the original shape of the time series (albeit with modified numerical axes labels), while dividing by the mean provides a final unit that is independent of the original measurement and thus provides a method to more readily visualize fractional changes to the original time series. It must be noted that detrending processes, regardless of which approach is selected, can help remove deterministic trends (e.g., time-dependent changes), but often struggle to alleviate stochastic processes from the resulting time series.

Fig. 6

An H$\alpha$ line core intensity time series (upper left; solid black line) extracted from a penumbral location of the HARDcam data described in Sect. 2.1.1. The intensities shown have been normalized by the time-averaged H$\alpha$ intensity established in a quiet Sun region within the field-of-view. A dashed red line shows a third-order polynomial fitted to the lightcurve, which is designed to detrend the data to provide a stationary time series. The upper-right panel displays the resulting time series once the third-order polynomial trend line has been subtracted from the raw intensities (black line). The solid red line depicts an apodization filter designed to preserve 90% of the original lightcurve, but gradually reduce intensities to zero towards the edges of the time series to help alleviate any spurious signals in the resulting FFT. The lower panel reveals the final lightcurve that is ready for FFT analyses, which has been both detrended and apodized to help ensure the resulting Fourier power is accurately constrained. The horizontal dashed red lines signify the new mean value of the data, which is equal to zero due to the detrending employed

The dashed red line in the upper-left panel of Fig. 6 displays a third-order polynomial trend line fitted to the raw H$\alpha$ time series. The line of best fit is relatively low order, yet still manages to trace the global time-dependent trend. Subtracting the trend line from the raw intensity lightcurve provides fluctuations about a constant mean equal to zero (upper-right panel of Fig. 6), helping to ensure the resulting time series is stationary. It can be seen that wave-like signatures are present in the lightcurve, particularly towards the start of the observing sequence, where fluctuations on the order of $\approx 8{\%}$ of the continuum intensity are visible. However, it can also be seen from the right panel of Fig. 6 that between times of approximately 300–1300 s there still appears to be a local increase in the mean (albeit no change to the global mean, which remains zero). To suppress this local change in the mean, higher order polynomial trend lines could be fitted to the data, but it must be remembered that such fitting runs the risk of manipulating the true wave signal. Hence, for the purposes of this example, we will continue to employ third-order polynomial detrending, and make use of the time series shown in the upper-right panel of Fig. 6.

For data sequences that are already close to being stationary, one may question why the removal of such background trends is even necessary since the Fourier decomposition with naturally put the trend components into low-frequency bins. Of course, the quality and/or dynamics of the input time series will have major implications regarding what degree of polynomial is required to accurately transform the data into a stationary time series. However, from the perspective of wave investigations, non-zero means and/or slowly evolving backgrounds will inappropriately apply Fourier power across low frequencies, even though these are not directly wave related, which may inadvertently skew any subsequent frequency-integrated wave energy calculations performed. The sources of such non-stationary processes can be far-reaching, and include aspects related to structural evolution of the feature being examined, local observing conditions (e.g., changes in light levels for intensity measurements), and/or instrumental effects (e.g., thermal impacts that can lead to time-dependent variances in the measured quantities). As such, some of these sources (e.g., structural evolution) are dependent on the precise location being studied, while other sources (e.g., local changes in the light level incident on the telescope) are global effects that can be mapped and removed from the entire data sequence simultaneously. Hence, detrending the input time series helps to ensure that the resulting Fourier power is predominantly related to the embedded wave activity.

Another step commonly taken to ensure the reliability of subsequent Fourier analyses is to apply an apodization filter to the processed time series (Norton and Beer 1976). An Fourier transform assumes an infinite, periodically repeating sequence, hence leading to a looping behavior at the ends of the time series. Hence, an apodization filter is a function employed to smoothly bring a measured signal down to zero towards the extreme edges (i.e., beginning and end) of the time series, thus mitigating against sharp discontinuities that may arise in the form of false power (edge effect) signatures in the resulting power spectrum.

Typically, the apodization filter is governed by the percentage over which the user wishes to preserve the original time series. For example, a 90% apodization filter will preserve the middle 90% of the overall time series, with the initial and final 5% of the lightcurve being gradually tapered to zero (Dame and Martic 1987). There are many different forms of the apodization filter shape that can be utilized, including tapered cosines, boxcar, triangular, Gaussian, Lorentzian, and trapezoidal profiles, many of which are benchmarked using solar time series in Louis et al. (2015). A tapered cosine is the most common form of apodization filter found in solar physics literature (e.g., Hoekzema et al. 1998), and this is what we will employ here for the purposes of our example dataset. The upper-right panel of Fig. 6 reveals a 90% tapered cosine apodization filter overplotted on top of the detrended H$\alpha$ lightcurve. Multiplying this apodization filter by the lightcurve results in the final detrended and apodized time series shown in the bottom panel of Fig. 6, where the stationary nature of this processed signal is now more suitable for Fourier analyses. It is worth noting that following successful detrending of the input time series, the apodization percentage chosen can often be reduced, since the detrending process will suppress any discontinuities arising at the edges of the data sequence (i.e., helps to alleviate spectral leakage; Nuttall 1981). As such, the apodization percentage employed may be refined based on the ratio between the amplitude of the (primary) oscillatory signal and the magnitude of the noise present within that signal (i.e., linked to the inherent signal-to-noise ratio; Stoica and Moses 2005; Carlson and Crilly 2010).

Performing a fast Fourier transform (FFT; Cooley and Tukey 1965) of the detrended time series provides a Fourier amplitude spectrum, which can be displayed as a function of frequency. An FFT is a computationally more efficient version of the discrete Fourier transform (DFT; Grünbaum 1982), which only requires $N\log {N}$ operations to complete compared with the $N^{2}$ operations needed for the DFT, where N is the number of data points in the time series, which can be calculated by dividing the time series duration by the acquisition cadence. Following a Fourier transform of the input data, the number of (non-negative) frequency bins, $N_{f}$, can be computed by adding one to the number of samples (to account for the zeroth frequency representing the time series mean; Oppenheim and Schafer 2009), $N+1$, dividing the result by a factor of two, before rounding to the nearest integer. The Nyquist frequency is the highest constituent frequency of an input time series that can be evaluated at a given sampling rate (Grenander 1959), and is defined as $f_{\text{Ny}} = {\mathrm {sampling~rate}}/2 = 1/(2 \times {\text{cadence}})$. To evaluate the frequency resolution, $\varDelta {f}$, of an input time series, one must divide the Nyquist frequency by the number of non-zero frequency bins (i.e., the number of steps between the zeroth and Nyquist frequencies, N/2), providing,

$$\begin{aligned} \varDelta {f} = \frac{f_{\text{Ny}}}{N/2} = \frac{\frac{1}{2 \times {\text{cadence}}}}{\frac{{\mathrm {time~series~duration}}}{2 \times {\text{cadence}}}} = \frac{1}{{\mathrm {time~series~duration}}} . \end{aligned}$$

(3)

As a result, it is clear to see that the observing duration plays a pivotal role in the corresponding frequency resolution (see, e.g., Harvey 1985; Duvall et al. 1997; Gizon et al. 2017, for considerations in the helioseismology community). It is also important to note that the frequency bins remain equally spaced across the lowest (zeroth frequency or mean) to highest (Nyquist) frequency that is resolved in the corresponding Fourier spectrum. See Sect. 2.2.1 for a more detailed comparison between the terms involved in Fourier decomposition.

The HARDcam dataset utilized has a cadence of 1.78 s, which results in a Nyquist frequency of $f_{{\text {Ny}}}\approx 280\,\text{mHz}$ $\left( \frac{1}{2\times 1.78}\right)$. It is worth noting that only the positive frequencies are displayed in this review for ease of visualization. Following the application of Fourier techniques, both negative and positive frequencies, which are identical except for their sign, will be generated for the corresponding Fourier amplitudes. This is a consequence of the Euler relationship that allows sinusoidal wave signatures to be reconstructed from a set of positive and negative complex exponentials (Smith 2007). Since input time series are real valued (e.g., velocities, intensities, spectral line widths, magnetic field strengths, etc.) with no associated imaginary terms, then the Fourier amplitudes associated with the negative and positive frequencies will be identical. This results in the output Fourier transform being Hermitian symmetric (Napolitano 2020). As a result, the output Fourier amplitudes are often converted into a power spectrum (a measure of the square of the Fourier wave amplitude), or following normalization by the frequency resolution, into a power spectral density. This approach is summarized by Stull (1988), where the power spectral density, PSD, can be calculated as,

$$\begin{aligned} \text{PSD}(n) = \frac{2 \cdot |{\mathcal {F}}_{A}(n)|^{2}}{\varDelta f} = \frac{2 \cdot \left( \left[ {\mathcal {F}}_{\mathrm {real~part}}(n) \right] ^{2} + \left[ {\mathcal {F}}_{\mathrm {imaginary~part}}(n) \right] ^{2} \right) }{\varDelta f} . \end{aligned}$$

(4)

In Eq. (4), ${\mathcal {F}}_{A}(n)$ is the Fourier amplitude for any given positive frequency, n, while $\varDelta f$ is the corresponding frequency resolution of the Fourier transform (see definition above and further discussion points in Sect. 2.2.1). Note that the factor of ‘2’ is required due to the wrapping of identical Fourier power at negative frequencies into the positive domain. The normalization of the power spectrum by the frequency resolution is a best practice to ensure that the subsequent plots can be readily compared against other data sequences that may be acquired across shorter or longer observing intervals, hence affecting the intrinsic frequency resolution (see Sect. 2.2.1). As an example, the power spectral density of an input velocity time series, with units of km/s, will have the associated units of km$^{2}$/s$^{2}$/mHz (e.g., Stangalini et al. 2021b). The power spectral density for the detrended HARDcam H$\alpha$ time series is depicted in the lower-middle panel of Fig. 7. Here, the intensity time series is calibrated into normalized data number (DN) units, which are often equally labeled as ‘counts’ in the literature. Hence, the resulting power spectral density has units of DN$^{2}$/mHz.

Fig. 7

Taking the raw HARDcam H$\alpha$ lightcurve shown in the upper-left panel of Fig. 6, the upper row displays the resultant detrended time series utilizing linear (left), third-order polynomial (middle), and nineth-order polynomial (right) fits to the data. In each panel the dashed red line highlights the line of best fit, while the dashed blue line indicates the resultant data mean that is equal to zero following detrending. The lower row displays the corresponding Fourier power spectral densities for each of the linear (left), third-order polynomial (middle), and nineth-order polynomial detrended time series. Changes to the power spectral densities are particularly evident at low frequencies

An additional step often employed following the calculation of the PSD of an input time series is to remove the Fourier components associated with noise. It can be seen in the lower panels of Fig. 7 that there is a flattening of power towards higher frequencies, which is often due to the white noise that dominates the signal at those frequencies (Hoyng 1976; Krishna Prasad et al. 2017). Here, white noise is defined as fluctuations in a time series that give rise to equal Fourier power across all frequencies, hence giving rise to a flat PSD (Bendat and Piersol 2011). Often, if white noise is believed to be the dominant source of noise in the data (i.e., the signal is well above the detector background noise, hence providing sufficient photon statistics so that photon noise is the dominant source of fluctuations), then its PSD can be estimated by applying Eq. (4) to a random light curve generated following a Poisson distribution, with an amplitude equivalent to the square root of the mean intensity of the time series (Fossum and Carlsson 2005a; Lawrence et al. 2011). Subtraction of the background noise is necessary when deriving, for example, the total power of an oscillation isolated in a specific frequency window (Vaseghi 2008). Other types of noise exist that have discernible power-law slopes associated with their PSDs as a function of frequency. For example, while white noise has a flat power-law slope, pink and red noise display 1/f and $1/f^{2}$ power-law slopes, respectively, resulting in larger amplitudes at lower frequencies (Kolotkov et al. 2016; Strekalova et al. 2018). The specific dominant noise profile must be understood before it is subtracted from the relevant data PSDs.

As a result of the detrending employed in Fig. 6, the absolute Fourier wave amplitude related to a frequency of 0 Hz (i.e., representing the time series mean; upper panel of Fig. 8) is very low; some 4 orders-of-magnitude lower than the power associated with white noise signatures at high frequencies. Of course, if the processed time series mean is exactly zero, then the Fourier wave amplitude at 0 Hz should also be zero. In the case of Fig. 8, the detrended time series does have a zero mean. However, because the time series is not antisymmetric about the central time value, it means that the application of the tapered cosine apodization function results in a very small shift in the time series mean away from the zero value. As a result, the subsequent Fourier amplitudes are fractionally (e.g., at the $10^{-8}$ level for the upper panel of Fig. 8) above the zero point. Once the processes of detrending and apodization are complete, it is possible to re-calculate the time series mean and subtract this value to ensure the processed mean remains zero before the application of Fourier analyses. However, for the purposes of Figs. 7 and 8, this additional mean subtraction has not been performed to better highlight this potential artifact at the lowest temporal frequencies.

Fig. 8

Fourier power spectrum of the HARDcam H$\alpha$ detrended lightcurve shown in the lower panel of Fig. 6 (top). For the purposes of wave filtering, a step function is shown on the Fourier spectrum using a dashed red line (middle left), where the step function equals unity between frequencies spanning 3.7–5.7 mHz (i.e., $4.7\pm 1.0$ mHz). Multiplying the Fourier power spectrum by this step function results in isolated power features, which are displayed in the middle-right panel. Alternatively, a Gaussian function centered on 4.7 mHz, with a FWHM of 2.0 mHz, is overplotted on top of the Fourier power spectrum using a red line in the lower-left panel. Multiplying the power spectrum by the Gaussian function results in similar isolated power features, shown in the lower-right panel, but with greater apodization of edge frequencies to help reduce aliasing upon reconstruction of the filtered time series

Note that Fig. 8 does not have the frequency axis displayed on a log-scale in order to reveal the 0 Hz component. As such, the upper frequency range is truncated to $\approx 28$ Hz to better reveal the signatures present at the lower frequencies synonymous with wave activity in the solar atmosphere. The suppression of Fourier wave amplitudes at the lowest frequencies suggests that the third-order polynomial trend line fitted to the raw H$\alpha$ intensities is useful at removing global trends in the visible time series. However, as discussed above, care must be taken when selecting the polynomial order to ensure that the line of best fit does not interfere with the real wave signatures present in the original lightcurve. To show the subtle, yet important impacts of choosing a suitable trend line, Fig. 7 displays the resultant detrended time series of the original HARDcam H$\alpha$ lightcurve for three different detrending methods, e.g., the subtraction of a linear, a third-order polynomial, and a nineth-order polynomial line of best fit. It can be seen from the upper panels of Fig. 7 that the resultant (detrended) lightcurves have different perturbations away from the new data mean of zero. This translates into different Fourier signatures in the corresponding power spectral densities (lower panels of Fig. 7), which are most apparent at the lowest frequencies (e.g., $<3$ mHz). Therefore, it is clear that care must be taken when selecting the chosen order of the line of best fit so that it doesn’t artificially suppress true wave signatures that reside in the time series. It can be seen in the lower-middle panel of Fig. 7 that the largest Fourier power signal is at a frequency of $\approx 4.7$ mHz, corresponding to a periodicity of $\approx 210$ s, which is consistent with previous studies of chromospheric wave activity in the vicinity of sunspots (e.g., Felipe et al. 2010; Jess et al. 2013; López Ariste et al. 2016, to name but a few examples).

2.2.1 Common misconceptions involving Fourier space

Translating a time series into the frequency-dependent domain through the application of a Fourier transform is a powerful diagnostic tool for analyzing the frequency content of (stationary) time series. However, when translating between the temporal and frequency domains it becomes easy to overlook the importance of the sampling cadence and the time series duration in the corresponding frequency axis. For example, one common misunderstanding is the belief that increasing the sampling rate of the data (e.g., increasing the frame rate of the observations from 10 frames per second to 100 frames per second) will improve the subsequent frequency resolution of the corresponding Fourier transform. Unfortunately, this is not the case, since increasing the frame rate raises the Nyquist frequency (highest frequency component that can be evaluated), but does not affect the frequency resolution of the Fourier transform. Instead, to improve the frequency resolution one must obtain a longer-duration time series or employ ‘padding’ of the utilized lightcurve to increase the number of data points spanning the frequency domain (Lyons 1996).

To put these aspects into better context, we will outline a worked example that conveys the importance of both time series cadence and duration. Let us consider two complementary data sequences, one from the Atmospheric Imaging Assembly (AIA; Lemen et al. 2012) onboard the SDO spacecraft, and one from the 4 m ground-based Daniel K. Inouye Solar Telescope (DKIST; Tritschler et al. 2016; Rimmele et al. 2020; Rast et al. 2021). Researchers undertaking a multi-wavelength investigation of wave activity in the solar atmosphere may choose to employ these types of complementary observations in order to address their science objectives. Here, the AIA/SDO observations consist of 3 h (10,800 s) of 304 Å images taken at a cadence of 12.0 s, while the DKIST observations comprise of 1 h (3600 s) of H$\alpha$ observations taken by the Visual Broadband Imager (VBI; Wöger 2014) at a cadence of 3.2 s.

The number of samples, N, for each of the time series can be calculated as $N_{\text{AIA}} = 10{,}800 / 12.0 = 900$ and $N_{\text{VBI}} = 3600 / 3.2 = 1125$. Therefore, it is clear that even though the AIA/SDO observations are obtained over a longer time duration, the higher cadence of the VBI/DKIST observations results in more samples associated with that data sequence. The number of frequency bins, $N_{f}$, can also be computed as $N_{f({\text{AIA}})} = (900+1) / 2 = 451$, while ${N_{f({\text{VBI}})} = (1125+1) / 2 = 563}$. Hence, the frequency axes of the corresponding Fourier transforms will be comprised of 451 and 563 positive real frequencies (i.e., $\ge 0$ Hz) for the AIA/SDO and VBI/DKIST data, respectively. The increased number of frequency bins for the higher cadence VBI/DKIST observations sometimes leads to the belief that this provides a higher frequency resolution. However, we have not yet considered the effect of the image cadence on the corresponding frequency axes.

In the case of the AIA/SDO and VBI/DKIST observations introduced above, the corresponding Nyquist frequencies can be computed as $f_{\mathrm {Ny(AIA)}} = 1/(2 \times 12.0) \approx 42$ mHz and ${f_{\mathrm {Ny(VBI)}} = 1/(2 \times 3.2) \approx 156}$ mHz, respectively. As a result, it should become clear that while the VBI/DKIST observations result in a larger number of corresponding frequency bins (i.e., $N_{f({\text{VBI}})} > N_{f({\text{AIA}})}$), these frequency bins are required to cover a larger frequency interval up to the calculated Nyquist value. Subsequently, for the case of the AIA/SDO and VBI/DKIST observations, the corresponding frequency resolutions can be calculated as ${\varDelta {f}_{\text{AIA}} = 1/10{,}800 = 0.0926}$ mHz and ${\varDelta {f}_{\text{VBI}} = 1/3600 = 0.2778}$ mHz, respectively. Note that while the frequency resolution is constant, the same cannot be said for the period resolution due to the reciprocal nature between these variables. For example, at a frequency of 3.3 mHz ($\approx 5$ min oscillation), the period resolution for VBI/DKIST is $\approx 25$ s (i.e., $\approx 303\pm 25$ s), while for AIA/SDO the period resolution is $\approx 8$ s (i.e., $\approx 303\pm 8$ s). Similarly, at a frequency of 5.6 mHz ($\approx 3$ min oscillation), the period resolutions for VBI/DKIST and AIA/SDO are $\approx 9$ s (i.e., $\approx 180\pm 9$ s) and $\approx 3$ s (i.e., $\approx 180\pm 3$ s), respectively.

Figure 9 depicts the Fourier frequencies (left panel), and their corresponding periodicities (right panel), as a function of the derived frequency bin. It can be seen from the left panel of Fig. 9 that the AIA/SDO observations produce a lower number of frequency bins (i.e., a result of less samples, $N_{\text{AIA}} < N_{\text{VBI}}$), alongside a smaller peak frequency value (i.e., a lower Nyquist frequency, ${f_{\mathrm {Ny(AIA)}}} < {f_{\mathrm {Ny(VBI)}}}$, caused by the lower temporal cadence). However, as a result of the longer duration observing sequence for the AIA/SDO time series (i.e., 3 h for AIA/SDO versus 1 h for VBI/DKIST), the resulting frequency resolution is better (i.e., ${\varDelta {f}_{\text{AIA}}} < {\varDelta {f}_{\text{VBI}}}$), allowing more precise frequency-dependent phenomena to be uncovered in the AIA/SDO observations. Of course, due to the AIA/SDO cadence being longer than that of VBI/DKIST (i.e., 12.0 s for AIA/SDO versus 3.2 s for VBI/DKIST), this results in the inability to examine the fastest wave fluctuations, which can be seen more clearly in the right panel of Fig. 9, whereby the VBI/DKIST observations are able to reach lower periodicities when compared to the complementary AIA/SDO data sequence. The above scenario is designed to highlight the important interplay between observing cadences and durations with regards to the quantitative parameters achievable through the application of Fourier transforms. For example, if obtaining the highest possible frequency resolution is of paramount importance to segregate closely matched wave frequencies, then it is the overall duration of the time series (not the observing cadence) that facilitates the necessary frequency resolution.

Fig. 9

The frequencies (left panel) and corresponding periodicities (right panel) that can be measured through the application of Fourier analysis to an input time series. Here, the solid blue lines depict AIA/SDO observations spanning a 3 h duration and acquired with a temporal cadence of 12.0 s, while the solid red lines highlight VBI/DKIST observations spanning a 1 h window and acquired with a temporal cadence of 3.2 s. It can be seen that both the cadence and observing duration play pivotal roles in the resulting frequencies/periodicities achievable, with the longer duration AIA/SDO observations providing a better frequency resolution, $\varDelta {f}$, while the higher cadence VBI/DKIST data results in a better Nyquist frequency that allows more rapid wave fluctuations to be studied. In the left and right panels, the dashed blue and red lines depict the Nyquist frequencies and corresponding periodicities for the AIA/SDO and VBI/DKIST data sequences, respectively (see text for more information)

Another important aspect to keep in mind is that the Fourier spectrum is only an estimate of the real power spectrum of the studied process. The finite-duration time series, noise, and distortions due to the intrinsic covariance within each frequency bin may lead to spurious peaks in the spectrum, which could be wrongly interpreted as real oscillations. As a result, one may believe that by considering longer time series the covariance of each frequency bin will reduce, but this is not true since the bin width itself becomes narrower. One way forward is to divide the time series into different segments and average the resulting Fourier spectra calculated from each sub-division—the so-called Welch method (Welch 1967), at the cost of reducing the resolution of frequencies explored. However, data from ground-based observatories are generally limited to 1–2 h each day, and it is not always possible to obtain such long time series. Therefore, special attention must be paid when interpreting the results.

It is also possible to artificially increase the duration of the input time series through the process known as ‘padding’ (Ransom et al. 2002), which has been employed across a wide range of solar studies incorporating the photosphere, chromosphere, and corona (e.g., Ballester et al. 2002; Auchère et al. 2016; Hartlep and Zhao 2021; Jafarzadeh et al. 2021). Here, the beginning and/or end of the input data sequence is appended with a large number of data points with values equal to the mean of the overall time series. The padding adds no additional power to the data, but it acts to increase the fine-scale structure present in the corresponding Fourier transform since the overall duration of the data has been artificially increased. Note that padding with the data mean is preferable to padding with zeros since this alleviates the introduction of low-frequency power into the subsequent Fourier transform. Of course, if the input time series had previously been detrended (see Sect. 2.2) so that the resulting mean of the data is zero, then zero-padding and padding with the time series mean are equivalent.

Note that the process of padding is often perceived to increase the usable Fourier frequency resolution of the dataset, which is unfortunately incorrect. The use of padded time series acts to reveal small-scale structure in the output Fourier transform, but as it does not add any real signal to the input data sequence, the frequency resolution remains governed by the original time series characteristics (Eriksson 1998). As such, padding cannot recover and/or recreate any missing information in the original data sequence. This effect can be visualized in Fig. 10. Here, a resultant wave consisting of two sinusoids with normalized frequencies 0.075 and 0.125 of the sampling frequency is cropped to 32 and 64 data points in length. Figure 10a shows the corresponding power spectral density (PSD) following Fourier transformation on both the raw 32 data samples array (solid black line with circular data points) and the original 32 data point array that has been padded to a total of 64 data points (dashed black line with crosses). In addition, Fig. 10b shows another PSD for the data array containing 64 input samples (solid black line with circular data points), alongside the same PSD for the original 32 data point array that has been padded to a total of 64 data points (dashed black line with crosses; same as Fig. 10a). From Fig. 10a it can be seen that while the padding increases the number of data points along the frequency axis (and therefore creates some additional small-scale fluctuations in the resulting PSD), it does not increase the frequency resolution to a value needed to accurately identify the two sinusoidal components. This is even more apparent in Fig. 10b, where the Fourier transform of the time series containing 64 data points now contains sufficient information and frequency resolution to begin to segregate the two sinusoidal components. The padded array (32 data points plus 32 padded samples) contains the same number of elements along the frequency axis, but does not increase the frequency resolution to allow the quantification of the two embedded wave frequencies. The use of padding is often employed to decrease the computational time. Indeed, FFT algorithms work more efficiently if the number of samples is an integer power of 2.

Fig. 10

Image reproduced with permission from Eriksson (1998)

Panels revealing the effect of padding an input time series on the resulting Fourier transform. For this example, two sinusoids are superimposed with normalized frequencies equal to 0.075 and 0.125 of the sampling frequency. Panels a, b show the resulting power spectral densities (PSDs) following the Fourier transforms of 32 input data points (solid black line with circular data points; left) and 64 input data points (solid black line with circular data points; right), respectively. In both panels, the dashed black lines with crosses represent the Fourier transforms of 32 input data points that have been padded to a total of 64 data points. It can be seen that the increased number of data points associated with the padded array results in more samples along the frequency axis, but this does not improve the frequency resolution to the level consistent with supplying 64 genuine input samples (solid black line in the right panel).

Of course, while data padding strictly does not add usable information into the original time series, it can be utilized to provide better visual segregation of closely spaced frequencies. To show an example of this application, Fig. 11 displays the effects of padding and time series duration in a similar format to Fig. 10. In Fig. 11, the upper-left panel shows an intensity time series that is created from the superposition of two closely spaced frequencies, here 5.0 mHz and 5.4 mHz. The resultant time series is $\approx 3275$ s ($\sim 55$ min) long, and constructed with a cadence of 3.2 s to remain consistent with the VBI/DKIST examples shown earlier in this section. The absolute extent of this 3275 s time series is bounded in the upper-left panel of Fig. 11 by the shaded orange background. In order to pad this lightcurve, a new time series is constructed that has twice as many data points in length, making the time series duration now $\approx 6550$ s ($\sim 110$ min). The original $\approx 3275$ s lightcurve is placed in the middle of the new (expanded) array, thus providing zero-padding at the start and end of the time series. The corresponding power spectral densities (PSDs) for both the original and padded time series are shown in the lower-left panel of Fig. 11 using black and red lines, respectively. Note that the frequency axis is cropped to the range of 1–10 mHz for better visual clarity. It is clear that the original input time series creates a broad spectral peak at $\approx 5$ mHz, but the individual 5.0 mHz and 5.4 mHz components are not visible in the corresponding PSD (solid black line in the lower-left panel of Fig. 11). On the other hand, the PSD from the padded array (solid red line in the lower-left panel of Fig. 11) does show a double peak corresponding to the 5.0 mHz and 5.4 mHz wave components, highlighting how such padding techniques can help segregate multi-frequency wave signatures.

Fig. 11

Upper left: Inside the shaded orange region is a synthetic lightcurve created from the superposition of 5.0 mHz and 5.4 mHz waves, which are generated with a 3.2 s cadence (i.e., from VBI/DKIST) over a duration of $\approx 3275$ s. This time series is zero-padded into a $\approx 6550$ s array, which is displayed in its entirety in the upper-left panel using a solid black line. Upper right: The same resultant waveform created from the superposition of 5.0 mHz and 5.4 mHz waves, only now generated for the full $\approx 6550$ s time series duration (i.e., no zero-padding required). Lower left: The power spectral density (PSD) of the original (un-padded) lightcurve is shown using a solid black line, while the solid red line reveals the PSD of the full zero-padded time series. It is clear that the padded array offers better visual segregation of the two embedded wave frequencies. Lower right: The PSDs for both the full $\approx 6550$ s time series (solid black line) and the zero-padded original lightcurve (solid red line; same as that depicted in the lower-left panel). It can be seen that while the padded array provides some segregation of the 5.0 mHz and 5.4 mHz wave components, there is no better substitute at achieving high frequency resolution than obtaining long-duration observing sequences. Note that both PSD panels have the frequency axis truncated between 1 and 10 mHz for better visual clarity

Of course, padding cannot be considered a universal substitute for a longer duration data sequence. The upper-right panel of Fig. 11 shows the same input wave frequencies (5.0 mHz and 5.4 mHz), only with the resultant wave now present throughout the full $\sim 110$ min time sequence. Here, the beat pattern created by the superposition of two closely spaced frequencies can be readily seen, which is a physical manifestion of wave interactions also studied in high-resolution observations of the lower solar atmosphere (e.g., Krishna Prasad et al. 2015). The resulting PSD of the full-duration time series is depicted in the lower-right panel of Fig. 11 using a solid black line. For comparison, the PSD constructed from the padded original lightcurve is also overplotted using a solid red line (same as shown using a solid red line in the lower-left panel of Fig. 11). It is clearly seen that the presence of the wave signal across the full time series provides the most prominent segregation of the 5.0 mHz and 5.4 mHz spectral peaks. While these peaks are also visible in the padded PSD (solid red line), they are less well defined, hence reiterating that while time series padding can help provide better isolation of closely spaced frequencies, there is no better candidate for high frequency resolution than long duration observing sequences.

On the other hand, if rapidly fluctuating waveforms are wanting to be studied, then achieving a high Nyquist frequency is necessary to achieve these objectives, which the duration of the observing sequence is unable to assist with. Hence, it is important to tailor the observing strategy to ensure the frequency requirements are met. This, of course, can present challenges for particular facilities. For example, if a frequency resolution of $\varDelta {f} \approx 35~\mu \text{Hz}$ is required (e.g., to probe the sub-second timescales of physical processes affecting frequency distributions in the aftermath of solar flares; Wiśniewska et al. 2019), this would require an observing duration of approximately 8 continuous hours, which may not be feasible from ground-based observatories that are impacted by variable weather and seeing conditions. Similarly, while space-borne satellites may be unaffected by weather and atmospheric seeing, these facilities may not possess a sufficiently large telescope aperture to probe the wave characteristics of small-scale magnetic elements (e.g., Chitta et al. 2012b; Van Kooten and Cranmer 2017; Keys et al. 2018) and naturally have reduced onboard storage and/or telemetry restrictions, thus creating difficulties obtaining 8 continuous hours of observations at maximum acquisition cadences. Hence, complementary data products, including ground-based observations at high cadence and overlapping space-borne data acquired over long time durations, are often a good compromise to help provide the frequency characteristics necessary to achieve the science goals. Of course, next-generation satellite facilities, including the recently commissioned Solar Orbiter (Müller et al. 2013, 2020) and the upcoming Solar-C (Shimizu et al. 2020) missions, will provide breakthrough technological advancements to enable longer duration and higher cadence observations of the lower solar atmosphere than previously obtained from space. Another alternative to achieve both long time-series and high-cadence observations is the use of balloon-borne observatories, including the Sunrise (Bello González et al. 2010b) and Flare Genesis (Murphy et al. 1996; Bernasconi et al. 2000) experiments, where the data are stored in onboard discs. Such missions, however, have their own challenges and are limited to only a couple of days of observations during each flight.

2.2.2 Calculating confidence levels

After displaying Fourier spectra, it is often difficult to pinpoint exactly what features are significant, and what power spikes may be the result of noise and/or spurious signals contained within the input time series. A robust method of determining the confidence level of individual power peaks is to compare the Fourier transform of the input time series with the Fourier transform of a large number (often exceeding 1000) of randomized lightcurves based on the original values (i.e., ensuring an identical distribution of intensities throughout the new randomized time series; O’Shea et al. 2001). Following the randomization and computation of FFTs of the new time series, the probability, p, of randomized fluctuations being able to reproduce a given Fourier power peak in the original spectrum can be calculated. To do this, the Fourier power at each frequency element is compared to the power value calculated for the original time series, with the proportion of permutations giving a Fourier power value greater than, or equal to, the power associated with the original time series providing an estimate of the probability, p. Here, a small value of p suggests that the original lightcurve contains real oscillatory phenomena, while a large value of p indicates that there is little (or no) real periodicities contained within the data (Banerjee et al. 2001; O’Shea et al. 2001). Indeed, it is worth bearing in mind that probability values of $p=0.5$ are consistent with noise fluctuations (i.e., the variance of a binomial distribution is greatest at $p=0.5$; Lyden et al. 2019), hence why the identification of real oscillations requires small values of p.

Following the calculation of the probability, p, the value can be reversed to provide a percentage probability that the detected oscillatory phenomenon is real, through the relationship,

$$\begin{aligned} p_{\text {real}} = (1 - p) \times 100 . \end{aligned}$$

(5)

Here, $p_{\text {real}}=100\%$ would suggest that the wave motion present in the original time series is real, since no (i.e., $p=0$) randomized time series provided similar (or greater) Fourier power. Contrarily, $p_{\text {real}}=0\%$ would indicate a real (i.e., statistically significant) power deficit at that frequency, since all (i.e., $p=1$) randomized time series provided higher Fourier power at that specific frequency. Finally, a value of $p_{\text {real}} = 50\%$ would indicate that the power peak is not due to actual oscillatory motions. A similar approach is to calculate the means and standard deviations of the Fourier power values for each independent frequency corresponding to the randomized time series. This provides a direct estimate of whether the original measured Fourier power is within some number of standard deviations of the mean randomized-data power density. As a result, probability estimations of the detected Fourier peaks can be estimated providing the variances and means of the randomized Fourier power values are independent (i.e., follow a normal distribution; Bell et al. 2018).

If a large number ($\ge 1000$) of randomized permutations are employed, then the fluctuation probabilities will tend to Gaussian statistics (Linnell Nemec and Nemec 1985; Delouille et al. 2008; Jess et al. 2019). In this case, the confidence level can be obtained using a standardized Gaussian distribution. For many solar applications (e.g., McAteer et al. 2002b, 2003; Andic 2008; Bello González et al. 2009; Stangalini et al. 2012; Dorotovič et al. 2014; Freij et al. 2016; Jafarzadeh et al. 2017d, to name but a few examples), a confidence level of 95% is typically employed as a threshold for reliable wave detection. In this case, $99\% \le p_{\text {real}} \le 100\%$ (or $0.00 \le p \le 0.01$) is required to satisfy the desired 95% confidence level.

To demonstrate a worked example, we utilize the HARDcam H$\alpha$ time series shown in the left panel of Fig. 6, which consists of 2528 individual time steps. This, combined with 1000 randomized permutations of the lightcurve, provides 1000 FFTs with 1000 different measures in each frequency bin; more than sufficient to allow the accurate use of Gaussian number statistics (Montgomery and Runger 2003). For each randomization, the resulting Fourier spectrum is compared to that depicted in the upper panel of Fig. 8, with the resulting percentage probabilities, $p_{\text {real}}$, calculated according to Eq. (5) for each of the temporal frequencies. The original Fourier power spectrum, along with the percentage probabilities for each corresponding frequency, are shown in the left panel of Fig. 12. It can be seen that the largest power signal at $\approx 4.7$ mHz ($\approx 210$ s) has a high probability, suggesting that this is a detection of a real oscillation. Furthermore, the neighboring frequencies also have probabilities above 99%, further strengthening the interpretation that wave motion is present in the input time series. It should be noted that with potentially thousands of frequency bins in the high-frequency regime of an FFT, having some fraction of points that exceed a 95% (or even 99%) confidence interval is to be expected. Therefore, many investigations also demand some degree of coherency in the frequency and/or spatial distributions to better verify the presence of a real wave signal (similar to the methods described by Durrant and Nesis 1982; Di Matteo and Villante 2018). To better highlight which frequencies demonstrate confidence levels exceeding 95%, the right panel of Fig. 12 overplots (using bold red crosses) those frequencies containing percentage probabilities in excess of 99%.

Fig. 12

The full frequency extent of the Fourier power spectral densities shown in the lower-middle panel of Fig. 7, displayed using a log–log scale for better visual clarity (left panel). Overplotted using a solid red line are the percentage probabilities, $p_{\text {real}}$, computed over 1000 randomized permutations of the input lightcurve. Here, any frequencies with $p_{\text {real}} \ge 99\%$ correspond to a statistical confidence level in excess of 95%. The same Fourier power spectral density is shown in the right panel, only now with red crossed symbols highlighting the locations where the Fourier power provides confidence levels greater than 95%

2.2.3 Lomb-scargle techniques

A requirement for the implementation of traditional Fourier-based analyses is that the input time series is regularly and evenly sampled. This means that each data point of the lightcurve used should be obtained using the same exposure time, and subsequent time steps should be acquired with a strict, uniform cadence. Many ground-based and space-borne instruments employ digital synchronization triggers for their camera systems that can bring timing uncertainties down to the order of $10^{-6}$ s (Jess et al. 2010b), which is often necessary in high-precision polarimetric studies (Kootz 2018). This helps to ensure the output measurements are sufficiently sampled for the application of Fourier techniques.

However, often it is not possible to obtain time series with strict and even temporal sampling. For example, raster scans using slit-based spectrographs can lead to irregularly sampled observations due to the physical times required to move the spectral slit.¹ Also, some observing strategies interrupt regularly sampled data series for the measurement of Stokes I/Q/U/V signals every few minutes, hence introducing data gaps during these times (e.g., Samanta et al. 2016). Furthermore, hardware requiring multiple clocks to control components of the same instrument (e.g., the mission data processor and the polarization modulator unit on board the Hinode spacecraft; Kosugi et al. 2007) may have a tendency to drift away from one another, hence effecting the regularity of long-duration data sequences (Sekii et al. 2007). In addition, some facilities including the Atacama Large Millimeter/submillimeter Array (ALMA; Wootten and Thompson 2009; Wedemeyer et al. 2016) require routine calibrations that must be performed approximately every 10 min (with each calibration taking $\sim 2.5$ min; Wedemeyer et al. 2020), hence introducing gaps in the final time series (Jafarzadeh et al. 2021). Finally, in the case of ground-based observations, a period of reduced seeing quality or the passing of a localized cloud will result in a number of compromised science frames, which require removal and subsequent interpolation (Krishna Prasad et al. 2016).

If the effect of data sampling irregularities is not believed to be significant (i.e., is a fraction of the wave periodicities expected), then it is common to interpolate the observations back on to a constant cadence grid (e.g., Jess et al. 2012c; Kontogiannis et al. 2016). Of course, how the data points are interpolated (e.g., linear or cubic fitting) may effect the final product, and as a result, care should be taken when interpolating time series so that artificial periodicities are not introduced to the data through inappropriate interpolation. This is particularly important when the data sequence requires subsequent processing, e.g., taking the derivative of a velocity time series to determine the acceleration characteristics of the plasma. Under these circumstances, inappropriate interpolation of the velocity information may have drastic implications for the derived acceleration data. For this form of analysis, the use of 3-point Lagrangian interpolation is often recommended to ensure the endpoints of the time series remain unaffected due to the use of error propagation formulae (Veronig et al. 2008). However, in the case for very low cadence data, 3-point Lagrangian interpolation may become untrustworthy due to the large temporal separation between successive time steps (Byrne et al. 2013). For these cases, a Savtizky–Golay (Savitzky and Golay 1964) smoothing filter can help alleviate sharp (and often misleading) kinematic values (Byrne 2015).

If interpolation of missing data points and subsequent Fourier analyses is not believed to be suitable, then Lomb–Scargle techniques (Lomb 1976; Scargle 1982) can be implemented. As overviewed by Zechmeister and Kürster (2009), the Lomb–Scargle algorithms are useful for characterizing periodicities present in unevenly sampled data products. Often, least-squares minimization processes assume that the data to be fitted are normally distributed (Barret and Vaughan 2012), which may be untrue since the spectrum of a linear, stationary stochastic process naturally follows a $\chi _{2}^{2}$ distribution (Groth 1975; Papadakis and Lawrence 1993). However, a benefit of implementing the Lomb–Scargle algorithms is that the noise at each individual frequency can be represented by a $\chi ^{2}$ distribution, which is equivalent to a spectrum being reliably derived from more simplistic least-squares analysis techniques (VanderPlas 2018).

Crucially, Lomb–Scargle techniques differ from conventional Fourier analyses by the way in which the corresponding spectra are computed. While Fourier-based algorithms compute the power spectrum by taking dot products of the input time series with pairs of sine- and cosine-based waveforms, Lomb–Scargle techniques attempt to first calculate a delay timescale so that the sinusoidal pairs are mutually orthogonal at discrete sample steps, hence providing better power estimates at each frequency without the strict requirement of evenly sampled data (Press et al. 2007). In the field of solar physics, Lomb–Scargle techniques tend to be more commonplace in investigations of long-duration periodicities spanning days to months (i.e., often coupled to the solar cycle; Ni et al. 2012; Deng et al. 2017), although they can be used effectively in shorter duration observations where interpolation is deemed inappropriate (e.g., Maurya et al. 2013).

2.2.4 One-dimensional Fourier filtering

Often, it is helpful to filter the time series in order to isolate specific wave signatures across a particular range of frequencies. This is useful for a variety of studies, including the identification of beat frequencies (Krishna Prasad et al. 2015), the more reliable measurement of phase variations between different wavelengths/filters (Krishna Prasad et al. 2017), and in the identification of various wave modes co-existing within single magnetic structures (Keys et al. 2018). From examination of the upper panel of Figs. 8 and 12, it is clear that the frequency associated with peak Fourier power is $\approx 4.7$ mHz, and is accompanied by high confidence levels exceeding 95%.

If we wish to reconstruct a filtered time series centered on this dominant frequency, then we have a number of options available. The dashed red line in the middle-left panel of Fig. 8 depicts a step function frequency range of $4.7\pm 1.0$ mHz, whereby the filter is assigned values of ‘1’ and ‘0’ for frequencies inside and outside, respectively, this chosen frequency range. Multiplying the Fourier power spectrum by this step function frequency filter results in the preserved power elements that are shown in the middle-right panel of Fig. 8, which can be passed through an inverse FFT to create a Fourier filtered time series in the range of $4.7\pm 1.0$ mHz. However, by employing a step function frequency filter, there is a sharp and distinct transition between elevated power signals and frequencies with zero Fourier power. This abrupt transition can create aliasing artifacts in the reconstructed time series (Gobbi et al. 2006). Alternatively, to help mitigate against aliasing (i.e., sharp Fourier power transitions at the boundaries of the chosen frequency range), the Fourier power spectrum can be multiplied by a filter that peaks at the desired frequency, before gradually reducing in transmission towards the edges of the frequency range. An example of such a smoothly varying filter is documented in the lower panels of Fig. 8, where a Gaussian centered at 4.7 mHz, with a full-width at half-maximum (FWHM) of 2 mHz, is overplotted on top of the Fourier spectrum using a solid red line, which can be multiplied by the original Fourier spectrum to gradually decrease the power down to zero at the edges of the desired frequency range (lower-right panel of Fig. 8). Performing an inverse FFT on this filtered Fourier power spectrum results in the reconstruction of an H$\alpha$ lightcurve containing dominant periodicities of $\approx 210$ s, which can be seen in Fig. 13. This process is identical to convolving the detrended intensity time series with the given Gaussian frequency filter, but we perform this process step-by-step here for the purposes of clarity.

Fig. 13

The original HARDcam time series (upper solid black line), normalized by the quiescent H$\alpha$ continuum intensity, and displayed as a function of time. The lower solid black line is a Fourier filtered lightcurve, which has been detrended using a third-order polynomial (right panel of Fig. 6), convolved with a Gaussian frequency filter centered on 4.7 mHz with a FWHM of 2.0 mHz (lower-right panel of Fig. 8), before applying an inverse FFT to reconstruct the filtered time series. For visual clarity, the filtered lightcurve has been offset to bring it closer to the original time series intensities

It must be noted that here we employ a Gaussian frequency filter to smoothly transition the Fourier power to values of zero outside of the desired frequency range. However, other filter shapes can also be chosen, including Lorentzian, Voigt, or even custom profile shapes depending upon the level of smoothing required by the investigators. At present, there is no firm consensus regarding which filter profile shape is best to use, so it may be necessary to choose the frequency filter based upon the specifics of the data being investigated, e.g., the frequency resolution, the amplitude of the spectral components wishing to be studied, the width of the documented Fourier peaks, etc. Of course, we must remind the reader that isolating a relatively limited range of frequencies in Fourier space and transforming these back into real (temporal) space will always result in the appearance of a periodic signal at the frequency of interest, even if the derived Fourier transform was originally noise dominated. Therefore, it is necessary to combine confidence interval analysis (see Sect. 2.2.2) with such Fourier filtering techniques to ensure that only statistically significant power is being considered in subsequent analyses.

2.2.5 Fourier phase lag analysis

Many observational datasets will be comprised of a combination of multi-wavelength and/or multi-component spectral measurements. For example, the Rapid Oscillations in the Solar Atmosphere (ROSA; Jess et al. 2010b) instrument at the DST is able to observe simultaneously in six separate bandpasses. It is common practice to acquire contemporaneous imaging observations through a combination of G-band, 3500 Å and 4170 Å broadband continuum filters, in addition to Ca ii K, Na i D$_{1}$, and H$\alpha$ narrowband filters, which allows wave signatures to be studied from the depths of the photosphere through to the base of the transition region (e.g., Morton et al. 2011, 2012; Jess et al. 2012a, b, c; Kuridze et al. 2012; Grant et al. 2015; Krishna Prasad et al. 2015, 2016, 2017; Keys et al. 2018). On the other hand, Fabry–Pérot spectral imaging systems such as the Crisp Imaging Spectropolarimeter (CRISP; Scharmer et al. 2008) and the Interferometric Bi-dimensional Spectrometer (IBIS; Cavallini 2006), are able to capture two-dimensional spatial information (often including spectropolarimetric Stokes I/Q/U/V measurements) across a single or multiple spectral lines. This allows a temporal comparison to be made between various spectral parameters of the same absorption line, such as the full-width at half-maximum (FWHM), intensity, Doppler velocity, and magnitudes of circular/linear polarization (providing spectropolarimetric measurements are made). As a result, harnessing multi-wavelength and/or multi-component observations provides the ability to further probe the coupling of wave activity in the lower solar atmosphere.

The upper panel of Fig. 14 displays two synthetic intensity time series generated with a cadence of 1.78 s (consistent with the HARDcam H$\alpha$ data products overviewed in Sect. 2.1.1), each with a frequency of 5.6 mHz ($\approx 180$ s periodicity) and a mean intensity equal to 2. However, the red lightcurve (LC2) is delayed by $45^{\circ }$, and hence lags behind the black lightcurve (LC1) by 0.785 radians. As part of the standard procedures prior to the implementation of Fourier analysis (see, e.g., Sect. 2.2), each of the time series are detrended (in this case by subtracting a linear line of best fit) and apodized using a 90% tapered cosine apodization filter. The final intensity time series are shown in the lower panel of Fig. 14, and are now suitable for subsequent Fourier analyses.

Fig. 14

Synthetic time series (upper panel), each with a cadence of 1.78 s, displaying a frequency of 5.6 mHz ($\approx 180$ s periodicity) and a mean intensity equal to 2. The red lightcurve is delayed by 45 $^{\circ }$ (0.785 radians) with respect to the black lightcurve. The lower panel displays the detrended and apodized time series, which are now suitable for subsequent FFT analyses

Following the approaches documented in Sect. 2.2.2, FFTs of the detrended and apodized time series are taken, with 95% confidence levels calculated. The resulting FFT power spectral densities are shown in Fig. 15, where the red crosses indicate frequencies where the associated power is in excess of the calculated 95% confidence levels for each respective time series. It can be seen in both the upper and lower panels of Fig. 15 that the input 5.6 mHz signal is above the 95% confidence threshold for both LC1 and LC2. Next, the cross-power spectrum, $\varGamma _{12}(\nu )$, between the FFTs of LC1 and LC2 is calculated following the methods described by Bendat and Piersol (2000) as;

$$\begin{aligned} \varGamma _{12}(\nu ) = F(LC1) * \overline{F(LC2)} , \end{aligned}$$

(6)

with F denoting an FFT and ${\overline{F}}$ the complex conjugate of the FFT. The cross-power spectrum is a complex array (just like the FFTs from which it is computed), and therefore has components representative of its co-spectrum ($d(\nu )$; real part of the cross-power spectrum) and quadrature spectrum ($c(\nu )$; imaginary part of the cross-power spectrum). The co-spectrum from the input time series LC1 and LC2 is shown in the upper panel of Fig. 16. The red cross signifies the frequency where the Fourier power exceeded the 95% confidence level in both FFTs, namely 5.6 mHz, which is consistent with the synthetic lightcurves shown in Fig. 14.

Fig. 15

FFT power spectral densities for LC1 (upper panel) and LC2 (lower panel), corresponding to the solid black and red lines in the lower panel of Fig. 14, respectively. The red crosses highlight frequencies where the calculated Fourier power is above the 95% confidence level. It can be seen that the synthetic 5.6 mHz input signal is accurately identified in both corresponding power spectra, with its associated Fourier power being in excess of the 95% confidence threshold. The oscillatory behavior at high frequencies is due to the selected apodization filter

Fig. 16

Co-spectrum (upper panel; real part of the cross-power spectrum) of the input time series LC1 and LC2 shown in the lower panel of Fig. 14. The lower panel displays the phase angle between the input time series LC1 and LC2, which corresponds to the phase of the complex cross-spectrum. Here, a positive phase angle indicates that LC1 leads LC2 (i.e., LC2 lags behind LC1), which can be seen visually through examination of the individual lightcurves depicted in Fig. 14. The red crosses indicate the frequency where the calculated Fourier power for LC1 and LC2 both exceed the 95% confidence levels (see Fig. 15). The horizontal dashed blue line in the lower panel highlights a phase angle of 0 $^{\circ }$

Finally, the co-spectrum and quadrature spectrum can be utilized to calculate the phase lag between the input lightcurves LC1 and LC2 as a function of frequency, defined by Penn et al. (2011) as,

$$\begin{aligned} \phi (\nu ) = atan\left( \frac{\langle c(\nu ) \rangle }{\langle d(\nu ) \rangle }\right) . \end{aligned}$$

(7)

Here, the phase angle, commonly chosen to span the interval $-180^{\circ } \rightarrow +180^{\circ }$, is simply the phase of the complex cross-spectrum (see the nomenclature of Vaughan and Nowak 1997). The lower panel of Fig. 16 displays the calculated phase angles, again with the red cross highlighting the phase value at the frequency where the Fourier power exceeds the 95% confidence level in both FFTs corresponding to LC1 and LC2. In this example, the phase angle corresponding to a frequency of $\approx 5.6$ mHz is equal to $45^{\circ }$, which is consistent with the input lightcurves depicted in Fig. 14. Here, a positive phase angle indicates that LC1 leads LC2 (i.e., LC2 lags behind LC1), which can be visually confirmed in Fig. 14 with LC1 (solid black line) leading LC2 (solid red line).

It must be noted that phase angles can be computed for all possible frequencies (see, e.g., the lower panel of Fig. 16). However, it is important to determine which of these phase values are reliable before they are used in subsequent scientific interpretations. For the purposes of the example shown here, we selected that frequency at which both times series LC1 and LC2 demonstrated Fourier power exceeding the 95% confidence levels in both of their corresponding FFTs. However, a common alternative is to calculate the coherence level for each constituent frequency, which can then be employed (independently of the confidence levels) to pinpoint reliable frequencies in the corresponding cross-power spectrum. The coherence level is estimated from the normalized square of the amplitude of the complex cross-spectrum (see, e.g., Storch and Zwiers 1999), providing a measure, ranging between ‘0’ and ‘1’, of the linear correlation between the two input time series. Under this regime, values of ‘0’ and ‘1’ indicate no and perfect correlation, respectively. For the purposes of solar physics research, it is common to adopt a coherence value $>0.80$ to signify robust and reliable phase measurements (McAteer et al. 2003; Bloomfield et al. 2004a, b; Stangalini et al. 2013b, 2018; Kontogiannis et al. 2016).

Therefore, the cross-power spectrum and coherence are both used to examine the relationship between two time series as a function of frequency. The cross spectrum identifies common large power (i.e., significant peaks) at the same frequencies in the power spectra of the two time series, and whether such frequencies are related to each other (the relationship is quantified by phase differences). Such correlations cannot, however, be revealed if one or both time series do not have significant power enhancements at particular frequencies, e.g., if the power spectra at those frequencies are indistinguishable from red noise. Nonetheless, there still may be coherent modes at such frequencies, that can be identified in the coherence spectrum, i.e., two time series can have a large coherence at a frequency even though both or one of the power spectra do not show large power at that frequency. Thus, the coherence is a measure of the degree of linear correlation between the two time series at each frequency. In solar physics, the coherence is particularly useful when the two signals are associated to, e.g., different solar atmospheric heights (with, e.g., different amplitudes) and/or two different physical parameters. An example, from real observations, where oscillatory power (at specific time-frequency locations) appears only in one of the signals is demonstrated in Fig. 25. Hence, no significant power is detected in the cross-power spectrum, whereas a large coherence level, exceeding 0.8, is identified. The significance of phase measurements for reliable coherence values can be evaluated by either introducing a coherence floor level (e.g., the 0.8 threshold mentioned above) or estimating confidence levels. To approximate a floor level, Bloomfield et al. (2004a) randomized both time series for a very large number of realizations and calculated the coherence for each, from which, the threshold was estimated as an average over all realizations plus some multiples of the standard deviation of the coherence values. For the confidence levels, the coherence values should be tested against the null hypothesis of zero population coherence, i.e., whether the coherence exceeds expected values from arbitrary colored (e.g., white or red) noise backgrounds. While various methods have been employed for this statistical test, one common approach is to estimate the confidence levels by means of Monte Carlo simulations (Torrence and Compo 1998; Grinsted et al. 2004; Björg Ólafsdóttir et al. 2016).

With reliable phase angles calculated, it then becomes possible to estimate a variety of key wave characteristics. If T is the period of the wave, then the phase lag, $\phi$ (in degrees), can be converted into a physical time delay through the relationship,

$$\begin{aligned} \text {time delay (s)}\ = \frac{\phi }{360} \times T . \end{aligned}$$

(8)

The time delay value (arising from the measured phase lag) corresponds to a wave propagating between the different atmospheric layers. Of course, phase angles deduced from the co-spectrum and quadrature spectrum (see Eq. (7)) inherently have phase wrapping at $\pm 180^{\circ }$, hence introducing a $360^{\circ }$ ambiguity associated with the precise phase angle (as discussed in Centeno et al. 2006b; Beck et al. 2008; Cally 2009). Hence, the true time delay may need to include multiples of the period to account for the $360^{\circ }$ ambiguity, hence transforming Eq. (8) into,

$$\begin{aligned} \text {time delay (s)}\ = \frac{\phi }{360} \times nT , \end{aligned}$$

(9)

where n is a non-zero integer. Many studies to date have examined the propagation of relatively long-period oscillations (e.g., 100–300 s), which permit the assumption of $n=1$ without violating theoretical considerations (e.g., sound speed restrictions; Jess et al. 2012c), hence allowing direct use of Eq. (8). However, as future studies examine higher-frequency (lower-period) wave propagation, then more careful consideration of Fourier phase wrapping will need to be taken into consideration to ensure the derived time delay is consistent with the observations. As part of a phase ‘unwrapping’ process, the identification of quasi-periodic waves and/or those with modulated amplitudes will allow phase ambiguities to be practically alleviated. For example, by tracking the commencement of a wave, and hence the time delay as it propagates between closely-spaced atmospheric layers, the phase angle can be computed without the $\pm 360^{\circ }$ phase wrapping uncertainty. Alternatively, a modulated waveform will provide secondary peaks associated with the propagating group, which supplies additional information to better establish the precise value of n in Eq. (9), hence assisting with the phase unwrapping of the data, which will enable much more precise tracking of wave energy flux through the solar atmosphere.

Finally, if the geometric height separation, d (in km), between the two layers is known or can be estimated (González Manrique et al. 2020), then the average phase velocity, $v_{\text {ph}}$, of the wave propagating between these two distinct layers can be deduced via,

$$\begin{aligned} v_{\text {ph}}~\text {(km/s)}\ = \frac{360d}{T\phi } . \end{aligned}$$

(10)

Similar estimations of the phase velocities of embedded waves have been made by Mein (1977), Athay and White (1979), White and Athay (1979), Centeno et al. (2006b), Bello González et al. (2010a), Jess et al. (2012c), Grant et al. (2015) and Jafarzadeh et al. (2017c), to name but a few examples. Importantly, Eq. (10) can also be rearranged to estimate the atmospheric height separation between two sets of observations. For example, the acoustic sound speed is approximately constant in the lower photosphere, hence this value, alongside the derived time lag, can be utilized to provide an estimate of the height separation, d (e.g., Deubner and Fleck 1989).

2.3 Three-dimensional Fourier analysis

Telescope facilities deployed in a space-borne environment, which benefit from a lack of day/night cycles and atmospheric aberrations, have long been able to harness three-dimensional Fourier analyses to examine the temporal ($t \leftrightarrow \omega$) and spatial ($[x,y] \leftrightarrow [k_{x},k_{y}]$) domains. Here, t and $\omega$ represent the coupled time and frequency domains, respectively, while the [x, y] and $[k_{x},k_{y}]$ terms represent the coupled spatial distances and spatial wavenumbers in orthogonal spatial directions, respectively. Such three-dimensional Fourier analyses has been closely coupled with the field of helioseismology, which is employed to study the composition and structure of the solar interior by examining large-scale wave patterns on the solar surface (Kosovichev et al. 1997; Turck-Chièze 2001; Braun and Lindsey 2000; Gizon et al. 2010; Kosovichev 2011; Buldgen et al. 2019), which often give rise to patterns consistent with ‘rings’ and ‘trumpets’ when viewed in Fourier space (Hill 1988).

Up until recently, it has been challenging to apply the same three-dimensional Fourier techniques to high-resolution datasets from ground- and space-based observatories (Leighton 1963; Spruit et al. 1990). These techniques have been applied with ground-based observations to study convective phenomena (Chou et al. 1991; Straus et al. 1992) and plage (Title et al. 1992). With the advent of high image pointing stability, brought to fruition through a combination of high-order AO, photometrically accurate image reconstruction algorithms, precise telescope control hardware, and sub-pixel cross-correlation image co-alignment software, it is now possible to achieve long-duration image and/or spectral sequences that are stable enough to allow Fourier analyses in both temporal and spatial domains. The benefit of using high-resolution facilities is that they offer unprecedented Nyquist temporal frequencies ($\omega$) and spatial wavenumbers ($[k_{x},k_{y}]$) due to their high temporal and spatial sampling, respectively. For example, the HARDcam H$\alpha$ dataset described in Sect. 2.1.1 has a temporal cadence of 1.78 s and a spatial sampling of $0{\,}.{\!\!}{''}138$ per pixel, providing a Nyquist frequency of $\omega _{{\text {Ny}}}\approx 280$ mHz $\left( \frac{1}{2\times 1.78}\right)$ and a Nyquist wavenumber of $k_{{\text {Ny}}}\approx 22.8$ arcsec$^{-1}$ $\left( \frac{2\pi }{2\times 0.138}\right)$. This allows for the examination of the smallest and most rapidly varying phenomena currently visible in such high-resolution datasets.

Applying an FFT to a three-dimensional dataset converts the spatial/temporal signals, [x, y, t], into its frequency counterparts, [$k_{x}, k_{y}, \omega$]. An example of this process can be seen in Fig. 17, whereby an FFT has been applied to the HARDcam H$\alpha$ dataset documented by Grant et al. (2018). It can be seen in the right panel of Fig. 17 that the Fourier power signatures are approximately symmetric in the $k_{x}/k_{y}$ plane. As a result, it is common for [$k_{x}, k_{y}$] cross-cuts at each frequency, $\omega$, to be azimuthally averaged providing a more straightforward two-dimensional representation of the Fourier power in the form of a $k{-}\omega$ diagram (Duvall et al. 1988; Krijger et al. 2001; Rutten and Krijger 2003; Kneer and Bello González 2011; Jess et al. 2012c, 2017).

Fig. 17

An example application of an FFT to a three-dimensional datacube, converting [x, y, t] (left) into its frequency counterparts [$k_{x}, k_{y}, \omega$] (right). The HARDcam H$\alpha$ dataset presented here is taken from the work of Grant et al. (2018)

An azimuthally averaged $k{-}\omega$ diagram for the HARDcam H$\alpha$ sunspot observations described in Sect. 2.1.1 is shown in the right panel of Fig. 18. A number of important features are present in this diagram, including consistency with many quiet-Sun and internetwork Fourier power peaks documented by Krijger et al. (2001), Kneer and Bello González (2011) and Jess et al. (2012c), whereby high power observed at larger spatial wavenumbers tends to be correlated with higher temporal frequencies. This can be visualized in the right panel of Fig. 18, whereby the dominant Fourier power is associated with the smallest spatial wavenumbers and temporal frequencies. However, as the wavenumber is increased to $>1$ arcsec$^{-1}$, the temporal frequencies corresponding to maximal Fourier power are concentrated within the $3{-}6$ mHz interval. This is consistent with the general trends observed in classical photospheric $k{-}\omega$ diagrams, such as that shown in Fig. 19. Here, two $k{-}\omega$ diagrams from the photospheric SDO/AIA 1700 Å time series that is co-spatial (and overlaps temporally) with the HARDcam H$\alpha$ observations (used to produce Fig. 18) are displayed. The information displayed in both panels of Fig. 19 is identical, however, the left panel is displayed on a linear wavenumber (k) and frequency ($\omega$) scales, while the right panel is displayed on log–log axes. In both panels, similar trends (e.g., heightened Fourier power with increasing temporal frequency in the interval of 3–6 mHz is linked to larger spatial wavenumbers) can be identified, which is consistent with the overall trends depicted in the right panel of Fig. 18. However, as discussed in Jess et al. (2017), within the region highlighted by the solid black box in the right panel of Fig. 18, there is evidence of elevated Fourier power that spans a large range of spatial scales, yet remains confined to a temporal frequency on the order of 5.9 mHz ($\approx 170$ s). This suggests that the embedded wave motion has strong coherency across a broad spectrum of spatial scales, yet can be represented by a very narrow range of temporal frequencies. Looking closely at the right panel of Fig. 18, it can be seen that elevated levels of Fourier power extend down to the smallest spatial wavenumbers allowable from the HARDcam dataset. This implies that the 5.9 mHz frequency is still significant on spatial scales much larger than the field of view captured by HARDcam.

Fig. 18

Image reproduced with permission from Jess et al. (2017), copyright by the authors

A two-dimensional $[k_{x},k_{y}]$ cross-cut for a single temporal frequency, $\omega$, corresponding to the HARDcam H$\alpha$ data acquired on 2011 December 10 and described in Sect. 2.1.1 (left panel). Due to the symmetries often found between $k_{x}$ and $k_{y}$, it is common to perform azimuthal averaging (e.g., along the solid green contour) to collapse this two-dimensional information into a single dimension, i.e. $[k_{x}, k_{y}] \rightarrow [k]$. This allows the three-dimensional FFT cube (see, e.g., the right panel of Fig. 17) to be simplified into a standardized two-dimensional image, forming a $k{-}\omega$ diagram (right panel). Here, the $k{-}\omega$ diagram is cropped between approximately $1<\omega <10$ mHz and $0.3<k<10.0$ arcsec$^{-1}$, and displayed on a log–log scale to assist visual clarity. The colors represent oscillatory power that is displayed on a log-scale, while the vertical dashed and dotted lines correspond to the spatial size of the umbral diameter ($\approx 13{\,}.{\!\!}{''}50$) and the radius of the umbra ($\approx 6{\,}.{\!\!}{''}75$), respectively. The solid black box indicates a region of excess wave power at $\approx 5.9$ mHz ($\approx 170$ s) over the entire spatial extent of the sunspot umbra.

Fig. 19

A set of $k{-}\omega$ diagrams, derived from the photospheric SDO/AIA 1700 Å time series of active region NOAA 11366, which is co-spatial (and overlaps temporally) with the chromospheric HARDcam measurements presented in Fig. 18. Both $k{-}\omega$ diagrams are identical, however, the left panel is displayed on linear wavenumber (k) and frequency ($\omega$) scales, while the right panel is displayed on log–log axes. It is clear from inspection of the two panels that each have their merit when presenting results, with the linear axes giving less visual emphasis to the lower wavenumbers/frequencies, while the log–log axes allowing power-law trends in the power spectral densities to be modeled more easily through straight-line fitting

However, there are a number of key points related to Figs. 18 and 19 that are worth discussing. First, Fig. 19 highlights the merits of utilizing either linear or log–log axes depending on the features being examined. For example, the use of a linear scale (left panel of Fig. 19) results in less visual emphasis being placed on the lowest spatial waveneumbers and temporal frequencies. This can help prevent (visual) over-estimations of the trends present in the $k{-}\omega$ diagram since all of the frequency bins occupy identical sizes within the corresponding figure. However, as spatial and temporal resolutions dramatically improve with next generation instrumentation, the corresponding spatial/temporal Nyquist frequencies continue to become elevated, often spanning multiple orders-of-magnitude. If these heightened Nyquist frequencies are plotted on a purely linear scale, then many of the features of interest may become visually lost within the vast interval occupied by the $k{-}\omega$ diagram. An option available to counter this would be to crop the $k{-}\omega$ diagram to simply display the spatial wavenumbers and temporal frequencies of interest, although this comes at the price of discarding information that may be important within the remainder of the frequency space. Alternatively, it is possible to use log–log axes for the $k{-}\omega$ diagram, which can be visualized in the right panels of Figs. 18 and 19. This type of log–log display also benefits the fitting of any power-law trends that may be present within the $k{-}\omega$ diagram, since they will manifest as more straightforward (to fit) linear slopes in the plot. Finally, the right panel of Fig. 18 reveals some horizontal banding of power that appears slightly different than the diagonal ‘arms’ of Fourier power visible in Fig. 19. This may be a consequence of the reduced spatial wavenumber and temporal frequency resolutions achievable with large-aperture ground-based observatories, which naturally have a reduced field-of-view size (causing a relatively low spatial wavenumber resolution when compared to large field-of-view observations from, e.g., SDO) and limited time series durations (creating relatively low temporal frequency resolutions when compared to space-borne satellite missions that are unaffected by day/night cycles and/or atmospheric effects). Therefore, it is imperative that the investigative team examines the merits of each type of $k{-}\omega$ display and selects the use of either linear or log–log axes to best represent the physical processes at work in their dataset.

2.3.1 Three-dimensional Fourier filtering

Taking the one-dimensional Fourier filtering methodology described in Sect. 2.2.4 a step further, it is often useful to filter an input three-dimensional dataset ([x, y, t]) in terms of both its temporal frequencies, $\omega$, and its spatial wavenumbers, k. While it is common for the frequency to be defined as the reciprocal of the period, i.e., $\omega = 1/T$, where T is the period of oscillation, the wavenumber is often defined as $k = 2\pi /\lambda$ (Krijger et al. 2001), where $\lambda$ is the wavelength of the oscillation in the spatial domain (i.e., [x, y]). Hence, it is often important to bear in mind this additional factor of $2\pi$ when translating between wavenumber, k, and spatial wavelength, $\lambda$. Figures 18 and 20 employ this form of frequency/wavenumber notation, meaning that the spatial wavelengths can be computed as $\lambda = 2\pi /k$, while the period is simply $T = 1/\omega$ (similar to that shown in Straus et al. 1992; Jess et al. 2012c). However, some research programs, particularly those adopting helioseismology nomenclature, utilize the factor of $2\pi$ in both the wavenumber and frequency domains (e.g., $T = 2\pi /\omega$; Mihalas and Toomre 1981). As a result, it is important to select an appropriate scaling to ensure consistency across a piece of work. An example code capable of doing three-dimensional Fourier filtering is the QUEEn’s university Fourier Filtering (QUEEFF; Jess et al. 2017) algorithm, which is based around the original techniques put forward by Tarbell et al. (1988), Title et al. (1989), Rutten and Krijger (2003), Roth et al. (2010) and Krijger et al. (2001), but now adapted into a publicly available Interactive Data Language (idl; Stern 2000) package.²$^{,}$³

Fig. 20

Outputs provided by a commonly available three-dimensional Fourier filtering code (QUEEFF; Jess et al. 2017), showing a frequency-averaged wavenumber spectrum (upper-left), a Gaussian (with $2<k<10$ arcsec$^{-1}$) wavenumber filter that resembles a torus shape when viewed in the $[k_{x}, k_{y}]$ plane (upper-middle), and the resulting transmitted wavenumber spectra once multiplied by the chosen filter (upper-right). The lower panel displays the wavenumber-averaged frequency spectrum (solid black line), where the Fourier power is displayed (using a log-scale) as a function of the temporal frequency, $\omega$. The dashed blue line highlights a chosen frequency filter, $20\pm 10$ mHz, with a Gaussian shape to more smoothly reduce Fourier power at the edges of the chosen spectral range to reduce aliasing. The solid red line displays the resulting transmitted frequency spectrum once multiplied by the chosen Gaussian filter. In each panel, dashed black or white lines highlight the $k_{x}/k_{y}=0$ arcsec$^{-1}$ or $\omega =0$ mHz locations

Importantly, the QUEEFF code provides the user with the ability to apply Gaussian smoothing windows to both frequency and wavenumber regions of interest in order to help mitigate against elements of aliasing during subsequent dataset reconstruction. Figure 20 shows an example figure provided by the QUEEFF code, which displays the frequency-averaged wavenumber power (upper-left panel), the chosen wavenumber filter (upper-middle panel) utilizing a Gaussian structure providing a torus-shaped filter spanning 2–10 arcsec$^{-1}$, alongside the resulting filtered wavenumber spectra (upper-right panel). The lower panel of Fig. 20 displays the spatially-averaged frequency spectrum of the HARDcam H$\alpha$ dataset, where the Fourier power is displayed as a function of the frequency, $\omega$, using a solid black line. A Gaussian frequency filter, spanning $20\pm 10$ mHz, is overplotted using a dashed blue line. The preserved temporal frequencies (i.e., once the original frequency spectrum has been multiplied by the chosen frequency filter) is shown using a solid red line. This filtered three-dimensional Fourier cube can then be passed through an inverse FFT to reconstruct an intensity image cube that contains the wavenumbers and frequencies of interest to the user.

Again, as discussed in Sect. 2.2.4, the QUEEFF three-dimensional Fourier filtering code constructs a Gaussian-shaped filter, which is applied in the Fourier domain. This ensures that the filter is symmetric about the chosen peak frequency (see, e.g., the black line in the left panel of Fig. 21). Of course, due to the oscillation period having a reciprocal relationship with the temporal frequency (i.e., $1/\omega$), this results in asymmetric sampling about the desired peak period (see, e.g., the solid black line in the right panel of Fig. 21). Depending upon the science requirements of the user, it may be more advantageous to apply a Gaussian-shaped filter in the period domain (e.g., the solid blue line in the right panel of Fig. 21), which ensures less inclusion of lower frequency (higher period) terms that may be undesirable in the final reconstructed time series. This is highlighted by the more rapid truncation of the filter (solid blue line in the left panel of Fig. 21) towards lower frequencies. Additionally, the user may select alternative frequency filters, such as a Voigt profile (Zaghloul 2007), which is shown in Fig. 21 using a solid red line. Furthermore, Fig. 21 shows possible filtering combinations that can be applied to the temporal domain, yet similar options are available when filtering the spatial wavenumbers ($[k_{x}, k_{y}]$) too. Ultimately, it is the science objectives that drive forward the wave filtering protocols, so possible options need to be carefully considered before applying to the input data.

Fig. 21

Different types of frequency ($\omega$) filter that can be applied to time-resolved data products. The left panel displays the filter transmission (as a percentage) in terms of the frequency, while the right panel displays the same filters as a function of the oscillatory period. Presented using a solid black line is a Gaussian-shaped filter in the frequency domain with a FWHM equal to 10 mHz, while the solid red line indicates a Voigt-shaped filter in the frequency domain, both centered on 20 mHz. Contrarily, a Gaussian-shaped filter in the period domain, with a FWHM equal to 10 s, is shown using a solid blue line, again centered on 50 s to remain consistent with the 20 mHz profiles shown using red and black lines. It is clearly evident that the filter profile shape changes dramatically between the time and frequency domains, and hence it is important to select the correct filter based upon the science requirements

Combination Fourier filters (i.e., that are functions of $k_{x}$, $k_{y}$ and $\omega$) have been utilized in previous studies to extract unique types of wave modes manifesting in the lower solar atmosphere. For example, specific Fourier filters may be employed to extract signatures of f- and p-mode oscillations manifesting in photospheric observations (e.g., Hill 1988; Schou et al. 1998; Gizon and Birch 2004; Bahauddin and Rast 2021). Another example of a well-used Fourier filter is the ‘sub-sonic filter’, which can be visualized as a cone in $k{-}\omega$ space (Title et al. 1989),

$$\begin{aligned} v_{\text{ph}} = \frac{\omega }{k} , \end{aligned}$$

(11)

where $v_{\text{ph}}$ is the phase velocity of the wave. Here, all Fourier components inside the cone, where propagation velocities are less than the typical sound speed (i.e., $v_{\text{ph}} < c_{s}$), are retained while velocities outside the cone are set to zero. An inverse Fourier transform of this filtered spectrum provides a dataset that is embodied by the convective part of the solar signal since the non-convective phenomena (e.g., solar p-modes) have been removed (Straus and Bonaccini 1997; Rutten and Krijger 2003). Alternatively, modification of the sub-sonic filter to include only those frequencies above the Lamb mode, $\omega = c_{s}k$ (Fleck et al. 2021), provides a reconstructed dataset containing oscillatory parts of the input signal. As highlighted above, it is the science objectives that define the filtering sequences required to extract the underlying time series of interest. However, well-proven examples of these exist for common phenomena (e.g., solar f- and p-modes), hence providing an excellent starting point for the community.

2.4 Wavelet analyses

While FFT analyses is very useful for identifying and characterizing persistent wave motion present in observational datasets, it begins to encounter difficulties when the time series contains weak signals and/or quasi-periodic signatures. Figure 22 shows example time series containing a persistent wave signal with a 180 s periodicity (5.56 mHz) and no embedded noise (top-left panel), a quasi-periodic 5.56 mHz wave signal with no noise (middle-left panel), and a qausi-periodic 5.56 mHz wave signal embedded in high-amplitude noise (lower-left panel). It can be seen for each of the corresponding right-hand panels, which reveal the respective Fourier power spectral densities, that the detected 5.56 mHz Fourier peak becomes progressively less apparent and swamped by noise, even becoming significantly broadened in the lower-right panel of Fig. 22. As a result, the application of Fourier analyses to solar time series often displaying quasi-periodic wave motion (e.g., spicules, fibrils, rapid blueshift excursions (RBEs), etc.; Beckers 1968; De Pontieu et al. 2004, 2007a, b; Zaqarashvili and Erdélyi 2009; Sekse et al. 2013a, b; Kuridze et al. 2015) may not be the most appropriate as a result of the limited lifetimes associated with these features.

Fig. 22

An example time series consisting of a pure 180 s periodicity (5.56 mHz) signal, which is sampled at a cadence of 1.44 s to remain consistent with modern instrument capabilities (upper left). The middle-left panel shows the same example time series, only now with the first three and last two complete wave cycles suppressed, hence making a quasi-periodic wave signal. The lower-left panel shows the same quasi-period wave signal shown in the middle-left panel (solid green line), only now with superimposed Poisson (shot) noise added on top of the signal. Each of the right panels display the corresponding FFT-generated Fourier spectra, with the frequency and Fourier power values plotted on log-scales for better visual clarity. The vertical dashed red lines highlight the input 5.56 mHz signal

Wavelet techniques, pioneered by Torrence and Compo (1998), employ a time-localized oscillatory function that is continuous in both time and frequency (Bloomfield et al. 2004b), which allows them to be applied in the search for dynamic transient oscillations. The time resolution of the input dataset is preserved through the modulation of a simple sinusoid (synonymous with standard FFT approaches) with a Gaussian envelope, providing the Morlet wavelet commonly used in studies of waves in the solar atmosphere (Bloomfield et al. 2004a; Jess et al. 2007; Stangalini et al. 2012; Kobanov et al. 2013, 2015; Jafarzadeh et al. 2017d). As a result, a wavelet transform is able to provide high frequency resolution at low frequencies and high time resolution at high frequencies, which is summarized by Kehtarnavaz (2008).

Figure 23 displays the wavelet power spectrum (lower panel) resulting from the application of a Morlet wavelet transform on the detrended and apodized HARDcam H$\alpha$ lightcurve (upper panel). Here, it is possible to see the effects of quasi-periodic wave phenomena, where there is clear evidence of a large-amplitude periodicity between times of $0{-}2200$ s at a period of $\approx 210$ s ($\approx 4.7$ mHz). This wave activity is highlighted in the wavelet transform by being bounded by the 95% confidence level isocontours across these times and periods, which is equivalent to the oscillatory behavior being significant at the 5% level (Torrence and Compo 1998). To calculate the wavelet power thresholds corresponding to the 95% confidence isocontours, the wavelet background spectrum (i.e., the output theoretical background spectrum that has been smoothed by the wavelet function) is multiplied by the 95$^{\text {th}}$ percentile value for a $\chi _{2}^{2}$ distribution (Gilman et al. 1963). Please note that for some considerations, including expensive computation times, the Monte Carlo randomization method is not preferred for wavelet transform (Lau and Weng 1995; Torrence and Compo 1998). The $\approx 210$ s wavelet power signatures shown in the lower panel of Fig. 23 are consistent with the standardized FFT approach documented in Sect. 2.2, although the quasi-periodic nature of the wave motion is likely a reason why the corresponding power in the traditional FFT spectrum (upper panel of Fig. 8) is not as apparent. Importantly, with the wavelet transform it is possible to identify more clearly the times when this periodicity appears and disappears from the time series, which is seen to correlate visibly with the clear sinusoidal fluctuations present at the start of the H$\alpha$ time series (upper panel of Fig. 23). Also, the lack of significant wavelet power at very long periods (low frequencies) suggests that the lightcurve detrending applied is working adequately.

Fig. 23

The detrended and apodized HARDcam H$\alpha$ lightcurve shown in the lower panel of Fig. 6 (top). The bottom panel shows the corresponding wavelet transform, where the wave power is displayed as a function of the oscillatory period (y-axis) and observing time (x-axis). The color bar displays the normalized wavelet power, while the cross-hatched region (bounded by a dashed white line) highlights locations of the wavelet transform that may be untrustworthy due to edge effects. Solid white lines contour regions where the wavelet power exceeds the 95% confidence level (i.e., significant at the 5% level)

Due to wavelet analyses preserving the time domain of the original input signal, care must be taken to ensure that any power visible in wavelet transforms is the result of wave motion and not an instantaneous spike in intensity. To achieve this, it is typical to exclude oscillations from subsequent analysis that last, in duration, for less than $\sqrt{2}$ wave cycles. This requirement is often referred to as the decorrelation time (Torrence and Compo 1998), which involves comparing the width of a peak in the wavelet power spectrum (defined as the time interval over which the wavelet power exceeds the 95% confidence level—see Sect. 2.2.2) with the period itself to determine the number of complete wave cycles (Ireland et al. 1999; McAteer et al. 2004). Oscillations that last for less time than $\sqrt{2}$ wave cycles are subsequently discarded as they may be the result of spikes and/or instrumental abnormalities in the data. In addition, periodicities manifesting towards the extreme edges of the lightcurve need to be considered carefully due to the possible presence of edge effects arising due to the finite duration of the time series (Meyers et al. 1993). This region where caution is required is highlighted in the lower panel of Fig. 23 using the cross-hatched solid white lines. Here, the “cone of influence” (COI) is defined as as the e-folding time for the autocorrelation of wavelet power at each scale, and for the traditional Morlet wavelet this is equal to $\sqrt{2}$ wave cycles (Torrence and Compo 1998), hence why longer periods are more heavily effected (in the time domain) than their shorter (high-frequency) counterparts.

Finally, many research studies employ the global wavelet spectrum to characterize the frequencies present in the input time series. Here, the global wavelet spectrum is defined as the average spectrum across all local wavelet spectra along the entire input time axis (Torrence and Compo 1998). Essentially, the global wavelet spectrum can be considered as an estimation of the true Fourier spectrum. For example, a time series comprised of mixed wave frequencies that are superimposed on top of a white noise background should produce Fourier spectral peaks equal to $2\sigma _{\epsilon }^{2} + NA_{i}^{2}/2$, where $A_{i}$ are the amplitudes of the oscillatory components, $\sigma _{\epsilon }^{2}$ is the variance of the noise, and N is the number of steps in the time series (Priestley 1981). However, the corresponding peaks in the global wavelet spectrum will usually be higher at larger scales when compared to smaller scales, which is a consequence of the wavelet transform having better frequency resolution at long periods, albeit with worse time localization (Marković and Koch 2005).

As such, the global wavelet spectrum is often considered a biased estimation of the true Fourier spectrum (Wu and Liu 2005). This effect can be clearly seen in Fig. 24, which displays both the Fourier and global wavelet power spectra for the same HARDcam H$\alpha$ time series shown in the lower panel of Fig. 6. In Fig. 24, the higher power at larger scales (lower frequencies) is visible in the global wavelet spectrum (red line), when compared to that derived through traditional Fourier techniques (black line). However, at smaller scales (higher frequencies), both the global wavelet and Fourier spectra are in close agreement with one another, with the global wavelet spectrum appearing as a smoothed Fourier spectrum. The reason for these effects is due to the width of the wavelet filter in Fourier space. At large scales (low frequencies), the wavelet is narrower in frequency, resulting in sharper peaks that have inherently larger amplitudes. Contrarily, at small scales (high frequencies), the wavelet is more broad in frequency, hence causing any peaks in the spectrum to become smoothed (Torrence and Compo 1998). As such, it is important to take such biases into consideration when interpreting any embedded wave motion. Indeed, Banerjee et al. (2001), Christopoulou et al. (2003), Samanta et al. (2016), Kayshap et al. (2018) and Chae et al. (2019) have discussed the implementation of global wavelet and Fourier power spectra in the context of solar oscillations.

Fig. 24

Fourier (black line) and global wavelet (red line) power spectra of the HARDcam H$\alpha$ detrended lightcurve shown in the lower panel of Fig. 6. It can be seen that at larger scales (lower frequencies) the global wavelet spectrum has increased power over that calculated from traditional Fourier techniques, due to the increased wavelet frequency resolution in this regime. Contrarily, at smaller scales (higher frequencies) the global wavelet spectrum appears as a smoothed Fourier spectrum due to the reduced frequency resolution at these smaller scales. While the global wavelet spectrum is a good estimation of the Fourier power spectrum, these biases need to be carefully considered when interpreting the embedded wave motion

2.4.1 Wavelet phase measurements

Similar to the Fourier phase lag analysis described in Sect. 2.2.5, it is also useful to obtain phase angles, cross-power spectrum, and coherence between wavelet power spectra at different wavelengths, spatial locations, and/or multi-component spectral measurements. Hence, the phase angles are determined not only as a function of frequency, but also as a function of time. These phase angles are usually demonstrated as small arrows on a wavelet co-spectrum (or wavelet coherence) map, where their directions indicate the phase angles at different time-frequency locations. The convention with which an arrow direction represents, e.g., zero and $90^{\circ }$ phase angles (and which lightcurve leads or lags behind) should be specified.

Reproduced from Jafarzadeh et al. (2017d), the lower- and upper-left panels of Fig. 25 display two wavelet power spectra (from a Morlet wavelet transform) of transverse oscillations in a small magnetic element (marked with circles in Fig. 5) at two atmospheric heights sampled by the SuFI/Sunrise 300 nm and Ca ii H bands (with an average height difference of $\approx 450$ km), respectively. Islands of high power, particularly those marked by the 95% confidence level contours, are evident in both wavelet power spectra. The wavelet co-spectrum and coherence maps of these two power spectra are shown in the upper- and lower-right panels of Fig. 25, respectively. The phase-lag arrows are over plotted on the entire cross-power spectrum, while the same arrows are depicted on the latter map only where the coherence exceeds 0.7. Here, the arrows pointing right represent in-phase oscillations and those pointing straight up identify 90 degrees phase lags where the oscillations in 300 nm lag behind those observed in the Ca ii H time series. Note here the changes of phase lags from one time-frequency region to another, particularly, in regions with confidence levels larger than 95%, and/or areas with coherence exceeding 0.7 (or 0.8). However, most of the arrows point upwards (with different angles) in this example, implying an upward wave propagation in the lower solar atmosphere (i.e., from the low photosphere, sampled by the 300 nm band, to the heights corresponding to the temperature minimum/low chromosphere, sampled by the Ca ii H images). A slight downward propagation is also observed in a small area. These may associate to various wave modes and/or oppositely propagating waves at different frequencies and times. We note that such phase changes with time could not be identified using a Fourier phase lag analysis (see Sect. 2.2.5), where phase angles are computed as a function of frequency only.

Fig. 25

Images reproduced with permission from Jafarzadeh et al. (2017d), copyright by AAS

Wavelet power spectra of transverse oscillations in a small magnetic element (marked with circles in Fig. 5), from time-series of images acquired in 300 nm (lower left) and in Ca ii H (upper left) bands from SuFI/Sunrise. The right panels display the wavelet co-spectrum power (on the top) and coherence map (on the bottom). The 95% confidence levels are identified with dashed/solid contours in all panels and the COIs are marked with the cross-hatched/shaded regions. The arrows on the right panels show the phase angles between oscillations at the two atmospheric heights, with in-phase oscillations depicted by arrows pointing right and fluctuations in Ca ii H leading those in 300 nm by $90^{\circ }$ marked by arrows pointing straight up.

Whether the cross-power spectrum or coherence should be used for the wave identification greatly depends on the science and the types of data employed. While the co-spectrum (which is obtained through multiplying the wavelet power spectrum of a time series by the complex conjugate of the other) identifies regions with large power in common (between the two time series), the coherence (i.e., square of the cross-spectrum normalized by the individual power spectra; Grinsted et al. 2004) highlights areas where the two time series co-move, but not necessarily sharing a high common power. An example is the area around the time and period of 70 and 47 s, respectively, that is associated to a coherence level exceeding 0.8 (and within the 95% confidence levels), but with no significant power in the co-spectrum (only one of the power spectra, i.e., that from the Ca ii H data, show large power at that time and period location).

As a working example, from the right panels of Fig. 25, the phase lag at the time and period of 75 and 41 s, respectively, reads about 140 $^{\circ }$, which is translated to a time lag of $\approx 16$ s. Given the average height difference of 450 km between the two atmospheric layers, it results in a wave propagation speed of $\approx 28$ km/s (due to the transverse oscillations in the small-scale magnetic element). A similar analysis for intensity oscillations in the same small-scale magnetic element has also been presented in Jafarzadeh et al. (2017d). Of course, as highlighted in Sect. 2.2.5, phase measurements are always subject to an associated uncertainty of $\pm 360 ^{\circ }$ ($\pm 2\pi$), which arises via phase wrapping. As a consequence, to alleviate ambiguities in phase angles, in addition to subsequently derived phase velocities, care must be taken to select observational time series where the atmospheric height separation is not too substantial (see Sect. 2.2.5 for more discussion), which helps to minimize the ambiguities associated with phase wrapping.

Depending on science objectives, it may be helpful to inspect the variation of phase lags with frequency (or period). To this end, a statistical phase diagram can be created, where all reliable phase angles (e.g., those associated to power significant at 5%, and/or with a coherence exceeding 0.8) are plotted as a function of frequency (Jess et al. 2012c). Such a phase diagram can provide information about the overall wave propagation in, e.g., similar magnetic structures. Figure 26 illustrates a phase diagram (i.e., a 2D histogram of phase angle versus period; from Jafarzadeh et al. 2017d) constructed from all the reliable phase angles obtained from the transverse oscillations in 7 small magnetic elements, similar to that discussed above. The background colors represent the occurrence frequency and the contours mark regions which are statistically significant (i.e., compared to the extreme outliers). From this phase diagram, it is evident that the upward propagating waves (i.e., the positive phase angles in the convention introduced here) appear preferential.

Fig. 26

Image reproduced with permission from Jafarzadeh et al. (2017d), copyright by AAS

Phase diagram of transverse oscillations in 7 small magnetic elements observed in two layers of the lower solar atmosphere (with $\approx 450$ km height difference) from SuFI/Sunrise.

2.5 Empirical mode decomposition

Empirical Mode Decomposition (EMD; Huang et al. 1998, 1999) is a statistical tool developed to decompose an input time series into a set of intrinsic timescales. Importantly, EMD is a contrasting approach to traditional FFT/wavelet analyses since it relies on an empirical approach rather than strict theoretical tools to decompose the input data. Due to the decomposition being based on the local characteristic timescales of the data, it may be applied to non-linear and non-stationary processes without the detrending often applied before the application of Fourier-based techniques (i.e., under the assumption that such detrending is able to accurately characterize any non-stationary and/or non-periodic fluctuations in the time series with a low-order polynomial). As such, it is possible for EMD to overcome some of the limitations of FFT/wavelet analyses, including aspects of wave energy leakage across multiple harmonic frequencies (Terradas et al. 2004). non-stationary/non-period fluctuations that can be characterized by a low-order polynomial

Following the methodology described by Terradas et al. (2004), we apply EMD techniques to the HARDcam H$\alpha$ time series depicted in the upper-left panel of Fig. 6. To begin, extrema in the lightcurve are identified, and are then connected by a cubic spline fit to provide an upper envelope of the positive intensity fluctuations (i.e., fluctuations above the mean). Next, the same process is applied to find the lower envelope corresponding to negative intensity fluctuations (i.e., fluctuations below the mean). The mean value between the upper and lower envelopes, at each time step, is denoted $m_{1}(t)$. The difference between the original input data and the mean function is called the first component, $h_{1}(t)$. Providing the input time series contains no undershoots, overshoots, and/or riding waves (Huang et al. 1998), then the first intrinsic mode function (IMF) is equal to $h_{1}(t)$.

Unfortunately, many input time series contain signal blemishes, and removal of the first component, $h_{1}(t)$, from the original lightcurve will generate additional (false) extrema. Hence, to mitigate against these potential issues, the above procedure is repeated numerous times until the first true IMF is constructed (see Huang et al. 1998, for more information). The first IMF constructed, $c_{1}(t)$, is comprised of the most rapid fluctuations of the signal. This can then be subtracted from the original time series, producing a residual lightcurve made up of longer duration fluctuations. The process can subsequently be repeated numerous times to extract additional IMFs until the amplitude of the residual lightcurve falls below a predetermined value, or becomes a function from which no more IMFs can be extracted (Terradas et al. 2004).

Figure 27 shows a collection of IMFs extracted from the HARDcam H$\alpha$ time series depicted in the upper-left panel of Fig. 6. It is clear that the most rapid fluctuations are present in IMF $c_{1}$, with IMF $c_{8}$ documenting the slowest evolving intensity variations. Plotted on top of IMF $c_{8}$ is the original H$\alpha$ time series, along with the polynomial best-fit line (dashed red line) used to detrend the lightcurve in Sect. 2.2 before the application of FFT/wavelet techniques. The global trends highlighted by IMF $c_{8}$ and the polynomial best-fit line are similar, again highlighting the appropriate use of detrending in Sect. 2.2, but now compared with generalized empirical methods. Figure 28 displays the 8 extracted IMFs in the form of a two-dimensional map, which can often be used to more readily display the corresponding interplay between the various amplitudes and variability timescales.

Fig. 27

IMFs $c_{1} \rightarrow c_{8}$, extracted from the original (non-detrended and non-apodized) HARDcam H$\alpha$ time series overplotted in the lower-right panel. In addition, the lower-right panel also shows the polynomial best-fit line (dashed red line) used to detrend the data prior to FFT/wavelet analyses. It can be seen that the longest period fluctuations making up IMF $c_{8}$ are similar to the global trend line calculated in Sect. 2.2.Note that a summation of IMFs $c_{1} \rightarrow c_{8}$ will return the original signal

Fig. 28

IMFs $c_{1} \rightarrow c_{8}$ displayed as a two-dimensional map (left), where yellow and blue colors represent the peaks and troughs, respectively, of the IMF intensity fluctuations. The horizontal dashed black lines represent cuts through each IMF, with the corresponding intensity time series displayed in the right panel. The horizontal dashed red lines represent the zero value corresponding to each IMF

Once the IMFs have been extracted from the input time series, it is possible to employ Hilbert spectral analysis (Huang et al. 1998; Oppenheim and Schafer 2009) to examine the instantaneous frequencies with time for each IMF. The combined application of EMD and Hilbert spectral analysis is often referred to as the Hilbert–Huang transformation (Huang and Wu 2008). From the outputs of the Hilbert–Huang transformation, it is possible to display the instantaneous frequencies for each of the extracted IMFs as a function of time. The left panel of Fig. 29 displays the instantaneous frequencies corresponding to IMFs $c_{2} \rightarrow c_{7}$ using the purple, blue, dark green, light green, orange, and red lines, respectively. IMFs $c_{1}$ and $c_{8}$ have been removed from the plot as these correspond to very high and low frequency fluctuations, respectively, which clutter the figure if included. The solid colored lines represent the running mean values of the instantaneous frequencies (calculated over a 30 s window), while the vertical colored error bars indicate the standard deviations of the frequency fluctuations found within the running mean sampling timescale. As already shown in Figs. 27 and 28, the frequencies associated with higher IMFs are naturally lower as a result of the residual time series containing less rapid fluctuations. It can be seen the in left panel of Fig. 29 that IMF $c_{2}$ contains frequencies in the range of 50–300 mHz (3–20 s), while IMF $c_{7}$ displays lower frequencies spanning 1–30 mHz (33–1000 s). We must note that the left panel of Fig. 29 is simply a representation of the instantaneous frequencies present in the time series as a function of time and does not contain information related to their relative importance (e.g., their respective amplitudes), although this information is indeed present in the overall Hilbert–Huang transform.

Fig. 29

Instantaneous frequencies computed from applying a Hilbert–Huang transform to the HARDcam H$\alpha$ lightcurve shown in the lower-right panel of Fig. 27 and displayed as a function of time (left panel). The solid purple, blue, dark green, light green, orange, and red lines correspond to moving average frequencies (computed over a 30 s window) for the IMFs ${c_{2} \rightarrow c_{7}}$, respectively. The vertical error bars correspond to the standard deviations of frequencies found within the 30 s moving average windows. The right panel displays the corresponding Hilbert–Huang power spectrum, calculated by integrating the instantaneous frequency spectra over time and normalized to the largest power value computed, hence providing a plot of relative changes in spectral energy as a function of frequency. Features within the power spectrum are consistent with the FFT and wavelet outputs shown in Figs. 12 and 23

Finally, it is possible to integrate the instantaneous frequency spectra (including their relative amplitudes) across time, producing the Hilbert–Huang power spectrum shown in the right panel of Fig. 29. The features of the Hilbert–Huang power spectrum are very similar to those depicted in the FFT spectrum shown in the right panel of Fig. 12. Notably, there is a pronounced power enhancement at $\approx 4.7$ mHz, which is consistent with both the FFT power peaks (right panel of Fig. 12) and the heightened wave amplitudes found at $\approx 210$ s in the wavelet transform shown in the bottom panel of Fig. 23. This shows the consistency between FFT, wavelet, and EMD approaches, especially when visible wave activity is evident.

2.6 Proper orthogonal decomposition and dynamic mode decomposition

Recently, for the first time, Albidah et al. (2021, 2022) applied the methods of Proper Orthogonal Decomposition (POD; see e.g., Pearson 1901; Lumley 1967) and Dynamic Mode Decomposition (DMD; see e.g., Schmid 2010) to identify MHD wave modes in a sunspot umbrae. The POD method defines the eigenfunctions to be orthogonal in space but places no constraints on their temporal behaviour. On the other hand, DMD puts no constraints on the spatial structure of the eigenfunctions but defines them to be orthogonal in time, meaning that each DMD mode has a distinct frequency. Hence, POD modes are permitted to have broadband frequency spectra but DMD modes are not. This is shown in the right panel of Fig. 30, which shows a broadband power spectral density (PSD) of 8 POD modes detected in a sunspot umbra by Albidah et al. (2021) using HARDcam H$\alpha$ intensity observations from Jess et al. (2017).

Fig. 30

The left panel shows a sunspot from Jess et al. (2017) in H$\alpha$ intensity using data from HARDcam (one pixel has a width of 0.138”). The middle panel shows the mean intensity of the time series, the colourbar displays the magnitude of the mean time series, the solid black line shows umbra/penumbra boundary (intensity threshold level 0.85) and the green box (101 $\times$ 101 pixels) shows the region where Albidah et al. (2021) applied POD and DMD. The right panel displays the PSD of the time coefficients of the first 20 POD modes (in log scale). The PSD shows peaks between frequencies 4.3–6.5 mHz (corresponding to periods of 153–232 s). Image adapted from Albidah et al. (2021)

Both the POD and DMD produce 2D eigenfunctions as shown in the left and middle columns of Fig. 31, however, they achieve this using different approaches. Essentially, DMD identifies the spatial modes which best fit a constant sinusoidal behavior in time, as with a Fourier transform. POD ranks the spatial modes in order of contribution to the total variance, which DMD cannot do.

Since POD can produce as many modes as there are time snapshots, the challenge is to identify which modes are physical and which are not. Similarly, not all DMD modes may be physical. For practical purposes a physical model, such as the magnetic cylinder model (see Sect. 2.9.2 for discussion of MHD wave modes of a magnetic cylinder) can be used to select POD and DMD modes which most closely correspond to predicted MHD wave modes. For the approximately circular sunspot shown in Fig. 31, the predicted MHD cylinder modes which are in the strongest agreement with the selected POD and DMD modes are shown in the right column. These are the fundamental slow body sausage (top row) and kink modes (bottom row).

Fig. 31

The top and bottom rows show snapshots of the slow body sausage and kink modes, respectively. From left to right, the columns show the POD and DMD modes from HARDcam H$\alpha$ intensity observations of a sunspot (Jess et al. 2017), then the corresponding magnetic cylinder model modes. As shown in the color bar, the intensity oscillations are normalized between $-1$ and 1, hence the blue and red regions are in anti-phase. The methods of POD and DMD provide a most promising approach to decompose MHD wave modes in pores and sunspots, even if their cross-sectional shapes are much more irregular than this example. Image adapted from Albidah et al. (2021)

In the case of the magnetic cylinder model, assuming a static background, the eigenmodes, e.g., kink, sausage and fluting, are orthogonal to each in other in space by definition. Furthermore, each mode can have a broadband signal in $\omega$ and k as shown for a real sunspot in the right panel of Fig. 30. Hence, POD can identify such modes in pores and sunspots, providing there is no significant background flow that will break the condition of orthogonality. Furthermore, if a mode has a dominant power frequency, this can be identified with DMD as well as POD. Indeed, this was done by Albidah et al. (2021) for the 8 POD modes shown in the PSD plot in the right panel of Fig. 31 which have distinct power peaks between 4.3 and 6.5 mHz. In such cases a combined POD/DMD approach is a most promising avenue for identifying physical modes. However, it must be highlighted, as initially introduced in Sect. 2.2, that the characterization of waves using POD and DMD techniques must be treated with the same caution as traditional FFT approaches. For example, it is essential that the relative amplitudes of each eigenfunction are compared to noise and/or background sources to establish its true significance.

As will be shown in Sect. 3.3, POD and DMD methods are especially useful for decomposing MHD wave modes in pores and sunspots of more irregular cross-sectional shapes than the example shown in Fig. 31. This is because POD and DMD do not have the limitation of having their eigenfunctions pre-defined as they are with Fourier decomposition, where the basis functions are simply fixed as sinusoids. Even in the standard cylinder model, the eigenfunctions in the radial direction are Bessel functions not sinusoids. Hence, when it comes to identifying the spatial structure of individual MHD wave modes in pores and sunspots, the methods of POD and DMD are more suited to the job than Fourier decomposition. However, Fourier transforming the time coefficient of a POD mode is still necessary to calculate its PSD as shown in the right panel of Fig. 30.

2.7 $B{-}\omega$ diagrams

Imaging spectropolarimetry offers the additional possibility to study the variations in the wave power spectrum as a function of magnetic flux. To this aim, Stangalini et al. (2021b) have proposed a new visualization technique, called a $B{-}\omega$ diagram (see Fig. 32), which combines the power spectrum of a particular quantity (e.g., Doppler velocities) with its corresponding magnetic information. In this diagram, each column represents the average power spectrum of pixels within a particular magnetic field strength interval as inferred from polarimetry (e.g., via spectropolarimetric inversions or center-of-gravity methods; Rees and Semel 1979). The $B{-}\omega$ diagram therefore has the capability to help visualize changes in the oscillatory field as one transitions from quiet Sun pixels outside the magnetic tube to the inner (more concentrated) magnetic region. In Fig. 32 we show an example of $B{-}\omega$ diagram taken from Stangalini et al. (2021b), which reveals unique wave information for a magnetic pore observed by IBIS in the photospheric Fe i 6173 Å spectral line. Here, we clearly see that the the amplitude of 5-min ($\approx 3$ mHz) oscillations in the quiet Sun is progressively reduced as one approaches the boundary of the magnetic pore (increasing B values). On the other hand, immediately inside the boundary of the pore (highlighted using a dashed vertical line), a set of spectral features is observed in both Doppler velocity and CP (circular polarization) oscillations (i.e., magnetic field oscillations), which are interpreted as specific eigenmodes of observed magnetic cylinder.

Fig. 32

Image reproduced with permission from Stangalini et al. (2021b)

Doppler velocity (left) and circular polarization (CP; center) $B{-}\omega$ diagrams of a magnetic pore observed in the photospheric Fe i 6173 Å spectral line. The vertical blue dashed lines represents the boundary of the umbral region as inferred from intensity images. The right panel shows the average spectra outside and inside the magnetic umbra of the pore. The 5-min (p-mode) oscillations dominate the quiet Sun, but their amplitude is progressively reduced absorbed as one approaches the concentrated magnetic fields of the pore, until a series of eigenmodes are excited within the magnetic tube itself.

2.8 Effects of spatial resolution

The solar atmosphere is highly structured, presenting features across a wide range of spatial scales down to the resolution limit of current instrumentation. Oscillations can be localized at particular spatial scales/features (see, e.g., the discussion in Sect. 3.3). This means that, for instance, the Doppler velocity, or indeed any other diagnostic, is the average within the resolution angle of the observations. For this reason, the signal itself and its inherent temporal oscillations associated with features below (or close to) the resolution limit can be underestimated (MacBride et al. 2021).

To illustrate this effect, we consider a case study based on CRISP observations acquired at the SST of a quiet Sun region, which were previously deconvolved using the MOMFBD code (van Noort et al. 2005) to reduce the effects of residual image aberrations. Here, for the seek of simplicity we consider the starting data as “perfect data” for the only purpose of illustrating the effects of spatial resolution of the final power spectra of the oscillations. In the left panel of Fig. 33 we show the original instantaneous Doppler velocity field obtained from the Fe i 6301.5 Å photospheric spectral line. In order to mimic the effect of a lower spatial resolution, we convolve this data using a point spread function (PSF), assumed here to be Gaussian, with a larger full-width at half-maximum (FWHM). In order to simplify the process, we do not consider the effects of residual seeing aberrations present in the original and convolved images. Therefore, our PSF model only considers the effect of the instrumental PSF, which can be represented by the Houston diffraction limited criterion (Houston 1927),

$$\begin{aligned} {\text{FWHM}} = \frac{1.03\lambda }{D} , \end{aligned}$$

(12)

where $\lambda$ is the observed wavelength and D is the diameter of the telescope. Local seeing effects in ground-based observations can further reduce the effective resolution, in addition to the seeing conditions themselves varying significantly throughout the observations, thus providing further (time varying) degradation to the data. In the left panel of Fig. 33, the photospheric velocity field is the result of two components: downflows in the intergranular lanes (red colors) and upflows in the granules (blue colors). Since the integranular lanes are much smaller and narrower with respect to the granules, the velocity signals associated with the integranular regions become more affected (i.e., reduced) by the lower spatial resolution induced by worsening seeing conditions. This effect is apparent in the middle and right panels of Fig. 33, where the progressively worsening seeing conditions (${\text{FWHM}}=0.2''$ middle panel; ${\text{FWHM}}=0.5''$ right panel) result in lost fine-scale velocity information.

Fig. 33

Estimated effects of the spatial resolution (i.e., different FWHMs of the instrumental PSF; see Eq. (12)) on the observed Doppler velocity field. The original Doppler velocity field observed by the CRISP instrument at the SST in the Fe i 6301.5 Å photospheric spectral line (left panel) is convolved with a Gaussian PSF with larger and larger FWHMs to mimic the effects of a lower spatial resolution (middle and right panels). The sign convention employed shows downflows (positive velocities) as red colors and upflows (negative velocities) as blue colors. It can be seen in the middle (${\text{FWHM}}=0.2''$) and right (${\text{FWHM}}=0.5''$) panels that progressively worsening seeing conditions results in lost velocity signals from primarily small-scale features (e.g., intergranular lanes)

If the resolution angle is smaller than the angular size of the feature being studied, then the measured signal will approach the true value. This is due to the ‘filling factor’ being equal to ‘1’, whereby the feature of interest occupies the entirety of the resolution element on the detector. On the contrary, if the resolution element is larger than a particular spatial feature, then the signal measured will be a combination of both the feature of interest and plasma in its immediate vicinity. Here, the filling factor of the desired structure is $<1$, resulting in a blended signal that produces the measured parameters. In the specific case of integranular lanes (see, e.g., Fig. 33), this means that if the resolution element is larger than their characteristic width, signal from the neighboring granules will be collected too. This effect is shown in Fig. 34, where the probability density functions (PDFs) of the instantaneous velocities for different spatial resolutions is shown. By lowering the spatial resolution, the original skewed distribution of the velocity, which is a consequence of the different spatial scales associated with the upflows (blueshifts) and downflows (redshifts), is transitioned into a more symmetric distribution that is characterized by smaller velocity amplitudes.

Fig. 34

Probability density functions (PDFs) of the instantaneous velocity fields shown in Fig. 33 as a function of spatial resolution. Here, the blue, orange, and green lines represent the PDFs for three different seeing conditions represented by a ${\text{FWHM}}=0.16''$, ${\text{FWHM}}=0.2''$, and ${\text{FWHM}}=0.5''$, respectively. It can be seen that worse seeing conditions (e.g., the green line) produce more symmetric distributions and smaller velocity amplitudes due to the spatial averaging of the measured signals

These effects, in turn, also translate into a reduction of the measured amplitudes of any oscillations present in the data. This effect can be seen in Fig. 35, where the suppression factor of the Doppler velocity amplitudes (upper panel) and the resulting power spectral densities in two distinct frequency bands, namely 3 mHz and 5 mHz (1 mHz bandwidth; lower panel), are shown as a function of the spatial resolution. The suppression factor gives an idea of the underestimation of the amplitudes of the embedded oscillations, and in the top panel of Fig. 35 it is normalized to the value associated with the original SST/CRISP data used here (i.e., ${\text{FWHM}}=0.16''$ provides a suppression factor equal to 1.0). From the upper panel of Fig. 35 we can also predict the amplitudes of the velocity oscillations captured in forthcoming observations from the new 4 m DKIST facility, which could be as large as 1.3–1.4 times that of the velocity amplitudes measured with a 1 m class telescope at the same wavelength (under similar local seeing conditions).

Fig. 35

Wave amplitude suppression factor (upper panel) and the resulting power spectral densities (lower panel) for observations acquired with different spatial resolutions. In the upper panel, the wave amplitude suppression factors (blue dots) are computed with respect to the velocity information displayed in Fig. 33, with the vertical dotted lines highlighting telescope aperture sizes of 4 m (DKIST), 1 m (SST), and 0.1 m. The dashed red line displays an exponential fit (using Eq. (13)), with the fit parameters shown in the figure legend. The lower panel displays the resulting power spectral densities, as a function of spatial resolution, for two key frequencies commonly found in observations of the solar atmosphere, notably 2.5–3.5 mHz (orange dots) and 4.5–5.5 mHz (blue dots). Again, the power spectral densities are fitted using Eq. (13), with the corresponding fit parameters shown in the figure legend. These panels document the importance of spatial resolution when attempting to measure weak oscillatory processes, since poor spatial resolution (either through small telescope aperture sizes or poor local seeing conditions) may result in complete suppression of the observable signal

Both the suppression factor and the resulting power reduction, as a function of spatial resolution, are well modeled by an exponential decay of the form,

$$\begin{aligned} A = A_{0} e^{-\frac{{\text{FWHM}}}{{s_{0}}}} + C , \end{aligned}$$

(13)

where $A_{0}$ is either the amplitude of the velocity signals or the wave power, $s_{0}$ is a characteristic spatial length, and C is a constant. Equation (13) characterizes very nicely the impact spatial resolution has on the visible wave characteristics, whereby when the resolution element is larger than the characteristic physical scale of the observed process in the solar atmosphere (i.e., ${\text{FWHM}}>s_{0}$), then the oscillatory signal is strongly suppressed. This may result in weak oscillatory amplitudes being lost from the final data products, a process that was recently discussed by Jess et al. (2021b) in the context of sunspot oscillations.

Such amplitude suppression effects imply that when estimating the energy flux of waves, one needs to consider the specific spatial resolution achieved and correct the resulting estimates by a factor depending on the FWHM of the instrumental PSF and the local seeing effects. We note that this effect strongly depends on the characteristic spatial length of the processes observed in the solar atmosphere. In order to illustrate the problem we have made use of photospheric observations (i.e., Figs. 33, 34, 35). However, due to the presence of narrow filamentary structures observed in the chromosphere, the power of the oscillations can be even more underestimated at those atmospheric heights.

2.9 Identification of MHD wave modes

In this section, we will not review MHD wave theory in any great detail since this has been covered previously in many books and reviews (see e.g., Aschwanden 2005; Nakariakov and Verwichte 2005; Priest 2014; Jess et al. 2015; Roberts 2019). Instead, we would like to highlight the particular challenges of identifying MHD wave modes from observational data given what is known from MHD wave theory.

2.9.1 Homogeneous and unbounded plasma

In most textbooks, for simplicity MHD waves are rightly introduced by assuming a homogeneous unbounded plasma with a straight and constant magnetic field. This highly idealized plasma configuration only permits propagating Alfvén, slow, and fast magnetoacoustic wave modes. In stark contrast, the Sun’s atmosphere is actually very inhomogeneous and the newest high resolution instrumentation reveal the solar plasma to be ever more finely structured. But let us assume the wavelengths are large enough so that these MHD wave modes do not “feel” the effect of any plasma fine structure, hence allowing us to apply the unbounded homogeneous plasma model, as a zeroth order approximation, to observational data. How can we actually identify the Alfvén, slow, and fast magnetoacoustic MHD wave modes? As we shall discuss, in practical terms, even in this simplest of plasma configurations, each MHD wave mode would actually be non-trivial to identify without ambiguity, even from excellent quality spectropolarimetric data.

First, let us consider the Alfvén wave (Alfvén 1942). The only restoring force of this wave is magnetic tension, but since this wave is incompressible the magnetic field lines remain equidistant from each other as they are oscillating. Hence, although the direction of the magnetic field vectors will change with time as the field lines oscillate the magnitude of the vectors will remain constant. Therefore, this wave will not reveal itself through variations in the magnetic field strength using the Zeeman or Hanle effects. Also, due to its incompressibility the Alfvén wave would not reveal itself in intensity oscillations since the density is not perturbed. This only leaves the velocity perturbations associated with this wave, which could in principle be detected in Doppler measurements. However, to truly identify an Alfvén wave it would have to be established that the velocity perturbations were perpendicular to the magnetic field lines and that the wave vector was not perpendicular to the direction of the magnetic field. To add even more difficulty to the challenge of identifying an Alfvén wave, it is only approximately anisotropic, in the sense that the fastest propagation is along the direction of the magnetic field and only completely perpendicular propagation is forbidden, i.e., the more perpendicular the wave vector becomes relative to the magnetic field the slower the propagation will be.

What about identifying the slow and fast magnetoacoustic modes? The allowed directions for the slow magnetoacoustic wave vector are very similar to that of the Alfvén wave, meaning that it is only approximately anisotropic and propagation perpendicular to the magnetic field direction is forbidden. However, unlike the Alfvén wave, the slow magnetoacoustic wave is compressible and should reveal itself in intensity oscillations if the amplitude of the perturbations are large enough relative to the background. However, to establish even more convincing evidence, a slow magnetoacoustic wave requires validation that the plasma and magnetic pressure perturbations are in anti-phase. Of course, this is not an easy task in observational data and would require both a fortuitous line-of-sight and an excellent signal-to-noise ratio to determine perturbations in both intensity and Zeeman/Hanle effect measurements. In contrast to the Alfvén and slow magnetoacoustic waves, the fast magnetoacoustic wave is more isotropic in nature since it can also propagate perpendicular to the magnetic field. A further key difference to the slow magnetoacoustic wave is that the plasma and magnetic pressure perturbations associated with a fast magnetoacoustic wave are in phase. To show this from observational data would provide compelling evidence that a fast magnetoacoustic wave mode has indeed been identified, but, as with showing the anti-phase behavior between plasma and magnetic pressures for a slow magnetoacoustic wave, this is not a trivial task, even with excellent quality spectropolarimetric data.

There are also more subtle points in distinguishing between the Alfvén, slow, and fast magnetoacoustic wave modes depending on the value of plasma-$\beta$, which itself is difficult to determine from observational data. Importantly, for MHD wave modes the value of plasma-$\beta$ also indicates the relative values of the sound and Alfvén speeds. Especially problematic is the case when the sound speed is close to the Alfvén speed, since here the propagation speeds of the Alfvén, slow and fast magnetoacoustic waves along the direction of the magnetic field are practically indistinguishable. This effect is clearly demonstrated in Fig. 36, which is based on the ‘NC5’ flux tube model presented by Bruls and Solanki (1993), and clearly shows how the localized velocities associated with different wave modes can become difficult to disentangle in the lower solar atmosphere, hence providing some ambiguity when attempting to diagnose the true wave mode from the propagation velocity alone. But remember, the nuanced discussion we have had here on wave mode identification assumed that the solar plasma was both homogeneous and unbounded. In practical terms, it is more likely that the analysis of waves in the lower solar atmosphere will be directly related to their excitation in, and propagation through, large scale magnetic structures such as sunspot and pores (see Sect. 3.2) or smaller scale structures such as spicules and fibrils (see Sect. 3.4). In such cases the most applied model is that of the magnetic cylinder (e.g., Wentzel 1979; Wilson 1979, 1980; Spruit 1982; Edwin and Roberts 1983), which we shall discuss next.

Fig. 36

Various wave speeds in a flux tube in the lower solar atmosphere, from the hot ‘NC5’ flux tube model put forward by Bruls and Solanki (1993), in combination with the surrounding cool VAL-A atmosphere (Vernazza et al. 1981)

2.9.2 Magnetic cylinder model

The advantage of the magnetic cylinder model is that it allows for the key plasma parameters, e.g., magnetic field strength and plasma density, to differ inside and outside of the flux tube, allowing us to introduce inhomogeneity in the direction perpendicular to the cylinder axis. In this model, relative to the cylindrical coordinates $(r, \theta , z)$, where r, $\theta$, and z are the radial, azimuthal, and axial directions, respectively, waves can either be standing or propagating in all three orthogonal directions (see the left panel of Fig. 37). If the wave is propagating in the radial direction this is a so-called “leaky” wave, which is not trapped by the cylindrical waveguide and damps due to MHD radiation. The so-called “trapped” modes are standing in the radial direction with the greatest wave energy density in the internal region of the cylinder. Outside of the cylinder the trapped mode is evanescent and decays with increasing distance from the tube.

Fig. 37

Images reproduced with permission from Arregui et al. (2005, left), copyright by ESO; and Morton et al. (2012, middle and right panels), copyright by Macmillan

A typical cylindrical flux tube model (left panel) represented by a straightened magnetic tube of length, L, and radius, R. The magnetic field, B, is uniform and parallel to the z-axis and the whole configuration is invariant in the azimuthal direction, $\theta$ (labeled as $\varphi$ in the diagram). In the schematic, the density varies in a non-uniform transitional layer of width, l, from a constant internal value, $\rho _{i}$, to a constant external value in the local plasma environment, $\rho _{e}$. The middle and right panels show the effects of $m=0$ (sausage) and $m=1$ (kink) wave perturbations, respectively, to the equilibrium flux tube. The sausage wave (middle) is characterized by an axi-symmetric contraction and expansion of the tube’s cross-section. This produces a periodic compression/rarefaction of both the plasma and magnetic field. The kink wave (right) causes a transverse displacement of the flux tube. In contrast to the sausage wave, the kink wave displacement/velocity field is not axi-symmetric about the flux tube axis. The red lines show the perturbed flux tube boundary and thick arrows show the corresponding displacement vectors. The thin arrows labelled B show the direction of the background magnetic field.

Beyond the basic descriptions of whether the mode is “leaky” or “trapped”, the azimuthal integer wave number, m, defines whether the waves are the so-called “sausage”, “kink”, or “fluting” modes. The sausage mode has $m=0$ and is azimuthally symmetric, the kink mode has $m=1$ and is azimuthally asymmetric (see the middle and right panels of Fig. 37). The fluting modes are higher order in the azimuthal direction with $m \ge 2$. A further classification of wave types in a magnetic cylinder is “body” or “surface” modes. A body wave is oscillatory in the radial direction inside the tube and evanescently decaying outside. Because the body wave is oscillatory inside the tube, it has a fundamental mode in the radial direction and also higher radial harmonics. In contrast, a surface wave is evanescent inside and outside of the tube with its maximum amplitude at the boundary between the internal and external plasma. Since it is strictly evanescent inside the tube, the surface mode cannot generate higher radial harmonics.

At this point it will be worth explaining why confusion has arisen over the years since the seminal publication by Edwin and Roberts (1983), who also introduced the terms “fast” and “slow” to classify the propagation speeds of MHD wave modes along the axis of the magnetic cylinder. In the dispersion diagrams of a magnetic cylinder, distinct bands appear for a particular wave mode where the axial phase speed is bounded by characteristic background speeds. As an example, we can model a photospheric waveguide as being less dense than the surrounding plasma and having a stronger magnetic field internally than externally. This would be a reasonable basic model for, e.g., a pore or sunspot, where the internal density depletion is a result of the increased magnetic pressure (Maltby et al. 1986; Low 1992; Cho et al. 2017; Gilchrist-Millar et al. 2021; Riedl et al. 2021). In this case, we can form the inequality of the characteristic background speeds as $v_A>c_e>c_0>v_{Ae}$, where $v_A$ is the internal Alfvén speed, $c_e$ is the external sound speed, $c_0$ is the internal sound speed, and $v_{Ae}$ is the external Alfvén speed. This results in a slower band with phase speeds between $[c_T, c_0]$, where the internal tube speed, $c_T$, is defined as,

$$\begin{aligned} c_T = \frac{c_0 v_A}{\sqrt{c_0^2+v_A^2}} . \end{aligned}$$

(14)

In addition, a faster band also exists with phase speeds between $[c_0, c_e]$. Wave modes with phase speeds below the “slow” band and above the “fast” band are not trapped modes (having real $\omega$ and $k_z$ values). The “slow” and “fast” bands for these chosen photospheric conditions are shown in the dispersion diagram in the left panel of Fig. 38.

Fig. 38

Left panel: A dispersion diagram is shown for a representative photospheric magnetic cylinder. It can be seen that there are two distinct horizontal bands with slower and faster phase speeds. The fast band is bounded between $[c_0,c_e]$ and the slow band between $[c_T,c_0]$. The adjectives “slow” and “fast” here have a quite distinct meaning from the terms slow and fast when referring to the magnetoacoustic wave modes of a homogeneous and unbounded plasma. Right panel: A cartoon of theoretically predicted MHD wave modes in a sunspot, and their possible sources, based on the magnetic cylinder model of Edwin and Roberts (1983). Images adapted from Edwin and Roberts (1983, left panel) and Evans and Roberts (1990, right panel)

Although Edwin and Roberts (1983) used the perfectly apt adjectives, “slow” and “fast”, to describe the phase speed bounds of these distinct bands of trapped MHD wave modes, they have quite a different physical meaning to the terms of the slow magnetoacoustic and fast magnetoacoustic waves from the homogeneous and unbounded plasma model. This is most clearly illustrated when comparing the same label “fast” in both scenarios. For a cylindrical waveguide any trapped fast MHD mode is strictly anisotropic since the propagating wave vector is restricted to being absolutely parallel to magnetic field direction, which is also aligned with the cylinder axis. However, a fast magnetoacoustic wave in a homogeneous plasma can propagate with any angle relative to the magnetic field orientation.

There is a special class of incompressible Alfvén modes that can exist in a magnetic cylinder with any azimuthal wave number, m, the so-called torsional Alfvén waves (see e.g., Spruit 1982). Like the Alfvén wave in a homogeneous plasma, the only restoring force is magnetic tension. However, torsional Alfvén waves are strictly anisotropic since they can only propagate along the direction of the tube axis, whereas their counterpart in a homogeneous plasma can propagate at any angle (with the exception of perpendicular) relative to the magnetic field. The torsional Alfvén wave can only be excited if the driver itself is incompressible, meaning that the tube boundary is not perturbed at all in the radial direction. However, in reality it likely that the boundary of solar magnetic flux tubes are perturbed to some degree in the radial direction. If the boundary is only slightly perturbed in the radial direction, and the dominant perturbations are in the axial direction, then this will excite a slow mode. If the radial perturbation dominates over the axial perturbation, resulting in a greater perturbation of the boundary, then this will excite a fast mode. The greater radial perturbation for a fast mode means that magnetic tension plays a larger role in the restoring force than for a slow mode, where the longitudinal compressive forces of plasma and magnetic pressure dominate.

Understanding the phase relations between the restoring forces for MHD wave modes in a magnetic cylinder is not as straightforward as it is for the three possible MHD modes in a homogeneous plasma. This is because the phase relations between plasma pressure, magnetic pressure, and magnetic tension restoring forces depend on whether the wave is propagating or standing in each of the three orthogonal directions, i.e., radial (r), azimuthal $(\theta )$, and axial (z). Also, the radial spatial structuring of the plasma in a magnetic cylinder means that perturbed MHD variables, such as the magnetic field $(B_r,B_{\theta },B_z)$ and velocity $(v_r,v_{\theta },v_z)$ components, are related, not only by time derivatives, but spatial derivatives dependent on the variation of the background plasma properties.

A simplified thin tube or “wave on a string” approximation was made by Fujimura and Tsuneta (2009) to derive the phase relations between $v_r$ and $B_r$ for a kink mode, and $v_z$ and $B_z$ for a sausage mode. This was done for both propagating and standing waves in the axial direction, but caution should be taken in applying these results to structures of finite width. A more detailed investigation into the phase relations of these MHD variables was done for the sausage mode by Moreels and Van Doorsselaere (2013), utilizing a magnetic cylinder of finite width under photospheric conditions. Like Fujimura and Tsuneta (2009), this model predicted the phase relations for both standing and propagating waves in the axial direction. A note of caution should be introduced here to state that both the models of Fujimura and Tsuneta (2009) and Moreels and Van Doorsselaere (2013) assume the kink and sausage modes are “free” oscillations of the structure and are not being driven. To correctly derive the phase relations between the MHD wave variables in a driven system demands that system is solved as an initial value problem. However, currently the exact spatial and temporal structures of the underlying drivers of the waves observed in pores and sunspots are not universally understood.

Although the phase relations between the perturbed variables for any MHD wave mode may be not simple to predict theoretically, at the least the spatial structure of these variables (independent of time), providing the cross-section of the wave guide is resolved (e.g., particularly in the case of larger magnetic structures such as pores and sunspots), should correlate in straightforward way. First, let us consider a fixed axial position, z, which for a vertical tube would correspond to a fixed height in the solar atmosphere. If the magnetic cylinder is oscillating with an eigenmode, then the variables related to compressible axial motion, i.e., $v_z$, $B_z$ and plasma pressure (also related to perturbations in temperature and plasma density), should have the same spatial structure in the radial (r) and azimuthal $(\theta )$ directions. Likewise, the spatial structure of variables related to radial perturbations of the magnetic field, i.e., $v_r$ and $B_r$, should be consistent. The same is also true for the variables that relate to the torsional motions of the magnetic field, i.e., $v_{\theta }$ and $B_{\theta }$. Again, all these theoretical predictions assume free oscillations of the entire magnetic structure, e.g., a pore or sunspot. If the oscillations are being driven, then this is a more complicated and computationally expensive modeling problem to solve. Also, the spatial scale of the driver relative to the size of the magnetic structure is crucial. To excite the global eigenmodes of magnetic structures the driver has to be at least as large as the structure itself. If the driver is much smaller than the magnetic structure, it will still excite localized MHD waves, but these will not be global eigenmodes of the entire magnetic structure. This too requires a different modeling approach, see e.g., Khomenko and Collados (2006), who modeled p-mode propagation and refraction through sunspots.

High resolution images of sunspots, pores, magnetic bright points, and fibrillar structures are continually telling us that modeling these features using cylindrical flux tube geometries, while more mathematically simplistic, is far from realistic. Even from basic membrane models, in which separation of variables is possible, the cross-sectional shape has a fundamental effect on the structure of the eigenfunctions. For elliptical magnetic flux tubes, Aldhafeeri et al. (2021) investigated the effect of eccentricity, $\epsilon =\sqrt{1-b^2/a^2}$, where a and b are the semi-major and semi-minor axes, respectively, on the spatial structure of eigenfunctions. See, for example Fig. 39, which shows two sunspot umbrae fitted with ellipses with eccentricities $\epsilon =0.58$ and $\epsilon =0.76$. These are not negligible values since a circle has $\epsilon =0$. Figure 40 shows $m=1$ (kink) and $m=2,3$ (fluting) fast body modes where the phase is odd with respect to the major axis as eccentricity increases, while Fig. 41 shows the same modes where the phase is even with respect to the major axis. Although all MHD wave modes in flux tubes of elliptical cross-section have their spatial structure distorted when compared to their equivalent versions in flux tubes of circular cross-section, it can be seen that the fluting modes that have even phase with respect to the major axis (shown in Fig. 41) become notably different in character as eccentricity increases, since previously distinct regions of phase or anti-phase end up coalescing. This advancement from the cylindrical flux tube model demonstrates that more sophisticated modeling of magnetic flux tubes with more realistic, and hence more irregular, cross-sectional shapes is required to more accurately interpret what type of wave modes are present in pores and sunspots. Recently this was done by Albidah et al. (2022) and Stangalini et al. (2022) to identify MHD wave modes in sunspot umbrae and this will be discussed in Sect. 3.2.

Fig. 39

Two active regions, NOAA AR12565 (left) and NOAA AR12149 (right), captured in the G-band by ROSA at the Dunn Solar Telescope. To show the departure from circular cross-sectional shape, ellipses are fitted to the sunspot umbrae. The eccentricity of the left umbra is $\epsilon =0.76$, while the right umbra is $\epsilon =0.58$. Image adapted from Aldhafeeri et al. (2021)

Fig. 40

The normalized density perturbations of fast body wave modes under representative coronal conditions for the different values of eccentricity $\epsilon$. Note that the eigenfunctions for slow body wave modes under photospheric conditions would have a very similar appearance. From top to bottom, the $m=1$ (kink) and $m=2,3$ (fluting) modes are shown which have an odd phase structure with respect to the major axis of the ellipse. Image adapted from Aldhafeeri et al. (2021)

Fig. 41

The same wave modes are shown as in Fig. 40 but their phase structure is even with respect to the major axis of the ellipse. Image adapted from Aldhafeeri et al. (2021)

In Sect. 2.8 the crucial issue of spatial resolution was discussed. In smaller scale magnetic structures, such as off-limb spicules or on disc fibrils, it is not possible to observe the true cross-section of the wave guide (as is possible for larger on-disc features such as pores and sunspots) in order to identify eigenmodes. However, fast sausage and kink modes can still be identified in these smaller structures if the amplitude of the radial motion (i.e., transverse to the magnetic field direction) is large enough. The kink mode is the only cylinder mode which causes a transverse displacement of the axis. For smaller magnetic structures, such as fibrils, the kink mode will appear as a “swaying” motion. If the radial motion of the fast sausage mode is large enough, then this causes periodic changes in the width of the structure, which can be resolved. Wave mode identification in smaller magnetic structures is addressed in detail in Sect. 3.4. As for larger scale magnetic waveguides, where the cross-section can be resolved fully, such as in pores or sunspots, the right panel of Fig. 38 shows the wide variety of theoretically predicted MHD wave modes, including slow/fast and body/surface, that can exist in such structures based on the magnetic cylinder model of Edwin and Roberts (1983). Recent progress in the identification of such wave modes from observations is discussed in Sect. 3.3.

Across Sect. 2, we have discussed the fundamental theoretical considerations of waves manifesting in the solar atmosphere (Sect. 2.9), we have provided an overview of the techniques used to characterize them (Sects. 2.2–2.7), and summarized the challenges faced in light of variable spatial resolution (Sect. 2.8). Regardless of these challenges, over the last number of decades the solar community has overcame many obstacles, which has allowed for the successful acquisition, extraction, and identification of many different types of wave modes across a wide variety of solar features. In the following section, we will overview recent discoveries in the field of waves in the lower solar atmosphere, as well as comment on the difficulties still facing the global community in the years ahead.

3 Recent studies of waves

In the past, review articles have traditionally segregated wave activity in the solar atmosphere into a number of sub-topics based on the specific wave types and structures demonstrating the observed behavior. For example, Jess et al. (2015) divided up the review content on a feature-by-feature basis, including sections related to compressible and incompressible waveforms, which were subsequently further sub-divided into quiet Sun, magnetic network, and active region locations. However, as modern observations and modeling approaches continue to produce data sequences with ever improving spatial resolutions, placing the physical boundary between two locations becomes even more challenging. Indeed, emerging (and temporally evolving) magnetic fields often blur the boundaries between magnetic network elements, pores, proto-sunspots, and fully developed active regions. Hence, it is clear that solar complexity continues to increase with each improvement in spatial resolution made. As a result, dividing the content between previously well-defined structures becomes inappropriate, which is even more apparent now that mixed MHD waves (e.g., compressible and incompressible modes; Morton et al. 2012) are being identified in a broad spectrum of magnetic features.

Hence, for this topical review we employ just three (deliberately imprecise) sub-section headings, notably related to ‘global wave modes’, as well as ‘large-scale’ and ‘small-scale’ structures. This is to avoid repetition and confusion, and to allow the overlap between many of the observables in the Sun’s atmosphere to be discussed in a more transparent manner. Importantly, while discussing the recent developments surrounding wave activity in the lower solar atmosphere, we will attempt to pinpoint open questions that naturally arise from the cited work. We must stress that closing one research door more often than not opens two (or more) further avenues of investigation. Therefore, discussion of the challenges posed is not to discredit the cited work, but instead highlight the difficult research stepping stones facing the solar physics community over the years and decades to come.

3.1 Global wave modes

The field of helioseismology has employed long-duration data sequences (some spanning multiple continuous solar cycles; Liang et al. 2018) to uncover the internal structure and dynamics of the Sun through its global oscillation properties. Pioneering observations by Frazier (1968) suggested, for the first time, the presence of dual oscillating modes in the solar atmosphere, something that contradicted previous interpretations where the observed oscillations were simply considered to be an atmospheric response to granular impacts. It was subsequently shown that a variety of global wavenumbers could be seen in the photospheric velocity field of the C i 538 nm line (Deubner 1975). Importantly, the pioneering work of Deubner (1975) revealed clear ridges in photospheric $k{-}\omega$ power spectra, which helped to highlight, for the first time, that the ubiquitous 5-min p-mode oscillations are in-fact resonant eigenmodes of the Sun. Novel observations acquired during austral summer at the South Pole discovered global 5-min global oscillations at a wide range of horizontal wavelengths, revealing the true extent of oscillation modes associated with global solar resonances (Grec et al. 1980; Duvall and Harvey 1983). Traditionally, in the field of helioseismolgy, the Sun is considered as an approximate spherically symmetric body of self-gravitating fluid that is suspended in hydrostatic equilibrium (Christensen-Dalsgaard et al. 2000). This results in the modes of solar oscillations being interpreted as resonant vibrations, which can be represented as the product of a function of radius and a spherical harmonic, $Y_{l}^{m}(\theta , \phi )$. Here, l relates to the horizontal scale on each spherical shell (commonly referred to as the ‘angular degree’), while m determines the number of nodes in solar longitude (commonly referred to as the ‘azimuthal order’).

The specific modes of oscillation can be divided up into three main categories: (1) Pressure modes (p-modes), which are essentially acoustic waves where the dominant restoring force is pressure, providing frequencies in the range of $\sim 1{-}5$ mHz and angular degrees spanning $0 \le l \le 10^{3}$ (Rhodes et al. 1997; Kosovichev 2011; Korzennik et al. 2013), (2) Internal gravity modes (g-modes), where the restoring force is predominantly buoyancy (hence linked to the magnitude of local gravitational forces), which typically manifest in convectively stable regions, such as the radiative interior and/or the solar atmosphere itself (Siegel and Roth 2014), and (3) Surface gravity modes (f-modes), which have high angular degrees and are analogous to surface waves in deep water since they obey a dispersion relation that is independent of the stratification in the solar atmosphere (Mole et al. 2007). In the limit that the wavelength is much smaller than the solar radius these wave are highly incompressible. The main restoring force for f-modes is gravity, which acts to resist wrinkling of the Sun’s surface.

The intricacies of helioseismology become even more complex once isolated magnetic features, such as sunspots, develop within the solar atmosphere, which impact the velocities and travel times of the embedded global wave modes (Braun and Lindsey 2000; Rajaguru et al. 2001, 2004, 2019; Kosovichev 2012; Schunker et al. 2013, 2016). A complete overview of the progress in helioseismology is beyond the scope of the present review. Instead, we refer the reader to the vast assortment of review articles that focus on the widespread development of helioseismology over the last few decades (e.g., Deubner 1983; Bonnet 1983; Christensen-Dalsgaard 2002; Gizon and Birch 2005; Gizon et al. 2010, 2017; Gough 2013; Basu 2016; Buldgen et al. 2019)

Importantly, the magnetic field in the solar photosphere is inhomogeneous and found in discrete concentrations across all spatial scales (Zwaan 1987). Outside of the magnetic concentrations, where plasma pressure and gravity are the dominant restoring forces, longitudinal acoustic waves (i.e., p-modes) are generated at the top of the convection zone from the turbulent motions constituting the convective motion (Stein 1967; Goldreich and Kumar 1990; Bogdan et al. 1993). The p-modes can propagate upwards and contribute to heating of the higher layers if their frequency is larger than the acoustic cut-off frequency (Ulmschneider 1971b; Wang et al. 1995). Thus, the acoustic waves can dissipate their energy in the solar chromosphere by forming shocks (as a result of gas-density decreases with height), which are manifested in intensity images as, e.g., intense brightenings (Rutten and Uitenbroek 1991; Carlsson and Stein 1997; Beck et al. 2008; Eklund et al. 2020, 2021), or drive, e.g., Type i spicules, in the so-called ‘magnetic portals’ (Jefferies et al. 2006). Moreover, Skogsrud et al. (2016) showed that these shocks are associated with dynamic fibrils in an active region they exploited from observations with SST and the Interface Region Imaging Spectrograph (IRIS; De Pontieu et al. 2014b) space telescope.

Properties of propagating acoustic/magnetoacoustic waves through the lower solar atmosphere have also been reported in a number of recent studies from both ground-based (e.g., Sobotka et al. 2016; Abbasvand et al. 2020a, b) and space-borne observations (e.g., Martínez-Sykora et al. 2015; Zhao et al. 2016; Abbasvand et al. 2021). From IRIS observations of a quiet-Sun area in Mn i 2801.25 Å, Mg ii k 2796.35 Å, and C ii 1334.53 Å spectral lines (sampling the photosphere, chromosphere, and transition region, respectively), Kayshap et al. (2018) found upwardly propagating p-modes with periods on the order of 1.6–4.0 min, and downward propagation in the higher period regime (i.e., periods larger than $\approx 4.5$ min). Furthermore, Kayshap et al. (2020) identified the propagation of slow magnetoacoustic waves (with 2–9 min periodicities), within a plage region, from the high-photosphere/low-chromosphere to the transition region, using SDO/AIA and IRIS observations.

In addition, g-modes (i.e., internal gravity waves) can be produced within turbulent convective flows (Mihalas and Toomre 1981, 1982) and propagate through the lower solar atmosphere, with frequencies shorter than $\approx 2$ mHz (Newington and Cally 2010; Kneer and Bello González 2011; Vigeesh et al. 2017; Jefferies et al. 2019; Vigeesh and Roth 2020). Their identifications in the solar atmosphere have, however, been in a challenging task since they become evanescent in the convection zone and their amplitudes at the surface are exceedingly small (Schunker et al. 2018; Calchetti et al. 2021). The internal gravity waves can potentially carry a large amount of energy flux (of $\approx 5$ kW/m$^2$; Straus et al. 2008) to the chromosphere, thus contributing to its radiative losses (Vigeesh et al. 2021).

Fortunately, unlike g-modes, f-modes (i.e., surface gravity modes) have been detected in abundance and have provided valuable diagnostic information about flows and magnetic field in the near surface region (Ghosh et al. 1995; Rosenthal and Christensen-Dalsgaard 1995; Christensen-Dalsgaard et al. 1996). Furthermore, f-modes have been exploited to quantify the Sun’s effective seismic (or acoustic) radius (Schou et al. 1997; Antia et al. 2000; Dziembowski et al. 2001; Dziembowski and Goode 2005), a relatively new concept driven by results from helioseismology, as opposed to the measuring the Sun’s physical (or true) radius. Such studies have shown that f-mode frequencies, as well as being sensitive to the seismic radius, are also modified by changes in the magnetic field during the solar cycle.

3.1.1 Global p-modes in the lower solar atmosphere

Acoustic waves (i.e., p-modes) can propagate both outside and inside magnetic concentrations. Through a number of studies prior to the turn of the century, properties of their ‘global’ oscillations (i.e., properties averaged over a relatively large field of view) became a “basic fact”, describing the characteristic periodicity of p-modes as 5 min in the solar photosphere (Leighton et al. 1962; Ulrich 1970; Ruiz Cobo et al. 1997; Schunker et al. 2009), and 3 min in the chromosphere (Evans et al. 1963; Orrall 1966; Cram 1978; Fleck and Schmitz 1991; Rutten and Uitenbroek 1991; Carlsson and Stein 1992). Standing acoustic waves have also been reported from multi-line observations in the solar chromosphere (Fleck et al. 1994a), though wave patterns (and power spectra) were found to be somewhat different in He i 1080 Å observations (Fleck et al. 1994b, 1995), compared to those in other chromospheric diagnostics (e.g., Ca ii H  & K, Ca ii 8542 Å, and H$\alpha$). While the global p-modes are more purely acoustic in nature in the photosphere, they are more likely to manifest as magnetoacoustic waves in the upper atmosphere, where the magnetic forces dominate (Khomenko and Calvo Santamaria 2013).

Many of the observations demonstrating the characteristic periodicities of the global p-modes have been based on wide-band imaging at low-spatial resolutions, with very large fields of view. A recent study by Fleck et al. (2021) highlighted the presence of ubiquitous 3-min characteristic periodicities by exploring several advanced state-of-the-art numerical models. Even so, considerable differences between the various simulations were also reported, including the height dependence of wave power, in particular for high-frequency waves, varying by up to two orders of magnitude between the models (Fleck et al. 2021). Thus, although the numerical simulations provide us with important information regarding the physical processes embedded within observational data, they should be interpreted with caution since the numerical domains are too small to resolve the true physics driving large-scale global eigenmodes.

Development of modern instruments in recent years, resulting in relatively narrow-band (often spectrally-resolved) observations at high resolution, have further explored the highly dynamic nature of the lower solar atmosphere. These novel observations reveal that the physical properties and structure of the lower solar atmosphere may significantly vary over different solar regions (with different levels of magnetic flux and/or topology), as well as through different atmospheric layers. Therefore, chromospheric wide-band filtergrams, that often integrate over a significant portion of strong chromospheric lines (hence, sampling across a large range of heights), may result in mixing (or averaging) of observable information (e.g., the oscillatory power), which can largely vary within a short distance in the lower solar atmosphere. A large variation in the height of formation can also cause a strong temporal modulation that may consequently destroy the oscillatory signal. Furthermore, the effect of spatial resolution can be crucial, as information may be lost in lower resolution observations due to, e.g., smearing (see Sect. 2.8 for more discussion related to resolution effects). Moreover, an average power spectrum over a very large field of view can predominantly be dominated by characteristics of quiet-Sun regions (which cover the majority of the solar surface at any given time).

An example of the influence of spatial resolution is the larger (total) energy flux of acoustic waves (larger by a factor of $\approx 2$) found in a quiet-Sun region by the 1 m Sunrise telescope (Bello González et al. 2010b), compared to that from the 0.7m VTT telescope (Bello González et al. 2009). However, the effect of seeing-free observations with Sunrise could also play a role in that difference, highlighting again the importance of spatial resolution (as discussed in Sect. 2.8). Such variations in atmospheric seeing (that directly affect the spatial resolution achievable) can influence the measured periodicities, in particular the global p-modes that are ubiquitously visible across the photosphere and chromosphere.

In the presence of strong magnetic fields (e.g., in network or plage regions, where a group of concentrated small-scale magnetic features reside), the global acoustic power is enhanced at photospheric and low chromospheric heights (known as a ‘power halo’; Brown et al. 1992; Kontogiannis et al. 2010; Rajaguru et al. 2013), while it is suppressed in the high chromosphere (so-called ‘magnetic shadows’; Leighton et al. 1962; Title et al. 1992; McIntosh and Judge 2001). While the exact mechanisms behind such power variation are not yet fully understood, a number of suggestions have been provided in recent years, from both observations and simulations. In particular, models have shown that the power enhancement at lower heights can be due to the reflection of fast waves at the magnetic canopy, as a result of a large Alfvén speed gradient (Khomenko and Cally 2012; Rijs et al. 2016). From observations, both magnetic-field strength and inclination have been found to play an important role, with greater power in the stronger and more horizontal fields (Schunker and Braun 2011; Rajaguru et al. 2019). The power suppression of the acoustic waves in the high chromosphere has been suggested to be due to the mode conversion at the plasma-$\beta \approx 1$ level (i.e., as a result of interactions between p-mode oscillations and the embedded magnetic fields; Moretti et al. 2007; Nutto et al. 2012a), less efficient wave propagation under the canopy, or the wave-energy dissipation before it reaches the canopy (Ulmschneider 1971a; Ulmschneider et al. 2005; Song 2017; Martínez-Sykora et al. 2020; Srivastava et al. 2021). The power suppression and its spatial scale has found to be directly correlated with the magnetic-field strength and/or geometric height (Chitta et al. 2012a; Jain et al. 2014; Krishna Prasad et al. 2016). From MHD simulations with the Bifrost code (Gudiksen et al. 2011), Heggland et al. (2011) showed that field inclination plays an important role in propagation of long-period waves (longer than 3 min) in the solar chromosphere. As such, they primarily found 3-min periodicities in regions with weak or vertical magnetic fields (including the center of strong flux tubes), whereas 5-min dominant waves in strong or inclined magnetic fields (such as the edges of flux tubes).

Power suppression of 3 min oscillations in the upper solar chromosphere has been reported by Samanta et al. (2016), where almost no oscillatory power at this period was observed in time series of H$\alpha$ line-core intensity images from SST/CRISP. The authors, however, found a slightly larger number of pixels demonstrating 3 min oscillations in Doppler velocity signatures of the same spectral line. In addition, they found power halos at lower atmospheric heights. In this study, the presence of ubiquitous chromospheric transient events (i.e., short-lived fibrillar structures) was speculated to be responsible for the power enhancements at lower heights. In addition, it was speculated that mode conversion was causing the magnetic shadows found around the 3 min periodicity in the upper chromosphere. Figure 42 illustrates the multi-height observations studied by Samanta et al. (2016) (on the left) along with their corresponding distribution of dominant periods of the oscillations (on the right), representing periods corresponding to the maximum power at each pixel. The lack of 3-min oscillations (i.e., the green color) on the top layer is evident. Using high-resolution H$\alpha$ line-core observations with SST, De Pontieu et al. (2007a) had found longer periods in regions where the field is supposedly more inclined. From spatial distribution of dominant periods (from a wavelet analysis) they showed that while sunspots and plage regions were dominated by 3-min global p-modes, 5-min and longer periodicities were found in adjacent to the dense plage regions and in more inclined-field areas, respectively. We should, however, note that such dominant-period maps demonstrated by Samanta et al. (2016) and De Pontieu et al. (2007a) should be interpreted with great caution, since multiple peaks with comparable (or even equal) power may co-exist in a power spectrum. As such, the period associated to one absolute maximum of the power may not solely be representative of the oscillations in that pixel. In addition, we should note that the global wavelet spectrum is often considered a biased estimation of the true Fourier spectrum, with variable frequency resolution through the entire spectrum, that can also depend on the choice of wavelet function (see Sect. 2.4 for more details).

Fig. 42

Images reproduced with permission from Samanta et al. (2016), copyright by AAS

Multi-layer observations of a quiet-Sun region from the low photosphere to the high chromosphere (left) whose dominant oscillatory periods are shown on the right. The green, red, and yellow colors in the dominant-period maps roughly represent periods around 3, 5, and 7 min, respectively. The bottom panels illustrate the corresponding line-of-sight magnetogram, from Stokes inversions of Fe i 630.2 nm spectral line.

A recent investigation of such global oscillations (in brightness temperature) from millimeter observations with ALMA also revealed the lack of 3-min oscillations in the solar chromosphere in datasets with relatively large amounts of magnetic flux (Jafarzadeh et al. 2021). Conversely, the same study showed the presence of dominating 3-min oscillations in the most magnetically quiescent datasets employed. However, due to the uncertain nature of those millimeter observations, particularly, their exact heights of formation, further investigations are required. Furthermore, Norton et al. (2021) reported on global oscillations in the photosphere, from SDO/HMI data, in various regions, namely, the quiet-Sun, plage, umbra and the polarity inversion line of an active region. While the 5-min periodicity, with a considerably large power, was found in all four areas in Doppler velocity perturbations, much smaller power enhancements could be observed in intensity and line-width observations of the quiet and plage regions.

Of particular importance is also the effect of magnetic topology in the chromosphere, with the multi-layer magnetic canopy whose strength and thickness depends on, e.g., the magnetic flux involved in their generation (Jafarzadeh et al. 2017a). By exploring the formation and properties of various chromospheric diagnostics, Rutten (2017) showed that the dense canopies of long opaque fibrils in the upper chromosphere, seen in H$\alpha$ line-core intensity images, could act as an ‘umbrella’, obscuring the dynamics underneath. Thus, this could perhaps explain the lack of 3-min oscillations in the high chromosphere (in addition to the magnetic shadows effect). In the case of ALMA observations, Rutten (2017) speculated that the same phenomena could also occur, though at those wavelengths the dense fibrillar structures might not be visible due to their reduced lateral contrast (i.e., an insensitivity to Doppler shifts; ALMA observes continuum emission, and as such cannot be used to derive Doppler velocities) We note that similarities between ALMA observations (at 3 mm) and H$\alpha$ line-width images have been shown by Molnar et al. (2019).

All in all, it is important to investigate the variation of the global p-modes with height, throughout the lower solar atmosphere, in greater details. This can hopefully clarify whether the characteristic periodicity reported in previous studies is constant through the photosphere and the chromosphere, or whether they vary with height and/or in various solar regions.

3.2 Large-scale magnetic structures

Large magnetic structures, in the form of sunspots and solar pores, are considered ideal laboratories for the study of the excitation and propagation of MHD waves. Modern high resolution observations have revealed an extremely complex physical scenario in which different wave modes simultaneously co-exist in the same magnetic structure, hampering an unambiguous identification of individual wave modes. This is even more the case for highly structured magnetic fields, where the wave propagation reflects the geometrical complexity of the field lines acting as waveguides. However, in recent years our understanding of MHD waves in large magnetic structures, and their corresponding role in the heating of the solar atmosphere, has dramatically changed thanks to the opportunity offered by high-resolution fast-cadence tomographic imaging and to the new spectropolarimetric diagnostic capabilities, which have progressively extended up to chromospheric heights thanks to the technological advances of modern instrumentation. In particular, the inference of the plasma and magnetic field parameters obtained by spectropolarimetric inversion techniques have enabled the investigation of the effects of the magnetic field geometry on the wave propagation itself. In addition, new spectropolarimetric diagnostics have started to provide additional information about the magnetic field fluctuations, which are expected from several MHD wave modes.

The first oscillatory phenomena in the umbra of sunspots were observed by Beckers and Tallant (1969) and Beckers and Schultz (1972), where the spatially localized brightenings (so-called ‘umbral flashes’) were immediately associated with locally excited magnetoacoustic waves propagating upwards along the field lines. Today, after more than 50 years from the first discovery of these oscillations, our view of wave excitation and propagation of MHD waves in large magnetic structures has changed dramatically. From the observational point of view, in addition to the localized wave phenomena and disturbances in sunspots and pores (e.g., umbral flashes), the aforementioned instrumental advances have also allowed the identification of global eigenmodes of the magnetic structure (e.g., sausage modes, kink modes). These are generally mixed with other local disturbances, requiring specific filtering techniques (e.g., $k{-}\omega$ filtering; see Sect. 2.3) for their identification. Although most of the literature on the subject mainly reflects this apparent dichotomy with the two classes of waves (i.e., global and local oscillations) addressed independently, most recent observations have started suggesting a superposition of locally excited magnetoacoustic waves resulting, for example, from p-mode absorption or residual local magneto-convection (Krishna Prasad et al. 2015), and global resonances of the magnetic structure. These two components can coexist, both of which contribute to the physical complexity of the observed velocity patterns; an aspect that was highlighted by Roberts (2019). In the following we wish to strike a balance between local disturbances and global resonances, and we will summarize the results from the most relevant studies in recent literature.

3.2.1 Magnetoacoustic waves in large-scale magnetic structures

Sunspots and other large magnetic structures, such as solar pores, typically display intensity and Doppler velocity power spectra that are dominated by 5-min ($\sim 3$ mHz) oscillations in the photosphere, and 3-min ($\sim 5$ mHz) oscillations in the chromosphere (see for instance, Centeno et al. 2006b, 2009; Felipe et al. 2010; Felipe 2020; Felipe and Sangeetha 2020, and references therein). Of course, it must be noted that the frequencies/periodicities found at photospheric and chromospheric heights are not universal values at precisely 3 mHz and 5 mHz, respectively. Indeed, windows of power are normally referred to when discussing the corresponding Fourier spectra (Centeno et al. 2006b; Heggland et al. 2011; Gupta et al. 2013; Khomenko and Collados 2015), for example, $5\pm 0.5$ min (3.0–3.7 mHz) and $3\pm 0.5$ min (4.8–6.7 mHz) for the photosphere and chromosphere, respectively. These spectral features are depicted in Fig. 43, which clearly shows the frequency transition of peak power between the photospheric and chromospheric layers of two sunspots. While some authors have interpreted this to be the combined action of an acoustic cut-off ($\omega _{c} \sim 5.3$ mHz, allowing the upward propagation of magnetoacoustic waves with $\omega > \omega _{c}$; Deubner and Gough 1984; Duvall et al. 1991; Fossat et al. 1992; Vorontsov et al. 1998) and the atmospheric density stratification resulting in the subsequent amplification of the wave amplitudes with height, others have explained the spectral features as the result of the presence of an acoustic resonator (Jess et al. 2020, 2021b; Felipe et al. 2020). Power spectra similar to those shown in Fig. 43 were also obtained by Kanoh et al. (2016) from observations of a sunspot with Hinode/SP (in Fe i 6301.5/6302.5 Å) and IRIS (in Si iv 1403 Å), corresponding to photospheric and transition-region heights, respectively. By comparing energy fluxes at the two atmospheric regions, Kanoh et al. (2016) speculated the three orders of magnitude energy decrease with height could suggest wave dissipation in the chromosphere.

Fig. 43

Image reproduced with permission from Centeno et al. (2006b), copyright by AAS

Average umbral power spectra for two different sunspots. Solid lines indicate the power spectra of the chromospheric velocity oscillations averaged over each entire umbra, demonstrating a peak around 6 mHz. Dashed lines reveal the photospheric velocity power spectra averaged over each entire umbra, with a peak around 3.3 mHz and secondary peaks around 6 mHz.

When considering the umbrae of sunspots, the aforementioned umbral flashes dominated early observations after their initial detection. Umbral flashes (UFs) initially manifested as intensity brightenings in the core of the Ca ii K spectral line, subtending multiple arcseconds across an umbra (Beckers and Tallant 1969). Consistent with the dominant chromospheric frequencies mentioned above, UFs exhibited a 3-min periodicity, however observed brightness increases of up to 150% (Bogdan 2000) and line-of-sight velocity excursions of 10 km s$^{-1}$ (Beckers and Schultz 1972; Phillis 1975) implied that these were not the signatures of linear MHD magneto-acoustic oscillations. Subsequent modelling efforts (see the seminal works of Carlsson and Stein 1997; Bard and Carlsson 2010) established that UFs were the signature of shocks formed from the steepening of slow magneto-acoustic waves as they propagate through the large negative density gradients of the low chromosphere. When these non-linear shock fronts are formed, the intensity brightenings correspond to the dissipation of wave energy into plasma heat. At this stage, the plasma is no longer frozen into the magnetic field, and can propagate isotropically, however as the shocked plasma radiatively cools, gravitational effects will cause this overdense plasma to infall. The observational signatures of this morphology can be seen in Fig. 44, with periodic 3-min intensity brightenings seen in concert with large velocity excursions consistent with the steepening of slow modes. The right panel of Fig. 44 details the development of the shocked plasma, with the impulsive shock formation process, characterized by a notable blue-shifted velocity, leading to a more gradual red-shifted signature due to the infall of the plasma as it cools. This spectral morophology, known as a ‘saw-tooth’ is distinctive in comparison to the sinusoidal behavior of linear MHD waves.

Fig. 44

Image reproduced with permission from Grant et al. (2018), copyright by Macmillan

A velocity–time graph extracted from IBIS observations of the sunspot umbral core on 24th August 2014 by Grant et al. (2018). The horizontal axis represents Doppler line-of-sight velocity shifts from the rest wavelength, with the brightnesses displayed correlating to the Ca ii 8542Å spectral profile of a single umbral pixel over the full time series. The left panel displays velocities of up to 30 km s$^{-1}$ (or 0.85Å), while the right panel zooms in to a smaller sub-set for a closer examination of the associated signatures. The red and green lines in the right panel denote the accelerations associated with the rising (blue-shifted) and falling (red-shifted) plasma, respectively.

The nature of shock development in the solar atmosphere is deserving of its own dedicated review, as the three characteristic MHD wave speeds lead to a plethora of potential shock configurations (Delmont and Keppens 2011). Wave activity is also not the only driver, with magnetic reconnection capable of generating a range of shocks (Petschek 1964; Yamada et al. 2010). The physical processes involved in shock dynamics entails that current modeling work still strives to replicate their behavior in realistic conditions (e.g., Snow and Hillier 2019; Snow et al. 2021) and there have only been initial detections of other modes in the magnetic solar atmosphere (Grant et al. 2018; Houston et al. 2020). In the context of UFs, it is more instructive to consider them as a dissipative process of waves, as opposed to wave behavior synonymous with the focus of this review. It is, however, useful to outline recent studies that characterize the effect of UFs on umbral plasma and their effectiveness as wave dissipators. It is also valuable to discuss the effect UFs have on observables, and their influence on attempts to extract linear MHD modes from sunspot umbrae.

At the turn of the 21st century, as instrumental capabilities took a leap forward, initial studies into UF atmospheres still could not resolve any influence on the umbral magnetic field from shock fronts (Rouppe van der Voort et al. 2003). Instead, the temperature enhancements of UFs were characterized by de la Cruz Rodríguez et al. (2013) by applying NICOLE inversion techniques on Ca ii 8542 Å data, with temperature excursions of up to 1000 K inferred. However, magnetic field perturbations were still unresolvable, likely due to the coarse spatial sampling of the data and small sample size of profiles inverted. The modification of the umbral magnetic field due to UFs was finally detected by Houston et al. (2018) using polarimetric He i 10830 Å observations. By sampling this high-chromospheric spectral line with the high spectral resolution of the FIRS instrument, HAZEL inversions (Asensio Ramos et al. 2008) revealed $\sim 200$ G fluctuations in the vector magnetic field, and incremental changes in the inclination and azimuth of the field, approximately $8^{\circ }$, implying that the magnetic field enhancement is predominantly along the direction of wave and shock propagation. Houston et al. (2018) also corroborated the temperature enhancements of shocks, though with smaller average values of $\sim 500$ K, consistent with the higher atmospheric height sampled by He i 10830 Å, leading to observation of the shocked plasma as it enters into its cooling stage. Subsequently, magnetic field changes in Ca ii 8542 Å were reported (Joshi and de la Cruz Rodríguez 2018), and the subsequent derivation of semi-empirical UF atmospheric models performed by Bose et al. (2019).

These results further reinforce the scenario where UFs perturb the magnetic field geometry of the umbra. However, the perturbation is always predominantly along the magnetic field vector, as the shock propagates, and the field returns to its unperturbed state once the shock has propagated through. Thus, shocks are not seen as candidates to further incline umbral fields to permit longer period waves to pass, or to greatly impact on adjacent waves as they propagate. Indeed, the perceived propagation of UF shock fronts horizontally across the umbra towards the penumbral boundary was instead interpreted as successive UFs developing along more inclined fields (Madsen et al. 2015), further implying that the shocked plasma does not play a notable role in umbral morphology. Despite this, UFs have proved valuable in uncovering fine-scale umbral features and waves. Henriques et al. (2017) utilized the brightenings of UFs to reveal small-scale horizontal magnetic fields across the umbra, revealing a complex ‘corrugated’ structure to the field geometry in the chromosphere. UFs have been shown to generate a number of plasma flows with wave implications, notably Henriques et al. (2020) detected downflows, upflows, and counter-flows before, after, and during the UFs, respectively. Recently, downflowing UFs have been found to be a signature of standing oscillations above sunspot umbrae (Felipe et al. 2021).

When considering the processes necessary to balance the chromospheric energy budget, shocks provide a macroscopic method for converting wave energy directly into local plasma heating. As discussed earlier, the intensity excursions of UFs are confirmed as signatures of heating, with between 500 and 1000 K temperature increases observed as a result. When UFs are judged in terms of sole heaters of the chromosphere, Anan et al. (2019) derived a UF shock heating energy per unit mass of plasma that was insufficient to balance radiative losses. It is unsurprising that UFs alone are incapable of heating the chromosphere, particularly as they only occur in localized umbrae. However, the identification that slow-mode shocks impart heating energy is notable, given that it is proposed that such wave-driven shocks occur on a variety of scales across the solar atmosphere (Snow et al. 2021), they present a viable method to potentially contribute to heating. In addition, as has been discussed, shock formation is not limited to slow-mode waves. Grant et al. (2018) observed shocks at the umbra-penumbra boundary of a sunspot that are inconsistent with the scenario of UF formation, as the inclined magnetic fields in this region does not produce a large density gradient to steepen magneto-acoustic waves. Instead, it was proposed that Alfvén waves were coupling to, and resonantly amplifying, magneto-acoustic waves in at the penumbral boundary to allow for shock formation, which was verified through the transverse velocity signatures in the shocks. These Alfvén-induced shocks produced local temperature enhancements of 5%, less than simultaneous UFs, but providing dissipation of Alfvén waves in the chromosphere. This further highlights that the range of possible shock configurations are capable of dissipating a wide assortment of waves, including the elusive incompressible Alfvén mode, inferring that there is an tapestry of shock heating across the atmosphere that is yet to be fully uncovered.

When seeking to observe waves in sunspot umbrae, UFs must always be taken into account. The intensity excursions associated with UFs have noticable effects on the spectral profiles sensitive to the density and temperature perturbations of UFs, such as Ca ii H/K, Ca ii 8542 Å, He i 10830 Å, and upper chromospheric/transition region channels from IRIS, such as C ii 1335.71 Å, Mg ii 2796.35 Å, and Si iv 1393.76 Å (Tian et al. 2014; Kayshap et al. 2021). The line-core emission from these brightenings causes non-trivial profile shapes to develop, and introduces opacity effects that inhibits velocity inference through profile fitting (as seen by classes 3 and 4 of Fig. 45). A method of recovering velocity information from these profiles is to use inversion methods such as NICOLE and HAZEL, however, these are computationally intensive, time consuming, and do not account for multi-component atmospheres. The seminal work of Socas-Navarro et al. (2000) showed that in the observation of a column of shocked plasma, there is always an intermixing of active and quiescent atmospheres below the resolution limit of the observations. Thus, an observed spectral profile will in-fact be the result of a two-component atmosphere, where the filling factor is unknown. It is therefore possible to extract the linear oscillations embedded in a bright pixel if the atmospheres can be separated. This was investigated by MacBride et al. (2021) and MacBride and Jess (2021), who employed machine learning techniques to develop a neural network capable of classifying Ca ii 8542 Å as a function of line-core emission (see Fig. 45). The authors were then able to separate the two individual atmospheres through model fitting to provide values for both the shocked plasma velocity and the associated quiescent component. For any observer looking into waves in the umbra of a sunspot, care must be taken to account for UF signatures, either through exclusion of these signatures, two-component fitting, or the use of spectral lines that are less sensitive to temperature changes, such as H$\alpha$ (Cauzzi et al. 2008, 2009).

Fig. 45

Image reproduced with permission from MacBride et al. (2021)

Plots of stacked Ca ii 8542 Å umbral line spectra grouped by the neural network classification of MacBride et al. (2021), where the intensity scale for each spectrum is normalized between ‘0’ and ‘1’ to aid visualization. A two-dimensional map (lower right) reveals the prominent neural network classifications for the Ca ii 8542 Å spectra present spatially across the umbra for a single IBIS spectral imaging scan.

The propagation of waves from the umbrae outward (i.e., along penumbral filaments) are known as running (penumbral) waves (RPWs; Giovanelli 1972; Zirin and Stein 1972; Brisken and Zirin 1997; Christopoulou et al. 2000, 2001; Georgakilas et al. 2000; Kobanov and Makarchik 2004; Bogdan and Judge 2006; Bloomfield et al. 2007; Sych et al. 2009; Jess et al. 2013; Madsen et al. 2015; Löhner-Böttcher and Bello González 2015; Löhner-Böttcher et al. 2016; Stangalini et al. 2018), which have also been attributed to magnetoacoustic wave modes (Brisken and Zirin 1997; Kobanov and Makarchik 2004). RPW phenomena are mostly prominent in the mid-to-upper chromosphere, though they have also been observed at photospheric heights (Löhner-Böttcher and Bello González 2015; Zhao et al. 2015).

The origin of RPWs has long been debated as either a chromospheric phenomenon visible as a result of trans-sunspot wave interactions (e.g., Alissandrakis et al. 1992; Tsiropoula et al. 1996, 2000; Tziotziou et al. 2006; Bogdan and Judge 2006; Sharma et al. 2017a; Zhou and Liang 2017; Priya et al. 2018), or as the chromospheric signature of upwardly propagating and magnetically guided p-mode waves from the sub-photospheric layers (e.g., Christopoulou et al. 2000, 2001; Georgakilas et al. 2000; Rouppe van der Voort et al. 2003; Bloomfield et al. 2007; Reznikova and Shibasaki 2012; Reznikova et al. 2012; Jess et al. 2013; Yuan et al. 2014; Madsen et al. 2015). When viewed as a function of radial distance from the umbral center, RPW signatures manifest with large apparent phase speeds ($\sim 40$ km/s) and relatively high frequencies ($\sim 5$ mHz) at the umbra/penumbra boundary, decreasing to lower apparent phase speeds ($\sim 10$ km/s) and reduced frequencies ($\sim 1$ mHz) towards the outer penumbral edge (Kobanov et al. 2006). This effect can be also seen in Figs. 46 and 47, which are reproduced from Jess et al. (2013). Here, both the amplitude and frequency of the captured wave modes is found to depend strongly on the magnetic field geometry at chromospheric heights. The dominant frequency of the waves progressively extends towards lower values (longer periods) as one moves from the center of the umbra and into the surrounding regions with more heavily inclined magnetic fields (see the discussions below involving the ramp effect).

Fig. 46

Image reproduced with permission from Jess et al. (2013), copyright by AAS

Simultaneous images of the blue continuum (photosphere; upper left) and H$\alpha$ core (chromosphere; upper middle) acquired by the DST using the ROSA imaging instrument. A white cross marks the center of the sunspot umbra, while a white dashed line in the continuum image displays the extent of the photospheric plasma-$\beta = 1$ isocontour. The white concentric circles overlaid on the chromospheric image depict a sample annulus used to extract wave characteristics as a function of distance from the center of the umbra, while the solid white line extending into the north quadrant reveals the slice position used for the time-distance analysis displayed in Fig. 47. The dashed white lines isolate the active region into four distinct regions, corresponding to the North (N), South (S), East (E), and West (W) quadrants. The scale is in heliocentric coordinates where $1'' \approx 725$ km. The remaining panels display a series of chromospheric power maps extracted through Fourier analysis of the H$\alpha$ time series, indicating the locations of high oscillatory power (white) with periodicities equal to 180, 300, 420, and 540 s. As the period of the wave becomes longer, it is clear that the location of peak power expands radially away from the center of the umbra. This effect is synonymous with the presence of running penumbral waves (RPWs), which were first identified in solar images by Giovanelli (1972) and Zirin and Stein (1972).

Fig. 47

Image reproduced with permission from Jess et al. (2013), copyright by AAS

Top: azimuthally averaged absolute Fourier power displayed as a function of radial distance from the center of the umbra. Middle: power spectra from the top panel normalized by the average power for that periodicity within the entire field of view. Thus, the vertical axis represents a factor of how much each period displays power above its spatially and temporally averaged background. Bottom: power spectra normalized to their own respective maxima. The vertical dashed lines represent the radial extent of the umbral and penumbral boundaries, while the graduated color spectrum, displayed in the color bar at the top, assigns display colors to a series of increasing periodicities between 45 and 1200 s.

Madsen et al. (2015) examined datasets from both the SDO and IRIS spacecrafts and concluded that the apparent trans-sunspot motion associated with RPWs is not a real effect, but instead is a result of the waves (originating from the photospheric p-modes) traveling along magnetic field lines of increasing inclination angle away from the umbral core. On the other hand, Priya et al. (2018) examined high resolution observations from the Goode Solar Telescope (GST; Cao et al. 2010). The authors found that oscillatory events in the sunspot umbra appeared to initiate from earlier occurring RPWs, which, in turn, caused the development of new RPW events. This was proposed to be evidence that many of the RPW signatures that are seen at high spatial resolutions may be entirely chromospheric in origin. However, the authors also suggest that complex, twisted magnetic field geometry can create a scenario where wave emergence seems to contradict Madsen et al. (2015). As a result, with next generation instrumentation and facilities imminent, close attention will need to be paid to multi-wavelength (i.e., multi-height) observations in order to compare the small- and large-scale characteristics of RPWs, which will help to unequivocally determine the underlying physics that underpins their visible signatures.

In Fig. 48 we show a typical phase diagram from Felipe et al. (2010) that was obtained by simultaneously measuring the Doppler velocity at both photospheric and chromospheric heights in a sunspot umbra. Here, we see the clear effect of an embedded acoustic cut-off, with frequencies above $\sim 5.3$ mHz having a positive phase lag, highlighting the upward propagation of these waveforms. Consequently, these waves experience a rapid density drop as they propagate into the chromosphere, thus resulting in strong amplification of their amplitudes (see the lower panel of Fig. 48), which eventually results into shock formation. Following the methodology put forward by Ferraro and Plumpton (1958) and Centeno et al. (2006b), the amplitude, A, of a monochromatic wave with frequency, $\omega$, in a plane-parallel isothermal atmosphere permeated by a uniform vertical magnetic field comes from the solution to the equation,

$$\begin{aligned} c_{s}^{2} \frac{d^{2}A(z)}{dz^{2}} - \gamma g \frac{dA(z)}{dz} + \omega ^{2}A(z)=0 , \end{aligned}$$

(15)

where z is the vertical coordinate, g is the acceleration due to gravity, ${c_{s}=\gamma g H_{0}}$ is the speed of sound, $H_{0}$ is the pressure scale height, and $\gamma$ is the ratio of specific heats, which equals 5/3 for a monoatomic gas demonstrating adiabatic propagation. A solution to Eq. (15) is given by the trial function,

$$\begin{aligned} A(z)=e^{i k_{z} z} , \end{aligned}$$

(16)

where $k_{z}$ represents the vertical wavenumber. Solving for the vertical wavenumber, $k_{z}$, provides a dispersion relation of the form,

$$\begin{aligned} k_{z}=\frac{1}{c_{s}} \left( -i \omega _{c} \pm \sqrt{\omega ^{2} -\omega ^{2}_{c}} \right) , \end{aligned}$$

(17)

where $\omega _{c}=\gamma g/2c_{s}$ is the cut-off frequency. For $\omega <\omega _{c}$, $k_{z}$ takes imaginary values and the wave is evanescent. In the opposing regime (i.e., $\omega >\omega _{c}$), waves are able to propagate. This is illustrated in Fig. 49, where the phase angle is displayed as a function of frequency for waves measured at two independent geometric heights for both non-stratified and stratified atmospheric models. The cut-off frequency appears as a natural consequence of the stratification itself. However, as also shown in Fig. 49, a distinct cut-off frequency only exists in the limit of negligible radiative losses. In more realistic models, which include aspects of radiative cooling, a sharp separation between the propagating and evanescent regimes does not exist (solid line in Fig. 49), with the resulting phase diagram displaying a smooth transition at around 3 mHz. Of course, the above equations that represent the cut-off frequency are only valid in the limit of an isothermal atmosphere.

Fig. 48

Image reproduced with permission from Felipe et al. (2010), copyright by AAS

Phase spectrum (upper), coherence (middle; see Sect. 2.2.2), and amplification (lower) of Doppler velocity oscillations observed in the photospheric (Si i) and chromospheric (He i) spectral lines in a sunspot atmosphere. The red line in the upper panel represents the best fit from a theoretical model. The horizontal dashed black line at a coherence value of 0.7 in the middle panel highlights the lower confidence threshold. The red line in the lower panel represents the best fit from a theoretical model of acoustic waves propagating in an isothermal and stratified atmosphere (see text for more details).

Fig. 49

Image reproduced with permission from Centeno et al. (2006b), copyright by AAS

Phase difference spectrum for acoustic oscillations sampled at two geometric heights. The dashed and dot-dashed lines represent the expected phase as a function of frequency for non-stratified and vertically stratified isothermal atmospheres, respectively. The solid line indicates the phase relationship in a vertically stratified atmosphere including radiative losses from Newton’s cooling law. For each of the three cases depicted, the same plasma parameters are utilized ($T = 9000$ K, $\varDelta {z} = 1600$ km, and $g = 274$ m s$^{-2}$).

The strong vertical stratification of the atmospheric parameters in sunspots result in orders-of-magnitude changes to the propagation speeds of the embedded magnetoacoustic waves; namely the Alfvén speed, $v_{A}$, and the sound speed, $c_{s}$. This, together with the vertical and horizontal gradients of the background magnetic field, alongside other variations in local plasma parameters, constitutes important ingredients in the propagation characteristics of magnetoacoustic waves in these magnetic structures (MacBride et al. 2022). In Fig. 50 we show the variation of the sound and Alfvén speeds in a typical small sunspot model. To relate these variations to real observations we must employ approximations, which leads to slightly different but significant changes (Felipe and Sangeetha 2020). The cut-off frequency has been found to significantly change as a function of atmospheric height, generating important implications for both the heating of the upper layers of the Sun’s atmosphere and the wave propagation itself (Wiśniewska et al. 2016; Felipe et al. 2018a).

Fig. 50

Image reproduced with permission from Khomenko and Collados (2006), copyright by AAS

Contours of constant density (upper left) and constant magnetic field strength (lower left) for a typical small sunspot. In the lower-left panel the labels indicate the magnetic field strength in units of kG, while the thick black line denotes the isosurface where $v_{A} = c_{s}$. The dotted lines indicate the geometries of the embedded magnetic field lines. The red box corresponds to the domain size displayed in the right hand panels, where contours of constant Alfvén speed ($v_{A}$) and constant sound speed ($c_{s}$), in units of km/s, are depicted in the upper-right and lower-right panels, respectively. Note the strong horizontal gradients of both $v_{A}$ and $c_{s}$ due to the Wilson depression. The dotted lines, as per the lower-left panel, indicate the geometries of the magnetic field lines.

Further to changes with atmospheric height, it has been shown that the cut-off frequency depends on the magnetic field inclination, with more inclined fields allowing the upward propagation of frequencies below the $\sim 5.3$ mHz threshold (Bloomfield et al. 2007; Jess et al. 2013)—the so-called ramp effect. These waves are generally interpreted as longitudinal slow magnetoacoustic waves, with the general consensus on their origin being photospheric, through the absorption of externally driven p-modes (see for instance, Spruit and Bogdan 1992; Crouch and Cally 2003; Jess et al. 2012a; Krishna Prasad et al. 2015), with trans-sunspot oscillations at chromospheric heights potentially influencing the behaviour of these wave trains (e.g., Chae et al. 2017; Sych et al. 2020). Although these observational results are in agreement with the theoretical scenario of propagating slow magnetoacoustic modes, which in the low plasma-$\beta$ regime correspond to acoustic-like waves propagating along field lines with the magnetic pressure as the dominant restoring force, other magnetoacoustic modes exist in spatially uniform plasmas: namely an incompressible wave with magnetic tension as the sole restoring force (Alfvén wave), and an intermediate wave mode that can be thought of as a generalization of an acoustic wave with contributions from magnetic pressure (fast wave in the presence of a low plasma-$\beta$). In other words, in low plasma-$\beta$ environments the fast mode is an acoustic wave modified by the magnetic tension, capable of propagating isotropically with respect to the magnetic field.

Interestingly, at the equipartition layer where the sound and Alfvén speeds are nearly equal, a fraction of the energy, C, can be either channeled from a fast mode in the high plasma-$\beta$ regime (which is mainly an acoustic-like wave) to a fast magnetoacoustic mode in the low plasma-$\beta$ regime, or converted into a slow mode, thus preserving its acoustic nature. If we consider the sound speed, $c_{s}=\sqrt{\gamma P_{0}/\rho _{0}}$, and the Alfvén speed, $v_{A}=B/(4 \pi \rho _{0})$, where $\gamma$ is the adiabatic index, $P_{0}$ is the gas pressure, $\rho _{0}$ is the density, and B is the magnitude of the magnetic field strength, then the ratio between the two speeds squared can be given by,

$$\begin{aligned} \frac{c_{s}^2}{v_{A}^2}=\gamma \frac{4 \pi P_{0}}{B^{2}} . \end{aligned}$$

(18)

Substituting the magnetic pressure, $P_{B}=B^{2}/(8 \pi )$, into Eq. (18) we obtain,

$$\begin{aligned} {\frac{c_{s}^2}{v_{A}^2} = \frac{\gamma }{2} \beta . } \end{aligned}$$

(19)

This means that the equipartition layer is in practice close to the plasma-$\beta =1$ surface, and although they are conceptually different, they are often difficult to segregate from one another in observational data sequences (e.g., see the discussion points raised by Grant et al. 2018).

The wave translation process from one form to another is generally referred to as one of two processes: mode conversion or mode transmission (Crouch and Cally 2005; Suzuki and Inutsuka 2005; Cally and Khomenko 2015; Pagano and De Moortel 2017). Here, ‘mode conversion’ refers to a wave that retains its original character (i.e., fast-to-fast or slow-to-slow), yet converts its general nature in the form of acoustic-to-magnetic or magnetic-to-acoustic. Contrarily, ‘mode transmission’ refers to a wave that preserves its general nature (i.e., remains a ‘magnetic’ or ‘acoustic’ mode), yet changes character from fast-to-slow or vice versa. The fraction of energy that can be converted from fast to slow modes depends on the attack angle of the wave with respect to the magnetic field lines. The precise transmission coefficient, T, is defined as the proportion of incident wave energy flux transmitted from fast to slow acoustic waves (Cally 2001, 2007), which is governed by,

$$\begin{aligned} T = e^{-\pi k h_{s} sin^{2}(\alpha )} , \end{aligned}$$

(20)

where k is the wavenumber, $h_{s}$ the thickness of the conversion layer, and $\alpha$ the attack angle itself. The fast-to-fast conversion coefficient, C, can then be obtained by invoking energy conservation: $T + |C|=1$, where C is a complex energy fraction to take into account possible phase changes during the process of mode conversion (Hansen and Cally 2009). We note that the conversion coefficient, C, is larger when the frequency of the incident waves is higher and the attack angle is larger (Kontogiannis et al. 2014).

Schunker and Cally (2006) have shown how the combination of mode conversion alongside the ramp effect can result in an acoustic flux which is strongly dependent on the magnetic field geometry. This was also confirmed by Stangalini et al. (2011), who found a strong dependence of the wave flux between the photosphere and chromosphere on the magnetic field geometry inferred from spectropolarimetric inversions.

Grant et al. (2018) have also shown, by exploiting unique high-resolution observations and magnetic field extrapolations, combined with thermal inversions and MHD wave theory, that magnetoacoustic waves can couple with Alfvén waves at the equipartition layer in a sunspot, resulting in Alfvén-driven shocks that can efficiently contribute to the overall energy budget of the chromosphere (see Fig. 51).

Fig. 51

A cartoon representation of a sunspot umbral atmosphere demonstrating a variety of shock phenomena. A side-on perspective of a typical sunspot atmosphere, showing magnetic field lines (orange cylinders) anchored into the photospheric umbra (bottom of image) and expanding laterally as a function of atmospheric height, into the upper atmospheric regions of the transition region (TR) and corona. The light blue annuli highlight the lower and upper extents of the mode-conversion region for the atmospheric heights of interest. The mode-conversion region on the left-hand side shows a schematic of non-linear Alfvén waves resonantly amplifying magnetoacoustic waves, increasing the shock formation efficiency in this location. The mode-conversion region on the right-hand side demonstrates the coupling of upwardly propagating magnetoacoustic oscillations (the sinusoidal motions) into Alfvén waves (the elliptical structures), which subsequently develop tangential blue- and red-shifted plasma during the creation of Alfvén shocks. The central portion represents the traditional creation of umbral flashes that result from the steepening of magnetoacoustic waves as they traverse multiple density scale heights in the lower solar atmosphere. The image is not to scale and is reproduced from Grant et al. (2018)

Mode conversion and propagation of magnetoacoustic wave modes in sunspot atmospheres has also been investigated through numerical two-dimensional simulations (e.g., Khomenko and Collados 2008) incorporating realistic sunspot atmospheres (see Fig. 50). In particular, the atmospheric response to both longitudinal and transverse pulses (with respect to the magnetic field lines) has been investigated. Figure 52 reveals the velocity, magnetic field, and pressure fluctuations for a wave with input frequency above the cut-off value (10 s periodicity or 100 mHz) following 100 s of simulation run time for both longitudinal and transverse pulses, respectively. It is clear from Fig. 52 that the specific input pulse type results in different mixtures of transverse and longitudinal wave modes, which undergo mode conversion at the Alfvén/acoustic equipartition layer, and may also be reflected or refracted by the vertical and horizontal gradients.

Fig. 52

Image reproduced with permission from Khomenko and Collados (2006), copyright by AAS

Variations of the magnetic field (upper-left), pressure (upper-right), and velocity (transverse, lower-left; longitudinal, lower-right) at an elapsed time, ${t = 100}$ s, after the beginning of the simulations for a vertical longitudinal driver with a 10 s (100 mHz) periodicity. The lower 4 panels are identical to the upper 4 panels, only now show the variations for a horizontal transverse driver with a 10 s (100 mHz) periodicity. In each panel, the horizontal axis represents the radial distance from the center of the sunspot, while the inclined black lines highlight the magnetic field orientation. The two red lines in each panel indicate contours of constant $c_{s}/v_{A}$, with the thicker red line corresponding to $v_{A} = c_{s}$, and the thinner red line to $c_{s}/v_{A} = 0.1$. The thick black lines inclined towards the left (with increasing atmospheric height) of each panel indicate the direction of $\nabla v_{A}$, starting at the location of the pulse. Here, the $\nabla v_{A}$ line represents the boundary that separates waves refracting to the right from those refracting to the left, which is perpendicular to the contours of constant $v_{A}$ at every geometric height.

From theoretical studies, it was suggested that the enhanced 3-min wave power observed at chromospheric heights in sunspot umbrae may come from the presence of an acoustic resonance cavity, which is established by the temperature gradients at both the photospheric and transition region boundaries (see for instance, Hollweg 1979; Botha et al. 2011; Snow et al. 2015; Felipe 2019). Recently, Jess et al. (2020) exploited multi-height high spatial and temporal resolution observations, spectropolarimetric inversions, and numerical modeling to provide an observational confirmation of this physical mechanism. The authors examined the Fourier power spectra originating within a sunspot umbra and compared this to high-precision simulations encompassing a variety of different atmospheric stratifications. It was found that once steep temperature gradients were introduced into the simulation, the resulting cavity produced resonant amplification of the 3-min oscillations (see Fig. 53). Following on from the study by Jess et al. (2020) and Felipe et al. (2020) independently confirmed the presence of an acoustic resonator for another sunspot structure, and highlighted the potential importance of such findings for future helioseismic investigations of the solar atmosphere.

Fig. 53

A vertical stack of narrowband images taken across the Ca ii 8542 Å spectral line by the IBIS Fabry–Pérot instrument at the DST, revealing the photospheric (blue) to chromospheric (red) stratification of the sunspot atmosphere (left panel). This sunspot was found to display characteristics consistent with the presence of a resonance cavity, which caused the manipulation of spectral energies across the frequency domain (right panel), in particular providing a resonant enhancement at $\approx 20$ mHz. Image adapted from Jess et al. (2020)

3.3 Eigenmodes of large-scale magnetic structures

MHD theory applied to cylindrical magnetic flux tubes predicts a series of wave modes (Edwin and Roberts 1983), which represent the intrinsic overall response of the magnetic structure to the external forcing. These wave modes manifest over all scales of magnetic flux tube, however, due to their irregular cross sectional shapes and corrugated boundaries with surrounding plasma, sunspots and pores may significantly alter the characteristics of their resonant modes and their azimuthal properties (Albidah et al. 2021). Although as in any natural system one may also expect resonant modes in large scale magnetic structures, such as sunspots and pores, they are still poorly investigated in these structures. One reason is likely coupled to the fact that the amplitudes of the oscillations belonging to these modes are inherently weaker than that of the intensity and velocity excursions associated with local magnetoacoustic fluctuations such as umbral flashes, requiring particular filtering techniques for their identification.

In this regard, a first attempt to disentangle the signal associated with resonant modes from the rest of the spatially incoherent oscillations was made by Jess et al. (2017), who applied three-dimensional Fourier (i.e., $k{-}\omega$; see Sect. 2.3.1) filtering to study the intensity oscillations of the sunspot (shown in Fig. 54) at chromospheric heights. The $k{-}\omega$ diagram of the intensity oscillations displays several horizontal ridges, which are evident in the right panel of Fig. 18. These spectral features are associated to the spatial scale of the entire sunspot (i.e., at wavelengths corresponding to the size of the overall sunspot umbra), highlighting the presence of a spatially coherent oscillations that encompass the entire sunspot umbra itself. It is worth stressing here that, in order for these modes to be readily identified, one requires sufficient spectral resolution in both the temporal and spatial frequency domains. As discussed in Sect. 2.2.1, to maximize the temporal frequency resolution requires the acquisition of long duration data sequences. In an analogous manner, large field-of-view sizes are required to provide sufficient frequency resolution in the spatial domain (i.e., a small wavenumber resolution, $\varDelta {k}$).

Fig. 54

Image reproduced with permission from Jess et al. (2017), copyright by AAS

Sample images of a near circularly symmetric sunspot, including the vertical component of the magnetic field ($B_{z}$; upper left), the photospheric continuum (upper middle), and the chromospheric H$\alpha$ line core (upper right). The color bar corresponding to the strength of the magnetic field is saturated at $\pm 1000$ G for visual clarity. The lower left panel displays a snapshot of H$\alpha$ intensities following $k{-}\omega$ filtering (i.e., in both temporal and spatial domains; see Sect. 2.3.1). The solid green contour outlines the time-averaged umbra/penumbra boundary, while the red annulus depicts the extent of the region used for examining azimuthal wave motion within the umbra, where the center of the annulus is placed at the center of the umbra. The lower right panel is a time-azimuth diagram following the polar transformation of the signals contained within the red annulus in the lower left panel, which allows the circular nature of the wave rotation to be investigated in a similar way to traditional time-distance diagrams. The horizontal dashed green line highlights the azimuthal intensity signal corresponding to the filtered image shown in the lower left panel, while the solid red line represents the fitted angular frequency (i.e., degrees per second) of the rotating wave amplitudes.

After filtering at the wavelengths and frequencies corresponding to the horizontal spectral ridges shown in the $k{-}\omega$ diagram in the right panel of Fig. 18, Jess et al. (2017) were able to detect a coherent rotational wave within the sunspot umbra. The observed rotational motion is visible in the lower panels of Fig. 54, where the lower-left panel shows the reconstructed intensities following $k{-}\omega$ filtering within the solid black box depicted in the right panel of Fig. 18. The lower-right panel of Fig. 54 reveals the time-azimuth map following the polar transformation of the circular sunspot wave patterns, where the straight diagonal trends highlight coherent bulk rotations of the MHD wave phenomena. Thanks to numerical MHD modeling, it was possible to interpret these results as the first detection of an $m = 1$ slow magnetoacoustic mode in the chromospheric umbra of a sunspot. Building upon the early work of Ulrich (1970) and Deubner (1975), which has subsequently been improved through the application of more sensitive instrumentation and modern techniques (e.g., Kosovichev et al. 1997; Rhodes et al. 1997; Haber et al. 1999; Christensen-Dalsgaard 2002; Howe et al. 2004; González Hernández et al. 2006, to name but a few), significant Fourier power at $\sim 5$ mHz (i.e., consistent with the generalized p-mode spectrum) has demonstrated coherency down to wavenumbers on the order of $k \sim 0.05$ arcsec$^{-1}$, corresponding to spatial wavelengths of approximately $140''$ ($\sim 100{,}000$ km), with radio observations from missions such as GOLF and BiSON observing global wavenumbers (e.g., Jiménez-Reyes et al. 2004; Chaplin et al. 2007). This may suggest that the elevated Fourier power bands shown in the right panel of Fig. 18 may be linked to large-scale, sub-surface drivers. While the HARDcam H$\alpha$ observations are chromospheric in nature, being formed at a geometric height of $\sim 1500$ km (Vernazza et al. 1981; Leenaarts et al. 2012), the highly magnetic composition of sunspot structures may enable direct and efficient coupling with the solar layers below (see, e.g., the recent review by Cally et al. 2016).

In addition, sausage modes (see Dorotovič et al. 2008; Morton et al. 2011; Dorotovič et al. 2014; Grant et al. 2015; Freij et al. 2016; Grant et al. 2022) have also been identified in magnetic pores through variations in their cross-sectional area of the magnetic structures and associated out-of-phase intensity oscillations. This highlights the interesting possibility to simultaneously exploit several diagnostics for the identification of global resonances in magnetic structures. Besides velocity and intensity perturbations, MHD waves are also expected to be characterized by magnetic perturbations. However, their identification has been long debated, as opacity effects can also mimic magnetic perturbations (for further details, see Khomenko and Collados 2015, and the references therein). However, the recent advances in multi-height spectropolarimetric imaging observations have enabled the use of phase lag analyses between different layers of the solar atmosphere, which can be used to robustly identify real magnetic oscillations and disentangle them from spurious effects (e.g., changes in opacity). By doing this, Stangalini et al. (2018) were able to identify propagating magnetic fluctuations at the umbra-penumbra boundaries of a large sunspot observed by IBIS, which constitutes a spatially coherent oscillation that was interpreted as the signature of a surface mode of the sunspot flux tube. By studying the phase relationship between circular polarization (CP) and intensity signals, it was also argued that the oscillations were not consistent with opacity effects (see Fig. 55).

Fig. 55

Left: Phase lag map of the amplitude of circular polarization (CP) fluctuations at 3 mHz (across a bandwidth of 0.7 mHz) between the photosphere and chromosphere. The solid white contour depicts the umbra/penumbra boundary. The positive phase lags towards the edges of the sunspot umbra reveal the presence of upwardly propagating magnetic waves. Right: Phase lag diagram obtained in umbra-penumbra boundary region showing the presence of a distinct positive peak (downward propagation) in the spectrum. Plots reproduced from Stangalini et al. (2018)

Studying sausage modes in numerous pore data sets obtained with ROSA, Keys et al. (2018) were able to identify signatures of surface and body waves associated with the oscillating pores. Surface and body modes can be identified by the spatial distribution of the amplitude across the flux tube (Rae and Roberts 1983; Erdélyi and Fedun 2010). For surface modes, all perturbations have a maximum at the tube boundary and will be zero at the center of the tube. For body modes, the kinetic gas pressure, $v_z$ and $B_z$ will have maximum amplitude at the center of the tube, decreasing to the tube boundary. In the case of body modes, higher harmonics in the radial direction may result in nodes between the axis of symmetry and the boundary of the flux tube. For both surface and body modes in a homogeneous ambient plasma, the external wave power should decrease exponentially as a function of distance from the flux tube boundary. Figure 56 shows a schematic diagram (adapted from Keys et al. 2018) of the expected spatial distribution of power for both body (left image in panel a) and surface (right image of panel a) waves in a cylindrical flux tube. To detect these signatures, the authors identified oscillations in the pore data sets by looking for oscillations in area and intensity in the pores, which would indicate the presence of a sausage mode. Periodicities were found to range from 90–700 s, with the most common periods in area and intensity occurring at $\sim 300 \pm 45$ s. To determine if the wave was a surface or body mode, the spatial distribution of the power was then analyzed. This was performed by employing Gaussian filtering of the data for the dominant oscillatory frequencies within the data. This was limited to sections of the time series where significant power was found from the wavelet and EMD analysis of the area and intensity signals for the pores.

Fig. 56

Panel a shows a schematic of the spatial structure of the pressure perturbation for the body mode (left) and the surface mode (right) in both two and three dimensions. Arrows at the bottom of each schematic show the sausage oscillation in the flux tube. Panel b shows an example of a body mode in an elliptical pore. The lower left image shows a G-band intensity image of the pore, while the image above shows the two-dimensional power plot for the pore filtered at a central frequency of 11.1 mHz. The white contour here indicates the location of the pore. The blue cross-cut on the intensity image indicates the region taken for the one-dimensional power plot across the pore shown on the right. Panel c shows a similar example for a surface mode. Again, the blue cross-cut on the intensity image shows the location for the corresponding one-dimensional power plot shown to the right in this panel. In this instance, the power plot is produced for the pore filtered at a central frequency of 2.2 mHz, consistent with the timescale of granular evolution. Red dashed lines in the one-dimensional power plots in both panels b and c indicate the pore boundary in both instances. The body mode is characterized by a central peak decaying to the pore boundary, while the surface mode is characterized by peaks in power at the pore boundary decaying to zero in the center of the pore. Image adapted from Keys et al. (2018)

By looking at the spatial distribution of the power, Keys et al. (2018) were able to identify surface and body modes associated with their data sets. The right panels of Fig. 56 shows examples of the power distribution observed by Keys et al. (2018) for both the body (upper plot) and surface (lower plot) waves. The surface mode was determined to be the most frequently occurring of the two. The authors suggest that this could be due to either the size or the magnetic field strength of the pore, as smaller weaker pores were more likely to display signatures of body modes. The authors suggest that it is possible that with a stronger field strength, there is a larger magnetic field gradient between the pore and the ambient plasma, possibly resulting in surface modes when the field strength is larger. The authors are clear to note though that the number of samples they have is small and, so, this relation is something that needs to be analyzed further. Estimates were also made for the energies associated with the observed surface and body modes using the framework of Moreels et al. (2015). Surface modes had an energy flux estimate of $22\pm 10$ kW m$^{-2}$, while the body modes had an observed energy flux of $11\pm 5$ kW m$^{-2}$.

Recently, Stangalini et al. (2022) and Albidah et al. (2022) have found multiple slow body modes in sunspot umbrae, from low to high order in both the photosphere and chromosphere (see Fig. 57, 58, 60). These works have shown that higher order modes (which have smaller spatial regions of phase and anti-phase) are particularly sensitive to the cross-sectional shape of the umbra. It can be seen for the $m=2$ slow body mode shown in Fig. 57 detected in an approximately elliptical sunspot by Albidah et al. (2022) that the Pearson correlation between the POD/DMD modes and the the modes predicted by the elliptical and exact shape models are very similar (the regions of red show a strong in phase correlation). However, the first kink overtone, which has smaller regions of phase and anti-phase compared with the $m=2$ fluting mode, shown in Fig. 58, demonstrates that model with the exact umbral shape has a much stronger correlation to the observed POD/DMD modes. This shows that higher order modes “feel” the irregularities in umbral cross-sectional shapes more than the lower order modes (see Albidah et al. 2022, for more examples).

Fig. 57

Image reproduced with permission from Albidah et al. (2022), copyright by the authors

The first row displays the spatial structure of the modes that were detected from the observational HARDcam data in an approximately elliptical sunspot by Albidah et al. (2022): the first POD mode (middle) and the DMD mode that corresponds to the frequency of 5.6 mHz. In the first column, the theoretical spatial structure of the fundamental slow body fluting mode ($m=2$), even with respect to the major axis, in the elliptical magnetic flux tube model (middle) and the corresponding model using the exact umbral shape (bottom). The rest of the panels show the Pearson correlation between the theoretical and POD/DMD modes. The positive/negative numbers in the color bars denote regions of phase/anti-phase. The dashed line show the boundary of the theoretical elliptical tube, and the solid black line shows the actual umbra/penumbra boundary.

Fig. 58

Image reproduced with permission from Albidah et al. (2022), copyright by the authors

This is the equivalent display for the same approximately elliptical sunspot as shown in Fig. 57 for the first overtone of the slow body kink mode which is odd with respect to the major axis. Here the DMD mode frequency is 5.3 mHz.

Furthermore, if the observed sunspot umbra is far from circular or elliptical such simple models cannot be applied and the actual cross-sectional shape must be used, even to get an accurate representation of low order modes. Such an example from Stangalini et al. (2022) is shown in Fig. 59 where first fifty eigenfunctions are modelled for the vertical velocity component of slow body modes in a large sunspot umbra (40 Mm across). The nine most dominant modes are highlighted in Fig. 59 and were used to reconstruct the observed Doppler signal to high accuracy as shown in Fig. 60. This would not have been possible by simply approximating the observed irregular umbra cross-sectional shape with circular or even elliptical cylinder MHD wave models.

Fig. 59

Image reproduced with permission from Stangalini et al. (2022), copyright by the authors

The first fifty slow body mode eigenfunctions calculated for the vertical component of velocity using the observed sunspot umbra shape by Stangalini et al. (2022). The dominant nine mode numbers are highlighted in red and these were used to reconstruct the observed Doppler velocity signal to high accuracy as shown in Fig. 60.

Fig. 60

Image reproduced with permission from Stangalini et al. (2022), copyright by the authors

Comparison between the spatially coherent wave pattern observed in a sunspot after properly filtering the data (left), and the one expected from numerical modeling (right, see Fig. 59 for the eigenfunctions that were used). The filtered wave pattern displays a high order oscillation in the umbra that agrees very well with numerical predictions of eigenmodes. B–$\omega$ diagram of the same sunspot showing a rich variety of frequencies excited within the umbra and consistent with resonant modes of the magnetic structure.

From theory, each individual MHD wave mode can be broadband in both $\omega$ and k, as is shown for the magnetic cylinder model dispersion diagram in the left panel of Fig. 38, i.e., each MHD mode on the dispersion diagram forms a continuous line, meaning that for every $\omega$ there is a unique k that will satisfy the phase speed relation for a particular mode. This broadband behaviour has now actually been observed, as can be seen in Fig. 60 taken from Stangalini et al. (2022), where the $B{-}\omega$ diagram of a sunspot in the photosphere shows that, within the umbra, there are multiple strong frequency peaks, other than the usual dominant p-mode frequency, excited at the same time. It must be emphasised that the actual fine structure of the frequency power spectrum inside the umbra was not predicted at all from any model and therefore potentially opens up a whole new field of lower solar atmospheric MHD waves research.

The limitation of studies by Keys et al. (2018), Albidah et al. (2022) and Stangalini et al. (2022) is that they only looked at either the photospheric or chromospheric signatures of MHD waves in pores and sunspots. Multi-height studies of MHD waves in pores and sunspots are needed to determine the variation of properties and energy flux with height to give a clearer understanding of their contributions and relative importance to heating.

It also must be noted that the identification of global eigenmodes in large-scale magnetic wave guides is synonymous with the identification of coherent wave motions in a particular flux tube, be it a pore or sunspot umbra. Nevertheless, the limited number of works in this direction (Yuan 2015; Jess et al. 2017; Keys et al. 2018; Albidah et al. 2021, 2022; Stangalini et al. 2022) testifies to the intrinsic difficulty of this task. As mentioned above, a potential reason for that could be due to the intrinsically small amplitudes of the associated oscillations, compared to the other omnipresent spatially incoherent fluctuations. However, filtering techniques represent a viable solution, but these require very high spatial and temporal resolution data with enough temporal and spatial coverage to reach the necessary frequency and wavelength resolutions in $k{-}\omega$ space.

3.4 Small-scale magnetic structures

Concentrations of intense magnetic fields at small scales are mostly found in intergranular regions (in the photosphere), where strong downflows occur (Stenflo 1973; Jess et al. 2010a; Riethmüller et al. 2014; Borrero et al. 2017). In such small magnetic elements, in addition to gravity and pressure, the magnetic field also acts as a restoring force (e.g., Steiner 2010). Therefore, in addition to the longitudinal (compressible) acoustic and gravity waves, other MHD wave modes may also propagate along the flux tubes (Roberts and Webb 1978; Edwin and Roberts 1983; Hasan and Sobouti 1987; Fleck et al. 1993; Steiner et al. 1998; Bogdan et al. 2003; Musielak and Ulmschneider 2003a; Khomenko et al. 2008b; Fedun et al. 2011a). Such waves can be generated as a result of, e.g., (1) buffeting of the flux tubes (i.e., bundles of magnetic-field lines; Spruit 1976, 1981a; Solanki et al. 1996) by the surrounding granules (and intergranular turbulence), i.e., kink modes (Musielak and Ulmschneider 2003b), (2) compression and contraction of the flux tubes by convective forces from opposite directions (sausage modes), or (3) twisting the flux tubes by rotating flows around the tubes (torsional Alfvén waves; Spruit 1982; Solanki 1993). It is, however, more complex in real situations where wave modes of various eigenmodes may co-exist in the same magnetic elements. Excitation of various wave modes in a magnetic cylinder has been reviewed in detail by Nakariakov and Verwichte (2005). Such waves may occur in either propagating or standing states (Rosenthal et al. 2002; Dorotovič et al. 2014). The propagating magnetoacoustic (or MHD) waves are channeled from the base of the photosphere to the chromosphere and beyond along the magnetic field lines (where the magnetic field acts as a guide; Khomenko et al. 2008a) whose strength and inclination plays a role in their leakage to the upper solar atmosphere (Michalitsanos 1973; Jefferies et al. 2006; Schunker and Cally 2006; Stangalini et al. 2011). In addition, the magnetoacoustic waves may propagate as ‘fast’ or ‘slow’, traveling at speeds faster or slower compared to the ratio of the Alfvén and sound speeds. The various wave modes may be converted from one form to another (i.e., mode conversion) and/or only switch the fast and slow labels (i.e. mode transmission) at the plasma-$\beta \approx 1$ level, where the Alfvén and sound speeds nearly coincide (Khomenko and Cally 2019; see Sects. 2.9 and 3.2.1 for greater details).

As a (thin) flux tube extends into the solar atmosphere, it expands with height (while the gas density and pressure decrease). In addition, many of such flux tubes bend over at lower atmospheric heights (compared to larger and stronger magnetic-field structures, thus producing the multi-height magnetic canopy; Giovanelli and Jones 1982; Solanki et al. 1991; Rosenthal et al. 2002; Jafarzadeh et al. 2017a). The (inclined) flux tubes are thought to be observed throughout the solar chromosphere as (dark or bright) thread-like structures in intensity images (De Pontieu et al. 2007b; Pietarila et al. 2009; Rouppe van der Voort et al. 2009; Pereira et al. 2012; Gafeira et al. 2017b; Jafarzadeh et al. 2017a; Kianfar et al. 2020). Only recently, has it been possible to identify the MHD wave modes in small-scale structures through the entire lower solar atmosphere, thanks to high-resolution observations provided by modern facilities. Such small structures are often very dynamic and short lived (with observational timescales on the order of a few seconds to a few minutes), therefore, not only are high-spatial resolution observations required to resolve them, but any study of their rapid evolution also needs a high temporal resolution. Furthermore, multi-line (i.e., multi height) observations are also essential in order to trace waves as they propagate through the atmosphere. In this regards, narrow-band observations in spectral lines (with a relatively high spectral resolution) would reduce mixing information from different atmospheric heights, but ambiguities may still exist.

In the following, we summarize the most recent advances on detection and analysis of various MHD wave modes in small-scale magnetic elements and fibrillar structures in the solar photosphere and the chromosphere.

3.4.1 Excitation, propagation, and dissipation of MHD waves in small-scale magnetic structures

Due to different (distinct) kinematics of various solar photospheric regions (with different levels of magnetic flux; Abramenko et al. 2011; Stangalini 2014; Keys et al. 2014; Jafarzadeh et al. 2017b), characteristics of waves and oscillations in small-scale magnetic elements may depend on the environment in which they are embedded. Small magnetic elements are able to laterally move within a supergranular body (i.e., the internetwork). In concert with surrounding granules interacting, due to expansions and explosions, and effects of intergranular turbulence, a variety of MHD wave modes can be generated across a range of frequencies. Of particular interest, from the detection point of view, is also a lower number density of these elements in the internetwork, hence, they are more isolated compared to those found in network/plage regions. In the latter, while such interactions also occur between the plasma and magnetic elements, they are often found in groups of concentrated field structures trapped in sinks (stagnation points) where inflows from surrounding supergranules prevent them from moving in a preferred direction, but rather they experience random walks within a relatively small area (van Ballegooijen et al. 1998; Nisenson et al. 2003; Utz et al. 2010; Manso Sainz et al. 2011; Chitta et al. 2012b; Jafarzadeh et al. 2014a; Giannattasio et al. 2014).

Oscillations (of different properties) in the low photosphere, in both network and internetwork magnetic elements as well as those in plage areas (in the vicinity of large magnetic structures), have been identified over the past decades (e.g., Noyes and Leighton 1963b; Simon and Leighton 1964; Howard 1967; Ulrich 1970; Giovanelli 1972; Canfield and Mehltretter 1973; Stein and Leibacher 1974; Goldreich and Keeley 1977; Christensen-Dalsgaard and Gough 1982; November and Simon 1988; Toutain and Froehlich 1992; Fontenla et al. 1993; Zaqarashvili 1999; Gizon et al. 2003; Vecchio 2006, to name but a few). Thanks to simultaneous multi-height (multi spectral-line) observations of the entire solar photosphere and chromosphere (although at different resolution/band width), propagation of the various types of magnetoacoustic waves have also been investigated. However, it is worth noting that the highest chromospheric layer to which such waves have been traced can depend on the properties of the lines in which the observations were made and/or the magnetic topology of the observed area. Thus, the relatively wide-band observations of the chromospheric lines include a wide range of chromospheric heights. In contrast, narrow-band observations have revealed more filamentary structures of the chromosphere, while their number density and thickness tend to increase with height through the entire solar chromosphere. Therefore, formation heights at which waves are identified/traced should be interpreted with caution. In addition, the presence or density number of fibrillar structures at chromospheric heights may depend on the level of magnetic flux within (and/or in the immediate vicinity of) the observed photospheric field of view. Using the Ca ii H filter (with a width of 0.1 nm) onboard the Sunrise balloon-borne solar observatory, Jafarzadeh et al. (2017a) illustrated a field of view of an active region filled with slender fibrils where almost no magnetic bright points could be observed. On the contrary, they presented a quiet-Sun region (taken with the same filter) where no fibrillar structures appeared.

Various kinds of MHD modes have been identified at small scales in high-resolution observations. Often they are identified in intensity images of features such as magnetic bright points (MBPs) or in fibrillar structures at chromospheric heights, utilizing various spectral lines, thus, sampling different atmospheric layers. As a result, their propagation from the solar photosphere to the chromosphere, or in a particular region within these layers, has been characterized in a number of studies (e.g., Lites and Chipman 1979; Kalkofen 1997; McAteer et al. 2002a; De Pontieu et al. 2005; Okamoto and De Pontieu 2011; Jess et al. 2012c; Kuridze et al. 2012; Jafarzadeh et al. 2017d; Bate et al. 2022, to name a few). Such waves include transverse oscillations (often interpreted as MHD kink or Alfvénic waves), oscillations in intensity and/or Doppler velocity (characterizing longitudinal magnetoacoustic waves), twist perturbations (describing torsional Alfvén waves), as well as fluctuations in the size/width of small magnetic structures, as a signature of (compressible) MHD sausage modes, particularly, when they are in anti-phase relationships with intensity oscillations (Edwin and Roberts 1983; Centeno et al. 2006a; Erdélyi and Morton 2009; Mathioudakis et al. 2013; Jess and Verth 2016). Kink modes and Alfvén waves are incompressible and their dissipation in the solar atmosphere requires a large gradient in the Alfvén speed (Goossens et al. 2009). Alfvén waves also result in Doppler velocity perturbations, while the longitudinal magnetoacoustic waves result in fluctuations in both intensity and Doppler velocity. Therefore, spectral observations with sufficiently high wavelength resolution for fine doppler studies should be interpreted by taking the nature of various wave modes, as well as mode-coupling/mixing, into consideration.

It is thought that rapid (greater than $\approx 2$ km/s) pulse-like kicks to small magnetic elements, as a result of, e.g., granular explosions, can excite transverse kink waves along the flux tubes (Spruit 1981b; Choudhuri et al. 1993a, b; Steiner et al. 1998; Hasan and Kalkofen 1999; Musielak and Ulmschneider 2003b; Möstl et al. 2006). Such impulsively excited waves, as a result of rapid continuous jostling of the flux tube by granules (Hasan et al. 2000) can upwardly propagate into the upper solar chromosphere, with their amplitudes increasing exponentially. The kink waves may become nonlinear in the upper chromosphere where their propagating speeds are comparable to tube speeds (Kalkofen 1997). Such waves may, however, couple to the longitudinal magnetoacoustic waves in the low-to-mid chromosphere and dissipate by forming shocks (Ulmschneider et al. 1991; Zhugzhda et al. 1995). Muller et al. (1994) examined the pulse-like excitation mechanism and estimated an energy flux of 2000 W/m$^2$ that could be carried by kink waves in network MBPs. A smaller energy flux on the order of 440 W/m$^2$ was later reported by Wellstein et al. (1998), based on horizontal motion of chromospheric Ca ii K bright points in network areas. However, the magnetic nature of such small-scale brightenings were not determined.

Thanks to the high spatial resolution provided by Sunrise, Jafarzadeh et al. (2013) were able to identify such jerky rapid (sometimes supersonic) pulse-like motions in small internetwork MBPs observed in the upper photosphere/low chromosphere. Such waves were found to be energetic enough (with a net energy flux of $\approx 300$ W/m$^2$) to potentially heat the outer solar atmosphere. The somewhat large difference between the energy fluxes found by different authors could be due to, e.g., different geometric heights, network versus internetwork (the latter hosts considerably smaller number of MBPs, that can move around more freely, compared to the former), spatial resolution (influencing the number of detected elements, their sizes, and horizontal-motion measurements), and their nature (magnetic versus non-magnetic features).

Incompressible horizontal convective motions, on the solar surface, are generally thought to be the prime excitation mechanism of the transverse MHD waves in magnetic flux tubes, by being either perpendicular or tangential to the surface of magnetic elements (resulting in kink or torsional Alfvén modes, respectively; e.g., Narain and Ulmschneider 1996; Erdélyi and Fedun 2007). In addition, turbulent convective downflows have been suggested to generate transverse displacements inside the magnetic concentrations (van Ballegooijen et al. 2011), occurring on smaller length and time scales, compared to those from the granular motions.

Morton et al. (2014) exploited observations from multiple instruments (i.e., from the Swedish Solar Telescope (SST; Scharmer et al. 2003), Hinode/SOT (Kosugi et al. 2007; Tsuneta et al. 2008), DST/ROSA, and Coronal Multi-channel Polarimeter, CoMP) to study the generation and transport of energy by kink waves in small scale structures through the entire quiescent solar atmosphere. They found similar power spectra for transverse oscillations of photospheric MBPs (and granular flows) and the chromospheric H$\alpha$ fibrils, suggesting the granular motions have excited the kink waves identified in the small structures. In addition, Morton et al. (2014) found that the higher-frequency wave energy was significantly diminished in the corona’s power spectra, thought to be a signature of energy dissipation at those frequencies. However, the authors give no consideration to the attenuation effects of the contribution function and transmission of spectral lines in the solar atmosphere at high frequencies. Deubner (1976) found notable reductions in wave power at evenly spaced frequencies above 10 mHz, and proposed that the spacing equated to wavelengths that were integers of the length of the atmospheric region contributing to the spectral line. Subsequently, it was conclusively shown that the extent of the atmospheric column that contributes to a spectral line directly impacts the transmission of high frequency waves in the atmosphere (Durrant 1979; Cram et al. 1979; Mein and Mein 1980). As a result, the transmission of frequencies is dependent on both the wavelength, and the spectral window that observations are integrated over (Fossum and Carlsson 2005b). In relation to energy transport and dissipation, differential transmission has been shown to impact the potential of acoustic waves in particular to propagate energy in higher frequency modes beyond the photosphere (e.g., Bello González et al. 2009, 2010b). In the case of Morton et al. (2014), a variety of different spectral lines are studied, including $H\alpha$ imaging with a spectral window of 0.25 Å and sampling a range of formation heights. Without consideration of the relative transmission functions of the multi-wavelength observations, it cannot be verified whether the damping observed in Morton et al. (2014) is dissipative, or due to transmission effects. This must also be taken into account whenever damping is observed across different spectral line observations, either in mitigation, or through transmission function analysis.

Nonlinear propagation of transverse waves to the solar chromosphere at small scales, previously predicted by theoretical models (Moreno-Insertis 1986; Ulmschneider et al. 1991), was identified by Stangalini et al. (2015), where the authors exploited transverse perturbations in several MBPs simultaneously recorded in both the photosphere and the low chromosphere at high resolution with SST. They found the identified kink waves to nonlinearly propagate upward above a cut-off frequency of $\approx 2.6$ mHz. The nonlinearity was concluded due to remarkable differences between the photospheric and chromospheric power spectra (i.e., considerably different patterns of peaks in the power spectra).

Using a relatively long time series (of $\approx 4$ h) from Hinode/NFI (Tsuneta et al. 2008), Stangalini et al. (2013a) provided the full power spectra of transverse (kink) oscillations in small photospheric magnetic elements (limited to the spatial resolution of the 0.5 m telescope). They found a wide range of frequencies of 1–12 mHz, of which, the lower frequencies would only be reliably identifiable by exploiting such a long image sequence, which is rare for ground-based observations. However, on the higher frequency end, the measurements were limited to the 30 s cadence of the observations, hence, a Nyquist frequency of 16.7 mHz. Therefore, detection of higher frequencies would only be possible with a higher temporal resolution. In addition, it is worth noting that spatial resolution is also an important factor for detecting power at high frequencies (Wedemeyer-Böhm et al. 2007). Therefore, the higher spatial and temporal resolutions are, the higher frequency oscillations can be detected (if they exist).

Following earlier works by Hasan et al. (2003, 2005) and Hasan and van Ballegooijen (2008) employed numerical simulations to propose that the excess brightness in the network Ca ii H & K MBPs (in the solar chromosphere) could be due to (1) high-frequency (higher than 10 mHz) transverse oscillations at the base of the magnetic flux tubes, and (2) absorption of acoustic waves from the surrounding medium. These result in temperature perturbations (of up to 900 K) due to shock dissipation at chromospheric heights, following the upwardly propagating (slow) magnetoacoustic waves along the flux tubes. Such propagating high-frequency transverse waves (up to 23 mHz) as well as (longitudinal) intensity oscillation (up to 30 mHz) in small MBPs were detected by Jafarzadeh et al. (2017d) from high temporal and spatial resolution observations with Sunrise. These authors studied the MBPs at two atmospheric heights, corresponding to the low photosphere and the low chromosphere, with an approximate height difference of 450 km on average (estimated using two independent approaches). Together with phase differences between the intensity oscillations at the two atmospheric layers, a wide range of propagating velocities were determined. Of which, phase speeds larger than 30 km/s could not satisfy expected propagating speeds (from theoretical models) at these heights. Uncertainties in $I-I$ phase analysis can be introduced through radiative damping, particularly in atmospheric regions where the radiative relaxation time is equivalent to, or less than, the wave period (e.g., Souffrin 1972; Schmieder 1978; Deubner et al. 1990). Estimates predict that 3–5 min oscillations are on the order of the relaxation time in the low chromosphere (Mihalas and Toomre 1982; Severino et al. 2013), thus the unexpectedly large phase speeds will be influenced by non-adiabatic atmospheric evolution. In addition, refraction of the propagating path of these waves may influence phase speed estimations. Using numerical simulations, Nutto et al. (2012b) studied such possible modifications of wave propagation (and wave travel time) in small magnetic features in the solar network atmosphere. They found that the travel time (and hence the propagating speed) is strongly influenced by mode conversion, sometimes at multiple plasma-$\beta =1$ levels which are placed on top of each other as a result of the highly dynamic atmosphere. In addition, they found that the measured wave travel-time could significantly be reduced as a result of the fast waves being refracted above the magnetic canopy due to the large gradient of the Alfvén speed. Thus, the two mechanisms, i.e., the fast waves due to (multiple) mode conversion inside the magnetic canopy and the refraction of the propagation path above the canopy may lead to observations of wave travel time that is too short (i.e., propagating speeds which are too large) between two atmospheric layers. Short time delays of 4 s and 29 s were also reported by Keys et al. (2013) between horizontal velocity variations within 500 km from simulated data and within heights sampled by the G-band and Ca ii K MBPs from DST/ROSA, respectively. The authors interpreted such short time intervals as the results of oblique granular shock waves in the simulations, and of a semi-rigid flux tube in the observations.

Incompressible kink and compressible sausage modes at high frequencies (of $\approx 12$ and 29 mHz, respectively) were also detected in slender Ca ii H fibrils (located in the low-to-mid chromosphere) from high-resolution observations with Sunrise (Jafarzadeh et al. 2017c; Gafeira et al. 2017a). Figure 61 shows such slender fibrillar structures, filling the entire field of view (left), with an example of the detected transverse oscillations at 5 locations along one fibril (right). The fibril and locations of the artificial slits (marked with a-e) are illustrated on the top of the right panel. The transverse oscillations are identified in space-time plots, where time variations of the location of fibrils has been inspected at each ‘cut’ perpendicular to the fibril’s axis. Slope of the lines (green), connecting the same peaks/troughs of the oscillations at different locations, indicate the propagation of the transverse (kink) waves from right to left in the top panel. To quantify the propagating speeds (and periods), Jafarzadeh et al. (2017c) employed a wavelet analysis to compute phase differences between oscillations at different locations (whose distances are known). The energy flux transported by the kink waves along these slender fibrils was found to be $\approx 15$ kW/m$^2$, on average.

Fig. 61

Images reproduced with permission from Jafarzadeh et al. (2017a, c), copyright by AAS

Small-scale (slender) Ca ii H fibrils in the low chromosphere (left). Transverse oscillations in one fibril are illustrated on the right panel in space-time plots, at multiple cuts in different locations along the fibril (shown on the top). The solid (green) lines connect the extrema of the fluctuations, indicating wave propagation from right to left along the fibril.

Furthermore, Gafeira et al. (2017a) identified sausage modes (with periods on the order of 32–35 s) by measuring intensity and size at various locations (cuts) along these small-scale fibrils. The left panel in Fig. 62 illustrates one example where the measurements of both intensity and width (by fitting a Gaussian function perpendicular to the fibril’s axis) at multiple locations are shown on the stack of one fibril (on top of each other) at different times. The vertical lines, depicted at the different spatial locations, indicate the width of the fibril at those locations (marked with a small circle). The fluctuations of intensity and width at one location, marked with the arrow on the left panel, is presented on the right, where a clear anti-correlation between the two oscillations is evident. The authors also measured the wave properties by means of wavelet analysis, resulting in a propagating speed of 11–15 km/s, interpreted as fast sausage modes. These phase speeds are considerably smaller than those found in chromospheric H$\alpha$ fibrils, propagating with $\approx 67$ km/s, on average, for the fast sausage modes (Morton et al. 2012), also with longer periods on the order of 197 s. Morton et al. (2012) also detected transverse oscillation in the same H$\alpha$ fibrils with periods and propagating speeds of 232 s and $\approx 80$ km/s, respectively. The increase in propagating speed with height is expected from theoretical models, due to the height stratification of the physical parameters or, possibly, these waves are close to the cut-off frequency and are close to being evanescent resulting in the observed high speeds.

Fig. 62

Images reproduced with permission from Gafeira et al. (2017a), copyright by AAS

Measurements of intensity and width in different locations along one slender Ca ii H fibril is illustrated on the stack of the fibril at different times (left). Oscillations of the two quantity at one location (marked with the arrow) is illustrated on the right panel.

Generation and propagation of both kink and sausage modes in chromospheric fibrillar structures (i.e., on-disc Type i spicules) were studied in detail by Jess et al. (2012b) from both observations (from DST/ROSA) and MHD simulations, where the mode conversion at the lower solar atmosphere was found to be the main driver of the MHD waves. Jess et al. (2012b) showed that the longitudinal waves in the photospheric MBPs (with periods on the order of $130-440$ s) could be converted to kink modes at higher frequencies (higher by a factor of $\approx 2$), which were concluded to be the result of a $90^{\circ }$ phase difference encompassing opposite sides of the photospheric driver. Indeed, they found these waves are energetic enough (with an energy flux of $\approx 300$ kW/m$^2$) to heat the outer solar atmosphere (or accelerate the solar wind). Stangalini et al. (2014) provided observational evidence for excitation of kink modes in small photospheric magnetic elements as a result of granular buffeting. Their relatively long time-series of images (of about 4 h) from Hinode/NFI, and the use of EMD technique, allowed Stangalini et al. (2014) to reveal hints about the mechanisms of excitation of low frequency kink oscillations in small-scale magnetic tubes through their sub-harmonic response. Indeed, the probability density function of the periodicities of horizontal oscillations in a large sample of magnetic tubes, revealed several peaks in the statistical distribution corresponding to sub-harmonic oscillations with periods multiple of a fundamental one of $\approx 7.6$ min, which is comparable with evolution time of granular cells. Furthermore, the application of EMD approach on horizontal-velocity fluctuations of small (low) chromospheric MBPs, seen in SST Ca ii H images, led Stangalini et al. (2017) to find an elliptic polarization of the velocity vector associated to the low-frequency (smaller than 5–6 mHz) kink oscillations. The Ca ii H MBPs were showed to more freely move around, in a helical motion (while fluctuating transversely) compared to their photospheric counterparts (bounded to the granular flows). The left panel of Fig. 63 schematically illustrates such a polarized kink wave where the superposition of both helical motion and transverse kink waves co-exist in the same flux tube. The power spectra of the x and y components of the horizontal velocity of a Ca ii H MBP, as well as coherence spectrum between the two components, are plotted on the top-right panel (with solid red, dashed blue, and open circles, respectively), indicating the presence of higher frequencies, in addition to the larger peaks in the lower end. The helical motion, which was characterized from a phase relationship between the two components of the horizontal velocity, can also be visualized through the plot of the vector velocity on the bottom-right panel of Fig. 63, where a rotation in the displacement direction of the MBP is observed.

Fig. 63

Images reproduced with permission from Stangalini et al. (2017), copyright by AAS

Left: A cartoon illustrating a low-frequency helical displacement superimposed on a high-frequency kink wave in the solar chromosphere. Top right: Power spectra of the x and y components of the horizontal velocity of a Ca ii H MBP (solid red and dashed blue lines, respectively). The coherence spectrum of the two components are also plotted with the red open circles. Bottom right: Vector horizontal-velocity (both direction and magnitude) as a function of time. The helical motion (as the rotation of the velocity direction) is evident.

Transverse kink waves have also been studied in small-scale structures in sunspots. Pietarila et al. (2011) identified kink waves in dynamic fibrils (in the immediate vicinity of a sunspot; from Ca ii 8542 Å observations with SST) with periods of $\approx 135$ s. More recently, Morton et al. (2021) used high-resolution observations in Ca ii 8542 Å spectral line from SST/CRISP to demonstrate that transverse waves also pervaded the sunspot super-penumbral fibrils (in the solar chromosphere). They interpreted the oscillations as MHD kink modes with periods and propagation speeds on the order of 754 s and 25 km/s, on average, respectively. The velocity amplitudes (with an average of $0.76\pm 0.47$ km/s) were found to increase with distance from the umbral center by about 80%, as illustrated in Fig. 64. Morton et al. (2021) speculated this variation as, possibly, a result of a density decrease along the fibrils as the super-penumbra is extending to higher atmospheric heights while moving away from the umbra until reaching its highest magnetic-canopy point and returning to the surface. Thus, considering the field topology (in the chromosphere) is an important key when interpreting the observations, particularly, in intensity images in which the projection effects cannot be directly realized. Morton et al. (2021) also discussed a number of possible excitation mechanisms (for the transverse oscillations), namely, convection driven, reconnection, and mode conversion, of which, they found the latter to be more convincing. Oscillations in these sunspot’s small-scale structures may be different when compared to other chromospheric features due to a number of reasons, e.g., the very strong magnetic fields of the sunspots.

Fig. 64

Images reproduced with permission from Morton et al. (2021), copyright by the authors

Left: a sample sunspot super-penumbral fibril, along which transverse kink waves have been identified. Right: average velocity amplitudes and periods of the transverse oscillations as a function of distance to the umbral center. The error bars indicate the standard deviation of distributions of the parameters.

As has been discussed at length, observed frequencies in the lower solar atmosphere center around a range between 2 and 10 mHz. Observations with the Vacuum Tower Telescope (VTT) led Volkmer et al. (1995) to find relatively high-frequency oscillations in horizontal motions (with frequencies of $\approx 10$ mHz), consistent with kink modes, and in Doppler velocity (with frequencies that peaked at $\approx 8$ mHz) in small-scale structures in a plage region in the photosphere. The Doppler velocities were computed from Stokes-V profiles of the Fe i 630.15 nm spectral line. High-resolution observations from SST led Lin et al. (2007) to identify the signature of propagating kink oscillations (in both intensity and Doppler velocity) along numerous thin, thread-like structures in a H$\alpha$ filament. The 3–9 min perturbations found to travel along the small-scale structures with an average phase speed of 12 km/s. Using high-spatial resolution with DST/ROSA, and in agreement with numerical simulations, Jess et al. (2012c) reported upwardly propagating longitudinal magnetoacoustic waves in photospheric MBPs with periods in the range 100–600 s. They also found standing waves at shorter periods in about 27% of their MBPs. By employing time-series of slit-jaw images from IRIS (in 279.6 nm, 133 nm, and 140 nm channels), Zeighami et al. (2020) found 2–5.5 min intensity oscillations in small MBPs, propagating from the chromosphere to transition region with phase speeds ranging from 30 to 200 km/s.

More recently, Guevara Gómez et al. (2021) identified both kink and sausage modes in small, likely magnetic, bright points from pioneering observations with ALMA (at 3 mm), with periods on the order of 60 s, on average, for the transverse oscillations, and periodicities of about 90 s for the brightness temperature and size fluctuations. Although the exact heights of formation of these observations are still unclear, there have been indications to suggest that the ALMA Band 3 observations represents a wide range of heights, mostly from the mid-to-high chromosphere, but there may also be contributions from the lower chromosphere, and possibly the upper atmosphere (Wedemeyer et al. 2020; Jafarzadeh et al. 2021). Thus, it is difficult to conclude at this moment where these small structures reside, although it is highly likely to be in the chromosphere. The high-frequency oscillation reported by Guevara Gómez et al. (2021) are comparable to those previously found in the low-to-mid chromosphere.

The high frequencies observed in the photosphere and the lower/middle chromosphere ($> 10$ mHz) have been reported less frequently in the upper chromospheric fibrillar structures (De Pontieu et al. 2007a; Kuridze et al. 2012; Morton 2012; Morton et al. 2013, 2014), which however, were observed at different resolutions (and with different properties) compared to those seen at lower heights. These could be speculated as the result of wave energy dissipation (associated to those high frequencies) through the chromosphere. However, no clear observational evidence for such energy release has been found to date. We note that frequencies higher than 10 mHz were also observed by Morton et al. (2013, 2014) in fibrillar structures, however their mean values lie in lower frequencies. One exception is the high-frequency (on the order of 22 mHz) transverse oscillations that Okamoto and De Pontieu (2011) found in type ii spicules, however, this was postulated by the authors to be a result of the method they employed in their study. Direct evidence of kink wave damping in the solar chromosphere (i.e., in a Ca ii H spicule from Hinode/SOT) was provided by Morton (2014) when an initial rapid increase in the oscillation’s amplitude with height was followed by an amplitude decrease in the upper chromosphere. The conclusion of wave damping was reached by combining the amplitude variations with changes in width (of the spicule) and phase speed with height, while also comparing to theoretical models.

Oscillations in the chromospheric thread-like structures, including the off-limb Type i and Type ii spicules, and the on-disk counterparts of the latter, so-called rapid blue-/red-shifted events (RBEs/RREs; Rouppe van der Voort et al. 2009), have also been reported in a number of studies from both ground-based and space-born observing facilities. By exploiting joint observations of the lower solar atmosphere with SST and IRIS (Rouppe van der Voort et al. 2020) and the help of MHD simulations, Martínez-Sykora et al. (2017) described the generation of spicules as a result of magnetic tension and ion-neutral interactions. These authors found that impulsive release of the magnetic tension to excite Alfvén waves in these small-scale thread-like structures. Sekse et al. (2013b) identified longitudinal, transversal, and torsional oscillations in numerous RBEs/RREs from H$\alpha$ and Ca ii 8542 Å observations with SST. The three types of oscillations were found to propagate with velocity amplitudes on the order of 50–100 km/s, 15–20 km/s, and 25–30 km/s, respectively. Later, Rouppe van der Voort et al. (2015) speculated that bright features around their extended network regions observed in IRIS 1330 Å and 1400 Å slit-jaw images could be heating signatures associated to (waves in) H$\alpha$ RBEs and/or RREs from their coordinated observations with SST.

Another important interaction between small-scale magnetic concentrations and the convective motion is in the form of vortices at the solar surface (Steiner and Rezaei 2012). The presence and and properties of vortex motions in the solar photosphere and in the chromosphere have been studied from both observations and numerical simulations (e.g., Bonet et al. 2008; Wedemeyer-Böhm and Rouppe van der Voort 2009; Steiner et al. 2010; Shelyag et al. 2011; Park et al. 2016; Shetye et al. 2019; Yadav et al. 2020; Silva et al. 2020; Khomenko et al. 2021). Of particular interest is that the vortex flows can excite a variety of MHD wave modes, including torsional (Alfvén) waves, at small scales, as the magnetic field lines are frozen in the plasma in the lower photosphere (Fedun et al. 2011b; Tziotziou et al. 2020).

Jess et al. (2009) provided the first observational evidence of the torsional (Alfvén) waves, detected as full-width half-maximum oscillations in a small MBP through the lower solar atmosphere (with periods on the order of 126–700 s). They estimated an energy flux of $\approx 15{,}000$ W/m$^2$ carried by these waves. Later, Morton et al. (2013) demonstrated the excitation of incompressible kink modes by vortex motions of strong photospheric magnetic concentrations whose chromospheric counterparts showed quasi-periodic torsional motions. In addition, they identified transverse waves in the chromospheric fibrillar structures, connected to the magnetic concentrations, to be driven by the torsional motion. Using MURaM radiation-MHD simulations (Vögler et al. 2005), Yadav et al. (2021) discussed the formation mechanism of the small-scale vortices and showed how they can heat the solar chromosphere, though propagation of torsional (Alfvén) waves. Particularly, they showed that small-scales vortices are produced as a result of cascading in the relatively larger scales (residing in the interganular lanes in the photosphere) due to the turbulent nature of the plasma. That is, the twisted flux tubes create turbulence in the chromosphere, where the magnetic-field pressure dominates that of gas, by co-rotating the surrounding plasma. It is worth noting that some of the small features they found in their simulations (with diameters of 50–100 km in the photosphere; 100–200 km in the chromosphere) cannot yet be resolved in observations from currently available instruments.

According to MHD wave theory, an infinite number of wave modes may co-exist in the same magnetic structure, where phase mixing and mode coupling may also occur (Verth and Jess 2016). However, many fundamental wave modes (specifically the higher order modes) and their coupling/interaction, particularly at small scales, have been difficult to identify in observations. Stangalini et al. (2013b) reported interaction between transverse and longitudinal waves in small magnetic elements from both observations (with Sunrise/IMaX) and MHD simulations (with MURaM). They particularly found a $90^{\circ }$ phase difference between transverse oscillations (with frequencies larger than 10 mHz) and longitudinal (velocity) perturbations characterized by frequencies smaller than 7–8 mHz. The interaction between the two type of MHD waves were, however, found to take place with a high confidence level at periods shorter than 200 s.

High spatial and temporal resolution H$\alpha$ observations from SST/CRISP provided Sharma et al. (2017b) three-dimensional velocity vectors to identify MHD kink modes in spicules from two approximately perpendicular angles. Furthermore, Sharma et al. (2018) used the same dataset to inspect coupling between various MHD wave modes in the off-limb thread-like structures. In this regard, they also explored time variations of longitudinal, cross-sectional width, photometric, and azimuthal shear/torsion parameters at selected spicules, that were concluded to be coupled over the period scale, supported by mutual phase relationships. In particular, they found that the nonlinear kink waves (identified as the displacement of the spicule’s axis in both the plane-of-sky and Doppler directions) were coupled with the longitudinal (field-aligned) flows. These led Sharma et al. (2018) to explain the coupling of the independent wave modes in the spicules as a result of a single pulse-like driver following a twist. Figure 65 visualizes a 3D structure of a spicule studied by Sharma et al. (2018), where the coupled transverse and width oscillations (top row) as well as transverse and azimuthal shear components (bottom row) are shown in four time steps.

Fig. 65

Image reproduced with permission from Sharma et al. (2018). An animation of this figure is also available

Visualization of the coupled MHD wave modes in a spicule, constructed from high-resolution observations with SST/CRISP. The four columns illustrate the 3D structure at different time steps indicated on the top. Top row: coupled transverse and width with intensity. Bottom row: transverse and azimuthal shear components.

While ubiquitous rapid (supersonic) and high-frequency intensity fluctuations in the chromosphere and transition region, observed with the rocket-borne Chromospheric Lyman-Alpha Spectropolarimeter (CLASP; Kano et al. 2012) instrument, were explained as MHD fast-mode waves (Kubo et al. 2016), they were later attributed to both waves and jets (i.e., small-scale transient features) from joint observations with CLASP and IRIS (Schmit et al. 2020). In the latter study, the authors found non-linear wave propagation in the core of plages and linear propagation of fluctuations as a result of non-recurrent jet-like features. Moreover, using an unprecedented high temporal cadence of 0.3 s with CLASP in hydrogen Ly$\alpha$ 1216 Å line, Yoshida et al. (2019) found high-frequency oscillations (of the Doppler velocity) in the early phase of a spicule evolution, with period and propagating speed on the order of 30 s and 470 km s$^{-1}$, respectively.

Other excitation mechanisms have been proposed for the observation of magnetoacoustic waves at small scales. Of particular interest is small magnetic reconnection which have been thought to be the driver of kink modes (He et al. 2009; Ebadi and Ghiassi 2014). Furthermore, through numerical studies, Kato et al. (2011) proposed a new mechanism called “magnetic pumping” to excite upwardly propagating slow modes in magnetic flux concentrations. They showed that the convective downdrafts around a flux tube can eventually result in pumping downflows inside the tube, hence, creating magnetoacoustic oscillations

Table 1 summarizes the average properties of some of the various MHD waves in small-scale magnetic structures, reviewed in this section. As evidenced, the mean periods and phase speeds measured in different studies have wide ranges on the order 34–754 s and 13–270 km s$^{-1}$, respectively. We note that comparison of these wave characteristics, from different studies, should be performed with great caution. These may include results obtained for the same wave types/modes and/or in similar structures (e.g., in magnetic bright points or fibrillar structures). Such values may not always be one-to-one comparable due to various reasons, such as structures with different spatial/temporal scales residing in different geometric heights and/or different solar environments, as well as some measurement effects. Estimating accurate formation heights is a challenging task, even when the observations are made at similar wavelengths. For instance, the spectral resolution or the width of filters employed, and/or the level of magnetic flux contained can play important roles in observations of different geometric heights, not only on average, but also across the field of view or along the magnetic structures. The choice of analysis approaches is another important factor in the reported (often) mean values (due to, e.g., some selection effects). Furthermore, the spatial and temporal resolutions of the observations as well as the length of the time series can limit, e.g., the identified structures (which are found in a variety of spatial and temporal scales, though with similar names), and the range of detectable frequencies. Hence, the wave characteristics reported in the literature (using different observations) may not necessarily represent waves in the same structures, same geometric heights, and/or same solar regions.

Table 1 Average (or median) period (T), phase speed ($v_{\text{ph}}$), and energy flux ($F_E$) of MHD waves in small-scale magnetic structures observed in the lower solar atmosphere

3.4.2 Magnetic-field perturbations in small-scale magnetic structures

Measurement of the magnetic fields at small scales has been a challenge due to various reasons, particularly, due to the fact that most of them are spatially unresolved. Furthermore, the photospheric linear polarization signals are often weak (i.e., on the order of the photon noise level) in the quiet Sun where the small-scale magnetic structures reside. Therefore, computations of the full vector magnetic field at these structures are rare (Bellot Rubio and Orozco Suárez 2019) and may lead to incorrect field parameters, such as field inclination angles, when, e.g., traditional Stokes inversions, are employed (Borrero and Kobel 2011, 2012; Jafarzadeh et al. 2014b). Such measurements at chromospheric heights are even more challenging, due to the smaller magnetic flux (Lagg et al. 2017).

As a result, it is often difficult, if not impossible, to accurately constrain the vector magnetic fields of small-scale magnetic elements as they weave their way from the base of the photosphere through to the chromosphere and beyond. A novel way of uncovering the magnetic field information associated with such small-scale features is to examine them off-limb, where the background polarimetric signals contaminate the Stokes profiles to a lesser extent. Hence, viewing structures, like spicules, against the black background of space is a compelling way of measuring their small-amplitude magnetic signals (e.g., López Ariste and Casini 2005; Socas-Navarro and Elmore 2005; Trujillo Bueno et al. 2005; Centeno et al. 2010). In particular, Utilizing high-resolution Ca ii 8542 Å limb observations acquired with the SST, Kriginsky et al. (2020) employed the weak field approximation to radiative transfer equations to infer the line-of-sight component of the magnetic field in a multitude of different spicule types. Magnetic field strengths on the order of 100 G were found ubiquitously along the spicule structures, with little difference found between spicules embedded close to active regions and those associated with the quiet Sun (Kriginsky et al. 2020).

Moving on from the weak field approximation, Kuridze et al. (2021) harnessed a new version of the non-LTE NICOLE inversion code to invert a prominent off-limb spicule captured in the Ca ii 8542 Å spectral line by the SST. By considering true geometry effects through the inclusion of vertical stratification, Kuridze et al. (2021) were able to provide a semi-empirical model for the specific spicule structure examined, which consisted of a uniform temperature of 9560 K, coupled with an exponential density decrease as a function of atmospheric height, providing a density scale height on the order of 1000–2000 km. These results are consistent with those previously deduced by Beckers (1968), Alissandrakis (1973) and Krall et al. (1976), albeit with more modern non-LTE considerations invoked. However, the spicule studied by Kuridze et al. (2021) is an interesting structure with atypical characteristics. Specifically, the spicule demonstrated a clear inverted Y-shaped base, consistent with anemone jets driven by magnetic reconnection in the lower solar atmosphere (e.g., Yokoyama and Shibata 1995; Shibata et al. 2007; He et al. 2009). In addition, the spicule structure reached atmospheric heights of $\sim 10$ Mm above the surface and was clearly visible in the far blue-wing of the Ca ii 8542 Å spectral line, suggesting the feature may be similar to dynamic type ii spicule events (Martínez-Sykora et al. 2018). However, the lifetime of the structure examined by Kuridze et al. (2021) was $> 20$ min, which is not consistent with the shorter duration lifetimes ($\sim 50{-}150$ s) of traditional type ii spicules (Pereira et al. 2012, 2016; Sekse et al. 2012, 2013a). As such, the spicule feature examined by Kuridze et al. (2021) may be more closely related to the macrospicules initially observed as cool plasma in the He ii 304 Å spectral line by Skylab (Bohlin et al. 1975), and later as vibrant O v emission by the Solar Ultraviolet Measurements of Emitted Radiation (SUMER; Wilhelm et al. 1995; Wilhelm 2000) spectrograph onboard the Solar and Heliospheric Obervatory (SOHO; Domingo et al. 1995) spacecraft.

Small-scale magnetism can also be studied in the upper-chromospheric He i 10830 Å spectral line. Here, it is possible to make use of the Hanle effect to deduce magnetic field information since it is more sensitive (compared to the traditionally employed Zeeman effect) to the weaker magnetic fields present in small magnetic elements (Stenflo 1994; Lin et al. 1998). Similar to the work of Kriginsky et al. (2020) and Kuridze et al. (2021), many studies have attempted spectropolarimetric inversions of off-limb spicular features to uncover the magnetic field variations away from the solar disc (e.g., Trujillo Bueno et al. 2005; Centeno et al. 2010; Orozco Suárez et al. 2015). Utilizing He i 10830 Å diagnostics, Trujillo Bueno et al. (2005) uncovered spicule magnetic fields as low as 10 G. This highlights the difficulties when attempting to uncover magnetic field perturbations arising from propagating wave phenomena. Even with a relatively large 10% amplitude variation, this only equates to a $\pm 1$ G fluctuation in the associated magnetic field strength of the spicule. As a result, it becomes statistically challenging to reliably quantify such minuscule variations in the magnetic field strength, especially with these plasma parameters inferred from a multitude of spectral lines, often with weak Stokes Q/U components. Hence, measurements of oscillations in individual components of the polarization signals has predominantly been limited to the (often dominant) circular polarization (Stokes V) component.

Only recently, small-scale kilogauss magnetic elements could spatially be fully resolved (Lagg et al. 2010), not only because the high-spatial resolution provided by the 1-m Sunrise balloon-borne solar telescope, but also due to the seeing-free observations (and the high precision of the IMaX spectropolairemeter; Martínez Pillet et al. 2011) which in turn resulted in higher polarization signal-to-noise compared to those normally achieved with similar ground-based instruments currently available. Thus, fluctuations of polarization signals could also be detected in small magnetic elements observed by Sunrise (Jafarzadeh et al. 2013). With a relatively lower spatial resolution from Hinode, but also from seeing-free data, Utz et al. (2013b) were also able to measure the kilogauss field strength at small magnetic elements. However, it should be noted that the definition of ‘small’ features may vary from one study to another, with various spatial sizes reported on.

Martínez González et al. (2011) presented magnetic-field oscillations (from Sunrise/IMaX observations) in a very quiet photospheric area whose field strength did not exceed 500 G. They concluded that the oscillations were not associated to oscillatory modes of magnetic concentrations, but rather, to buffeting of the magnetic-field lines by granular flows. It is worth noting that they found two different prominent period ranges, one corresponding to magnetic flux density patches of $10^{16}$–$10^{17}$ Mx (with a period range of 4–11 min) and one associated with more intense patches of $10^{18}$ Mx (with periods on the order of 3–5 min). Thus, the small magnetic patches studied by Martínez González et al. (2011) were representative of relatively weak magnetic flux concentrations that preferentially emerge within granules (Martínez González et al. 2007; Centeno et al. 2007; Orozco Suárez et al. 2008; Martínez González and Bellot Rubio 2009; Danilovic et al. 2010; Kianfar et al. 2018). On the other hand, Martínez González et al. (2011) found clear evidence of p-modes oscillations (with a 5 min periodicity) in the Doppler velocity, as well as both continuum and line-core intensities (of the Fe i 525.02 nm spectral line). Interestingly, they found a 180 degree phase difference between fluctuations in the continuum and line-core intensities, as well as an anti-phase between Doppler velocity and continuum intensity perturbations. Such anti-correlations between Doppler velocity and intensity in small-scale magnetic structures were also found by Campbell et al. (2021) in high-resolution observations with GREGOR (Schmidt et al. 2012).

Recently, Norton et al. (2021) examined the identification of magnetic-field perturbations from SDO/HMI observations, in various solar regions. They found that except in the umbra, almost no oscillatory power could be detected in the field-strength oscillations in, e.g., the quiet-Sun and plage regions. No oscillatory signature could also be identified in the field inclination and azimuth.

Using high spatial and temporal resolution observations with SST/CRISP, Keys et al. (2020) detected rapid (within 33–99 s) magnetic-field amplification (by a factor of $\approx 2$ on average) in numerous MBPs whose field strengths followed a bimodal distribution (Keys et al. 2019). With the help of numerical simulations, Keys et al. (2020) found that the field amplification could possibly be explained as a result of convective collapse (Spruit 1979), as the most frequent process, thus the decrease in size of the MBPs is accompanied by amplification of the field strength. This process is similar to that responsible for excitation of sausage modes. Furthermore, the authors find evidence, albeit less frequently, for granular compression of the inter-granular lanes leading to magnetic field amplification in the MBPs. Although the ,mechanisms leading to field amplification are somewhat similar to convective collapse with both displaying a decrease in size, the authors pick out differences that imply a distinct process occurring in this case for amplification due to granular expansion into the lanes. Like the case with convective collapse, the authors note that this process is similar to that responsible for excitation of sausage modes.

In addition to the level of polarization signals (compared to the noise level) which is necessary for a reliable detection, inferring physical parameters from Stokes inversions seems to be challenging when waves pass through the solar atmosphere. Keys et al. (2021) did an experiment with a synthesized dataset for the Fe i 6301 Å and Fe i 6302 Å line pair and found that, e.g., the Doppler velocities, could not be returned accurately by the inversion codes, after synthesis with NICOLE (Socas-Navarro et al. 2015) and inversion with SIR (Ruiz Cobo and del Toro Iniesta 1992) (compared to their initial values in the numerical simulations), at the presence of an upwardly propagating wave in a thin flux tube. This was explained as waves perturb the atmosphere over a smaller height range compared to that sampled by the spectral lines. Thus, development of inversion codes in this regard is crucial, otherwise, the prevalence of waves and oscillations in the solar atmosphere may largely influence the inferred parameters from inversions. It is worth nothing that some recent advancements in development of powerful (multi-line) non-LTE inversion codes (e.g., Riethmüller et al. 2017; Milić and van Noort 2018; de la Cruz Rodríguez et al. 2019; Ruiz Cobo et al. 2022), and other similar approaches (Centeno 2018; Asensio Ramos and Díaz Baso 2019), will significantly improve these measurements, particularly in the quiet regions through the entire lower solar atmosphere.

Although, the current pioneering instruments, with enhanced spatial resolutions and spectropolarimetric sensitivities over the past decade, have advanced our understanding around the magnetic-field oscillations at small scales (and in the quiet-Sun regions), a clear detection of such perturbations remains difficult to date. It is foreseen that the next generation spectropolarimetric instruments, including those on DKIST and Sunrise III, together with new generations of Stokes inversion codes, will revolutionize our vantage points of magnetic-field oscillations in the lower solar atmosphere.

4 Future directions and progress

Section 3 highlights key advancements and accomplishments made by the global solar physics community over the last number of years. Importantly, these achievements also reveal areas where our collective understanding is currently lacking. Here, we will pinpoint specific topics that we believe are within scientific reach following the realization of next-generation observatories, modeling techniques and analysis tools.

It has been shown in a multitude of studies that wave signatures manifesting in the lower solar atmosphere may be the result of upward wave propagation from sub-photospheric layers (e.g., Löhner-Böttcher and Bello González 2015; Chae et al. 2017, 2019; Griffiths et al. 2018; Kayshap et al. 2018; Yurchyshyn et al. 2020; Zeighami et al. 2020, to name but a few recent examples). However, while it seems that global eigenmodes may be responsible for many of the observed wave signatures, they are unable (by themselves) to account for all of the perceived signals found in lower atmospheric wave guides, which often require additional (often unknown) perturbations, pulses, and/or excitation sources to explain (e.g., Zhao et al. 2016; Stangalini et al. 2018, 2021a). Therefore, what other mechanisms are responsible for the complex wave signatures captured in high-resolution photospheric and chromospheric data sequences? This is where Fourier filtering techniques (see, e.g., Sect. 2.3.1) can play an important role, since they have the ability to disentangle small-scale wave perturbations from macroscopic flows and dominant MHD modes. Indeed, wavelet filtering and pixelized wavelet techniques can also provide useful diagnostic potential, especially if transient and/or rapidly developing wave signatures are present (Antoine et al. 2002; Sych and Nakariakov 2008; Sych et al. 2021). Therefore, tailored Fourier and/or wavelet filtering algorithms will be of paramount importance to extract the smallest amplitude fluctuations, particularly those associated with higher order MHD modes, from the dominant (with often orders-of-magnitude larger power) wave signals.

The power of novel visualization tools, such as $B{-}\omega$ diagrams (see Sect. 2.7), cannot be overstated. In Fig. 30, the suppression of p-modes towards the boundary of a magnetic structure is clearly evident, as is the growth of distinct eigenmodes once magnetic fields synonymous with the solar pore are isolated. Such techniques go well beyond traditional one-dimensional Fourier power spectra, and allow not only the wave behavior to be more closely scrutinized as a function of other local plasma properties, but also provide impressive visualizations that allow the work to be more readily disseminated. We expect such powerful techniques to be widely exploited by the solar physics community in years to come. In particular, $B{-}\omega$ diagrams may be harnessed in longer-duration synoptic studies to see whether driving mechanisms in magnetic features (sunspots, pores, etc.) vary over the course of a solar cycle. As such, long-duration data sequences will provide excellent frequency resolutions that may help isolate true eigenmodes trapped within the boundaries of such magnetic structures.

Observations have shown systematic differences in the chemical abundances found in the corona and in the photosphere. In particular, in closed-loop systems, low ($<10$ eV) first ionization potential (FIP) elements are occasionally more abundant in the corona than in the photosphere, by a factor ranging from 2 to 4. In contrast, the plasma composition in open magnetic field regions remains practically unfractionated between the lower and upper solar atmosphere (Sheeley 1995, 1996; Baker et al. 2013, to mention but a few examples). Interestingly, the FIP bias is also observed in the solar wind (Schwadron et al. 1999; Laming et al. 2019) and could, therefore, be used to infer the magnetic connectivity in the heliosphere, and possibly help to identify the sources of the solar wind itself (Brooks et al. 2015). This is of primary importance to investigate and identify the solar wind acceleration mechanisms, which is one of the fundamental questions that the Solar Orbiter mission is being specifically designed to address through a combination of both remote sensing and in-situ instruments. It was predicted from theoretical modeling that the FIP bias in the solar corona could be due to the presence of magnetic perturbations which, through the ponderomotive force, are responsible for the plasma fractionation (Laming 2015). This was recently confirmed thanks to a combination of ground-based and space-borne observations, which identified magnetic oscillations at chromospheric heights magnetically linked to the coronal locations where the FIP bias was observed (Baker et al. 2021; Stangalini et al. 2021a).

Theoretical and numerical modeling work has continually speculated that the lower solar atmosphere should be replete with the entire assortment of MHD wave modes: slow and fast magnetoacoustic modes, plus Alfvén waves (e.g., Cally 1983, 2017; Cally et al. 1994, 2016; Khomenko and Collados 2006; Cally and Goossens 2008; Terradas et al. 2011; Hansen and Cally 2012; Khomenko and Cally 2012; Cally and Moradi 2013; Moreels et al. 2015; Arber et al. 2016; Leonard et al. 2018; Cally and Khomenko 2019; González-Morales et al. 2019; Pennicott and Cally 2019; Raboonik and Cally 2019, to name but a few examples). Indeed, observations have shown evidence for ubiquitous slow mode waves (Grant et al. 2015; Kanoh et al. 2016; Tsap et al. 2016; Jess et al. 2017; Chen et al. 2018; Kang et al. 2019), and indeed even (more challenging to detect) Alfvén modes (Jess et al. 2009; De Pontieu et al. 2014a; Srivastava et al. 2017). Unfortunately, the universal existence of fast magnetoacoustic modes in the lower solar atmosphere is less well documented. Morton et al. (2012) revealed evidence for fast incompressible MHD wave modes in chromospheric fibrils and mottles, which may be driven by photospheric flows and vortices (Morton et al. 2013; Murawski et al. 2018; Liu et al. 2019; Murabito et al. 2020). However, the unequivocal ubiquity of fast mode waves in the lower solar atmosphere has still not been verified.

Spectropolarimetric inversion processes are powerful techniques that allow key plasma parameters (magnetic fields, densities, temperatures, velocities, etc.) to be established as a function of optical depth. Examples of such software have allowed crucial plasma conditions associated with solar structures exhibiting wave activity to be uncovered (e.g., Beck et al. 2013b; Rouppe van der Voort and de la Cruz Rodríguez 2013; Houston et al. 2018, 2020; Grant et al. 2018; Felipe et al. 2019; Jess et al. 2020). Next-generation inversion routines, including the Stockholm inversion code (STiC; de la Cruz Rodríguez et al. 2019), the Spectropolarimetic NLTE Analytically Powered Inversion (SNAPI; Milić and van Noort 2018) code, and the Departure coefficients Stokes Inversion based on Response functions (DeSIRe; Ruiz Cobo et al. 2022) code offer powerful new approaches for accurately constraining the derived plasma profiles, including the simultaneous use of multiple spectral lines. This is especially important when lines may not strictly be formed under LTE conditions, for example, the Fe i photospheric absorption lines in the presence of UV overionization (Smitha et al. 2020). Recently, Riethmüller and Solanki (2019) showed that many-line inversions of photospheric spectropolarimetric data, particularly at short wavelengths (where the photon noise is considerably higher than that at longer wavelengths) can significantly improve the outputs compared to inversions of a few spectral lines only. In addition, van Noort (2012) and Asensio Ramos and de la Cruz Rodríguez (2015) have recently been developing spatially-coupled inversion routines, whereby the authors found that inclusion of the point spread function of the telescope, alongside the degree of spatial correlation between neighboring pixels, respectively, helps to minimize associated errors of the inversion outputs. As such, utilizing numerous spectral lines with differing magnetic field sensitivities (i.e., different Landé g-factors) spanning the base of the photosphere through to the upper extremities of the chromosphere will be required to converge spectropolarimetric inversion outputs to trustworthy values, which will be of paramount importance when attempting to benchmark small-scale seismological fluctuations throughout the lower solar atmosphere.

There are a myriad of complexities involved in inversion studies involving oscillatory phenomena. One such issue is that most inversion codes return the plasma parameters as a function of optical depth, which is distinctly different from the geometric atmospheric height. It is a challenging endeavor to convert optical depth to atmospheric height, since this conversion hinges upon the reliability of the input model atmosphere. Indeed, Ishikawa et al. (2018) recently showed the influence of different input atmospheric models on the inverted diagnostic outputs. This highlights the importance of ensuring the stationary (background) model is representative of the structure being investigated. Indeed, if wave-based perturbations are superimposed on top of this stationary input model, then it becomes difficult to disentangle what is a true fluctuation from what is just a simple deviation from the simplified input model. This complexity is even more pronounced when systematic effects associated with spectral fitting routines contaminate the inversion outputs (Asensio Ramos et al. 2016). The next logical step is to try and generate background model atmospheres that already contain stratified periodic fluctuations in, e.g., density, temperature, velocity, etc., which are synonymous with the wave features wishing to be evaluated. Of course, such a-priori knowledge of the dominant embedded wave modes may not always be available to the researcher before the inversions are performed. As such, making use of inversion processes based on machine learning and neural networks (e.g., Asensio Ramos and Díaz Baso 2019; Milić and Gafeira 2020; Socas-Navarro and Asensio Ramos 2021) can help expedite the inversion process, especially when dealing with large degrees of freedom. Ultimately, this will allow Stokes profiles that contain wave perturbations to be reliably inverted, hence minimizing uncertainties and providing the first step in converting optical depths into true geometric heights once compared to the accurate background model that is reflective of the observations acquired.

Another possible solution to this issue is the MHD-assisted Stokes Inversion (MASI; Riethmüller et al. 2017) code, which is based on the spectral syntheses of state-of-the-art MHD simulations and is a first step to being able to directly output plasma parameters as a function of true geometric scales. This work builds upon the legacy codes provided and documented by Molowny-Horas et al. (1999), Tziotziou et al. (2001), Berlicki et al. (2005), Carroll and Kopf (2008) and Beck et al. (2013a, 2015). As highlighted by Riethmüller et al. (2017), the MASI code provides a platform for inversions that give results consistent with the underlying MHD equations, which is likely to be of huge benefit to researchers working exclusively in wave perturbations in the lower solar atmosphere.

Yet another solution for inferring the magnetic-field vector throughout the lower solar atmosphere as a function of geometric height is the non-force-free fields extrapolations introduced by Wiegelmann et al. (2015, 2017). These extrapolations use high-spatial resolution photospheric field vector as the boundary condition for a magnetohydrostatic (MHS) model. This extrapolation code was specifically designed for accurate approximation of magnetic-field vector at heights below $\approx 2000$ km (i.e., the solar photosphere and the chromosphere) where the non-vanishing Lorentz force exists. Thus, to account for such a mixed plasma-$\beta$ environment, the MHS model self-consistently considers the pressure gradients and gravity forces through these regions of the atmosphere. We note that such MHS extrapolations are different from the traditional force-free field extrapolations which are mostly accurate in the solar corona where the Lorentz force vanishes (Wiegelmann et al. 2006, 2008; Wiegelmann and Sakurai 2012). By combining MHS constraints with Stokes inversions, Borrero et al. (2021) developed a new inversion code which is capable of retrieving the physical parameters in the solar photosphere as a function of geometric height. The new code makes use of an MHS solver and the Firtez-DZ code (a solver of the polarized radiative transfer equation in geometrical scale; Pastor Yabar et al. 2019; Borrero et al. 2019) which are based on three-dimensional MHD simulations of sunspots (Rempel 2012).

Recent work by Keys et al. (2021) showed the difficulties in constraining atmospheric parameters in an oscillating waveguide using commonly used inversion techniques on a simulation with propagating MHD waves. The authors utilized simple two-dimensional MHD models with known driver and atmospheric parameters. From this, the authors synthesized Stokes profiles for the 6301 Å and 6302 Å line pair and then employed the Stokes Inversion based on Response functions (SIR; Ruiz Cobo and del Toro Iniesta 1992) code to establish if the atmospheric parameters of the oscillation could be accurately returned. Results from a study by Vigeesh et al. (2011) on the same dataset showed the necessity of high-spatial resolution in resolving asymmetries in the Stokes parameters, therefore, Keys et al. (2021) degraded the spatial resolution to that typical of DKIST and the upcoming 4-m European Solar Telescope (EST; Collados et al. 2010, 2013; Quintero Noda et al. 2022). The results highlighted that inversion codes such as SIR, could return the atmospheric parameters fairly accurately within typical height of formation regions for the lines inverted. However, the authors note that the Doppler velocity is less accurately fit for perturbed spectra. This is an issue that could possibly be improved upon in the future.

This is a potential difficulty for inversion algorithms that utilize a nodal approach in minimizing the differences between the post radiative transfer profiles to the input models. Due to the fact that a spline interpolation is generally used to calculate the atmospheric parameters at grid points between nodes, there is a trade-off in selecting nodes (free parameters) in the inversion. Too few and it is unlikely that the perturbation will be accurately fit in the atmosphere. Too many nodal points could lead to potentially over-fitting the atmosphere when the spline interpolation is calculated (Bellot Rubio et al. 2000). Therefore, there is a risk in inversion algorithms that employ the nodal approach that an over-fit atmosphere is misinterpreted as a perturbation in the atmosphere. Furthermore, inversion algorithms generally assume the atmosphere is in hydrostatic equilibrium to make the computation manageable. Perturbations and oscillatory phenomena in the atmosphere will affect the validity of this assumption could affect the resultant output from the inversion. However, initial studies have shown that inversion methods can extract reliable wave behaviour. For instance, Nelson et al. (2021) estimated the magnitude of the magnetic field of a pore with SIR inversions and the strong-field approximation method. Both methods detected magnetic perturbations of similar amplitude and frequency, consistent with sausage modes. Despite this, it is an uncertainty that deserves more comprehensive study, particularly when impulsive, non-linear wave modes are under consideration (e.g., Houston et al. 2020).

To an extent, these issues could be alleviated by multi-line studies using inversion codes such as STiC and DeSIRe to invert the atmosphere at multiple heights. Studies ascertaining the success of such approaches in the context of wave propagation need to be performed, however, to establish any potential issues with this approach. Another solution to this problem was suggested by Keys et al. (2021), whereby the authors suggest adapting an inversion code such as CAISAR (or MASI) to take spectra from MHD simulations archives with propagating wave phenomena with known wave properties to return the atmospheric parameters in observations with wave phenomena present. A similar idea has already been implemented with IRIS data by Sainz Dalda et al. (2019). Here, the authors employed machine and deep learning techniques using a database of representative IRIS profiles of various solar features to ascertain the thermodynamic state in the upper photosphere and chromosphere. This method was found to be both accurate and significantly faster to process over more traditional techniques. A similar approach utilising MHD simulations as suggested by Keys et al. (2021) is not without issues, as it is constrained by the accuracy and diversity of profiles generated by the MHD simulations. This is highlighted further by Fleck et al. (2021), who find substantial differences in wave properties between the four commonly used MHD codes. Even considering these limitations, and keeping in mind all inversion algorithms will have limitations, it may be desirable to explore the possibilities of adapting existing tools, to perhaps provide a somewhat robust method of inverting data with potential wave phenomena present. Regardless, the issue of accurately constraining atmospheric parameters from inversions in an oscillating atmosphere is still an open question in the field. It is highly likely that either a dedicated inversion code is needed or an existing code is adapted to invert datasets with suspected wave activity for a more accurate inversion of the observations. This is an even greater requirement for small-scale dynamic features, such as MBPs.

Further issues may arise due to the dataset being inverted as well. That is, instrumental limitations may make it harder to retrieve accurate inversion results in the presence of upwardly propagating waves. One such example is with the use of Fabry–Pérot interferometers. With these instruments there will be a time lag as the instrument scans across a given line, which can result in a scan time of at least around 10 s for an individual line after data reduction. Obviously this will have implications on wave studies. One study by Felipe et al. (2018b) looked at the effect of scanning time in umbral flashes appearing in simulations obtained with the MANCHA code. The authors synthesized (and inverted) the 8542 Å line with the NLTE code NICOLE (Stangalini et al. 2015). The authors performed the inversions on the instantaneous synthesized Stokes parameters (with a single time step) as well as ‘synthetic scanned’ Stokes parameters, where line scanning similar to instrument scans was introduced to the profiles to simulate actual observations. The authors report that, for profiles simulating a line scan, the inversions reasonably infer the atmosphere prior to and after the umbral flash is fully developed. However, during the early stage of the umbral flashes, issues appeared due to the short time for changes in the Stokes profiles in comparison to the temporal cadence of the scans. The Stokes-V profiles in this case were similar to those for the instantaneous profiles, but with a remarkably different intensity profile associated to them. Such a difference can complicate the interpretation of the spectropolarimetric data, with the authors estimating that approximately 15% of profiles in flashing regions are affected by this issue.

A subsequent study by Felipe and Esteban Pozuelo (2019) looked at the effect of wavelength sampling on inversion results for wave studies. Using a similar approach to previous work, the authors synthesized the 8542 Å line from a simulated umbral flash event, and degraded the spectral resolution (i.e., the wavelength sampling) before inverting the data with NICOLE and analyzing. The authors find that the vertical magnetic field inferred from the inversions is more accurate with finer wavelength sampling. However, finer sampling means an increase in scan time, which will miss sudden changes in the profiles given the dynamic nature of the umbral flash. The authors conclude that sampling positions should be selected with observing goal in mind: studies of the magnetic field should prioritize spectral resolution while studies of the temperature and Doppler velocity should prioritize faster scan times.

In essence, the work of Felipe and Esteban Pozuelo (2019) is intrinsically linked to the line scan time study (Felipe et al. 2018b), in that spectral resolution will affect the scan time for a line. However, the study highlights some key challenges facing the wave community in considering and interpreting the results from inversions. Key information can be lost and/or misinterpreted depending on the dataset and instrument employed. In terms of the limitations investigated in these studies linked to Fabry–Pérot interferometers, the next generation of instruments and facilities should alleviate the issues in obtaining reliable inversion results from the data. For example, instruments employing detectors with lower noise levels should attain faster scan times. Likewise, instruments such as integral field units can circumvent some of these issues as they render scanning superfluous. These solutions are not without their issues, however, as faster scan times from better detectors will likely lead to more complex, multi-line scans being used by observers (thus making the scan time issue emerge again). Also, integral field units have a limited field of view, which may impact certain wave studies.

It should be noted as well, that these studies focused on umbral flashes in the 8542 Å line. Therefore, the focus is on a specific phenomena in the chromosphere. Further study is likely warranted for the effect of scan time and spectral resolution on inversion results for perturbations in the photosphere as well for other waveguides and oscillations. It is likely that similar issues will be reported, however, such studies would still be beneficial to determine optimal observing parameters for different features within the context of inversions. Realistically, a combination of instrument improvements and observing sequences with more modern detectors may be needed to fully understand the results from inversions in the context of wave studies in the lower solar atmosphere. Dedicated inversion algorithms for wave studies would likely complement such spectropolarimteric datasets and, therefore, give a deeper understanding and clearer understanding of the atmospheric parameters in the presence of propagating disturbances across a range of structures.

Of course, in order to obtain high-precision full Stokes spectropolarimetric signals that are suitable for the study of wave activity, instrumentation is required to capture such spectra with high cadences and low noise levels. Traditional instrumentation, including slit-based spectropolarimeters (e.g., the Facility Infrared Spectropolarimeter [FIRS] at the DST, Jaeggli et al. 2010; the Visible Spectropolarimeter [ViSP] on DKIST, de Wijn et al. 2012; and the TRI-Port Polarimetric Echelle-Littrow [TRIPPEL] spectrograph at the SST, Kiselman et al. 2011) and Fabry–Pérot spectral imagers (e.g., CRISP, IBIS, and the GREGOR Fabry–Pérot Interferometer [GFPI], Puschmann et al. 2012, 2013; and the Near Infrared Imaging Spectropolarimeter [NIRIS] at Big Bear Solar Observatory’s New Solar Telescope, Cao et al. 2012), have long been the main acquisition systems for such observations. However, these instruments often run into difficulties when attempting to achieve the ‘trifecta’ of high acquisition cadences, spectral precision, and spatial resolutions, all at the same time. For example, it is difficult to maintain high cadences for slit-based spectrographs when scanning their one-dimensional slit across a large field-of-view. Similarly, the simultaneity of spectral profile shapes can be lost if the wavelength scanning time is not fast enough in the case of Fabry–Pérot interferometers (Felipe et al. 2018b).

The next generation of solar instrumentation is attempting to address the inherent weaknesses of previous observing systems. Many are turning their attention to Integral Field Units (IFUs), which can help provide simultaneous high-precision spectra across a two-dimensional field of view (Iglesias and Feller 2019). IFUs typically come under three distinct guises: (1) fiber-fed bundles, (2) image slicers, and (3) microlens arrays, which are summarized in Fig. 66. At present, one drawback of IFUs is the reduced field-of-view sizes, which is a natural consequence of requiring the imaging detector to record spectra across two dimensions (i.e., x and y domains) simultaneously. Multi-slit variations (e.g., operating the FIRS instrument at the DST in a multiple-slit configuration), which are also shown in Fig. 66, can be viewed as an in-between solution that offers better cadences that single slit spectrographs, but lacks the true simultaneity of IFUs.

Fig. 66

Image reproduced with permission from Iglesias and Feller (2019), copyright by the authors

Six different spectral mapping techniques used in solar spectropolarimetry. The input cube ($x,y,\lambda$) at the top can be mapped on to the detector (bottom row) through a number of different mechanisms. Towards the left-hand side of the figure are more traditional methods, including scanning one-dimensional slit-based spectrographs and Fabry–Pérot interferometers, while the right-hand side of the panel depicts the three different types of IFU technology currently available. A multi-slit configuration (middle of the figure) can be considered an in-between solution.

Much progress has been achieved in recent years with regards to IFU development. The fiber-fed Diffraction Limited Near Infrared Spectropolarimeter (DL-NIRSP; Elmore et al. 2014) is currently in the commissioning stage at the DKIST, providing spatially resolved spectra of infra-red features such as Ca ii 8542 Å and He i 10830 Å. Prototypes of image slicers (e.g., the Multi-Slit Image Slicer based on collimator-Camera, MuSICa; Calcines et al. 2013), and microlens-fed spectrographs (MiHi; Jurčák et al. 2019) have also been trialed with great success. In addition, a new hyper-spectropolarimetric imager is being developed for the Indian National Large Solar Telescope (NLST; Hasan et al. 2010; Dhananjay 2014), which is due to be built close to the Chinese border in the Merak region of northern India. This prototype instrument, named the Fibre-Resolved opticAl and Near-infrared Czerny-Turner Imaging Spectropolarimeter (francis), is specifically designed to probe the lower solar atmosphere within the wavelength range of 350–700 nm. A plethora of photospheric and chromospheric lines are available, including the Ca i 423 nm and Ca ii 393 nm spectral lines, so is optimally designed to contribute to the future goals of lower atmospheric wave studies through high-cadence polarimetry across multiple ionization states to help with subsequent spectropolarimetric inversions. francis utilizes 400 optical fibers in a $20\times 20$ configuration that are each only 40 $\upmu$m in diameter. By focusing on the blue portion of the electromagnetic spectrum, the cladding thickness can be significantly reduced over its infrared counterparts (e.g., DL-NIRSP), giving each fiber a total overall diameter of 55 $\upmu$m (40 $\upmu$m diameter core plus 7.5 $\upmu$m cladding). When arranged in an offset pattern to maximize the filling factor of the fiber inlet aperture, this results in the $20\times 20$ fiber array occupying a $1.10\times 1.13$ mm$^{2}$ surface area on the entrance ferrule, which can be seen in the upper-left panel of Fig. 67. Each of the 400 fibers are mapped on to a one-dimensional linear array that is approximately 22 mm in length (upper-right panel of Fig. 67), before being passed into the spectrograph entrance slit. The lower panels of Fig. 67 show a sample active region captured in the blue continuum at 390 nm (left) and in the core of the Ca ii K spectral line at 393 nm (right), with the yellow circles depicting an example sampling provided by the $20\times 20$ input fiber array. Here, each fiber covers an approximate $2''$ diameter on the solar surface, providing a total field of view spanning $40\times 40$ arcsec$^{2}$. In this configuration, approximately 35 fibers cover the sunspot umbra, 170 fibers cover the penumbra, with the remaining 195 fibers sampling the surrounding quiet Sun. By adjusting the size of the imaged beam incident on the fiber entrance ferrule, it becomes possible to easily adjust the spatial field-of-view achievable with this instrument. Furthermore, since the inlet fiber ferrule is polished, it becomes possible to utilize a slitjaw camera to accurately coalign the resulting hyper-spectropolarimetric images with other instruments in operation at the same time. Installation and commissioning of francis prototype instrument took place at the DST during the summer of 2022, where it will remain as a common-user instrument until such times as the Indian NLST becomes close to operations.

Fig. 67

The inlet fiber ferrule (upper left) of the hyper-spectropolarimetric imager being developed for the Indian NLST facility currently planned for construction. The dark rectangle ($1.10\times 1.13$ mm$^{2}$) towards the center of the ferrule contains a $20\times 20$ array of optical fibers, which are mapped into a linear configuration (upper right) allowing the light to be passed into the spectrograph housing. The $20\times 20$ fiber grid is arranged in an alternating offset configuration to maximize the filling factor of the fibers in the inlet aperture, which can be visualized by the yellow circles in the lower two panels. The lower left and right panels show images of a solar active region captured in the blue continuum (390 nm) and core of the Ca ii K (393 nm) spectral line, respectively. The yellow circles depict an example spacing (approximately $2''$ across each fiber), providing two-dimensional spectropolarimetric information at rapid cadences across a relatively large field of view. Here, approximately 35, 170, and 195 fibers sample the sunspot umbra, penumbra, and surrounding quiet Sun, respectively. An image of the Earth is provided in the lower-right corner of each solar image to provide a sense of scale. Images courtesy of the Astrophysics Research Centre at Queen’s University Belfast

Of course, one cause for concern with any instrument attempting to acquire Stokes I/Q/U/V modulation states for the purposes of spectropolarimetric inversions is the introduction of time delays while the modulators re-orientate themselves in order to cycle between the various polarization states. Common hardware to perform this task include the use of liquid crystal variable retarders (LCVRs) and rotating waveplates, which can often operate at rates exceeding 10 Hz to try and mitigate any timing delays (e.g., Lites 1987; Terrier et al. 2010; Hanaoka 2012; Alvarez-Herrero et al. 2017; Anan et al. 2018; Antonucci et al. 2020). Even with the ability to progress rapidly through modulation states, the non-simultaneity of the resulting exposures means that the resulting Stokes I/Q/U/V spectra will not have identical noise values at each wavelength, making inversions challenging in the limit of low photon counts. One method to combat this involves the use of a charge-caching camera, which is able to utilize fast charge transfer in the detector’s silicon chip to allow the use of a modulator running at kHz rates. By operating at (or above) the atmospheric turbulence timescales, such technology allows residual seeing effects to be largely eliminated. As a result, all Stokes modulation states can be accumulated at kHz rates before readout is performed, thus ensuring each Stokes modulation state has identical noise properties. An early development in this technique revolved around the Zurich Imaging Polarimeters (ZIMPOL; Povel 1995), in which a CCD detector array was alternately divided into photosensitive rows and storage rows that were shielded from light by a mask. By shifting the accumulated photons into neighboring storage rows, the acquisition process could be repeated over many modulation cycles until sufficient charge had been accumulated, hence ensuring high signal-to-noise values. Building upon the earlier work of Povel (1995) and harnessing the more versatile and photosensitive back-illuminated CMOS technology, Keller (2004) discusses the Charge Caching CMOS Detector for Polarimetry (C$^{3}$Po) concept, whereby a typical multiplexed CMOS pixel would provide a well depth exceeding $6\times 10^{6}$ electrons, providing approximately 500,000 electrons per modulation state, helping to achieve polarimetric precision down to the $10^{-5}$ level (Povel 1995; Stenflo 2013).

By employing a kHz custom-made camera (using a frame-transfer, back-illuminated pn-type CCD sensor), Iglesias et al. (2016) introduced a novel Fast SpectroPolarimeter (FSP) based on a polarization modulator with ferroelectric liquid crystals. The FSP could result in high temporal- and spatial-resolution full-Stokes measurements with a high sensitivity in the visible wavelength range (i.e., 400–700 nm), providing almost 100% duty cycle. The $1024\times 1024$ detector was shown to reach a high frame rate of 400 fps and low readout noise on the order of 5 e$^-$ rms (i.e., considerably smaller than that of, e.g., CRISP and IBIS, by a factor of 4). This implies that the FSP can capture the full Stokes parameters in 10 ms, hence, resulting in a high signal-to-noise polarimetric measurements with a high-cadence of, e.g., 2 s, after image restoration. While the FPS, which can spatially resolve scattering polarization signals (Zeuner et al. 2018), is an optimal instrument for studying high frequency waves through the solar photosphere and chromosphere, it has not yet been trialed with suitable hardware (e.g., a Fabry–Pérot) to allow sequential tuning and recording across multiple spectral lines. The latter is necessary for accurately constraining the plasma parameters through the multi-line inversions codes discussed above.

Next-generation instruments onboard the next (third) flight of the Sunrise balloon-borne solar observatory have been designed for such multi-line spectropolarimetric observations. In particular, two grating-based spectropolarimeters (with slit scanning and context imaging with slit-jaw cameras) are being employed (Barthol et al. 2018; Suematsu et al. 2018). The Sunrise UV Spectropolarimeter and Imager (SUSI) will sample the rich near-UV wavelength range of $300-430$ nm, which includes thousands of photospheric and over 150 chromospheric lines, poorly observable from the ground. Simultaneously, the Sunrise Chromospheric Infrared spectroPolarimeter (SCIP) will explore two spectral windows in the near-infrared, within the 765–855 nm wavelength range, which contains many magnetically sensitive lines sampling various heights in the lower solar atmosphere. The compromise for such many (full-Stokes) lines measurements with the slit-scanning instruments is, however, a trade-off between size of the field of view and the cadence of observations, all of which are important for tracing wave dynamics in the 3D atmosphere. Finally, of particular interest is the Polarimetric and Helioseismic Imager (PHI; Solanki et al. 2020) instrument on Solar Orbiter, which examines the Zeeman and Doppler effects linked to the photospheric Fe i 617.3 nm spectral line. PHI employs two telescopes: the first is a full-disk view of the Sun designed to capture information across all phases of the orbit, while the second is a high-resolution telescope that will allow structures as small as 200 km ($\approx 0{\,}.{\!\!}{''}27$) on the surface of the Sun to be examined at closest perihelion. This capability has the potential to provide long-duration, seeing-free polarimetric observations of the photosphere, which is particularly exciting for the examination of waves in developing magnetic features, such as network bright points, pores, and sunspots.

Importantly, as can be seen from the text above, new types of instruments (e.g., IFUs) and acquisition techniques (e.g., charge caching cameras) are paving the way for higher precision spectropolarimetry with improved spatial resolutions and temporal cadences, i.e., converging towards the ‘trifecta’ once thought impossible. Over the coming years, we expect the instrumentation described above to play a pivotal role in investigations of waves and oscillations in the lower solar atmosphere by providing the solar physics community with high precision data sequences necessary to take advantage of the cutting-edge inversion software concurrently being developed. What is clearly required at this stage is a ‘Level 2’ data repository for observing sequences that have already been processed with a specific inversion scheme. This would be similar to what is commonly available from space-borne observatories, such as the Joint Science Operations Center⁴ (JSOC), which provides calibrated Level 1 observations, complete with Level 2 data products, such as disambiguated vector magnetograms from the HMI/SDO instrument (Martens et al. 2012; Couvidat et al. 2016). Providing a similar data repository scheme for current- and next-generation ground-based facilities would unquestionably boost the accessibility of these key scientific data products for the global solar physics community. We therefore look forward to additional proposals to increase the compute and storage capabilities of such a Level 2 data center, particularly if it can incorporate observations from other ground-based solar facilities, including the DST, SST, NST, GREGOR, ALMA, and once commissioned, the NLST and EST.

5 Conclusions

In this review, we have attempted to overview (see Sect. 2) the main wave analyses techniques currently harnessed by the global solar physics community. The reason for this is two-fold: (1) We would like to provide the next generation of researchers studying oscillatory phenomena a concise, yet helpful introduction to the techniques that underpin current research efforts, and (2) We wish to establish an analysis tool framework from which researchers can build upon. All too often analyses techniques are reinvented by researchers who do not have easy access to existing algorithms, which unfortunately delays time frames associated with the dissemination and publication of results. Hence, as part of this review, we wish to encourage researchers to visit the Waves in the Lower Solar Atmosphere⁵ (WaLSA) dedicated website repository, where a collection of the most readily used codes are available to download and employ. As a natural part of the feedback process, if researchers update and improve the existing wave analyses software, then their updated codes can be hosted on the WaLSA platform for others to avail of, with appropriate references provided to document the underlying improvements.

Sections 3 and 4 document a plethora of high-impact wave studies that have come to fruition over the last number of years. While these studies have unquestionably improved our understanding of the tenuous solar atmosphere, they naturally pose yet more unanswered questions. Thankfully, we are entering an era of discovery with the current and upcoming commissioning of high value (monetary and scientifically) ground-based and space-borne facilities such as DKIST, Solar Orbiter, and Solar-C. This places the community in an ideal position to explore new high-resolution observations with unrivaled accuracy. Many of the instruments associated with such observing facilities are well suited to wave studies, with multi-wavelength capabilities, rapid cadences, high polarimetric precisions, and unprecedented spatial resolutions planned from the very beginning.

In addition to the imminent next-generation observational facilities we will have at our disposal, researchers will also be able to capitalize on continually updated high performance computing (HPC) infrastructures to model and simulate the observed wave signatures with unprecedented resolution. Globally, the Distributed Research utilising Advanced Computing⁶ (DiRAC), ARCHER,⁷ the Norwegian academic e-infrastructure,⁸ La Red Española de Supercomputación⁹ (RES, Spanish Supercomputing Network), and the NASA Pleiades¹⁰ HPC supercomputing systems (to name but a few examples) provide tens of petaflops of compute capability. Such evolving HPC facilities are crucial for the accurate replication of physics, particularly down to the spatial and temporal scales imminently visible by the newest observing facilities.

As a community, we therefore look forward in anticipation to the new era of understanding that will be brought to fruition by the newest researchers, observatories, and computing facilities we will have at our disposal over the decades to come.