Maximum likelihood fitting of X-ray power density spectra: Application to high-frequency quasi-periodic oscillations from the neutron star X-ray binary 4U1608-522

High frequency quasi-periodic oscillations (QPOs) from weakly magnetized neutron stars display rapid frequency variability and high coherence with quality factors up to at least 200 at frequencies around 850 Hz. Their parameters have been estimated so far from standard min(chi2) fitting techniques, after combining a large number of Power Density Spectra (PDS), as to have the powers normally distributed. Accounting for the statistical properties of PDS, we apply a maximum likelihood method to derive the QPO parameters in the non Gaussian regime. The method presented is general, easy to implement and can be applied to fitting individual PDS, several PDS simultaneously or their average, and is obviously not specific to the analysis of kHz QPO data. It applies to the analysis of any PDS optimized in frequency resolution and for low frequency variability or PDS containing features whose parameters vary on short timescales, as is the case for kHz QPOs. It is equivalent to the standard chi^2 minimization fitting when the number of PDS fitted is large. The accuracy, reliability and superiority of the method is demonstrated with simulations of synthetic PDS. We show that the maximum likelihood estimates of the QPO parameters are asymptotically unbiased, and have negligible bias when the QPO is reasonably well detected. By contrast, we show that the standard min(chi2) fitting method gives biased parameters with larger uncertainties. The maximum likelihood fitting method is applied to a subset of archival Rossi X-ray Timing Explorer (RXTE) data of the neutron star X-ray binary 4U1608-522. We show that the kHz QPO parameters can be measured on 8 second timescales and that the time evolution of the frequency is consistent with a random walk. This enables us to estimate the intrinsic quality factor of the QPO to be around 260, whereas previous analysis indicated a maximum value around 200 (abridged).

💡 Research Summary

The paper addresses a fundamental methodological limitation in the analysis of X‑ray timing data, specifically the fitting of power density spectra (PDS) that contain high‑frequency quasi‑periodic oscillations (kHz QPOs) from weakly magnetized neutron stars. Traditional approaches rely on averaging a large number of PDS so that the resulting powers become approximately Gaussian, allowing the use of a minimum‑χ² fitting routine. This averaging, however, erases short‑timescale variability of the QPO parameters (frequency, width, amplitude) and introduces biases when the number of spectra is modest.

To overcome these issues, the authors develop a maximum‑likelihood estimation (MLE) framework that directly incorporates the true statistical distribution of PDS powers, namely the exponential (χ² with two degrees of freedom) distribution derived from the underlying Fourier amplitudes. The model for a single PDS consists of a flat white‑noise level plus a Lorentzian component describing the QPO, characterized by four parameters: centroid frequency ν₀, full‑width at half‑maximum (FWHM) Δν, amplitude (integrated power) A, and background level C. The likelihood for each frequency bin is the exponential probability density evaluated at the observed power given the model expectation; the total log‑likelihood is the sum over all bins. Maximization is performed with standard numerical optimizers (restricted BFGS, Newton‑Raphson), and parameter uncertainties are obtained from the inverse Fisher information matrix, providing asymptotically correct error estimates.

The authors validate the method through extensive Monte‑Carlo simulations. Synthetic PDS are generated with a range of signal‑to‑noise ratios, observation lengths (from 8 s to 256 s), and numbers of averaged spectra. The results demonstrate that MLE yields essentially unbiased estimates of ν₀, Δν, and A, even when the QPO is weak or when only a few spectra are available. In contrast, the conventional χ² fitting systematically underestimates the width and overestimates the uncertainties, leading to a biased quality factor Q = ν₀/Δν. The bias of the MLE becomes negligible once the QPO is detected at a modest significance (≈5σ). Moreover, the MLE can be applied to individual PDS without any averaging, preserving the intrinsic rapid frequency drift of the oscillation.

The method is then applied to archival Rossi X‑ray Timing Explorer (RXTE) observations of the neutron‑star low‑mass X‑ray binary 4U 1608‑522. The data set comprises several hundred kiloseconds of high‑time‑resolution event mode recordings. The authors segment the light curves into 8‑second intervals, compute a PDS for each segment, and fit each spectrum with the MLE procedure. This yields a time series of QPO frequencies that exhibits a random‑walk behaviour on timescales of seconds to minutes, consistent with previous studies that used longer averaging windows. By correcting for the observed frequency drift, the intrinsic quality factor of the QPO is inferred to be Q ≈ 260, significantly higher than the previously reported maximum of Q ≈ 200 obtained with χ² fitting. The higher Q implies that the oscillation is more coherent than earlier believed, placing tighter constraints on theoretical models of kHz QPO generation (e.g., relativistic precession, disk‑oscillation modes, or beat‑frequency mechanisms).

Beyond the specific case of kHz QPOs, the authors emphasize the generality of the MLE approach. It can be employed for any PDS analysis where the frequency resolution is optimized and where spectral features evolve on short timescales, such as low‑frequency noise components, broad Lorentzians, or transient burst oscillations. The technique is computationally efficient, easy to implement, and yields statistically optimal parameter estimates without the need for large ensembles of spectra. Consequently, it is poised to become a standard tool for forthcoming X‑ray timing missions (e.g., eXTP, STROBE‑X) that will deliver high‑throughput, high‑resolution data where preserving short‑timescale variability is essential.

In summary, the paper introduces a robust maximum‑likelihood fitting methodology for X‑ray power spectra, demonstrates its superiority over traditional χ² fitting through simulations, and applies it to real RXTE data to reveal a higher intrinsic coherence of the kHz QPO in 4U 1608‑522. The work represents a significant advance in timing analysis techniques and opens new avenues for probing the physics of neutron‑star accretion flows.