MatBYIB: A Matlab-based code for Bayesian inference of extreme mass-ratio inspiral binary with arbitrary eccentricity

TL;DR

MatBYIB tackles the challenge of Bayesian parameter estimation for gravitational waves from EMRIs with arbitrary eccentricity by marrying Computationally efficient Analytical Kludge waveforms with both Fisher Information Matrix forecasts and full Metropolis–Hastings MCMC posterior sampling. The framework enables rapid, near-analytic error estimates and robust, full posterior exploration, using Gelman–Rubin convergence across parallel chains to guarantee sampling adequacy. Its modular MATLAB implementation (waveform generation, detector response, FIM, and MCMC) is open-source and demonstrated on representative EMRI-like cases, showing strong agreement between FIM predictions and MCMC posteriors, with MatBYIB achieving convergence on standard hardware. This approach provides a practical, scalable tool for parameter estimation of eccentric EMRIs for current and future space-based detectors such as LISA and Taiji, and lays groundwork for incorporating more physics and higher-order waveform models.

Abstract

Accurate parameter estimation(PE) of gravitational waves(GW) is essential for GW data analysis. In extreme mass-ratio inspiral binary(EMRI) systems, orbital eccentricity is a critical parameter for PE. However, current software for for PE of GW often neglects the direct estimation of orbital eccentricity. To fill this gap, we have developed the MatBYIB, a MATLAB-based software package for PE of GW with arbitrary eccentricity. The MatBYIB employs the Analytical Kludge (AK) waveform as a computationally efficient signal generator and computes parameter uncertainties via the Fisher Information Matrix (FIM) and the Markov Chain Monte Carlo (MCMC). For Bayesian inference, we implement the Metropolis-Hastings (M-H) algorithm to derive posterior distributions. To guarantee convergence, the Gelman-Rubin convergence criterion (the Potential Scale Reduction Factor R) is used to determine sampling adequacy, with MatBYIB dynamically increasing the sample size until R < 1.05 for all parameters. Our results demonstrate strong agreement between FIM- based predictions and full MCMC sampling. This program is user-friendly and allows for estimation of gravitational wave parameters with arbitrary eccentricity on standard personal computers. Code availability:The implementation is open-source at https://github.com/GenliangLi/MatBYIB.

Figures (6)

• Figure 1: Schematic Diagram of the Software Architecture
• Figure 2: The evolution of orbital frequency (Left) and eccentricity (Right) with time. The system parameters are as follows: the mass of the central BH is $M = 10^6 M_{\odot}$, the mass of the orbiting object is $10 M_{\odot}$, the spin parameter of the central BH is $S/M^2 = 0.4$, the eccentricity at the last stable orbit $e_\mathrm{LSO} = 0.0,0.3,0.6$, and the frequency at the last stable orbit is calculated by ${f_\mathrm{LSO}}=c^{3} /(2 \pi G M)\left(\left(1-e_{\mathrm{LSO}}^{2}\right) /\left(6+2 e_{\mathrm{LSO}}\right)\right)^{3 / 2}$. The evolution time is $t_c = 3.14 \times 10^6$ seconds.
• Figure 3: The time-domain gravitational waveforms h(t) are shown for three distinct time intervals: ${0 \sim 3000}$ (left panel), ${1.109\times10^6 \sim 1.204\times10^6}$ (middle panel) and ${3.137\times10^6 \sim 3.140\times10^6}$ (right panel). The system parameters are configured as follows: CO's mass: ${m_2=10 M_{\odot}}$; MBH's mass: ${M=10^{6} M_{\odot}}$; MBH's spin magnitude: ${S=0.01 M^{2}}$; Angle between MBH's spin and orbital angular momentum: ${\lambda=60{^\circ}}$; We set $\phi_\mathrm{LSO} = \gamma_\mathrm{LSO} = \alpha_\mathrm{LSO} = 0$; Sky angle $\theta_S = \phi_S = \theta_k =\phi_k = 60 ^\circ$.
• Figure 4: GW signal on frequency domain. The frequency-domain GW waveforms by FFT on time-domain GW in Fig. \ref{['fig:ht']} (Left); The characteristic strain of GW signal by $h_c = 2f|h(f)|$ (Right) . The parameters are the same as in Fig. \ref{['fig:ht']}.
• Figure 5: The corner plot from FIM exploration with fiducial/injected values $m_1=10^{6} M_{\odot}$ (central mass), $m_2=10 M_{\odot}$ (orbiting mass), $S/M^{2}=0.01$ (dimensionless spin of the central BH), $e_\mathrm{LSO}=0.3$. We have assumed an observation of $t_c = 3.14\times10^6 s$. redshift $z=0.01$, Median and $68 \%$ confidence interval are $m_1 =1_{-1.09\times10^{+2}}^{+1.08\times10^{+2}} \times 10^{6} M_{\odot}$, $m_2=10_{-1.70\times10^{-3}}^{+1.69\times10^{-3}} M_{\odot}$, $e_\mathrm{LSO} =0.3_{-7.39\times10^{-5}}^{+7.49\times{10^{-5}}}$, and $S/M^{2}=0.01_{-2.49\times10^{-4}}^{+2.47\times10^{-4}}$.