.epspdf.pdfconvert #1

11institutetext: Aberystwyth University, Ceredigion, Cymru, SY23 3BZ, UK
11email<EMAIL_ADDRESS>22institutetext: Centre for Fusion, Space and
Astrophysics, Department of Physics, University of Warwick, Coventry CV4 7AL,
UK 22email<EMAIL_ADDRESS>
# Automated analysis of oscillations in coronal bright points

B. Ramsey 1122 E. Verwichte 22 H. Morgan 11

(Received date / Accepted date)

###### Abstract

Context. Coronal bright points (BPs) are numerous, bright, small-scale
dynamical features found in the solar corona. Bright points have been observed
to exhibit intensity oscillations across a wide range of periodicities and are
likely an important signature of plasma heating and/or transport mechanisms.

Aims. We present a novel and efficient wavelet-based method that automatically
detects and tracks the intensity evolution of BPs using images from the
Atmospheric Imaging Assembly (AIA) on board the Solar Dynamics Observatory
(SDO) in the 193Å bandpass. Through the study of a large, statistically
significant set of BPs, we attempt to place constraints on the underlying
physical mechanisms.

Methods. We used a continuous wavelet transform (CWT) in 2D to detect the BPs
within images. One-dimensional CWTs were used to analyse the individual BP
time series to detect significant periodicities.

Results. We find significant periodicity at 4, 8-10, 17, 28, and 65 minutes.
Bright point lifetimes are shown to follow a power law with exponent $-1.13\pm
0.07$. The relationship between the BP lifetime and maximum diameter similarly
follows a power law with exponent $0.129\pm 0.011$.

Conclusions. Our wavelet-based method successfully detects and extracts BPs
and analyses their intensity oscillations. Future work will expand upon these
methods, using larger datasets and simultaneous multi-instrument observations.

###### Key Words.:

Sun: corona - Sun: oscillations - Sun: atmosphere

##  1 Introduction

Coronal bright points (BPs) are ubiquitous forms of activity seen as areas of
point-like emission in extreme-ultraviolet (EUV) to X-ray wavelengths, in both
quiet-Sun and coronal hole regions of the solar corona (Madjarska 2019). They
have been the subject of intense interest since their discovery in the 1970s
(Vaiana et al. 1973; Sheeley & Golub 1979). A BP is a collection of small
coronal loops (a mini active region) that forms an area of diffuse emission
10-60 arcsecs in size with a bright core of about 5 to 10 arcsecs (Golub et
al. 1977; Hirzberger et al. 2008). Bright points are associated with small
magnetic bipolar regions with typical photospheric magnetic fluxes of
$10^{-19}$–$10^{-20}\mathrm{Mx}$ (Golub et al. 1976).

On average, there are 400–800 BPs per day on the solar disk (Sattarov et al.
2002; McIntosh & Gurman 2005; Alipour & Safari 2015). The BP frequency outside
the active region belt does not show much variation with the solar cycle,
though the number decreases with temperature since BP production occurs at
temperatures well below the temperatures that soft X-ray detectors are
sensitive to (Hara & Nakakubo-Morimoto 2003; McIntosh & Gurman 2005). The
lifetimes of BPs seen in X-ray are found to exhibit a statistical
distribution, with a mean value of around 8 hours (Golub et al. 1974). Alipour
& Safari (2015) show that BP size approximately follows a power-law
distribution with exponent 0.14 with respect to the lifetime. Bright coronal
features with smaller scales, around $4\mathrm{Mm}$, and lifetimes of less
than 1 hour have been observed. These features are more commonly known as
brightenings, or transient coronal brightenings (Chen et al. 2021; Berghmans
et al. 2021), and are not considered BPs in this work.

Bright points have been observed to host intensity oscillations in both X-ray
and EUV with a broad range of periodicities, from a few minutes to hours
(Sheeley & Golub 1979; Christensen-Dalsgaard & Frandsen 1983; Strong et al.
1992; Ugarte-Urra et al. 2004; Kariyappa & Varghese 2008; Tian et al. 2008;
Kumar et al. 2011; Chandrashekhar et al. 2013). Early observations by Sheeley
& Golub (1979) showed a morphological evolution of about 6 minutes. Long-
period oscillations from 8 to 64 minutes were observed by Tian et al. (2008).
Longer-period oscillations were found by Kariyappa & Varghese (2008) in X-ray
BPs as seen by Hinode using power spectra analysis, with periods between 9 and
133 minutes. Wavelet analysis of time series by Ugarte-Urra et al. (2004)
showed BP oscillation periods as short as 236 seconds, with dominant
periodicities at 8 and 13 minutes. It is unknown if these oscillations are a
result of propagating magneto-acoustic waves or recurrent magnetic
reconnection.

In addition to intensity oscillations, decayless kink oscillations have also
been observed in BPs. These physical oscillations have been seen with periods
between 1 and 8 minutes, with an average of 5 (Gao et al. 2022).

Intensity oscillations have also been observed in other coronal structures.
Longer-period oscillations of 8–27 hours have been detected in coronal
filaments (Foullon, C. et al. 2004). Auchère et al. (2014) detected long-
period intensity oscillations of 3–16 hours. These oscillations have been seen
across active regions, visually associated with coronal loops, and in the
quiet Sun.

The physical cause of these long-period oscillations, which can last for
several days, is still uncertain. Numerical thermal non-equilibrium models by
Müller et al. (2005), used to explain coronal rain, can produce periods within
the range observed by Auchère et al. (2014). Thermal non-equilibrium has been
suggested as a proposed mechanism for coronal loop formation (Mok et al.
2008). Froment et al. (2015) show that intensity oscillations in loops are
linked to loop heating, with evaporation and condensation cycles, with
simulations to support this (Froment et al. 2017).

Various physical mechanisms have been proposed to explain intensity
oscillations, such as ‘leaky p-modes’ propagating along inclined magnetic flux
tubes (De Pontieu et al. 2005), standing waves within chromospheric cavities
(Leibacher & Stein 1981), recurrent small-scale reconnection events (Parker
1988; Chandrashekhar et al. 2013; Chandrashekhar & Sarkar 2015), and cyclic
loop heating (Habbal & Withbroe 1981; Verwichte & Kohutova 2017). The study of
acoustic waves provides physical insight into the mechanism of deposition
energy transport in the solar atmosphere (e.g. De Moortel & Hood 2003; Wang et
al. 2003). The simpler geometry of BPs compared to active regions makes it
easier to disentangle the observed signatures.

Figure 1:  Mexican hat wavelet applied to a 2D Gaussian with increasing
scales. The left panel shows a 2D Gaussian shape with width $\sigma_{0}=3$.
The remaining panels show, from left to right, the application of a Mexican
hat CWT with scales of $a$=1, 3, and 6, respectively. The solid curve shows
the 1D profile along a direction intersecting with the position of the
Gaussian.

Figure 2: Profile of the CWT ($n$=2) at the location of a 2D Gaussian of
amplitude $A$=1 and width $\sigma_{0}$=3, as a function of scale, $a$.

Many methods have been developed to automate the detection and
characterisation of a statistically meaningful set of BPs. Brajša et al.
(2001) used the regions of interest segmentation package in Interactive Data
Language (IDL). Hara & Nakakubo-Morimoto (2003) identified BPs in X-ray images
through thresholding based on the estimated noise level and other criteria,
such as the size and shape of the candidate regions. McIntosh & Gurman (2005)
developed a filtering and thresholding technique. Humphries et al. (2021)
developed a method for detecting and characterising small-scale brightenings
in EUV imagery using spatio-temporal bandpass filtering and adaptive
thresholding. Alipour & Safari (2015) used an automated machine learning
method, developed by Alipour et al. (2012), to study the statistical
properties of BPs. This method uses a training dataset of coronal BPs that is
collated by eye. Coronal BPs between 2 Mm and 20 Mm — co-located in
Atmospheric Imaging Assembly (AIA) 171Å, 193Å, and 211Å bands (Lemen et al.
2012) and detected by the Helioseismic and Magnetic Imager (HMI; Schou et al.
2012) on board the Solar Dynamics Observatory (SDO; Pesnell et al. 2012) —
were used to inform the machine learning algorithm. Alipour & Safari (2015)
implemented their method using data at a reduced 45-second cadence.

We wished to expand upon these investigations by applying automated analysis
techniques to determine the statistical properties of oscillations in
spatially resolved BPs. We focused on persistent BPs that have minimum
lifetimes of 1 hour (i.e. not the small, transient, brightenings studied by
Chen et al. 2021 and Berghmans et al. 2021). This allowed us to investigate
oscillations across four magnitudes of temporal scales. We analysed imaging
data from SDO/AIA, which has been collecting data almost continuously since
2010, at the full time cadence of 12s for the EUV channels of AIA. This
provides an unprecedented quantity of data, which allows for the systematic
analysis of BP oscillations across a large range of periods and across the
solar cycle.

The paper is structured as follows. Section 2 lays out the methodology of the
automated detection and tracking of BPs, with a focus on the following areas:
the acquisition of the data; detection using 2D continuous wavelet transforms
(CWTs); BP tracking; BP morphology; and BP time series analysis. The main
results and their implications are discussed in Sect. 3. Section 4 concludes
the work.

##  2 Automated analysis procedure

###  2.1 Data acquisition

We focused on analysing imaging data from AIA in the 193Å channel, which
covers the Fe XII and Fe XXIV emission lines and is sensitive to temperatures
around 1.5 and 20 MK. We chose this channel as a compromise between (i)
detecting many small-scale features, such as loop foot points, in lower
temperature bandpasses and only the hottest BPs in the hotter bandpasses and
(ii) maintaining a good signal-to-noise ratio. For this study, imaging data
from three days, 01 January 2020 to 03 January 2020 inclusive, were used at
the full time cadence of 12s. Level 1 data were used, that is, images have CCD
read–out noise (the noise of the on–chip amplifier) removed and rotated so as
to align with solar north. The images have a pixel size of 0.6 arcsec.

To find a suitable compromise between image detail and processing time during
the detection phase, images were reduced in resolution by a factor of 4, from
$4096\times 4096$ to $1024\times 1024$, using linear interpolation. The first
image in a time sequence was used as the initial reference image, and the
pointings of all subsequent images were aligned to it. This reduces offsets
and helps keep BPs centred within sub-images during tracking. Consistent image
headers were maintained. The near and off-limb areas of the image were
excluded from the region of interest (ROI) by use of a circular image mask
centred on the solar disk and with a radius equivalent to 0.7 $R_{\odot}$.

###  2.2 Detection using 2D continuous wavelets

For the automated detection of BPs in the images, we applied a 2D CWT to the
imaging data (Antoine et al. 2002; Hochedez et al. 2002; Delouille et al.
2005; Mallat 2008; White et al. 2012). The CWT of a 2D image, $I(\@vec{r})$,
is defined as

$\mathrm{CWT}(I)(\@vec{b},a,\theta)=\frac{1}{a^{n}}\iint\limits_{-\infty}^{+\infty}\\!I(\@vec{r})\,\,\psi\left(\frac{1}{a}R_{-\theta}(\@vec{r}-\@vec{b})\right)\mathrm{d}^{2}\@vec{r}\,\,,$
(1)

where $\psi\@vec{(r)}$ is called the mother wavelet. In the transform it is
translated by a 2D displacement vector, $\@vec{b}$, scaled by the finite scale
parameter $a$, and rotated through the rotation matrix $R_{\theta}$ at an
angle $\theta$ (Wang & Lu 2010). The transform can be defined depending on the
chosen norm, which is expressed through the index $n$. For the typical
$\text{L}_{1}$, norm $n$=2, which we use throughout this paper. In an
$\text{L}_{1}$ norm, magnitudes of vectors are calculated as the sum of the
absolute values of their complex components. The CWT corresponds to a scale-
space convolution of an image with a mother wavelet.

Figure 3: Full disk AIA and CWT images illustrating the detection of BPs.
Panel (a): Full disk AIA 193Å on 01 January 2020 at 01:00:05 UT. Panel (b): 2D
CWT at the active region scale, $a_{\mathrm{AR}}$. The active region is
detected by applying a threshold value at the 97th percentile of the CWT
value. This threshold area is the area within the red contour. Panel (c): 2D
CWT at the BP scale, $a_{\mathrm{BP}}$. The active region is masked by
removing the area denoted by the red contour, as in panel (b). The larger
black circle is applied at 0.7$R_{\odot}$, obscuring limb, off-limb, and edge
effects. In this case, the active region mask is outside the limb mask. Panel
(d): Full disk AIA 193Å on 01 January 2020 at 01:00:05 UT. Candidate BPs are
shown as white crosses. Masked areas are denoted by red and black contours.

Figure 4: Sub-image of BP #149 and the corresponding CWT image. Panel (a): AIA
193Å example BP #149 on 01 January 2020 at 01:59:59 UT. Panel (b): Application
of a CWT at scale $a_{\mathrm{BP}}$. White crosses highlight the maxima in the
CWT, which exceed the threshold value after the application of a 2D weighted
Gaussian to the CWT image. The red cross highlights the maxima closest to the
centre of the image.

Figure 5: Illustration of the process of extracting a BP region in four steps.
Panel (a): AIA 193Å example BP #149 on 01 January 2020 at 01:59:59 UT. Panel
(b): Binary mask showing the region that exceeds the median of the whole
weighted Gaussian data cube by 6 sigmas. Panel (c): Binary mask showing the
previous region but expanded to include pixels down to 3 sigmas. Panel (d):
Binary mask with holes filled using morphological filtering. This is the final
BP area. The temporal evolution of the BP in panel (a) is available as an
online movie.

The $\psi\@vec{(r)}$ is a wavelet if it is localised both in space and in its
reciprocal space, that is, $\iint\\!\psi(\@vec{r})\,\mathrm{d}^{2}\@vec{r}=0$,
and if it fulfils the admissibility condition in reciprocal space
$\iint|\hat{\psi}(\@vec{k})|^{2}\mathrm{d}^{2}\@vec{k}/k^{2}<0$, where
$\hat{\psi}$ and $\@vec{k}$ are the Fourier-transform of $\psi$ and the 2D
wave vector, respectively. The Mexican hat wavelet, which is proportional to
the Laplacian of a 2D Gaussian profile, satisfies these conditions:

$\psi(\@vec{r})=-\@vec{\nabla}^{2}\exp\left(-\frac{|\@vec{r}|^{2}}{2}\right)=(2-|\@vec{r}|^{2})\exp\left(-\frac{|\@vec{r}|^{2}}{2}\right)\,\,.$
(2)

The Mexican hat wavelet is isotropic, so the transform is independent of
rotation angle $\theta$. To understand how this 2D CWT can be employed for the
detection of BPs, we calculated the transform of a 2D Gaussian profile,
$I(\@vec{r})=A\,\exp(-(\@vec{r}-\@vec{c})^{2}/2\sigma_{0}^{2})$, where $A$ is
the amplitude and $\sigma_{0}$ is the width, centred on position $\@vec{c}$:

$\displaystyle\mathrm{CWT}(I)(\@vec{b},a)$ $\displaystyle=$
$\displaystyle\frac{1}{a^{2}}\iint\limits_{-\infty}^{+\infty}A\exp\left(-\frac{|\@vec{r}\\!-\\!\@vec{c}|^{2}}{2\sigma_{0}^{2}}\right)\,\psi\left(\frac{\@vec{r}\\!-\\!\@vec{b}}{a}\right)\,\mathrm{d}^{2}\@vec{r}\,\,,$
(3) $\displaystyle=$ $\displaystyle\frac{2\pi
A\left(\frac{a}{\sigma_{0}}\right)^{2}}{\left(1+\left(\frac{a}{\sigma_{0}}\right)^{2}\right)^{2}}\,\,\psi\left(\frac{\@vec{b}\\!-\\!\@vec{c}}{\sqrt{a^{2}\\!+\\!\sigma_{0}^{2}}}\right)\,\,.$

The result is a Mexican hat wavelet located at $\@vec{c}$ with a width and
amplitude that depend explicitly on scale $a$ (illustrated in Fig. 1). The
transform is maximal at $\@vec{b}=\@vec{c}$. At that position, the transform
has a clear profile as a function of $a$ with a maximum value of $\pi A$ at
$a=\sigma_{0}$. For scales smaller or larger than $\sigma_{0}$, the transform
drops away sharply. This is illustrated in Fig. 2. Therefore, this wavelet can
be used to differentiate point-like features across scales. Furthermore, the
Mexican hat CWT of a constant or linear trend in intensity is zero due to its
characteristic as a Laplacian, that is,
$\mathrm{CWT}(Axy\\!+\\!Bx\\!+\\!Cy\\!+\\!D)=0$. With the removal of
background signals, the CWT enhances image contrast (White et al. 2012). All
these characteristics make this wavelet particularly well suited for the
detection of point-like intensity features in solar images

The CWT was applied in two stages, illustrated in Fig. 3. First, we applied a
simple circular mask to the image with a diameter that is 1% shorter than that
of the Sun. This obscures the bight limb detail. Then we performed a 2D CWT on
that image with a scale of approximately $a_{\mathrm{AR}}$=90 arcsec. This
scale is much larger than the typical BP, and of the order of active regions.
Therefore, regions of high CWT intensity at that scale correspond to active
regions. As it is easy for an automated algorithm to confuse BPs and the foot
points of active region loops, we excluded active regions from the ROI by
subtracting the detected regions from the image mask. An active region was
designated where the intensity of the real part of the CWT is greater than the
97th percentile of the intensity. A mask of this defined area was then
applied.

Secondly, a 2D CWT was performed on the image with the active regions and limb
removed. Another simple circular mask with a diameter of 70% of the Sun-disk
radius was applied. This prevents detections close to the limb that may be
affected by the line of sight reduction and geometric projection. This means
that BPs close to the $0.7\text{R}_{\odot}$ boundary appear approximately 30%
smaller than at disk centre. The 2D CWT was applied at the scale of
$a_{\mathrm{BP}}$=7 arcsec, the typical scale of BPs (Golub et al. 1977;
Hirzberger et al. 2008). Candidate BPs were initially detected within the
second circular mask as regions at the 99th percentile of the real part of the
CWT.

###  2.3 BP tracking

Figure 6: Histogram showing the distribution of lifetimes. A power-law fitting
with exponent $-1.13\pm 0.07$ is shown in red. The grey background shows the
$2\sigma$ confidence level

Figure 7: Histogram showing the distribution of the average BP diameter. We
show a Gaussian fit with $\mu=24.06\pm 0.19,\leavevmode\nobreak\
\sigma=4.93\pm 0.19$, and $fwhm=11.60\pm 0.45.$

Figure 8: Lifetime versus maximum BP diameter. A total of 656 BPs are used. A
power law with the form $y=\alpha{}\tau^{\beta}{}$, where $\alpha{}=7.72\pm
0.87$ and $\beta{}=0.129\pm 0.011,$ is shown in red, and the grey background
shows the $3\sigma{}$ confidence level. The root–mean–square error for the fit
is 9.7. The average maximum BP diameter for each lifetime bin is shown as a
white circle, with the standard deviation as the error.

For each candidate BP ROI, properties such as location, total intensity, semi-
minor and semi-major radii, elliptical shape, and orientation were extracted.
First, we used this to further eliminate false detections. Regions with an
eccentricity greater than 0.6 were removed from the ROI, thereby eliminating
likely elongated loop structures and other structures that do not fit the
general morphological shape of a BP, although this does not prevent a
candidate BP from changing shape over its lifetime. Furthermore, regions with
total pixel areas of less than $30\leavevmode\nobreak\ \text{arcsec}^{2}$ were
also removed as they could be short-lived, small-scale solar transients or
cosmic rays. The above selection steps reduced the number of detections by
about 50%. We call the remainder of the detections simply BPs.

In order to ascertain the time range over which a BP is visible, a detection
procedure was repeated at 1-hour time intervals. At each hour, the found BPs
were compared with those found an hour earlier. This was to distinguish
between newly formed and pre-existing BPs and was achieved in the following
way. A binary image of the current image’s detections was subtracted from a
binary image of the previous image’s detections. This created a third image
with the following characteristics. Newly formed BPs are designated with a
value of -1, pre-existing BPs a value of 0, and BPs that disappear a value of
1. We could use these values to extract the newly formed BPs, for which we
could use the average position coordinates to generate a sub-image centred on
BP, with a size approximately eight times larger than the size of the detected
BP at its birth hour in pixels. This was a compromise intended to avoid the
computational expense of an unnecessarily large sub-image but include the
largest BP size found in the literature (Golub et al. 1977; Madjarska 2019),
approximately 60 arcsec. This sub-image sizing allows the BP to grow over its
lifetime and remain within the sub-image. The corresponding heliographic
coordinates were then used to track the BP position with time by rotating them
according to the local synodic solar differential rotation rate.

We then used a detection procedure and similarity test to determine the last
image in which a BP exists, as follows. The 2D CWT was applied to the first
sub-image at a scale length of 7 arcsecs. As with the initial detection, a
threshold mask was applied to the CWT image, initially at the 95th percentile.
If more than one area was detected within the CWT image, the process was
repeated in increasing 0.01% percentile increments until only one area
remained. This area should be the brightest point of the BP sub-image, but it
might not be in the centre of the sub-image; therefore, the sub-image was re-
centred on this point and the CWT re-applied. This process was then repeated
on the sub-image 1 hour later, using the coordinates of the re-centred BP.

Two tests were performed to determine if there is a BP present in a sub-image.
First, we determined the difference between the CWT value at the centre of the
real CWT image and the minimum CWT value in the image. If this value is small,
less than 100, then the CWT value at the centre and the minimum are close
together and therefore a BP is unlikely to be in the centre of that image.
Next, the centre value of the new BP sub-image and the standard deviation of
the average of the previous and current BP sub-image were found. If the
quotient of these two numbers was greater than 5 and the first criterion was
met, then a BP is said to exist. If either of these conditions was not met,
then the BP has disappeared and so we took the time of the last known image
containing the BP to be the BP’s death hour. This death hour was further
refined during the tracking process.

The birth hours were determined during the detection and comparison with
previous images. The birth hours for the BPs detected in the first image of
the full dataset are not known as they lie outside the bounds of the dataset
and were eliminated from the statistics. With birth and death hours
established, we created for each BP a 3D data cube from the 193Å dataset, at
the full spatial and temporal resolution (0.6 arcsec per pixel and 12s), of a
restricted field of view of between 72x72 and 115x115 arcsec, at the
heliographic rotating coordinates centred on the BP.

###  2.4 BP morphology

Figure 9: 1D CWT applied to BP#149. The top panel shows the time series of the
average BP intensity normalised by the standard deviation of the time series.
The bottom-left panel shows the wavelet power, with the COI in red and the
global confidence shown within the white contours. The bottom-right panel
shows the normalised Fourier spectrum in grey, the global wavelet spectrum in
black, the global significance level in red, and the local significance level
in orange. The noise model components are shown as follows: power law in
dashed orange, the kappa function in dashed blue, and the white noise in
dashed green.

Figure 10: From top to bottom: 1D CWT plots for BPs #195, #685, #735, #840,
and #1375. The left panels show the average BP intensity time series, with the
wavelet power below. The right panels show the Fourier spectrum in grey and
the global wavelet power spectrum in black. The solid red and orange lines
show the global and local wavelet significance levels, respectively. The noise
model components are shown as follows: power law in dashed orange, the kappa
function in dashed blue, and the white noise in dashed green. Figure 11:
Illustration of the process of bandpass-filtering the BP region. The leftmost
panel shows the AIA 193Å image at the beginning of the time range of interest.
The remaining three panels show images at three times separated by 2 minutes.
At each pixel, the time series has been bandpass-filtered around the period of
4 minutes with a Hann filter with a typical width equal to the mode frequency.
The temporal evolution of the bandpass-filtered images is available as an
online movie Figure 12: Illustration of the process of bandpass-filtering the
BP region. The leftmost panel shows the AIA 193Å image at the beginning of the
time range of interest. The remaining three panels show images at three times
separated by 4 minutes. At each pixel, the time series has been bandpass-
filtered around the period of 8 minutes with a Hann filter with a typical
width equal to the mode frequency. The temporal evolution of the bandpass-
filtered images is available as an online movie Figure 13: Illustration of the
process of bandpass-filtering the BP region. The leftmost panel shows the AIA
193Å image at the beginning of the time range of interest. The remaining three
panels show images at three times separated by 8 minutes. At each pixel, the
time series has been bandpass-filtered around the period of 16 minutes with a
Hann filter with a typical width equal to the mode frequency. The temporal
evolution of the bandpass-filtered images is available as an online movie

Figure 14: Normalised and weighted by BP lifetime, the total significant power
for oscillations that last for at least three periods for all BPs whose
lifetime is greater than 1 hour (in black). Blue shows BPs with lifetimes of 1
hour; the power has been scaled down by a factor of 7.

For each BP and at each time step, we extracted the relevant morphological
characteristics. To achieve this, we needed to more closely identify which
portion of the sub-image is identified as being part of the BP. We again
applied the 2D CWT to the image with a scale $a$=7 arcsecs. We then further
weighted the CWT signal by multiplying it by a 2D Gaussian with unity
amplitude centred on the field of view. We then identified the BP maximum as
the maximum in the image that is within the 95th percentile of the total
maximum and closest to the centre of the field of view (see Fig. 4). The 2D
Gaussian weight was then centred on the location of the found BP maximum. The
BP itself was then identified as the region overlapping the maximum and that
exceeds the sub-image median by 6$\sigma$. The region is represented by a
binary mask. This region was grown to encompass all neighbouring pixels down
to $3\sigma$. Any remaining holes in the found region were closed using
morphological filtering. This BP extraction method is illustrated in Fig. 5.
From the binary mask and original 193Å sub-image, statistical properties could
then be extracted for the BP, such as average and maximum intensity, size, and
shape. By repeating this procedure for all time steps, time series of BP
properties were created.

###  2.5 BP time series analysis

The method of extracting a BP from the background, as defined in Sect. 2.4,
may not successfully detect and extract the BP across the whole time series.
In order to prepare the morphological statistics of the BPs, some further
processing steps must be considered. In these instances, two actions were
performed on the morphological data prior to extracting statistical results.
First, instances in the time series were found wherein a BP is not extracted
consecutively for eight images; this is approximately 96 seconds and is about
half the period of the frequently observed 3-minute oscillations. Any
instances of missing data that last for longer than or equal to eight images
remove potentially important data; therefore, the time series were cut off at
the beginning of the gap. The remaining gaps in the data, fewer than eight
images, were then filled in using linear interpolation. Second, if the number
of remaining gaps equated to more than 15% of the number of data points, the
whole time series was dismissed. This removes the statistical unreliability
that comes with interpolating too many gaps in a dataset. Additionally, any BP
time series that has non-physical values, such as negative intensity or a
zero-valued area, were dismissed. Some time series show discontinuity in the
form of intensity jumps, which will result in the 1D CWT (Torrence & Compo
1998; Verwichte et al. 2004; Auchère et al. 2016) power being dominated by the
discontinuous jumps if present; therefore, before these time series were
analysed, a point filter was applied to smooth sporadic jumps. The point
filter identified outliers by comparing values to the local standard deviation
and replaced the value with a new value closer to the local mean. Lastly, each
time series was manually checked to ensure that any anomalous time series were
not included in the analysis results. Each time series of average BP intensity
was analysed using a 1D CWT and a custom noise model given by Auchère et al.
(2016) in the form $\sigma{}(v)=Av^{s}+BK_{\rho}{}(v)+C$, where the first term
represents the power-law dependence given by background stochastic
fluctuations. The second term is a kappa function that is related to pulses in
the time series. The final term corresponds to high-frequency white noise.

##  3 Results and discussion

This detection method found 3308 BPs over the three-day period: 1191 on 01
January 2020, 1141 on 02 January 2020, and 976 on 03 January 2020. The number
of BPs used in the analysis was eventually reduced to 656.

###  3.1 General BP characteristics

Figure 6 shows the distribution of BP lifetimes, with a mean lifetime of 6.8
hours and a range from 1 hour to 22 hours. We fitted a power law of the form
$y\/=\/ax^{\beta}$ with an exponent equal to $-1.13\pm 0.05$. The lifetime
shows an almost power-law distribution across all lifetimes, with a good fit
up to approximately 700 minutes, after which the function begins to diverge.
Previous work by McIntosh & Gurman (2005) found a power-law behaviour with
exponentials at longer lifetimes. The exponents for the power law by McIntosh
& Gurman (2005) vary with temperature and time over the solar cycle, but fall
between -1.24 and -2.00 for the 195Å channel of the Extreme ultraviolet
Imaging Telescope (EIT) on board the Solar and Heliospheric Imager (SOHO).
Alipour & Safari (2015) find an exponent of -1.6. These previous exponents
were found at varying stages of the solar cycle, and, considering our very
small dataset, we cannot comment on the significance of our exponent in
relation to the solar cycle. Our exponent is consistently different from
literature values, perhaps due to our power-law-only fit with a much shorter
maximum lifetime. McIntosh & Gurman (2005), for example, show EUV BPs with
lifetimes in excess of 100 hours. In contrast, the maximum lifetime we find is
22 hours. This is most likely the source of the discrepancy between the power-
law fits. Additionally, short-lifetime detections are not considered BPs
within our data as BPs with lifetimes of less than an hour and scales of
around 4 Mm could be considered coronal brightenings (Chen et al. 2021;
Berghmans et al. 2021).

The average BP diameter was obtained by taking twice the mean semi-major
radius for each BP across its lifetime. Similarly, the maximum BP diameter is
defined as twice the maximum semi-major radius that a BP achieves across its
lifetime. The mean diameter is $24.06\pm 0.19$ arcsec, with a normal
distribution of width $\sigma=4.93\pm 0.19$ and a range from 10 arcsec to 39
arcsec, as shown in Fig. 7. When we compared the mean diameter against a BP’s
lifetime, we find no clear relationship. However, there is a relationship
between the maximum diameter and lifetime, as shown in Fig. 8. The
distribution in this figure shows a general increase in the diameter with
lifetime. We binned the lifetimes in one-hour increments and calculated the
mean and standard deviation of the maximum diameters for each bin. These are
illustrated in Fig. 8 and show an increase from the shortest to longest
lifetimes. At the very longest lifetimes, we see a decrease in the average
maximum diameter, but there are only a few BPs with such long lifetimes, so
the statistics are not reliable. To further characterise this relationship, we
fitted the maximum diameters to a power law of the lifetimes. We note that
this fit was done with all the individual BP data points, not the bin
averages. Our power-law fit has the form
$D_{\mathrm{max}}=\alpha\tau^{\beta}$, where

$D_{\mathrm{max}}\,=\,(7.72\pm 0.87)\,\tau^{0.129\pm 0.011}\,\,.$ (4)

Here, $D_{\mathrm{max}}$ and $\tau$ are the maximum BP diameter in mega-metres
and lifetime in seconds, respectively. The parameter errors here were obtained
as the standard deviation of the non-linear least mean square fitting method
used to fit the data, with an RMS error of 9.7. This power law confirms the
clear relationship between the maximum diameter and lifetime. That there is no
clear relationship between mean diameter and lifetime is surprising,
particularly given the relationship shown by Fig. 11 of Alipour & Safari
(2015). Their study, however, looks at considerably smaller lifetime and
spatial scales than ours. The lack of a relationship between lifetime and mean
diameter and the clear relationship between lifetime and maximum diameter is
interesting and requires further study. Our detection method does have a bias,
in that we discard BPs with some minimum diameter ( 6.2”) and minimum lifetime
(1 hour). Furthermore, at long lifetimes we have far fewer BPs. A larger study
would help with the statistics at longer lifetimes.

###  3.2 Example BPs

Using BPs #149, #195, #685, #735, #840, and #1378 as examples, we applied 1D
CWTs to the average BP intensity, as shown in Fig. 9. We can see, in the
bottom-left panel, several periodicities with significant power appearing
regularly across the lifetime of the BP, namely at periods between 1 and 10
minutes. We can see, in the right panel, a highly structured Fourier spectrum
shown along with the wavelet power. We can see regions of wavelet power above
the significance levels; these regions are more concentrated at the beginning
of the BP’s lifetime, especially at shorter periods.

The global and local significance levels are shown in Fig. 9. The global
significance level represents the wavelet power that lies above a global
confidence level when compared with the noise model. The local significance
level represents the probability that power in a single bin is significant
when compared to the noise model. Auchère et al. (2016) provide a detailed
explanation of the noise model and accompanying significance levels.

At shorter periods of between 1 and 10 minutes the wavelet power is above both
the local and global significance levels. This is in contrast to longer
periods of approximately 30 and 70 minutes, where the wavelet power is only
above the local significance level.

We can see further examples of these wavelet power spectra in Fig. 10. Here we
show the wavelet spectra and powers for BPs #195, #685, #735, #840, and #1378.
These BPs show a range of lifetimes, from 2 to over 16 hours. We can see some
common periodicities, for example for BPs #685 and #735, fairly distinct peaks
in both the Fourier and wavelet power spectra that lie above the local and
global significance level, at approximately 4 minutes. On the other hand, BP
#840 has a suggestion of periodicity at 4 minutes above the local significance
only. At longer periods we can start to see some common periodicity in BPs
#195, #685, and #735, at approximately 30-40 minutes.

The periods of interest seen in Fig. 9 are shown for BP #149 in Figs. 11, 12,
and 13. We chose to look at periods of 4, 8, and 16 minutes as they can be
seen in the main wavelet plot as areas above the significance level. We
applied a temporal bandpass filter around the period of interest with a Hann
filter with a typical width equal to the mode frequency. This shows a change
over the oscillation cycle, with the spatial structure of the three
oscillation periods showing clear differences. For the 4-minute oscillation in
Fig. 11, the BP shows hints of anti-phase behaviour between its two sides,
which could be an m=1 mode. Such modes have also been observed in sunspots
(Jess et al. 2017) and chromospheric vortices (Murabito et al. 2020). While
the structure of a BP is significantly different from that of a sunspot or
chromospheric vortex, this behaviour could suggest leaky p modes on the
eastern and western sides, which would constitute an apparent m=1 mode
structure. The 8- and 16-minute oscillations show a phase coherence that maps
out the coronal loops in the BP. The acoustic cut-off frequencies in the quiet
Sun vary in the range 4-6 mHz (3-4 minutes; Felipe & Sangeetha 2020).
Intensity oscillations with similar or shorter periods are interpreted to be
acoustic in nature. Those with periods substantially longer, 10 minutes and
more, are unlikely to be acoustic waves propagating up from below and instead
may be evanescent acoustic tails or be associated with thermal limit cycles
(e.g. Foullon, C. et al. (2004); Auchère et al. (2014); Froment et al. (2017);
Verwichte & Kohutova (2017))

###  3.3 General periodicity

There are, more generally, many periodicities potentially detected across all
of the analysed BPs. To better visualise the significance of the detected
periodicities, we took the significant normalised wavelet power – that is, the
power within the cone of influence (COI; see Fig. 9) – as an example; it lasts
for at least three complete periods. These powers were then summed together
for all BPs and normalised. All BPs can contain periodicities of less than 1
hour; conversely, periodicities greater than 1 hour are not possible for BPs
whose lifetimes are shorter than 1 hour. We therefore applied a weighting to
the total normalised power. This result can be seen in Fig. 14, where we show
the normalised and weighted total significant power for BPs with lifetimes
greater or less than 1 hour. Combining these BPs results in a discontinuity at
10 minutes as this is the longest periodicity that the wavelet can detect;
periods greater than 10 minutes fall outside the COI for BPs with lifetimes of
1 hour.

We can see from these two plots a peak below and above 4 minutes (average
time, 4.11 minutes). For BPs of 1 hour, we see another peak between 8 and 9
minutes, whereas for the remaining BPs we see this peak at 10 minutes. There
is a final noticeable peak at approximately 17 minutes. For longer
periodicities, there are slight humps within the plot at approximately 28, 40,
and 49 minutes. At longer periods we see hints of additional periodicity at
approximately 65 to 75 minutes.

Ugarte-Urra et al. (2004) showed dominant periods of between 8 and 13 minutes
in BP oscillations. Our 10-minute peak falls within this range; its physical
nature is still uncertain it but could be evanescent acoustic modes moving
through the transition regions resonant in loops or thermal-limit cycles
(Habbal & Withbroe 1981; Verwichte & Kohutova 2017). The peaks at 17, 28, and
65 minutes seen here have been seen by Tian et al. (2008) (16, 28, and 64
minutes). We do not see the clear peak at 32 minutes noted by Tian et al.
(2008); however, this potential peak in our case could be too broad and
therefore not discernable in Fig. 14. Zhang et al. (2012) found a periodicity
of about 1 hour, which they suggested was due to quasi-periodic recurrent
flashes. The peak at 17 minutes falls within the range of 15–25 minutes seen
by Chandrashekhar & Sarkar (2015) within their simulated loop and nanoflare
model, in addition to the observed values by Chandrashekhar et al. (2013). The
most dominant periods are the 4- and 10-minute ones seen in Fig. 14. The
4-minute period is around the 5-minute photospheric p-mode, which could
contribute to oscillations in the corona under certain circumstances, namely
the ‘leaky’ p-mode, (De Pontieu et al. 2005). Srivastava & Dwivedi (2010)
found periodicity of $241\pm 60s$ ($4.0\pm 1$ minutes). Our 4-minute period
falls within that range. Future work will focus on the multi-bandpass (multi-
thermal) and spatial mode structure aspects of the BP oscillations. A
comparison of various analysis methodologies employed by various studies may
also be undertaken.

##  4 Conclusion

The aim of this work was to analyse a large set of coronal BPs using CWTs in
2D and 1D. We present a novel method for the detection and tracking of BPs
using 2D CWTs. We analysed the morphology of these BPs and investigated
intensity oscillations using 1D CWTs.

We find that BPs have a lifetime distribution that follows a power law, with
exponent $-1.13\pm 0.07$. We find that the relationship between a BP’s
lifetime and maximum diameter roughly follows a power law with exponent
$0.129\pm 0.011$. These statistical results compare well with previous studies
of BPs (Alipour & Safari 2015; McIntosh & Gurman 2005). The analysis of
intensity oscillations within BPs shows a broad range of significant
periodicity between 1 and 100 minutes, with notable peaks at 4, 10, and 17
minutes and the suggestion of peaks at longer periods, namely 28 and 65
minutes. Further in-depth analysis is required to study the spatial mode
structure of these oscillations and place constraints on their physical
nature.

There is a clear relationship between a BP’s area and intensity. However, the
effects of limb projection are not considered here, and future work will
endeavour to address this. The hope is the automated methods described here
will allow for a much larger statistical study of BP intensity oscillations
and their morphological characteristics in the future. Additionally, further
work will endeavour to expand the automated methods into different SDO/AIA
passbands and different instruments.

###### Acknowledgements.

The wavelet transform has been performed using the Python wavelet module by
Erwin Verwichte (University of Warwick) and was supported by a UK STFC grant
ST/L006324/1. We acknowledge STFC studenship ST/V506527/1, and STFC grant
ST/S000518/1 to Aberystwyth University.

## References

  * Alipour & Safari (2015) Alipour, N. & Safari, H. 2015, ApJ, 807, 175
  * Alipour et al. (2012) Alipour, N., Safari, H., & Innes, D. E. 2012, The Astrophysical Journal, 746, 12
  * Antoine et al. (2002) Antoine, J. P., Demanet, J., Hochedez, J.-F., et al. 2002, Physicalia Magazine, 24, 93
  * Auchère et al. (2014) Auchère, F., Bocchialini, K., Solomon, J., & Tison, E. 2014, A&A, 563, A8
  * Auchère et al. (2016) Auchère, F., Froment, C., Bocchialini, K., Buchlin, E., & Solomon, J. 2016, The Astrophysical Journal, 825, 110
  * Berghmans et al. (2021) Berghmans, D., Auchère, F., Long, D. M., et al. 2021, A&A, 656, L4
  * Brajša et al. (2001) Brajša, R., Wöhl, H., Vršnak, B., et al. 2001, A&A, 374, 309
  * Chandrashekhar et al. (2013) Chandrashekhar, K., Krishna Prasad, S., Banerjee, D., Ravindra, B., & Seaton, D. B. 2013, Sol. Phys., 286, 125
  * Chandrashekhar & Sarkar (2015) Chandrashekhar, K. & Sarkar, A. 2015, ApJ, 810, 163
  * Chen et al. (2021) Chen, Y., Przybylski, D., Peter, H., et al. 2021, A&A, 656, L7
  * Christensen-Dalsgaard & Frandsen (1983) Christensen-Dalsgaard, J. & Frandsen, S. 1983, Sol. Phys., 82, 165
  * De Moortel & Hood (2003) De Moortel, I. & Hood, A. W. 2003, A&A, 408, 755
  * De Pontieu et al. (2005) De Pontieu, B., Erdélyi, R., & De Moortel, I. 2005, ApJ, 624, L61
  * Delouille et al. (2005) Delouille, V., De Patoul, J., Hochedez, J. F., Jacques, L., & Antoine, J. P. 2005, Sol. Phys., 228, 301
  * Felipe & Sangeetha (2020) Felipe, T. & Sangeetha, C. R. 2020, A&A, 640, A4
  * Foullon, C. et al. (2004) Foullon, C., Verwichte, E., & Nakariakov, V. M. 2004, A&A, 427, L5
  * Froment et al. (2017) Froment, C., Auchère, F., Aulanier, G., et al. 2017, ApJ, 835, 272
  * Froment et al. (2015) Froment, C., Auchère, F., Bocchialini, K., et al. 2015, ApJ, 807, 158
  * Gao et al. (2022) Gao, Y., Tian, H., Doorsselaere, T. V., & Chen, Y. 2022, The Astrophysical Journal, 930, 55
  * Golub et al. (1977) Golub, L., Krieger, A. S., Harvey, J. W., & Vaiana, G. S. 1977, Sol. Phys., 53, 111
  * Golub et al. (1974) Golub, L., Krieger, A. S., Silk, J. K., Timothy, A. F., & Vaiana, G. S. 1974, ApJ, 189, L93
  * Golub et al. (1976) Golub, L., Krieger, A. S., & Vaiana, G. S. 1976, Sol. Phys., 50, 311
  * Habbal & Withbroe (1981) Habbal, S. R. & Withbroe, G. L. 1981, Sol. Phys., 69, 77
  * Hara & Nakakubo-Morimoto (2003) Hara, H. & Nakakubo-Morimoto, K. 2003, ApJ, 589, 1062
  * Hirzberger et al. (2008) Hirzberger, J., Gizon, L., Solanki, S. K., & Duvall, T. L. 2008, Sol. Phys., 251, 417
  * Hochedez et al. (2002) Hochedez, J. F., Jacques, L., Verwichte, E., et al. 2002, in ESA Special Publication, Vol. 477, Solspa 2001, Proceedings of the Second Solar Cycle and Space Weather Euroconference, ed. H. Sawaya-Lacoste, 115–118
  * Humphries et al. (2021) Humphries, L. D., Morgan, H., & Kuridze, D. 2021, Sol. Phys., 296, 140
  * Jess et al. (2017) Jess, D. B., Van Doorsselaere, T., Verth, G., et al. 2017, ApJ, 842, 59
  * Kariyappa & Varghese (2008) Kariyappa, R. & Varghese, B. A. 2008, A&A, 485, 289
  * Kumar et al. (2011) Kumar, M., Srivastava, A. K., & Dwivedi, B. N. 2011, MNRAS, 415, 1419
  * Leibacher & Stein (1981) Leibacher, J. W. & Stein, R. F. 1981, Oscillations and pulsations., ed. S. Jordan, Vol. 450, 263–287
  * Lemen et al. (2012) Lemen, J. R., Title, A. M., Akin, D. J., et al. 2012, Sol. Phys., 275, 17
  * Madjarska (2019) Madjarska, M. S. 2019, Living Reviews in Solar Physics, 16, 2
  * Mallat (2008) Mallat, S. 2008, A Wavelet Tour of Signal Processing: The Sparse Way
  * McIntosh & Gurman (2005) McIntosh, S. W. & Gurman, J. B. 2005, Sol. Phys., 228, 285
  * Mok et al. (2008) Mok, Y., Mikić, Z., Lionello, R., & Linker, J. A. 2008, ApJ, 679, L161
  * Müller et al. (2005) Müller, D. A. N., De Groof, A., Hansteen, V. H., & Peter, H. 2005, A&A, 436, 1067
  * Murabito et al. (2020) Murabito, M., Shetye, J., Stangalini, M., et al. 2020, A&A, 639, A59
  * Parker (1988) Parker, E. N. 1988, ApJ, 330, 474
  * Pesnell et al. (2012) Pesnell, W. D., Thompson, B. J., & Chamberlin, P. C. 2012, Sol. Phys., 275, 3
  * Sattarov et al. (2002) Sattarov, I., Pevtsov, A. A., Hojaev, A. S., & Sherdonov, C. T. 2002, ApJ, 564, 1042
  * Schou et al. (2012) Schou, J., Scherrer, P. H., Bush, R. I., et al. 2012, Sol. Phys., 275, 229
  * Sheeley & Golub (1979) Sheeley, N. R., J. & Golub, L. 1979, Sol. Phys., 63, 119
  * Srivastava & Dwivedi (2010) Srivastava, A. K. & Dwivedi, B. N. 2010, Monthly Notices of the Royal Astronomical Society, 405, 2317
  * Strong et al. (1992) Strong, K. T., Harvey, K., Hirayama, T., et al. 1992, PASJ, 44, L161
  * Tian et al. (2008) Tian, H., Xia, L. D., & Li, S. 2008, A&A, 489, 741
  * Torrence & Compo (1998) Torrence, C. & Compo, G. P. 1998, Bulletin of the American Meteorological Society, 79, 61
  * Ugarte-Urra et al. (2004) Ugarte-Urra, I., Doyle, J. G., Nakariakov, V. M., & Foley, C. R. 2004, A&A, 425, 1083
  * Vaiana et al. (1973) Vaiana, G. S., Davis, J. M., Giacconi, R., et al. 1973, ApJ, 185, L47
  * Verwichte & Kohutova (2017) Verwichte, E. & Kohutova, P. 2017, A&A, 601, L2
  * Verwichte et al. (2004) Verwichte, E., Nakariakov, V. M., Ofman, L., & Deluca, E. E. 2004, Sol. Phys., 223, 77
  * Wang & Lu (2010) Wang, N. & Lu, C. 2010, Journal of Atmospheric and Oceanic Technology, 27, 652
  * Wang et al. (2003) Wang, T. J., Solanki, S. K., Innes, D. E., Curdt, W., & Marsch, E. 2003, A&A, 402, L17
  * White et al. (2012) White, R. S., Verwichte, E., & Foullon, C. 2012, A&A, 545, A129
  * Zhang et al. (2012) Zhang, Q. M., Chen, P. F., Guo, Y., Fang, C., & Ding, M. D. 2012, ApJ, 746, 19