11institutetext: Institute of Theoretical Astrophysics, University of Oslo,
P.O. Box 1029 Blindern, NO-0315 Oslo, Norway 22institutetext: Rosseland Centre
for Solar Physics, University of Oslo, P.O. Box 1029 Blindern, NO-0315 Oslo,
Norway
22email<EMAIL_ADDRESS>
# Spicules and downflows in the solar chromosphere

Souvik Bose , 1122 Jayant Joshi, 1122 Vasco M.J. Henriques, 1122 Luc Rouppe
van der Voort, 1122

(XXXX; accepted XXXX)

###### Abstract

Context. High speed downflows have been observed in the solar transition
region (TR) and lower corona for many decades. Despite their abundance, it has
been hard to find signatures of such downflows in the solar chromosphere.

Aims. In this work, we target an enhanced network region which shows ample
occurrences of rapid spicular downflows in the H$\mathrm{\alpha}$ spectral
line that could potentially be linked to high-speed TR downflowing
counterparts.

Methods. We used the $k$-means algorithm to classify the spectral profiles of
on-disk spicules in H$\mathrm{\alpha}$ and Ca ii K data observed from the
Swedish 1-m Solar Telescope (SST) and employed an automated detection method
based on advanced morphological image processing operations to detect such
downflowing features, in conjunction with rapid blue shifted and red shifted
excursions (RBEs and RREs).

Results. We report the existence of a new category of RREs (termed as
downflowing RRE) for the first time that, contrary to earlier interpretation,
are associated with chromospheric field aligned downflows moving towards the
strong magnetic field regions. Statistical analysis performed on nearly 20,000
RBEs and 15,000 RREs (including the downflowing counterparts), detected in our
97 min long dataset, shows that the downflowing RREs are very similar to RBEs
and RREs except for their oppositely directed plane-of-sky motion.
Furthermore, we also find that RBEs, RREs and downflowing RREs can be
represented by a wide range of spectral profiles with varying Doppler offsets,
and H$\mathrm{\alpha}$ line core widths, both along and perpendicular to the
spicule axis, that causes them to be associated with multiple substructures
that evolve together.

Conclusions. We speculate that these rapid plasma downflows could well be the
chromospheric counterparts of the commonly observed TR downflows.

###### Key Words.:

Sun: chromosphere Sun: atmosphere line: profiles methods: statistical analysis
techniques: image processing proper motions

##  1 Introduction

Figure 1: Full FOV of the SST dataset observed on 25 May 2017 at 09:16:58 UTC.
Panels (a) and (b) show the H$\mathrm{\alpha}$ red wing and blue wing images
observed with CRISP at a Doppler offset of +40 km s-1and $-$40 km s-1,
respectively. Spicules are seen in the two panels as dark and elongated
thread-like structures. Panel (c) shows a dense chromospheric canopy of
fibrils in the corresponding CHROMIS Ca ii K line core image and (d) shows the
map of the LOS magnetic field (BLOS) saturated between $\pm$100 G derived from
inversions of the Fe I 6301 and 6302 Å spectral profiles. The direction to
solar North is pointing upwards. An animation of this figure is available at
https://www.dropbox.com/s/nxnz5twzczjwaac/Context_movie1.mp4?dl=0

Downflows are known to commonly occur in the solar atmosphere. High speed
downflows have been observed in the transition region (TR) that can sometimes
last from several hours to even several days (e.g. Gebbie et al., 1981; Dere,
1982). Moreover, stronger plasma downflows with speeds ranging from $60$–$200$
km s-1 have been observed over or in the close vicinity of active regions and
quiet Sun alike, with the latter revealing relatively weaker downflows
(Doschek et al., 1976; Kjeldseth-Moe et al., 1988; Brekke et al., 1997; Peter
& Judge, 1999). Decades of observations have revealed the prevalence of
predominant downflows (or red shifts) in the spectral lines of the TR (Doschek
et al., 1976; Scharmer, 1981; Dere, 1982; Peter & Judge, 1999; Dadashi et al.,
2011). These downflows are considered to play an important role in the mass
and energy balance of the solar atmosphere and are crucial to further
understand the physics of the lower solar corona.

These downflows are predominantly seen mainly in the lower coronal and TR
passbands. Brekke et al. (1997); Peter & Judge (1999) reported that most of
the high speed quiet Sun TR downflows usually vanish at chromospheric
temperatures. Chitta et al. (2016) have found supersonic downflows in the TR
of sunspots, which are highly intermittent in time and location (Nelson et
al., 2020). Statistical analysis of such downflows in $48$ sunspots by Samanta
et al. (2018) show that at most half of them show signatures in the upper
chromospheric spectral lines such as Mg ii k, whereas the rest are limited to
the TR. Sunspots (and active regions in general) occupy only a very small
fraction of the solar surface compared to the non active regions. Therefore,
it has remained a mystery as to what happens to the remaining high speed
downflows in the chromosphere that are so predominant in the TR of non active
regions. Why have their signatures been so elusive in the chromosphere? A
possible reason for such lack of evidence could be because the plasma density
in the chromosphere is several orders of magnitude higher compared to the
lower corona or the TR that causes these downflows to simply breakdown as soon
as they reach chromospheric heights.

The chromosphere is a layer of the Sun’s atmosphere that is sandwiched between
the visible photosphere and the million degree corona. It is dominated by a
myriad of different features such as spicules, dynamic fibrils, and filaments
to name a few. All these features contribute towards the vigorous dynamical
changes that the chromosphere experiences on time scales ranging from seconds
to minutes. The chromosphere also plays a crucial role in mass loading and
heating the solar corona since all the non-thermal energy that is responsible
for these mechanisms, propagates through the chromosphere. Only a small
fraction of this energy escapes into the corona or the solar wind while the
majority remains trapped within (Carlsson et al., 2019).

One of the most abundant and ubiquitous features that appear in the
chromosphere are spicules. Spicules are thin, elongated, thread-like, and
highly dynamic structures that are omnipresent in the solar atmosphere both in
active and non active regions. To our knowledge, the first lithographed
drawings resembling spicules were performed by Pietro Tachinni (see, for
historic references and a modern reproduction, Chinnici, 2006) and Angelo
Secchi (Secchi, 1871), both in April 1871, with the drawings in the second
edition of Le Soleil (Secchi, 1877) being better known examples of these early
observations. Spicules are visible all over the solar surface when observed in
chromospheric spectral lines and in some cases are known to reach coronal
temperatures. They have been of great interest to the solar physics community
for a long time. An overview of some of such work can be found at Beckers
(1968); Sterling (2000); Tsiropoula et al. (2012); Hinode Review Team et al.
(2019).

The discovery of a more elusive but energetic category of spicules by De
Pontieu et al. (2007) divided them into two different classes: type-I and
type-II, and stirred a fierce debate in the community regarding their
existence, formation mechanisms and role in coronal heating. The type-Is are
mainly used to refer to the dynamic fibrils found in active regions or mottles
in the quiet Sun, and they usually appear in the close vicinity of strong
magnetic fields (De Pontieu et al., 2004). Further, they display
characteristic parabolic paths in H$\mathrm{\alpha}$ line core space-time
diagrams with typical up and down motions of the order of $10$–$40$ km s-1,
lifetimes between $3$–$5$ min and quasi-periodicities of roughly the same time
periods (Rouppe van der Voort et al., 2007). Advanced numerical simulations by
Hansteen et al. (2006); Heggland et al. (2007) show remarkable similarities
with the observations described and they firmly established that type-I
spicules are formed due to the leakage of photospheric magnetoacoustic
oscillations into the solar chromosphere. The type-II spicules, on the other
hand, are comparatively more dynamic with vigorous sideways motions, high
apparent speeds (80–300 km s-1) and shorter lifetimes (~$1$–$3$ minutes) (De
Pontieu et al., 2007; Pereira et al., 2012, 2016). Moreover, they are known to
be heated as they propagate beyond the chromosphere and become visible in TR
(Pereira et al., 2014; Rouppe van der Voort et al., 2015) and coronal
passbands in both active regions and the quiet Sun (De Pontieu et al., 2011;
Henriques et al., 2016; Kuridze et al., 2016; Samanta et al., 2019). These
qualities render the importance of attributing type-II spicules towards mass
loading and heating the solar corona (see, e.g., Martínez-Sykora et al., 2017;
Kontogiannis et al., 2018; Martínez-Sykora et al., 2018), however, their
detailed physical processes remain far from known.

The spectral signatures of the on-disk counterparts of type-II spicules were
observed for the first time by Langangen et al. (2008). They reported sudden
rapid excursions in the blue wing of Ca II 8542 Å spectral line that led them
to be termed as rapid blue shifted excursions (RBEs). Later, Rouppe van der
Voort et al. (2009); Sekse et al. (2012); Pereira et al. (2016); Bose et al.
(2019a) observed high resolution on-disk images and spectra in the
chromospheric H$\mathrm{\alpha}$, Ca ii 8542 Å, and Ca ii K lines associated
with them and established beyond doubt that RBEs are indeed the on-disk
counterparts of the earlier known type-II spicules.

De Pontieu et al. (2012) reported the presence of torsional motions in type-II
spicules that explained the occurrence of spicular bushes in the red wing
images of chromospheric spectral lines, that were morphologically similar to
their blue counterparts. Later, Sekse et al. (2013b) also described the
existence of red wing counterparts of RBEs in H$\mathrm{\alpha}$ and Ca ii
8542Å with very similar properties as the former and termed them rapid red
shifted excursions (RREs). Next, Kuridze et al. (2015a); Rouppe van der Voort
et al. (2015) and more recently Bose et al. (2019a) also confirmed their
existence, both in the chromosphere and the TR, with high resolution on-disk
observations.

The complex twisting and swaying found in type-II spicules, in addition to the
flows along the chromospheric magnetic field lines, are important
characteristics of spicules that represents outward propagating Alfvénic waves
which can be of the order of several hundred km s-1 and can cause heating in
the hotter TR lines as they propagate (De Pontieu et al., 2014). Like De
Pontieu et al. (2012), Sekse et al. (2013b) and later Kuridze et al. (2015a)
used the transverse motion of spicules to argue for the existence and
appearance of RREs. They interpreted that RREs (though less abundant), like
RBEs, are a manifestation of the same phenomenon with very similar statistical
properties and occur when the latter harbor such complex motions that can
sometimes result in a net red shift when observed on the disk. Therefore, RBEs
can transition to RREs and vice-versa depending on the orientation between the
line-of-sight (LOS) and transverse motion of the structures. Moreover, the
torsional motions can also cause RBEs and RREs to be in close association with
each other (Sekse et al., 2013b; De Pontieu et al., 2014). Both RBEs and RREs
are generally seen to originate in the vicinity of strong network regions with
enhanced magnetic fields, and appear to propagate away as they evolve (Rouppe
van der Voort et al., 2009; McIntosh & De Pontieu, 2009; Sekse et al., 2013b;
De Pontieu et al., 2014; Kuridze et al., 2015a)

Movies of high-resolution H$\mathrm{\alpha}$ blue wing filtergrams at Doppler
offsets $\gtrsim$30 km s-1 are generally dominated by RBEs that appear to move
away from the magnetic network. Corresponding red wing movies at equivalent
Doppler offset, are not nearly as much dominated by outward moving RREs, but
also show elongated absorption features that appear to be moving downward to
the magnetic network. In this paper, with the help of high resolution on-disk
observations from the Swedish 1-m Solar Telescope (SST, Scharmer et al.,
2003), we aim to characterize these returning flows that appear
morphologically similar to RBEs and RREs, but seem to be downflowing in
nature. In the following sections, we describe the methods undertaken to
detect and further investigate their dynamical characteristics and spatio-
temporal evolution. We also discuss their interpretation and the possible
physical processes that could be responsible for their appearance and suggest
that these downflows could be a representative of the chromospheric
counterparts of the lower coronal and TR downflows.

##  2 Observations and data reduction

We observed an enhanced network region close to disk center in a coordinated
SST and IRIS campaign on 25 May 2017 shown in Fig. 1. The heliocentric
coordinates were $(x,y)=(45\arcsec,-93\arcsec)$ with corresponding observing
angle $\mu=\cos\theta=0.99$ ($\theta$ being the heliocentric angle). The
temporal duration was close to 97 min starting from 09:12 UT until 10:49 UT.
We acquired imaging spectroscopic data in H$\mathrm{\alpha}$ and
spectropolarimetric data in Fe I 6302 Å from the CRisp Imaging
SpectroPolarimeter (CRISP, Scharmer et al., 2008). The CHROmospheric Imaging
Spectrometer (CHROMIS) was used to obtain imaging spectroscopic data in Ca ii
K. Both CRISP and CHROMIS are tunable Fabry-Pérot instruments installed at the
SST with the first light of the latter being in 2016. CRISP sampled
H$\mathrm{\alpha}$ at 32 wavelength positions between $\pm$1.85 Å and the Fe I
6301 and 6302 Å line pair at 16 wavelength positions respectively, with a
temporal cadence of 19.6 s and a spatial sampling of
0$\aas@@fstack{\prime\prime}$058\. The Fe I spectral profiles were subjected
to an early version of a robust Milne-Eddington (ME) inversion scheme based on
a parallel `C++`/`Python` implementation111
https://github.com/jaimedelacruz/pyMilne (de la Cruz Rodríguez, 2019). We use
the LOS magnetic field component derived from these inversions in our analysis
as shown in panel (d) of Fig. 1. CHROMIS sampled Ca ii K at 41 wavelength
positions within $\pm$1.28 Å with 63.5 mÅ steps. Furthermore, a continuum
position was sampled at 4000 Å. The CRISP H$\mathrm{\alpha}$ red-wing and
blue-wing images for the full field-of-view (FOV) at +40 km s-1 and $-$40 km
s-1 are respectively shown in panels (a) and (b) of Fig. 1. The corresponding
CHROMIS Ca ii K line core image is shown in panel (c). The temporal cadence
and the spatial scale of this data were 13.6 s and
0$\aas@@fstack{\prime\prime}$038\. High spatial resolution down to the
telescopic diffraction limit (given by $\lambda/\mathrm{D}$ =
0$\aas@@fstack{\prime\prime}$08 at 3934Å) was achieved through excellent
seeing conditions, the SST adaptive optics system (Scharmer et al., 2019), and
the Multi-Object Multi-Frame Blind Deconvolution (MOMFBD, van Noort et al.,
2005) image restoration technique.

We used the CRISPRED data reduction pipeline (de la Cruz Rodríguez et al.,
2015) and an early version of the CHROMIS pipeline (Löfdahl et al., 2018) for
further data reduction while including the spectral consistency method
described in Henriques (2012). The images from both the instruments were de-
rotated to account for diurnal field rotation, aligned and de-stretched to
remove warping due to seeing effects before making it ready for scientific
analyses. Later, both the CRISP and the CHROMIS datasets were co-aligned by
cross-correlating the corresponding photospheric wideband channels and were
rotated to align the direction to solar North along the $y$-axis. The wideband
channel for CRISP has a full-width at half maximum (FWHM)=4.9 Å centered at
H$\mathrm{\alpha}$, whereas for CHROMIS the FWHM equals 13.2 Å centered at
3950 Å between the Ca II H and Ca II K lines.

These observations were first presented by Bose et al. (2019a) where the
cotemporal IRIS observations were included in the analysis. The orientation
and total area covered by the dataset matches the IRIS observations. This data
will be made available in the future as a part of the SST and IRIS database
(Rouppe van der Voort et al., 2020). A brief overview of the dataset and
targetted features of interest are given in the following section.

###  2.1 Overview of the dataset

Figure 1 shows an overview of the FOV and the qualification, enhanced network,
is illustrated by extended patches of strong magnetic field that are
predominately of negative magnetic polarity in the $B_{\mathrm{LOS}}$ map. The
Ca ii K line core image is dominated by a dense canopy of chromospheric
fibrils that covers most of the FOV. The H$\mathrm{\alpha}$ far wing images at
$\pm$40 km s-1 Doppler offset show a large number of elongated absorption
features, and a close look at the animation associated with this figure
clearly reveals that many display the rapid and complex dynamical evolution
that is characteristic for the on-disk counterparts of type-II spicules. More
interestingly, when comparing the evolution in the red and blue wing images,
we see that a majority of the dark threads in the red wing move towards the
strong network field concentrations, opposite to the familiar RBE dynamics
that is predominantly outward in the blue wing images.

These structures, which we term as downflowing RREs, are described here for
the first time and form the major theme of this paper. They are characterized
and compared with the traditional RREs and RBEs which have been widely been
observed in the recent past.

Figure 2: All 50 RPs of H$\mathrm{\alpha}$ spectral lines considered in the
final configuration of the $k$-means algorithm. The RP number is indicated in
the lower right corner. RPs 0–3 are colored in blue signifying blue-ward
excursions whereas 4–7 are colored in red indicating red-ward excursions. The
remaining RPs (8–49) are drawn in grey. The dashed gray spectral line in each
of the panels shows the mean H$\mathrm{\alpha}$ profile averaged over the
entire data set and serves as reference. The magenta profiles show the
positive part of the difference between the mean H$\mathrm{\alpha}$ profile
and the respective RP in each panel. The dashed horizontal magenta line marks
the level of the peak (IP-diff) of the differential profiles, and the vertical
solid magenta line marks the COG Doppler offset of the positive value of the
differential profiles ($\lambda_{\mathrm{COG-diff}}$). The width $W$ of the
H$\mathrm{\alpha}$ line core is given for each RP in Å.

##  3 Method of analysis

###  3.1 Characterizing RBEs and RREs/downflowing RREs from their spectral
profiles

The first and foremost task in our analyses was to identify the rapid blue and
red excursions in our data. We followed the $k$-means clustering method as
described in Bose et al. (2019a) to identify the different H$\mathrm{\alpha}$
and Ca ii K profiles that occur during these rapid excursions.

The method works by partitioning similar observations into $k$ number of
groups or clusters. Each observation is assigned to the cluster with the
nearest mean (also called cluster center). It is an iterative algorithm whose
main objective is to minimize the sum of distances between the observed points
and their respective cluster centers. In our data, each observation (pixel) is
a spectral profile which is compared with the mean spectral profiles of each
of the $k$ clusters before assigning it to one. An advantage of using this
technique is that it relies on the complete spectral signature of the
different features on the FOV rather than their appearance at a particular
wavelength position. This enables efficient detection and characterization of
the features as discussed below.

Thus, we applied this algorithm to cluster intensity profiles of each pixel on
the FOV on the combined H$\mathrm{\alpha}$ and Ca ii K data. The detailed
algorithmic steps, data pre-processing methods, followed by finding the
optimum number of clusters, among which the data can be divided, has been
discussed in detail in Bose et al. (2019a). In this paper, we leverage the
analysis performed in Bose et al. (2019a) and proceed with the 50 clusters
among which the data had been grouped. Each cluster has a representative
profile (RP) that corresponds to the mean over all profiles in that cluster.
Naturally, each and every pixel in the FOV is uniquely assigned to a
particular cluster that can be represented in the form of a 2D map, like the
one shown in panel (c) of Fig. 1 in Bose et al. (2019a) (also see panel (d) of
Fig. 5). We base further analyses, performed in this study, solely on the
H$\mathrm{\alpha}$ RPs and the quantities extracted from them. Figure 2 shows
all 50 H$\mathrm{\alpha}$ RPs which resulted from the $k$-means clustering.
With the exception of RPs-$13$ and $21$, that coincide with the gray borders
outside the common CRISP and CHROMIS FOV in Fig. 1, all the remaining RPs
correspond to various features observed in the FOV. In the following
subsections, we describe the methods undertaken to first identify the RBE and
RRE/downflowing RRE RPs and use them further to detect on-disk spicules in our
dataset.

####  3.1.1 Identifying the RPs

RBEs and RREs have typically been observed in H$\mathrm{\alpha}$ and their
spectral profiles show significant absorption asymmetries in the blue and red
wing positions compared to an average quiet Sun profile (Rouppe van der Voort
et al., 2009). These asymmetries are a sign of velocity gradients that are
commonly observed in spicules. Preliminary analysis with `CRISPEX` (Vissers &
Rouppe van der Voort, 2012), a widget-based analysis tool written in
Interactive Data Language (IDL), suggests that downflowing RREs have spectral
signatures similar to the traditional RREs with an absorption asymmetry in the
red wing of H$\mathrm{\alpha}$. Therefore, the distinction between RREs and
downflowing RREs was not made while identifying their characteristic RPs.

We begin our identification strategy by computing the differential profiles,
i.e. the difference between average H$\mathrm{\alpha}$ profile and the RP,
similar to Rouppe van der Voort et al. (2009) for each cluster. By considering
only the positive values of these differential profiles, we have a measure to
determine the enhanced absorption part. The positive part of the differential
profiles are shown in magenta in Fig. 2. This is followed by determining the
center-of-gravity (COG) and the peak intensity of all these magenta
differential profiles (their values are indicated by the solid vertical and
horizontal dotted magenta lines, respectively). We used the method described
in Uitenbroek (2003) for computing the COG, and its position was used as a
measure of the Doppler offset of an absorption feature in the RPs. The
combination of large values for COG and positive peak of the difference
profile serves as an effective selection criterion to determine enhanced
absorption in the H$\mathrm{\alpha}$ wings.

Pereira et al. (2016) highlighted the importance of H$\mathrm{\alpha}$ line
width in the detection and analysis of spicules. They reported that the
statistical properties measured solely from the H$\mathrm{\alpha}$ line wing
images at selected wavelength positions would result in an underestimate of
the extracted physical quantities because spicules show a large range of
Doppler offsets during their evolution. For such cases line width maps provide
a more robust tool for spicule detection. Typically, RBEs and RREs show
enhanced line widths in addition to large Doppler offsets, with the fastest
ones having even broader line profiles. This indicates that both Doppler
offset and the line width are important properties that should be considered
while detecting spicules. Therefore, we leverage the detection of RBEs, and
RREs/downflowing RREs by combining the two factors. We choose the method
described in Cauzzi et al. (2009) to calculate the widths of the
H$\mathrm{\alpha}$ line core over the full-width at half maximum (FWHM) method
utilized by Pereira et al. (2016), because FWHM tends to mix the photospheric
and the chromospheric signals, which is effectively avoided in the former.
Nevertheless, spectral profiles of spicules show enhancements in both FWHM and
line core width of H$\mathrm{\alpha}$ rendering the validity of both
approaches.

Figure 3: Scatter plot between the COG Doppler offsets and the peak intensity
(IP-diff) of the differential profiles for all RPs. The points are annotated
with the respective RP index and color coded based on the H$\mathrm{\alpha}$
line-core width. The horizontal dashed red and blue lines indicate a Doppler
offset of $\pm$20 km s-1 as a reference.

Figure 3 shows a scatter plot between the COG Doppler offsets and peak
intensities (IP-diff) of the magenta difference profiles, the color of the
data points follow their respective line core widths. Such a representation
takes into account all possible spectral characteristics important for RBE/RRE
qualification, thereby leading to an efficient and a robust characterization.
For visualization purpose, we choose to restrict the lower limit of the
colormap to a line core width of $1.1$ Å. Thereafter, we imposed the following
three criteria to consider RPs belonging to RBEs and RREs/downflowing RREs.
They should: (1) have a minimum COG Doppler offset of 20 km s-1, (2) have IP-
diff $\geq$ 0.2 and (3) have line core widths $\geq$ 1.2 Å. Evidently, RPs
$2$, $3$, $6$, and $7$ stand out quite distinctly in Fig. 3, as they not only
have the highest Doppler offsets, but also high values of peak intensity and
line core widths making them clear RBE and RRE/downflowing RRE RPs. Moreover,
the above criteria are also satisfied by RPs-0, 1, 4 and 5, though not as
distinctly, whereas the rest of the profiles clearly do not make the cut.
Therefore, we chose RPs 0–7 for the analyses presented in the rest of the
paper and assign RPs 0–3 to RBEs and 4–7 to RREs/downflowing RREs in
increasing order of the strength of their spectral properties and they are
accordingly arranged in the first 7 panels of Fig. 2. The remaining RPs,
numbered from 8–49, are indicated in gray.

According to Fig. 3, the computed Doppler offsets of RP-0 and RP-4 are
comparable to RPs-2 and 6 respectively, but, the latter show stronger
absorption in their line wings and a stronger shift in their spectral lines
that in turn also contributes to an increased line core width making them
stronger than the rest. Therefore, RPs-0 and 4 are assigned lower in the order
of labeling the RPs in Fig. 2. Additionally, we also exploited the difference
between the spectral properties of RPs 6 and 7 (2 and 3) and RPs 4 and 5 (0
and 1) described earlier, and segregated the detected RBEs and
RREs/downflowing RREs into two compartments with the purpose of comparing
them. Basically, we grouped RPs 6 and 7 (2 and 3) for the red shifted (blue
shifted) events under stronger red (blue) excursions, whereas RPs 4 and 5 (0
and 1) were grouped together under weaker red (blue) excursions in a way such
that the events in one group were unique with respect to the other. Such a
segregation formed the basis of the investigations carried out in Sects. 4.1
and 4.2, where we discuss them further. In addition to the above properties,
we also note from Fig. 3 that the line-core widths of downflowing RRE/RRE-like
RPs are slightly enhanced in comparison to the RBE-like RPs. The difference is
more pronounced for the RPs belonging to the strong excursions with respect to
the weaker ones. We would like to note that the estimates based on the COG
technique serves well to identify RBE/RRE like profiles, but it fails for some
of the other RPs (such as RP-12 and 31) for which this technique is rather
meaningless. Nevertheless, the method adopted in this paper provides one of
the most comprehensive approach so far in characterizing spicular spectral
profiles.

Now, once the RPs have been identified, a question arises as to how well are
they able to represent individual profiles in a particular cluster? Figure 4
sheds light in this direction and shows the RBE and RRE/downflowing RRE-like
spectral profiles in H$\mathrm{\alpha}$ (top two rows) and Ca ii K (bottom two
rows) for each RP in the form of density distributions, with darker shades
indicating a higher number density of spectral lines. The colored solid line
in the different panels shows the mean over all profiles for a particular
cluster, which, in our case, is equivalent to an RP. We clearly see that the
distribution of the spectral profiles belonging to each cluster is narrow and
is mostly concentrated near their respective mean profiles. This indicates
that the identified RPs are able to efficiently describe the profiles in the
respective clusters for both H$\mathrm{\alpha}$ and Ca ii K.

Figure 4: Density plots of the H$\mathrm{\alpha}$ (top two rows) and Ca ii K
(bottom two rows) spectra for RBE and RRE/downflowing RRE RPs (0–7). The
density (darker meaning higher concentration of spectra) corresponding to each
RPs shows the distribution of profiles over the whole time series. The solid
lines overplotted correspond to the RPs discussed in the text. The blue (red)
color coding indicates the excursions in the blue ward (red ward) side of
H$\mathrm{\alpha}$ respectively.

From the distribution of the spectral profiles in Ca ii K, we see that for a
vast majority of cases, the strongest Doppler shifted K3 (in both the red and
blue excursions) has a significantly stronger opposite K2 intensity
enhancement relative to its line core. In other words, we see that the
absolute difference between the intensities $I$(K2) - $I$(K3) is correlated
with the shift of K3, as was also shown in Bose et al. (2019a). This is the
case exclusively for the RBE and RRE-like RPs, which also formed an additional
basis for the identification of the RPs in Ca ii K. The rest of the profiles
show no such characteristic behavior (refer to Fig. B.2 in the appendix of
Bose et al., 2019a).

The intensity enhancement in the K2 peaks is inherently linked to lower layer
photons, observed due to the presence of strong velocity gradients in spicules
which remove top-layer opacities at those wavelengths. Such lower layers
probably feature enhanced emission as it is due to increased temperature or
other local source function enhancing effects, but the relation between the K2
and the K3 features, set by velocity and reproducible by a simple model,
reveals the role of the velocity gradients. In Bose et al. (2019a) this was
described as an opacity window effect due to the two-dimensional presence of
background features imaged in K2 wavelengths and to distinguish it from other
effects such as the reflector effect (Scharmer, 1981, 1984). Similar
enhancements are generally observed whenever there is a gradient in the LOS
velocity in strong resonant lines such as Ca ii K or Mg ii k (Carlsson &
Stein, 1997; de la Cruz Rodríguez et al., 2015; Kuridze et al., 2015b). A
further example in a very different solar feature is found in umbral flashes
in the Mg ii k spectral line (Bose et al., 2019b). For a recent discussion on
this topic see Sect. 4.1 of Henriques et al. (2020).

###  3.2 Spicule detection and Morphological operations

####  3.2.1 Halos around spicules and their substructures

Figure 5: Overview of the detected RBEs and RREs/downflowing RREs in
H$\mathrm{\alpha}$ wing and Ca ii K. Panels (a), (b) and (c) show RBEs and
their associated substructures (halos) in blue and RREs/downflowing RREs and
their halos in red colored contours against a background of H$\mathrm{\alpha}$
wing images at +40 km s-1, $-$40 km s-1 and Ca ii K line core, respectively at
time t=09:17UT. Panel (d) shows the corresponding RP index map with gradients
in the color indicating the contribution from the different RPs as shown in
the colorbar. An animation of this figure is available at
https://www.dropbox.com/s/wav35xo38vo6s9d/shadow_movie1.mp4?dl=0.  Figure 6:
Overview of the morphological processing techniques described in the text.
Panel (a) shows a portion of the FOV in H$\mathrm{\alpha}$ wing ($-38$ km s-1)
at time t=09:24 UT (scan:50). Blue contours on this image indicate the RBEs
detected with RPs 0–3 as discussed in the text, panel (b) shows the binary
mask of the detected RBEs, and panel (c) indicates the result after performing
connected component labeling of the binary mask in panel (b). The labels are
color coded as seen in the adjoining colorbar. The total number of labels for
RBEs over the full FOV and full duration corresponds to 19,643. Panel (d)
shows the same FOV as in (a) but at time t=10:43 UT (scan:400). Cobalt blue
colored contours in panel (e) show the masks of detected RBEs obtained after
performing morphological labeling, and (f) shows the fitted ellipses with
their centers (green), semi-major axis and semi-minor axis (red) for each
label shown in (e).

Spicules are known to display complex dynamical behavior owing to their: (a)
flows aligned along the chromospheric magnetic field lines, (b) swaying
motions and (c) torsional motions (De Pontieu et al., 2012). Furthermore,
Sekse et al. (2013b) established that many RBEs and RREs exhibit a variation
in Doppler offsets along their lengths and breadths due to which the widths of
H$\mathrm{\alpha}$ spectral lines are enhanced, as was shown by Pereira et al.
(2016). This, in turn is reflected in their spectra which essentially
translates to having multiple spectral profiles (RPs) for the same spicular
body. In this section, we attempt to fully capture this complex behavior by
showing how the different RPs contribute in detecting the complete morphology
of spicules (including the downflowing RREs) and their implications.

Figure 5 and its associated animation shows an overview of the spicules
detected by including multiple RPs (0–7). Panels (a) and (b) respectively show
the contours of the detected RBEs and RREs/downflowing RREs in blue and red
colors against H$\mathrm{\alpha}$ +40 km s-1 and $-$40 km s-1 images. A close
look at panel (b), for example, clearly reveals that in many cases the
overlaid contours enclose a region that appears to be larger than the visible
dark thread like structures, such as the ones around
($\mathrm{X}$,$\mathrm{Y}$) = (27″,5″) or (15″,17″). Similar examples can also
be found in panel (a) for the RREs/downflowing RREs. These lighter shades
around the central dark regions are here termed as spicule halos. They
represent structures with weaker Doppler offsets surrounding the centrally
stronger offset regions. The halos become prominent in the images observed
closer to the H$\mathrm{\alpha}$ core, implying that though they have weaker
Doppler offsets they form a part of the same morphology and evolve together.
Panel (c) shows the same contours as in (a) and (b) but against a background
of a thick chromospheric canopy imaged in Ca ii K line core. As shown in Bose
et al. (2019a), spicules do not have a sharp intensity contrast in Ca ii K
(due to the opacity window effect) unlike H$\mathrm{\alpha}$, which makes it
nearly impossible to observe them as we do in the H$\mathrm{\alpha}$ wing
images. However, looking at the animation of this panel and closely following
the overlaid contours does provide a better impression of the evolution of
spicules in such band-pass.

Further justification of including multiple RPs in the analysis and
visualization become clearer when we look at the RP index map shown in panel
(d) of Fig. 5. It shows the combined RBE and RRE/downflowing RRE RP indices
with varying strengths as set by their their associated Doppler offsets, line
core widths and peak intensities of the respective enhanced differential
profiles. The darker blue (red) shades are an indicator of higher values of
the above three parameters compared to the lighter ones. Clearly, we see that
in several cases the darker shades are accompanied or surrounded by lighter
shades or the previously described halos, in both blue and red shifted
structures. The animation of this figure distinctly indicates that these halos
evolve in conjunction with their corresponding darker cores making it clear
they are a part of the same, but bigger, morphological structure. As seen
earlier for panels (a) and (b), this again strongly indicates widespread group
behavior among spicules. Such behavior was also pointed out by Skogsrud et al.
(2014) among the off-limb spicules. The animation also shows a large number of
downflowing RREs (in red color), co-evolving with their halos, with an
opposite apparent motion compared to the RBEs/RREs (in blue/red color).

The discussion presented in the preceding paragraphs demonstrates that
spicules have multiple substructures that are possible to be detected only
after the inclusion of multiple RPs. It also further evinces that it is
inaccurate to infer about their nature based on detections from a single
wavelength position since their morphology, length and lifetimes could also be
very different in reality. As an example, we show the multi-structural
features among the red and blue excursions belonging to the stronger category
in Appendix A. Such a representation makes the variation in their spectral
properties abundantly clear and prominent.

####  3.2.2 Morphological processing techniques

The discussion presented in the preceding section justifies the importance of
including multiple RPs in spicule detection. This section further advances and
describes the details of the image processing algorithms employed in detecting
the RBEs and RREs/downflowing RREs. The remainder of this section describes
the detailed steps followed in the detection of RBEs in our data. The exact
same procedure was also employed in the detection of RREs/downflowing RREs.

We started off by creating a 3D binary mask (in spatial and temporal domains)
containing all the pixels in the FOV belonging to RPs 0–3 and assigning them a
value of 1 (bright) and the rest 0 (dark). A morphological opening, followed
by a closing operation, with a 3$\times$3 diamond shaped structuring element
was applied to each of the binary masks on a per time step basis. The opening
operation is analogous to an erosion followed by a dilation, and is useful to
remove tiny bright dots in the binary mask. This helps to get rid of small
(1-pixel) connectivity that might be present between different morphological
structures in the 2D space. A morphological opening operation however, can
also create small dark holes in between bright structures which can be
effectively closed by using a closing operation that is the reverse of
opening. It is however important to keep the same 3$\times$3 structuring
element as before. These operations were performed on a per time step basis,
i.e. in 2D, because we intended to preserve the connections in the temporal
domain that would enable us to effectively label them.

The next step was to perform a 3D connected component labeling (Fiorio &
Gustedt, 1996) so that we could identify components uniquely based on a given
heuristic. This technique is widely used in computer vision and image
processing technology. Two pixels are said to be connected when they are
neighbors and have the same numerical value. In this case, we aim to label
pixels in 3D space that are connected and are similar in spectra as set by the
$k$-mean procedure. To not bias for direction we selected an 8-neighbourhood
connectivity in 3D. The top row of Fig. 6 (panels (a)–(c)) provides an
overview of the steps undertaken for the detection of RBEs at the indicated
scan.

When applied to the complete dataset, the method described above led us to
identify 19,643 RBEs and 14,650 RREs (including the downflowing RREs). The two
panels of Fig. 13 in Appendix B, shows the location of all the detected red
and blue excursions respectively, over the co-spatial CRISP and CHROMIS FOV
for the entire time series. On average, most of the spicular activities are
seen to exist in the vicinity of the strong field regions or the network
fields (shown in black colored contours). The difference between the number of
RBE and RRE/downflowing RRE detection is consistent with Sekse et al. (2013b),
with the RRE:RBE detection ratio being ¡ 1. However, it is important to remark
that our red excursions also consists of the downflowing RREs along with the
traditional RREs, and we report far more number of on-disk red and blue
excursions than any of the preceding works cited in this paper, mainly because
of the unique detection method followed by advanced morphological operations
discussed above. Therefore, this allowed us to perform much more exhaustive
data analysis than many such earlier works.

###  3.3 Dimensional analyses and lifetime statistics

One of the major goals of this study is to statistically compare the
properties of the newly reported downflowing RREs to the traditionally known
RBEs and RREs. The detection method undertaken in Sect. 3.2 yielded a large
number of on-disk spicules which provides a perfectly fertile ground to
explore further in this direction. In this section, we focus on the technique
employed to compute their dimensions and lifetimes.

We began by fitting an ellipse to each of the labels obtained after the 3D
connected component labeling, and then computing their lengths of the major
axis, eccentricities ($e$, that is given by $e=\sqrt{1-{b^{2}}/{a^{2}}}$, with
$a$ and $b$ being the semi-major and semi-minor axis of a standard ellipse,
respectively), and the occupied area. Panels (e) and (f) of Fig. 6 provide an
illustrative example where we show the identified spicules in cobalt blue
contours for a given scan (shown in panel (d)), along with their fitted
ellipses–together with their centers and semi-major axis, for a small area on
the FOV. To avoid erroneous detections, a lower limit for the length of the
major axis of detected blue shifted and red shifted features was set at ~100
km or 4 CHROMIS pixels, such that any label with length below the threshold is
not included. Strictly speaking, spicules are not perfect ellipses. However,
they share a common morphology where their lengths (generally) far exceed
their widths (Beckers, 1968; Pereira et al., 2012), making them appear as dark
elongated streaks when observed in H$\mathrm{\alpha}$ on-disk wing images. The
elliptical fitting is performed solely for the purpose of determining the
length occupied by RBEs/RREs by measuring the length of their respective major
axes which ensures that the ends are located at the widest points of the
perimeter of the rapid excursions. Furthermore, elliptical shapes allow a
certain degree of freedom even for those spicules that are not highly
elongated but rather have their lengths only slightly greater than their
widths. These shapes can readily be fitted by ellipses with $0\leq e\leq 1$,
which would mean that the length of the major axis of the ellipses provide a
good approximation of the length of RREs/downflowing RREs and RBEs.

Earlier studies, such as Rouppe van der Voort et al. (2009) and Sekse et al.
(2012, 2013a), relied on obtaining the morphological skeleton of the detected
RBEs and RREs in order to compute their lengths. In such cases, the skeleton
is a thin version of that shape that is equidistant to its boundaries. The
major axis of an ellipse is very similar to the skeleton because it passes
through the COG of the ellipse and therefore their lengths will be comparable.
Both these techniques are basically approximations and skeletonizing certainly
has its own benefits as it preserves the shape of any given structure.
However, since we are interested in determining the maximum extent from both
ends of a feature, the major axis of an ellipse can prove to be advantageous
for complex morphological structures that are often associated with RBEs and
RREs.

It is however important to recall that the 3D connected component labeling
produces a chain of events that are attached both in space and time.
Therefore, if an event (or a label) lasts for multiple time frames, we only
consider the length, area and $e$ when it is at its maximum extent. A similar
approach is followed for spicules with multiple structures or halos around
them. This is justified because inclusion of multiple dimensions for the same
spicule would lead to incorrect statistics. The results obtained after
performing the analysis are shown and discussed in Sect. 4.3.

The lifetimes of spicules were determined on a per event basis, where we
considered the difference between the first and the last occurrence of the
same event in the temporal dimension. The lower limit of the measured
lifetimes is set by the CRISP cadence of ~19.6 s. The results are further
discussed in the latter half of Sect. 4.3.

###  3.4 Apparent motion of spicules

Evidently, the major difference between a downflowing and a traditional RRE is
their plane-of-sky (apparent) motion with respect to the strong magnetic
network areas. In this paper, we investigate the apparent motion of 19,643
RBEs and 14,650 RREs/downflowing RREs on a morphological event-by-event basis,
so as to statistically analyze and describe their trajectories in the plane-
of-sky. Such an analysis would help us to get an idea as to what extent the
red shifted excursions show an opposite trajectory in our datatset, thereby
revealing the abundance of such events. We computed the area weighted COG of
each morphological label (spicule) and followed their evolution in both space
and time. The weighting by area allows the algorithm to follow the trajectory
of the larger substructures in a label. It starts tracking the COG from the
first occurrence of a label and continues until the last. The $X$ and $Y$
coordinates of the COG were stored for each time step and and were then
plotted to display their apparent motion. The results are described further in
Sect. 4.1.

##  4 Results

###  4.1 Downflowing rapid red shifted excursions

Figure 7: Examples of RBEs, RREs and downflowing RREs (labelled as DRREs). The
top panels show their spatio-temporal evolution in H$\mathrm{\alpha}$ blue and
red wing images, whereas the bottom row shows the spectral evolution in
$\lambda t$-diagrams corresponding to the locations marked with white crosses.
Animations are available at
https://www.dropbox.com/sh/3x3kvx158u57xho/AAAA5Xm0sH8RSxPyxe8vqA5Wa?dl=0
Figure 8: Plane-of-sky trajectories of RBEs, RREs and downflowing RREs mapped
according to the strength of their excursions towards the blue (top row) and
red (bottom row) wings. The background shows a BLOS map saturated between
$\pm$ 500 G. The rainbow colored trajectories map the travel direction from
the origin in blue to the termination in red. Here, only a zoom in on a
network patch is shown, Fig. 15 shows the full FOV.

The animations linked to the H$\mathrm{\alpha}$ red wing images in Figs. 1 and
5 seem to convey the impression that in many cases the apparent motion of the
events is rather inward moving instead of the characteristic outward movement
(with respect to the background network areas) associated with traditional
RREs. Except for their apparent motion, we find that these events are very
similar in their appearance and morphology of RBEs and RREs. In this section,
we report clear observational and statistical evidences of these new class of
events termed as downflowing RREs.

Figure 7 shows two examples each of RBEs, RREs and downflowing RREs (labelled
DRREs in the figures). The temporal evolution of these events are shown along
the column for each case and they immediately suggest that, morphologically,
the downflowing RREs are very similar to both RBEs and RREs, and like RREs
they appear in the far red wings of the H$\mathrm{\alpha}$ line core at +40 km
s-1. However unlike the former, the downflowing RREs clearly move towards the
bright network structures as they evolve, whereas the traditional RBEs and
RREs move outwards in the opposite direction. The online animation associated
with each of the examples shown in Fig. 7 establishes the scenario quite
convincingly.

The temporal evolution of their corresponding H$\mathrm{\alpha}$ spectra shown
in $\lambda t$-diagrams also display typical type-II spicule-like behavior for
the downflowing RREs, with a sudden development of a highly asymmetric line
profile towards the red side of the line core. This is very similar to the
characteristic RRE $\lambda t$ evolution first reported by Sekse et al.
(2013b). For the sake of completeness, we also show extended $\lambda
t$-diagrams in Appendix 14 for the two downflowing RREs stretching well before
(~550 s) the occurrence of the red excursions. They clearly show no signs of
preceding blue shift that are typical for type-I spicules. Moreover, we also
note that the lifetime of their excursions in the red wing are similar to RREs
and RBEs. Furthermore, the redshifts associated with the downflowing RREs, in
addition to their inward apparent motion, strongly suggests that these are
real plasma flows (moving away from the observer) and not simply apparent
motions that are often associated with the TR network jets (Tian et al.,
2014). The high apparent speeds in the network jets (of the order of 100–300
km s-1) are most likely not caused by real mass flows. Instead, they are most
likely due to heating fronts propagating at Alfvenic speeds (De Pontieu et
al., 2017). Therefore, appertaining to their spectral and morphological
similarities, we suggest that the downflowing RREs are basically like RREs
simply with an opposite plane of sky motion.

To rule out the possibility that these are singular events, we performed a
statistical study based on the COG tracking method described in Sect. 3.4,
where each and every event belonging to the blue and red side of the
H$\mathrm{\alpha}$ line core was tracked individually. Both the blue and red
excursions were first segregated into two compartments, according to the
strength of their spectral features, by grouping them in the manner described
in Sect. 3.1.1. Consequently, after morphological processing operations,
events that have at least one stronger RP (belonging to either 2 and 3 or 6
and 7) over their whole lifetime are considered under the stronger excursion
category, whereas the rest are grouped into weaker excursions. Panels (a)–(d)
of Figs. 8 and 15 show the apparent trajectories of these excursions as
indicated in their title. Figure 8 shows a zoom in to the central network
patch and Fig. 15 shows the full FOV. The direction of motion is shown by
drawing the trajectories of the COG in a rainbow colormap, where the blue
marks the origin and red the final destination. For the sake of clarity, we
display only those events that have sufficient apparent displacement: $\geq$
0$\aas@@fstack{\prime\prime}$5\. ($\geq$ 1″for the full FOV in Fig. 15). The
$B_{\mathrm{LOS}}$ map is shown as background in order to enhance the
visibility of the trajectories and to facilitate better understanding of the
apparent motion of the events with respect to the strong magnetic field
regions.

We immediately notice that almost all the events in panels (a) and (b)
originate close to the network region and move further outwards as they
evolve–a behavior that is typical to RBEs. The bottom row, on the other hand,
shows that the paths traversed by the red excursions are predominantly
opposite to their blue shifted counterparts. Panel (c), for example, shows an
inward apparent motion for all the red excursions. Moreover, the origin of the
majority of such events can be seen to lie outside the enhanced magnetic field
region and they tend to terminate within or in close proximity to the boundary
of the strong network regions. The weaker red excursions in panel (d) mostly
shows a mixed-bag scenario where a large number of events show the inward
apparent motion, but many among them trace the traditional RRE trajectories,
such as the ones around ($X$,$Y$) = (22″,32″) or (35″,45″). A similar trait is
also observed for a vast majority of events detected in the full FOV, shown in
panels (c) and (d) of Fig. 15, which strongly suggests that the newly reported
downflowing RREs are widely prevalent in our dataset in conjunction with RBEs
and RREs.

###  4.2 Spatial distribution of the stronger and weaker excursions

Figure 9: Spatial occurrence of the stronger and weaker blue and red
excursions with respect to the strong field network regions. Panels (a) and
(b) show the stronger and weaker blue shifted excursions (RBEs) whereas panels
(c) and (d) shows the distribution of the stronger and weaker red shifted
excursions (RREs/downflowing RREs), respectively. The colors represent the
number density of the events shown. The events are mutually exclusive, meaning
that the excursions shown in any one panel are unique and are not related to
the events in the other. The black contour indicates the regions with an
absolute LOS magnetic field $\geq$ 100 G.

The trajectories of the strong and weak excursions presented in the preceding
section sparks interest in investigating their detailed spatial occurrences
with respect to the background network areas over the full FOV. Traditionally,
RBEs and RREs are known to appear in the close vicinity of strong magnetic
field network regions which are also thought to be their foot-points (Rouppe
van der Voort et al., 2009). Therefore, it is worthwhile to explore if such a
behavior is also seen among the stronger and weaker excursions detected in
this study. Figure 9 shows the occurrence of these excursions in the form of
2D density maps. Once again, the events were grouped in exactly the same way
as in Sect. 4.1, with the stronger excursions belonging to the group with RPs
2 and 3 (6 and 7), shown in panels (a) and (c), and the weaker ones belonging
to RPs 0 and 1 (4 and 5) panels (b) and (d), in such a way that events
belonging to these categories are mutually exclusive with respect to one
another.

A closer examination of Fig. 9 reveals that stronger excursions are located
closest to the enhanced network regions (indicated by black contours), whereas
their weaker category counterparts are located further outwards. Moreover,
panels (b) and (d) also suggest that the weaker spicular excursions appear to
exist all over the FOV, but their number density is, on average, roughly
10–15% lesser than their stronger counterparts. Panels (a) and (c) also
highlight an important difference between the stronger blue and red
excursions. In panel (a) it appears that the density of blue excursions are
mostly concentrated outside the network regions and appear to spread outwards
and away from them, whereas the red excursions (including the downflowing
RREs) in panel (c) are mostly located on or within the boundaries of the
strong network regions.

###  4.3 Statistical properties

Figure 10: Dimensional analysis and lifetime statistics of RREs/downflowing
RREs in our dataset. 1D histograms of length, area covered, eccentricity, and
lifetime are shown in (a), (b), (c) and (d) respectively. Moreover, panel (d)
also shows the ECDF of the lifetime distribution in solid green and the dashed
black horizontal line indicates the 98% mark. The vertical dashed lines in
panels (a), (b) and (d) indicates the spatial and temporal resolution limits
of our data. Panels (e)–(h) show the multivariate JPDFs between various
quantities as indicated in a rainbow colormap with red (blue) indicating
highest (lowest) density regions. The magenta contour overlaid on each of the
density distribution indicates the region within which 70% of the events lie.
Figure 11: Same format as Fig. 10 above but for RBEs.

Following the method discussed in Sect. 3.3, we computed the maximum lengths,
maximum areas, the corresponding eccentricities of the fitted ellipses and the
lifetimes of 14,650 downflowing RREs/RREs and 19,643 RBEs. This section is
dedicated to the description and comparison of the statistical properties of
the newly reported downflowing RREs with respect to the traditional RBEs and
RREs.

From the 1D histograms of length and area shown in panels (a) and (b)
respectively, of Figs. 10 and 11, we find that the maximum length of the RBEs
varies between 0.102 Mm to 7.8 Mm whereas for their red shifted counterparts,
an upper limit of 7.75 Mm was found. In both cases a lower threshold length of
4 CHROMIS pixels was imposed which roughly translates to 0.1 Mm. The area
occupied by RBEs range from 0.003 Mm2 to 5.53 Mm2, while the RREs/downflowing
RREs occupied an area that have a maximum value of 7.38 Mm2, with roughly the
same minimum. Both the area and length distributions appear to be skewed with
a large number of data points clustered around 1 Mm.

Panel (c) shows the 1D histogram of eccentricities of the fitted ellipses to
the red and blue shifted excursions. A first glance at them shows that the
events detected in both cases are rather elongated with $e$ $\geq$ 0.5. This
is well aligned to the known morphology of on-disk spicules that has been
described in many studies in the past. Moreover, we also find that highly
eccentric events tend to be more abundant. A close look at the eccentricity
histograms, however, reveal that the distribution appears to be flatter for
the RREs/downflowing RREs in Fig. 10, compared to the RBEs in Fig. 11. This is
further supported by the fact that the median value for the eccentricities
corresponds to 0.85 for the former whereas for the latter it equates to 0.92.
This indicates that, on average the red shifted excursions are slightly less
elongated than their blue shifted counterparts.

Panels (e)–(g) shows the combined joint probability distribution functions
(JPDFs) between various quantities as indicated. The magenta contours specify
the regions within which 70 % of the events lie. The motivation behind these
JPDFs was to highlight the bi-variate relationships between the extracted
quantities in the form of a 2D probability density function. The JPDF between
maximum length vs. maximum area of the blue shifted and red shifted events in
panel (e) show a strong correlation indicating that spicules that occupy
larger area are also longer. Further investigation from panels (f) and (g)
reveal that the bulk (in this case 70 % or more) of the events are highly
eccentric and are less than 1 Mm in length and 0.5 Mm2 in area. Also in
general, events that are more eccentric also tend to be larger in size. This
is true for both the RBEs and RREs/downflowing RREs.

Panels (d) and (h) in Figs. 10 and 11 show the 1D distributions and the JPDFs
between lifetime and length of downflowing RREs/RREs and RBEs, respectively.
We find RBEs lasting from 13.6 s to 2140 s, and RREs/downflowing RREs lasting
from 13.6 s to 1618 s with a median lifetime of about 27.2 s in both cases.
Moreover, we also show an empirical cumulative distribution function (ECDF) of
the lifetime that clearly indicates that 98% of the events are shorter than
200 s. The events lasting longer than $\geq$ 300 s in the above distributions
make up less than 2% in our data and can be considered as outliers. Part of
these are related to a few surge-like events that are morphologically quite
different from type-II spicules but can have similar spectral properties. It
can further not be excluded that some of these are actually multiple recurring
spicules that the method cannot separate and classifies as single long
duration events. The JPDFs in panel (h) also show that more than 70% of the
events that last shorter than 100 s, are also shorter than 0.7 Mm. The results
described above complies well with the characteristic properties of type-II
spicules and are well in agreement with the values reported in earlier works
such as Rouppe van der Voort et al. (2009); Sekse et al. (2012) and Pereira et
al. (2012). Further analyses showing the relationship between the eccentricity
of spicules vs. their lifetimes are shown in Fig. 16 in appendix B, which
reinforces that the majority of the detected events are both highly eccentric
and rapid.

The results presented above do not explicitly distinguish between the
downflowing and traditional RREs, but they strongly suggest that the newly
reported downflowing RREs have properties similar to their upflowing
counterparts observed in both the blue ward and red ward side of the
H$\mathrm{\alpha}$ line core. This further reinforces our proposition that the
downflowing RREs belong to the family of type-II spicules with an opposite
apparent motion.

##  5 Discussion

###  5.1 Interpretation of the downflowing RREs

The examples shown in Fig. 7 and the statistical analysis of the proper motion
of stronger and weaker excursions described in Fig. 8, indicate that there
exists a new category of RREs that behave contrary to the traditional
interpretation of RREs. De Pontieu et al. (2012); Sekse et al. (2013b); De
Pontieu et al. (2014) and Kuridze et al. (2015a) explained that RREs, like
RBEs, are a manifestation of the same physical phenomenon and appear either in
the blue or red wing of H$\mathrm{\alpha}$ depending on their transverse
motion along the LOS of chromospheric magnetic field lines. In other words,
RREs are generally found in close association with RBEs. The downflowing RREs,
discussed in this paper, do not satisfy this explanation since they are not
often found associated with their blue wing counterparts. Moreover, their
opposite apparent direction of motion indicates that they do not occur due to
the swaying and torsional motions often associated with type-II spicules (De
Pontieu et al., 2012).

The motion associated with chromospheric spicules is quite complex and often
all the three, i.e., the field aligned flows, the torsional motions and the
transverse swaying motions, are at play. So far, the appearance of isolated
RBEs are interpreted to be the result of upflows along the field lines
parallel to the spicule axis and the RREs the result of the combination of the
three. The appearance, morphology and their statistical properties discussed
in this paper suggests that downflowing RREs are similar to RBEs and they are
the result of the field aligned downflows in the solar chromosphere.

These downflows, however, are quite different from the type-I spicules such as
active region dynamic fibrils or quiet sun mottles that are commonly observed
near the H$\mathrm{\alpha}$ line core images (Rouppe van der Voort et al.,
2007). The $\lambda$t slices shown at the bottom of Fig. 7 show that the
downflowing RREs have no a priori blue shift associated with them, a
characteristic that is often linked with type-I spicules. Moreover, various
studies such as (De Pontieu et al., 2007; Pereira et al., 2012) show that
type-Is have much lesser LOS velocities (Doppler offsets) and seldom appear at
wing positions so far from the line core. Furthermore, the morphological
appearance of the downflowing RREs are very similar to RBEs and RREs, which
also provides a strong evidence that the former category is not associated
with the red shifts of the fibrils but are more like the downflowing
counterparts of RREs.

###  5.2 What might cause downflowing RREs?

One of the most important questions that remains to be addressed is what might
be the origin of these downflowing RREs? In this section we discuss some of
the possible mechanisms that could be responsible for such ubiquitous
downflows.

Type-II spicules sometimes exhibit parabolic up-down motion much like the
dynamic fibrils (Pereira et al., 2014). However, unlike the latter, the type-
II spicules are much faster, shorter lived (in the chromosphere) and show
signatures in the transition region (TR) and even in the solar corona (De
Pontieu et al., 2011; Henriques et al., 2016; Samanta et al., 2019). During
their ascending phase they get rapidly heated to TR with signatures in
passbands sampling coronal (1 MK) temperatures, in both active regions and the
quiet Sun, and eventually fall back in the chromosphere. These downflowing
plasma could, in principle, be responsible for the observed downflowing RREs
in the H$\mathrm{\alpha}$ red wing. The combined upflowing and downflowing
lifetimes in such cases can well be above 600–700 s (Pereira et al., 2014;
Samanta et al., 2019) which makes them typically last longer than their purely
chromospheric type-I counterparts (lifetime ~3–5 min).

A supporting proposition in favor of this returning plasma in the spicular
form can stem from the fact that observations of the spectral lines formed in
the temperature regimes between ~15,000 K and 2.5 $\times$ 105 K, reveal the
prevalence of an average red shift or downflowing motion of the order of
$10$–$15$ km s-1 in the TR (Doschek et al., 1976; Gebbie et al., 1981; Peter &
Judge, 1999; Dadashi et al., 2011) as discussed in Sect. 1. These studies
imply the existence of plasma flows or wave motions in the quiet Sun with
amplitudes that are significant fractions of the sound speed in the TR.

Different mechanisms have been proposed in the past to explain the net
downflowing structure of the TR. Peter et al. (2006) synthesized spectra of
the coronal and TR lines with the help of a 3D numerical simulation spanning
the photosphere to the corona and showed that the persistent red shifts in the
TR could possibly be explained by the effect of the flows caused due to
heating by magnetic braiding. Hansteen et al. (2010) expanded on this model
and injected an emerging flux to the 3D model of Peter et al. (2006) and
concluded that these downflows are a result of the rapid episodical heating
between the upper chromosphere and lower corona. Zacharias et al. (2018)
further complemented these earlier studies by establishing that pressure
driven downflows along the magnetic field lines could be identified as one of
the key mechanisms responsible for these observed red shifts in the TR.

Perhaps the most likely mechanism applicable to the observations of
downflowing RREs is the proposition that TR red shifts are caused due to the
emission from the return flows, that had been formerly heated and injected
into the solar corona by spicules (Pneuman & Kopp, 1977; Athay & Holzer, 1982;
Athay, 1984). Based on their modeling efforts, the studies by the above sets
of authors estimated that spicular material heated to coronal temperatures,
carry a large upward flux (upflows) that is almost 100 times the flux measured
due to solar wind at 1 AU. Consequently, they reasoned that the rest of the
mass must return to solar atmosphere and it should happen in the form of
spicular downflows which are then in turn responsible for the observed net red
shifts in the TR. Older numerical calculations such as the one by Mariska
(1987), however, could not obtain these observed red shifts. Moreover, many
earlier numerical models could not even account for any blue shift (upflows)
in the coronal or TR passbands. Therefore, the community at large remained
skeptical about the contribution of spicules in these downflows. Recent
observations by McIntosh & De Pontieu (2009); Rouppe van der Voort et al.
(2015) and state-of-the-art numerical modeling efforts by Martínez-Sykora et
al. (2017, 2018) addressed these concerns and proved beyond doubt that
spicules harbor a clear blue shift when observed in hotter TR and coronal
spectral channels.

The observations of these downflowing features have mainly been limited to the
TR till date. None of the studies in the past found evidences of the
chromospheric counterparts of these TR downflows. In some studies (such as
Mariska, 1987) the non existence of the chromospheric signatures of the TR
downflows raised serious doubts on their impact in spicular plasma in the
past. However, the examples of the ubiquitous downflows and the detailed
analysis presented in this paper seem to suggest that the downflowing RREs
could well be the signatures of the TR downflows in the chromosphere. Whether
they follow a typical parabolic trajectory or not warrants further
investigation but our results strongly revives the possibility that returning
spicular materials can be one of major driving mechanisms for the observed TR
red shifts and downflowing RREs.

Alternate possible explanation for the origin of the downflowing RREs could
lie in the investigation carried out by Rutten et al. (2019) where they
conjectured that RBEs could display return flows in the form of cool plasma
that traces the trails of the preceding type-II spicules. Their statistical
analysis strongly suggests that there is a tight correlation between the
occurrences of RBEs and subsequent H$\mathrm{\alpha}$ fibrils within a certain
time delay. Consequently, they attribute the dark fibrilar appearance around
the chromospheric network regions (seen predominantly in H$\mathrm{\alpha}$)
mainly to preceding type-II spicules. The downflowing RREs, like the ones
presented in this paper, could well be the immediate following state of these
firbils as they continue to evolve. However, so far there are no studies that
could firmly establish this relation.

Finally, the chromospheric counterparts of the low lying loops found in the TR
(Hansteen et al., 2014; Pereira et al., 2018) could also be attributed to the
appearance of these downflowing RREs in H$\mathrm{\alpha}$. The low lying
loops are not found to show prominent signals in the H$\mathrm{\alpha}$ line
wing images, except in the far wings, including the red wing. If such red-wing
signatures are true flows then a footpoint could appear as a downflowing RRE
whereas most of the loop would be ”hiding” in TR passbands. The footpoints of
nearly all the loops observed by Pereira et al. (2018) lie in close proximity
or share the footpoints of chromospheric spicules predominantly rooted in the
network regions which makes them possible to have a relation with the
downflowing RREs as described above. Furthermore, Martínez-Sykora et al.
(2020) also discussed the possibility of spicules forming along the loops in
numerical simulations which could be further compared with the observations.
However, given the fact that these loops are not as ubiquitous as spicules on
the solar surface it is very unlikely that downflowing RREs are always
associated with the chromospheric signatures of these low lying TR loops.

###  5.3 Comparison with flocculent flows

Flocculent flows in the solar chromosphere were first described by Vissers &
Rouppe van der Voort (2012) as distinct small-scale features that move
intermittently towards and away from a sunspot. They bear morphological
resemblance to coronal rain, both qualitatively and quantitatively (Antolin &
Rouppe van der Voort, 2012), but their sizes are somewhat smaller and they
move at much lower average velocities and over shorter distances. Flocculent
flows have been suggested to be a result of a siphon flow driven by pressure
difference between the footpoints in a loop. Compared to downflowing RREs,
flocculent flows typically travel over larger distances and consist more of
distinct and isolated blob-like features. We observe downflowing RREs
exclusively in close vicinity to magnetic network areas while flocculent flows
were observed to travel along extended parts of the sunspot superpenumbra and
long active region fibrils.

###  5.4 Spatial distribution of spicules and their substructures

The results presented in Fig. 9 suggests that, on average, the stronger
excursions in both the redward and blueward side of H$\mathrm{\alpha}$ line
core are located closer to the strong magnetic field regions, whereas the
weaker counterparts appear to be scattered across the FOV. This seems to be
the case for both the red and blue shifted excursions. Moreover, the weaker
excursions occupy lesser area on the FOV on an individual event-by-event
basis. Since most of the events in the stronger red excursion category have
apparent inward motion (refer to Figs. 8 and 15), we conclude that on average
the downflowing RREs are predominantly located in the regions close to the
network areas, whereas their traditional counterparts appear to exist slightly
farther away.

The spicule halos, described in Sect. 3.2.1, are a consequence of the
substructures seen in RBEs and RREs/downflowing RREs. These substructures are
mainly due to the fact that spicules often have a large distribution of
spectral properties such as line core width, Doppler shift (Pereira et al.,
2016) or absorption in the wings of the H$\mathrm{\alpha}$ spectral line, not
only along the spicule axis but also in the transverse direction, and are a
part of the same morphological structure that evolves in a collective fashion.
Figure 5 and its associated animation strongly suggests that many times
spicules do not evolve independently but in groups that maybe difficult to
discern when observed at one wavelength position. This reinforces the fact
that group behavior is common among type-II spicules (Skogsrud et al., 2014).

Another important feature commonly observed in spicules is the rapid
morphological transformations that they undergo both in space and time,
including multiple splitting and branching (Yurchyshyn et al., 2020). We refer
to the examples shown in Fig. 7 where we clearly see that the DRRE-1 initially
starts as one structure but towards the end of its evolution we see it clearly
split into two as they terminate at the network bright points. RBE-1 also
shows clear signs of morphological branching with a ”Y” pattern just before it
disappears. Often these branches have different LOS velocities but as we see
in the examples they are clearly associated with the original parent structure
evolving together. Therefore it implies that, downflowing RREs, like RBEs and
RREs, can have multiple structures either in situ, due to a range of spectral
properties, or due to the morphological transformations they undergo during
their evolution.

###  5.5 Significance of the detection technique and their limitations

The detection method based on $k$-means clustering and morphological image
processing technique presented in this paper enabled us to study and analyze
spicules in an unprecedented detail. Identifying the RBE and RRE/downflowing
RRE RPs on the basis of the strength of their Doppler offset, line core width
and enhanced absorption measure in the H$\mathrm{\alpha}$ spectral line
provides one of the most compendious approach in their characterization.

Nearly 20,000 RBEs and 15,000 RREs (including downflowing RREs) have been
identified in our 97 minute long dataset that allowed us to perform varied
statistical analysis and compare the properties of the downflowing RREs in the
context of traditional RBEs and RREs. Though no distinction was possible
between the downflowing RREs and RREs in terms of their spectral signatures,
we see that out of the 15,000 red shifted events a substantial number of
strongest excursions are downflowing in nature, with their apparent motion
directed towards the strong network areas. Furthermore, the statistical
analysis presented in this paper reveals that downflowing RREs have similar
dimensions and lifetimes when compared against RBEs and traditional RREs. This
makes them more likely to be a part of the family of type-II spicules.

Previous reports (Rouppe van der Voort et al., 2009; Sekse et al., 2012;
Pereira et al., 2012; Kuridze et al., 2015a) found comparatively longer
lengths for RBEs and RREs. This could possibly be due to the fact that, except
for Pereira et al. (2012), most of these earlier works have focused on
performing the statistical analyses based on images at single wavelength
positions far in the blue or the red wing of the H$\mathrm{\alpha}$ line
profile. Moreover, on most occasions they employed a lower length detection
threshold of ~$0.730$ Mm that strictly limited the lower estimate of their
analysis, thereby facilitating the detections of relatively longer events. On
the other hand, the method employed in this paper is compendious both in terms
of spectra as well as in size thereby enabling smaller detections.

The technique presented in this paper exploits the complete spectral profile
(RPs) instead of relying on the morphology of spicular features at single
wavelength positions. To get rid of erroneous detections we also impose a
lower limit on the length of the detected events but it is far lower (~0.1 Mm)
than the ones chosen in the earlier works. Rouppe van der Voort et al. (2009)
reported the presence of special RBEs in the form of ”black beads” in
H$\mathrm{\alpha}$ that appeared as tiny roundish darkenings and were
interpreted as spicules that were aligned closely along the LOS. The dimension
of such features were reportedly between 0.15 and 0.3 Mm. Interestingly, Anan
et al. (2010) reported the shortest mean lengths of spicules compared to all
other works in the recent past. Pereira et al. (2012) interpreted this
discrepancy owing it to the exponential drop in intensity along the body of
the spicule which can possibly render the tops too faint to be picked up in
their detection when viewed against the disk. This explanation is also valid
in our case because, despite the spatial resolution, the top portion of the
spicules most likely has too little opacity to be considered as a part of
either an RBE or an RRE/downflowing RRE. Therefore, there could be an
intrinsic bias in determination of the lengths of spicules as was also
reported by Sekse et al. (2012). The lifetimes, on the other hand, are well in
agreement with most of the former studies mentioned above, despite the
increased capture of events over the FOV and increased sensitivity to smaller
scales.

Despite undertaking advanced methods to characterize and detect RBEs and RREs,
we do encounter a few limitations. In this paper, we have referred to the
detected events as RBEs or RREs/downflowing RREs which are on-disk
counterparts of type-II spicules; but are we clearly detecting the type-IIs?
One of the limitations of $k$-means clustering is that it only accounts for
the spectral lines in our case. Therefore, both type-I and type-II spicules
could very well be included in our detections as they have similar spectral
profiles in H$\mathrm{\alpha}$. However, the results from the detailed
statistical analysis presented in this paper strongly suggest that we are in
the domain of type-II spicules since, most importantly, 98% of the events last
under ~$3$ minutes and the median lifetime is about $27$ s. Furthermore the
morphology and the evolution of the detected events, like the ones indicated
in Figs. 6 and 5, clearly suggest that we are detecting the rapid type-II
spicules.

##  6 Concluding remarks

We report the first observation and characterization of rapid downflows in the
solar chromosphere in the form of spicules. Their rapidity and apparent motion
in the far red wing images of H$\mathrm{\alpha}$ strongly suggests that they
are the downflowing counterparts of the traditional RREs first reported by
Sekse et al. (2013b), unlike coronal rain (Antolin & Rouppe van der Voort,
2012) and flocculent flows (Vissers & Rouppe van der Voort, 2012). We
therefore term them downflowing RREs. In depth statistical analyses performed
on 14,650 RREs and 19,643 RBEs imply that downflowing RREs are similar to the
already known RREs and RBEs. Moreover, they also undergo rapid morphological
transformations during their evolution in the same way as RBEs and RREs. The
only evident difference of this new class of RREs lies in their apparent
motion where they seem to originate away and terminate close to the strong
field network regions. This suggests that they are a result of the field
aligned downflows in the solar chromosphere. Furthermore, the downflowing RREs
could also undergo the transverse and torsional motions that are often
associated with type-II spicules which hints at the possibility of finding
downflowing RBEs in the solar chromosphere. Future work in this direction
could shed more light in this context.

The downflowing RREs could possibly be linked to the return flows of type-II
spicules which ascend rapidly from the chromosphere to the TR and even coronal
heights, and eventually fall back. Moreover, we present arguments suggesting
that these downflowing RREs could well be responsible for the observed TR
downflows as was first hypothesized by Pneuman & Kopp (1977). It would be
interesting to find direct signatures of these downflowing RREs in the corona
and the TR with coordinated observations from the Interface Region Imaging
Spectrograph and Solar Dynamics Observatory because it will enable better
understanding of the mass and energy cycle in the solar atmosphere. Such a
study is currently underway and will be the subject of a forthcoming paper.

###### Acknowledgements.

We thank Ainar Drews for his help with the observations. The Swedish 1-m Solar
Telescope is operated on the island of La Palma by the Institute for Solar
Physics of Stockholm University in the Spanish Observatorio del Roque de los
Muchachos of the Instituto de Astrofísica de Canarias. The Institute for Solar
Physics is supported by a grant for research infrastructures of national
importance from the Swedish Research Council (registration number 2017-00625).
This research is supported by the Research Council of Norway, project number
250810, and through its Centers of Excellence scheme, project number 262622.
VMJH is also funded by the European Research Council (ERC) under the European
Union’s Horizon 2020 research and innovation programme (SolarALMA, grant
agreement No. 682462).

## References

  * Anan et al. (2010) Anan, T., Kitai, R., Kawate, T., et al. 2010, PASJ, 62, 871
  * Antolin & Rouppe van der Voort (2012) Antolin, P. & Rouppe van der Voort, L. 2012, ApJ, 745, 152
  * Athay (1984) Athay, R. G. 1984, ApJ, 287, 412
  * Athay & Holzer (1982) Athay, R. G. & Holzer, T. E. 1982, ApJ, 255, 743
  * Beckers (1968) Beckers, J. M. 1968, Sol. Phys., 3, 367
  * Bose et al. (2019a) Bose, S., Henriques, V. M. J., Joshi, J., & Rouppe van der Voort, L. 2019a, A&A, 631, L5
  * Bose et al. (2019b) Bose, S., Henriques, V. M. J., Rouppe van der Voort, L., & Pereira, T. M. D. 2019b, A&A, 627, A46
  * Brekke et al. (1997) Brekke, P., Hassler, D. M., & Wilhelm, K. 1997, Sol. Phys., 175, 349
  * Carlsson et al. (2019) Carlsson, M., De Pontieu, B., & Hansteen, V. H. 2019, ARA&A, 57, 189
  * Carlsson & Stein (1997) Carlsson, M. & Stein, R. F. 1997, ApJ, 481, 500
  * Cauzzi et al. (2009) Cauzzi, G., Reardon, K., Rutten, R. J., Tritschler, A., & Uitenbroek, H. 2009, A&A, 503, 577
  * Chinnici (2006) Chinnici, I. 2006, Memorie della Societa Astronomica Italiana Supplementi, 9, 28
  * Chitta et al. (2016) Chitta, L. P., Peter, H., & Young, P. R. 2016, A&A, 587, A20
  * Dadashi et al. (2011) Dadashi, N., Teriaca, L., & Solanki, S. K. 2011, A&A, 534, A90
  * de la Cruz Rodríguez (2019) de la Cruz Rodríguez, J. 2019, A&A, 631, A153
  * de la Cruz Rodríguez et al. (2015) de la Cruz Rodríguez, J., Löfdahl, M. G., Sütterlin, P., Hillberg, T., & Rouppe van der Voort, L. 2015, A&A, 573, A40
  * De Pontieu et al. (2012) De Pontieu, B., Carlsson, M., Rouppe van der Voort, L. H. M., et al. 2012, ApJ, 752, L12
  * De Pontieu et al. (2004) De Pontieu, B., Erdélyi, R., & James, S. P. 2004, Nature, 430, 536
  * De Pontieu et al. (2017) De Pontieu, B., Martínez-Sykora, J., & Chintzoglou, G. 2017, ApJ, 849, L7
  * De Pontieu et al. (2007) De Pontieu, B., McIntosh, S., Hansteen, V. H., et al. 2007, PASJ, 59, S655
  * De Pontieu et al. (2011) De Pontieu, B., McIntosh, S. W., Carlsson, M., et al. 2011, Science, 331, 55
  * De Pontieu et al. (2014) De Pontieu, B., Rouppe van der Voort, L., McIntosh, S. W., et al. 2014, Science, 346, 1255732
  * Dere (1982) Dere, K. P. 1982, Sol. Phys., 77, 77
  * Doschek et al. (1976) Doschek, G. A., Feldman, U., & Bohlin, J. D. 1976, ApJ, 205, L177
  * Fiorio & Gustedt (1996) Fiorio, C. & Gustedt, J. 1996, Theoretical Computer Science, 154, 165
  * Gebbie et al. (1981) Gebbie, K. B., Hill, F., November, L. J., et al. 1981, ApJ, 251, L115
  * Hansteen et al. (2014) Hansteen, V., De Pontieu, B., Carlsson, M., et al. 2014, Science, 346, 1255757
  * Hansteen et al. (2006) Hansteen, V. H., De Pontieu, B., Rouppe van der Voort, L., van Noort, M., & Carlsson, M. 2006, ApJ, 647, L73
  * Hansteen et al. (2010) Hansteen, V. H., Hara, H., De Pontieu, B., & Carlsson, M. 2010, ApJ, 718, 1070
  * Heggland et al. (2007) Heggland, L., De Pontieu, B., & Hansteen, V. H. 2007, ApJ, 666, 1277
  * Henriques (2012) Henriques, V. M. J. 2012, A&A, 548, A114
  * Henriques et al. (2016) Henriques, V. M. J., Kuridze, D., Mathioudakis, M., & Keenan, F. P. 2016, ApJ, 820, 124
  * Henriques et al. (2020) Henriques, V. M. J., Nelson, C. J., Rouppe van der Voort, L. H. M., & Mathioudakis, M. 2020, A&A, 642, A215
  * Hinode Review Team et al. (2019) Hinode Review Team, Al-Janabi, K., Antolin, P., et al. 2019, PASJ, 71, R1
  * Kjeldseth-Moe et al. (1988) Kjeldseth-Moe, O., Brynildsen, N., Brekke, P., et al. 1988, ApJ, 334, 1066
  * Kontogiannis et al. (2018) Kontogiannis, I., Gontikakis, C., Tsiropoula, G., & Tziotziou, K. 2018, Sol. Phys., 293, 56
  * Kuridze et al. (2015a) Kuridze, D., Henriques, V., Mathioudakis, M., et al. 2015a, ApJ, 802, 26
  * Kuridze et al. (2015b) Kuridze, D., Mathioudakis, M., Simões, P. J. A., et al. 2015b, ApJ, 813, 125
  * Kuridze et al. (2016) Kuridze, D., Zaqarashvili, T. V., Henriques, V., et al. 2016, ApJ, 830, 133
  * Langangen et al. (2008) Langangen, Ø., De Pontieu, B., Carlsson, M., et al. 2008, ApJ, 679, L167
  * Löfdahl et al. (2018) Löfdahl, M. G., Hillberg, T., de la Cruz Rodriguez, J., et al. 2018, arXiv e-prints, arXiv:1804.03030
  * Mariska (1987) Mariska, J. T. 1987, ApJ, 319, 465
  * Martínez-Sykora et al. (2018) Martínez-Sykora, J., De Pontieu, B., De Moortel, I., Hansteen, V. H., & Carlsson, M. 2018, ApJ, 860, 116
  * Martínez-Sykora et al. (2017) Martínez-Sykora, J., De Pontieu, B., Hansteen, V. H., et al. 2017, Science, 356, 1269
  * Martínez-Sykora et al. (2020) Martínez-Sykora, J., Leenaarts, J., De Pontieu, B., et al. 2020, ApJ, 889, 95
  * McIntosh & De Pontieu (2009) McIntosh, S. W. & De Pontieu, B. 2009, ApJ, 707, 524
  * Nelson et al. (2020) Nelson, C. J., Krishna Prasad, S., & Mathioudakis, M. 2020, A&A, 636, A35
  * Pereira et al. (2012) Pereira, T. M. D., De Pontieu, B., & Carlsson, M. 2012, ApJ, 759, 18
  * Pereira et al. (2014) Pereira, T. M. D., De Pontieu, B., Carlsson, M., et al. 2014, ApJ, 792, L15
  * Pereira et al. (2016) Pereira, T. M. D., Rouppe van der Voort, L., & Carlsson, M. 2016, ApJ, 824, 65
  * Pereira et al. (2018) Pereira, T. M. D., Rouppe van der Voort, L., Hansteen, V. H., & De Pontieu, B. 2018, A&A, 611, L6
  * Peter et al. (2006) Peter, H., Gudiksen, B. V., & Nordlund, Å. 2006, ApJ, 638, 1086
  * Peter & Judge (1999) Peter, H. & Judge, P. G. 1999, ApJ, 522, 1148
  * Pneuman & Kopp (1977) Pneuman, G. W. & Kopp, R. A. 1977, A&A, 55, 305
  * Rouppe van der Voort et al. (2015) Rouppe van der Voort, L., De Pontieu, B., Pereira, T. M. D., Carlsson, M., & Hansteen, V. 2015, ApJ, 799, L3
  * Rouppe van der Voort et al. (2009) Rouppe van der Voort, L., Leenaarts, J., De Pontieu, B., Carlsson, M., & Vissers, G. 2009, ApJ, 705, 272
  * Rouppe van der Voort et al. (2020) Rouppe van der Voort, L. H. M., De Pontieu, B., Carlsson, M., et al. 2020, A&A, 641, A146
  * Rouppe van der Voort et al. (2007) Rouppe van der Voort, L. H. M., De Pontieu, B., Hansteen, V. H., Carlsson, M., & van Noort, M. 2007, ApJ, 660, L169
  * Rutten et al. (2019) Rutten, R. J., Rouppe van der Voort, L. H. M., & De Pontieu, B. 2019, A&A, 632, A96
  * Samanta et al. (2018) Samanta, T., Tian, H., & Prasad Choudhary, D. 2018, ApJ, 859, 158
  * Samanta et al. (2019) Samanta, T., Tian, H., Yurchyshyn, V., et al. 2019, Science, 366, 890
  * Scharmer (1981) Scharmer, G. B. 1981, ApJ, 249, 720
  * Scharmer (1984) Scharmer, G. B. 1984, Accurate solutions to non-LTE problems using approximate lambda operators, ed. W. Kalkofen, 173–210
  * Scharmer et al. (2003) Scharmer, G. B., Bjelksjö, K., Korhonen, T. K., Lindberg, B., & Petterson, B. 2003, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 4853, Innovative Telescopes and Instrumentation for Solar Astrophysics, ed. S. L. Keil & S. V. Avakyan, 341–350
  * Scharmer et al. (2019) Scharmer, G. B., Löfdahl, M. G., Sliepen, G., & de la Cruz Rodríguez, J. 2019, A&A, 626, A55
  * Scharmer et al. (2008) Scharmer, G. B., Narayan, G., Hillberg, T., et al. 2008, ApJ, 689, L69
  * Secchi (1871) Secchi, A. 1871, Atti dell’Accademia Pontificia de’ Nuovi Lincei. Anno 24., tomo 24. 1871
  * Secchi (1877) Secchi, A. 1877, L’astronomia in Roma nel pontificato DI Pio IX.
  * Sekse et al. (2012) Sekse, D. H., Rouppe van der Voort, L., & De Pontieu, B. 2012, ApJ, 752, 108
  * Sekse et al. (2013a) Sekse, D. H., Rouppe van der Voort, L., & De Pontieu, B. 2013a, ApJ, 764, 164
  * Sekse et al. (2013b) Sekse, D. H., Rouppe van der Voort, L., De Pontieu, B., & Scullion, E. 2013b, ApJ, 769, 44
  * Skogsrud et al. (2014) Skogsrud, H., Rouppe van der Voort, L., & De Pontieu, B. 2014, ApJ, 795, L23
  * Sterling (2000) Sterling, A. C. 2000, Sol. Phys., 196, 79
  * Tian et al. (2014) Tian, H., DeLuca, E. E., Cranmer, S. R., et al. 2014, Science, 346, 1255711
  * Tsiropoula et al. (2012) Tsiropoula, G., Tziotziou, K., Kontogiannis, I., et al. 2012, Space Sci. Rev., 169, 181
  * Uitenbroek (2003) Uitenbroek, H. 2003, ApJ, 592, 1225
  * van Noort et al. (2005) van Noort, M., Rouppe van der Voort, L., & Löfdahl, M. G. 2005, Sol. Phys., 228, 191
  * Vissers & Rouppe van der Voort (2012) Vissers, G. & Rouppe van der Voort, L. 2012, ApJ, 750, 22
  * Yurchyshyn et al. (2020) Yurchyshyn, V., Cao, W., Abramenko, V., Yang, X., & Cho, K.-S. 2020, ApJ, 891, L21
  * Zacharias et al. (2018) Zacharias, P., Hansteen, V. H., Leenaarts, J., Carlsson, M., & Gudiksen, B. V. 2018, A&A, 614, A110

##  Appendix A Spicule substructures

Figure 12: Discerning the substructures of the RBEs, RREs and, downflowing
RREs but only for the spicules harboring at least one of strong excursion RPs
(2, 3, 6 and 7). Panels (a) and (b) show the strong (RPs 2 and 3) and weak
(RPs 0 and 1) offset counterparts of the stronger excursion RBEs,
respectively; whereas panels (c) and (d) shows the strong (RPs 6 and 7) and
weak (RPs 4 and 5) offset counterparts belonging to the stronger
RREs/downflowing RREs. The color gradient indicates the density of the events
detected over the entire duration of the data. The black contour indicates the
regions that have an absolute value of the LOS magnetic field $\geq$ $100$ G.

We describe the substructures of spicules belonging to the stronger excursions
in blueward and redward side of H$\mathrm{\alpha}$ line core. We reiterate
that both the stronger and weaker excursion spicules can have spatial
substructures due to variation in their spectral properties which causes them
to have different RPs (stronger excursion RPs along with weaker excursion
RPs). As described in Sect 3.1.1 and 4.1, we considered those spicules as
stronger excursion that have at least one of stronger excursion RPs, i.e., RPs
2, 3, 6, and, 7, at least once during their entire lifetime. However, as
discussed earlier in the paper, even the strongest excursions have a spatial
variation in their Doppler offset along and perpendicular to the spicule axis,
hence weaker RPs, i.e., RPs 0, 1, 4 and, 5, could also be present within their
morphological structure (see Sect. 3.2.1). The major motivation of the
analysis presented in this section is to discern the substructures of strong
RBEs, RREs, and downflowing RREs depending on RPs and see if stronger and
weaker excursion RPs have any spatial preference. These variations add to the
structural complexity of spicules as shown in Fig. 12. Panel (b) shows the
”weaker” Doppler offset counterparts of the stronger blue excursion events in
panel (a). The arrangement is identical for the red excursions in panels (c)
and (d). Therefore, unlike Sect. 4.2, only the events common to the stronger
excursion category are identified and represented in panels (b) and (d) of
Fig. 12, while all others are excluded. Consequently, Fig. 9 shows the
distribution of spicules belonging uniquely to the strong and weak excursion
categories whereas Fig. 12 shows the spatial distribution of the spicules
belonging only to the stronger excursion category, but with different
excursion RPs. Naturally, the events in the different velocity compartments
are now no longer unique.

Upon closer investigation, we find that the number density of events increases
roughly by 45% for the strong blue excursions (panels (a)–(b)) and 27% (panels
(c)–(d)) for the strong red excursions. Furthermore, this increase is
predominantly reflected in the regions away from the strong network areas in
the former, whereas for the red excursions this is also linked to an increase
within and inside the strong field contours such as within the network patch
located at ($X$,$Y$) =(45″,33″) or close to the bottom of the FOV at $X$=20″in
panel (d). The increase in the number density can be attributed to the
variation in the Doppler offsets, and line widths of spicules in different
parts of their body.

##  Appendix B Supplementary figures

This section provides supplementary resources for the results presented in
Sects. 3.2, 4.1 and 4.3. The details of the figures shown in this section are
discussed in the main text.

Figure 13: Overview of the location and density distribution of the total
detected events over the whole FOV. Panel (a) shows the occurrence of
RREs/downflowing RREs over the whole time series against an H$\mathrm{\alpha}$
red-wing image at 93 km s-1; and panel (b) shows the occurrence of RBEs
against a background of CHROMIS continuum image at 4000 Å. Figure 14: Extended
$\lambda t$ diagrams corresponding to the two downflowing RREs shown in Fig. 7
displaying the lack of any blue shift preceding the downflowing RREs. The
maximum excursion towards the red-wing of H$\mathrm{\alpha}$ is indicated by
blue horizontal markers which also corresponds to the time shown in the
$\lambda t$ diagram of Fig. 7. Figure 15: Apparent motion of spicules in the
same format as Fig. 8 but for the full FOV. The behavior explained in the text
is also seen to be followed all over the FOV. Displacements smaller than
1″have been removed in order to improve the visibility.  Figure 16: Bi-variate
joint probability density between lifetime and eccentricity (of the fitted
ellipses) of (a) RREs/downflowing RREs and (b) RBEs in rainbow colormap as in
Fig. 11. The magenta contour indicates the region within which 70 % of the
events lie.