# High Frequency Waves in Chromospheric Spicules

W. Bate Astrophysics Research Centre, School of Mathematics and Physics,
Queen’s University Belfast, Belfast, BT7 1NN, UK  D. B. Jess Astrophysics
Research Centre, School of Mathematics and Physics, Queen’s University
Belfast, Belfast, BT7 1NN, UK  Department of Physics and Astronomy, California
State University Northridge, 18111 Nordhoff Street, Northridge, CA 91330, USA
V. M. Nakariakov Centre for Fusion, Space and Astrophysics, Physics
Department, University of Warwick, Coventry, CV4 7AL, UK  St. Petersburg
Branch, Special Astrophysical Observatory, Russian Academy of Sciences,
196140, St. Petersburg, Russia  S. D. T. Grant Astrophysics Research Centre,
School of Mathematics and Physics, Queen’s University Belfast, Belfast, BT7
1NN, UK  S. Jafarzadeh Leibniz Institute for Solar Physics (KIS), Schöneckstr.
6, 79104 Freiburg, Germany  Rosseland Centre for Solar Physics, University of
Oslo, P.O. Box 1029 Blindern, NO-0315 Oslo, Norway  M. Stangalini ASI, Italian
Space Agency, Via del Politecnico snc, 00133, Rome, Italy  P. H. Keys
Astrophysics Research Centre, School of Mathematics and Physics, Queen’s
University Belfast, Belfast, BT7 1NN, UK  D. J. Christian Department of
Physics and Astronomy, California State University Northridge, 18111 Nordhoff
Street, Northridge, CA 91330, USA  F. P. Keenan Astrophysics Research Centre,
School of Mathematics and Physics, Queen’s University Belfast, Belfast, BT7
1NN, UK

###### Abstract

Using high cadence observations from the Hydrogen-alpha Rapid Dynamics camera
imaging system on the Dunn Solar Telescope, we present an investigation of the
statistical properties of transverse oscillations in spicules captured above
the solar limb. At five equally separated atmospheric heights, spanning
approximately $4900-7500$ km, we have detected a total of $15{\,}959$
individual wave events, with a mean displacement amplitude of $151\pm 124$ km,
a mean period of $54\pm 45$ s, and a mean projected velocity amplitude of
$21\pm 13$ km s-1. We find that both the displacement and velocity amplitudes
increase with height above the solar limb, ranging from $132\pm 111$ km and
$17.7\pm 10.6$ km s-1 at $\approx 4900$ km, and $168\pm 125$ km and $26.3\pm
14.1$ km s-1 at $\approx 7500$ km, respectively. Following the examination of
neighboring oscillations in time and space, we find 45% of the waves to be
upwardly propagating, 49% to be downwardly propagating, and 6% to be standing,
with mean absolute phase velocities for the propagating waves on the order of
$75-150$ km s-1. While the energy flux of the waves propagating downwards does
not appear to depend on height, we find the energy flux of the upwardly
propagating waves decreases with atmospheric height at a rate of
$-13{\,}200\pm 6500$ W m-2/Mm. As a result, this decrease in energy flux as
the waves propagate upwards may provide significant thermal input into the
local plasma.

Sun: atmosphere — Sun: chromosphere — Sun: oscillations — Sun: solar spicules

††facilities: DST (HARDcam; Jess et al., 2012a)††software: MGN (Morgan &
Druckmüller, 2014), Auto-NUWT (Morton et al., 2013; Weberg et al., 2018)

##  1 Introduction

Spicules are dynamic plasma jets that are prevalent within the solar
chromosphere, and which generally have diameters on the order of hundreds of
km. They are relatively short-lived features, typically having a lifetime of
less than $10$ minutes (Pereira et al., 2012). When viewed in the visible and
UV bands at the solar limb, spicules appear ubiquitously as a dense forest of
narrow, straw-like features (Sterling, 1998, 2000).

Figure 1: Contextual SDO/HMI continuum (left), ROSA G-band (middle), and
H$\alpha$ line core (right) images acquired at 13:59:09 UT. The area imaged by
ROSA and HARDcam is marked by the red square in the full disk SDO/HMI
continuum image. Numerous spicules are clearly visible above the solar limb as
narrow, straw-like structures in the corresponding H$\alpha$ image.

Secchi (1877) was the first to observe solar spicules, and they have remained
a focal point of solar physics research over the last 140 years. Transverse
oscillations in spicules were first identified in the $1960$s (Pasachoff et
al., 1968), utilizing ground-based observations obtained at the Sacramento
Peak Observatory (for a comprehensive review see Zaqarashvili & Erdélyi,
2009). The magnetic cylinder model is generally accepted as being the most
applicable to spicules, allowing their oscillatory behavior to be interpreted
and modeled as magnetohydrodynamic (MHD; Alfvén, 1942) modes (Edwin & Roberts,
1983). Sterling (2000) highlighted that high-resolution observations, due to
the small width of the structures, are vital for a complete description of the
spicule wave phenomena. Wedemeyer-Böhm et al. (2007) also note that the
ability to detect oscillatory power at higher frequencies is influenced by the
spatial resolution of the observations (see also the discussions provided by
Jess et al., 2020, 2021).

One of the major focuses of current solar physics research is the so-called
‘coronal heating paradox’. One of the proposed theoretical mechanisms to
explain the source of this heating is linked to the propagation and
dissipation of wave phenomena, commonly referred to as the ‘AC’ heating
mechanism (Alfvén & Lindblad, 1947). Spicules are of particular interest when
attempting to explain the heating of the solar atmosphere due to their
potential to facilitate the transfer of mass and energy between the
photosphere and corona. They are often categorized by their properties into
two types, type i and type ii (De Pontieu et al., 2007a; Pereira et al.,
2012), although such distinct classifications are still under debate (e.g.,
Zhang et al., 2012; Pereira et al., 2013).

Observations of transverse oscillations of spicules, fibrils, and mottles in
the upper chromosphere typically find mean periods on the order of $80-300$ s,
often with $\sim 50-1000$ examples found during the course of an individual
data sequence (Nikol’Skii & Sazanov, 1967; Nikolsky & Platova, 1971; De
Pontieu et al., 2007b; Kuridze et al., 2012; Morton et al., 2012, 2013, 2014).
A major exception are the short period ($45\pm 30$ s) transverse oscillations
found in $89$ type ii spicules within a coronal hole identified by Okamoto &
De Pontieu (2011). However, Okamoto & De Pontieu (2011) suggest that this
short average period is likely due to them utilizing a methodology which did
not allow for the measurement of the properties of the longer period waves
($>100$ s). Another exception is the more recent work of Shetye et al. (2021),
which found 30 examples of transverse spicule oscillations with periods
ranging from $16-100$ s. These authors also note a selection effect due to
only choosing an event if its oscillation period is less than its lifetime,
resulting in longer period waves not being considered.

Finding the energy flux of spicule oscillations is an important step in
investigating their contribution to the heating of the chromosphere and
corona. It is estimated that an energy flux of $10^{3}-10^{4}$ W m-2 is
required to heat the chromosphere. The energy required to heat the corona is
around an order of magnitude less than that required to heat the chromosphere
(Withbroe & Noyes, 1977). This suggests that accounting for chromospheric
heating is a challenge of equal or greater magnitude than for coronal heating
when investigating solar atmospheric heating mechanisms.

Energy flux estimations are based on the interpretation of these transverse
oscillations as MHD wave modes. De Pontieu et al. (2007b) interpreted the
transverse oscillations of spicules as bulk Alfvén waves and assumed a filling
factor of unity. Using this interpretation, an energy flux estimate of
$4000-7000$ W m-2 was suggested. However, this bulk Alfvén interpretation has
attracted criticism, with Erdélyi & Fedun (2007) and Van Doorsselaere et al.
(2008) pointing out that Alfvén waves do not result in the bulk transverse
motions observed, and instead proposing that these transverse oscillations are
best interpreted as kink modes. The filling factor, a measure of what fraction
of the total volume is occupied by oscillating spicules, is another extremely
important consideration for energy flux calculations (Van Doorsselaere et al.,
2014). An equivalent interpretation, assuming that spicules have an
approximately constant width across varying heights, would be the ratio of the
area of the solar surface covered by spicules to the total solar surface area.

Makita (2003) found a spicule filling factor of $5\%$ at a height of $4000$ km
using Ca ii H & K line observations taken during a solar eclipse. This
suggests that a filling factor of $0.05$ is more appropriate than unity. Using
the revised interpretations of the most realistic MHD mode and associated
filling factor, Van Doorsselaere et al. (2014) found that the energy flux
estimates by De Pontieu et al. (2007b) were reduced from $4000-7000$ W m-2 to
$200-700$ W m-2, a difference exceeding one order of magnitude. Furthermore,
Morton et al. (2012) used high-resolution H$\alpha$ observations taken by the
Dunn Solar Telescope to find a similar upper limit for the filling factor
($4-5$%) for open chromospheric structures that connected to higher layers of
the solar atmosphere. By interpreting the transverse oscillations of fibrils
as kink modes, the authors estimated the energy flux as $170\pm 110$ W m-2,
similar to that derived by Van Doorsselaere et al. (2014). However, in
addition to the $4-5$% filling factors commonly used in modern literature,
lower estimates have also been put forward, with Beckers (1972) suggesting a
filling factor of $0.6\%$. As a result, it is generally believed that the
spicule filling factor spans an approximate order-of-magnitude (ranging
between $\approx 0.5-5$%), with differing values being applicable depending on
factors such as the atmospheric height sampled and the degree of solar
activity (i.e., it is not a quantity that can be applied universally across
all observations). The influence of the chosen filling factor on energy flux
calculations is discussed further in Section 3.

An important caveat when interpreting these energy flux estimates is that they
are only based on resolved transverse oscillations. Waves with amplitudes too
small to be spatially resolved or periods too short to be temporally resolved
are not included in these estimations, leaving the possibility that a
significant amount of wave energy may be unaccounted for (Verth & Jess, 2016).
Another aspect contributing to the underestimation of the total energy flux
may be the presence of kink motions along the observer’s line-of-sight, which
will not manifest as visible transverse oscillations. Examples of this have
been documented by Sharma et al. (2018) and Shetye et al. (2021), who measured
helical motions of spicules through Doppler measurements (see also the
modeling work by Zaqarashvili & Skhirtladze, 2008). If these line-of-sight
motions are not taken into account when calculating the energy flux, it may
result in an underestimation of the true value.

The aim of the current study is to identify the properties of spicule
oscillations across a statistically significant sample that is extracted from
different chromospheric heights. With oscillation characteristics measured
across a range of atmospheric layers, we calculate the energy flux carried by
these waves as a function of geometric height. To achieve this objective, we
utilize ground-based instrumentation with high spatial and temporal
resolutions, providing unprecedented data products that are ideally suited for
this study.

##  2 Observations

Our analysis employs data collected on 2015 July 27 from 13:52 – 15:29 UT
using the Dunn Solar Telescope (DST; Dunn, 1969) at the National Solar
Observatory in New Mexico, USA. The Rapid Oscillations in the Solar Atmosphere
(ROSA; Jess et al., 2010b) and Hydrogen-alpha Rapid Dynamics camera (HARDcam;
Jess et al., 2012a) imaging systems were used to observe a large sunspot,
which was part of NOAA AR12391, close to the solar limb at N07.8E73.6 in the
conventional heliographic coordinate system. Seeing conditions remained
excellent throughout the first hour of the observing period, gradually
worsening towards the latter stages of the observing window.

HARDcam observations employed a narrow ($0.25$Å FWHM) bandpass filter centered
on the H$\alpha$ line core ($6562.8$Å), while the ROSA camera system observed
the same region through G-band (10Å FWHM centered at 4305Å) and broadband
4170Å continuum filters. The HARDcam data have a pixel size of
$0.092^{\prime\prime}$ ($66.5$ km), providing a $180^{\prime\prime}\times
180^{\prime\prime}$ field-of-view, while the ROSA system was slightly
undersampled ($0.180^{\prime\prime}$ per pixel) to provide an identical field-
of-view size to that of the HARDcam observations. To correct for wavefront
deformations in real time, higher-order adaptive optics (AO) were used during
the observations (Rimmele, 2004; Rimmele & Marino, 2011). Original data from
both ROSA and HARDcam were taken at a frame rate of $30.3$ s-1, with the
images synchronized by way of a master trigger with microsecond precision
(Jess et al., 2010b). The resulting HARDcam H$\alpha$ images were then
improved using speckle reconstruction algorithms (Wöger et al., 2008),
utilizing a $30\rightarrow 1$ restoration, resulting in a final reconstructed
cadence of $0.990$ s.

Figure 2: An H$\alpha$ core sub-field ($67\times 16$ Mm2) image acquired using
HARDcam at 14:49:45 UT. Numerous spicules are clearly visible above the solar
limb as narrow, straw-like structures. The two most extreme slits used to take
the time distance diagrams are shown by the white lines, at heights of $4890$
and $7500$ km. The axes are shown using different scales to aid with visual
clarity.

ROSA continuum observations were coaligned using cross-correlation techniques
(see, e.g., Jess et al., 2010a) with contemporaneous continuum images from the
Helioseismic and Magnetic Imager (HMI; Schou et al., 2012) on board the Solar
Dynamics Observatory (SDO; Pesnell et al., 2012), providing sub-arcsecond
pointing accuracy for the field-of-view covered by the DST. Following this,
the HARDcam field was aligned with the master ROSA images using sequences of
targets acquired during the calibration procedures at the DST, resulting in
H$\alpha$ observations that have precise pointing metadata that is consistent
with modern space-based observatories. Contextual images from SDO/HMI, ROSA,
and HARDcam, following the processing steps outlined above, are shown in
Figure 1.

##  3 Analysis and Discussion

During the course of the observations, the DST’s AO system was locked onto the
high-contrast sunspot structure that was very close to the limb. As a result,
limb spicules close to the central portion of the field-of-view were
accurately corrected from atmospheric seeing effects by the AO. Hence, the
current HARDcam H$\alpha$ dataset offers an unprecedented opportunity to
examine limb spicules at extremely high time cadence (0.990 s) and spatial
resolution (133 km two-pixel resolution).

Figure 3: A time-distance diagram captured using a curved slit at a height of
$6850$ km above the solar limb. Each bright streak is a feature passing
through the slit, with the clear oscillatory features representative of
transverse motions displaying a range of amplitudes and periods.

A sub-field, spanning approximately 70 Mm along the central portion of the
field-of-view, where the AO corrections were operating optimally, was isolated
for further study. As the DST was tracking the sunspot contained within the
field-of-view, over time the pixel coordinates corresponding to the limb
position change as a result of the sunspot rotating on to the disk. The image
sequence was hence stabilized with respect to the limb, which was achieved by
first choosing a reference frame towards the beginning of the dataset. Next, a
small area of the limb image with high contrast was selected, with subsequent
images compared and shifted using two-dimensional cross correlation
techniques. Pixel shift values that produced the highest cross correlation
coefficients were selected and applied to each image in the time series
iteratively. The resulting shifted images lead to the limb remaining
stationary at the same pixel location throughout the dataset, providing a
robust baseline from which to examine spicule oscillations above the fixed
limb.

Figure 4: Histograms of the wave properties identified at a height of $6850$
km above the solar limb. The upper, middle, and lower panels display
information related to the displacement amplitudes, periods, and velocity
amplitudes, respectively. Measurements of the corresponding averages and
deviations are displayed on the right of each histogram. MAD denotes the
median absolute deviation.

Multiscale Gaussian Normalization (MGN; Morgan & Druckmüller, 2014) was
applied to each image in the dataset in order to more easily identify each
spicule and its associated motion. It must be noted that MGN does not preserve
photometric accuracy. However, this is not an issue when mapping the
transverse oscillations of features since we are not concerned with
comparisons of relative intensities. For the application of MGN, we employed
the convolution of HARDcam images with Gaussian kernels with one-sigma widths
of $w=1.25,2.5,5,10,20,40$ pixels, followed by the production of gamma-
transformed images with a $\gamma$ value of $3.2$ (Poynton, 2003).

Five slits were placed at equally spaced, constant radial heights above the
limb, spanning approximately $4900$ km to $7500$ km in steps of $\approx 650$
km. These slits were curved in nature in order to maintain a constant radial
height above the limb, and the highest and lowest slits are shown by the white
lines in Figure 2. When taking this approach, it is important to note that
superposition along the line of sight of spicules anchored behind the limb, in
front of the limb, and on the limb is unavoidable. As a result of the slit
heights being based on a geometric distance above the limb, this will result
in the foreground/background spicules being sampled further along their
lengths than those precisely located on the limb. We have carefully selected
the minimum and maximum heights of the slits to be in the range of $4890-7500$
km (see Figure 2), which is towards the upper end of the ‘dense forest’ of
spicules, hence minimizing the degree of feature superposition. Due to the
(minimized) spicule superposition affecting each of the slits in a similar
way, and considering the large numbers of spicules observed at each height,
comparisons between wave properties taken with different slits will still be
valid. However, it is still important to consider this effect when examining
wave properties taken from a single slit in isolation, since the chosen slit
height will be a minimum value of the distance sampled along the spicule due
to these geometric considerations. Time-distance diagrams were then produced
from each of the slits, with an example shown in Figure 3.

The Automatic Northumbria University Wave Tracking (Auto-NUWT; Morton et al.,
2013; Weberg et al., 2018) code was utilized in order to identify the location
of the spicules as a function of time, track their transverse motion, and
extract the properties of their oscillations. Features are identified by
fitting a sum of Gaussians to each time slice in the time-distance diagrams,
enabling the determination of sub-pixel values for the location of the center
of the feature. The transverse oscillatory behavior of these features is
probed through the application of Fourier analysis to the position of the
center of the feature as a function of time.

Figure 5: The displacement curves, corresponding to two of the waves
identified in our dataset, are shown (in their raw form) using the solid black
lines. The dashed red lines highlight the fitted waves with properties derived
using Fourier analysis. The top panel shows a wave with a period of $\approx
138$ s and a displacement amplitude of $\approx 358$ km, while the bottom
panel shows a wave with a period of $\approx 20$ s and a displacement
amplitude of $\approx 79$ km.

At each of the five heights considered, over $3000$ spicule features are
detected in the time-distance diagrams. Employing Fourier analysis, the
properties of the waves present in the transverse motions of these features
were determined. As a representative example, the averages and deviations of
wave properties found at a height of $6850$ km above the limb are displayed in
Figure 4, where the distributions of the displacement amplitudes, periods, and
calculated velocity amplitudes of these waves are plotted as histograms. These
properties follow approximate log-normal distributions, which are shown by the
solid green lines in Figure 4. Log-normal distributions for these properties
are consistent with those found in previous studies (e.g., De Pontieu et al.,
2007b; Okamoto & De Pontieu, 2011; Pereira et al., 2012).

Figure 6: Mean values of spicule displacement amplitudes (top panel) and velocity amplitudes (lower panel) plotted against height above the solar limb. Linear lines of best fit are shown in both panels using a dashed red line. In the bottom panel a quadratic fit is shown using a dashed blue line. Errors associated with each data point represent the standard error on the mean. Height | Number of Waves | Displacement Amplitude | Period | Velocity Amplitude
---|---|---|---|---
(km) |  | (km) | (s) | (km s-1)
4890 | 4880 | 132.4$\pm 111.2$ | 55.1$\pm 45.0$ | 17.7$\pm 10.6$
5550 | 4920 | 139.6$\pm 118.8$ | 57.0$\pm 48.1$ | 18.3$\pm 10.7$
6200 | 5022 | 152.4$\pm 128.5$ | 55.8$\pm 48.1$ | 20.6$\pm 11.6$
6850 | 5209 | 163.7$\pm 128.6$ | 53.4$\pm 45.4$ | 23.4$\pm 13.1$
7500 | 5298 | 168.1$\pm 125.4$ | 48.9$\pm 39.9$ | 26.3$\pm 14.1$
Table 1: Mean wave properties and their standard deviations at each sampled
height.

Average wave properties for each height are presented in Table 1. Importantly,
the averages of the displacement and velocity amplitudes appear to be
consistent with those found in previous studies of transverse waves in
spicules. However, the average period of the waves in the current study are
shifted to lower values than those found previously (see the summary provided
by Jess et al., 2015). Specifically, the majority of earlier studies found
average periods on the order of $80-300$ s (e.g., Nikolsky & Platova, 1971; De
Pontieu et al., 2007b; Kuridze et al., 2012; Morton et al., 2012, 2013, 2014),
while we find the average period to be $53\pm 45$ s (middle panel of Figure
4). It should be noted that the mean has been chosen for comparison with
previous studies here. However, as the wave properties approximately follow
log-normal distributions, the modal value represents a more useful statistic
in understanding the peak of this distribution. As a point of comparison,
Nikolsky & Platova (1971) found modal and mean periods of $60$ s and $85$ s,
respectively, for their observed spicule oscillations, whereas example modal
and mean periods found at a height of $6850$ km above the solar limb in this
study are $\approx$24 s and $\approx$54 s, respectively, as shown in Figure 4.
We consider the detection of these shorter period waves likely due, at least
in part, to the unprecedented $\sim 1$ s time cadence of the dataset utilized.
For comparison, previous investigations using data from the Swedish Solar
Telescope provided cadences on the order of $5$ s, which would make it very
difficult for the lowest oscillation periods ($<10$ s) identified here to be
detected.

Across all five defined slits, over $16{\,}600$ spicular threads were
identified, of which $15{\,}959$ (95.9%) exhibit at least one complete wave
cycle. Of these examples, $8568$ (51.5%) threads exhibit a single wave, $5770$
(34.7%) consist of two superposed waves, and $1621$ (9.2%) have three (or
more) superposed waves. These proportions are similar to those found in
transverse oscillations of coronal plumes using Auto-NUWT by Weberg et al.
(2018).

Two examples of the identified waves are shown in Figure 5. These are chosen
as they have radically different periods and displacement amplitudes,
consisting of $\approx$138 s and $\approx$358 km, respectively, for the top
panel, with $\approx$20 s and $\approx$79 km, respectively, for the lower
panel. Both waves are observed for longer than one full period. The wave
identified in the top panel of Figure 5 has properties consistent with those
found in previous studies of transverse spicule oscillations (see the review
by Jess et al., 2015), highlighting that these longer period ($>50$ s) waves
are also present within our data and are fitted well using our techniques.
However, due to the high spatial and temporal resolutions provided by HARDcam,
much shorter period waves are able to be identified, including the example
shown in the lower panel of Figure 5.

As the wave properties have been determined for each of the five equally-
spaced slits above the solar limb, we are able to compare and study
characteristics as a function of atmospheric height. The mean values for
displacement amplitude and velocity amplitude are shown in Figure 6, where
both parameters can be seen to increase with height. By fitting a linear line
of best fit through the corresponding data points (see the dashed red lines in
Figure 6), the displacement amplitude increases at a rate of $14.6\pm 0.8$
km$/$Mm, and the velocity amplitude at $3.33\pm 0.08$ km s${}^{-1}/$Mm. The
conservation of energy flux requires a reciprocal relation between density and
velocity amplitude (see, e.g., the discussions in Stein & Leibacher, 1974;
Ebadi et al., 2012; de la Cruz Rodríguez et al., 2013; Khomenko & Collados,
2015; Grant et al., 2018; Houston et al., 2018, 2020; Riedl et al., 2021).
Employing spectropolarimetric inversions of the Ca ii spectral line, Kuridze
et al. (2021) revealed evidence that the mass density of spicules decreases
exponentially with height, requiring the velocity amplitude to similarly
increase to conserve energy flux. Hence, a quadratic fit is presented in the
lower panel of Figure 6 using a blue dashed line to show the potential synergy
between expected mass density and velocity amplitude. However, we note that
this is presented only for completeness, since it is difficult to infer the
true nature of the relationship with only five data points. The average
periods do not show a similar trend with atmospheric height, instead ranging
within the same interval of $48.8-57.0$ s for the five heights considered.

In order to measure the phase velocity of these waves, it was necessary to
identify the same feature across different heights. This was achieved by
extracting individual wave properties from a certain atmospheric height and
searching through the wave catalog for waves with similar properties
identified at an adjacent height. The properties considered for this study
were the equilibrium $x$-position of the spicule ($\pm$5 pixels or $\pm$330
km), the midpoint time (whether-or-not the next spicule feature lay between
the start and end times of the wave being compared), duration of the
oscillation ($\pm 50\%$), and the frequency ($\pm 10\%$). Based on these
criteria, around $140$ waves were found to be suitably similar between each
set of adjacent heights, providing large number statistics with a similar
proportion (tracked waves in relation to total identified waves) to that
documented by Jafarzadeh et al. (2017).

Figure 7: A histogram showing the calculated phase velocities of the $135$
propagating waves identified between the heights of $6200$ and $6850$ km above
the solar limb. The solid red line denotes the cumulative probability
function. A bin width of $20$ km s-1 was used for the creation of this
histogram.

The phase difference between all sets of waves identified at adjacent heights
was calculated using Fourier phase lag analysis. The cross-power spectrum was
calculated using the representative Fourier spectra of the two waves found at
adjacent heights (Bendat & Piersol, 2000). The real part of the cross-power
spectrum (co-spectrum) was used to verify that each original Fourier spectrum
had a peak at the same frequency. The phase of the cross-power spectrum was
then computed at the same frequency to determine the phase lag between the two
heights (Vaughan & Nowak, 1997). Finding this phase lag, $\phi$ (in degrees),
allows for the calculation of the phase velocity, $v_{ph}$ (in km$/$s),

$v_{ph}=\frac{360d}{\mathcal{T}\phi}\ ,$ (1)

where $d$ is the height difference between the two slits in km, and
$\mathcal{T}$ is the period of the wave in seconds (Jess et al., 2012b;
González Manrique et al., 2020).

Importantly, calculation of the phase velocities of the waves embedded within
the spicules allows for the eventual calculation of their energy fluxes. The
distribution of phase velocities for the $135$ propagating waves identified
traveling between the heights of $6200\rightarrow 6850$ km are shown in Figure
7, where waves propagating in both the upward and downward directions are
identified. Any waves displaying zero phase lags (i.e., providing infinite
phase speeds in Equation 1) were classified as standing modes. Due to the
relatively small number of standing modes present in our dataset, this type of
wave is not included in the histogram depicted in Figure 7. We must highlight
that the distribution of upwardly and downwardly propagating waves shown in
Figure 7 appear to originate from the same population, with an approximately
Gaussian distribution encompassing waves propagating both upwards and
downwards. However, when examining the energy flux carried by these waves it
is important to examine the direction of energy propagation, which is
determined by the sign of their associated phase speed. Thus, distinctions are
made between upwardly and downwardly propagating waves for the sake of further
analysis, but it must be emphasized that there do not seem to be two distinct
populations present in Figure 7.

It should be noted that any plasma flows within the spicules will affect the
apparent phase velocities of the measured kink oscillations. In the case of
upflowing plasma, the apparent (i.e., measured) phase velocities of the
upwardly propagating waves will be related to $(v_{ph}+U)$, while the
downwardly propagating waves will have apparent phase velocities equal to
$(-v_{ph}+U)$, where $v_{ph}$ is the true phase velocity and $U$ is the
velocity of the upflow (Nakariakov & Roberts, 1995). This is similar to
observations put forward by Grant et al. (2015), who examined the bulk plasma
upflow within a magnetic pore and the subsequent effect this had on the
apparent wave speeds of sausage mode oscillations. Strong upflows are
typically associated with type ii spicules. However, the spicules observed in
this study are likely not best characterized by this classification (De
Pontieu et al., 2007a). This effectively means that the apparent phase
velocities of the upwardly propagating waves can be considered an upper limit
to their true phase velocities. Conversely, in the case of the downwardly
propagating waves, this can be considered as a lower limit. As the velocity of
any possible upflows are not known, the measured phase velocities have been
used in all further calculations, but it is important to note that this will
result in the calculated energy fluxes being upper/lower limits for the
upwardly/downwardly propagating kink waves.

Across the four sets of adjacent heights, the occurrence rates of upwardly
propagating, downwardly propagating, and standing mode waves were found to be
$45\%$, $49\%$, and $6\%$, respectively. This is in contrast to the occurrence
rates found by Okamoto & De Pontieu (2011) of $59\%$, $21\%$, and $20\%$,
respectively. However, the spicules examined by Okamoto & De Pontieu (2011)
were observed within a coronal hole, so may have different properties to those
examined here. Importantly, our present study highlights a more equal balance
of upward/downward propagation, with fewer examples categorized as standing
modes. The lack of standing mode detections may also be a consequence of the
incredibly high spatial and temporal resolutions of the HARDcam dataset, since
phase precision is drastically improved as a result of the sub-1 second
cadence.

It might initially be assumed that a roughly equal balance of upward/downward
propagation should be expected, due to the high reflection coefficient of the
transition region (Hollweg et al., 1982). Liu et al. (2014) also observed
downwardly propagating transverse waves within solar spicules and note that
low-frequency (periods of $\approx 100$ s) are expected to reflect strongly in
the transition region (Suzuki & Inutsuka, 2005). However, Okamoto & De Pontieu
(2011) argue that such reflection would result in more standing modes being
observed due to the superposition of upwardly and downwardly propagating
waves. This would create an imbalance in observations, with more upward than
downward propagations detected. This superposition is, however, heavily
dependent on the height of the reflecting boundary, the phase velocity of the
upward wave, the lifetime of the spicule, and the time that the wave persists
for. If there is insufficient time for the reflected wave to interact with the
upward wave, due to any combination of the aforementioned criteria, then wave
superposition (and hence a standing wave) will not be observed. Although a
full characterization of the driving mechanisms behind the downwardly
propagating waves, as well as the clear domination of these waves with a phase
speed around zero (see Figure 7), is beyond the scope of the present work,
these are important questions to be investigated in future studies.

Height | Upward | Downward
---|---|---
| Phase Velocity | Phase Velocity
(km) | (km s-1) | (km s-1)
$4890\rightarrow 5550$ | $128\pm 23$ | $75\pm 12$
$5550\rightarrow 6200$ | $131\pm 23$ | $82\pm 23$
$6200\rightarrow 6850$ | $139\pm 25$ | $101\pm 15$
$6850\rightarrow 7500$ | $147\pm 23$ | $128\pm 23$
Table 2: Mean phase velocities for each set of adjacent heights that are
defined in Table 1.

For each of the four sets of adjacent heights, waves identified as upwardly
propagating were segregated from their downwardly propagating counterparts. It
was hence possible to calculate the rate of change of phase velocity as a
function of atmospheric height for both the upwardly and downwardly
propagating waves. Average phase velocities for each set of adjacent heights
are shown in Table 2 and plotted in Figure 8, where the upper panel
corresponds to the average phase velocity of the upwardly propagating waves,
while the lower panel depicts the average phase velocity of the downward
propagation. The uncertainties shown in Figure 8 have been calculated
following the ‘bootstrap’ methodologies described by Efron et al. (1979). Due
to the combined presence of traditional (periods $\geq 50$ s) and high-
frequency (periods $<50$ s) spicule oscillations, it is challenging to assign
basic standard errors to the derived phase velocities, especially since the
equivalence (or lack thereof) of the driving mechanisms responsible for these
characteristics have yet to be observationally and/or theoretically verified.
As such, we apply bootstrapping techniques to better constrain the confidence
intervals of data following non-normal or unknown distributions (similar to
that presented by Simpson & Mayer-Hasselwander, 1986; Desmars et al., 2009;
Yao et al., 2017).

The upward phase velocities appear to increase with atmospheric height at a
rate of approximately ${10\pm 15}$ km s${}^{-1}/$Mm. However, due to the size
of the associated uncertainties (see the error bars in the upper panel of
Figure 8), it is difficult to unequivocally stipulate the precise
relationship. A more pronounced trend is present in the downward phase
velocities (lower panel of Figure 8), which appear to decrease (as the height
sampled decreases) at a rate of approximately $24\pm 11$ km s${}^{-1}/$Mm,
implying that the wave slows as it travels down the spicule and encounters
more dense layers of the lower solar atmosphere.

Figure 8: Mean values of phase velocity shown with height above the solar
limb. The values for upwardly and downwardly propagating waves are shown in
the top and bottom graphs respectively. Errors are calculated using
bootstrapping.

With the velocity amplitudes and phase velocities of the oscillations
measured, it was possible to estimate the energy flux associated with both the
upwardly and downwardly propagating waves. In order to calculate the energy
flux, a model for the density of the spicules with height is required. Kuridze
et al. (2021) observed a limb spicule and derived a model of its density using
the Non-LTE Inversion Code using the Lorien Engine (NICOLE; Socas-Navarro et
al., 2015) inversion code. The final density model takes the form,

$\rho(y)=\rho_{0}e^{(y-h_{0})/\Lambda}\ ,$ (2)

where $y$ is the height above the solar limb, $\rho(y)$ is the spicule density
as a function of height, $h_{0}$ is the base height, $\rho_{0}$ is the density
at the base height, and $\Lambda$ is the density scale height. Values for our
energy flux calculations were taken directly from Kuridze et al. (2021), where
$\rho_{0}\approx 6\times 10^{-7}$ kg m-3, $h_{0}=2000$ km, and $\Lambda=1500$
km.

The energy flux, $F$, from transverse waves in a multiple flux tube system can
be calculated as,

$F\approx f\frac{1}{2}(\rho_{i}+\rho_{e})v^{2}v_{gr}\ ,$ (3)

where $f$ is the density filling factor, $\rho_{i}$ is the density inside the
flux tube filled in by the spicule, $\rho_{e}$ is the density outside the
spicule, $v$ is the velocity amplitude, and $v_{gr}$ is the group speed (Van
Doorsselaere et al., 2014). For propagating kink waves, the group velocity can
be approximated by the phase speed, $v_{ph}$, as they are only weakly
dispersive (Terradas et al., 2010; Nakariakov et al., 2021). The internal
density for spicules can be assumed to be much larger than the external
density, i.e., $\rho_{i}\gg\rho_{e}$ (Uchida, 1961), providing a simplified
equation for the energy flux,

$F\approx f\frac{1}{2}\rho_{i}v^{2}v_{ph}\ .$ (4)

Taking the upper limit of the spicule density filling factor as $5\%$ (Morton
et al., 2012) allowed the energy fluxes to be calculated for each adjacent set
of heights, which are displayed individually for all propagating waves (top
panel) alongside upwardly (middle) and downwardly (bottom) propagating waves
in Figure 9. For all waves examined, it can clearly be seen that there is a
decrease in energy flux with height, indicated using solid black data points
in the upper panel of Figure 9. A linear line of best fit is presented using a
dashed black line in the upper panel of Figure 9, with a gradient of
$-12{\,}600$ W m${}^{-2}/$Mm. However, an exponential fit would perhaps be
more appropriate, since the main factor for the energy flux decrease is
expected to be density stratification, which is typically represented by a
decaying exponential profile with height (e.g., Verth et al., 2011). Due to
the relatively small number of data points under consideration, a linear fit
has been chosen for simplicity. Regardless of the fitting function employed,
the important message is that the energy flux of the propagating transverse
waves clearly decreases with atmospheric height, hinting at some sort of
damping and/or dissipation process.

It is important to consider the effect of using a filling factor of $5\%$.
This means that Equation 4 estimates the energy flux under the assumption that
the waves are omnipresent, i.e., does not take into account the sporadic
nature of the observed wave motion. In addition, as the waves are not seen to
exist in all spicules, the actual filling factor, $f$, should be reduced to
account for this effect. Thus, the estimation based on Equation 4 gives us the
upper limit of the energy flux in the waves. However, as the filling factor is
a multiplicative term, this only affects the magnitude of any energy flux
estimations. The trends in energy flux examined with respect to height are
independent of any adjustment to the filling factor. For example, using the
relatively low filling factor of $0.6\%$ suggested by Beckers (1972) will
simply lower all energy flux and rate of change of flux values by a linear
factor of $0.12$ when compared to those values calculated with a filling
factor of $5\%$. The values presented in the text and within Figure 9 utilize
a filling factor of $5\%$, unless stated otherwise, and should therefore be
taken as upper limits.

For all upwardly propagating waves, we observe the energy flux to decrease as
a function of height at a rate of $-13{\,}200\pm 6500$ W m${}^{-2}/$Mm, which
is indicated in the middle panel of Figure 9 using a dashed black line derived
from a linear least-squares fit. For completeness, it is estimated that energy
fluxes in the range of $10^{3}-10^{4}$ W m-2 are required to heat the
chromosphere (Withbroe & Noyes, 1977). Hence, the total energy flux, in
addition to the measured rate of energy flux decay with height, are on the
same order as the total energy input required to provide basal heating to the
solar chromosphere. Even considering the relatively low filling factor of
$0.6\%$, as suggested by Beckers (1972), the rate of energy flux decrease
would be $-1580\pm 780$ W m${}^{-2}/$Mm, still within the range that is needed
to balance the radiative losses of the chromosphere. By contrast, the energy
flux for all of the waves propagating downwards does not appear to depend on
the height sampled (black data points in the lower panel of Figure 9), with
the energy flux estimates remaining consistent (${\sim 4\times 10^{4}}$ W m-2)
across the height range of approximately $7500\rightarrow 4900$ km above the
solar limb. The similarity in the rate of energy flux drop off in height is
consistent between the full set of waves (upper panel of Figure 9) and the
upwardly propagating ones (middle panel of Figure 9). This is to be expected,
since the downward energy flux remains approximately constant with atmospheric
height.

Figure 9: Energy flux estimations as a function of atmospheric height for all
propagating waves (upper panel), upwardly propagating waves (middle panel),
and downwardly propagating waves (lower panel). The total energy flux provided
by short/long period waves is shown in black, while the energy fluxes for
short- ($<50$ s) and long-period ($>50$ s) waves are shown in red and blue,
respectively. The energy fluxes provided by the full set of waves (including
upwardly and downwardly propagating) and for all upwardly propagating waves
are depicted, using a linear line of best fit, as a dashed black line in the
upper and middle panels, with gradients equal to $-12{\,}600$ W m${}^{-2}/$Mm
and $-13{\,}200$ W m${}^{-2}/$Mm, respectively.

The decrease in upward energy flux with atmospheric height may be due to at
least three different factors: (1) physical thermalization of wave energy into
localized heat via dissipation mechanisms (e.g., Hollweg, 1986; He et al.,
2009; Antolin et al., 2015, 2018; Okamoto et al., 2015, to name but a few
examples), (2) damping of detectable transverse waves through the process of
mode conversion, where kink mode amplitudes decay as a result of the transfer
of energy from transverse kink oscillations to azimuthal Alfvén motions
(Pascoe et al., 2010, 2012, 2013), and/or (3) reflection of the waves downward
at varying heights above the solar limb (Hollweg et al., 1982; Suzuki &
Inutsuka, 2005). Tentative observational evidence has shown that torsional
Alfvén and kink waves may exist concurrently in spicules, providing credence
for the applicability of mode conversion processes (De Pontieu et al., 2012).
Previous modeling work by Sterling & Hollweg (1984) has shown that Alfvén
waves within spicules can produce high-frequency signatures, including
periodicities of $112$, $37$, and $22$ s for the fundamental, first, and
second harmonic resonant periods, respectively, which are similar to the
periodicities found in our current work. Employing simultaneous plane-of-sky
imaging and line-of-sight Doppler measurements will allow more precise
definitions of the embedded spicule wave modes, which will allow the high-
frequency Alfvén modes to be examined and compared to the models put forward
by Sterling & Hollweg (1984).

In order to establish if the wave energy is dissipated in the form of
localized heating, measurements of thermal processes in the vicinity of these
spicules are necessary. This may be achieved using differential emission
measures of optically thin coronal EUV observations directly above the
spicules (McIntosh, 2012; Vanninathan et al., 2012). An alternative approach
would be to use the Atacama Large Millimeter/Submillimeter Array (ALMA;
Wootten & Thompson, 2009; Wedemeyer et al., 2016) to find the temperature of
the spicules and the surrounding plasma (Chintzoglou et al., 2021; Henriques
et al., Preprint; Jafarzadeh et al., 2021). Importantly, the timing
information related to the decay of the spicule oscillations would need to be
harnessed to provide both spatial and temporal information to examine
localized temperature fluctuations that may be a result of thermalization
mechanisms. While this is beyond the scope of the present work, it will form
the basis of a follow-up study over the coming months.

The downwardly propagating waves maintain an approximately constant energy
flux through a reduction in both velocity amplitude and phase velocity as they
travel down the spicule, visible in Figures 6 and 8, respectively. It is
likely that this is due to the wave interacting with the denser plasma at
lower heights above the solar limb, resulting in a slower Alfvén speed in
these regions (Okamoto & De Pontieu, 2011). This is not unexpected, as the
theoretical modeling of propagating kink waves in longitudinally stratified
waveguides found that phase velocities and velocity amplitudes decrease with
height (Soler et al., 2011).

It has been proposed that in order to supply the quasi-steady effects needed
to heat the solar atmosphere, the dissipation of short period waves is of
paramount importance (Hasan et al., 2005; Hasan & Van Ballegooijen, 2008; Van
Ballegooijen et al., 2011). The energy flux carried by both short period
($<50$ s) and long period ($\geq 50$ s) waves between each set of adjacent
heights is shown in Figure 9 using red and blue data points, respectively. In
order to calculate the associated energy flux for the propagating wave modes,
new filling factors were calculated by combining the previously used spicule
density filling factor ($5\%$; Morton et al., 2012) with the fraction of waves
which were found to fall into each relevant category (i.e., $<50$ s or $\geq
50$ s). The new filling factors were approximately $2.5\%$, which is a result
of the $50$ s boundary being very close to the average period found at each
height (see Table 1).

It can be seen from Figure 9 that the energy flux of the short period waves is
greater than that of the long period waves for the full set of propagating
waves (upper panel), and both the upwardly propagating (middle panel) and
downwardly propagating (lower panel) waves. For the full set of propagating
waves and the upwardly propagating waves, both the short and long period waves
show a similar energy flux decrease with height as that for the total energy
flux values. The energy flux of both the short and long period downwardly
propagating waves show a similar lack of dependence on atmospheric height,
which is consistent with the total energy flux measurements. This suggests
that both short and long period upwardly propagating waves have the potential
to heat the solar atmosphere, although the short period waves have a larger
energy flux across all heights, giving them a greater potential capacity for
heating.

##  4 Conclusions

The results presented here represent a sizable increase in the statistical
population of examined transverse spicule oscillations. Our use of data with a
time cadence of $\sim 1$ s also allowed for the identification of high
frequency waves, similar to those found by Okamoto & De Pontieu (2011), with
periods as short as $10-20$ s, only now with a significant increase in the
examined population size. Observations with even higher spatial and temporal
resolutions may allow for the detection of even shorter period and smaller-
scale oscillations, and further extend the statistical distributions (see,
e.g., Figure 4) down to even smaller values.

We examined the wave properties of spicule oscillations across multiple
atmospheric heights, which facilitated the calculation of associated phase
speeds, hence allowing us to categorize the waves as either being
upwardly/downwardly propagating or standing. Almost an equal balance was found
between upwardly (45%) and downwardly (49%) propagating waves, in contrast to
the earlier study by Okamoto & De Pontieu (2011), who found that upwardly
propagating waves were dominant in their time series. However, the
observations presented here are in close proximity to the solar active region
NOAA AR12391 and may therefore have distinctly different properties to the
coronal hole observations examined by Okamoto & De Pontieu (2011).

Directional information for the spicule waves allowed the calculation of their
associated energy flux as a function of upwardly and downwardly propagating
waves across a number of atmospheric heights. Energy flux estimates are
relatively consistent across all heights for the waves propagating in a
downwards direction. However, for the upwardly propagating waves, a negative
correlation with height is demonstrated, with the overall energy flux
decreasing at a rate of $-13{\,}200\pm 6500$ W m-2/Mm calculated with a
spicule filling factor of $5\%$ (or at a rate of $-1580\pm 780$ W
m${}^{-2}/$Mm using a lower-limit filling factor of $0.6\%$). The mechanism
responsible may either be due to thermalization of the mechanical wave energy
or mode coupling, although investigation of the proportional contributions of
each mechanism are beyond the scope of this study. If even a small fraction of
the wave energy carried in the transverse waves of the spicules examined is
deposited as thermal energy, then it may significantly contribute to the
$10^{3}-10^{4}$ W m-2 requirements needed to balance the radiative losses of
the chromosphere (Withbroe & Noyes, 1977).

WB, DBJ, and FPK acknowledge support from the Leverhulme Trust via the
Research Project Grant RPG-2019-371. DBJ and SDTG wish to thank Invest NI and
Randox Laboratories Ltd. for the award of a Research & Development Grant
(059RDEN-1), in addition to the UK Science and Technology Facilities Council
(STFC) for the consolidated grant ST/T00021X/1. DBJ also acknowledges funding
from the UK Space Agency via the National Space Technology Programme (grant
SSc-009). VMN acknowledges support from STFC grant ST/T000252/1 and Russian
Foundation for Basic Research (RFBR) grant (grant No. 18-29-21016). SJ
acknowledges support from the European Research Council under the European
Union Horizon 2020 research and innovation program (grant agreement No.
682462) and from the Research Council of Norway through its Centres of
Excellence scheme (project No. 262622). The Dunn Solar Telescope at Sacramento
Peak/NM was operated by the National Solar Observatory (NSO). NSO is operated
by the Association of Universities for Research in Astronomy (AURA), Inc.,
under cooperative agreement with the National Science Foundation (NSF). The
authors wish to acknowledge scientific discussions with the Waves in the Lower
Solar Atmosphere (WaLSA; www.WaLSA.team) team, which is supported by the
Research Council of Norway (project number 262622), and The Royal Society
through the award of funding to host the Theo Murphy Discussion Meeting “High-
resolution wave dynamics in the lower solar atmosphere” (grant Hooke18b/SCTM).

## References

  * Alfvén (1942) Alfvén, H. 1942, Nature, 150, doi: 10.1038/150405d0
  * Alfvén & Lindblad (1947) Alfvén, H., & Lindblad, B. 1947, Monthly Notices of the Royal Astronomical Society, 107, 211, doi: 10.1093/mnras/107.2.211
  * Antolin et al. (2015) Antolin, P., Okamoto, T. J., De Pontieu, B., et al. 2015, Astrophysical Journal, 809, doi: 10.1088/0004-637X/809/1/72
  * Antolin et al. (2018) Antolin, P., Schmit, D., Pereira, T. M. D., De Pontieu, B., & De Moortel, I. 2018, The Astrophysical Journal, 856, 44, doi: 10.3847/1538-4357/aab34f
  * Beckers (1972) Beckers, J. M. 1972, ARA&A, 10, 73, doi: 10.1146/annurev.aa.10.090172.000445
  * Bendat & Piersol (2000) Bendat, J. S., & Piersol, A. G. 2000, Measurement Science and Technology, 11, doi: 10.1088/0957-0233/11/12/702
  * Chintzoglou et al. (2021) Chintzoglou, G., De Pontieu, B., Martínez-Sykora, J., et al. 2021, The Astrophysical Journal, 906, 82, doi: 10.3847/1538-4357/abc9b1
  * de la Cruz Rodríguez et al. (2013) de la Cruz Rodríguez, J., Rouppe van der Voort, L., Socas-Navarro, H., & van Noort, M. 2013, A&A, 556, A115, doi: 10.1051/0004-6361/201321629
  * De Pontieu et al. (2012) De Pontieu, B., Carlsson, M., Rouppe Van Der Voort, L. H. M., et al. 2012, Astrophysical Journal Letters, 752, doi: 10.1088/2041-8205/752/1/L12
  * De Pontieu et al. (2007a) De Pontieu, B., Mcintosh, S., Hansteen, V. H., et al. 2007a, Publications of the Astronomical Society of Japan, doi: 10.1093/pasj/59.sp3.S655
  * De Pontieu et al. (2007b) De Pontieu, B., McIntosh, S. W., Carlsson, M., et al. 2007b, Science, 318, doi: 10.1126/science.1151747
  * Desmars et al. (2009) Desmars, J., Arlot, S., Arlot, J. E., Lainey, V., & Vienne, A. 2009, A&A, 499, 321, doi: 10.1051/0004-6361/200811509
  * Dunn (1969) Dunn, R. B. 1969, S&T, 38, 368
  * Ebadi et al. (2012) Ebadi, H., Zaqarashvili, T. V., & Zhelyazkov, I. 2012, Ap&SS, 337, 33, doi: 10.1007/s10509-011-0811-x
  * Edwin & Roberts (1983) Edwin, P. M., & Roberts, B. 1983, Solar Physics, 88, 179, doi: 10.1007/BF00196186
  * Efron et al. (1979) Efron, B., et al. 1979, Annals of Statistics, 7, 1
  * Erdélyi & Fedun (2007) Erdélyi, R., & Fedun, V. 2007, Science, 318, 1572, doi: 10.1126/science.1153006
  * González Manrique et al. (2020) González Manrique, S. J., Noda, C. Q., Kuckein, C., Cobo, B. R., & Carlsson, M. 2020, Astronomy and Astrophysics, 634, doi: 10.1051/0004-6361/201937274
  * Grant et al. (2015) Grant, S. D. T., Jess, D. B., Moreels, M. G., et al. 2015, ApJ, 806, 132, doi: 10.1088/0004-637X/806/1/132
  * Grant et al. (2018) Grant, S. D. T., Jess, D. B., Zaqarashvili, T. V., et al. 2018, Nature Physics, 14, 480, doi: 10.1038/s41567-018-0058-3
  * Hasan & Van Ballegooijen (2008) Hasan, S. S., & Van Ballegooijen, A. A. 2008, ApJ, 680, 1542, doi: 10.1086/587773
  * Hasan et al. (2005) Hasan, S. S., Van Ballegooijen, A. A., Kalkofen, W., & Steiner, O. 2005, ApJ, 631, 1270, doi: 10.1086/432655
  * He et al. (2009) He, J. S., Tu, C. Y., Marsch, E., et al. 2009, A&A, 497, 525, doi: 10.1051/0004-6361/200810777
  * Henriques et al. (Preprint) Henriques, V. M. J., Jafarzadeh, S., Guevara Gómez, J. C., et al. Preprint.  http://arxiv.org/abs/2109.02374
  * Hollweg (1986) Hollweg, J. V. 1986, Journal of Geophysical Research, 91, doi: 10.1029/JA091iA04p04111
  * Hollweg et al. (1982) Hollweg, J. V., Jackson, S., & Galloway, D. 1982, Sol. Phys., 75, 35, doi: 10.1007/BF00153458
  * Houston et al. (2018) Houston, S. J., Jess, D. B., Asensio Ramos, A., et al. 2018, ApJ, 860, 28, doi: 10.3847/1538-4357/aab366
  * Houston et al. (2020) Houston, S. J., Jess, D. B., Keppens, R., et al. 2020, ApJ, 892, 49, doi: 10.3847/1538-4357/ab7a90
  * Jafarzadeh et al. (2017) Jafarzadeh, S., Solanki, S. K., Gafeira, R., et al. 2017, The Astrophysical Journal Supplement Series, 229, 9, doi: 10.3847/1538-4365/229/1/9
  * Jafarzadeh et al. (2021) Jafarzadeh, S., Wedemeyer, S., Fleck, B., et al. 2021, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 379, doi: 10.1098/rsta.2020.0174
  * Jess et al. (2012a) Jess, D. B., De Moortel, I., Mathioudakis, M., et al. 2012a, ApJ, 757, 160, doi: 10.1088/0004-637X/757/2/160
  * Jess et al. (2010a) Jess, D. B., Mathioudakis, M., Christian, D. J., Crockett, P. J., & Keenan, F. P. 2010a, ApJ, 719, L134, doi: 10.1088/2041-8205/719/2/L134
  * Jess et al. (2010b) Jess, D. B., Mathioudakis, M., Christian, D. J., et al. 2010b, Sol. Phys., 261, 363, doi: 10.1007/s11207-009-9500-0
  * Jess et al. (2015) Jess, D. B., Morton, R. J., Verth, G., et al. 2015, Space Science Reviews, 190, 103, doi: 10.1007/s11214-015-0141-3
  * Jess et al. (2012b) Jess, D. B., Shelyag, S., Mathioudakis, M., et al. 2012b, ApJ, 746, 183, doi: 10.1088/0004-637X/746/2/183
  * Jess et al. (2021) Jess, D. B., Snow, B., Fleck, B., Stangalini, M., & Jafarzadeh, S. 2021, Nature Astronomy, 5, 5, doi: 10.1038/s41550-020-1158-4
  * Jess et al. (2020) Jess, D. B., Snow, B., Houston, S. J., et al. 2020, Nature Astronomy, 4, 220, doi: 10.1038/s41550-019-0945-2
  * Khomenko & Collados (2015) Khomenko, E., & Collados, M. 2015, Living Reviews in Solar Physics, 12, 6, doi: 10.1007/lrsp-2015-6
  * Kuridze et al. (2012) Kuridze, D., Morton, R. J., Erdélyi, R., et al. 2012, The Astrophysical Journal, 750, doi: 10.1088/0004-637X/750/1/51
  * Kuridze et al. (2021) Kuridze, D., Socas-Navarro, H., Koza, J., & Oliver, R. 2021, ApJ, 908, 168, doi: 10.3847/1538-4357/abd100
  * Kuridze et al. (2021) Kuridze, D., Socas-Navarro, H., Koza, J., & Oliver, R. 2021, The Astrophysical Journal, 908, 168, doi: 10.3847/1538-4357/abd100
  * Liu et al. (2014) Liu, Z. X., He, J. S., & Yan, L. M. 2014, Research in Astronomy and Astrophysics, 14, 299, doi: 10.1088/1674-4527/14/3/004
  * Makita (2003) Makita, M. 2003, Publications of the National Astronomical Observatory of Japan, 7, 1
  * McIntosh (2012) McIntosh, S. W. 2012, Space Science Reviews, 172, 69, doi: 10.1007/s11214-012-9889-x
  * Morgan & Druckmüller (2014) Morgan, H., & Druckmüller, M. 2014, Solar Physics, 289, 2945, doi: 10.1007/s11207-014-0523-9
  * Morton et al. (2013) Morton, R. J., Verth, G., Fedun, V., Shelyag, S., & Erdélyi, R. 2013, Astrophysical Journal, 768, doi: 10.1088/0004-637X/768/1/17
  * Morton et al. (2014) Morton, R. J., Verth, G., Hillier, A., & Erdélyi, R. 2014, Astrophysical Journal, 784, doi: 10.1088/0004-637X/784/1/29
  * Morton et al. (2012) Morton, R. J., Verth, G., Jess, D. B., et al. 2012, Nature Communications, 3, doi: 10.1038/ncomms2324
  * Nakariakov & Roberts (1995) Nakariakov, V. M., & Roberts, B. 1995, Sol. Phys., 159, 213, doi: 10.1007/BF00686530
  * Nakariakov et al. (2021) Nakariakov, V. M., Anfinogentov, S. A., Antolin, P., et al. 2021, Space Sci. Rev., 217, 73, doi: 10.1007/s11214-021-00847-2
  * Nikol’Skii & Sazanov (1967) Nikol’Skii, G. M., & Sazanov, A. A. 1967, Soviet Ast., 10, 744
  * Nikolsky & Platova (1971) Nikolsky, G. M., & Platova, A. G. 1971, Sol. Phys., 18, 403, doi: 10.1007/BF00149062
  * Okamoto et al. (2015) Okamoto, T. J., Antolin, P., De Pontieu, B., et al. 2015, Astrophysical Journal, 809, doi: 10.1088/0004-637X/809/1/71
  * Okamoto & De Pontieu (2011) Okamoto, T. J., & De Pontieu, B. 2011, Astrophysical Journal Letters, 736, doi: 10.1088/2041-8205/736/2/L24
  * Pasachoff et al. (1968) Pasachoff, J. M., Noyes, R. W., & Beckers, J. M. 1968, Solar Physics, 5, 131
  * Pascoe et al. (2012) Pascoe, D. J., Hood, A. W., De Moortel, I., & Wright, A. N. 2012, Astronomy and Astrophysics, 539, doi: 10.1051/0004-6361/201117979
  * Pascoe et al. (2013) —. 2013, Astronomy and Astrophysics, 551, doi: 10.1051/0004-6361/201220620
  * Pascoe et al. (2010) Pascoe, D. J., Wright, A. N., & De Moortel, I. 2010, Astrophysical Journal, 711, 990, doi: 10.1088/0004-637X/711/2/990
  * Pereira et al. (2012) Pereira, T. M., De Pontieu, B., & Carlsson, M. 2012, Astrophysical Journal, 759, doi: 10.1088/0004-637X/759/1/18
  * Pereira et al. (2013) Pereira, T. M. D., De Pontieu, B., & Carlsson, M. 2013, ApJ, 764, 69, doi: 10.1088/0004-637X/764/1/69
  * Pesnell et al. (2012) Pesnell, W. D., Thompson, B. J., & Chamberlin, P. C. 2012, Sol. Phys., 275, 3, doi: 10.1007/s11207-011-9841-3
  * Poynton (2003) Poynton, C. 2003, Digital Video and HD: Algorithms and Interfaces (Morgan Kaufmann)
  * Riedl et al. (2021) Riedl, J. M., Gilchrist-Millar, C. A., Van Doorsselaere, T., Jess, D. B., & Grant, S. D. T. 2021, A&A, 648, A77, doi: 10.1051/0004-6361/202040163
  * Rimmele (2004) Rimmele, T. R. 2004, in Advancements in Adaptive Optics, doi: 10.1117/12.551764
  * Rimmele & Marino (2011) Rimmele, T. R., & Marino, J. 2011, Living Rev. Solar Phys, 8.  http://www.livingreviews.org/lrsp-2011-2http://www.nso.edu/
  * Schou et al. (2012) Schou, J., Scherrer, P. H., Bush, R. I., et al. 2012, Sol. Phys., 275, 229, doi: 10.1007/s11207-011-9842-2
  * Secchi (1877) Secchi, A. 1877, Le Soleil, Vol. 2 (Gauthier-Villars)
  * Sharma et al. (2018) Sharma, R., Verth, G., & Erdélyi, R. 2018, ApJ, 853, 61, doi: 10.3847/1538-4357/aaa07f
  * Shetye et al. (2021) Shetye, J., Verwichte, E., Stangalini, M., & Doyle, J. G. 2021, ApJ, 921, 30, doi: 10.3847/1538-4357/ac1a12
  * Simpson & Mayer-Hasselwander (1986) Simpson, G., & Mayer-Hasselwander, H. 1986, A&A, 162, 340
  * Socas-Navarro et al. (2015) Socas-Navarro, H., de la Cruz Rodríguez, J., Asensio Ramos, A., Trujillo Bueno, J., & Ruiz Cobo, B. 2015, Astronomy and Astrophysics, 577, doi: 10.1051/0004-6361/201424860
  * Soler et al. (2011) Soler, R., Terradas, J., Verth, G., & Goossens, M. 2011, ApJ, 736, 10, doi: 10.1088/0004-637X/736/1/10
  * Stein & Leibacher (1974) Stein, R. F., & Leibacher, J. 1974, ARA&A, 12, 407, doi: 10.1146/annurev.aa.12.090174.002203
  * Sterling (1998) Sterling, A. 1998, in ESA Special Publication, Vol. 421, Solar Jets and Coronal Plumes, ed. T.-D. Guyenne, 35
  * Sterling (2000) Sterling, A. C. 2000, Solar Physics, 196, 79
  * Sterling & Hollweg (1984) Sterling, A. C., & Hollweg, J. V. 1984, ApJ, 285, 843, doi: 10.1086/162563
  * Suzuki & Inutsuka (2005) Suzuki, T. K., & Inutsuka, S. 2005, The Astrophysical Journal, 632, L49, doi: 10.1086/497536
  * Terradas et al. (2010) Terradas, J., Goossens, M., & Verth, G. 2010, Astronomy and Astrophysics, 524, doi: 10.1051/0004-6361/201014845
  * Uchida (1961) Uchida, Y. 1961, Publications of the Astronomical Society of Japan, 13, 321
  * Van Ballegooijen et al. (2011) Van Ballegooijen, A. A., Asgari-Targhi, M., Cranmer, S. R., & DeLuca, E. E. 2011, ApJ, 736, 3, doi: 10.1088/0004-637X/736/1/3
  * Van Doorsselaere et al. (2014) Van Doorsselaere, T., Gijsen, S. E., Andries, J., & Verth, G. 2014, Astrophysical Journal, 795, doi: 10.1088/0004-637X/795/1/18
  * Van Doorsselaere et al. (2008) Van Doorsselaere, T., Nakariakov, V. M., & Verwichte, E. 2008, ApJ, 676, L73, doi: 10.1086/587029
  * Vanninathan et al. (2012) Vanninathan, K., Madjarska, M. S., Scullion, E., & Doyle, J. G. 2012, Sol. Phys., 280, 425, doi: 10.1007/s11207-012-9986-8
  * Vaughan & Nowak (1997) Vaughan, B. A., & Nowak, M. A. 1997, The Astrophysical Journal, 474, doi: 10.1086/310430
  * Verth et al. (2011) Verth, G., Goossens, M., & He, J. S. 2011, Astrophysical Journal Letters, 733, doi: 10.1088/2041-8205/733/1/L15
  * Verth & Jess (2016) Verth, G., & Jess, D. B. 2016, Washington DC American Geophysical Union Geophysical Monograph Series, 216, 431, doi: 10.1002/9781119055006.ch25
  * Weberg et al. (2018) Weberg, M. J., Morton, R. J., & McLaughlin, J. A. 2018, The Astrophysical Journal, 852, 57, doi: 10.3847/1538-4357/aa9e4a
  * Wedemeyer et al. (2016) Wedemeyer, S., Bastian, T., Brajša, R., et al. 2016, Space Sci. Rev., 200, 1, doi: 10.1007/s11214-015-0229-9
  * Wedemeyer-Böhm et al. (2007) Wedemeyer-Böhm, S., Steiner, O., Bruls, J., & Rammacher, W. 2007, in The Physics of Chromospheric Plasmas
  * Withbroe & Noyes (1977) Withbroe, G. L., & Noyes, R. W. 1977, Annual Review of Astronomy and Astrophysics, 15, doi: 10.1146/annurev.aa.15.090177.002051
  * Wootten & Thompson (2009) Wootten, A., & Thompson, A. R. 2009, Proceedings of the IEEE, 97, 1463, doi: 10.1109/JPROC.2009.2020572
  * Wöger et al. (2008) Wöger, F., Lühe, O. V. D., & Reardon, K. 2008, Astronomy and Astrophysics, 488, 375, doi: 10.1051/0004-6361:200809894
  * Yao et al. (2017) Yao, J. M., Manchester, R. N., & Wang, N. 2017, MNRAS, 468, 3289, doi: 10.1093/mnras/stx729
  * Zaqarashvili & Erdélyi (2009) Zaqarashvili, T. V., & Erdélyi, R. 2009, Space Sci. Rev., 149, 355, doi: 10.1007/s11214-009-9549-y
  * Zaqarashvili & Skhirtladze (2008) Zaqarashvili, T. V., & Skhirtladze, N. 2008, ApJ, 683, L91, doi: 10.1086/591524
  * Zhang et al. (2012) Zhang, Y. Z., Shibata, K., Wang, J. X., et al. 2012, Astrophysical Journal, 750, doi: 10.1088/0004-637X/750/1/16