Quasi-periodic Oscillations of the X-Ray Burst from the Magnetar SGR J1935–2154 and Associated with the Fast Radio Burst FRB 200428

As the first X-ray astronomy satellite of China, Insight-HXMT was launched on 2017 June 15 and carries three main instruments on board (Zhang et al. 2020b): the High Energy X-ray telescope (HE; Liu et al. 2020), the Medium Energy X-ray telescope (ME; Cao et al. 2020), and the Low Energy X-ray telescope (LE; Chen et al. 2020). The high flux of the X-ray burst associated with FRB 200428 caused the data saturation and loss in HE and LE detectors. The raw data in some time intervals were discarded on board, and their light curves have several gaps as described in Li et al. (2021). Therefore, only the data of the ME instrument are used to search for the QPO, and ME also provides the high time resolution (6 us) required for timing analysis (Tuo et al. 2021).

The data reduction of the X-ray burst associated with FRB 200428 is performed with the Insight-HXMT Data Analysis Software package (HXMTDAS) version 2.04. The main steps of ME data reduction are as follows: (1) Use the commands mepical in HXMTDAS to calibrate the photon events from the 1L data product according to the Calibration Database (CALDB) of Insight-HXMT. (2) Select the good time intervals (GTIs) from T₀ to T₀+1 s, where T₀ is 2020-04-28 14:34:24 UTC (satellite time). (3) Identify the grade of each event according to their arrival time and calculate the dead time of each Field Programmable Gate Array (FPGA) at the specified interval. (4) Extract the good events based on the GTIs using the commands mescreen. (5) Produce the light curve of ME with the dead-time correction using the commands, melcgen. Although ME detectors do not have the problems of the data saturation like HE and LE, the dead-time effect is still significant for bright X-ray burst detection.

We elaborate the dead-time effect of the ME instrument here. ME consists of three detector boxes. Each box has three independent FPGAs, and each FPGA manages six Application Specified Integrated Circuits (ASICs). Each ASIC is responsible for the readout of 32 pixels (Li et al. 2020). The dead time of ME can be calculated by each FPGA independently. When 1 pixel in the FPGA is triggered, a window with 62.5 ns is opened and the pixels triggered in other ASICs of the same FPGA will be recorded. If there are N_hit ASICs in the same FPGA that are triggered within the time window, the electronics will take 163 + 93 · N_hit μs to process these triggered pixels. The typical timescale of dead time is 256 μs because only 1 pixel is triggered at most situations (Cao et al. 2020). The dead-time process of ME is similar to that of nonparalyzed detectors, and the hardware algorithm is implemented in the HXMTDAS to calculate the dead time of ME. The dead time causes the loss of photons and the distribution of photon arrivals to deviate from Poisson distribution and alters the distribution of power in the periodogram as well. Note that the dead time strongly depends on the count rate. As a consequence, dead-time corrections are especially important for bright bursts.

Special types of events in ME data could be an indicator to determine whether the dead-time correction is performed well. There are two kinds of thresholds for ME. One is the triggered threshold; the other is the readout threshold. The special events are those where the pulse height is higher than the triggered threshold and lower than the readout threshold. These events are mainly the noise events that are independent of the source brightness, so the light curve of these events is expected to be constant after dead-time correction as shown in Figure 1.

Zoom In Zoom Out Reset image size

KW consists of two identical NaI(Tl) scintillation detectors, each with a 2π sr field of view. One detector (S1) points toward the south ecliptic pole, thereby observing the south ecliptic hemisphere, while the other (S2) observes the north ecliptic hemisphere. Here the light curves from KW are obtained from Ridnaia et al. (2021).

2.3.1. Monte Carlo Simulations of Pulse Profile

Monte Carlo simulation of light curves is a standard tool in timing analysis. One can fit an empirical function to the burst profile and then generate a large number of realizations of that burst profile by adding the Poisson photon counting noise. The periodograms computed from these simulated light curves form a distribution against which to compare the periodogram of the real data. The X-ray burst profiles of FRB 200428 are fitted with three Gaussian functions, which can be found in Li et al. (2021):

$$\begin{eqnarray}&&{R}_{3}=\displaystyle \sum _{i=1}^{3}{N}_{p,i}G(t,{t}_{p,i},{\sigma }_{p,i})+l,\end{eqnarray} \tag{ 1 }$$

where $G(t,{t}_{p,i},{\sigma }_{p,i})=\frac{1}{\sqrt{2\pi }{\sigma }_{p,i}}\exp \left(-,\tfrac{{\left(t-{t}_{p,i}\right)}^{2}}{2{\sigma }_{p,i}^{2}}\right)$; N_p,i , t_p,i , and σ_p,i are the normalization, arrival time, and Gaussian width of the ith peak, respectively; and l is the background level of the light curve.

Actually, the true pulse profile of the X-ray burst associated with FRB 200428 cannot be known precisely owing to the degeneracy between the overall burst profile and the red noise component. If the pulse profile can be presented by Equation (1), Monte Carlo simulations of the pulse profile caused by the Poisson noise of photon counting can be used to estimate the significance of the QPO candidate.

2.3.2. Bayesian Method

In the following, we briefly introduce the Fourier-based periodograms analyzed by the Bayesian method commonly used to derive the significance of QPOs (Vaughan 2010; Huppenkothen et al. 2013). A periodogram is a sample of the power spectral density (PSD) of the signal based on a given time series x_k (t) (k = 0, 1, ..., N − 1), N = 2^m (m is an integer). We calculate the periodogram as the absolute square of the discrete Fourier transform of the signal, and it follows a ${\chi }_{2}^{2}$ distribution around the PSD S_j = S(f_j ), j = 0, 1, N/2 at Fourier frequency f_j (Van der Klis 1989; Timmer & Koenig 1995; Vaughan 2010). Therefore, the periodogram I_j follows an exponential distribution about S_j and can be described as follows:

$$\begin{eqnarray}&&p({I}_{j}| {S}_{j})=\displaystyle \frac{1}{{S}_{j}}{e}^{-{I}_{j}/{S}_{j}}.\end{eqnarray} \tag{ 2 }$$

Many astrophysical transients show excess power at low frequencies, and it is often assumed that this can be modeled as a red-noise process. One basic noise model S_p (f) that is commonly used is a combination of a red-noise power law and a white-noise amplitude,

$$\begin{eqnarray}&&{S}_{p}(f)={{Af}}^{-\alpha }+C,\end{eqnarray} \tag{ 3 }$$

where A is the amplitude and α is the index. C is the constant representing the white-noise component. Another noise model S_b (f) is composed of a broken power law and a constant white noise,

$$\begin{eqnarray}{S}_{b}(f)=\left\{\begin{array}{ll}{{Bf}}^{-{\alpha }_{1}}+C & \mathrm{if}\,f\leqslant {f}_{\mathrm{break}},\\ {{Bf}}_{\mathrm{break}}^{{\alpha }_{2}-{\alpha }_{1}}{\left(f\right)}^{-{\alpha }_{2}}+C & \mathrm{if}\,f\gt {f}_{\mathrm{break}},\end{array}\right.\end{eqnarray} \tag{ 4 }$$

where α₁ and α₂ are indices of the two power-law components, f_break is the break point for the frequency, and B is the normalization value. When α₁ equals α₂, the model of S_b (f) will be the same as S_p (f), and they satisfy the requirement of nested models. Following Bayes's rule, the posterior probability of a set of model parameters θ, like A, α, and C in noise model S_p (f), can be defined given the uniform distribution of the prior.

With the Bayesian method, we obtain the maximum a posteriori (MAP) estimates of the model parameters by maximum likelihood estimation (MLE) of observed Leahy-normalized periodogram. As introduced in Vaughan (2010), maximizing the joint likelihood function $p({\boldsymbol{I}}| {\boldsymbol{\theta }},H)=\prod _{j}^{N/2}p({I}_{j}| {S}_{j})$ is equivalent to minimizing the deviance function:

$$\begin{eqnarray}&&D({\boldsymbol{I}},{\boldsymbol{\theta }},H)=-2\mathrm{log}p({\boldsymbol{I}}| {\boldsymbol{\theta }},H)=2\sum _{j=1}^{N/2}\left(\displaystyle \frac{{I}_{j}}{{S}_{j}}+\mathrm{log}{S}_{j}\right).\end{eqnarray} \tag{ 5 }$$

In other words, finding the MLEs of the parameters of a model S( θ ) can be expressed as

$$\begin{eqnarray}&&{\hat{{\boldsymbol{\theta }}}}_{\mathrm{MLE}}=\arg \mathop{\min }\limits_{\theta }D({{\boldsymbol{I}}}^{\mathrm{obs}},{\boldsymbol{\theta }},H)\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\,=\arg \mathop{\max }\limits_{\theta }p({\boldsymbol{I}}={{\boldsymbol{I}}}^{\mathrm{obs}}| {\boldsymbol{\theta }},H),\end{eqnarray} \tag{ 7 }$$

where I ^obs is the observed periodogram, $\mathrm{argmin}$ is to find the set of parameters θ to obtain the minimum value of function D, and $\mathrm{argmax}$ of distribution p is to find the maximum likelihood.

To verify the possible QPO in the periodogram, a model selection task is taken to assess whether the periodogram can be represented by a simple power law or requires a more complex model like a broken power law. The likelihood ratio test (LRT) is employed to distinguish two nested models, null hypothesis H₀ and the alternative hypothesis H₁. H₀ is the power law, and H₁ is the broken power law. The LRT statistic is twice the minimum log likelihood of the two models:

$$\begin{eqnarray}&&{T}_{\mathrm{LRT}}=-2\mathrm{log}\displaystyle \frac{p\left({\boldsymbol{I}}| {\hat{{\boldsymbol{\theta }}}}_{\mathrm{MLE}}^{0},{H}_{0}\right)}{p\left({\boldsymbol{I}}| {\hat{{\boldsymbol{\theta }}}}_{\mathrm{MLE}}^{1},{H}_{1}\right)}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\,={D}_{\min }({H}_{0})-{D}_{\min }({H}_{1}).\end{eqnarray} \tag{ 9 }$$

Once the noise model is selected using the LRT statistic, the QPO search can also be considered as a model selection problem. The new H₁ model is the sum of the selected noise model and a Lorentz model. When the amplitude of the Lorentz is zero, the H₁ model becomes the noise model. The asymptotic theory suggests that the LRT statistic should be distributed as a chi-square variable under certain regularity conditions, ${T}_{\mathrm{LRT}}\sim {\chi }_{\nu }^{2}$, where the number of degrees of freedom ν is the difference between the number of free parameters in H₁ and H₀. When the regularity conditions are not satisfied, one cannot use the distribution to be that of the asymptotic theory (Protassov et al. 2002). Nevertheless, LRT is still a powerful statistic for comparing models and can be calibrated by posterior predictive simulations, as shown by Protassov et al. (2002).

In order to investigate narrow features, another statistic T _R to use is the maximum ratio of observed to model power (Vaughan 2010; Huppenkothen et al. 2013),

$$\begin{eqnarray}&&{T}_{R}=\max ({R}_{j}),\end{eqnarray} \tag{ 10 }$$

where

$$\begin{eqnarray}&&{R}_{j}=2{I}_{j}/{S}_{j},\end{eqnarray} \tag{ 11 }$$

with I_j and S_j the observed and model powers. This is the candidate for single-bin periodicity when searching the periodogram for the highest data/model outlier. After simulating a large number of periodograms from the posterior parameter sets and fitting each simulated periodogram with the preferred noise model, we can also compute the data/model and find the maximum data/model outlier. Then, the distributions of the highest outlier can be calibrated under the selected noise model. The posterior p-value for the observed outlier can be obtained when calculating the tail areas from the simulated outlier distributions.

Although Monte Carlo simulation of the pulse profile is an important tool in timing analysis, it must know the physically motivated burst profile function in advance. We assume that the burst profile of the X-ray burst associated with FRB 200428 can be fitted with three Gaussian functions, and the fit parameters can be found in Table 1. It should be noted that we just take the burst as an example to show that an assumption for the burst profile without taking the red noise into account will produce an uncorrected significance estimation of the QPO candidate.

The procedures of the Monte Carlo simulations of the burst profile are described below. As introduced in Section 2.3.1, the burst profile can be fitted with Equation (1) and is shown in Figure 2 with a green solid line. The Leahy-normalized periodograms of the light curve can be derived using the FFT code in the Python software stingray as shown in Figure 3. A potential QPO candidate with moderate significance at a frequency of about 40 Hz can be found. In order to establish the p-value of the QPO, Monte Carlo simulations of the burst profile are used to validate whether the QPO signal is generated from the fluctuations of the light curve or not. We then generate 10,000 simulated light curves sampled from the fitted three-Gaussian burst profile by adding the Poisson photon-counting noise. The periodogram is computed for each simulated light curve to form a distribution at each frequency bin. The maximum, mean, and minimum values of power distribution at each frequency bin are calculated and plotted in Figure 3. Comparing the observed power in each frequency bin with the distribution of simulated powers allows us to make a statement that the observed power at 40 Hz is difficult to produce by the Poisson noise process alone according to the distribution of simulated powers. From the simulated distribution at each frequency, the distribution of Leahy power above 10 Hz satisfies a ${\chi }_{2}^{2}$ distribution. Consequently, the p-value for the QPO at 40 Hz under the assumption of the three-Gaussian burst profile is 1.7e–7, indicating that the Poisson noise of the profile with only three broad Gaussians cannot explain the observed excess in power at 40 Hz.

Zoom In Zoom Out Reset image size

The QPO candidate in the light curve as shown in Figure 2 can also be fitted with

$$\begin{eqnarray}&&Q={N}_{g}G(t,{t}_{g},{\sigma }_{g})\sin (2\pi {f}_{q}(t-{t}_{q}))+{R}_{3},\end{eqnarray} \tag{ 12 }$$

where N _g , t _g , and σ _g are the normalization, arrival time, and width of the broad Gaussian profile of the X-ray burst; f _q is QPO frequency; and t _q is the time offset of the sinusoidal function. The fitting results are N _g = 1.4 ± 0.2, t _g = 0.452 ± 0.004, σ _g = 0.024 ± 0.004, f _q = 40.4 ± 1.2, t _q = 0.09 ± 0.01. The frequency f _q is consistent with the result from the periodogram analysis presented below. The lower and upper envelopes of the amplitude modulation of the QPO candidate can thus be expressed by

$$\begin{eqnarray}&&{Q}_{\mathrm{en}}={R}_{3}\pm {N}_{g}G(t,{t}_{g},{\sigma }_{g}).\end{eqnarray} \tag{ 13 }$$

The fitted ME light curve with Equation (12) is shown in Figure 2(b), where the lower and upper envelopes of the amplitude modulation of the QPO described by Equation (13) are also plotted.

The Monte Carlo method outlined above is versatile and powerful, but it has limitations. In this section, we will explore the Bayesian method, which is a conservative method based on the assumption that the red-noise component dominates the periodogram. This approach is to model the observed periodogram rather than the light curve. The data between 0.1 and 0.74 s are used to generate the dead-time-corrected light curve with time resolution of 2.5 ms and the corresponding periodogram. The 0.64 s duration basically covers the entire X-ray burst of FRB 200428 as shown in Figure 2. In order to study the nonstationary effect on this burst, we also compare the results of different duration times such as the main peak of the burst (0.32 s) and more background data besides the burst such as 1.28 s.

3.2.1. Selection of a Noise Model

In the following, we briefly illustrate the analysis procedure with the X-ray burst associated with FRB 200428. This burst has a duration of T₉₀ with 0.65 s, and the periodogram of 0.64 s (from 0.1 to 0.74 s) length with time resolution of 2.5 ms for this burst is presented in Figure 4. After fitting both a simple power law (H₀ model) and the broken power law (H₁ model), the likelihood ratio between the two models using Equation (8) is calculated as LRT = 5.36. The fits to the periodogram and the ratios (data/model) are presented in Figure 4 using the MAP estimate of the parameters.

Zoom In Zoom Out Reset image size

We utilize the Markov Chain Monte Carlo (MCMC) simulations implemented in the Python package emcee (Foreman-Mackey et al. 2013) to generate sets of parameters of posterior probabilities. We then sample from the posterior distribution of the power-law model via MCMC and simulate the periodograms according to Equation (2) from chains of posterior distribution as shown in Figure 5. These simulated periodograms are again fitted with models H₀ and H₁, respectively, so that we can build the likelihood ratio distribution of the power-law model. Comparing the distribution of simulated likelihood ratios with that of the observed LRT, we can compute the tail area probability or the posterior p-value as 0.03, which is shown in Figure 6. Thus, we can reject the power-law model and accept the broken power-law model for the noise model.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.2.2. Searching for the QPO

The noise model has been selected above as a broken power law, and then we can search for the QPO as a model selection problem. We compare the selected noise model to a more complex model combining both the noise model and a Lorentz function to account for the QPO. The burst periodograms are fitted by bknpow+constant as a simple model (H₀) and can be also fitted by bknpow+Lorentz+constant as a nested model (H₁). bknpow and Lorentz are the power-law function with a break frequency and the Lorentz function in Xspec. After fitting the periodogram with the two models above, the likelihood ratio between the two models can be obtained, denoted as ${T}_{\mathrm{LRT}}^{\mathrm{obs}}$. The ${D}_{\min }({H}_{0})$ and ${D}_{\min }({H}_{1})$ of the observed burst periodogram are 743.75 and 712.76, respectively. The ratio of the periodogram and H₁ model is also displayed in the bottom panel of Figure 4. Thus, the observed LRT of the burst is ${T}_{\mathrm{LRT}}^{\mathrm{obs}}=30.99$.

Next, we will clarify how we calibrate the posterior distribution of LRT and calculate the p-value. Using the H₀ model, we pick n sets of parameters for the H₀ model from the posterior MCMC chains. Each set of parameters is simulated to produce a periodogram, for which an LRT value for H₀ and H₁ is computed. We then have n LRT values for n periodograms. The distribution of the n LRT values is the posterior distribution of LRT, which is shown in Figure 7 with a red solid line. The observed LRT value is also shown in Figure 7, with a red vertical line. As a result, the p-value or the ratio of the tail areas of ${T}_{\mathrm{LRT}}^{\mathrm{obs}}$ is 1e–6, which corresponds to 4.75σ. We find a posterior centroid frequency for the QPO of ${39.3}_{-1.20}^{+1.22}$ Hz, with a width of ${3.7}_{-2.01}^{+3.74}$ Hz. We use the posterior distribution of QPO parameters to derive the fractional rms amplitude A_frac and find A_frac = 0.36 ± 0.096. The results of 0.64 s duration are also displayed in Table 2.

Zoom In Zoom Out Reset image size

Table 2. Results for the Significance of QPO and QPO Parameters at Different Duration Times Using the LRT Statistic

(s)(Hz)(Hz)

Duration Time	H0 Model	H1 Model	${D}_{\min }({H}_{0})$	${D}_{\min }({H}_{1})$	${T}_{\mathrm{LRT}}^{\mathrm{obs}}$	p-value	Frequency	Width	Norm
0.32	PL	PL+Lorentz	378.91	359.07	19.84	2.9e–4	${39.14}_{-1.94}^{+1.88}$	${4.88}_{-2.73}^{+5.67}$	${389.15}_{-167.17}^{+260.55}$
0.64	BPL	BPL+Lorentz	743.75	712.76	30.99	1e–6	${39.28}_{-1.20}^{+1.22}$	${3.67}_{-2.01}^{+3.74}$	${315.27}_{-116.74}^{+218.63}$
1.28	BPL	BPL+Lorentz	1479.71	1421.63	58.08	<5e–7	${39.26}_{-0.8}^{+0.78}$	${3.24}_{-1.41}^{+2.02}$	${267.95}_{-76.85}^{+135.87}$

Download table as:  ASCIITypeset image

The analysis presented above makes a strong assumption about the data: our choice of a ${\chi }_{2}^{2}$ distributed likelihood around the model power spectrum implies that the periodogram is the result of a stationary noise process. This assumption is not strictly the case for burst signals, as shown in Huppenkothen et al. (2013); it is a conservative assumption that holds for all but the lowest frequencies in the periodogram. In this paper, we also explore the effects of the nonstationary nature of a transient light curve and a QPO within it. We show that selecting more background signals before or after the burst signal leads to greatly overestimating the significance of potential QPOs. We also show that the significance will be overestimated if the whole burst signal is selected for analysis but the QPOs occur in only part of the burst.

To search for the duration of the 40 Hz QPO, we utilize wavelet analysis and the two most widely used methods: the generalized Lomb–Scargle periodogram (LSP; Lomb 1976; Scargle 1982) and weighted wavelet Z-transform (WWZ; Foster 1996) are used to obtain the power spectra of the 0.64 s light curve with time resolution of 2.5 ms. In this work the power spectra of the LSP method are checked with the independent results of the WWZ approach. As shown in Figure 8, the QPO signal appears in the main peak of the time range from 0.3 to 0.6 s. Then, we choose the duration of the main peak to be about 0.32 s and more background before and after the main peak, such as 1.28 s, to compare their differences.

Zoom In Zoom Out Reset image size

Using the same method for 0.64 s, we select the noise model for 0.32 and 1.28 s. The H₀ model is a power law for 0.32 s and a broken power law for 1.28 s. As shown in Table 2, for each duration time, after selecting the noise model H₀, we use the H₀ model to generate a large set of parameters (for power law, three parameters as shown in Equation (3), and for broken power law, five parameters as shown in Equation (4)) from their posterior distributions. Each parameter set is used to simulate a periodogram using the exponential distribution in Equation (2). For the simulated periodogram we also use the H₀ model and H₁ model to get the ${D}_{\min }$ values as defined in Equation (5) using MCMC. The T_LRT can be derived from Equation (8). Then, the distribution of the simulated likelihood ratios can be derived from multiple sets of parameters of the H₀ model. Comparing the distribution of simulated likelihood ratios with that of the observed T_LRT as shown in Figure 7, we can compute the tail area probability, and this is also the normalized posterior p-value. A total of 150,000 and 2,000,000 sets of parameters for the H₀ model from the posterior Markov chain are used to simulate the periodogram for 0.32 and 1.28 s, respectively. The likelihood ratio distributions and the observed LRT values are also displayed in Figure 7. The posterior parameters of the QPO, such as the centroid frequency, width, and normalization, are also shown in Table 2.

The smaller p-value shown in Table 2 indicates a higher significance of the QPO. Although the parameters of the QPO are consistent within the uncertainties in the three durations, however, the uncertainties of the QPO parameters are smaller for the longer duration. Our explanation about why the significance is higher for a long duration used is as follows. We find that the QPO is only present in the main peak of the X-ray burst. However, the noise level varies through the time series, i.e., the noise level is higher during the main peak, where the count rate is higher (and also varies) than that outside the main peak. This is a typical nonstationary process. Extending the selected duration beyond the main peak is equivalent to using the noise level beyond the main peak to estimate the noise level during the main peak, which does not change the strength of the QPO but decreases the estimated noise. Consequently, the signal-to-noise ratio and detection significance of the QPO are artificially increased.

To test this explanation, we have also done simulations for 0.64 and 1.28 s and assume that the light curve before and after the main peak (0.32 s) are only Poisson processes, i.e., by replacing the observed light curve with the Poisson sampled light curve of the same count rate. We have obtained the same results: a longer duration outside the main peak leads to a higher significance of the QPO. Therefore, a reliable way to estimate the significance of the QPO is to use the segment of data in which the QPO is present.

As a consequence, we report the final results of QPO from 0.32 s, which are shown in Table 2.

As described in Section 2.3.2, another statistic for investigating narrow QPOs is T _R , which is defined in Equation (10). We also search the periodogram for the highest data/model outlier and compare this outlier to those distributed by pure noise to find narrow features that may be candidates for a possible QPO. One shortcoming of the T _R statistic is that it optimally detects periodic signals confined to one frequency bin. Since the QPO power spreads over several bins, this is not an optimal way of detecting broad signals. There are several ways to reduce this restriction. One is to bin (or smooth) the data in some way and compute T _R for the binned data. If we bin the simulated periodograms in the same way, then the test statistic T _R for the binned data is comparable to the distribution approximated by our simulations, and the latter can be used to derive posterior predictive p-values.

We use the same chosen noise model (H₀ model) to fit the periodogram and the same samples of simulated periodograms in the previous step to search for the highest data/model outlier. We bin the periodogram geometrically, where the bin size grows with frequency and the bin factor is shown in Table 3. Figure 9 also presents as an example for the bin of the 1.28 s periodogram, and we choose the bin factor as 0.25. We also compare the observed value of the T _R statistic and compute the corresponding p-values for the three durations of the burst. The results are summarized in Table 3. The QPO with central frequency at about 40 Hz can be found for different durations. The logarithmic rebin of the periodogram cannot ensure that all the QPO frequency bins are in the same bin; therefore, this method is just used as a verification of the QPO.

Zoom In Zoom Out Reset image size

To further verify the robustness of the ∼40 Hz QPO detection, cross-correlations between the detrended light curves of Insight-HXMT/ME and KW are made and shown in Figure 10, along with cross-correlations between each light curve and white noise. The QPOs in both light curves match each other in both phase and period of 24.5 ms, while no statistically meaningful correlation is found with white noise.

Zoom In Zoom Out Reset image size

A cross-correlation function (CCF) is used here to understand whether the QPO signals found in the light curves of Insight-HXMT/ME and KW are synchronized. The CCF between the detrended light curves (after removing the three broad Gaussian peaks by fitting with Equation (1)) of ME and KW shows strong oscillations of ∼24.5 ms from τ = − 0.1 s to τ = 0.05 s as plotted in Figure 10. The errors for CCF are estimated from the Monte Carlo method, by sampling the light curves 1000 times. We also make the cross-correlation between the detrended light curves and white noise with the same average counting rate; as expected, there are no oscillation structures in CCFs. However, there are some oscillation signals around −0.15 and 0.15 s, which are also generated by the oscillation structures in the real light curves.

To further understand the above results of cross-correlation analysis, we use sinusoidal signals to simulate the structures in CCFs obtained above. As shown in Figure 11(a) (background + sinusoidal signals), two sinusoidal signals of amplitudes of 8 and 24 are simulated with background 24, which have frequency 10/2π Hz and only appear for 3 s in the light curves. The Poisson fluctuations are added on the light curves. The background light curves only including white noise are generated from Poisson sampling with an expected count of 24 in each time bin, as shown in Figure 11(b) (background only) for two of the simulated background light curves of white noise.

Zoom In Zoom Out Reset image size

We first make cross-correlations between the two simulated light curves (including their backgrounds) with different amplitudes of sinusoidal oscillations. Very significant and periodic structures are shown in Figure 12, demonstrating the power of cross-correlation analysis in identifying weakly oscillating signals ("Seq 2"; blue line in Figure 11(a)), given a template of the oscillations ("Seq 1"; red line in Figure 11(b)). Therefore, the cross-correlation analysis presented in Figure 10 is equivalent to using ME's light curve, which contains strong QPO signals, as a template in identifying the relatively weak QPOs in the light curve of KW.

Zoom In Zoom Out Reset image size

The CCFs between the white noise and both light curves with weak and strong oscillation signals are shown in Figure 13. Panel (a) shows the result from one single simulation run, where many oscillatory structures exist. The processes are repeated 1000 times to obtain the mean CCFs, which show clear oscillation structures around ±5 (Figure 13(b)). These structures are similar and should have the same origin as the structures near both ends of the correlation functions of white noise with the light curves of ME and KW, as shown in Figure 10.

Zoom In Zoom Out Reset image size

The 40 Hz QPO reported in this paper poses the following constraints on the available models for FRB 200428:

• 1.  
  The fact that this burst is special with respect to the majority of SGR bursts suggests that the rarity of FRB–SGR burst associations (Lin et al. 2020c) is likely due to intrinsic rather than extrinsic (beaming or narrow spectra; Lin et al. 2020c) factors.
• 2.  
  Since QPOs are directly related to crust oscillations, and since the two spikes on the X-ray light curve that are associated with the two pulses of FRB 200428 are roughly consistent with QPO peaks, it is almost certain that the X-ray emission and the associated FRB pulses originate from the magnetosphere of the neutron star (Katz 2020; Lu et al. 2020; Wang et al. 2020; Wadiasingh & Chirenti 2020; Wang 2020; Yang et al. 2020). Models invoking emission outside the magnetosphere (Yu et al. 2020; Margalit et al. 2020) are thus disfavored, unless the ejections of the emitter are associated with crust activities (Yuan et al. 2020).
• 3.  
  The trigger of QPOs in the neutron star crust is likely internal, but an external trigger powered by infall of an asteroid is also possible (Dai 2020; Geng et al. 2020), which might explain the rarity of the FRB-associated bursts. Such models, however, need to interpret the existence of QPOs already before the spikes associated with FRB pulses.

The identified QPO at 40 Hz may be explained as the toroidal mode of the crust (Duncan 1998), or the magnetoelastic modes if the coupling effect between the crust and core is considered (van Hoven & Levin 2011). For the crustal toroidal mode, the 40 Hz QPO can be explained as the ₃ t₀ mode if the tangled magnetic field in the crust is B_t ≪ B_μ , or the fundamental toroidal mode ₂ t₀ if B_t ≈ B_μ , where _ℓ t_n is the notation for the general toroidal modes, ℓ is the angular quantum number, n is the number of radial nodes in the eigenfunctions, and ${B}_{\mu }={\left(4\pi \mu \right)}^{1/2}\approx 4\times {10}^{15}{\rho }_{14}^{0.4}$ G, with μ being the crust shear modulus and ρ₁₄ being crust density in units of 10¹⁴ g cm⁻³ (Duncan 1998). Note that the internal magnetic field is usually expected to be much higher than the surface dipole magnetic field, which is 2.2 × 10¹⁴ G for SGR J1935+2154 (Israel et al. 2016).