© ESO 2018

1 Introduction

The discovery of quasi-periodic oscillations (QPOs) in the giant flare of the magnetar SGR 1806-20 by Israel et al. (2005) may have been the first detection of neutron star oscillations and triggered a wealth of theoretical work explaining the reported frequencies. The giant flare was likely caused by a large-scale reconnection or an interchange instability of the magnetic field (Thompson & Duncan 1995). Large amounts of energy are released as an expanding e^±-pair plasma, observable as the initial spike of the giant flare. Parts of this plasma are trapped by the ultra-strong magnetic field and form a so-called trapped fireball (Thompson & Duncan 1995), which then slowly evaporates on a timescale of up to a few 100 s. The QPOs were detected in this decaying tail of the giant flare. Other groups have not only confirmed this detection, they even found additional oscillation frequencies in different magnetars: 18, 26, 29, 92, 150 , 625, and 1840 Hz in the giant flare of SGR 1806-20, and 28, 53, 84, and 155 Hz in the giant flare of SGR 1900+14 (e.g. Strohmayer & Watts 2005; 2006; Watts & Strohmayer 2006; Huppenkothen et al. 2014c). With different methods, more oscillation frequencies were found in the giant flare of SGR 1806-20 by Hambaryan et al. (2011) 17, 21, 36, 59, and 116 Hz. The number of giant flares at the time of writing is limited to three events. The more frequent but less energetic bursts of several magnetars have therefore also been investigated for frequencies, and some candidates were found: 57 Hz in SGR 1806-20 (Huppenkothen et al. 2014b), and at 93, 127, and 260 Hz in SGR J1550–5418 (Huppenkothen et al. 2014a).

It was soon realized that these frequencies are probably related to oscillations of the neutron star, and several groups tried to identify them as elastic oscillations of the crust (Duncan 1998; Messios et al. 2001; Strohmayer & Watts 2005; Piro 2005; Sotani et al. 2007; 2013, 2016; Samuelsson & Andersson 2007, 2009; Steiner & Watts 2009; Deibel et al. 2014), Alfvén oscillations (Cerdá-Durán et al. 2009; Sotani et al. 2008; Colaiuda et al. 2009), or coupled magneto-elastic oscillations (Levin 2006; 2007; Glampedakis et al. 2006; Gabler et al. 2011; 2012; Colaiuda & Kokkotas 2011; van Hoven & Levin 2011, 2012). The theoretical models based on the observed frequencies are very elaborate and may be able to constrain properties of high-density matter as found in the interior of neutron stars. Some of the models, for instance, require a superfluid component in the core of the star (van Hoven & Levin 2011, 2012; Glampedakis et al. 2011; Passamonti & Lander 2013, 2014; Gabler et al. 2013; 2016, and in prep.). Different models depend sensitively on the identification of the fundamental oscillation frequency, and may not explain all of the observed frequencies. Even when the fundamental frequency is identified, the interpretation and parameter estimation is not yet straightforward because of degeneracies in the parameter space. However, keeping other stellar parameters fixed, some general trends of the fundamental oscillation frequency can be summarized as follows (see Gabler et al. 2016, and in prep. for a detailed discussion): i) The frequency scales linearly with the magnetic field strength. ii) It decreases with increasing compactness (Sotani et al. 2008). The compactness is related to the hardness of the equation of state (EOS): Material with a stiff equation of state is harder to compress, leading to larger radii and hence lower compactnesses. iii) It can only reach the surface for significantly strong magnetic fields , whose thresholds depend on the square root of the shear velocity (Gabler et al., in prep.).

It is of great importance to understand which of the frequencies are in the signal. In previous attempts at identifying possible frequencies in the light curve, the statistical noise of the detectors was not modelled consistently. To improve on this, we employ a Bayesian method that properly takes the photon shot noise as the largest contributor into account.

In this work, we re-analyse the data for the two giant flares of SGR1806-20 and SGR1900+14, which were obtained with the proportional counter array (PCA) of the Rossi X-rayTiming Explorer (RXTE). The data are available online at the High Energy Astrophysics Science Archive Research Center (HEASARC). In Sect. 2 we briefly discuss our numerical approach to estimating the influence of the noise and to reconstructing the likely signal with a reduced contribution from the noise. The next section is devoted to the reconstruction and investigation of the light curves of the two giant flares of SGR 1806-20 and SGR 1900+14, which are then discussed in Sect. 4. We conclude our analysis in Sect. 5.

2 Inference of photon observations

Our goal is to reconstruct the light curve and possible frequencies of QPOs in giant X-ray flares for the neutron stars SGR 1806-20 and SGR 1900+14 using the data taken by RXTE. Owing to experimental constraints such as limited data storage, finite time resolution of each photon count, finite detector size, and sensitivity, RXTE cannot record all physical relevant and available information of a continuous photon flux. Most importantly, the recorded photon counts contain significant photon shot noise. Hence, we have to use probabilisitic data analysis methods to obtain an estimate of the time-dependent, continuous photon flux ϕ(t), including its uncertainties. Since the flux varies on a logarithmic scale, we reconstructed the logarithm of the flux s(t) = log(ϕ(t)/ϕ₀ and its temporal power spectrum directly from the data.

In this context, it is natural to build the inference upon the Bayes theorem, allowing us to investigate the a posteriori probability distribution P(ϕ|d), stating how likely a given photon flux field ϕ is given the data set d . For the two data sets of SGR 1806-20 and SGR 1900+14, we assumed the data d to be the result of a Poisson process whose expectation value λ is given by the photon flux ϕ(t) of the photon burst convolved with the response operator R , (1)

The response operator encodes all instrument specifications and states how ϕ imprints itself on λ. Here, we only considered the readout deadtime periods of the instrument and neglected all other instrument responses. In Eq. (1) the constant ϕ₀ absorbs numerical constants and the physical units of the time dependent log- flux signal field s(t) = s. As the signal describes the logarithmic flux, we naturally ensured the positive definiteness of the photon flux. Since we analysed QPOs, it is expected that ϕ exhibits unknown but spatial correlations. Hence, we did not enforce any concrete spatial correlations. We only assumed ϕ to follow a multivariate log-normal statistic, with an a priori unknown covariance. Assuming stationary statistics, the underlying covariance is fully determined by a power spectrum P_s(v), described by P_s(v)= P₀e^τ(v), to ensure positivity of the power spectrum. We inferred this power spectrum as well as ϕ from the data themselves by setting up a hierarchical prior model (Selig et al. 2015; Junklewitz et al. 2016; Pumpe et al. 2016; Enßlin & Frommert 2011; Enßlin & Knollmüller 2016). Hence the posterior of our Bayesian inference is given by (2) where P(d|ϕ,τ) is the likelihood describing how a potential photon flux ϕ including a certain covariance structure described by τ imprints itself in a potential data set. In addition to strong spectral peaks in the power spectrum induced by the almost discrete frequencies of QPOs, we enforce some kind of spectral smoothness on the logarithmic scale, favouring power-law spectra. This behaviour is enforced and tuned by the prior term P(τ|σ_sm) and its parameter σ_sm (Oppermann et al. 2013). If σ_sm → 0, infinite spectral smoothness is enforced, while σ_sm →∞ does not enforce any spectral smoothness. As we wish to be sensitive to spectral lines in the power spectrum, the influence of σ_sm on the inference is discussed in greater detail in Appendix A.

For the detailed mathematical rigorous derivation and discussion of the used D³PO- algorithm, we refer to Selig & Enßlin (2015) and Enßlin & Knollmüller (2016). The D³ PO- algorithm was successfully applied on the Fermi-LAT data (Selig et al. 2015). In order to handle high power spectrum resolutions as is needed to infer spectral lines in the power spectrum, we reimplemented the algorithm in NIFTy 3 (Steininger et al. 2017). The capabilities of the algorithm to analyse QPOs are demonstrated in the Appendix A to C.

3 Results

3.1 SGR 1806-20

We re-analysed the archival RXTE data of the giant flare of SGR 1806-20, which occurred on 2004 December 27. Owing to the high photon flux, the instrument telemetry was saturated, causing several data gaps in the first seconds of the flare. We therefore neglected the very beginning of the flare in the current analysis and started our investigation roughly 4.5 s after the initial rise, immediately after the third deadtime interval. However, we accounted for the only remaining operational down time of the instrument in our data between 7.3865 s ≤t ≤ 7.9975 s. For the analysis, we binned the data into pixels with a volume of 1∕800 s. To overcome the periodic boundary conditions introduced by the fast-Fourier transformation, which we used to switch between signal space and its harmonic space, we performed the signal inference on a regular grid with 2¹⁹ pixel, each also with a volume of 1∕800 s, by adding sufficient buffer time.

In Fig. 1 we plot the inferred ϕ for the entire duration of the giant flare and for a selected period of time for a smoothness parameter σ_sm = 5 × 10⁵. Our algorithm is able to significantly reduce the scatter of the light curve that is caused by the photon shot noise. Thus the reconstructed light curves can now be further analysed for potential periodic signals.

Figure 2 shows the reconstructed profiles of one pulse rotation period. There, the mean pulse profile is given as a thick black line, and all individual pulses are plotted in different colours. Dark blue are the first pulses that have a significantly lower second maxima around 5.5 s than the mean pulse. The major maximum (time ~ 3 s) is very close to the mean maximum. At intermediate times (green lines), both maxima of the pulse take their maximum values, roughly 40% and 130% more than at the beginning for the first and second maximum, respectively.At late times (red lines) the main peak declines by 40%, while the second maximum almost stays constant and has an amplitude similar to the main peak. This behaviour indicates a complex evolution of the fireball, or of the fireballs, if there are more than one.

The power spectrum of the entire flare is plotted in Fig. 3. The rotation period of the magnetar is recovered with the frequency of the first main peak ν₀ = 0.1323 Hz, which is very close to ν₀ = 0.13249 Hz, as given in Olausen & Kaspi (2014). We are able to find up to the 31st overtone of this frequency at ν = 4.245 Hz. In Fig. 3b we show the reconstructed power spectrum P_rec from Fig. 3a in black together with the power spectrum obtained from the reconstructed light curve (Fig. 1a) P_s(t)(v) in red. At low frequencies, that is, ν ≲ 3 Hz, the two spectra are in good agreement, as the reconstruction is well constrained by the data on large scales. Between 3 Hz ≲ ν ≲ 20 Hz, the inference algorithm enters the regime of a lower signal-to-noise ratio (S/N), which in principle leads to noisier P_rec . However, this natural behaviour is counteracted by the smoothness-enforcing prior. At higher ν ≳ 20 Hz, the shapes of both spectra start to deviate significantly. The reason for this is that in the noise-dominated frequency regime, D³ PO filters out the photon shot noise. From a naive perspective, small-scale features in the signal therefore need to be significantly strong in order to be detectable after a pure photon shot noise filtering operation on the data set. However, D³ PO accounts for the power loss of this filtering when it reconstructs the power spectrum from the data themselves. Thus, Fig. 3b indicates that above 20 Hzthe data are noise dominated, and spectral features there have to be very strong to be recognisable. To test the dependence of our method on the chosen smoothness prior σ_sm as discussed also in Appendix A, we additionally calculated the reconstructed light curve and its corresponding power spectrum for σ_sm = 10⁵. The latter is given in Fig. 3c. Obviously, a smaller σ_sm leads to a smoothing of the spectrum and the algorithm suppresses the detection of periodic signals at higher frequencies. We found σ_sm = 5 × 10⁵ to be the optimal value to still observe power in the Fourier transform at higher frequencies. For higher values of σ_sm, we qualitatively obtain similar but more noisy results for the reconstructed light curve and the corresponding power spectra.

In addition to the obvious peaks that are related to the rotation period and the corresponding overtones, there are still other features in the reconstructed power spectrum of Fig. 3a that seem to have higher powers than the noise. To estimate the significance of these spectral peaks, we calculated a residual χ between the inferred log-spectrum τ and its local median τ̄ weighted with the local variance σ, (3)

The local median and local variance were calculated over a window of 401 pixels, corresponding to a frequency window of approximately 1 Hz. In the top panel of Fig. 5 we plot the histogram of χ₀, where the index 0 refers to the fundamental frequency, that is, χ at the respective frequency. The resulting distribution deviates significantly from a Gaussian, as there is a significant excess for large χ₀. These counts can easily be identified with the highest spectral peaks in Fig. 3a as integer multiples of the neutron star rotation frequency of ν₀ = 0.1323 Hz. The fat tails of the distribution make it hard to identify whether a peak sticks out of the tail, in particular for χ₀ ≳ 10. For χ₀ ≳ 15 all peaks can be identified and are related to the neutron star rotation period, except for one at ν ~ 9.187 Hz. In Table 1 we show all frequencies that have χ₀ > 11 to have a selection of possible oscillation candidates. If there are two or more neighbouring frequencies with χ₀ > 11, we list the highest value.

All other candidates in Table 1 have significantly lower χ than the oscillation at 9.187 Hz and are consistent with being in the tail of the distribution that is shown in the top panel of Fig. 5. This indicates that these are artefacts of the Poisson noise.

We also checked the χ₀ values of previously reported frequencies. None of them reach more than χ₀ ≳ 5. We therefore extended our search in a ± 5% interval around the frequencies in Table 2. The only frequency higher than χ₀ = 11 is at ν = 16.27 Hz. Thus all reported lines are consistent with noise. However, we already see an interesting pattern emerging: In Table 2 we locally (within the ± 5% interval) find the highest powers at 18.265 and 36.412 Hz. These are almost twice and four times the only significant frequency at ν= 9.187 Hz that our method detects beyond the rotational frequency and its first 31 harmonics. In Fig. 4 we plot the reconstructed power spectrum around the ν = 9.187 Hz candidate oscillation (black line) and its first overtone (red line). The amplitude at ν = 9.187 Hz is significantly larger (factor 2) than other amplitudes in the shown frequency range, while the amplitude at ν = 18.265 Hz is comparable with other spectral peaks. As discrete frequencies in a power spectrum are likely to show spectral peaks at integer multiplies of a ground frequency, we also display the two-dimensional histogram of the calculated weighted residual χ₀ at some ground frequency on the x-axis and its first harmonics χ₁ on the y-axis in the bottom panel of Fig. 5. We marked all counts with χ₁ ≥ 5 and χ₀ ≥ 10 with their corresponding frequency in Hz. We find about ten frequencies that satisfy this criterium. Obviously, all but the frequencies around ν ~ 9.186 Hz are integer multiples of the rotation frequency ν₀ .

This further increases our confidence that ν ~ 9.186 Hz is a candidate for an additional periodic signal in the data. We do not find any significant features for frequencies higher than the corresponding overtones that the algorithm recovered at ν ~ 18.265 Hz and ν ~ 36.412 Hz.

Fig. 1 Reconstructed light curve of the giant flare of SGR 1806-2 using a smoothness-enforcing prior, with σsm = 5 × 105. The grey narrow rectangle indicates the operational down time of the instrument. For better visibility, we plot the light curve between ≈164 s and ≈174 s in Figure 1b. In addition to the raw photon counts (black dots), the black line indicates the reconstruction of the expected photon counts, i.e. λ as well as its one-σ confidence interval. Owing to the high resolution of the photon flux, only every fourth point of the regular grid is plotted. In both plots each pixel has a duration of 1∕800 s.

Fig. 2 Pulse profiles of different rotation periods overplotted as given in Fig. 1a. The temporal evolution of the pulse profile is visible from blue to green to red. The mean pulse is shown as a thick black line. All pulses are shifted such that their log-means are zero.

Fig. 3 Reconstructed power spectra of the giant flare of SGR 1806-20: For Figure 3a and Figure 3c we used smoothness-enforcing priors with σsm = 5 × 105 and σsm = 105, respectively. The uncertainty intervals are given as grey shaded areas. In Figure 3b we show the reconstructed power spectrum again from the data themselves as in Figure 3a, along with the power spectrum of the logarithmic reconstructed light curve of Fig. 1a. Because of the high resolution of the reconstructed power spectra, only every fourth point of the regular grid is plotted. Each pixel has a volume of 1∕655 Hz.

Fig. 4 Zoomed-in view of the reconstructed power spectra of giant flare SGR 1806-20 around ν ~ 9.186 Hz in black and its first overtone around ν ~ 18.265 Hz in red. The spectra correspond to σsm = 5 × 105 as in Fig. 3a.

Fig. 5 Top panel: histogram of the residuals between the logarithmic power spectrum of Fig. 3a and its local median, weighted with its local variance σ. Bottom panel: two-dimensional histogram of the same quantity at some frequency at the x-axis and its first harmonic at the y-axis. In addition, we indicate for all counts with χ1 ≥ 5 and χ0 ≥ 10 the corresponding frequency in Hertz and its multiple of the neutron star frequency, ν0 = 0.1323 Hz. To generate the histograms, we used the power spectrum shown in Fig. 3a, i.e. σ = 5 × 105.

3.2 SGR 1900+14

Analogouslyto SGR 1806-20, we also re-analysed the archival RXTE data of the giant flare of SGR 1900+14 that occurred on 1998 August 27. The resolution of the signal field is again 1∕800 s, but here we only used2¹⁸ pixels, as the giant flare did not last as long as that of SGR 1806-20. The operational down times²of RXTE during the observation were taken into account for the inference.

As for SGR1806-20, we were able to recover the frequency corresponding to the rotation as ν₀ = 0.1938 Hz, which is consistent with the reported ν₀ = 0.1923 Hz (Olausen & Kaspi 2014). For SGR 1900+14, the corresponding spectrum with the same smoothness prior σ_sm = 5 × 10⁵ is more noisy because of the various operational down times of the instrument, which are marked grey in Fig. 6. In Fig. 6b, we show a snapshot of the entire reconstructed light curve for a more detailed view. Even though no data were recorded between 88.245 s ≤ t ≤ 95.1375 s, the algorithm was able to infer a light curve by extrapolating from the data-constrained regions. This extrapolation is mainly driven by periodic features appearing in the data. As a further natural consequence, the one-σ confidence interval increases drastically during down times. In total, these multiple operational down times of the instrument lead to fewer observed photons and in consequence to a weaker constrained reconstruction of the light curve as well as its power spectrum. Since a smoothness-enforcing prior σ_sm = 5 × 10⁵ does not sufficiently denoise the light curve, we used a slightly stronger prior, σ_sm = 10⁵ (Fig. 6d), for our further analysis.

With the same detection threshold as before, χ₀ > 11, we find twocandidate frequencies at 7.693 and 11.768 Hz, see Table 3.

In Table 4 we give the maximum χ₀ around previously observed frequencies. As for SGR 1806-20, we cannot confirm any of the previously detected frequencies, the highest significance we see is χ₀ > 5.4 for 28.0 Hz. We further investigated our combined criterium, but now reduced to the giant flare of SGR 1900+14 compared to the flare of SGR 1806-20. Figure 9 shows only one additional frequency at ν = 7.695 Hz, which strengthens our confidence in the previous finding in Table 3 at ν = 7.693 Hz. In Fig. 7 we overplot the reconstructed power spectra around ν = 7.693 Hz (black line) and its first overtone (red line). These two peaks are the largest in the given frequency range.

We plot the time evolution of the different pulses of SGR 1900+14 in Fig. 8. As before, blue indicates the beginning of the flare, green the middle, and red the end. Here, the data are of inferior quality compared to the SGR 1806-20 and the parts of the curves without significant short time-variability are purely reconstructed by our algorithm. For SGR 1900+14 we observed four maxima, which also evolved differently compared to each other. For example, the weakest maximum (time ~ 1 s) declines almost monotonically with time, while the others remain rather constant with far less variability in time.

In orderto facilitate secondary studies based on our reconstructions and spectra, we provide these data online at the CDS.At www.mpa-garching.mpg.de/ift/data/QPO we also provides interactive plots of Figs. 6, 1, and 3 for more detailed views.

Fig. 6 Figure 6a shows the reconstructed light curve ϕ for SGR 1900+14. The grey rectangles indicate the operational down times of the instrument. Figure 6b: snapshot of the entire reconstruction between ≈80 s and ≈91 s. For t ≳ 88 the instrument had an operational down time, leading to zero counts in this time interval. Figure 6c and Figure 6d: reconstructed power spectra as well as their uncertainties, using a smoothness-enforcing prior with σsm = 5 × 105 and σsm = 105, respectively. The plots have the same volumes and sampling rate as in Fig. 1.

Fig. 7 Zoomed-in view of the reconstructed power spectra of the giant flare SGR 1900+14 around ν ~ 7.693 Hz in black and its first overtone around ν ~ 15.386 Hz in red. The plot corresponds to σsm = 1 × 105 as in Fig. 6d.

Fig. 8 Pulse profiles of different rotation periods overplotted for SGR 1900+14. The temporal evolution of the pulse profile is visible from blue to green to red. Smooth curves without significant short time-variability are reconstructed during instrument deadtimes.

Fig. 9 Same as in Fig. 5 for the SGR 1900+14 event, with σ = 105 . We indicate for all counts with χ0 ≥ 5 and χ1 ≥ 5 the corresponding frequency in Hertz and its multiple of the neutron star frequency ν0 = 0.1938 Hz.

4 Discussion

We recovered the rotation period of the two magnetars ν¹⁸⁰⁶ = 0.1323 Hz and ν¹⁹⁰⁰ = 0.1938 Hz, and up to the 31th overtone for SGR 1806-20 (χ₀ > 11) and to the 18th for SGR 1900+14 (χ₀ > 7.5). In addition, we only found potential periodic signals at 9.2 Hz (SGR1806-20) and 7.7 Hz (SGR1900+14), which are significant in χ₀ and in the combination of χ₀ and χ₁ according to our respective criterium. For SGR 1900+14 and with the single criterium on χ₀ , we found another candidate frequency at ν = 11.8 Hz, with a similar χ₀ as the signal at 7.7 Hz, but without detected overtones. Therefore we only consider ν = 7.7 Hz as a potential signal. The most robust feature is at 9.2 Hz for SGR 1806-20, it has the highest χ₀ , and we detect some of its overtones as local maxima in χ₀ at 18.3 and 36.4 Hz, respectively. The two latter frequencies are similar to previously reported frequencies at 17.9 (Watts & Strohmayer 2006) and 36.8 Hz (Hambaryan et al. 2011).

Although the reimplemented Bayesian inference algorithm D³ PO proved (see Appendix A–C) its basic ability to reconstruct QPO in photon bursts in many tests, even in the low S/N regime, we cannot confirm the previously reported frequencies at 17, 21, 26, 29, 59, 92.5, 116, and 150 Hz for SGR 1806-20 and 28, 53, 84, and 155 Hz for SGR 1900+14. In particular, we show in Appendix C that our algorithm is able to recover sufficiently strong quasi-periodic signals around 90 Hz with a width of ~0.6 Hz.

Our analysis and therefore all reconstructed power spectra and reconstructed photon fluxes depend on the particular chosen smoothness-enforcing prior σ_sm. Selecting a smaller σ_sm allows the algorithm to denoise the power spectrum even at low S/N. However, at high frequencies, the small σ_sm leads to a lower sensitivity for spectral lines. Therefore, a trade-off needs to be found for σ_sm between better denoising (small σ_sm) and better sensitivity (large σ_sm). Owing to the sharply deteriorating S/N for frequency ranges around 625 Hz and 1840 Hz, which would require a very small σ_sm, we are currently not able to investigate these frequencies for QPOs.

We assume that our findings hold and determine the effect that the results have on the interpretation of the theoretical model of magnetar oscillations. First of all, the new candidate frequencies at 9.2 and 7.7 Hz for SGR1806-20 and SGR1900+14, respectively, are much lower than the frequencies reported so far. Here, our method has a clear advantage over previous studies, as we can analyse the entire data set at once, meaning that no parts of the data set are left out. In principle, this leads to a more robust reconstruction and also counteracts statistical artefacts that may be introduced by a windowed analysis. Our method is fully noise-aware and has in principle no problems with observational dead times. If these two frequencies are related to neutron star oscillations, the parameters of current models change significantly. When we assume that the oscillations are pure crustal shear modes, the identification of the low frequency with the n = 0, l = 2 mode would favour models with rather low shear speeds and therefore EoS with a fast increase of symmetry energy (Steiner & Watts 2009). However, pure shear models are not very likely because of the strong interaction with the magnetic field in the interior of the star (Levin 2006; 2007; Gabler et al. 2011; 2012). If the magnetic field is not neglected, coupled magneto-elastic oscillations need to be considered. In this case, the much lower fundamental oscillation frequency indicates a magnetic field of the order of B̅ ~ 6× 10¹³–3 × 10¹⁴ G, lower by a factor of ~3 than our estimates in Gabler et al. (2013). With such weak magnetic fields, the oscillations would also remain confined to the core for models with very strong shear modulus, and there would be no chance to observe the oscillations exterior to the star. Therefore, similarly to the case of pure crustalshear oscillations, the low frequencies of 9.2 or 7.7 Hz favour EoSs that give a small shear modulus. However, the estimation of the magnetic field strength has some serious problems : i) The spin-down estimate is only accurate to a factor of a few. ii) The particular magnetic field configuration inside a magnetar is unknown, meaning that even for dipolar-like fields in the exterior, there are different possible realization in the interior. iii) In addition to the degeneracies of the dipolar magnetic field strength and the magnetic field configuration, which lead to comparable frequencies, the compactness of the neutron star and the superfluid properties of the core matter also influence the frequencies significantly (Galber et al., in prep.). In order to advance and to lift some of these degeneracies, we need more observationsof QPOs, and if possible, observations for different sources.

In respect of the potential of improving our method by modelling the spectra as an independent combination of a continuum and lines, we postpone a more thorough discussion of these consequences to future work. We also plan to improve our algorithm by taking the deadtime after each photon detection into account and include further instrument responses. However, we expect these improvements to increase ourdetection limits only by about a few per cent.

We have furthermore shown the capability of our code D³ PO to denoise and recover the shape of the light curve. This may have interesting applications in the field of modelling the pulses of X-ray bursters.

5 Conclusion

We have applied a new Bayesian method, D³ PO to analyse time-series data of the giant flares of SGR 1806-20 and SGR 1900+14. Thereby, we ensured that the photon shot noise, operational down times of the instrument, and the a priori unknown power spectrum of the photon flux were taken into account. In order to denoise and reconstruct the logarithmic photon flux and its power spectrum simultaneously from the data, we had to assume some kind of spectral smoothness on the power spectrum to counteract data variance. Our analysis cannot confirm previous findings of QPOs in the giant flare data, but we found new candidates for periodic signals at 9.2 Hz for SGR 1806-20 and 7.7 Hz for SGR1900+14. If these are real and related to the lowest frequency oscillation of the magnetar, our results favour high-compactness, weaker magnetic fields than were assumed before, and low shear moduli. The necessary application of a smoothness-enforcing prior results in a reduced sensitivity to the spectral line in the power spectrum. Hence, we propose for future work to decompose the power spectrum into two components, one smooth background spectrum, and one consisting of spectral lines. In this way, the previously reportedlines may be confirmed or refuted with higher confidence, or additional and as yet unknown frequencies in QPOs might be discovered.

Acknowledgements

We would like to thank Martin Reinecke forhis computational and numerical support, and Ewald Müller and the anonymous referee for their fruitful discussions and comments. Furthermore, we would like to acknowledge Anna Watts, Phil Uttley, and Daniela Huppenkothen for a controversial debate. Work supported by the ERC Advanced Grant no. 341157-COCO2CASA. This research hasmade use of data and/or software provided by the High Energy Astrophysics Science Archive Research Center (HEASARC), which is a service of the Astrophysics Science Division at NASA/GSFC and the High Energy Astrophysics Division of the Smithsonian Astrophysical Observatory.

Appendix A Influence of the smoothness- enforcing parameter σ

To demonstrate the performance of the Bayesian inference algorithm, we challenged the new implementation of the algorithm with mock data. To do so, we drew Poisson samples from λ = Re^s. For simplicity, we assumed R to be the identity operator. The signal field s was drawn from a Gaussian random field s↩"G" (s, S), with (A.1) with ν₀ = 1 absorbing the physical units. To test whether the algorithm can reconstruct discrete frequencies as well, we further injected four δ-peaks into s_v′′ at randomly chosen positions. The peak heights decrease at higher frequencies to be as realistic as possible. For the test we used a one-dimensional regular grid with with 2¹⁶ pixels, each with a volume of 1∕1000 s.

In Fig. A.1 we tested the principle capabilities of the algorithm under the above-mentioned setting while varying the strength of the smoothness enforcing parameter σ_sm. In the scenario of σ_sm = 10⁴, shown in Figs. A.1 and A.1, the algorithm was able to recover the signal as well as the power spectrum well within the one-σ confidence interval. However, the two highest injected frequencies could not be recovered, as their injected strength and therefore their effect on the data was too weak to be distinguished from the Poisson shot noise.

Fig. A.1 Reconstruction, i.e. map and power spectra from mock data, created according to Appendix A. For clarity, Figure A.1a  only shows a snapshot of the events between 32.8 s ≲ t ≲ 33 s. In addition to the raw photon counts, the black line shows the original signal s , as well as the reconstruction in red, including its uncertainty. Figure A.1b shows the reconstructed power spectra as well as its uncertainty, as in Fig. A.1a, using a smoothness-enforcing prior with σsm = 104.  Fig. A.1c, shows a reconstructed power spectrum with a weaker smoothness-enforcing prior with σsm = 105.

In comparison to Figs. A.1b, A.1c show the reconstructed power spectrum using a significantly weaker smoothness enforcing σ_sm = 10⁵. As a first obvious consequence, the local variance of the spectrum increases by multiple orders of magnitude. However, the same two frequencies as with the stronger σ_sm were clearly recovered. As a further consequence of the smaller σ_sm , the power spectrum does not fully recover the shape of s_v′′ as it picks up more noise features from the recovered map and data. For a detailed discussion and mathematical derivation of this smoothness-enforcing parameter, we refer to Oppermann et al. (2013).

Appendix B Recovery of spectral lines at high frequencies

Now we estimate how strong periodic oscillations must be in order to be detectable in the reconstructed power spectrum. To be as close as possible to a realistic scenario, we manipulated the reconstructed photon flux of event SGR 1806-20. We added a periodic signal with discrete frequencies at 18.0, 26.0, 30.0, and 90.0 Hz to the reconstructed ϕ shown in Fig. 1a. The relative strength of this additional photon flux was varied, between 10² and 10⁵ times stronger than the local power of the injected frequencies. From these manipulated fluxes, we again drew Poisson samples and let the algorithm recover the power spectrum and the flux.

Figure B.1 shows the reconstructed power spectra. In the setting of the top panel, we were only able to recover the injected frequency at 18 Hz, the other three could not be inferred as their strength of 10² compared tothe local power was too weak to be identified as part of the signaland not the Poisson shot noise. However, as the one-σ confidence interval of the reconstructed spectrum at 18 Hz ≲ ν ≲ 40 Hz increases, the algorithm is just at the S/N threshold at which it is able to reconstruct discrete frequencies. If the injected strength of the frequencies, that is, the amplitude of the injected signal, becomes larger, as in Fig B.1, the algorithm can infer all frequencies.

Fig. B.1 Reconstructed power spectrum, according to the test scenario described in Appendix B. The strengths of the injected frequencies are 102 (top) and 105 (bottom) times stronger than the local power at these frequencies. The displayed uncertainties indicate the one-σ confidence interval.

Hence, we may conclude that if the data show periodic signals at high frequencies with sufficient strength, they can be reconstructed by the algorithm. The strength of the periodic signal must therefore be at least 10² times stronger than the local power of the signal in order to give a significant imprint on the reconstructed power spectrum.

Appendix C Recovery of quasi-periodic oscillations

Here, we estimate how strong QPOs have to be for our algorithm to still find a significant signal in the reconstructed power spectrum. Similar to Appendix B, we used the reconstructed photon flux, Fig. 1a of event SGR 1806-20, and added a quasi-periodic signal at around 20 Hz and 90 Hz. We assumed that the power spectrum of a QPO has an approximately Gaussian shape with variance σ ≈0.6. From this manipulated flux, we drew Poisson samples and let the algorithm recover the power spectrum and the flux. Figure C.1 shows the power spectrum of the Poissonian photon counts for a pure QPO-like injected signal, smoothed with a Gaussian convolution kernel whose variance is the same as that of the injected QPO. The reconstructed power spectrum is shown in Fig. C.2. The S/N around 20 Hz is much higher than the S/N around 90 Hz. Therefore, the QPO at 20 Hz has a significantly larger amplitude in the reconstructed power spectrum than that around 90 Hz. The strength of the injected QPO at 90 Hz is just at the threshold to give a significant signal in the reconstructed power spectrum. Thus, we may conclude that at this frequency, a QPO has to have a Leahy power of the order of 3 in order be detectable by D³PO.

Fig. C.1 Smoothed Leahy power of the injected pure QPO-like signal. The two QPOs at 20 Hz and 90 Hz were reconstructed by D3 PO, Fig. C.2.

Fig. C.2 Reconstructed power spectrum according to the QPO injection test described in Appendix C. The displayed uncertainties indicate the one-σ confidence interval.

References

All Tables

All Figures

Fig. 1 Reconstructed light curve of the giant flare of SGR 1806-2 using a smoothness-enforcing prior, with σsm = 5 × 105. The grey narrow rectangle indicates the operational down time of the instrument. For better visibility, we plot the light curve between ≈164 s and ≈174 s in Figure 1b. In addition to the raw photon counts (black dots), the black line indicates the reconstruction of the expected photon counts, i.e. λ as well as its one-σ confidence interval. Owing to the high resolution of the photon flux, only every fourth point of the regular grid is plotted. In both plots each pixel has a duration of 1∕800 s.

In the text

	Fig. 2 Pulse profiles of different rotation periods overplotted as given in Fig. 1a. The temporal evolution of the pulse profile is visible from blue to green to red. The mean pulse is shown as a thick black line. All pulses are shifted such that their log-means are zero.
In the text

Fig. 3 Reconstructed power spectra of the giant flare of SGR 1806-20: For Figure 3a and Figure 3c we used smoothness-enforcing priors with σsm = 5 × 105 and σsm = 105, respectively. The uncertainty intervals are given as grey shaded areas. In Figure 3b we show the reconstructed power spectrum again from the data themselves as in Figure 3a, along with the power spectrum of the logarithmic reconstructed light curve of Fig. 1a. Because of the high resolution of the reconstructed power spectra, only every fourth point of the regular grid is plotted. Each pixel has a volume of 1∕655 Hz.

In the text

	Fig. 4 Zoomed-in view of the reconstructed power spectra of giant flare SGR 1806-20 around ν ~ 9.186 Hz in black and its first overtone around ν ~ 18.265 Hz in red. The spectra correspond to σsm = 5 × 105 as in Fig. 3a.
In the text

Fig. 5 Top panel: histogram of the residuals between the logarithmic power spectrum of Fig. 3a and its local median, weighted with its local variance σ. Bottom panel: two-dimensional histogram of the same quantity at some frequency at the x-axis and its first harmonic at the y-axis. In addition, we indicate for all counts with χ1 ≥ 5 and χ0 ≥ 10 the corresponding frequency in Hertz and its multiple of the neutron star frequency, ν0 = 0.1323 Hz. To generate the histograms, we used the power spectrum shown in Fig. 3a, i.e. σ = 5 × 105.

In the text

Fig. 6 Figure 6a shows the reconstructed light curve ϕ for SGR 1900+14. The grey rectangles indicate the operational down times of the instrument. Figure 6b: snapshot of the entire reconstruction between ≈80 s and ≈91 s. For t ≳ 88 the instrument had an operational down time, leading to zero counts in this time interval. Figure 6c and Figure 6d: reconstructed power spectra as well as their uncertainties, using a smoothness-enforcing prior with σsm = 5 × 105 and σsm = 105, respectively. The plots have the same volumes and sampling rate as in Fig. 1.

In the text

	Fig. 7 Zoomed-in view of the reconstructed power spectra of the giant flare SGR 1900+14 around ν ~ 7.693 Hz in black and its first overtone around ν ~ 15.386 Hz in red. The plot corresponds to σsm = 1 × 105 as in Fig. 6d.
In the text

	Fig. 8 Pulse profiles of different rotation periods overplotted for SGR 1900+14. The temporal evolution of the pulse profile is visible from blue to green to red. Smooth curves without significant short time-variability are reconstructed during instrument deadtimes.
In the text

	Fig. 9 Same as in Fig. 5 for the SGR 1900+14 event, with σ = 105 . We indicate for all counts with χ0 ≥ 5 and χ1 ≥ 5 the corresponding frequency in Hertz and its multiple of the neutron star frequency ν0 = 0.1938 Hz.
In the text

Fig. A.1 Reconstruction, i.e. map and power spectra from mock data, created according to Appendix A. For clarity, Figure A.1a  only shows a snapshot of the events between 32.8 s ≲ t ≲ 33 s. In addition to the raw photon counts, the black line shows the original signal s , as well as the reconstruction in red, including its uncertainty. Figure A.1b shows the reconstructed power spectra as well as its uncertainty, as in Fig. A.1a, using a smoothness-enforcing prior with σsm = 104.  Fig. A.1c, shows a reconstructed power spectrum with a weaker smoothness-enforcing prior with σsm = 105.

In the text

	Fig. B.1 Reconstructed power spectrum, according to the test scenario described in Appendix B. The strengths of the injected frequencies are 102 (top) and 105 (bottom) times stronger than the local power at these frequencies. The displayed uncertainties indicate the one-σ confidence interval.
In the text

	Fig. C.1 Smoothed Leahy power of the injected pure QPO-like signal. The two QPOs at 20 Hz and 90 Hz were reconstructed by D3 PO, Fig. C.2.
In the text

	Fig. C.2 Reconstructed power spectrum according to the QPO injection test described in Appendix C. The displayed uncertainties indicate the one-σ confidence interval.
In the text