© The Authors 2025

1. Introduction

The prompt emission of gamma-ray bursts (GRBs) in X-rays and γ-rays shows a complex time evolution (see e.g. Pe’er 2015, for an overview). Unlike supernovae with their well-ordered rise and decay patterns, GRBs’ prompt light curves are extremely variable on short timescales, with a broad diversity in their temporal structure among the burst population. The origin of this variability is still largely unknown. There have been numerous studies of the minimum variability timescales, t_var, min, by various groups, including rapid rise times. In studies with the Compton Gamma Ray Observatory’s (GGRO) Burst And Transient Source Experiment (BATSE) and with Fermi’s Gamma-Ray Burst Monitor (GBM), t_var, min values around or shorter than a millisecond have been reported in the 100 keV–1 MeV band (Walker et al. 2000; MacLachlan et al. 2013), though Golkhou et al. (2014) suggests that these timescales may be underestimated. At the higher energy, > 100 MeV band accessed by Fermi’s Large Area Telescope (LAT), the variability can still be on timescales of less than a tenth of a second (Aldrich & Nemiroff 2024). Such t_var, min timescales yield estimates of the physical size of the GRB emission regions typically on the scale of light minutes (i.e. ∼1 AU), and have provided clues about the presence of a structured, weakly magnetised jet (Camisasca et al. 2023).

Searches for periodic or quasi-periodic signals in prompt GRB emission are particularly interesting, because the specific timescales implied by a quasi-periodic oscillation (QPO) imposes strong constraints on the possible underlying emission mechanism and can potentially constrain the central engine. Some short GRBs have recently been identified as extragalactic giant flares from magnetars, strongly magnetised neutron stars (NSs) known for their extraordinary bursting behaviour in X-rays and γ-rays (Roberts et al. 2021; Trigg et al. 2024). Galactic giant flares exhibit both periodic oscillations associated with the NS’s rotation period, and QPOs associated with torsional vibrations of the magnetar (Strohmayer & Watts 2005; Israel et al. 2005). Kilohertz QPOs were seen in GRB 200415 (Castro-Tirado et al. 2021), indicating the nature of a magnetar giant flare; these QPOs were subsequently identified as possible overtones of crustal oscillations using numerical simulations (Sotani et al. 2023).

Various QPOs have been claimed and contested in the literature in both long and short GRBs, using a range of different statistical methods and tests, sometimes with contradictory results. Of particular note is the contested detection in GRB 090709A, where a QPO was reported with a period of around 8 s in Swift, Konus, Suzaku, and INTEGRAL data (Markwardt et al. 2009; Golenetskii et al. 2009; Gotz et al. 2009), but where subsequently Cenko et al. (2010) showed that removing the overall trend from the light curve can lead to significant false positive detections, and similarly a reanalysis by de Luca et al. (2010) also excluded a significant oscillation. Many studies have focused on the very large catalogues of GRBs observed with the CGRO/BATSE and Swift/BAT instruments. A systematic search for QPOs in 44 bright short GRBs observed with Fermi/GBM, Swift/BAT, and CGRO/BATSE provided no detections (Dichiara et al. 2013) and, similarly, Liu & Zou (2024) searched for QPOs in 532 short GRBs observed with BATSE using Fourier-based techniques and reported no compelling detections. Guidorzi et al. (2016) analysed 215 bright long GRBs observed with Swift/BAT, and report three potential candidates, but at a low significance given the size of the sample. Conversely, Tarnopolski & Marchenko (2021) report on the detection of 34 QPOs in the prompt emission of GRBs observed with Swift/BAT using an approach based on wavelets, and Chirenti et al. (2023) present two kilohertz QPOs in short bursts, GRB 910711 and GRB 931101B, observed with BATSE. A systematic search for QPOs in precursors of both short and long GRBs yielded no significant detections so far (Xiao et al. 2022), except for a candidate in the Swift/BAT and Fermi/GBM observations of the precursor of the kilonova-associated long-duration GRB 211211A (Xiao et al. 2024). Chirenti et al. (2024) report a 19.5 Hz oscillation in GRB 211211A, thought to be the result of a merger between a NS and a black hole. Overall, the view of QPOs in GRBs is somewhat unclear: different methods come to different conclusions on the same datasets, and some detections have been contested on statistical grounds, leaving behind an uncertainty about the robustness and reliability of the existing detections.

Searches for QPOs in GRBs most often employ one of two approaches (or sometimes both): Fourier analysis and wavelets. Fourier periodograms have a long history across astronomy (e.g. van der Klis 1989) and other areas of science in the detection and characterisation of periodic and quasi-periodic signals, and are thus well tested and well understood. As we detail more below, however, Fourier periodograms struggle with non-stationary signals; that is, time series whose statistical properties change as a function of time. In such cases, wavelets are often considered the prime alternative for their ability to characterise time-variable signals (e.g. Foster 1996). Recently, Gaussian processes (GPs) have been presented as a possible alternative in QPO searches in short transients (Hübner et al. 2022b), and have been explored specifically for GRBs by Song & Mao (2024).

On March 7 2023 at 15:44:06.67 UTC, multiple space-based γ-ray telescopes detected GRB 230307A (see e.g. Fermi GBM Team 2023; Dichiara & Fermi GBM Team 2023b) originating from the direction of the Magellanic Bridge. The event’s T₉₀ duration (Kouveliotou et al. 1993) was 41.5 s (Sun et al. 2023); its gamma-ray fluence (10–1000 keV) reached the level of 3 × 10⁻³ erg cm⁻² (Sun et al. 2023), making it the second-brightest GRB ever observed. Further, a lanthanide-rich kilonova coincident with the position and time of the GRB was identified (Yang et al. 2024; Levan et al. 2024), indicating that this event was due to a merger of two compact objects, rather than the core collapse of an evolved star. Sun et al. (2023) performed a comprehensive broadband spectral analysis and identified emission likely related to a magnetar¹ central engine, suggesting that a QPO search might reveal the magnetar’s spin period.

In this paper, we present a thorough QPO search² of the prompt emission of GRB 230307A observed with both Fermi/GBM and INTEGRAL’s SPectrometer of INTEGRAL AntiCoincidence Shield (SPI-ACS). We applied all three major QPO search methods: Fourier analysis, wavelets, and GPs. We report on a candidate detection with a period of 0.82 s present across instruments and detection methods, and a less significant candidate with a period of ∼0.34 s. Beyond this specific search for QPOs in this specific event, we also critically evaluated the assumptions that each method makes of the data, and whether these assumptions are fulfilled by GRBs more generally and this event in particular. As we show, all of the currently available statistical methods for QPO searches in GRBs and other short transients make strong assumptions that are not supported by the data: as a result, the differing assumptions made by each method can lead to strong disparities in the significance of a candidate signal and no clear, robust conclusion about the presence of that signal.

In Section 2, we present the data used in this study. Section 3 presents our results using Fourier-based methods, Section 4 the analysis of the same data using wavelet transforms, and Section 5 a QPO search using GPs. For all three techniques, we critically evaluate the robustness of the method and the significance of the candidate signals. Finally, Section 6 connects these candidates with potential physical mechanisms, and considers the ability of current methods to robustly detect and characterise quasi-periodic signals in short transients.

2. Data

We analysed data collected with two γ-ray instruments: the Anti-Coincidence Shield (ACS) of the spectrometer (SPI) on board the INTErnational Gamma-Ray Astrophysics Laboratory (INTEGRAL) and the Fermi/Gamma-ray Burst Monitor (GBM). SPI-ACS provides GRB light curves above 80 keV with 50 ms time resolution (von Kienlin et al. 2003)³. GBM detectors provide Time-Tagged Events (TTE) with a minimum readout capability of 2 μs in the 8–700 keV, 200–10 000 keV bands for the GBM/NaI and GBM/BGO detectors, respectively. We analysed data of two Fermi/GBM detectors (NaI 10 and BGO 1). All remaining GBM detectors either had detector zenith to source angles exceeding 60^o or were blocked by other parts of the spacecraft (Dichiara & Fermi GBM Team 2023a), thus not suitable for analysis.

For the GBM detectors, we used TTE data binned to 50 ms time resolution. When binning the TTE data we corrected for the detector deadtime (τ = 2.6 μs; Meegan et al. 2009) using the nonparalysable formula:

$$\begin{array}{c}\hfill n=\frac{m}{1-m\tau },\end{array}$$(1)

where n is the corrected count rate and m is the recorded count rate (Knoll 2010). Deadtime introduces complex features into the periodogram, but is unlikely to strongly affect our timing analysis on very long timescales compared to the deadtime interval. At long timescales, nonparalysable deadtime such as that observed from Fermi/GBM effectively reduces the amplitude of any present signal (see, e.g. Bachetti et al. 2015; Huppenkothen & Bachetti 2022).

INTEGRAL/SPI-ACS saturates above ∼1.8 × 10⁶ counts/s (e.g. Savchenko et al. 2024). We checked the INTEGRAL light curve for bins that reach or exceed that saturation limit, and found only two bins. To check the effect of these bins, we perform the analysis in Section 3 both with the unaltered data and with a light curve where the two bins above the saturation limit are replaced by the average of its neighbouring bins. This allows us to perform a periodogram analysis, which requires evenly sampled data without gaps. Similarly, we performed the wavelet analysis in Section 4 both with the unaltered data, and with a light curve that has the two time bins in question removed. We find no appreciable differences between the results of the analysis, and thus in this paper solely report the results derived using the original light curve.

Due to the high count rates of the event, GBM TTE data suffered data loss between 2.5 s and 7.5 s after the trigger time (bad time interval, BTI; Dichiara & Fermi GBM Team 2023a), where the data packets are lost due to the bandwidth limit between the instrument and the spacecraft. With this caveat in mind, we proceeded by analysing the full light curves for both Fermi/GBM and INTEGRAL/SPI-ACS, as well as light curves where the BTI has been excised (where the methodological approaches described below allow for light curves with gaps). We converted all photon arrival times to the solar system barycenter to align the observations and remove any potential effects due to spacecraft orbital motion.

To determine the background level for each instrument, we fitted a linear function to data segments in the pre-burst (−25 to −5 s) and post-burst (100 to 120 s) intervals using the entire energy range. Then, we subtracted the background level from each corresponding time series to obtain the background-free burst data. While the T₉₀ duration of this burst is 41.5 s (Sun et al. 2023), we performed the search in the 0–60 s interval to capture any variability in the remainder of its gamma-ray emission phase. Figure 1 (left panels) shows the light curve of the kilonova-associated GRB, which exhibits high degree of flux variability, motivating us to search for possible periodic or quasi-periodic modulations in the light curves of this event.

Fig. 1. Top: Light curves of GRB 230307A in 50 ms temporal resolution as seen with SPI-ACS (left panel), and the brightest GBM NaI and BGO detectors; NaI 10 (middle panel) and BGO 1 (right panel). The time is in seconds since the trigger time. The orange area within the lower plots marks the time interval of 2.5 to 7.5 s, for which the GBM team issued a warning for possible data problems. Bottom: Fourier periodograms corresponding to the GRB light curves on the top. We show both the unbinned periodogram (black) and the log-binned periodogram (blue). Note that for the Fermi/GBM data, these do include the segment for which a warning was issued. All three periodograms contain strong variability above the instrumental noise limit at all frequencies considered here, and show peaks on top of the broadband variability present across all frequencies in the periodogram.

3. Fourier-based methods and analysis

The most common standard method of searching for periodic and quasi-periodic signals in astronomical data uses the Fourier transform. An introduction into the formalism and statistical background can be found in van der Klis (1989) and Bachetti & Huppenkothen (2023). We produced periodograms for the full GRB light curve for each instrument, normalised using the formalism by Leahy et al. (1983), and show the periodograms in Figure 1 (bottom row of panels). The periodograms show a high amount of excess power above the Poisson noise level across the entire frequency range we probed, and in particular a set of broad, prominent features in all three instruments peaking near 0.45, 0.6, and 1.2 Hz.

3.1. Search for (quasi-)periodic oscillations

We followed the approach introduced by Vaughan (2010), adapted to transients by Huppenkothen et al. (2013), and implemented in the stingray Python library (Huppenkothen et al. 2019; Bachetti et al. 2024) to search for periodicities and narrow QPOs. First, we fitted a model for the broadband timing variability at frequencies below 10 Hz. We chose between a simple power law and a broken power law by applying a likelihood ratio test (LRT). The LRT was calibrated using simulations from the simpler model, drawn from the posterior for the parameters of the power law model. We chose uniform priors between 0 and 5 for the power law indices, and wide, log-uniform priors for the normalisation of the model, and for the break frequency. We sampled the power law model using Markov chain Monte Carlo, implemented in the package emcee (Foreman-Mackey et al. 2013), and checked visually and through computation of the autocorrelation time whether the chains had converged. We then simulated 1000 periodograms from the posterior, and computed the LRT for each, in order to compute a tail probability (p value) for the observed LRT. For all three datasets, p > 0.05, suggesting that the null hypothesis (a power-law model) cannot be rejected.

We subsequently simulated another 1000 periodograms from the posterior for this model, performed a maximum likelihood fit of the power-law model for each, and computed the highest outlier in each periodogram. The null hypothesis that the maximum power in the observed periodogram can be explained by intrinsic, non-periodic variability can be rejected if the tail probability for the observation based on the sample of simulated highest outliers is small. For all three instruments, this null hypothesis cannot be rejected (p > 0.05).

Outlier detection methods are most powerful when the putative signal is concentrated in a single frequency bin, as is the case for a very narrow QPO or a strictly periodic signal. As is shown in Figure 1, the peaks visible below 1.5 Hz are distributed across multiple frequencies. To search for these signals, we performed another LRT, but comparing the power-law model for the broadband noise with one that also includes a Lorentzian component to model a putative QPO. As for the power law, we used a wide, log-uniform prior for the Lorentzian’s normalisation and a log-uniform prior across the entire frequency range of the spectrum for the Lorentzian centroid frequency. We parametrised the Lorentzian’s width in terms of the quality factor, q = ν_c/Δν, where the quality factor is defined as the ratio between the centroid frequency and the full width at half maximum (FWHM), and applied a uniform prior between 2 and 100. We compared that LRT with those calculated for periodograms simulated from the null hypothesis (only the power law).

For all three instruments, we confidently reject the null hypothesis with p < 0.001, with a putative QPO detection at 1.2 Hz (see Figure 2 for details for the INTEGRAL data; additional figures for the NaI and BGO data are in Appendix A, Figures A.1 and A.2). This is somewhat unsurprising: the null hypothesis involves a stochastic process where powers in neighbouring bins are statistically independent. It is obvious from the periodograms that this is not the case here: for the peaked structures in the periodogram, where neighbouring bins are clearly correlated in some way, it would make sense that the LRT highly favours a model with a Lorentzian. Adding a second Lorentzian component to the model and comparing that model with the power law and single Lorentzian also yields a significant rejection of the latter model in the INTEGRAL data (p = 0.006), but not in the Fermi/GBM NaI (p = 0.033) or Fermi/GBM BGO detectors (p = 0.199). We found including additional Lorentzians beyond two challenging in practice, even when exploring a wide range of starting parameters. Often, these components would optimise to local minima and very broad, flat features as part of the broadband variability.

Fig. 2. Left: Fourier periodogram of the INTEGRAL data with posterior draws from the three models compared via LRTs: in green, the power law model; in blue, a power law model with a Lorentzian component for a single QPO; in orange, a model comprising a power law and two Lorentzians. Middle: Distribution of the likelihood ratios from 1000 simulated periodograms. The likelihood ratio for the observed periodogram is a clear outlier. Right: Same as middle panel, but for the model with two QPOs. Again, the observed likelihood ratio is a clear outlier compared with the null hypothesis (a single QPO).

Overall, Fourier analysis suggests that there exist at least one, and possibly two QPOs in the data. We sampled the posterior for both QPO components along with the power law model for all three instruments in order to obtain credible intervals on the QPO properties. For INTEGRAL, we find a centroid frequency of the first QPO components of ${\nu }_1=1.{217}_{-0.027}^{+0.0281}$ Hz and a FWHM of $\mathrm{\Delta }{\nu }_1=0.{094}_{-0.059}^{+0.093}$ Hz, corresponding to a quality factor of ν₁/Δν₁ = 12.9. Parameter estimates for the BGO data are consistent, while the results for the NaI data suggest a larger FWHM, $\mathrm{\Delta }{\nu }_1=0.{32}_{0.16}^{+0.21}$ Hz (though we note the substantially larger errors on the NaI result). For the higher-frequency QPO in the INTEGRAL data, we find a centroid frequency of ${\nu }_2=2.{987}_{-0.035}^{+0.083}$ Hz and fairly large credible intervals for the QPO width, $\mathrm{\Delta }{\nu }_2=0.{120}_{-0.095}^{+0.133}$ Hz.

To check for any energy dependence of the QPOs, we generated Fermi light curves in different energy ranges: 8 − 40 keV, 40 − 200 keV, and 200 − 700 keV for the NaI detectors, and 200 − 700 keV, 700 − 3000 keV and 3000 − 10 000 keV for the BGO detectors. The shared energy band between the two GBM detectors serves as a useful cross-check for consistency. Since the SPI-ACS detector does not record the energy of incoming photons, we cannot perform a similar analysis on the Integral data (von Kienlin et al. 2003). For the Fermi/GBM light curves, we find that the fractional root mean square (rms) amplitude of the QPO at 1.2 Hz is strongly energy-dependent (Figure 3, increasing from ∼16% to ∼88% from the lowest to the highest-energy band). This behaviour persists across the detectors, and the fractional rms amplitude for the overlapping band is consistent within statistical uncertainties. Because of the lower significance, and the wide credible intervals on the QPO width for the second QPO in the Fermi/GBM data, we did not perform a similar analysis for the potential QPO candidate at 2.9 Hz.

Fig. 3. Fractional rms amplitude as a function of photon energy for the two Fermi/GBM detectors.

The increase in the rms amplitude with energy band might well be expected for a wide array of physical emission mechanisms. The time variability of the flux at different energies will naturally couple to spectral variability: the spectral shape will likely not be preserved as the flux goes up and down, and at energies where the spectrum is steeper, fluxes are likely to vary more. This effect is most pronounced if the break energy varies. When this happens, the fractional changes in fluxes are likely to be greater at a fixed photon energy above the break than at a fixed energy below the break where the spectrum is flatter. For this burst, Band model fits (see e.g. Table 2 of Dichiara et al. 2023) indicate that the break energy E_br ∼ 635 − 970 keV and the spectral indices α (below E_br) and β (above) all vary substantially over the time interval pertinent to the QPO analysis. Thus, fluctuations in E_br are likely to lead to a higher RMS amplitude for QPOs above this energy than below it.

3.2. Assumptions and limitations

The Fourier-based methods used in this section are well tested and well understood in the context of periodicity and QPO searches in astronomical light curves. As has been described in van der Klis (1989) and Bachetti & Huppenkothen (2023), we understand the statistical properties of periodograms of stochastic processes observed with X-ray and gamma-ray instruments very well. However, these methods make specific assumptions that are a challenge in QPO detections in fast transients such as GRBs. In particular, one assumption that underpins the analysis above is that of weak stationarity. Weak stationarity (or wide-sense stationarity) assumes that the mean of the process and the autocovariance of the process do not change as a function of time. This is generally true in the context of X-ray binaries and active galactic nuclei over the timescales of interest for QPO searches, where many of these methods are routinely applied successfully. However, it is obviously not true for short transients: by definition, they have a beginning and an end, and as in the case of this GRB, may show very significant changes in the variability as a function of time. As is shown in Hübner et al. (2022a), a key consequence of this non-stationarity is that powers in neighbouring bins are no longer statistically independent. Any method that relies on that independence, as the ones above do, will tend to overestimate the significance of any candidate signal.

The mismatch between the model above and the data considered here is easily illustrated by generating periodograms from the posterior for the parameters of the model with two QPOs, and simulating light curves from those periodograms using the method in Timmer & König (1995). We show the results in Figure 4. The mismatch between the GRB and the process assumed for both the overall variability and the possible QPO is immediately obvious: while a stationary stochastic process can still display large-amplitude variability, the properties of this variability do not change as a function of time, whereas for the GRB they do.

Fig. 4. Left: INTEGRAL light curve of GRB230307A (black), with three random light curves generated from the stochastic process used in the QPO detection methods outlined in this section (orange). The GRB has a well-defined beginning and an end, in between which there exists rapidly changing variability. The simulated light curves also contain variability at a high amplitude, but the overall process does not change throughout the light curve. This is expected for a stationary stochastic process. Right: Periodogram of the INTEGRAL data (black) and of the simulated light curves (orange). While the periodogram of the GRB exhibits peaks formed by excess power in correlated neighbouring frequency bins, the periodograms of the stochastic process contain – by design and construction – powers that are statistically independent.

One possible solution to the above problem could be to de-trend the light curve. This would entail fitting a model to the overall shape of the burst, subtracting the best-fit model, and then continuing with the Fourier analysis assuming the residuals follow a stationary stochastic process as assumed by the statistical methods above. However, this approach also comes with a range of challenges. First and foremost, it assumes that the observed light curve can be neatly separated into an overall burst shape and a stationary stochastic process. This may be the case, but the presence of variability from the stochastic process will significantly bias the fit to the overall burst, and thus also introduce biases into the residuals that again lead to a departure from the statistical assumptions of the method. Additionally, removing an overall trend generally removes power from the smallest frequencies in the periodogram. As Cenko et al. (2010) showed, doing so may lead to artefacts in the periodogram of the residuals that can mimic a QPO, and thus lead to spurious detections, a result recently confirmed by Song & Mao (2024). Finally, in this particular case, it is unlikely that detrending will solve the underlying problem. In all three instruments, the amplitude of the variability appears to be time-dependent; that is, the amplitude of the variability in the first ten seconds or so is significantly larger than twenty seconds after the burst. Thus, removing an overall trend will not yield residuals that can be modelled with a stationary stochastic process. We thus did not proceed with detrending.

4. Wavelets

A common method applied to short bursts and other astronomical transients is the wavelet transform (e.g. Hurley et al. 1998; Morris et al. 2010; Golkhou et al. 2014), a generalisation of the Fourier transform to a wider range of basis functions beyond sines and cosines. Crucially, most common wavelet functions are localised in both time and frequency. As a result, they do not require stationarity and are often used, for example, to detect time-dependent periodic and quasi-periodic signals.

4.1. Search for (quasi-)periodicities

We performed a wavelet transform with the Python library PyWavelets (Lee et al. 2019) using a complex Morlet wavelet, defined as

$$\begin{array}{c}\hfill \psi (t)=\frac{1}{\sqrt{\pi \beta }}exp(-\frac{t^2}{\beta })exp(2\pi i{\nu }_{c}t),\end{array}$$(2)

where β is the bandwidth and ν_c is the centre frequency. The bandwidth, β, describes the spread of power in the frequency domain (corresponding to a time decay in the time domain). The centre frequency, ν_c, should be chosen near the frequencies of interest to be explored using the wavelet transform. The complex Morlet equation above describes a complex exponential windowed by a zero-centred Gaussian, with a width set by $\sqrt{\beta /2}$. We chose a centre frequency of ν_c = 1.5 Hz, broadly in the range where we expect to see signals based on the periodogram, and a bandwidth of β = 10, corresponding to around 20 rotational cycles. We show the two-dimensional wavelet transform and the one-dimensional wavelet periodogram for the INTEGRAL data in Figure 5.

Fig. 5. Left: 2D wavelet transform (spectrogram) of the INTEGRAL observation of GRB230307A. The transform shows transient power at low frequencies in the first ∼20 seconds or so of the burst, where the variability is particularly strong. It also shows a short, transient signal between 0 − 10 s at 1.2 Hz, similarly to what was identified in the Fourier periodogram. The candidate signal at 0.32 Hz identified in the Fourier periodogram is less apparent here. Right: Fourier periodogram (black) and wavelet periodogram (blue) with the candidate signals found in the Fourier analysis noted as dashed orange lines. The Fourier and wavelet periodograms largely match.

The two-dimensional wavelet transform (or spectrogram) shows strong power at the lowest frequencies, and also clearly shows the candidate signal at 1.2 Hz. The candidate signal at 2.9 Hz, however, is less visible, though there are some faint structures at higher frequencies. All power is concentrated early in the light curve, suggesting that any signal is transient and only present in the first ∼10–20 seconds. The wavelet periodogram (Figure 5, right panel) is very consistent with the Fourier periodogram, as expected for this kind of wavelet.

Estimating significance in wavelet spectrograms is generally a challenging task. The standard approach in the literature (e.g. Ghosh et al. 2023) takes a variability model (often a power-law-like red noise model or equivalent) and generates simulated light curves from this model, and significance is then determined using per-bin outlier tests comparing the observed wavelet power in a bin to the distribution of simulated wavelet power in the same bin. We followed this approach and simulated 20 000 light curves from the posterior for the power law model generated in Section 3 using the approach in Timmer & König (1995), and subsequently produced wavelet spectrograms for each simulated light curve in the same way as for the real data. We then computed 99.99% percentiles for each time-frequency bin, and identified those bins where the observed wavelet power exceeds the 99.99% percentile derived from the simulations. We performed this analysis independently for each dataset (i.e. INTEGRAL, Fermi/GBM NaI and BGO detectors), using their individual posterior distributions for the power law model. The results are presented in Figure 6. This approach clearly identifies regions in the spectrogram where the observed wavelet power is an outlier compared to the simulations with p < 0.0001. In all three instruments, the candidate signal at 1.2 Hz is clearly detected, and present only in the first ∼10 s or so of the burst. In Fermi/GBM, this overlaps with much of the interval flagged by the Fermi/GBM team for potential data issues (BTI). However, we also note the similarity of the period, width and temporal extent of the signal in all three instruments, which may suggest that these data issues do not significantly affect the detection of the QPO. We find a weaker area of significance at the candidate signal around 2.9 Hz identified in the Fourier analysis. This, too, is highly localised to a small segment of the overall burst within the first ∼5 s after the trigger, though we also note that the signal evolves temporally to a higher frequency with time in all three instruments. There are additional small regions where bins exceed the significance threshold, but we find that these are not consistent in time or frequency across the three instrument, and thus consider them likely spurious candidates.

Fig. 6. Wavelet spectrograms for GRB 230307A for all three instruments. In all three, we find significant power at the lowest frequencies, as well as power at 1.2 Hz. In white contours, we overplot the 99.99% percentiles. In the Fermi/GBM data, we mark the segment flagged by the Fermi team for potential data issues as the shaded region.

The results presented in Figure 6 constitute single-trial p values, and are thus not corrected for the number of trials. Estimating the number of trials for a wavelet spectrogram is challenging because for all but very specific choices of frequency and time resolution, neighbouring bins in the spectrogram are not statistically independent. In addition, a trial-corrected p value calculated using the approach above would require millions of simulations given that we would need to correct for 76 680 frequency bins (making the conservative but incorrect assumption of statistical independence), which is computationally prohibitively expensive. The fact that there are patches of significance rather than individual pixels does not necessarily constitute additional evidence for the existence of a QPO: because neighbouring bins are not statistically independent, observing patches of high significance can be a result of that lack of independence, rather than evidence for the existence of a real signal.

We also performed the most conservative possible estimation of significance, and compared the single bin with the highest wavelet power to the distribution of the highest wavelet power found across all simulations, and found that, using this metric, none of the candidate QPOs are significant (p = 0.1 for both INTEGRAL and Fermi/GBM NaI data, respectively; p = 0.075 for the Fermi/GBM BGO data). However, we note that this is likely too conservative, because it does assume statistical independence of spectrogram bins.

Finally, we used the simulations to compute single-trial significances in the wavelet periodogram integrated over time (Figure 7). This analysis is similar to the Fourier periodogram in that integrating over the time dimension means we lose the advantage of time-dependent modelling of the wavelet spectrogram. Calculating single-trial significances at a significance threshold of p < 10⁻⁴ corresponds to a trial-corrected significance threshold of p < 0.012, or approximately 3σ. We find that in agreement with the results from the Fourier analysis, the candidate QPO at 1.2 Hz is above this 3σ threshold in all three instruments, whereas the 2.9 Hz signal only reaches 3σ in the Fermi/GBM BGO data.

Fig. 7. Wavelet periodograms (black) and red noise simulations (orange) for all three instruments (top: INTEGRAL; middle: Fermi/GBM NaI; bottom: Fermi/GBM BGO). The wavelet periodogram corresponds to the 2D wavelet spectrogram integrated over the time axis. In orange, the posterior mean derived from 1000 simulated wavelet periodograms, along with 10 posterior draws from the power-law stochastic model sampled in Section 3. In blue, we show the 99.99% single-trial detection limit. The candidate QPO at 1.2 Hz exceeds that limit for all three instruments, while the candidate QPO at 2.9 Hz exceeds this limit only for the Fermi/GBM BGO data. Note, however, that integrating over the time axis will necessarily yield lower significances given the short-lived nature of both candidate signals in the wavelet spectrogram.

4.2. Assumptions and limitations

A key limitation of the wavelet-based method is the challenging determination of significance and the number of trial correction in the presence of correlated frequency- and time-bins. In the analysis above, we have stated both single-trial significances as well as results for the most conservative possible assumption (that all bins are statistically independent), which will vastly decrease sensitivity of the QPO detection. The truth likely lies somewhere in the middle. Given the persistence of the signal at 1.2 Hz in the wavelet periodogram and its consistency with the Fourier analysis, we consider this to be a strong candidate for a QPO. The putative signal at 2.9 Hz is a little less clear: it, too, is present in all three wavelet spectrograms, but with less consistency than the 1.2 Hz signal. Additional power in the wavelet periodogram is inconsistent between the three instruments, and we consider these likely to be statistical artefacts. We note that a fraction of the 1.2 Hz signal falls into the segment flagged as subject to data quality issues in Fermi/GBM. However, given the consistency of the signal in Fermi/GBM and INTEGRAL, and that a part of the signal is present before the BTI, we suggest that the relevant data issues have not significantly impacted the detection of the QPO.

There is, however, a somewhat more fundamental open challenge with the analysis above. While wavelet transforms are better suited to the detection of transient signals, and especially of transient periodic signals, we still must compare them to a model parametrising the null hypothesis. Here, we followed the standard analysis for QPOs with wavelets in the literature and used the same stochastic power-law-shaped process as defined in Section 3. This means that irrespective of the transform applied to the data (and simulations), the same caveats nevertheless apply: a stationary stochastic process is not a good representation of the underlying data (see also Figure 4), which in turn will affect the trustworthiness of the significances derived using simulations generated from this process.

This is also the reason why we do not combine p values for the three instruments for any of the analysis methods considered in this paper: doing so is only permissible if the assumption of independence holds. Here, this assumption holds only for instrumental noise, which will be generated independently in each detector, but not for stochastic variability in the GRB itself, which is produced at the source. To derive more reliable significances requires a more realistic, non-stationary model for the data, which we define in the next section.

5. Gaussian processes

We searched for QPOs in the time domain using the method developed in Hübner et al. (2022b) based on GPs (e.g. Williams & Rasmussen 2006; for an introduction to GPs in astronomy, see Aigrain & Foreman-Mackey 2023). In short, the method directly models the light curve as a combination of a non-stationary trend function parametrising the overall shape of the burst and combines this with a stochastic process to model variability on top of this trend function. Through Bayesian model comparison, different classes of models (e.g. different types of trend functions) can be compared. This model can take into account the non-stationary nature of a burst in a more principled way through the trend function, though we note that we can currently not yet implement a nonstationary QPO except in some simple forms described in Hübner et al. (2022b). Due to the exceptionally bright nature of this GRB, the Gaussian measurement uncertainties for the data assumed in GPs are broadly applicable.

The approach chosen here appears similar to detrending, but has the advantage that it can take into account uncertainties in the parameters of the detrending function. It can also correctly account for correlations between the parameters of that function and the variability not modelled by the trend function. By simultaneously considering both the trend function and the variability on top of it, we can derive appropriately unbiased estimates of the trend function, and take into account uncertainties in that estimate, as well as correlations between the parameters of the trend function and the stochastic process.

We chose a skew-Gaussian function⁴ as a trend function, defined as

$$\begin{array}{c}\hfill f(t)=A\{\begin{array}{cc}exp\left(\frac{-{(t-t_c)}^2}{2{\sigma }_1^2}\right)\hfill & t<t_c\hfill \\ exp\left(\frac{-{(t-t_c)}^2}{2{\sigma }_2^2}\right)\hfill & t\ge t_c.\hfill \end{array}\end{array}$$(3)

The skew-Gaussian function provides a flexible model for the asymmetric, approximately exponential rise and decay of the burst. We compared two stochastic processes: one model parametrises a damped random walk (DRW), the second that same DRW combined with a QPO parametrised as a stochastically driven damped harmonic oscillator. The DRW, also known as an Ornstein-Uhlenbeck process or an autoregressive (AR) process of order 1 (AR(1)) is a fairly simple stochastic process that parametrises the flux at time t + 1 in terms of the flux at the previous time, t, and a random component. In Section 5.4, we consider higher-order AR moving-average processes as an alternative to this process.

We implemented wide, uninformative priors reflecting our lack of prior knowledge in most of the model parameters (see Table 1). Many relevant parameters depend on the properties of the data: amplitudes depend on the sensitivity of the instrument, and the frequency range we can search depends on the length of the GRB. For those parameters, we set priors that reflect the ranges we can reasonably expect to see in our data. Priors on parameters for the trend function and the DRW that exist in both the model with and without QPO are the same for both.

We constructed each model using the GP library tinygp (Foreman-Mackey et al. 2023) and sampled the posterior using nested sampling, implemented in the Python package jaxns (Albert 2020), with 2000 live points. We compared the two models using the Bayes factor, ℬ₂₁, for the model with a QPO (model ℳ₂) versus a model without (model ℳ₁). In logarithmic form, a positive value for log₁₀(ℬ₂₁) can be interpreted as evidence for the presence of a QPO component, whereas a negative Bayes factor can be interpreted as evidence against. The significance of the Bayes factor was calibrated using the common scale by Kass & Raftery (1995), which considers log(ℬ₂₁) > 2 as decisive evidence for model ℳ₂ (more information on the method, including simulations to calibrate its ability to detect QPOs, can be found in Hübner et al. (2022b)).

5.1. SPI-ACS GP QPO search

Comparing GP models with and without QPO yields a highly significant signal centred at the frequency of 1.21 ± 0.01 Hz with a Bayes factor of log₁₀(ℬ₂₁) = 3.70, indicating decisive evidence for the presence of a QPO component in the data (Kass & Raftery 1995). The posterior distribution for the QPO’s centroid frequency is narrow and unimodal, and the signal is fairly coherent, with a quality factor $log(q)=2.{63}_{-0.82}^{+1.18}$, corresponding to q ≃ 13, highly consistent with the estimate from the Fourier analysis. We present the corresponding posterior period distribution in the right panel of Figure 8, which peaks at 0.82 ± 0.01 s. Overall, the posterior distributions are well-constrained for both ℳ₁ and ℳ₂ (see the corner plot in Figure A.3 in Appendix A, generated with the Python package corner (Foreman-Mackey 2016)), and uncertainties in the individual evidences used to compute the Bayes factor are small (∼0.1). The sampling results are stable across multiple runs and the Bayes factor remains significant even when changing the uninformative priors significantly.

Fig. 8. Left: SPI-ACS light curve (black points) in units of counts per 0.05 s bin, the predicted rates from the maximum a posteriori combined model consisting of the mean function, the stochastic process and a QPO (blue), and posterior draws from the mean function (orange). Right: Posterior probability density for the QPO centroid period.

Even though the SPI-ACS data are not affected by the same instrumental issues during the brightest intervals of the GRB as Fermi/GBM is, we repeat the analysis with the Fermi/GBM BTI excised, in order to directly compare results across instruments, under otherwise identical conditions (i.e. with the same priors and nested sampling settings). We find evidence of a strong QPO in this light curve as well (Figure 9), but a different frequency broadly consistent with the second candidate signal found in the Fourier and wavelet representations: ν₀ = 2.880 ± 0.001 Hz, corresponding to a period of 0.347 s. The signal’s significance is ℬ₂₁ = 3.99 ± 0.1, indicating a decisive preference for the model with a QPO. As with the full dataset, this result is stable across nested sampling runs and when varying the already very wide priors. The signal has a high quality factor, $log(q)=4.{77}_{-0.52}^{+0.57}$, indicating that the signal is quite coherent. We note that while this period is much shorter, it does not correspond to a harmonic of the signal detected in the full data.

Fig. 9. Re-analysis of the SPI-ACS data excising the brightest part of the GRB (i.e. BTI), which is affected by data loss in Fermi/GBM. In the left and middle panel, we show observations in black, the maximum a posteriori combined model in blue, and draws from the posterior probability density for the mean function in orange. The left-most panel shows the results for ℳ1 (skew-Gaussian mean function and DRW), the middle panel for ℳ2 (skew-Gaussian mean function, DRW and QPO). In the right-hand panel, we show the posterior probability density for the period of the QPO, which is narrowly constrained around P = 0.347 s.

The wavelet representation suggests that much of the 1.2 Hz signal is covered by the BTI, and thus excised from the light curve along with that BTI. On the other hand, the candidate QPO at 2.9 Hz in the Fourier and wavelet representations is concentrated into a short interval before the BTI. When considering a model containing only one QPO, it makes sense that the posterior will be concentrated on the stronger of the two signals. When considering the whole light curve, this will be the 1.2 Hz QPO. However, removing the BTI, and much of the 1.2 Hz QPO with it, yields a light curve requiring a model with a QPO at 2.88 Hz.

5.2. Fermi/GBM GP QPO search

We performed a similar analysis with the data of both NaI (na) and BGO (b1) detectors. For each, we modelled the full light curve – including BTIs – with both ℳ₁ and ℳ₂, and found results that are highly consistent with the SPI-ACS data (see Figure 10). Both light curves show highly significant QPOs, with Bayes factors of log₁₀(ℬ₂₁) = 2.12 (NaI detector) and log₁₀(ℬ₂₁) = 4.68 (BGO detector). In both, the posterior probability density for the period is very constrained, $P_{\mathrm{NaI}}=0.{825}_{-0.02}^{+0.07}\mathrm{s}$ and P_BGO = 0.826 ± 0.009 s. The distribution for quality factors for the NaI data is somewhat broader, $log(q)=2.{37}_{-3.47}^{+1.52}$, reflecting a broader period posterior, but again we see very high coherence of the signal in the BGO detector, $log(q)=3.{34}_{-1.16}^{+1.51}$.

Fig. 10. Left panels: Light curve (black), maximum a posteriori model (blue) and posterior draws from the mean model (orange) for the Fermi/GBM NaI detector (na; top) and the BGO detector (b1; bottom). Right: Corresponding posterior probability densities for the period in the model including a QPO.

The strong consistency between the results for SPI-ACS and Fermi/GBM leads us to conclude that the data loss likely did not significantly impact the timing results. We note that the period posterior for the NaI detector shows a minor mode at ∼0.34, where we found a significant signal when excluding the Fermi/GBM BTI in the SPI-ACS data. We also analysed the Fermi/GBM data with the BTI excluded in order to obtain a light curve unaffected by data loss and less affected by deadtime. As with the SPI-ACS data, we find that the QPO period has shifted to P = 0.3476 ± 0.008 s, consistent with the SPI-ACS results (Figure A.4 in Appendix A). For both datasets, we find high Bayes factors, log(ℬ_21, NaI) = 5.25 and log(ℬ_21, BGO) = 5.57, respectively, and strongly constrained, high quality factors for the QPO, $log(q_{\mathrm{NaI}})=6.{75}_{-2.73}^{+1.35}$ and $log(q_{\mathrm{BGO}})=4.{69}_{-0.78}^{+0.72}$.

5.3. Detection sensitivity

Bayes factors are notoriously sensitive to prior choices. The empirical model we employed means that we largely chose wide, uninformative priors, which in turn may affect the sensitivity to QPOs in the data. To calibrate the detection sensitivity for the Bayes factor, we simulated fake GRBs drawing only from the joint posterior for ℳ₁; that is, the model without a QPO. For each draw, we generated a trend function and sampled from a GP combining that trend function with a realisation of the DRW process. We added Poisson photon counting noise to this simulated GRB to generate a realistic light curve. We modelled this simulated light curve in the same way we do for the real SPI-ACS data, using the same models and prior assumptions. We repeated this procedure for a 100 different simulations⁵ drawn from the posterior, and generated a distribution of Bayes factors expected under the model without a QPO. This allowed us to explore what range of Bayes factors we would expect under model ℳ₁.

In Figure 11, we present the resulting distribution; compared to the observation, we find much smaller Bayes factors for the simulated observations, indicating that we should not have observed the high Bayes factor recorded for the SPI-ACS data if the latter had been generated by a simple stochastic process and the trend function.

Fig. 11. Distribution of the Bayes factors obtained by modelling 100 simulated GRBs based on a model consisting of a skew-Gaussian trend function and a DRW, but without a QPO (grey). The Bayes factor obtained for the SPI-ACS data is a strong outlier compared to the simulations, indicating that our analysis captures additional variability not present in the simpler model.

5.4. Alternative stochastic models

The DRW model considered above to explain the variability in the GRB is a fairly simple model with a power spectrum constrained to Lorentzian centred at zero. While our goodness-of-fit test has shown that this model can explain the data, we nevertheless implement a somewhat more complex model: a continuous autoregressive, moving-average (CARMA) model of order p = 2 and q = 1. CARMA models consider both an AR process and a moving-average (MA) process simultaneously. Here, the AR describes the future of a system based on its current state and a random perturbation. The MA process parametrises the time series of a system in terms of a signal and its convolution with an impulse response function. The orders, p for the AR process and q for the MA process, describe the time lags and length of the impulse response, respectively (Moreno et al. 2019). The CARMA(2,1) process we implement here is capable of a wider range of power spectral shapes. We limit ourselves to this process, since higher-order CARMA processes can intrisically include QPO-like behaviour, making them not practicable as an alternative model in the context of QPO searches. We use the Python package tinygp for the implementation of an CARMA(2,1) process, and combine it with the same prescription for the mean function above. The priors for the parameters α and β of the CARMA process are similarly wide and uninformative as for the DRW (see Table 1); the amplitude of the process is in this parametrisation folded into the MA parameter, β.

We modelled the full light curves for all three instruments with both a mean function and a CARMA(2,1) process, and a model that additionally includes a QPO. For all three light curves, we once again find Bayes factors that strongly favour the model including a QPO (log₁₀(ℬ₂₁) > 2 for all three datasets). The QPO period is similarly constrained to 0.83 s in with this model as with the model containing a DRW (see Figure 12 for an example). The independence of our results of the chosen model for the underlying variability strongly suggest that the QPO is real.

Fig. 12. Left: INTEGRAL SPI-ACS observation (black) with the maximum a posteriori model comprising a CARMA(2,1) process, a QPO and a skew-Gaussian mean function (blue), with posterior draws from the mean function in orange. Right: Posterior probability density of the period parameter shows a clear peak at 0.82 s, consistent with results obtained using the simpler DRW model.

When we exclude the flagged segment (BTI) of the burst, we find strong signals at 2.9 Hz (0.34 s) in the data of all three instruments with similarly significant Bayes factors (see Figure 13) and results consistent with the analysis using the DRW model.

Fig. 13. Period posteriors for all three instruments for the full light curves (left) and for the light curves without the Fermi/GBM BTIs of 2.5–7.5 s after trigger (right). The posteriors show a remarkable degree of agreement between the three instruments, at frequencies corresponding to ∼1.2 Hz (0.83 s period) for the full light curves and at ∼2.9 Hz (0.34 s) without BTI.

5.5. Assumptions and limitations

Gaussian processes provide two key advantages over standard Fourier analysis and wavelet methods: they enable robust modelling of unevenly sampled light curves directly in the time domain, minimising aliasing and windowing effects, and they enable to joint modelling of a stochastic process together with an overall, deterministic trend, as we do here. This comes at considerable computational cost to calculate and calibrate the Bayes factors.

A main assumption of GP modelling is that data uncertainties are normally distributed. While not strictly true for the data considered here, GRB 230307A is bright enough for the Gaussian approximation to be justified. A second key assumption is that the GRB can be decomposed into a linear combination of a trend function parametrising the global rise and fall of the burst, and a stochastic process parametrising the variability on shorter timescales. Empirically, this appears to be not a bad assumption, but we also note that none of the models implemented here – GPs or in any of the previous methods – are physically motivated. While the DRW is a simple choice for a covariance function, increasing the flexibility by considering a higher-order CARMA process did not substantially alter our conclusions. We note, however, that visually the amplitude of the variability over the course of the GRB appears to change as a function of time. This is at odds with the models considered here, and implementation of a model including some form of non-linear variability component is beyond the scope of this paper.

6. Discussion

In this paper, we have presented a thorough analysis of the light curves of GRB 230307A taken with both INTEGRAL/SPI-ACS and Fermi/GBM’s NaI and BGO detectors. Across multiple detection methods, instruments and wavelengths, we consistently identify a short-lived QPO at 1.2 Hz (0.82s) in the first ten seconds of the GRB, when emission is at its brightest. Given that the signal is very strongly present in the SPI-ACS data, we exclude the possibility that its presence in the GBM data could be related to data issues concurrent with the signal’s presence.

A second, short-lived QPO appears to exist simultaneously in all three instruments at 2.8 Hz (0.34s), though not all tests return a confident detection. In the wavelet spectrogram, this signal appears to show an upward frequency trend over its lifetime. This signal is especially strongly present in QPO detection tests that excise the Fermi/GBM BTI. Given that most QPO detection methods are designed to find a single, strongest QPO candidate, this is unsurprising: removing most of the interval containing the 1.2 Hz QPO naturally leaves this second candidate as the strongest signal.

All three methods used have limitations and, as we have shown, make assumptions that are not met by the light curves analysed here. In particular, the assumption of pure stochastic (red) noise made by standard analyses in Fourier and wavelet domains hampers our ability to make robust detections in counter-intuitive ways: Hübner et al. (2022a) showed that many approaches may overconfidently detect QPOs in this context, because non-stationary light curves break the stationarity assumption of most methods. Wavelets – while excellent at detecting non-stationary, short-lived QPOs – cannot free us from this challenge if the assumed null hypothesis remains a stationary process. As theoretical modelling of GRBs improves, physically motivated predictions for the variability expected in prompt GRB emission would dramatically improve our robustness and sensitivity in finding these signals.

We quote significances for all three detections independently: combining detection probabilities across instruments or methods is generally not applicable in analyses containing substantial amounts of intrinsic source variability beyond detector noise, because the multiplication of probabilities relies on strict statistical independence of the underlying tests. This is given in the case of pure instrumental noise, but much more complicated when the underlying variability is intrinsic to the source emission, and thus shared across instruments (modulo some energy-dependence of the variability). However, we consistently find 3σ detections (or equivalent) in multiple methods for both QPOs, and thus conclude that both are very strong candidates.

Physically, the signature evidence of two QPOs at frequencies in the Hertz range can be explained in terms of a jet launching and evolution scenario. We note that different values for these frequencies might be obtained if different time windows could be probed with sufficient statistical precision in ν space, and the main information is that they are not in the kilo-Hertz range. The jet birth picture will focus on a binary merger progenitor, though the scenario described can also apply to a core collapse hypernova, should kilonova associations with long-duration bursts like GRB 230307A become more widely established.

As a binary NS merger proceeds, the tidal disruption that extracts plasma from one or both stars that eventually becomes the GRB jet only arises just before coalescence. This corresponds to an orbital semi-major axis of a few NS radii R_NS at most, and natural orbital frequencies of Ω = 2πν ∼ c/(3R_NS)∼10⁴ Hz for jet birth. As the merger proceeds to smaller radii, r, the Ω value increases modestly before the plasma shedding abruptly terminates at coalescence in presumably forming a black hole. Millisecond QPO periods from the merger proximity are likely mostly obscured by the plasma shroud, which thins out only at photospheric radii, R_ph ∼ 10¹² cm, for long-duration GRBs. Yet we note the recent report (Chirenti et al. 2023) of kilo-Hertz QPOs present in two short duration GRBs from the BATSE archive. This is an interesting result, albeit indicating the rarity of QPOs among the GRB population.

As the jet is launched, due to pressure from the surrounding medium comprising the disc and cocoon, it nominally develops a quasi-parabolic morphology (Tchekhovskoy et al. 2008; Komissarov et al. 2009) with an extraction of significant angular momentum, J, in the form of magnetic field helicity and plasma vorticity in a Poynting-flux dominated jet. The field maintains causal connection to the rotation (i.e. Ω) at the jet’s r ∼ R_NS base during the acceleration phase, wherein γ ∝ r^1/2 approximately (Tchekhovskoy et al. 2008). Eventually, the cocoon pressure declines and becomes insufficient to control the jet shape and dynamics, so that jet enters a coasting phase prior to its prime prompt GRB emission epoch. While the cocoon shapes the jet and controls its dynamics, it acts as the boundary to an acoustic cavity of circumference ∝r^1/2 transverse to the jet axis. Once the jet coasts, this boundary has effectively dissolved and rotational plasma fluctuations in the lateral dimension have no preferred timescale.

For the gas, at a light cylinder radius, typically quite close to the resultant compact object, the plasma vorticity saturates as the circular speed of the gas about the jet axis nudges c. Thereafter, as the jet expands to larger radii, r, its lateral (sheath) extent scales as r^1/2, while the circular speed is still c, so that the plasma rotation period, P, in the observer frame increases as P ∼ 2π(r R_NS)^1/2/c ∝ r^1/2. A ‘freeze-out’ of the plasma, P, occurs approximately when it enters the coasting epoch, which for fiducial jet launching models (see Fig. 8 of Tchekhovskoy et al. 2008) is around a resultant jet, γ = 10³, at around r/R_NS ∼ 10⁶, leading to P ∼ 0.1 − 1 s. This is consistent with the QPO periods observed herein, with the rotational freeze-out arising at distances of r ∼ 10¹² cm from the merger product, around the photospheric radius. More than one period may be sampled as the chaotic driver of the central engine interfaces with the cocoon sheath at the onset of the coasting and optically thin epoch.

In a quasi-acoustic phenomenon controlled by the cavity extent lateral to the jet direction, colliding plasma structures riding the gas vorticity in the jet would subsequently sample Fourier power at this rotational period as the jet becomes optically thin. Higher-frequency Fourier power from earlier epochs would generally be muted due to high photospheric densities. After dynamic decoupling of the jet from the cocoon, plasma fluctuations have no natural acoustic driver, and so their timescales decouple from further lateral expansion of the jet. Throughout, longitudinal fluctuations exist, but are not bounded geometrically, and so possess a chaos associated with the activity of the central driving engine, reflected in the light curves we see.

The QPOs observed thus constitute approximate images of the freeze-out plasma vorticity in a radially structured jet acquired at the larger distances, r_rad ∼ 10¹⁴ − 10¹⁶ cm, associated with prompt GRB radiation. As such, they enable jet archaeology by providing a window into the cessation of the jet launch and acceleration phase that is approximately contemporaneous with the jet’s exit from the photosphere. Accordingly, a core goal of jet launching simulations should be the reproduction of these rotational periods at the onset of coasting. The two QPO frequencies could be an imprint of the immediate pre-merger NS-NS binary evolution, carried forth to r_rad: higher-frequency signals correspond to jet plasma preparation deeper in the pre-merger gravitational potential where the Keplerian period is shorter. The observed values of P = 2π/Ω at ∼0.8 and ∼0.35 seconds do not distinguish whether the product of the merger is a black hole, a NS, or even a magnetar. They could well be analogous to year-timescale periods in flux levels observed from γ-ray blazars (Peñil et al. 2024) with their jets emanating from supermassive black holes.

Acknowledgments

We thank the anonymous reviewer for their helpful comments, Matteo Lucchini for advice on rms-energy spectra, and Volodymyr Savchenko for advice on dead time and saturation in INTEGRAL/SPI-ACS. DH is supported by NWO’s Women In Science Excel (WISE) fellowship. OJR gratefully acknowledges funding support from NASA under grant 80NSSC24M0035.

References

Appendix A: Supplementary figures

Fig. A.1. Left: Fourier periodogram of the Fermi/GBM NaI data with posterior draws from the three models compared via LRTs: in green, the power law model; in blue, a power law model with a Lorentzian component for a single QPO; in orange, a model comprising a power law and two Lorentzians. Middle: Distribution of the likelihood ratios from 1000 simulated periodograms. The likelihood ratio for the observed periodogram is a clear outlier. Right: Same as middle panel, but for the model with two QPOs. Again, the observed likelihood ratio is a clear outlier compared with the null hypothesis (a single QPO).

Fig. A.2. Left: Fourier periodogram of the Fermi/GBM BGO data with posterior draws from the three models compared via LRTs: in green, the power law model; in blue, a power law model with a Lorentzian component for a single QPO; in orange, a model comprising a power law and two Lorentzians. Middle: Distribution of the likelihood ratios from 1000 simulated periodograms. The likelihood ratio for the observed periodogram is a clear outlier. Right: Same as middle panel, but for the model with two QPOs. Again, the observed likelihood ratio is a clear outlier compared with the null hypothesis (a single QPO).

Fig. A.3. Posterior probability densities for all parameters in the DRW+QPO model for the SPI-ACS data described in Section 5.1.

Fig. A.4. Left panels: light curve (black) with BTI removed, maximum a posteriori model (blue) and posterior draws from the mean model (orange) for the Fermi/GBM NaI detector (na; top) and the BGO detector (b1; bottom). Right: corresponding posterior probability densities for the period in the model including a QPO. This figure suppements the results in Section 5.2.

All Tables

All Figures

Fig. 1. Top: Light curves of GRB 230307A in 50 ms temporal resolution as seen with SPI-ACS (left panel), and the brightest GBM NaI and BGO detectors; NaI 10 (middle panel) and BGO 1 (right panel). The time is in seconds since the trigger time. The orange area within the lower plots marks the time interval of 2.5 to 7.5 s, for which the GBM team issued a warning for possible data problems. Bottom: Fourier periodograms corresponding to the GRB light curves on the top. We show both the unbinned periodogram (black) and the log-binned periodogram (blue). Note that for the Fermi/GBM data, these do include the segment for which a warning was issued. All three periodograms contain strong variability above the instrumental noise limit at all frequencies considered here, and show peaks on top of the broadband variability present across all frequencies in the periodogram.

In the text

Fig. 2. Left: Fourier periodogram of the INTEGRAL data with posterior draws from the three models compared via LRTs: in green, the power law model; in blue, a power law model with a Lorentzian component for a single QPO; in orange, a model comprising a power law and two Lorentzians. Middle: Distribution of the likelihood ratios from 1000 simulated periodograms. The likelihood ratio for the observed periodogram is a clear outlier. Right: Same as middle panel, but for the model with two QPOs. Again, the observed likelihood ratio is a clear outlier compared with the null hypothesis (a single QPO).

In the text

	Fig. 3. Fractional rms amplitude as a function of photon energy for the two Fermi/GBM detectors.
In the text

Fig. 4. Left: INTEGRAL light curve of GRB230307A (black), with three random light curves generated from the stochastic process used in the QPO detection methods outlined in this section (orange). The GRB has a well-defined beginning and an end, in between which there exists rapidly changing variability. The simulated light curves also contain variability at a high amplitude, but the overall process does not change throughout the light curve. This is expected for a stationary stochastic process. Right: Periodogram of the INTEGRAL data (black) and of the simulated light curves (orange). While the periodogram of the GRB exhibits peaks formed by excess power in correlated neighbouring frequency bins, the periodograms of the stochastic process contain – by design and construction – powers that are statistically independent.

In the text

Fig. 5. Left: 2D wavelet transform (spectrogram) of the INTEGRAL observation of GRB230307A. The transform shows transient power at low frequencies in the first ∼20 seconds or so of the burst, where the variability is particularly strong. It also shows a short, transient signal between 0 − 10 s at 1.2 Hz, similarly to what was identified in the Fourier periodogram. The candidate signal at 0.32 Hz identified in the Fourier periodogram is less apparent here. Right: Fourier periodogram (black) and wavelet periodogram (blue) with the candidate signals found in the Fourier analysis noted as dashed orange lines. The Fourier and wavelet periodograms largely match.

In the text

	Fig. 6. Wavelet spectrograms for GRB 230307A for all three instruments. In all three, we find significant power at the lowest frequencies, as well as power at 1.2 Hz. In white contours, we overplot the 99.99% percentiles. In the Fermi/GBM data, we mark the segment flagged by the Fermi team for potential data issues as the shaded region.
In the text

Fig. 7. Wavelet periodograms (black) and red noise simulations (orange) for all three instruments (top: INTEGRAL; middle: Fermi/GBM NaI; bottom: Fermi/GBM BGO). The wavelet periodogram corresponds to the 2D wavelet spectrogram integrated over the time axis. In orange, the posterior mean derived from 1000 simulated wavelet periodograms, along with 10 posterior draws from the power-law stochastic model sampled in Section 3. In blue, we show the 99.99% single-trial detection limit. The candidate QPO at 1.2 Hz exceeds that limit for all three instruments, while the candidate QPO at 2.9 Hz exceeds this limit only for the Fermi/GBM BGO data. Note, however, that integrating over the time axis will necessarily yield lower significances given the short-lived nature of both candidate signals in the wavelet spectrogram.

In the text

	Fig. 8. Left: SPI-ACS light curve (black points) in units of counts per 0.05 s bin, the predicted rates from the maximum a posteriori combined model consisting of the mean function, the stochastic process and a QPO (blue), and posterior draws from the mean function (orange). Right: Posterior probability density for the QPO centroid period.
In the text

Fig. 9. Re-analysis of the SPI-ACS data excising the brightest part of the GRB (i.e. BTI), which is affected by data loss in Fermi/GBM. In the left and middle panel, we show observations in black, the maximum a posteriori combined model in blue, and draws from the posterior probability density for the mean function in orange. The left-most panel shows the results for ℳ1 (skew-Gaussian mean function and DRW), the middle panel for ℳ2 (skew-Gaussian mean function, DRW and QPO). In the right-hand panel, we show the posterior probability density for the period of the QPO, which is narrowly constrained around P = 0.347 s.

In the text

	Fig. 10. Left panels: Light curve (black), maximum a posteriori model (blue) and posterior draws from the mean model (orange) for the Fermi/GBM NaI detector (na; top) and the BGO detector (b1; bottom). Right: Corresponding posterior probability densities for the period in the model including a QPO.
In the text

	Fig. 11. Distribution of the Bayes factors obtained by modelling 100 simulated GRBs based on a model consisting of a skew-Gaussian trend function and a DRW, but without a QPO (grey). The Bayes factor obtained for the SPI-ACS data is a strong outlier compared to the simulations, indicating that our analysis captures additional variability not present in the simpler model.
In the text

	Fig. 12. Left: INTEGRAL SPI-ACS observation (black) with the maximum a posteriori model comprising a CARMA(2,1) process, a QPO and a skew-Gaussian mean function (blue), with posterior draws from the mean function in orange. Right: Posterior probability density of the period parameter shows a clear peak at 0.82 s, consistent with results obtained using the simpler DRW model.
In the text

	Fig. 13. Period posteriors for all three instruments for the full light curves (left) and for the light curves without the Fermi/GBM BTIs of 2.5–7.5 s after trigger (right). The posteriors show a remarkable degree of agreement between the three instruments, at frequencies corresponding to ∼1.2 Hz (0.83 s period) for the full light curves and at ∼2.9 Hz (0.34 s) without BTI.
In the text

Fig. A.1. Left: Fourier periodogram of the Fermi/GBM NaI data with posterior draws from the three models compared via LRTs: in green, the power law model; in blue, a power law model with a Lorentzian component for a single QPO; in orange, a model comprising a power law and two Lorentzians. Middle: Distribution of the likelihood ratios from 1000 simulated periodograms. The likelihood ratio for the observed periodogram is a clear outlier. Right: Same as middle panel, but for the model with two QPOs. Again, the observed likelihood ratio is a clear outlier compared with the null hypothesis (a single QPO).

In the text

Fig. A.2. Left: Fourier periodogram of the Fermi/GBM BGO data with posterior draws from the three models compared via LRTs: in green, the power law model; in blue, a power law model with a Lorentzian component for a single QPO; in orange, a model comprising a power law and two Lorentzians. Middle: Distribution of the likelihood ratios from 1000 simulated periodograms. The likelihood ratio for the observed periodogram is a clear outlier. Right: Same as middle panel, but for the model with two QPOs. Again, the observed likelihood ratio is a clear outlier compared with the null hypothesis (a single QPO).

In the text

	Fig. A.3. Posterior probability densities for all parameters in the DRW+QPO model for the SPI-ACS data described in Section 5.1.
In the text

	Fig. A.4. Left panels: light curve (black) with BTI removed, maximum a posteriori model (blue) and posterior draws from the mean model (orange) for the Fermi/GBM NaI detector (na; top) and the BGO detector (b1; bottom). Right: corresponding posterior probability densities for the period in the model including a QPO. This figure suppements the results in Section 5.2.
In the text