Significance testing for quasi-periodic pulsations in solar and stellar flares

TL;DR

This paper adapts a statistical method to detect quasi-periodic pulsations in solar and stellar flare data, effectively accounting for red noise and uncertainties, and demonstrates its application on synthetic and real data.

Contribution

It introduces an improved significance testing method for QPP detection that handles red noise and is applicable to rebinned spectra without detrending.

Findings

Method successfully detects QPPs in synthetic data.

Application to real flare data reveals significant QPP peaks.

No detrending needed for effective QPP detection.

Abstract

The robust detection of quasi-periodic pulsations (QPPs) in solar and stellar flares has been the topic of recent debate. In light of this, we have adapted a method described by Vaughan (2005) to aid with the search for QPPs in flare time series data. The method identifies statistically significant periodic signals in power spectra, and properly accounts for red noise as well as the uncertainties associated with the data. We show how the method can be further developed to be used with rebinned power spectra, allowing us to detect QPPs whose signal is spread over more than one frequency bin. An advantage of these methods is that there is no need to detrend the data prior to creating the power spectrum. Examples are given where the methods have been applied to synthetic data, as well as real flare time series data with candidate QPPs from the Nobeyama Radioheliograph. These show that, despite the transient nature of QPPs, peaks corresponding to the QPPs can be seen at a significant level in the power spectrum without any form of detrending or other processing of the original time series data, providing the background trends are not too steep.

Equations

$$p_{{\chi }_2^2}(x)=\frac{1}{2}{\mathrm{e}}^{-x/2}.$$

$$\text{Pr}\left\{X>x^{\prime }\right\}={\int }_{x^{\prime }}^{\infty }\frac{1}{2}{\mathrm{e}}^{-x/2}\mathrm{d}x=e^{-x^{\prime }/2}.$$

$$\text{Pr}\left\{X>x^{\prime }\right\}=1-{\left(1-{ϵ}_{N^{\prime }}\right)}^{1/N^{\prime }}\approx {ϵ}_{N^{\prime }}/N^{\prime },$$

$$I(f_j)=\mathcal{P}(f_j){\chi }_2^2/2,$$

$$\log \left[\hat{\mathcal{P}}(f)\right]=\{\begin{array}{ll}-\alpha \log \left[f\right]+c& \text{if }f<f_{break}\\ -\left(\alpha -\beta \right)\log \left[f_{break}\right]-\beta \log \left[f\right]+c& \text{if }f>f_{break},\end{array}$$

$$p_{2I_j}(x)=\frac{1}{2{\mathcal{P}}_j}{\mathrm{e}}^{-x/2{\mathcal{P}}_j}.$$

$$p_{{\hat{\mathcal{P}}}_j}(y)=\frac{1}{\sqrt{2\pi }y{S}_j}\exp \left\{-\frac{{\left(\ln [y]-\ln [{\mathcal{P}}_j]\right)}^2}{2S_j^2}\right\},$$

$$S_j=\text{err}\left\{\log \left[\hat{\mathcal{P}}(f_j)\right]\right\}\times \ln [10],$$

$$p_{{\gamma }_j}(z)={\int }_0^{+\infty }|y|p_{2I_j}(zy)p_{{\mathcal{P}}_j}(y)\mathrm{d}y,$$

$$\text{Pr}\left\{\widehat{{\gamma }_j}>{\gamma }_{{ϵ}_j}\right\}=\frac{1}{\sqrt{8\pi }{S}_j{\mathcal{P}}_j}{\int }_{{\gamma }_{{ϵ}_j}}^{\infty }{\int }_0^{\infty }\exp \left\{-\frac{{\left(\ln [y]-\ln [{\mathcal{P}}_j]\right)}^2}{2S_j^2}-\frac{yz}{2{\mathcal{P}}_j}\right\}\mathrm{d}y\mathrm{d}z.$$

$$\text{Pr}\left\{\widehat{{\gamma }_j}>{\gamma }_{{ϵ}_j}\right\}=\frac{1}{\sqrt{8\pi }{S}_j}{\int }_{{\gamma }_{{ϵ}_j}}^{\infty }{\int }_0^{\infty }\exp \left\{-\frac{{\left(\ln w\right)}^2}{2S_j^2}-\frac{wz}{2}\right\}\mathrm{d}w\mathrm{d}z.$$

$$\text{Pr}\left\{\widehat{{\gamma }_j}>{\gamma }_{{ϵ}_j}\right\}={\int }_0^{\infty }\frac{1}{\sqrt{2\pi }{S}_{j}w}\exp \left\{-\frac{{\left(\ln w\right)}^2}{2S_j^2}-\frac{{\gamma }_{{ϵ}_j}w}{2}\right\}\mathrm{d}w,$$

$$\text{Pr}\left\{\widehat{{\gamma }_j}>{\gamma }_{{ϵ}_j}\right\}\approx \frac{{ϵ}_{N^{\prime }}}{N^{\prime }},$$

$$p_{{\chi }_{2n}^2}(x)=\frac{x^{n-1}{\mathrm{e}}^{-x/2}}{2^n\Gamma (n)},$$

$$p_{2n{I}_j}(x)=\frac{x^{n-1}{\mathrm{e}}^{-x/2{\mathcal{P}}_j}}{2^n{\mathcal{P}}_j^n\Gamma (n)},$$

$$p_{{\gamma }_j}(z)={\int }_0^{\infty }\frac{(yz{)}^{n-1}}{2^n{\mathcal{P}}_j^n\Gamma (n)\sqrt{2\pi }{S}_j}\exp \left\{-\frac{{\left(\ln [y]-\ln [{\mathcal{P}}_j]\right)}^2}{2S_j^2}-\frac{yz}{2{\mathcal{P}}_j}\right\}\mathrm{d}y.$$

$$\text{Pr}\left\{\widehat{{\gamma }_j}>{\gamma }_{{ϵ}_j}\right\}={\int }_{{\gamma }_{{ϵ}_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_j^2}-\frac{wz}{2}\right\}\mathrm{d}w\mathrm{d}z.$$

$$\text{Pr}\left\{\widehat{{\gamma }_j}>{\gamma }_{{ϵ}_j}\right\}={\int }_0^{\infty }\frac{2}{\sqrt{8\pi }{S}_j\Gamma (n)w}\exp \left\{-\frac{{\left(\ln w\right)}^2}{2S_j^2}\right\}\left\{{\int }_{w{\gamma }_{{ϵ}_j}/2}^{\infty }\exp (-u)u^{n-1}\mathrm{d}u\right\}\mathrm{d}w,$$

$$\text{Pr}\left\{\widehat{{\gamma }_j}>{\gamma }_{{ϵ}_j}\right\}={\int }_0^{\infty }\frac{1}{\sqrt{2\pi }{S}_{j}w}\exp \left\{-\frac{{\left(\ln w\right)}^2}{2S_j^2}\right\}\frac{\Gamma (n,w{\gamma }_{{ϵ}_j}/2)}{\Gamma (n)}\mathrm{d}w,$$

Full text

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\propto f^{-\alpha}$). 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:

[TABLE]

The probability of having a value $X$ that is greater than some threshold $x^{\prime}$ is therefore:

[TABLE]

Since the power spectrum is positive everywhere, $x$ and $x^{\prime}$ in the above equation are always positive. When considering a power spectrum sampled at $N^{\prime}=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:

[TABLE]

where $\epsilon_{N^{\prime}}$ is the false alarm probability, and the approximation holds when is $\epsilon_{N^{\prime}}$ 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^{\prime}$. 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 Section 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 Section 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 Section 4, and in Section 5 a few examples are given where the methods have been applied to real solar flare data. A summary is given in Section 6.

2 Confidence levels on power-law power spectra

From Vaughan (2005), if $I(f_{j})$ is the periodogram power at a particular frequency, $f_{j}$, we have:

[TABLE]

where $\mathcal{P}(f_{j})$ is the ‘true’ power, and $\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:

[TABLE]

where $f_{break}$ is the frequency at which the power law break occurs, $\alpha$ and $\beta$ are power law indices and $c$ is a constant. The probability density for ${2I_{j}}$ (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:

[TABLE]

When fitting the ‘true’ power spectrum, $\mathcal{P}$, with a broken power law model, $\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:

[TABLE]

with

[TABLE]

where $\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 Section 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 Section 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 Section 5 parameters of $\alpha=2$, $\beta=1$, $c=-3$ and $f_{break}=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.

The probability density of the ratio $\hat{\gamma}_{j}=2I_{j}/\hat{\mathcal{P}}_{j}$ (the power spectrum with the red noise component removed) can be found from (Curtiss, 1941):

[TABLE]

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 $\gamma_{\epsilon_{j}}$ and infinity gives the probability that a value $\hat{\gamma_{j}}$ is greater than $\gamma_{\epsilon_{j}}$:

[TABLE]

Substituting in another dummy variable, $w=y/\mathcal{P}_{j}$ (with $\mathrm{d}y=\mathcal{P}_{j}\;\mathrm{d}w$), simplifies this to:

[TABLE]

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 $S_{j}$). Then the function can be integrated with respect to $z$, to get:

[TABLE]

which can be equated to (see Eq. (3) in Section 1):

[TABLE]

and solved numerically in order to find a $\gamma_{\epsilon_{j}}$ corresponding to each value of the fitted model power spectrum, $\hat{\mathcal{P}}_{j}$. Figure 2 shows plots of Eq. (12) and Eq. (13) as a function of $\gamma_{\epsilon_{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. $\langle\log\mathcal{P}_{j}\rangle\neq\log\langle\mathcal{P}_{j}\rangle$). 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 $\langle\mathcal{P}_{j}-\hat{\mathcal{P}}_{j}\rangle$) can be found. Therefore the confidence level in log space is then equal to $\log[\hat{\mathcal{P}}_{j}]+\log[\gamma_{\epsilon_{j}}\langle\mathcal{P}_{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):

[TABLE]

where $\Gamma$ is the gamma function. Hence the probability distribution followed by the rebinned power spectrum is:

[TABLE]

where $\mathcal{P}_{j}$ is the ‘true’ rebinned power spectrum which can be fitted by a power law model, $\hat{\mathcal{P}}_{j}$. Plugging this equation into Eq. (9), along with Eq. (7) gives:

[TABLE]

Integrating this probability density from $\gamma_{\epsilon_{j}}$ up to infinity and substituting $w=y/\mathcal{P}_{j}$, as before, gives:

[TABLE]

By swapping the order of integration and letting $u=wz/2$ (hence $\mathrm{d}z=2\mathrm{d}u/w$), this equation becomes:

[TABLE]

which, after writing the internal integral in gamma function notation, becomes:

[TABLE]

where $\Gamma(n,w\gamma_{\epsilon_{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 $\log[\hat{\mathcal{P}}_{j}]+\log[\gamma_{\epsilon_{j}}\langle\mathcal{P}_{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 Section 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 Section 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.

A scenario where the method described in Section 3 results in a peak above the 99% level in the power spectrum, while the method from Section 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.

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 hour section of a light curve was taken which contained a flare with QPPs that persisted for 30 minutes, then the significance of the peak in the power spectrum corresponding to the QPP signal would be lower than if a 30 minute 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 Section 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 Section 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\times 10^{-5}$. In the periodogram there is a peak with a period of $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.

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 Section 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 Section 3 shows that there is a peak above the 95% confidence level at a period of $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.

Acknowledgements.

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. & 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.