These authors contributed equally to this work.

These authors contributed equally to this work.

[2,5]\fnmLijing \surShao

1]\orgdivDepartment of Astronomy, School of Physics, \orgnamePeking University, \orgaddress\streetNo. 5 Yiheyuan Road, \cityBeijing, \postcode100871, \stateBeijing, \countryChina

2]\orgdivKavli Institute for Astronomy and Astrophysics, \orgnamePeking University, \orgaddress\streetNo. 5 Yiheyuan Road, \cityBeijing, \postcode100871, \stateBeijing, \countryChina

3]\orgdivSchool of Physics, \orgnameDalian University of Technology, \orgaddress\streetNo. 2 Linggong Road, \cityDalian, \postcode116024, \stateLiaoning, \countryChina

4]\orgdivInstitute for Gravitational Wave Astronomy, \orgnameHenan Academy of Sciences, \orgaddress\streetNo. 228 Chongshili Road, \cityZhengzhou, \postcode450046, \stateHenan, \countryChina

5]\orgdivNational Astronomical Observatories, \orgnameChinese Academy of Sciences, \orgaddress\street20A Datun Road, \cityBeijing, \postcode100012, \stateBeijing, \countryChina

Since the detection of the first gravitational wave (GW) event [1], we have opened new avenues for advancing our understanding of gravity [2, 3], cosmology [4], and astrophysics of compact stars [5, 6]. The final stage of binary black hole (BH) mergers, known as ringdown, describes the progression of remnants to the stable phase, offering a unique opportunity to study BH theories, such as the no-hair theorem [7], the BH area law [8], and properties of the event horizon [9]. A ringdown signal is usually modelled as the superposition of quasi-normal modes (QNMs) of the remnant BH, which can be further divided into spin-weighted spheroidal harmonics with angular indices $(\ell,m)$. Each set of $(\ell,m)$ contains a series of overtone modes, indexed by $n$ [10, 11]. The extraction and measurement of multiple QNMs provide unparalleled opportunities for studying BHs, known as the BH spectroscopy [12, 11].

The parameter inference of GW signals is mainly based on the Bayes’ theorem,

where in the second equality, we have introduced the terminologies in the Bayesian inference: $\pi(x)\equiv P(x)$ is called the prior, representing the knowledge of the model parameters $x$ before observing the data $d$; the likelihood ${\mathcal{L}}(x)\equiv P(d|x)$ is determined by the model, encoding the statistical relation between the parameter and the data; the denominator ${\mathcal{Z}}\equiv P(d)$ is called evidence, which can also be regarded as the normalization factor of the numerator according to the law of total probability, ${\mathcal{Z}}=\int\pi(x){\mathcal{L}}(x){\rm{d}}x$. The posterior distribution $P(x|d)$ is the conditional probability of parameters given the data, representing the measurement of parameters after the observation. Therefore, given the prior and the likelihood, the Bayesian inference calculates the posterior distribution according to Eq. (1). However, obtaining an analytical expression (or an efficient numerical formula) of the posterior is usually unfeasible, and a more practical approach is to sample from the posterior distribution and estimate the evidence as an alternative. There are already some techniques to carry out this, such as the Markov-Chain Monte Carlo (MC) [13, 14] and nested sampling methods [15, 16]. These methods partially reduce the difficulty of posterior calculations, but the computational cost is still high when sampling in high-dimensional parameter spaces.

Currently, a typical Bayesian analysis of GW signals can take several hours to days, depending on the number of parameters and the adopted waveform templates. In the future, the next-generation (XG) ground-based GW detectors, such as the Cosmic Explorer (CE) [17, 18] and the Einstein Telescope (ET) [19, 20], will be observing with one order of magnitude higher sensitivity than the current ones, while the near-future space-based detectors, such as the Laser Interferometer Space Antenna (LISA) [21], Taiji [22] and TianQin [23], are expected to detect GWs in the millihertz band. It will become feasible to extract multiple QNMs from ringdown signals [24, 25, 26]. Each QNM component is described by a damping harmonic oscillator consisting of four parameters: amplitude, phase, damping time, and oscillation frequency [27]. In the General Relativity, the latter two parameters are determined by the final mass and spin of the remnant BH. Therefore, adding one QNM mode to ringdown signals increases the dimension of the parameter space by two, and extracting multiple QNMs leads to a challenge in the data analysis. The computational costs are even more prohibitive in tests of no-hair theorem with multiple QNMs, where each mode introduces four free parameters [24, 25, 26].

In this work, we propose an algorithm, named FIREFLY ($\mathcal{F}$-statistic Inspired REsampling For anaLYzing GW ringdown signals), to accelerate the Bayesian analysis of parameter estimation. This algorithm is inspired by the $\mathcal{F}$-statistic method in the continuous GW searches, which reduces the computational cost by maximizing the likelihood over the extrinsic parameters [28]. In the FIREFLY, the amplitude and phase parameters in the QNMs are analytically marginalized under a specifically chosen prior, reducing the dimensionality of parameter space [29, 30]. Using this auxiliary step, the inference under the target prior can be efficiently performed by importance sampling. In our work, the FIREFLY algorithm reduces the computational time from several hours to a few minutes, while producing consistent posterior distribution and evidence as in the traditional full-parameter Bayesian inference.

FIREFLY does not refine the sampling technique itself, instead, it utilizes the special form of the likelihood and designs more efficient sampling strategies to accelerate the posterior calculation. Compared to the pure $\mathcal{F}$-statistic-based methods, which implicitly assume some specific priors [31, 32, 29, 30], FIREFLY is flexible in prior choices because of its efficient importance-sampling design. Compared to current machine-learning-based accelerating algorithms for Bayesian GW data analysis [33, 34, 35, 36], FIREFLY only utilizes the marginalization and resampling techniques in Bayesian analysis, offering a clear and intuitive statistical interpretability. Moreover, any future improvements in the stochastic sampling technique can be incorporated into FIREFLY, further accelerating the inference.

Now we explain the logic of FIREFLY. The recorded GW strain of a ringdown signal is expressed as [27, 30],

where $B^{\ell mn,1}=A_{\ell mn}\cos\phi_{\ell mn}$ and $B^{\ell mn,2}=A_{\ell mn}\sin\phi_{\ell mn}$ are reparameterization of the QNM amplitude $A_{\ell mn}$ and phase $\phi_{\ell mn}$, while $g_{\ell mn,j}(t)$ depends on other source parameters, denoted as ${\bm{\theta}}$. Assuming a stationary and Gaussian noise, one finds that the likelihood has a Gaussian form with respect to ${\bm{B}}$ [30]

where $B^{\mu}$ represents the unified notation for $B^{\ell mn,j}$, $\langle d|d\rangle$ is noise-weighted inner product of the observed data $d$, while ${\mathcal{F}}$, $\hat{B}^{\mu}$, and $M_{\mu\nu}$ are functions of ${\bm{\theta}}$ but independent of ${\bm{B}}$ (see Methods for more details). Based on this Gaussian-form likelihood, FIREFLY aims to accelerate the Bayesian inference, by computing the posterior distribution and evidence given the likelihood (3) and the prior $\pi({\bm{\theta}},{\bm{B}}|{\otimes})$ in a target inference problem.

The FIREFLY algorithm mainly composes of two steps: auxiliary inference and importance sampling. In the auxiliary inference, the prior for ${\bm{\theta}}$ is chosen as the same as the target inference, denoted as $\pi({\bm{\theta}})$, while the prior for ${\bm{B}}$ is chosen to be independent of ${\bm{\theta}}$ and flat in ${\bm{B}}$ with a large-enough range, $\pi({\bm{B}}|{\bm{\theta}},{\varnothing})=\pi_{B}$. Note that we use ${\otimes}$ and ${\varnothing}$ to distinguish the target and auxiliary priors. Under the auxiliary prior, the QNM parameters can be analytically marginalized, leaving an inference problem only involving ${\bm{\theta}}$. The marginal likelihood for ${\bm{\theta}}$ reads

where $N$ is the number of QNMs considered in the ringdown waveform. FIREFLY firstly draws samples of ${\bm{\theta}}$ with sampling techniques (such as the nested sampling) based on the marginal likelihood and the prior for ${\bm{\theta}}$. These samples are denoted by $\{{\bm{\theta}}_{i}|\varnothing\}$. The evidence ${\mathcal{Z}}_{\varnothing}$ is also obtained by the adopted sampling technique. The auxiliary inference is completed by resampling ${\bm{B}}$, which only needs to draw one sample ${\bm{B}}_{i}$ for the $i$-th sample ${\bm{\theta}}_{i}$ from the Gaussian distribution, ${\cal N}\big{(}\hat{{\bm{B}}}({\bm{\theta}}_{i}),M^{-1}({\bm{\theta}}_{i})\big{)}$; see Methods for details.

In the second step, FIREFLY performs importance sampling based on the posterior samples, $\{{\bm{\theta}}_{i},{\bm{B}}_{i}|{\varnothing}\}$, and the evidence ${\mathcal{Z}}_{\varnothing}$ produced by the auxiliary inference. The evidence under the target prior can be estimated by a MC integration, while for the posterior, a direct importance sampling weighted by the prior ratio $\pi({\bm{\theta}},{\bm{B}}|{\otimes})/\pi({\bm{\theta}},{\bm{B}}|{\varnothing})$ may lead to large variance [37], especially when the target and auxiliary priors are significantly different (see Methods). Therefore, FIREFLY adopts a two-step importance sampling, dealing with ${\bm{\theta}}$ and ${\bm{B}}$ separately. First, the sample weights for $\{{\bm{\theta}}_{i}|{\varnothing}\}$ are calculated according to the marginal-posterior ratio for ${\bm{\theta}}$ between the target and auxiliary inferences. After the importance resampling for ${\bm{\theta}}$, the sample ${\bm{B}}_{i}$, assigned to each ${\bm{\theta}}_{i}$, is drawn from the conditional posterior,

where the importance-sampling technique is applied once again with a Gaussian proposer.

The output of FIREFLY is the posterior samples $\{{\bm{\theta}}_{i},{\bm{B}}_{i}|{\otimes}\}$ and the evidence ${\mathcal{Z}}_{\otimes}$ under the target prior, which is the same as in the full Bayesian inference with all parameters. As we can see, FIREFLY takes advantage of the Gaussian-form likelihood in ${\bm{B}}$, and eliminates $2N$ parameters analytically before applying the time-consuming stochastic sampling techniques. We summarize the workflow of FIREFLY in Fig. 1.

Now we apply FIREFLY to analyze realistic ringdown signals, and compare the results with those from the full Bayesian inference with all parameters, which we refer as the “full-parameter sampling”. To test the performance of FIREFLY in accelerating the extraction of multiple QNMs, we assume the signals are recorded in the triangular-shape ET which consists of three detectors [20]. The data are generated from Numerical Relativity (NR) waveforms of Simulating eXtreme Spacetimes (SXS) catalog [38] and injected into a Gaussian stationary noise corresponding to the ET-D sensitivity.

Specifically, we inject the NR signal SXS:BBH:0305, which is generated by a GW150914-like, non-precessing binary BH (BBH) with a mass ratio of $1.2$, a detector-frame final mass of $68.2\,{M_{\odot}}$, and a dimensionless spin of $0.69$ [38]. The luminosity distance, inclination angle, and reference phase are set to $390\,{\rm Mpc}$, $3\uppi/4$, and $0$, respectively. The signal-to-noise ratio of ringdown is approximately $312$ in the ET, assuming that the signal’s starting time is at its peak amplitude [30]. Usually, for ringdown signals from nearly equal-mass BBHs, one focuses solely on the extraction of QNMs with $l=|m|=2$, while considering different number of overtones [39]. In this work, we consider three scenarios: (i) only the fundamental mode ($N=1$), (ii) one additional overtone mode ($N=2$), and (iii) two additional overtone modes ($N=3$). Since higher-order overtones decay faster and are only potentially significant at earlier times [39], we use the ringdown part from the whole signal at starting times of $18\,\rm{M}$, $23\,\rm{M}$ and $28\,\rm{M}$ (with the final mass $\rm{M}=68.2\,{M_{\odot}}$) after coalescence for $N=1$, $2$ and $3$, respectively.¹¹1Here we adopt the geometric units $G=c=1$, and $68.2\,{M}_{\odot}$ corresponds to $0.33\,{\rm ms}$. In each scenario, it was verified that reasonable adjustments of the starting time do not significantly change our results [30].

When extracting QNMs from the ringdown signal, the variable parameters in the waveform include the right ascension $\alpha$, declination $\delta$, polarization angle $\psi$, inclination angle $\iota$, coalescence time $t_{c}$, final mass $M_{f}$, and final dimensionless spin $\chi_{f}$, in addition to the QNM parameters ${\bm{B}}$. It is also customary to fix the first five parameters based on a more comprehensive study, such as the inspiral-merger-ringdown analysis. For simplicity, we adopt the sky-averaged antenna pattern functions over $\alpha$, $\delta$, and $\psi$, and fix $\iota$ and $t_{c}$ to their injected values. Therefore, in this work ${\bm{\theta}}$ only contains the final mass $M_{f}$ and final spin $\chi_{f}$ of the remnant BH, and the full set of variables is $\big{\{}M_{f},\chi_{f},A_{\ell mn},\phi_{\ell mn}\big{\}}$, where $(\ell,m,n)$ runs over all QNMs considered in the ringdown signal. Notably, expanding the set of free parameters does not affect the applicability of FIREFLY. In the target inference, we choose flat priors for all parameters, $M_{f}/M_{\odot}\sim{\cal U}(50,100)$, $\chi_{f}\sim{\cal U}(0,0.99)$, $\phi_{\ell mn}\sim{\cal U}(0,2\uppi)$ and $A_{\ell mn}\sim{\cal U}(0,A_{\max})$ with $A_{\max}=5\times 10^{-20}$. In the auxiliary inference, the prior for $M_{f}$ and $\chi_{f}$ remains the same, while the flat prior in ${\bm{B}}$ corresponds to $\pi\big{(}A_{\ell mn},\phi_{\ell mn}\big{|}{\bm{\theta}},{\varnothing}\big{)}\propto A_{\ell mn}$. We adopt this prior shape for all $A_{\ell mn}$ and $\phi_{\ell mn}$, and normalize it in the region of $A_{\ell mn}\leq A_{\max}$ and $0\leq\phi_{\ell mn}<2\uppi$.

Now we compare the posteriors obtained by FIREFLY and the full-parameter sampling, which directly samples all parameters $\{{\bm{\theta}},{\bm{B}}\}$ under the target prior. In both the auxiliary inference and the full-parameter sampling, we use the dynesty [40] sampler to draw samples from the posterior. In all three scenarios, we find that the FIREFLY results closely agree with those from the full-parameter sampling. For concreteness, we show the $N=3$ case in Fig. 2, where the green and blue contours represent the results from the FIREFLY and the full-parameter sampling, respectively. The P-P plots between the one-dimensional marginalized posteriors of the two methods are shown in Fig. 3. We also use the Wasserstein distance to quantitatively verify the consistency between the two methods [41], and find that the average difference of one-dimensional marginal distributions is less than 0.1 standard deviation. The computational times of FIREFLY and the full-parameter sampling are listed in Table 1. In the $N=3$ case, FIREFLY reduces the computational time from several hours to just a few minutes, speeding up by two orders of magnitude.

The main reason why FIREFLY significantly accelerates the inference is the analytical marginalization of the QNM amplitude and phase parameters in the auxiliary inference, which eliminates $2N$ free parameters in the stochastic sampling process. The computational advantage of FIREFLY becomes more pronounced when involving more QNMs in the analysis. Due to the same reason, FIREFLY is fully statistically interpretable, and also compatible with various sampling techniques. Since the majority of the computational time in FIREFLY is spent on the auxiliary inference (see Table 1), one can incorporate the state-of-the-art sampling techniques into FIREFLY, further reducing the computational cost.

FIREFLY also allows flexible prior choices, which is based on the effective two-step importance sampling. In the ringdown signal analysis, different prior choices can significantly affect posteriors, especially when considering multiple QNMs or analyzing weak signals (see Methods). In contrast, the pure $\mathcal{F}$-statistic-based methods reduce the number of parameters by analytical maximization, which inherently assumes some specific priors [31, 32, 29, 30]. The flexibility of priors in FIREFLY makes it possible to conduct a more comprehensive analysis when extracting multiple QNMs from ringdown signals. It is also worth noting that changing to a new target prior only requires one to rerun an importance sampling that takes only a few seconds, while the results of auxiliary inference can be reused.

Because of its structure, FIREFLY also gives the evidence estimates, and we compare them with those from the full-parameter sampling in Fig. 4. FIREFLY estimates the evidence with a smaller variance than in the full-parameter sampling. This is because that the uncertainties of the evidence estimation mainly come from the stochastic sampling of the auxiliary inference in FIREFLY, which has $2N$ fewer free parameters than in the full-parameter sampling. Nevertheless, the evidence estimated by FIREFLY tends to be slightly larger than that of the full-parameter sampling. Similar discrepancies have been reported by other approaches designed to accelerate evidence calculation in GW Bayesian inference [35]. To address this issue, we conduct a comprehensive investigation into the potential uncertainty source and confirm the consistency of the results produced by two methods after considering the additional uncertainties from the insufficient sampling in the full-parameter sampling (see Appendix). Having said that, the relative differences between the two methods are at the level of $10^{-3}$, so the evidence given by FIREFLY is still largely consistent with the full-parameter sampling.

After validating FIREFLY in the parameter estimation of ringdown signals, we now present the extraction of higher modes in the ringdown signal from an asymmetric-mass-ratio BBH system. We inject the NR signal SXS:BBH:0065, which is generated from a non-precessing binary with a mass ratio of $8.0$, a final redshifted mass of $122.45\,{M_{\odot}}$, and a final spin of $0.66$ [38]. The luminosity distance, inclination angle, and reference phase are set to $390\,{\rm Mpc}$, $3\uppi/4$, and $0$, respectively. We extract four QNMs $(\ell,m,n)=(2,2,1)$, $(2,1,1)$, $(3,3,1)$, and $(4,4,1)$ from the ringdown signal, at different start times $\Delta t=16\,{\rm M}$, $20\,{\rm M}$, $24\,{\rm M}$, and $28\,{\rm M}$ with ${\rm M}=122.45\,{M_{\odot}}$.²²2In the geometric units, ${\rm M}=122.45\,{M_{\odot}}$ corresponds to $0.60\,{\rm ms}$. The prior for the final mass is chosen to be $M_{f}/M_{\odot}\sim{\cal U}(100,150)$, while priors of other parameters are the same as those in the previous overtone analysis. The joint posteriors of the final mass and the final spin from FIREFLY and the full-parameter sampling are shown in Fig. 5. FIREFLY achieves close agreement with the full-parameter sampling, while reducing the computational time by about sixty times (see the last row in Table 1).

In this work, we develop the FIREFLY algorithm, which significantly accelerates the analysis of ringdown signals while maintaining nearly indistinguishable posterior and evidence estimates. In the spirit of ${\cal F}$-statistic, by constructing an auxiliary inference where the QNM amplitudes and phases are analytically marginalized, FIREFLY reduces the computational time by several orders of magnitude, demonstrating its potential to address the challenge of dimensionality curse when extracting multiple QNMs. Also, based on the delicately designed importance-sampling strategy, FIREFLY is highly flexible in prior choices, allowing for a more comprehensive analysis for identifying overtones and higher modes in the ringdown signals. FIREFLY presents a computational paradigm that only leverages the intrinsic structure of the likelihood to accelerate the Bayesian inference, which is fully statistically interpretable, and compatible with any improvements in stochastic sampling techniques. Similar structures are commonly encountered in GW analyses, such as for extreme-mass-ratio inspirals (EMRIs) [42] and GW localization [43]. We hope that the novelty of the FIREFLY algorithm will provide a viable solution to address many more data challenges in future GW data analysis and the beyond.