EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 602 (June 2017) A&A, 602 (2017) A47 Full HTML
Free Access
Issue
A&A
Volume 602, June 2017
Article Number A47
Number of page(s) 8
Section The Sun
DOI https://doi.org/10.1051/0004-6361/201730595
Published online 07 June 2017

© ESO, 2017

1. Introduction

Quasi-periodic pulsations (QPPs) have been widely observed in solar flares (e.g. Simões et al. 2015; Hayes et al. 2016; Myagkova et al. 2016) after they were first discovered by Parks & Winckler (1969), and they are also occasionally observed in stellar flares (e.g. Mathioudakis et al. 2003; Pugh et al. 2016). Since QPPs are a common feature of flares, the nature of them should be understood in order to fully understand flares. It is thought that QPPs could be the result of magnetohydrodynamic oscillations or a regime of repetitive magnetic reconnection: also referred to as “load/unload” mechanisms or “magnetic dripping”, where free magnetic energy continuously builds up but is released repetitively each time some threshold energy is surpassed (Nakariakov & Melnikov 2009; Nakariakov et al. 2016; Van Doorsselaere et al. 2016). The “magnetic dripping” term arises because an analogy can be made between this mechanism and leaking water accumulating at a steady rate at the bottom of a surface, which drips each time the weight of the water is great enough to overcome the surface tension. Although there is no strict definition of QPPs, it is generally accepted for stationary QPPs that the impulsive and/or decay phase of the flare should contain, at the very least, three cycles of oscillation, or pulses with approximately equal time spacing, visible above the noise level. There may also be non-stationary QPPs, where the time spacing between pulses increases or decreases in a non-random fashion (e.g. Kupriyanova et al. 2010; Huang et al. 2014).

The most commonly used methods to test for the presence of QPPs in a flare involve looking for a significant peak in the periodogram or wavelet spectrum. Early examples of these methods are shown by Lipa (1978), Aschwanden et al. (1998), and for more recent examples see Reznikova & Shibasaki (2011), Chowdhury et al. (2015), Kumar et al. (2016), Tian et al. (2016). The periodogram or wavelet spectrum of the unmodified flare light curve will, however, have power that is dependent on the frequency, following a power law relationship (P ∝ f− α). This spectral behaviour is the result of red noise in the light curve, which is intrinsic to flare time series data (McAteer et al. 2007; Gruber et al. 2011). In order to remove this red noise behaviour and make peaks due to a periodic component of the signal more prominent in the power spectrum, the flare light curve is often detrended. This is equivalent to the separation of the different physical phenomena of different time scales, and hence seems to be well justified. By removing the longer time-scale flare profile from the time series data, only the shorter time-scale QPPs (if present) and noise should be left. This is usually done by subtracting either a model fitted to the flare profile (e.g. Anfinogentov et al. 2013), a boxcar smoothed version of the time series with a pre-selected boxcar width (e.g. Kupriyanova et al. 2013), or an aperiodic trend determined by empirical mode decomposition (e.g. Cho et al. 2016). It has been noted, however, that detrending can lead to overestimating the significance of peaks in the power spectrum (Inglis et al. 2015) due to the artificial suppression of other spectral components (Gruber et al. 2011). A more serious potential consequence of detrending with an inappropriate model or a boxcar smooth is the introduction of a signal that looks periodic, but does not exist in the original data (Auchère et al. 2016).

In terms of assessing the significance of a peak in a power spectrum, Scargle (1982) showed how the false alarm probability can be found for data with white noise. For evenly time-spaced data, the periodogram is equivalent to the Fourier power spectrum (with additional normalisation), which is equal to the sum of the squares of the real and imaginary parts of the Fourier transform. For a white noise time series, where each value is drawn at random from a Gaussian distribution, the real and imaginary parts of the Fourier transform should also be Gaussian distributed random variables. Squaring a Gaussian distributed random variable results in a chi-squared, one degree of freedom (d.o.f.) distributed random variable, and adding together two chi-squared, one d.o.f. distributed variables results in a chi-squared, two d.o.f. distributed random variable. The probability density of a chi-squared, two d.o.f. distributed variable, x, which has a mean value of two is: pχ22(x)=12e−x/2.\begin{equation} p_{\chi ^2_2}(x) = \frac{1}{2}\mathrm{e}^{-x/2}. \end{equation}(1)The probability of having a value X that is greater than some threshold x′ is therefore: Pr{X>x′}=∫x′∞12e−x/2dx=e−x′/2.\begin{equation} \textit{Pr}\left\{X>x'\right\} = \int_{x'}^{\infty} \frac{1}{2}\mathrm{e}^{-x/2} \mathrm{d}x = {\rm e}^{-x'/2}. \end{equation}(2)Since the power spectrum is positive everywhere, x and x′ in the above equation are always positive. When considering a power spectrum sampled at N′ = N−1 independent frequencies (the Nyquist frequency is neglected as it follows a chi-squared one d.o.f. distribution; Vaughan 2005) the probability is equivalent to: Pr{X>x′}=1−(1−ϵN′)1/N′≈ϵN′/N′,\begin{equation} \textit{Pr}\left\{X>x'\right\} = 1 - \left(1 - \epsilon _{N'}\right)^{1/N'} \approx \epsilon _{N'}/N', \end{equation}(3)where ϵN′ is the false alarm probability, and the approximation holds when is ϵN′ small (Chaplin et al. 2002). The false alarm probability is defined as the probability of observing a peak in the power spectrum above some threshold power (Scargle 1982; Horne & Baliunas 1986), and this threshold power can be calculated by equating Eqs. (2) and (3) and solving for x′. For example, the 99% confidence level in the power spectrum is the power threshold corresponding to a false alarm probability of 0.01, and refers to the level above which there is only a 1% chance of observing a peak in the power spectrum of a Gaussian distributed random (white noise) time series. The above expressions are only valid, however, when considering data with white noise that is chi-squared, two d.o.f. distributed in the power spectrum. When considering solar and stellar flare time series data, the removal of an imperfect approximation of the underlying flare trend could alter the underlying noise distribution, leading to the calculation of a misleading confidence level. Therefore extreme care should be taken when assessing the significance of peaks in the power spectrum of detrended data.

Another point to consider when detrending, is that each flare must be treated separately in the analysis. For example, a model that gives a good fit to the flare profile must be chosen, along with suitable initial estimate parameters, but some flares are very complex in shape (e.g. Davenport et al. 2014), which makes finding a general model that fits all flares to a satisfactory standard more difficult. In addition to this, the same flare can look very different when viewed in different energy bands. On the other hand, if a boxcar smooth is used to detrend the most suitable boxcar width must be chosen for each flare. When undertaking a large-scale study of a number of events a minimal amount of manual intervention is preferable, hence methods that avoid detrending are more appropriate, such as that used by Inglis et al. (2016).

Vaughan (2005) demonstrates a method to assess the significance of a peak in a power spectrum, and shows how it can be used to test for periodic signals in X-ray light curves of active galaxies. The method avoids detrending and takes full account of red noise and data uncertainties. We build on this method and show in detail how it can be applied to solar flare data, so that peaks in the power spectrum found above a certain confidence level may be considered as candidate QPPs. For this study we address only stationary QPPs, with constant periods and without phase modulation. In Sect. 2 we summarise the method of Vaughan (2005), derive a more simple form of the equation to be solved to determine the confidence level, and describe how this method can be used to search for QPPs in flare light curve data. In Sect. 3 we adapt the method to be used with rebinned power spectra, which can help detect QPPs with a broad peak in the power spectrum. The results of testing the methods on simulated data are shown in Sect. 4, and in Sect. 5 a few examples are given where the methods have been applied to real solar flare data. A summary is given in Sect. 6.

2. Confidence levels on power-law power spectra

From Vaughan (2005), if I(fj) is the periodogram power at a particular frequency, fj, we have: I(fj)=𝒫(fj)χ22/2,\begin{equation} I\left(f_j\right) = \mathcal{P}\left(f_j\right)\chi ^2_2/2, \end{equation}(4)where \hbox{$\mathcal{P}(f_j)$} is the “true” power, and χ22/2\hbox{$\chi ^2_2/2$} is the chi-squared two d.o.f. distributed noise. For data with red noise, the power spectrum will follow a power law, which can be fitted with a straight line when working in log space. Solar flare power spectra often follow a broken power law, since the red noise component can fall below the white noise level at high frequencies (Gruber et al. 2011; McAteer et al. 2016). The broken power law model can be written as: log[𝒫̂(f)]={iff<fbreakiff>fbreak,\begin{equation} \log\left[\hat{\mathcal{P}}(f)\right] = \begin{cases} -\alpha\log\left[f\right] + c & \text{if } f < f_{\rm break} \\ -\left(\alpha \!-\! \beta\right)\log\left[f_{\rm break}\right] \!-\! \beta\log\left[f\right] \!+\! c & \text{if } f > f_{\rm break}, \end{cases} \end{equation}(5)where fbreak is the frequency at which the power law break occurs, α and β are power law indices and c is a constant. The probability density for 2Ij (the factor of two appears because the chi-squared distribution is conventionally defined assuming the values following that distribution have a mean equal to the number of degrees of freedom) can then be written as: p2Ij(x)=12𝒫je−x/2𝒫j.\begin{equation} p_{2I_j}(x) = \frac{1}{2\mathcal{P}_j}\mathrm{e}^{-x/2\mathcal{P}_j}. \end{equation}(6)When fitting the “true” power spectrum, \hbox{$\mathcal{P}$}, with a broken power law model, \hbox{$\hat{\mathcal{P}}$}, the uncertainties on this fitted model will follow a Gaussian distribution in log space, and hence a log-normal distribution in linear space: p𝒫̂j(y)=12πySjexp⎧⎪⎪⎪⎨⎪⎪⎪⎩−(ln[y]−ln[𝒫j])22Sj2⎫⎪⎪⎪⎬⎪⎪⎪⎭,\begin{equation} p_{\hat{\mathcal{P}}_j}(y) = \frac{1}{\sqrt{2\pi\;}yS_j}\exp{\left\{-\frac{\left(\ln[y] - \ln[\mathcal{P}_j]\right)^2}{2S^2_j}\right\}}, \end{equation}(7)with Sj=err{log[𝒫̂(fj)]}×ln[10],\begin{equation} S_j = \text{err}\left\{\log\left[\hat{\mathcal{P}}\left(f_j\right)\right]\right\} \times \ln[10], \end{equation}(8)where \hbox{$\text{err}\left\{\log\left[\hat{\mathcal{P}}(f_j)\right]\right\}$} is the uncertainty on the fitted model in log space, and the ln [10] factor accounts for the fact that the uncertainty is defined in terms of log base ten, whereas the log-normal distribution above is defined (by convention) in terms of log base e. In order to find the uncertainties on the model fitted to flare power spectra, uncertainties on the flare light curve data were used in Monte Carlo simulations. While some instruments that observe the Sun include uncertainties for the data provided, many do not. Fortunately in most cases reasonable estimates of the uncertainties can be made. For example, some X-ray observations of the Sun follow Poisson counting statistics, so the uncertainty on each count rate value can be found by taking the square root of the value. For the radioheliograph data used in Sect. 5, an estimate of the uncertainty can be found by calculating the standard deviation of several hours of flat data, in which no flares occur and there are no other features. For the Monte Carlo simulations, random numbers with a mean of zero and standard deviations equal to the uncertainties on the light curve values are added to each of the light curve values. The periodogram is then found, converted to log space, and a broken power law model is fitted using a least-squares method. This is repeated many times (10 000 times for the examples shown in Sect. 5), and for each iteration the initial guess parameters used in the model fit are allowed to vary, in order to prevent a local rather than a global minimum being found by the least-squares fit (for the examples in Sect. 5 parameters of α = 2, β = 1, c = −3 and fbreak = 0.1 were allowed to randomly vary with a standard deviation equal to 10% of the parameter value for each iteration). The distributions of the parameters from the repeated power spectrum fits should be approximately Gaussian, and so for each frequency bin of the power spectrum, the uncertainty on the fitted model value will be Gaussian distributed. Hence the distribution of fitted powers at each frequency index can be fitted by a Gaussian model, and the standard deviation of the Gaussian model can be used as an estimate of the uncertainty of the broken power law model at that index.

thumbnail Fig. 1

Plots of the integrand in Eq. (12) as a function of w. The solid black, dotted red and dashed blue lines show the function when Sj is equal to 0.4, 0.2 and 0.02 respectively. The value of γϵj has arbitrarily been chosen to be equal to 20, which is a typical value.

thumbnail Fig. 2

Plots of Eqs. (12) and (13) as a function of γϵj are shown by the solid red and dashed blue lines respectively. The values of Sj and ϵN′/N′ have arbitrarily been chosen to be equal to 0.2 and 0.01/100 respectively. The solution we require is where Eq. (12) is equal to Eq. (13), which corresponds to γϵj = 21.467.

The probability density of the ratio \hbox{$\hat{\gamma}_j = 2I_j/\hat{\mathcal{P}}_j$} (the power spectrum with the red noise component removed) can be found from (Curtiss 1941): pγj(z)=∫0+∞|y|p2Ij(zy)p𝒫j(y)dy,\begin{equation} p_{\gamma_j}(z) = \int_{0}^{+\infty} |y|p_{2I_j}(zy)p_{\mathcal{P}_j}(y) \mathrm{d}y, \end{equation}(9)where y and z are dummy variables representing different power levels in the power spectrum. The lower limit of this integral is zero rather than negative infinity because the power spectrum is always positive. Integrating this probability density between γϵj and infinity gives the probability that a value γjˆ\hbox{$\hat{\gamma _j}$} is greater than γϵj: Pr{γjˆ>γϵj}=18πSj𝒫j  ×∫γϵj∞∫0∞exp⎧⎪⎪⎪⎨⎪⎪⎪⎩−(ln[y]−ln[𝒫j])22Sj2−yz2𝒫j⎫⎪⎪⎪⎬⎪⎪⎪⎭dydz.\begin{eqnarray} &&\textit{Pr}\left\{\hat{\gamma _j} > \gamma _{\epsilon _j}\right\} = \frac{1}{\sqrt{8\pi\;}S_j\mathcal{P}_j} \nonumber\\ &&\qquad\times \int_{\gamma _{\epsilon _j}}^{\infty}\int_{0}^{\infty}\exp \left\{-\frac{\left(\ln[y] - \ln[\mathcal{P}_j]\right)^2}{2S^2_j} - \frac{yz}{2\mathcal{P}_j}\right\} \mathrm{d}y\mathrm{d}z. \end{eqnarray}(10)Substituting in another dummy variable, \hbox{$w = y/\mathcal{P}_j$} (with \hbox{$\mathrm{d}y = \mathcal{P}_j\;\mathrm{d}w$}), simplifies this to: Pr{γjˆ>γϵj}=18πSj∫γϵj∞∫0∞exp⎧⎪⎪⎨⎪⎪⎩−(ln w)22Sj2−wz2⎫⎪⎪⎬⎪⎪⎭dwdz.\begin{equation} \textit{Pr}\left\{\hat{\gamma _j} > \gamma _{\epsilon _j}\right\} = \frac{1}{\sqrt{8\pi\;}S_j}\int_{\gamma _{\epsilon _j}}^{\infty}\int_{0}^{\infty}\exp \left\{-\frac{\left(\ln\,w\right)^2}{2S^2_j} - \frac{wz}{2}\right\} \mathrm{d}w\mathrm{d}z. \end{equation}(11)Since the integrand is well-behaved and contains no discontinuities, the order of integration can be swapped (see Fig. 1 for plots of the function for different values of Sj). Then the function can be integrated with respect to z, to get: Pr{γjˆ>γϵj}=∫0∞12πSjwexp⎧⎪⎪⎨⎪⎪⎩−(ln w)22Sj2−γϵjw2⎫⎪⎪⎬⎪⎪⎭dw ,\begin{eqnarray} \textit{Pr}\left\{\hat{\gamma _j} > \gamma _{\epsilon _j}\right\} = \int_{0}^{\infty}\frac{1}{\sqrt{2\pi\;}S_jw} \exp \left\{-\frac{\left(\ln\,w\right)^2}{2S^2_j} - \frac{\gamma _{\epsilon _j}w}{2}\right\} \mathrm{d}w\,, \end{eqnarray}(12)which can be equated to (see Eq. (3) in Sect. 1): Pr{γjˆ>γϵj}≈ϵN′N′,\begin{equation} \textit{Pr}\left\{\hat{\gamma _j} > \gamma _{\epsilon _j}\right\} \approx \frac{\epsilon _{N'}}{N'}, \end{equation}(13)and solved numerically in order to find a γϵj corresponding to each value of the fitted model power spectrum, \hbox{$\hat{\mathcal{P}}_j$}. Figure 2 shows plots of Eqs. (12) and (13) as a function of γϵj, and the solution is where the two lines cross. In practise it is helpful to subtract Eq. (13) from Eq. (12), then the solution will be where this function is equal to zero and a root finding algorithm can be used. The final step is to ensure correct normalisation, since the conventional form of the probability density function for a chi-squared, two d.o.f. distribution (shown in Eq. (1)) assumes the values following that distribution have a mean equal to two, which is not necessarily the case. An important point to account for is that the mean calculated in log space is not the same as the log of the mean calculated in linear space (i.e. ⟨ log Ij ⟩ ≠ log  ⟨ Ij ⟩). Hence the power spectrum as well as the fit must be converted into linear space so that the mean, or expectation value, of the flattened power spectrum (denoted \hbox{$\langle I_j - \hat{\mathcal{P}}_j\rangle$}) can be found. Therefore the confidence level in log space is then equal to \hbox{$\log[\hat{\mathcal{P}}_j] + \log[\gamma _{\epsilon _j}\langle I_j - \hat{\mathcal{P}}_j\rangle /2]$}.

3. Confidence levels on rebinned power spectra

Appourchaux (2004) shows how rebinning the power spectrum can improve the detection of short-lived solar acoustic modes, such as wave trains with a highly modulated amplitude, which have power spread across several frequency bins. van der Klis (1989) and Papadakis & Lawrence (1993) also describe the use of binned or smoothed power spectra in the analysis of X-ray binaries and active galaxies exhibiting quasi-periodic oscillations. A similar method can be applied to candidate QPPs, for example those with exponential or Gaussian damping typical for solar and stellar flares (Pugh et al. 2016; Cho et al. 2016), or those with a small variation of the period. These QPP signals may appear as a broad peak in the power spectrum, with a power spanning more than one frequency bin, and hence considering all of the power contained within the peak (by rebinning the spectrum) rather than separately considering the power in each of the frequency bins will give a better assessment of the significance of the peak.

When summing every n frequency bins, the probability density follows a chi-squared 2n d.o.f. distribution with a mean equal to 2n (Appourchaux 2003): pχ2n2(x)=xn−1e−x/22nΓ(n),\begin{equation} p_{\chi ^2_{2n}}(x) = \frac{x^{n-1}\mathrm{e}^{-x/2}}{2^n\Gamma (n)}, \end{equation}(14)where Γ is the gamma function. Hence the probability distribution followed by the rebinned power spectrum is: p2nIj(x)=xn−1e−x/2𝒫j2n𝒫jnΓ(n),\begin{equation} p_{2nI_j}(x) = \frac{x^{n-1}\mathrm{e}^{-x/2\mathcal{P}_j}}{2^n\mathcal{P}^n_j\Gamma (n)}, \end{equation}(15)where \hbox{$\mathcal{P}_j$} is the “true” rebinned power spectrum which can be fitted by a power law model, \hbox{$\hat{\mathcal{P}}_j$}. Plugging this equation into Eq. (9), along with Eq. (7) gives: pγj(z)=∫∞0(yz)n−12n𝒫jnΓ(n)2πSj ×exp⎧⎪⎪⎪⎨⎪⎪⎪⎩−(ln[y]−ln[𝒫j])22Sj2−yz2𝒫j⎫⎪⎪⎪⎬⎪⎪⎪⎭dy.\begin{eqnarray} p_{\gamma_j}(z)& = &\int_{0}^{\infty}\frac{(yz)^{n-1}}{2^n\mathcal{P}_j^n\Gamma (n)\sqrt{2\pi\;}S_j} \nonumber\\[1.5mm] &&\quad \times\exp \left\{-\frac{\left(\ln[y] - \ln[\mathcal{P}_j]\right)^2}{2S^2_j} - \frac{yz}{2\mathcal{P}_j}\right\} \mathrm{d}y. \end{eqnarray}(16)Integrating this probability density from γϵj up to infinity and substituting \hbox{$w = y/\mathcal{P}_j$}, as before, gives: Pr{γjˆ>γϵj}=∫γϵj∞∫0∞(wz/2)n−18πSjΓ(n)exp⎧⎪⎪⎨⎪⎪⎩−(ln w)22Sj2−wz2⎫⎪⎪⎬⎪⎪⎭dwdz.\begin{equation} \textit{Pr}\left\{\hat{\gamma _j} \!>\! \gamma _{\epsilon _j}\right\} = \int_{\gamma _{\epsilon _j}}^{\infty}\int_{0}^{\infty}\frac{(wz/2)^{n-1}}{\sqrt{8\pi\;}S_j\Gamma (n)}\exp \left\{-\frac{\left(\ln\,w\right)^2}{2S^2_j} \!-\! \frac{wz}{2}\right\} \mathrm{d}w\mathrm{d}z. \end{equation}(17)By swapping the order of integration and letting u = wz/ 2 (hence dz = 2du/w), this equation becomes: Pr{γjˆ>γϵj}=∫∞028πSjΓ(n)wexp⎧⎪⎪⎨⎪⎪⎩−(ln w)22Sj2⎫⎪⎪⎬⎪⎪⎭ ×⎧⎪⎪⎨⎪⎪⎩∫wγϵj/2∞exp(−u)un−1du⎫⎪⎪⎬⎪⎪⎭dw,\begin{eqnarray} \textit{Pr}\left\{\hat{\gamma _j} > \gamma _{\epsilon _j}\right\} &=& \int_{0}^{\infty}\frac{2}{\sqrt{8\pi\;}S_j\Gamma (n)w}\exp \left\{-\frac{\left(\ln\,w\right)^2}{2S^2_j}\right\} \nonumber\\[1.5mm] & &\quad \times \left\{\int_{w\gamma _{\epsilon _j}/2}^{\infty} \exp (-u)u^{n-1} \mathrm{d}u\right\}\mathrm{d}w, \end{eqnarray}(18)which, after writing the internal integral in gamma function notation, becomes: Pr{γjˆ>γϵj}=∫0∞12πSjwexp⎧⎪⎪⎨⎪⎪⎩−(ln w)22Sj2⎫⎪⎪⎬⎪⎪⎭Γ(n,wγϵj/2)Γ(n)dw,\begin{equation} \textit{Pr}\left\{\hat{\gamma _j} > \gamma _{\epsilon _j}\right\} = \int_{0}^{\infty}\frac{1}{\sqrt{2\pi\;}S_jw} \exp \left\{-\frac{\left(\ln\,w\right)^2}{2S^2_j}\right\} \frac{\Gamma (n, w\gamma _{\epsilon _j}/2)}{\Gamma (n)}\mathrm{d}w, \end{equation}(19)where Γ(n,wγϵj/ 2) is the upper incomplete gamma function. Like before, this can be solved numerically by equating to Eq. (13), and the confidence level in log space is equal to \hbox{$\log[\hat{\mathcal{P}}_j] + \log[\gamma _{\epsilon _j}\langle I_j - \hat{\mathcal{P}}_j\rangle /2n]$}.

4. Testing the methods on simulated data

Figure 3 shows examples where confidence levels have been calculated for synthetic flare time series with QPP signals, and shows how different background trends (which are unknown for real flare data) affect the appearance of a QPP signal in the power spectrum. To create the time series, a polynomial background trend was added to an exponentially damped sinusoid, white noise, and additional red noise. The additional red noise was generated by a random walk, where each value in the time series is equal to a random number summed with the preceding value. The parameters were chosen to be comparable to flare time series data from the Nobeyama Radioheliograph (see Sect. 5) and, for all of these time series, the sinusoid, white noise and random walk noise terms were kept identical. The top two rows of Fig. 3 show time series with different polynomial background trends on the left, and the corresponding power spectra on the right. Despite the different background trends, peaks corresponding to the sinusoidal signals can be seen above the 99% confidence level in the power spectra. The bottom two rows show the same signals as the top two rows, but instead they have higher amplitude background trends. The steeper background trends mean that the sinusoidal signals are no longer seen at a significant level in the power spectra. Therefore, although the method described in Sect. 2 is useful for testing for the presence of a QPP signal when there is some unknown background trend in the data, when the amplitude of the background trend is sufficiently greater than the amplitude of a QPP signal, the QPP signal will be hidden in the power spectrum.

thumbnail Fig. 3

Examples of synthetic flare time series with QPPs are given on the left, and on the right are the corresponding power spectra, where the red solid line is a power law fit, and the red dotted and dashed lines correspond to the 95% and 99% confidence levels respectively. The top two rows show two signals with different background trends, both with a peak above the 99% level in the power spectrum. The bottom two rows show the same signals but with steeper background trends, the result of which is that the peaks no longer reach significant levels in the power spectra.

A scenario where the method described in Sect. 3 results in a peak above the 99% level in the power spectrum, while the method from Sect. 2 does not, is demonstrated in Fig. 4. The left panel shows a signal with the same background and noise as the signal in the top left panel of Fig. 3, but instead it has a sinusoidal term with a frequency that has a constant mean but a small amount of random variation with time. The result of this is that in the power spectrum, shown in the middle panel of Fig. 4, the peak corresponding to the sinusoidal signal is spread across more than one frequency bin, and hence no longer reaches a significant power level. The rebinned power spectrum is shown in the right panel, where the powers in every two frequency bins of the original spectrum have been summed together. Here the peak corresponding to the sinusoidal signal has a power above the 99% confidence level.

thumbnail Fig. 4

Example of how rebinning can help spectral peaks corresponding to certain kinds of periodic signals reach a significant power level. The left panel shows a synthetic time series signal, similar to that in the top left panel of Fig. 3, but with a sinusoidal component that has a frequency that fluctuates slightly with time. The middle panel shows the corresponding power spectrum, and the right panel shows the rebinned power spectrum (after summing the powers in every two frequency bins). As before, the red solid line is a power law fit, and the red dotted and dashed lines correspond to the 95% and 99% confidence levels respectively.

thumbnail Fig. 5

Left: a section of a GOES C7.1 class flare observed by Nobeyama Radioheliograph. Right: the corresponding power spectrum, where the red solid line is a power law fit to the spectrum, the red dotted line represents the 95% confidence level, and the red dashed line the 99% level. One peak is above the 99% level, at a period of 10.1-0.5+0.6\hbox{$10.1^{+0.6}_{-0.5}$} s.

thumbnail Fig. 6

As in Fig. 5, but for a different flare with a class of M8.7. Although this flare appears to have pulsations there is no peak close to the 95% level in the power spectrum.

thumbnail Fig. 7

As in Fig. 5, but for a different flare with a class of C3.6. Here the power spectrum contains a broad peak, which does not reach the 95% confidence level.

5. Examples of application to solar flare data

An important consideration when performing a periodogram analysis of time series data to search for a periodic signal that is localised in time (such as a wave train), is the choice of start and end times. For example, if a 5 h section of a light curve was taken which contained a flare with QPPs that persisted for 30 min, then the significance of the peak in the power spectrum corresponding to the QPP signal would be lower than if a 30 min section of the same light curve, centred around the QPP signal, was used instead. For the following solar flare light curves, the start and end times were chosen manually to best show off candidate QPPs: this was the only aspect of the analysis that was handled manually. An additional issue with the spectral analysis of time series data is the finite duration of the data. For a chosen section of a flare time series, the start and end values are unlikely to be equal to zero, hence there will be spectral leakage into side-lobes when performing some form of discrete Fourier transform, and this reduces the power of a peak in the power spectrum. One way to avoid this effect is to apply a window function, for example a Hann window, to the time series data before calculating the power spectrum. The application of such a window function is, however, not always helpful when searching for low-amplitude transient signals such as QPPs, since any QPP signals will be suppressed near the start and end of the time series, making detection even more challenging. In addition, the application of any window function other than a rectangular window will alter the distribution of noise in the data, and therefore this would need to be taken into account when using the methods described in this paper. Hence for the following examples, no window function has been applied.

An example of the method described in Sect. 2 being used to confirm candidate QPPs in a GOES C7.1 class flare, observed between 2014 October 29 23:40 and October 30 00:34 UT, is shown in Fig. 5. A section of 17 GHz microwave correlation signal from the Nobeyama Radioheliograph (NoRH; Nakajima et al. 1994) is shown in the left hand panel, and the corresponding Lomb-Scargle periodogram power spectrum is shown on the right. As mentioned in Sect. 2, the correlation data uncertainty was estimated by taking the standard deviation of a flat section of data, taken from 2016 October 27 00:00 until 05:00 UT. This gave an uncertainty of 1.1911749 × 10-5. In the periodogram there is a peak with a period of 10.1-0.5+0.6\hbox{$10.1^{+0.6}_{-0.5}$} s above the 99% confidence level, where the upper and lower uncertainties are taken to be plus or minus half a frequency bin, respectively, on either side of the peak frequency. Visual inspection of the light curve on the left confirms that the pulses have a time spacing approximately equal to this period, hence this can be considered a strong QPP candidate. An example where the method fails to support the possible presence of QPPs in a M8.7 class flare, observed between 2014 October 21 08:09 and 08:15 UT, is shown in Fig. 6, where there is no peak in the power spectrum above the 95% level. Although pulsations can be seen in the light curve, these are small in amplitude when compared to the underlying trend in the data, meaning that the trend dominates in the power spectrum and even though the pulsations may be periodic they are not detectable at a significant level.

Another point to note is that although a broken power law model was used to fit these NoRH power spectra in Figs. 5 and 6, the break cannot be seen. This can be explained by considering the white noise amplitude in the NoRH time series data, which is very small. Hence the frequency at which the white noise would start dominating over the red noise in the power spectrum is likely beyond the range of frequencies included in the spectrum.

thumbnail Fig. 8

Rebinned power spectrum for the flare shown in Fig. 7. The peak at a period of 15-3+5\hbox{$15^{+5}_{-3}$} s now surpasses the 95% confidence level, which is shown by the red dotted line.

For the section of the C3.6 class flare observed between 2014 October 24 03:56 and 04:30 UT, shown in Fig. 7, again the method described in Sect. 2 fails to show a peak in the power spectrum above the 95% level, however a broad peak can be seen. The rebinned power spectrum is shown in Fig. 8, where every three frequency bins have been summed over. Applying the method described in Sect. 3 shows that there is a peak above the 95% confidence level at a period of 15-3+5\hbox{$15^{+5}_{-3}$} s, hence this flare can be considered to have candidate QPPs.

6. Summary

We have demonstrated how the method of Vaughan (2005) can be applied to flare time series data in order to test for the presence of QPPs, subject to the careful choice of start and end times. The method has been adapted to be used with rebinned power spectra, which can aid the detection of QPP signals that have a period that varies slightly, or have a modulated amplitude. These methods avoid detrending the data, an approach which has been shown to have the potential to lead to false detections when the detrending is done by smoothing (e.g. Gruber et al. 2011; Auchère et al. 2016). An alternative method which also avoids detrending has been proposed by Inglis et al. (2015). This method involves doing a model comparison; different models, such as a power law and a power law plus Gaussian peak, are fitted to the power spectrum and are compared by calculating the Bayesian information criterion (BIC) for each. The model with the smaller BIC value will then be the preferred model. We therefore suggest that a thorough search for QPPs in flares could make use of both approaches. These methods, however, may not be suitable for searching for non-stationary QPP signals if the change in period is too great, in which case either some form of wavelet analysis (e.g. Kupriyanova et al. 2010) or empirical mode decomposition (e.g. Kolotkov et al. 2015) could be used.

Acknowledgments

We would like to thank the anonymous referee for their helpful comments. C.E.P. would like to thank Alex Seaton and Mark Hollands for useful discussions. C.E.P. and V.M.N.: This work was supported by the European Research Council under the SeismoSun Research Project No. 321141. A.-M.B. thanks the Institute of Advanced Study, University of Warwick for their support. The authors would like to acknowledge the ISSI International Team led by A.-M.B. for many productive discussions, and are also grateful to the Nobeyama team for providing the data used.

References

  1. Anfinogentov, S., Nakariakov, V. M., Mathioudakis, M., Van Doorsselaere, T., & Kowalski, A. F. 2013, ApJ, 773, 156 [NASA ADS] [CrossRef] [Google Scholar]
  2. Appourchaux, T. 2003, A&A, 412, 903 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  3. Appourchaux, T. 2004, A&A, 428, 1039 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Aschwanden, M. J., Kliem, B., Schwarz, U., et al. 1998, ApJ, 505, 941 [NASA ADS] [CrossRef] [Google Scholar]
  5. Auchère, F., Froment, C., Bocchialini, K., Buchlin, E., & Solomon, J. 2016, ApJ, 825, 110 [NASA ADS] [CrossRef] [Google Scholar]
  6. Chaplin, W. J., Elsworth, Y., Isaak, G. R., et al. 2002, MNRAS, 336, 979 [NASA ADS] [CrossRef] [Google Scholar]
  7. Cho, I.-H., Cho, K.-S., Nakariakov, V. M., Kim, S., & Kumar, P. 2016, ApJ, 830, 110 [NASA ADS] [CrossRef] [Google Scholar]
  8. Chowdhury, P., Srivastava, A. K., Dwivedi, B. N., Sych, R., & Moon, Y.-J. 2015, Adv. Space Res., 56, 2769 [NASA ADS] [CrossRef] [Google Scholar]
  9. Curtiss, J. H. 1941, Ann. Math. Stat., 12, 409 [CrossRef] [Google Scholar]
  10. Davenport, J. R. A., Hawley, S. L., Hebb, L., et al. 2014, ApJ, 797, 122 [NASA ADS] [CrossRef] [Google Scholar]
  11. Gruber, D., Lachowicz, P., Bissaldi, E., et al. 2011, A&A, 533, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Hayes, L. A., Gallagher, P. T., Dennis, B. R., et al. 2016, ApJ, 827, L30 [NASA ADS] [CrossRef] [Google Scholar]
  13. Horne, J. H., & Baliunas, S. L. 1986, ApJ, 302, 757 [NASA ADS] [CrossRef] [Google Scholar]
  14. Huang, J., Tan, B., Zhang, Y., Karlický, M., & Mészárosová, H. 2014, ApJ, 791, 44 [NASA ADS] [CrossRef] [Google Scholar]
  15. Inglis, A. R., Ireland, J., & Dominique, M. 2015, ApJ, 798, 108 [NASA ADS] [CrossRef] [Google Scholar]
  16. Inglis, A. R., Ireland, J., Dennis, B. R., Hayes, L., & Gallagher, P. 2016, ApJ, 833, 284 [NASA ADS] [CrossRef] [Google Scholar]
  17. Kolotkov, D. Y., Nakariakov, V. M., Kupriyanova, E. G., Ratcliffe, H., & Shibasaki, K. 2015, A&A, 574, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Kumar, P., Nakariakov, V. M., & Cho, K.-S. 2016, ApJ, 822, 7 [NASA ADS] [CrossRef] [Google Scholar]
  19. Kupriyanova, E. G., Melnikov, V. F., Nakariakov, V. M., & Shibasaki, K. 2010, Sol. Phys., 267, 329 [NASA ADS] [CrossRef] [Google Scholar]
  20. Kupriyanova, E. G., Melnikov, V. F., & Shibasaki, K. 2013, Sol. Phys., 284, 559 [NASA ADS] [CrossRef] [Google Scholar]
  21. Lipa, B. 1978, Sol. Phys., 57, 191 [NASA ADS] [CrossRef] [Google Scholar]
  22. Mathioudakis, M., Seiradakis, J. H., Williams, D. R., et al. 2003, A&A, 403, 1101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. McAteer, R. T. J., Young, C. A., Ireland, J., & Gallagher, P. T. 2007, ApJ, 662, 691 [NASA ADS] [CrossRef] [Google Scholar]
  24. McAteer, R. T. J., Aschwanden, M. J., Dimitropoulou, M., et al. 2016, Space Sci. Rev., 198, 217 [NASA ADS] [CrossRef] [Google Scholar]
  25. Myagkova, I. N., Bogomolov, A. V., Kashapova, L. K., et al. 2016, Sol. Phys., 291, 3439 [NASA ADS] [CrossRef] [Google Scholar]
  26. Nakajima, H., Nishio, M., Enome, S., et al. 1994, IEEE Proc., 82, 705 [NASA ADS] [CrossRef] [Google Scholar]
  27. Nakariakov, V. M., & Melnikov, V. F. 2009, Space Sci. Rev., 149, 119 [NASA ADS] [CrossRef] [Google Scholar]
  28. Nakariakov, V. M., Pilipenko, V., Heilig, B., et al. 2016, Space Sci. Rev., 200, 75 [NASA ADS] [CrossRef] [Google Scholar]
  29. Papadakis, I. E., & Lawrence, A. 1993, MNRAS, 261, 612 [NASA ADS] [CrossRef] [Google Scholar]
  30. Parks, G. K., & Winckler, J. R. 1969, ApJ, 155, L117 [NASA ADS] [CrossRef] [Google Scholar]
  31. Pugh, C. E., Armstrong, D. J., Nakariakov, V. M., & Broomhall, A.-M. 2016, MNRAS, 459, 3659 [NASA ADS] [CrossRef] [Google Scholar]
  32. Reznikova, V. E., & Shibasaki, K. 2011, A&A, 525, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Scargle, J. D. 1982, ApJ, 263, 835 [NASA ADS] [CrossRef] [Google Scholar]
  34. Simões, P. J. A., Hudson, H. S., & Fletcher, L. 2015, Sol. Phys., 290, 3625 [NASA ADS] [CrossRef] [Google Scholar]
  35. Tian, H., Young, P. R., Reeves, K. K., et al. 2016, ApJ, 823, L16 [NASA ADS] [CrossRef] [Google Scholar]
  36. van der Klis, M. 1989, ARA&A, 27, 517 [NASA ADS] [CrossRef] [Google Scholar]
  37. Van Doorsselaere, T., Kupriyanova, E. G., & Yuan, D. 2016, Sol. Phys., 291, 3143 [NASA ADS] [CrossRef] [Google Scholar]
  38. Vaughan, S. 2005, A&A, 431, 391 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

All Figures

thumbnail Fig. 1

Plots of the integrand in Eq. (12) as a function of w. The solid black, dotted red and dashed blue lines show the function when Sj is equal to 0.4, 0.2 and 0.02 respectively. The value of γϵj has arbitrarily been chosen to be equal to 20, which is a typical value.
In the text
thumbnail Fig. 2

Plots of Eqs. (12) and (13) as a function of γϵj are shown by the solid red and dashed blue lines respectively. The values of Sj and ϵN′/N′ have arbitrarily been chosen to be equal to 0.2 and 0.01/100 respectively. The solution we require is where Eq. (12) is equal to Eq. (13), which corresponds to γϵj = 21.467.
In the text
thumbnail Fig. 3

Examples of synthetic flare time series with QPPs are given on the left, and on the right are the corresponding power spectra, where the red solid line is a power law fit, and the red dotted and dashed lines correspond to the 95% and 99% confidence levels respectively. The top two rows show two signals with different background trends, both with a peak above the 99% level in the power spectrum. The bottom two rows show the same signals but with steeper background trends, the result of which is that the peaks no longer reach significant levels in the power spectra.
In the text
thumbnail Fig. 4

Example of how rebinning can help spectral peaks corresponding to certain kinds of periodic signals reach a significant power level. The left panel shows a synthetic time series signal, similar to that in the top left panel of Fig. 3, but with a sinusoidal component that has a frequency that fluctuates slightly with time. The middle panel shows the corresponding power spectrum, and the right panel shows the rebinned power spectrum (after summing the powers in every two frequency bins). As before, the red solid line is a power law fit, and the red dotted and dashed lines correspond to the 95% and 99% confidence levels respectively.
In the text
thumbnail Fig. 5

Left: a section of a GOES C7.1 class flare observed by Nobeyama Radioheliograph. Right: the corresponding power spectrum, where the red solid line is a power law fit to the spectrum, the red dotted line represents the 95% confidence level, and the red dashed line the 99% level. One peak is above the 99% level, at a period of 10.1-0.5+0.6\hbox{$10.1^{+0.6}_{-0.5}$} s.
In the text
thumbnail Fig. 6

As in Fig. 5, but for a different flare with a class of M8.7. Although this flare appears to have pulsations there is no peak close to the 95% level in the power spectrum.
In the text
thumbnail Fig. 7

As in Fig. 5, but for a different flare with a class of C3.6. Here the power spectrum contains a broad peak, which does not reach the 95% confidence level.
In the text
thumbnail Fig. 8

Rebinned power spectrum for the flare shown in Fig. 7. The peak at a period of 15-3+5\hbox{$15^{+5}_{-3}$} s now surpasses the 95% confidence level, which is shown by the red dotted line.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.