QUASI-PERIODIC PULSATIONS IN SOLAR AND STELLAR FLARES: RE-EVALUATING THEIR NATURE IN THE CONTEXT OF POWER-LAW FLARE FOURIER SPECTRA

A common feature of solar flare emission is the appearance of quasi-periodic pulsations (QPPs). Alternately known—particularly in stellar and astrophysical contexts—as quasi-periodic oscillations (QPOs), these phenomena have been observed in a wide range of wavelengths over several decades (see Nakariakov & Melnikov 2009, for a recent review). Although not precisely defined in the literature, the term QPP is most often used to describe variations in the flux from a flare or other astrophysical object as a function of time, which appear to include periodic components with characteristic timescales ranging from one second up to several minutes, although shorter (Tan & Tan 2012) and longer (Foullon et al. 2010) timescales are sometimes studied. They are typically observed during the impulsive phase of solar flares and have been observed over a wide range of wavelengths, from radio waves and microwaves to hard X-rays and gamma-rays. Similar signatures have also been observed from stellar flares (e.g., Mathioudakis et al. 2006; Kowalski et al. 2010). Since QPPs are directly linked to the properties of the flare reconnection region and flare acceleration sites, a full description of QPPs remains crucial for our understanding of solar flares.

There are two main theories that are currently being pursued as possible mechanisms for generating QPPs (see Nakariakov & Melnikov 2009), both of which assume the presence of a periodic driver. These are that the observed flux variations are driven either by (1) magnetohydrodynamic (MHD) wave behavior in the corona and in flare sites, or are instead (2) a result of a regime of periodic or bursty reconnection. As a consequence, recent studies have focused on the assumed periodic nature of QPPs; such studies are usually motivated by the concept of searching for a periodic signal obscured by random noise and long-term trends in the flare signal.

However, it has recently become clear (e.g., Gruber et al. 2011; Vaughan 2010) that flare time series are often dominated by a power law in the Fourier domain, rather than random white noise. In an astrophysical context this is often referred to as "red noise." Objects that are known to exhibit time series with power-law-like behavior in the Fourier power spectrum include XMM-Newton observations of Seyfert galaxies (Vaughan 2005, 2010), gamma-ray bursts (Cenko et al. 2010), active galactic nuclei (McHardy et al. 2006), and magnetars (Huppenkothen et al. 2013). The Fourier power spectra of all of these objects are well-described using power-law models with a negative slope, i.e., P(f) ≈ f^−α for α ⩾ 0 where f is frequency. Strictly speaking, "red noise" refers to a specific slope of the Fourier power-law spectrum, α = 2, however, the term is often more loosely used to refer to any power-law-dominated Fourier power spectrum where α ⩾ 0. The term "noise" in this context is also something of a misnomer; the observed emission is not noise in the conventional sense of detector noise or measurement error. Rather, the power-law shape of the Fourier power spectrum is an intrinsic property of the physical system and must, therefore, be taken into account when examining for other effects such as oscillations.

That solar flares exhibit Fourier power-law-like properties has in fact been known for some time (e.g., Ryabov et al. 1997; Aschwanden et al. 1998; Schwarz et al. 1998). Such properties are clearly present in the Fourier power spectra of RHESSI X-ray flare observations. However, until recently the implications for QPPs have not been considered. Gruber et al. (2011) were the first to demonstrate the importance of considering this property of flaring emission in QPP studies by analyzing Fermi/Gamma-ray Burst Monitor (GBM) and RHESSI solar flare data from the same events in two different ways. They demonstrated that the consequence of applying a white-noise assumption to data dominated by a power law in the Fourier power spectrum is a drastic overestimation of the significance of peaks detected in the Fourier domain. Additionally, the nature of the Fourier power spectrum in such signals means that the empirical subtraction of "background components" of the signal may lead to misleading results and should be avoided; the entire Fourier spectrum of the signal should always be considered (Vaughan 2010).

Hence the results of a number of previous studies, and the general prevalence of flare signals with oscillatory content, are both called into question. However, as Gruber et al. (2011) also point out, their study does not mean that observed variations in flare signals are not real, and we reiterate that here. Rather, it means that the physical mechanisms causing the observed flux variations need not be exclusively due to an oscillatory mechanism.

In this paper, we address this problem and the question of the true prevalence of statistically significant narrowband oscillations in solar flares. We achieve this by adopting an approach based on Vaughan (2010), which takes into account the true statistical properties of flare time series. In this paper, we introduce the method and demonstrate its application to selected solar and stellar events where the presence of QPPs has been suggested.

For this study, we concentrate on a joint data set consisting of X-ray data from the GBM (Meegan et al. 2009) on board Fermi, and soft X-ray and EUV data from the Large Yield Radiometer (LYRA) instrument (Dominique et al. 2013) on board PROBA2. Fermi was launched in 2008 and GBM, which consists of 12 NaI and 2 BGO detectors, has provided high-cadence solar observations in the 5 keV–40 MeV energy range since then. Up to six of the NaI detectors and one BGO detector may be sunward facing at any given time.

PROBA2 was launched in 2009 with a science payload including the LYRA, which observes the Sun in four channels with a nominal time resolution of 0.05 s, although the instrument is capable of 0.01 s resolution. Of these channels, two continue to take science quality data, namely the aluminum (Al) filter (17–80 nm + SXR) and the zirconium (Zr) filter (6–20 nm + SXR). The availability of these simultaneous data sets provides us with the necessary multi-wavelength information for analyzing QPPs effectively in different regimes, as well as the time cadence required for the detection of pulsations in the 1–300 s range.

Two recent solar flares considered to be QPP events are the GOES-class X2.2 flare of 2011 February 15 (e.g., Dolla et al. 2012), and the GOES-class M2.5 flare of 2011 June 7 (e.g., Inglis & Gilbert 2013). Both of these events were fully observed by the LYRA and GBM instruments, making them ideal candidates for a multi-wavelength investigation of their power spectral properties.

For each of these events we apply the analysis method described in Section 4 to a selection of GBM data channels. For this analysis, we utilize the CTIME data product, which provides X-ray count information at a 0.256 s temporal cadence in the nominal observing mode, and ≈0.064 s cadence during flare times. GBM CTIME data is available in eight energy channels, covering the range 4 keV–2 MeV. For the analysis presented in this paper, we select the following intervals: Channel 1 (≈ 12–27 keV), Channel 2 (≈27–50 keV), and Channel 3 (≈50–100 keV). The lowest energy emission observed in Channel 0 (4–12 keV) is not used. These intervals are suited to capturing both the soft and hard X-ray emission from flares. To ensure a consistent observational cadence and to improve signal-to-noise, all the analyzed GBM data are re-binned into 1 s intervals prior to analysis. Similarly, for LYRA the data are accumulated into 1 s intervals throughout this paper. Subsequently, identical analysis is applied to data obtained from the GBM channels to capture hard X-ray behavior, and to the Al and Zr filters, capturing EUV and soft X-ray emission. Hence the result will be a thorough understanding of the Fourier spectral properties of these flares across EUV, soft X-ray and hard X-ray wavelengths.

Additionally, we return to the flare of 1998 May 8, previously investigated in the microwave and X-ray regimes by Stepanov et al. (2004) and Inglis et al. (2008), perhaps one of the most pronounced examples of QPPs in the prior literature. During this flare, pronounced pulsations were observed in the 17 GHz data from the Nobeyama Radioheliograph (NoRH; Nakajima et al. 1994), as well as in X-ray flux in the 13–93 keV range as observed by the Yohkoh satellite (Ogawara et al. 1991). In Stepanov et al. (2004) and Inglis et al. (2008) the observed pulsations were analyzed using a white-noise assumption, resulting in an observed 16 s period that was interpreted as a signature of either an MHD ballooning mode or a magnetoacoustic sausage mode. However, it is timely to revisit this event and re-analyze the NoRH radio data and the Yohkoh X-ray data in the context of power-law Fourier spectra.

Finally, for a stellar flare perspective, we analyze near-UV U-band (∼3200–3900 Å) data from the New Mexico State University (NMSU) 1m telescope (Holtzman et al. 2010) for the "megaflare" of 2009 January 16, previously studied by Kowalski et al. (2010), Kowalski et al. (2013), and Anfinogentov et al. (2013). This flare exhibited interesting temporal variations that were recently investigated by Anfinogentov et al. (2013), who suggested that the observed emission was the manifestation of a magnetoacoustic wave. However, as with the solar flares above, the nature of the Fourier power spectrum was not considered in the analysis of this event.

3.1. LYRA Large-angle Rotations (LARs)

In order to apply our analysis method to solar flare data, the GBM and LYRA data require some preparation. In particular, a design feature of PROBA2 is the incidence of large-angle rotations (LARs), in which the spacecraft is rotated approximately 90° in order to accommodate the on-board star tracker. These maneuvers are performed four times per orbit, or approximately every 25 minutes. In some cases, a consequence of this is the appearance of temporal artifacts in the LYRA data during LARs, due to the movement of the Sun within the instrument's field of view.

A list of LARs and other spacecraft events (for example, passes through the South Atlantic Anomaly) is maintained by the PROBA2 instrument team, known as the LYRA timeline annotation file (LYTAF). In this work, we remove the effects of LARs in the simplest way, by removing the affected times from the data entirely. This is achieved using routines within the SunPy (Mumford et al. 2013) data analysis package that query the LYTAF and, for a given data set, automatically separate an LYRA time series into components, excluding the affected portions of the data (see Figure 1).

Zoom In Zoom Out Reset image size

3.2. Data Normalization and Apodizing

An additional step in data preparation for all data sources is to normalize and apodize the input signal F. This normalization is done according to,

$$\begin{equation} F_{{\rm norm}} = \frac{F - \bar{F}}{\bar{F}}, \end{equation} \tag{ 1 }$$

where $\bar{F}$ is the mean of the signal. Finally, the normalized signal is prepared by multiplying it by a window function. This step is crucial, otherwise the Fourier power spectral density (PSD) of the data is dominated by the discontinuity caused by the mismatch between the beginning and end of the signal. For this purpose, we choose the well-established Hann window (Blackman & Tukey 1959). For verification, other window functions (e.g., Blackman–Harris, Hamming) were tested and found to produce very similar results; for brevity these tests are not shown in this paper.

In all of the results shown in Section 5, the analysis method has been applied to the prepared, apodized data as described here.

A powerful and flexible methodology designed for testing the presence of oscillations in the presence of frequency-dependent noise is described in Vaughan (2010), where it is applied to determine the presence of QPOs in XMM-Newton observations of highly variable Seyfert 1 galaxies, objects that also exhibit frequency-dependent Fourier power spectra. A key advantage of this approach is that we avoid techniques subtracting empirically defined background trends, which inevitably modify the frequency content of the analyzed data (Cenko et al. 2010; Auchère et al. 2014). Additionally, one of the main motivations for adopting this Bayesian methodology is the enablement of posterior predictive checking, as described in Section 4.6. Finally, Bayesian methodology also makes it straightforward to incorporate any prior information we may have into the procedure, and to easily test different models of the observed power spectrum.

The Vaughan (2010) approach is a multi-stage process. In the first stage (see Section 4.1), Bayes' theorem is used to fit parameterized models S (e.g., a power law) to the Fourier spectral power density. This is achieved via a Markov-Chain Monte-Carlo (MCMC) procedure, which explores the parameter space of the chosen model to find the posterior probability distribution. From this the model that best fits the data may be extracted.

The next stage is to perform a model comparison, to establish which of the tested models S is most favored. Numerous methods are available to perform this model selection procedure. Having selected the most favored model, appropriate test statistics are calculated (Section 4.5). Then, (see Section 4.6) posterior predictive checking is performed to establish the distribution of the relevant test statistics. This allows us to establish whether the measured test statistics were extreme, and therefore whether the chosen model captures the statistics of the data. Hence, by choosing relevant models, this procedure can test whether a model that includes an oscillatory component is needed in order to explain the observed power spectra of solar and stellar flares.

In the following sections, we describe in brief the methodology used in this paper, which is similar to that used in Vaughan (2010). We begin with Bayes' theorem.

4.1. Determining the Model Parameters

Given that p(a|b) denotes the probability of a given b then Bayes' theorem states that

$$\begin{equation} p(H|D,I) = \frac{p(H|I) p(D|H,I)}{p(D|I)}, \end{equation} \tag{ 2 }$$

where H represents the hypothesis to be tested (in our context a model of the PSD), D is the data (in our case the measured PSD), and I represents the prior information available. The quantity p(D|H, I) is called the likelihood, and is the probability of the data D given the hypothesis H. The quantity p(H|I) is called the prior information, and relates what we know about the hypothesis from our previous information before we take any data. The denominator p(D|I) represents the prior probability of the data, and is a constant. The left-hand side—p(H|D, I)—is the posterior probability of the hypothesis after we have acquired some data, and expresses all of what we know about the hypothesis H (see Gregory 2010).

Vaughan (2010) uses Bayes' theorem to compare hypotheses—or models—of red-noise spectra to the observed PSD of Seyfert 1 galaxies. The same treatment may be applied to solar flares. A model spectral power density S_j = S(f_j, Θ) is chosen, where f_j is the Fourier frequency, and Θ = (θ₁, ....θ_m) represents the m parameters of the model (for example, the parameters of a simple power-law model consists of a power-law index and a normalization constant). The periodogram of any observed stochastic time series of length N has Fourier power D^obs = ($D^{{\rm obs}}_{1},...,D^{{\rm obs}}_{N/2}$) at Fourier frequency f_j = j/NΔT (with j = 1, ..., N/2) and is exponentially distributed (Press et al. 1992; Chatfield 2003) about the true spectral density S_j = S(f_j). This allows us to write the likelihood as

$$\begin{equation} p(D^{{\rm obs}}|\Theta,I) = \prod ^{N/2}_{j=1} \frac{1}{S_j} \exp \left(- \frac{D_j^{{\rm obs}}}{S_j} \right), \end{equation} \tag{ 3 }$$

where we have replaced the hypothesis H with its equivalent representation Θ.

The posterior p(H|D, I) can be calculated as a function of its parameters Θ using an MCMC method. In this work, the posterior distributions are calculated using PyMC, a Python analysis package that implements such an MCMC procedure. MCMC methods allow for the efficient mapping of Bayesian posterior probability density functions in multi-dimensional parameter space. After some initial period (known as "burn-in"), the Markov chain returns samples of the parameters Θ directly proportional to their probability density as defined by the Bayesian posterior; that is, the equilibrium distribution of the Markov chain is the same as the posterior probability density function (Gregory 2010). Therefore, volumes of the parameter search space containing probable solutions are sampled more often in preference to volumes containing less probable solutions. This Bayesian/MCMC approach calculates the probability density function (PDF) for Θ. This posterior PDF contains a lot of information: for example, we can find the best fit by finding Θ_mode, the value of Θ that has the maximum posterior probability. The posterior PDF can also be used to find the best-fit and probability distribution of each fit parameter (e.g., Ireland et al. 2013).

4.2. Candidate Models

In this work, we focus on two candidate models that may describe the PSD of flare time series data. First, we consider a single power-law model plus a constant, i.e.,

$$\begin{equation} S_A = A f^{-\alpha } + C. \end{equation} \tag{ 4 }$$

This model is physically motivated by the understanding that power-law Fourier power spectra are a common feature of astrophysical objects, as discussed in Section 1. The addition of the constant C accounts for the transition between a power-law ("red") regime to a white-noise regime or a Poisson regime, as observed by Cenko et al. (2010) among others.

The second model includes an additional component that accounts for the appearance of an oscillatory component in the data, in addition to any observed power laws. This model may be written,

$$\begin{equation} S_B = S_A + B \exp \left(\frac{-(\ln f - \beta)^2}{2\sigma ^2} \right), \end{equation} \tag{ 5 }$$

where $B = B_0 / \sqrt{2\pi \sigma ^2}$, such that the integral over the Gaussian component is equal to 1. Hence model S_B is model S_A with a function added that is equivalent to a Gaussian in log-frequency space. The width of this Gaussian is then given by σ, its location in log-frequency by β, and its amplitude by B. This function represents excess power in the signal that may arise due to an oscillation.

4.3. Selection of Priors

In Bayesian analysis, the choice of prior probabilities affects the final posterior probability distribution (Equation (2)). These priors represent the information we already possess about the model characteristics. In this work, we adopt uniform prior probabilities for the model parameters:

$$\begin{eqnarray} &&-10 < \log A < 10, \nonumber\\ &&-6 < \log B_0 < 5, \nonumber \\ &&-20 < \log C < 10, \\ &&-1 < \alpha < 6, \nonumber\\ &&-6.0 < \beta < {-}2.0, \nonumber\\ &&0.05 < \sigma < 0.25, \nonumber \end{eqnarray} \tag{ 6 }$$

Hence we define the prior probability of each parameter to be uniform within a set range, i.e., any value within that range is equally likely. This intentionally simple choice reflects our lack of precise prior knowledge about the signal and ensures that we do not accidentally exclude the optimal values of the model parameters from the analysis.

4.4. Model Selection

In order to determine which of the candidate models is a preferred fit to the data, we use the Bayesian Information Criterion (BIC; Schwarz 1978; Burnham & Anderson 2004). Similarly to log-likelihood measurements, the BIC is minimized in the process of finding the best-fit of the model parameters to the data. The BIC criterion is defined as

$$\begin{equation} {\rm BIC} = -2 \ln (L) + k \ln (n) \end{equation} \tag{ 7 }$$

for large n, where n is the number of data points, k is the number of parameters in the model, and L is the maximum likelihood.

A feature of the BIC is that, as Equation (7) shows, it includes built-in consideration of the number of parameters k in the model. Models with more parameters are penalized in comparison to those with fewer parameters, to compensate for that fact that the more complex model should always fit at least as well as the simple one.

In this paper, we will use the BIC to perform the model comparison. By comparing the BIC of two models one may find the extent to which one model is preferred over the other, e.g.,

$$\begin{equation} \Delta {\rm BIC} = {\rm BIC}_{A} - {\rm BIC}_{B}. \end{equation} \tag{ 8 }$$

Since lower values of BIC are preferred, a negative value of ΔBIC indicates that model S_A is preferred over model S_B, whereas a positive value indicates a preference for S_B. In general, a value of ±10 in ΔBIC is considered highly significant (see, e.g., Burnham & Anderson 2004).

It is important to note, however, that the BIC merely measures which of the tested models is more appropriate, or alternatively which model minimizes information loss. However, it does not explicitly test whether either model is a good choice in absolute terms—both models may be poor fits to the observed data. In order to determine this, additional test statistics are required, as described in Section 4.5 below.

4.5. Test Statistics

In addition to model selection via the BIC, we include other measures of the appropriateness of the chosen models. We consider two test statistics, defined as follows (Vaughan 2010).

$$\begin{equation} T_{{\rm SSE}} = \sum _j \left(\frac{D_j - S_j}{S_j} \right)^2 \end{equation} \tag{ 9 }$$

and,

$$\begin{equation} T_R = \max _j \left(2 D_j / S_j \right). \end{equation} \tag{ 10 }$$

The statistic T_SSE corresponds to a global goodness-of-fit measure of the model S to the observations D, similar to χ². T_R measures the maximum deviation of the observations D from the model S. Hence T_R is a valuable statistic for finding local anomalies between I and S, such as might occur due to a narrowband oscillation.

The distribution of test statistic values, and hence to derive p-values corresponding to the values of T_R and T_SSE, is achieved via posterior predictive checking, described in Section 4.6. This allows a determination of whether the measured values of T_R and T_SSE are extreme compared to this distribution, which could indicate a poor choice of model.

4.6. Posterior Predictive Checking: Estimating Test Statistic Distributions

4.7. Method Summary

5.1. The Flare of 2011 February 15

Here we apply the method described above to the GOES-class X2.2 flare of 2011 February 15, which originated from NOAA active region (AR) 11158. This flare was previously analyzed by Dolla et al. (2012), who concentrated on the correlation between short-term variations in the signal at various wavelengths using data from a variety of instruments. An explicit significance test of the observed variations in the signal was not performed. This flare has also been studied by various other authors, including Schrijver et al. (2011), Tan et al. (2012), Kerr & Fletcher (2014), and Wang et al. (2014).

For this event, the most sunward facing detectors of GBM show evidence of pulse pile-up effects due to high count rates. Therefore, we utilize detector NaI-2 for this analysis, which was pointed at a greater angle from the Sun, resulting in manageable counts. We apply the analysis method to the Fourier power spectra obtained from three X-ray energy ranges observed by GBM at 1 s cadence (12–27 keV, 27–50 keV, 50–100 keV), and to the Al and Zr filter data from LYRA, also at 1 s cadence. Figure 2 shows an example of the method applied to the 12–27 keV data. The left panel shows the original signal, while the second panel shows the Fourier power spectrum. In this example, the simple model S_A was strongly preferred over S_B, with ΔBIC = −14.4. The remaining two panels illustrate the T_R and T_SSE distributions obtained by posterior predictive checking as described in Section 4.6, where the actual observed values of these statistics are denoted by the red lines. From this it can be seen that the observed values of T_R and T_SSE are not extreme within the context of the power-law model.

Zoom In Zoom Out Reset image size

As Table 1 shows, the model comparison results for each X-ray band observed by GBM are similar, resulting in a negative value of ΔBIC, with values ranging from −14.4 at 12–27 keV to −5.4 in the 50–100 keV range. Similarly, for both LYRA filters the value of ΔBIC is < − 10. Hence the simple power-law model S_A is preferred in all cases over the more complex model S_B.

Table 1. Summary of Model Comparison Results and Statistics in Selected Energy Ranges for Four Flare Events

Event	Fig. Ref.	Obs. Instrument	Obs. Band	Obs. Interval	ΔBIC	Preferred Model	−α	p-value	p-value
				(UT)				(TR)	(TSSE)
2011 Feb 15	3(b)	GBM	12–27 keV	01:44–02:22	−14.4	SA	−2.89 ± 0.07	0.714	0.242
	3(c)	GBM	27–50 keV	01:44–02:22	−10.7	SA	−2.30 ± 0.06	0.261	0.077
	3(d)	GBM	50–100 keV	01:44–02:22	−5.4	SA	−2.03 ± 0.07	0.124	0.082
	3(e)	LYRA	Al (17–80 nm)	01:44–01:55	−14.3	SA	−3.52 ± 0.12	0.177	0.021
	3(f)	LYRA	Al (17–80 nm)	01:58–02:06	−14.2	SA	−4.90 ± 0.11	0.688	0.663
	3(g)	LYRA	Zr (6–20 nm)	01:44–01:55	−13.8	SA	−3.38 ± 0.13	0.649	0.368
	3(h)	LYRA	Zr (6–20 nm)	01:58–02:06	−18.8	SA	−4.95 ± 0.13	0.345	0.186
2011 Jun 7	4(b)	GBM	12–27 keV	06:16–06:58	−18.2	SA	−2.90 ± 0.07	0.990	0.688
	4(c)	GBM	27–50 keV	06:16–06:58	−9.6	SA	−3.12 ± 0.07	0.928	0.932
	4(d)	GBM	50–100 keV	06:16–06:58	−4.2	SA	−2.81 ± 0.08	0.844	0.871
	4(e)	LYRA	Al (17–80 nm)	06:16–06:27	−19.7	SA	−3.48 ± 0.1	0.035	0.046
	4(f)	LYRA	Al (17–80 nm)	06:30–06:50	−3.7	SA	−2.96 ± 0.11	0.110	0.742
	4(g)	LYRA	Zr (6–20 nm)	06:16–06:27	−9.7	SA	−3.46 ± 0.09	0.565	0.068
	4(h)	LYRA	Zr (6–20 nm)	06:30–06:50	−2.3	SA	− 2.80 ± 0.10	0.007	0.014
1998 May 8	5	NoRH	17 GHz	01:54–02:02	80.4	SB	−3.06 ± 0.12	0.824	0.798
	5	Yohkoh/HXT	13–23 keV	01:56–02:00	16.1	SB	−3.70 ± 0.66	0.146	0.222
	5	Yohkoh/HXT	23–33 keV	01:56–02:00	−3.4	SA	−2.27 ± 0.17	0.441	0.171
	5	Yohkoh/HXT	33–53 keV	01:56–02:00	−19.0	SA	−1.86 ± 0.13	0.641	0.929
2009 Jan 16	6	NMSU	U-band	04:03–07:35	−15.6	SA	−3.55 ± 0.12	0.999	0.998
"megaflare"

Download table as:  ASCIITypeset image

Figure 3 shows the lightcurves and best fits of model S_A to the Fourier power spectra each of the investigated energy bands. For brevity the T_R and T_SSE distributions are not shown, but for reference each best-fit power law is accompanied by a 99% significance line, which is obtained by finding the value of T_R at each f where the p-value p_TR < 0.01. It can therefore be seen that the variations in the Fourier power spectra are not extreme in the context of the power-law model expectations. The p-values for T_R and T_SSE are listed in Table 1.

Zoom In Zoom Out Reset image size

5.2. The Flare of 2011 June 7

This GOES-class M2.5 flare originated from AR 11226 and has been studied by various authors, including Innes et al. (2012), Inglis & Gilbert (2013), Reale et al. (2013), Gilbert et al. (2013), and Carlyle et al. (2014), primarily due to the spectacular prominence eruption associated with this event. However, this flare also exhibited strong pulsations in UV and X-ray wavelengths, which were investigated in the Inglis & Gilbert (2013) study. As with the event studied in Section 5.1, a statistical significance test was not performed on these pulsations, and their exact nature was left for debate in favor of timing studies and correlations with other flare parameters.

Figure 4 shows the result of considering power-law Fourier power spectra by applying the analysis method described in Section 4. For each of the six wavelengths studied, model S_B, which includes an additional spike that would be consistent with a narrowband oscillation, is never preferred over the pure power-law model S_A (see Table 1 for the associated ΔBIC values), although in the late phase of the Al and Zr channel data there is little difference between the two models. Additionally, the test statistics T_R and T_SSE describing the goodness of fit for S_A, listed in Table 1, are generally not extreme, indicating an overall good fit of the model S_A to the data D in most cases. An exception to this is the 12–27 keV data, where p_TR is 0.01, suggesting that although S_A is preferred over S_B, it is not a complete description of the observed data. Overall however, despite the visually striking pulses in the flare lightcurves, the observed PSDs at both X-ray and EUV wavelengths are consistent with and well-described by a power-law model in the Fourier power spectrum, without the need for explicit oscillations in the model.

Zoom In Zoom Out Reset image size

5.3. The Flare of 1998 May 8

The GOES-class M3.1 flare of 1998 May 8 originated from AR 8210 and began at approximately 01:49 UT in soft X-rays. It was also observed in radio and hard X-rays by the NoRH and the Yohkoh satellite, respectively, and is a pronounced example of a QPP event. Analysis of this data in two previous studies (Inglis et al. 2008; Stepanov et al. 2004) suggested the presence of a statistically significant oscillation with P ≈ 16 s, which was interpreted by Inglis et al. (2008) as the manifestation of a magnetoacoustic sausage mode. However, these studies did not account for the power-law-like nature of the signal Fourier power spectrum, and hence may have overestimated the statistical significance of their results.

Revisiting this event, the method described in Section 4.7 is applied to both 1 s cadence NoRH data at 17 GHz, and to 0.5 s cadence Yohkoh X-ray data. The results of the model comparison are shown in Figure 5. We find that, for the 17 GHz emission, ΔBIC ≈ 80, strongly favoring model S_B which includes a Gaussian bump in addition to a power law. The best-fitlocation for this bump is at f₀ = 0.062 ± 0.003 Hz, equivalent to P = 16.2 ± 0.8 s. Although the position of the Gaussian bump is determined very accurately, in this context the width of the Gaussian is also of interest, as it encapsulates information about the range, or variation, of the observed "quasi"-period. Hence, when the Gaussian width is included in the period uncertainty we obtain P = $16.2^{+3.8}_{-3.1}$ s. This is consistent with the results of Inglis et al. (2008) and Stepanov et al. (2004), who both estimated that P ≈ 16 s.

Zoom In Zoom Out Reset image size

In X-rays, analysis of the Yohkoh data in the 13–23 keV channel yields ΔBIC ≈ 16, also indicating a strong preference for model S_B. The best-fit location of the frequency bump with a 1σ uncertainty is f₀ = 0.071 ± 0.004 Hz, corresponding to P = 14.0 ± 0.75 s. As before, when we utilize the Gaussian width parameter to estimate the bounds of the quasi-period, we obtain P = $14.0^{+3.4}_{-2.8}$ s, consistent with the observations at radio wavelengths. In the higher energy X-ray channels, the preference for model S_B is not replicated—for the 23–33 keV channel, model S_A is marginally preferred with ΔBIC ≈ −3.4, while for the 33–53 keV channel model S_A is strongly preferred with ΔBIC ≈ −19. As Figure 5 shows however, the fitting of model S_B to the 23–33 keV channel, where the simple power law is marginally preferred, yields the same best-fit period of P ≈ 14 s as the 13–23 keV channel. This suggests that an oscillatory signature may indeed be present in the 23–33 keV data, but at insufficient strength to be unambiguously detected.

Examination of the test statistic values for this event (see Table 1) reveals that the measured T_R and T_SSE values for the 17 GHz data are smaller than the mean of the distribution obtained via posterior predictive checking (see Section 4.6), leading to p-values of 0.842 and 0.798, respectively. For the 13–23 keV data we find the reverse effect, but in neither case are the test statistic values extreme. However, although a good fit was obtained here, it should be noted that our choice of a Gaussian function is a purely empirical one; as pointed out by Vaughan (2010), the true statistics of the Fourier power spectrum in the region of the narrowband peak are unclear. Hence, although model S_B is strongly preferred for both the 17 GHz and 13–23 keV data, indicating the presence of an oscillation with P ≈ 14–16 s, the Gaussian function is likely not a statistically complete choice for modeling this oscillation.

The critical finding is that, for this event, a simple power-law Fourier power spectrum is not sufficient to explain the flare data at radio or soft X-ray wavelengths. Instead, there is strong evidence at both radio and X-ray energies for the presence of a narrow-band oscillation. This is an extremely important result, which confirms the existence of oscillation-like signatures in solar flare signals across multiple energies and emission regimes. This observation is consistent with classical ideas of QPPs, such as bursty reconnection processes or MHD wave modes.

5.4. The Stellar Megaflare of 2009 January 16

As a final example, we consider the stellar megaflare of 2009 January 16, previously studied by Kowalski et al. (2010) and Anfinogentov et al. (2013). This flare originated from the dM4.5e star YZ CMi, and exhibited several emission peaks in addition to an exponential-like decay. These variations were interpreted by Anfinogentov et al. (2013) as a decaying long-period oscillation with P ≈ 32 minutes, suggesting a longitudinal magnetoacoustic mode as a possible mechanism. Here, we revisit the sub-minute cadence near-UV data from the NMSU 1m telescope studied by Kowalski et al. (2010) and Anfinogentov et al. (2013), and re-evaluate this flare signature accounting for power-law effects in the Fourier power spectrum.

Figure 6 illustrates the U-band lightcurve and the best-fit Fourier power spectrum for this event. Here, model S_A is preferred over S_B by a substantial margin, with ΔBIC ≈ −16. Hence, despite the intriguing variations in the observed time profile, the emission from this stellar flare is better explained in terms of the power-law model than the model containing a narrowband oscillation. The observed T_R and T_SSE values are rather low, however, compared to their obtained posterior predictive distributions. This indicates that another untested model may better reproduce the statistics of the data than either S_A or S_B.

Zoom In Zoom Out Reset image size

Having studied 19 time series from 4 flares, we have found strong evidence for the existence of an explicit oscillation only in the 17 GHz radio data and the 13–23 keV X-ray data from the 1998 May 8 flare. For the remaining events, the simple power-law model S_A was preferred as a description of the data over the model S_B including an extra peak. However, in some cases the models were considered almost equally likely. For example, the Al and Zr channel during the latter part of the 2011 June 7 flare (Figure 4) showed ΔBIC ≈ −3.7 and ΔBIC ≈ −2.3, respectively, only marginally favoring S_A over S_B. Hence, for these signals, we do not rule out the presence of an additional component—we can only say that we do not have strong evidence in favor of one.

It is also important to emphasize that the adequate description of flares containing substantial temporal variations in terms of a power-law-like Fourier power spectrum does not mean that those variations are not real or do not exist. These variations are observed by multiple instruments and are of solar origin. Hence, studies of the timing of these variations, their variations across wavelengths and energy regimes (e.g., Dolla et al. 2012), and their correlation with other measurable flare parameters (Inglis & Gilbert 2013) remain crucial toward our understanding of flares. The results of Section 5 indicate that, when the frequency-dependent properties of flares are properly considered, an explicit oscillatory signal is not required in order to explain these observations. Instead, we must consider these variations as an intrinsic property of the flare system.

One possible mechanism that would give rise to a power-law structure in solar flares was suggested by Aschwanden (2011). In this model, a flare is considered as a superposition of many exponentially decaying energy deposition events. If the number of events is $N(E) \propto E^{-\alpha _{E}}$, and the energy in each event is E(T)∝T^{1 + γ}, where T is exponential decay time, then the resulting Fourier spectrum should have a power-law index

$$\begin{equation} p = -(2 - \alpha _E)(1+\gamma). \end{equation} \tag{ 11 }$$

Hence the observed power-law index can in principle provide information about the energy release process in a flare.

McAteer et al. (2007) also investigated the nature of bursty solar X-ray emission in the context of self-organized criticality and fractal behavior. In that work, the Hölder exponent was determined—using RHESSI data—as a function of X-ray energy for the X4.8-class flare of 2002 July 23. The Hölder exponent is related to the Fourier power spectral index and is a measure of the persistence—i.e., how later observed values in a signal depend on earlier values—of the physical system. In a persistent walk regime, a generally increasing signal means that the next observed value in time is also likely to show an increase ("smooth behavior"), whereas in an anti-persistent regime the opposite occurs, and a decrease is likely ("bursty behavior"). In general, it was found that, for the 2002 July 23 flare, the persistence of the signal decreased as a function of energy; low energy bands (e.g., 6–12 keV) were consistent with a persistent-walk regime, while at higher energies (e.g., 50–100 keV) anti-persistent behavior became more pronounced. This was also reflected in the measured Fourier power-law exponents; the power-law index α ranged from ≈−2.6 at 3–6 keV down to ≈−1.8 at 300–800 keV, indicating different signal characteristics.

Figure 7 shows the distribution of measured α values for the 2011 February 15, 2011 June 7, and 1998 May 8 events as a function of energy. The 2009 January 16 flare measurement is also shown. For 2011 February 15, the results of McAteer et al. (2007) are reproduced; the measured spectral index in hard X-rays monotonically decreases as a function of energy from ≈−2.9 in the 12–27 keV range, down to ≈−2.0 at 50–100 keV. Similarly, the three datapoints available from Yohkoh/HXT data on 1998 May 8 show a monotonic decrease versus energy, although it should be noted that the spectral index is very poorly constrained in the 13–23 keV channel. The 2011 June 7 flare however displays inconsistent results: the measured exponents remain approximately constant at ≈−3.0 as a function of energy in X-rays (see Figure 7). This suggests that flares may not universally follow the monotonic trend in spectral index and fractal dimension observed by McAteer et al. (2007); different flares may exhibit different emission characteristics as a function of energy. This is not unexpected given the earlier results of Aschwanden et al. (1998) who, in a large study of events in >25 keV X-rays observed with CGRO/BATSE, found substantial variation of the power-law index between flares (−1.5 <α < −3.2 for strong events).

Zoom In Zoom Out Reset image size

In interpreting these results, we should be wary of potential biases that may affect the measurement of α. One potential issue is a dependence of α on the total fluence of the data, however, this possibility is mitigated by the treatment of the input signal described by Equation (1). Additionally, the choice of GBM detector for the 2011 February 15 and 2011 June 7 events ensures that pile-up effects should not be an issue in the X-ray data. The other main factor in determining α, the length of the power-law component in the Fourier power spectrum, is captured by the uncertainties listed in Table 1. In general, we find that α is well and consistently constrained, with the notable exception of the 13–23 keV X-ray data from Yohkoh in the 1998 May 8 flare. In this case, the power-law slope is very poorly constrained due to the relatively early transition to a white-noise regime (see Figure 5) and the location of the Gaussian bump.

An additional feature of these results is that not only do the observed Fourier power-law indices vary between flares (e.g., Aschwanden et al. 1998); the observed power-law index of both LYRA channels changes substantially between the two sub-intervals of both the 2011 June 7 and 2011 February 15 flares, as highlighted in Figure 7. This could be related to the different physical phases of solar flares, e.g., the early "rise phase" that includes impulsive X-ray emission, and the long decay phase often observed at EUV wavelengths. For the 2011 February 15 event, the power-law index is much steeper in both the aluminum and zirconium filter channels during the observed decay phase of the LYRA emission compared with the earlier rise phase. In the context of Equation (11), this may indicate a change in either the number distribution of heating events, or a change in the energy distribution of these events, which would be expected during a transition to the flare decay phase. For the 2011 June 7 event, the difference in observed Fourier power-law index is much less pronounced (Figure 7) in the LYRA channels. One difference is that in this case, the later sub-intervals shown in Figures 4(e) and (g) are associated with a peak of observed EUV emission, rather than a strong decay phase, which is not observed due to an LAR of the LYRA spacecraft (see Section 3.1).

It is clear that power laws in the Fourier power spectrum are an intrinsic property of both solar and stellar flare time series, and one that must be taken into account when searching for pulsations or oscillations in flares. As has been shown here, events that may have been previously considered as containing signatures of oscillations can in fact be adequately described via a power-law model of the Fourier power spectrum. Hence, the prevalence of explicit oscillations in flares is likely to have been substantially overestimated in past literature, although few large-sample studies exist. However, the analysis of the 1998 May 8 flare shows that, for at least some events, there is strong evidence of oscillatory signals in the data over and above the background frequency-dependent spectrum, in both radio and X-ray emission regimes. In these cases, classical interpretations of QPPs, such as regimes of bursty reconnection or signatures of MHD wave modes, remain valid. Additionally, as discussed in Section 6, the detection and measurement of frequency-dependent Fourier spectra may provide us with a diagnostic tool for understanding the fundamental processes of flare energy release. The measured power-law index as a function of energy provides us with information regarding the energy release process in different energy regimes.

These results are relevant not just to flare studies, but throughout solar physics, and coronal seismology in particular, where many searches for waves and oscillations implicitly assume a Gaussian noise regime. For example, Auchère et al. (2014) recently observed that, for a number of regions of interest, the Fourier power spectra of long-duration EUV data obtained by SOHO/EIT also obey a power law. Similarly, Ireland et al. (2014) reveal Fourier power-law properties in the time series of ARs obtained in EUV by SDO/AIA. Hence, power-law Fourier power spectra must be accounted for in many solar features over many timescales and in multiple emission regimes.

A future large-scale study of this kind is required in order to understand how the Fourier power-law index may vary as a function of energy between flares, during the flares themselves, and how this information may be used to diagnose conditions at the energy release sites in flares. Such a study will also establish the true prevalence of explicit oscillations in solar and stellar flares.

A.R.I. and J.I. both gratefully acknowledge the support of the PROBA2 Guest Investigator Programme, which provided the opportunity and funds for both authors to collaborate with the PROBA2 team at the Royal Observatory of Belgium. The authors are also grateful to Dr. Adam Kowalski for providing the U-band data from the NMSU 1 m telescope. This work was carried out using the SunPy, SciPy, PyMC, and Matplotlib software packages. We also acknowledge the anonymous referee for useful comments that resulted in a number of improvements to this paper.