Persistence of post-Newtonian amplitude structure in binary black hole mergers

TL;DR

This study examines the persistence of post-Newtonian amplitude structures in binary black hole mergers, using numerical relativity data to develop polynomial fits that model mode amplitudes from inspiral to postmerger.

Contribution

It introduces polynomial fitting methods that extend post-Newtonian amplitude models through merger, improving waveform modeling accuracy in strong gravitational fields.

Findings

PN dependence persists in certain modes throughout merger.

Polynomial fits of degree ≤ 3 effectively model amplitude behavior near merger.

Discrepancies across catalogs highlight resolution effects.

Abstract

We analyze the spherical harmonic mode amplitudes of quasicircular, nonprecessing binary black hole mergers using 275 numerical relativity simulations from the SXS, RIT, and MAYA catalogs. We construct fits using the leading-order post-Newtonian (PN) dependence on intrinsic parameters, replacing the PN velocity with fit coefficients. We compare these to polynomial fits in symmetric mass ratio and spin. We analyze $(\ell, m)$ modes with $\ell \leq 4$ from late inspiral [$t = -500M$ relative to the $(2,2)$ peak] to postmerger ($t = 40M$). For nonspinning systems, the $(2,2)$, $(2,1)$, and $(3,3)$ modes retain the leading-order PN dependence on mass ratio throughout the merger. Higher-order modes deviate from the PN dependence only near and after the merger, where polynomial fits of degree $N \leq 3$ can capture the amplitude behavior up to $40M$. For aligned-spin systems at fixed mass ratio, the $(2,1)$ mode retains its PN spin dependence, while the $(3,2)$ and $(4,3)$ modes exhibit a quadratic spin dependence near merger. The PN-inspired fits lose accuracy with increasing mass ratio, particularly near merger. Results broadly agree across catalogs, though discrepancies appear in the $(3,1)$, $(4,2)$, and $(4,1)$ modes, likely from resolution differences. Our results clarify the extent to which PN structure persists in mode amplitudes. Although the fits cannot be fully interpreted within the PN formalism near merger, low-degree polynomial corrections to the PN amplitude Ans\"atze can capture strong-field behavior, enabling closed-form and efficient modeling of waveform amplitudes in this regime.

Tables

Table 1. Table 1: Optimal fit coefficients for the fits in Eqs. ( 7 ) at t = − 8000 ​ M t=-8000M .

Mode $(\ell ,m)$	$a_{\ell m}\left(t=-8000M\right)$
$(2,2)$	$0.23017$
$(2,1)$	$0.27610$
$(3,3)$	$0.24277$
$(4,4)$	$0.24192$

Table 2. Table 2: List of 103 nonspinning simulations from the SXS , RIT , and MAYA catalogs used in this work, grouped by mass ratio q q .

$q$	ID(s)		$q$	ID(s)		$q$	ID(s)
$1.00$	SXS:BBH:0389	SXS:BBH:2377	$2.60$	GT0651		$7.00$	SXS:BBH:2491	RIT:BBH:0416
	SXS:BBH:2378	SXS:BBH:3624
	SXS:BBH:4434	SXS:BBH:2496
	SXS:BBH:1132	SXS:BBH:3864
	SXS:BBH:1154	SXS:BBH:1155
	SXS:BBH:1153	SXS:BBH:3617
	SXS:BBH:2375	SXS:BBH:3634
	SXS:BBH:2376	SXS:BBH:2325
	SXS:BBH:2326	RIT:BBH:0112
	RIT:eBBH:1090
$1.08$	RIT:BBH:0458		$2.70$	GT0650		$7.19$	SXS:BBH:0188
$1.18$	RIT:BBH:0084	RIT:BBH:0113	$2.80$	GT0653		$7.50$	SXS:BBH:2490
$1.20$	SXS:BBH:0198		$2.90$	GT0654		$7.76$	SXS:BBH:0195
$1.33$	RIT:BBH:0045	RIT:BBH:0114	$3.00$	SXS:BBH:1179	SXS:BBH:2498	$8.00$	SXS:BBH:2707	SXS:BBH:3622
	RIT:eBBH:1241			SXS:BBH:2265	RIT:BBH:0102
				RIT:BBH:0138
$1.50$	SXS:BBH:2331	SXS:BBH:3984	$3.10$	GT0758		$8.27$	SXS:BBH:0186
	SXS:BBH:0593	RIT:BBH:0085
	RIT:BBH:0115
$1.67$	RIT:BBH:0086	RIT:BBH:0116	$3.50$	SXS:BBH:2483		$8.73$	SXS:BBH:0199
$1.75$	GT0727		$4.00$	SXS:BBH:2485	SXS:BBH:3631	$9.00$	SXS:BBH:2495
				SXS:BBH:2499	SXS:BBH:1220
				RIT:BBH:0088	RIT:BBH:0119
				RIT:eBBH:1133	GT0454
				MAYA1040
$1.82$	RIT:BBH:1020		$4.50$	SXS:BBH:2484	SXS:BBH:3144	$9.17$	SXS:BBH:0189
$2.00$	SXS:BBH:2497	SXS:BBH:1166	$5.00$	SXS:BBH:2374	SXS:BBH:2487	$9.20$	SXS:BBH:1108
	SXS:BBH:1167	RIT:BBH:0077		SXS:BBH:3619	RIT:BBH:0089
	RIT:BBH:0117	RIT:eBBH:1200		RIT:BBH:0120	RIT:BBH:0152
	GT0446
$2.25$	GT0757		$5.50$	SXS:BBH:2486		$9.99$	SXS:BBH:0185
$2.32$	SXS:BBH:0201		$6.00$	SXS:BBH:2164	SXS:BBH:2489	$10.00$	RIT:BBH:0978
				SXS:BBH:3630	RIT:BBH:0090
				RIT:BBH:0121
$2.41$	RIT:BBH:0139		$6.50$	SXS:BBH:2488		$14.00$	SXS:BBH:2480
$2.50$	RIT:BBH:0087	RIT:BBH:0118	$6.58$	SXS:BBH:0192		$15.00$	SXS:BBH:2477	RIT:BBH:0373
							RIT:BBH:0942	RIT:BBH:0943
							RIT:BBH:0957

Table 3. Table 3: List of 172 aligned-spin simulations from the SXS , RIT , and MAYA catalogs used in this work. Only SXS simulations were used to construct the fits in Figs. 10 – 12 , but RIT and MAYA simulations are included for visual comparison.

$q$	ID	${\chi }_{1z}$	${\chi }_{2z}$	$q$	ID	${\chi }_{1z}$	${\chi }_{2z}$	$q$	ID	${\chi }_{1z}$	${\chi }_{2z}$
1.0	SXS:BBH:2100	0.5	-0.9	1.0	SXS:BBH:1497	0.68	0.67	2.0	SXS:BBH:2131	0.85	0.85
1.0	SXS:BBH:2093	-0.5	0.9	1.0	SXS:BBH:3926	0.46	-0.32	2.0	SXS:BBH:2132	0.87	-0.85
1.0	SXS:BBH:1492	-0.47	-0.79	1.0	SXS:BBH:0330	-0.8	0.8	2.0	SXS:BBH:2109	-0.6	-0.6
1.0	SXS:BBH:2099	0.4	0.8	1.0	SXS:BBH:2085	-0.9	0.9	2.0	SXS:BBH:2127	0.5	0.5
1.0	SXS:BBH:2098	0.4	-0.8	1.0	SXS:BBH:2091	-0.6	0.6	2.0	SXS:BBH:2114	-0.3	-0.3
1.0	SXS:BBH:2095	-0.4	0.8	1.0	SXS:BBH:1476	-0.8	0.8	2.0	SXS:BBH:2116	-0.3	0.3
1.0	SXS:BBH:2094	-0.4	-0.8	1.0	SXS:BBH:2087	-0.8	0.8	2.0	SXS:BBH:2123	0.3	-0.3
1.0	SXS:BBH:3912	-0.47	-0.79	1.0	SXS:BBH:3927	0.51	0.29	2.0	SXS:BBH:0448	0.4	-0.4
1.0	SXS:BBH:3982	-0.18	0.72	1.0	SXS:BBH:1507	0.51	0.29	2.0	SXS:BBH:2130	0.6	0.6
1.0	SXS:BBH:3923	0.14	0.73	1.0	SXS:BBH:0376	0.6	-0.4	2.0	SXS:BBH:0554	0.2	-0.4
1.0	SXS:BBH:1502	-0.42	0.7	1.0	SXS:BBH:0394	0.6	0.4	2.0	SXS:BBH:2111	-0.6	0.6
1.0	SXS:BBH:3922	-0.42	0.7	1.0	SXS:BBH:0462	-0.6	-0.4	2.0	SXS:BBH:2112	-0.5	-0.5
1.0	SXS:BBH:4066	-0.95	0.95	1.0	SXS:BBH:0447	-0.6	0.4	2.0	SXS:BBH:2510	0.6	0.6
1.0	SXS:BBH:3977	0.95	-0.95	1.0	SXS:BBH:2088	-0.62	-0.25	2.0	SXS:BBH:2506	0.5	0.5
1.0	SXS:BBH:0436	-0.2	-0.4	1.0	SXS:BBH:2103	0.65	0.25	2.0	SXS:BBH:2501	0.25	0.25
1.0	SXS:BBH:3534	-0.8	0.8	1.0	SXS:BBH:3919	-0.75	0.34	2.0	SXS:BBH:0513	0.6	-0.4
1.0	SXS:BBH:0304	0.5	-0.5	1.0	SXS:BBH:3915	0.78	0.53	2.0	SXS:BBH:0387	-0.6	-0.4
1.0	SXS:BBH:0459	0.2	-0.4	1.0	SXS:BBH:3920	-0.77	-0.2	2.0	SXS:BBH:0614	0.75	-0.5
1.0	SXS:BBH:0370	-0.2	0.4	1.0	SXS:BBH:1495	0.78	0.53	2.0	SXS:BBH:2125	0.3	0.3
1.0	SXS:BBH:4029	0.95	-0.95	1.0	SXS:BBH:3553	0.8	0.4	2.0	SXS:BBH:0412	-0.2	-0.4
1.0	SXS:BBH:0327	0.8	-0.8	1.0	SXS:BBH:3582	-0.8	-0.4	2.0	SXS:BBH:3981	-0.12	-0.65
1.0	SXS:BBH:2092	-0.5	0.5	1.0	SXS:BBH:3578	-0.8	0.4	1.0	RIT:BBH:0016	0.8	-0.8
1.0	SXS:BBH:2422	0.6	0.6	1.0	SXS:BBH:1499	-0.75	0.34	1.0	RIT:BBH:0228	-0.4	0.85
1.0	SXS:BBH:2419	-0.2	-0.2	1.0	SXS:BBH:3601	0.8	-0.4	1.0	RIT:BBH:0261	0.4	-0.85
1.0	SXS:BBH:2418	-0.44	-0.44	1.0	SXS:BBH:1500	-0.77	-0.2	1.0	RIT:BBH:0324	0.9	0.9
1.0	SXS:BBH:4072	-0.95	-0.95	1.0	SXS:BBH:2106	0.9	0.5	1.0	RIT:BBH:0499	0.95	-0.95
1.0	SXS:BBH:3629	-0.9	-0.9	1.0	SXS:BBH:2083	-0.9	-0.5	1.0	RIT:BBH:0573	0.95	0.0
1.0	SXS:BBH:1475	-0.8	-0.8	1.5	SXS:BBH:0437	-0.2	-0.8	1.5	RIT:BBH:0289	0.8	-0.8
1.0	SXS:BBH:2102	0.6	0.6	1.5	SXS:BBH:0579	0.4	-0.8	1.5	RIT:BBH:0616	0.95	-0.95
1.0	SXS:BBH:2500	0.25	0.25	1.5	SXS:BBH:2342	-0.5	0.5	1.5	RIT:BBH:0687	0.95	0.95
1.0	SXS:BBH:0160	0.9	0.9	1.5	SXS:BBH:0369	0.6	-0.8	1.5	RIT:BBH:0788	0.95	0.0
1.0	SXS:BBH:1122	0.44	0.44	1.5	SXS:BBH:0392	-0.2	0.8	2.0	RIT:BBH:0083	0.5	-0.5
1.0	SXS:BBH:1123	0.5	0.5	1.5	SXS:BBH:2348	0.5	-0.5	2.0	RIT:BBH:0245	-0.85	0.5
1.0	SXS:BBH:2505	0.5	0.5	1.5	SXS:BBH:0441	0.6	0.8	2.0	RIT:BBH:0247	0.85	-0.5
1.0	SXS:BBH:2089	-0.6	-0.6	1.5	SXS:BBH:0397	-0.8	-0.4	2.0	RIT:BBH:0296	-0.85	-0.5
1.0	SXS:BBH:2512	0.7	0.7	1.5	SXS:BBH:1415	0.5	0.5	2.0	RIT:BBH:0314	0.85	0.5
1.0	SXS:BBH:0178	0.99	0.99	1.5	SXS:BBH:0372	0.8	-0.4	2.0	RIT:BBH:0336	0.0	-0.8
1.0	SXS:BBH:3897	0.8	0.8	2.0	SXS:BBH:0619	0.9	0.9	2.0	RIT:BBH:0338	0.0	0.8
1.0	SXS:BBH:0176	0.96	0.96	2.0	SXS:BBH:2113	-0.37	0.85	2.0	RIT:BBH:0446	0.5	-0.8
1.0	SXS:BBH:3895	-0.8	-0.8	2.0	SXS:BBH:2126	0.37	-0.85	2.0	RIT:BBH:0447	-0.25	-0.25
1.0	SXS:BBH:3627	0.95	0.95	2.0	SXS:BBH:2122	0.13	0.85	2.0	RIT:BBH:0451	-0.5	0.8
1.0	SXS:BBH:2086	-0.8	-0.8	2.0	SXS:BBH:0617	0.5	0.75	2.0	RIT:BBH:0452	0.25	-0.25
1.0	SXS:BBH:2509	0.6	0.6	2.0	SXS:BBH:2117	-0.13	-0.85	2.0	RIT:BBH:0454	-0.25	0.25
1.0	SXS:BBH:3628	0.97	0.97	2.0	SXS:BBH:0399	0.2	0.4	2.0	RIT:BBH:0457	0.25	0.25
1.0	SXS:BBH:3518	-0.8	-0.8	2.0	SXS:BBH:0333	0.8	0.8	2.0	RIT:BBH:0462	0.85	-0.25
1.0	SXS:BBH:3625	-0.8	-0.8	2.0	SXS:BBH:0618	0.8	0.8	2.0	RIT:BBH:0465	-0.85	-0.25
1.0	SXS:BBH:3976	-0.95	-0.95	2.0	SXS:BBH:4073	0.9	0.9	2.0	RIT:BBH:0466	-0.85	0.25
1.0	SXS:BBH:2420	0.2	0.2	2.0	SXS:BBH:4078	0.9	-0.9	2.0	RIT:BBH:0477	0.85	0.25
1.0	SXS:BBH:2421	-0.6	-0.6	2.0	SXS:BBH:4120	-0.9	-0.9	2.0	RIT:BBH:0679	0.95	-0.95
1.0	SXS:BBH:3896	-0.8	0.8	2.0	SXS:BBH:0584	-0.4	-0.4	2.0	RIT:BBH:0767	0.95	0.95
1.0	SXS:BBH:2423	0.85	0.85	2.0	SXS:BBH:4115	-0.9	0.9	2.0	RIT:BBH:0782	0.0	0.95
1.0	SXS:BBH:0172	0.98	0.98	2.0	SXS:BBH:2107	-0.87	0.85	2.0	RIT:BBH:0796	0.95	0.0
1.0	SXS:BBH:3978	0.95	0.95	2.0	SXS:BBH:0354	-0.2	0.4	1.0	GT0420	0.2	0.2
1.0	SXS:BBH:0329	-0.8	-0.8	2.0	SXS:BBH:2108	-0.85	-0.85	1.0	GT0546	0.0	0.2
1.0	SXS:BBH:3929	-0.24	-0.1	2.0	SXS:BBH:0334	-0.8	-0.8	1.0	GT0574	0.1	0.1
1.0	SXS:BBH:1509	-0.24	-0.1	2.0	SXS:BBH:0335	-0.8	0.8	1.5	GT0479	0.2	0.2
1.0	SXS:BBH:1506	0.46	-0.32	2.0	SXS:BBH:2128	0.6	-0.6
1.0	SXS:BBH:3917	0.68	0.67	2.0	SXS:BBH:0574	0.4	0.4

Equations

$$h(t;\theta ,\varphi )=\sum _{\ell =2}^{\infty }\sum _{m=-\ell }^{m=\ell }{h}_{\ell m}(t)Y_{\ell m}^{-2}(\theta ,\varphi ),$$

$$h_{\ell m}(t)=A_{\ell m}(t)e^{i{\varphi }_{\ell m}(t)},$$

$${\hat{A}}_{\ell m}(t)=\{\begin{array}{ll}{A}_{\ell m}(t)& \text{if }(\ell ,m)=(2,2)\\ \frac{A_{\ell m}(t)}{A_{22}(t)}& \text{otherwise}.\end{array}$$

$$\eta \equiv \frac{m_1{m}_2}{(m_1+m_2{)}^2},$$

$$\delta \equiv \frac{m_1-m_2}{m_1+m_2}=\sqrt{1-4\eta }.$$

$${\chi }_s\equiv \frac{1}{2}\left({\chi }_1+{\chi }_2\right),$$

$${\chi }_a\equiv \frac{1}{2}\left({\chi }_1-{\chi }_2\right).$$

$${\hat{A}}_{22}^{(\mathrm{LO})}$$

$${\hat{A}}_{21}^{(\mathrm{LO})}$$

$${\hat{A}}_{33}^{(\mathrm{LO})}$$

$${\hat{A}}_{32}^{(\mathrm{LO})}$$

$${\hat{A}}_{31}^{(\mathrm{LO})}$$

$${\hat{A}}_{44}^{(\mathrm{LO})}$$

$${\hat{A}}_{43}^{(\mathrm{LO})}$$

$${\hat{A}}_{42}^{(\mathrm{LO})}$$

$${\hat{A}}_{41}^{(\mathrm{LO})}$$

$$8\sqrt{\frac{\pi }{5}}\eta {\left[v(t;\eta )\right]}^2+\dots$$

$${\hat{A}}_{22}^{(\mathrm{LO},\mathrm{S})}$$

$${\hat{A}}_{21}^{(\mathrm{LO},\mathrm{S})}$$

$${\hat{A}}_{32}^{(\mathrm{LO},\mathrm{S})}$$

$${\hat{A}}_{43}^{(\mathrm{LO},\mathrm{S})}$$

$${\hat{A}}_{41}^{(\mathrm{LO},\mathrm{S})}$$

$${\hat{A}}_{\ell m}^{(N)}$$

$${\hat{A}}_{\ell m}^{(1,\text{S})}$$

$${\hat{A}}_{\ell m}^{(2,\text{S})}$$

$$C_{\ell m}=\frac{({\hat{A}}_{\ell m}^{(\mathrm{fit})}-{\bar{A}}_{\ell m}^{(\mathrm{fit})})\cdot ({\hat{A}}_{\ell m}^{(\mathrm{NR})}-{\bar{A}}_{\ell m}^{(\mathrm{NR})})}{\sqrt{({\hat{A}}_{\ell m}^{(\mathrm{fit})}-{\bar{A}}_{\ell m}^{(\mathrm{fit})}{)}^2({\hat{A}}_{\ell m}^{(\mathrm{NR})}-{\bar{A}}_{\ell m}^{(\mathrm{NR})}{)}^2}},$$

$$p(\vec{c}|{\hat{A}}_{\ell m}^{(\mathrm{NR})})=\frac{\mathcal{L}({\hat{A}}_{\ell m}^{(\mathrm{NR})}|\vec{c})\pi (\vec{c})}{\mathcal{Z}({\hat{A}}_{\ell m}^{(\mathrm{NR})})},$$

$$\mathcal{Z}({\hat{A}}_{\ell m}^{(\mathrm{NR})})=\int \mathcal{L}({\hat{A}}_{\ell m}^{(\mathrm{NR})}|\vec{c})\pi (\vec{c})d\vec{c}$$

$$\log B{F}_{12}=\log {\mathcal{Z}}_1-\log {\mathcal{Z}}_2$$

$${\hat{A}}_{22}^{(\text{PN})}$$

$${\hat{A}}_{21}^{(\text{PN})}$$

$${\hat{A}}_{33}^{(\text{PN})}$$

$${\hat{A}}_{32}^{(\text{PN})}$$

$${\hat{A}}_{31}^{(\text{PN})}$$

$${\hat{A}}_{44}^{(\text{PN})}$$

$${\hat{A}}_{43}^{(\text{PN})}$$

$${\hat{A}}_{42}^{(\text{PN})}$$

$${\hat{A}}_{41}^{(\text{PN})}$$

$${\hat{A}}_{22}^{(\text{PN, S})}$$

$${\hat{A}}_{21}^{(\text{PN, S})}$$

$$+\frac{1}{6}{v}^6({\chi }_a+\delta {\chi }_s)\left(3\pi -4\delta {\chi }_a+4(-1+\eta ){\chi }_s\right)+\mathcal{O}(v^6)$$

$${\hat{A}}_{32}^{(\text{PN, S})}$$

$${\hat{A}}_{43}^{(\text{PN, S})}$$

$${\hat{A}}_{41}^{(\text{PN, S})}$$

Full text

I Introduction

Recent advances in gravitational-wave (GW) detections [1] have dramatically improved our ability to observe binary black hole (BBH) mergers, opening opportunities to study their dynamics, test general relativity, and describe the population of black hole binaries. The essential tool for these analyses is the “waveform,” which describes the time evolution of the strain that arrives at gravitational-wave detectors, such as The Laser Interferometer Gravitational-Wave Observatory [2], Virgo [3], and The Kamioka Gravitational Wave Detector [4].

The waveform of a BBH system in a quasicircular orbit is characterized by eight intrinsic parameters: the two black hole masses and the six spin components. Efficiently extracting these parameters from an observed gravitational waveform and, conversely, constructing accurate waveforms for different intrinsic parameters helps maximize the scientific return of GW observations.

The most precise method for constructing gravitational waveforms is by numerically solving Einstein’s field equations [5]. However, numerical relativity (NR) is computationally expensive and can only provide solutions for discrete points in the intrinsic parameter space. To overcome this limitation, approximation methods are necessary.

During the early inspiral phase, when the two compact objects are widely separated and moving slowly, the post-Newtonian (PN) approximation provides an effective description of the waveform [6]. However, as the binary approaches the merger, the PN expansion breaks down due to the onset of strong-field effects and nonlinearities. In this regime, NR is required to accurately model the dynamics and GW emission.

Spherical harmonics with spin weight $-2$ provide a basis for describing GW radiation of inspiraling BBHs at large distances from the source. Consequently, both PN and NR waveforms for binary signals are conventionally decomposed into modes using this basis [6]. Searches and Bayesian inference of signal parameters can then selectively retain modes, balancing desired precision with computational cost.

In the ringdown phase, the GW signal can be approximated by a perturbed Kerr solution as a superposition of quasinormal modes [7]. The symmetries of the system during ringdown are best described by a decomposition into spheroidal harmonics, which correspond to the well-known quasinormal mode frequencies of a perturbed Kerr black hole. The angular indices specifying spheroidal harmonics are different from those specifying the spherical harmonics used for NR waveforms, a phenomenon known as mode mixing [8].

In principle, the amplitudes and phases of spherical or spheroidal harmonic modes can be mapped back to the intrinsic parameters of the source. While the PN expansion provides analytical expressions for this mapping in the inspiral regime, extracting similar relationships near the merger is more challenging. Nevertheless, BBH detections dominated by these late-time regimes, like GW190521 [9] and GW231123 [10], are already part of the observational record. Understanding and modeling the structure of mode amplitudes and phases in the strong-field regime across the vast BBH parameter space will become increasingly important as more of these signals are detected.

Historically, the phase has received more attention in waveform modeling, primarily because matched filtering [11], the main technique used for GW detection and parameter estimation, is more sensitive to phase variations than amplitude. Consequently, phase evolution has been extensively studied in the literature. At the same time, several works have highlighted the importance of amplitude modeling, especially for tests of general relativity, and the inclusion of higher-order modes [12, 13, 14]. The goal of this study is therefore to contribute to this line of work by focusing exclusively on mode amplitudes.

Previous studies [15, 16, 17, 18, 19, 20, 21] have explored how mode amplitudes near merger and during the ringdown depend on the intrinsic parameters of the binary, both for spheroidal quasinormal modes in the ringdown and for spherical harmonic modes near merger. For instance, [16] and [17] investigated the dependence of the quasinormal mode spectrum on the masses and spins of the progenitor system. Using NR simulations from the bam NR code [22], they obtained fitting functions that can describe the amplitudes of the dominant and loudest subdominant quasinormal modes in the ringdown phase. With a similar aim, the authors in [19, 20] used a Gaussian process regression to obtain mappings between intrinsic parameters and the ringdown quasinormal mode amplitudes for noneccentric precessing simulations from the Simulating Extreme Spacetimes (SXS) catalog [23, 24]. Eccentric simulations were considered in [25], where the author used simulations from the Rochester Institute of Technology (RIT) catalog [26, 27] to obtain expressions for the ringdown mode amplitudes of nonspinning unequal-mass BBHs.

On the other hand, Ref. [15] used simulations from the SXS catalog to study the mode structure of spin-weighted spherical harmonics near the merger, using fits inspired by the leading-order dependence of the PN mode amplitudes on binary parameters. Surprisingly, they found that PN-inspired fits can describe mode amplitudes at fixed times near the merger, despite the breakdown of the PN approximation in this region. The relationship between PN and NR waveforms was also studied in Ref. [18] for eccentric BBH mergers using SXS, RIT, and MAYA [28, 29] simulations. The author similarly found that PN approximations remain valuable for constructing waveform models even close to the merger. However, discrepancies between the catalogs emerged as the mass ratio or eccentricity increased.

Inspired by previous studies suggesting that specific modes retain the leading-order PN dependence on intrinsic parameters close to the merger, we first investigate the stability of leading-order PN fits for spherical harmonic mode amplitudes from the late inspiral through the postmerger regime. We connect the early-time behavior of the recovered fit coefficients to the PN velocity to aid in the interpretation of the recovered values.

Second, we construct polynomial fits in mass ratio and spins to assess the role of beyond-leading-order dependencies in improving amplitude fits, particularly for modes where leading-order fits break down near merger. We aim to identify the simplest model required to capture the dominant physical dependencies for each mode.

Finally, for the first time, we use these methods to carry out a detailed comparison of NR simulations from the SXS, RIT, and MAYA NR catalogs. Cross-catalog comparisons provide a critical consistency check and help flag simulations where systematic differences may emerge.

The rest of this paper is organized as follows. Section II.1 provides an overview of the spherical harmonic decomposition and introduces relevant quantities and notation. Section II.2 briefly describes the different NR catalogs used in this study and outlines the simulation selection criteria. Sections II.3 and II.4 present the validity of the PN expansion and the construction of PN-inspired fits, along with the methods used to assess their performance and stability. Section II.5 introduces fits with beyond-leading-order parameter dependencies. Section III.2 discusses the stability and evolution of the fit coefficients of the PN-inspired fits, evaluates their performance, interprets their connection to the true PN expansion at earlier times, and compares results across NR catalogs. Section III.3 presents the results of fits with higher-order dependencies and evaluates their improvements over the leading-order models. Finally, Sec. IV summarizes our findings and outlines directions for future work.

We use a geometrized system of units in which Newton’s gravitational constant $G$ and the speed of light $c$ are $c=G=1$.

II Methods

II.1 Decomposition of gravitational-wave strain into spherical harmonic modes

The multipolar expansion of GW radiation from BBH mergers can be expressed in a basis of $-2$ spin-weighted spherical harmonics $Y_{\ell m}^{-2}(\theta,\phi)$ as follows:

[TABLE]

where $h$ is the complex gravitational-wave strain, $h_{\ell m}$ are the wave modes, $(\theta,\phi)$ are the angular coordinates specifying the direction of wave propagation, and $t$ is the retarded time [30].

Each complex $(\ell,m)$ mode can be written in polar form as

[TABLE]

where $A_{\ell m}(t)$ and $\phi_{\ell m}(t)$ are its amplitude and phase, respectively. The $(2,2)$ mode dominates the gravitational radiation emitted by a binary system [31, 32]. In this study, we analyze the amplitudes of different modes relative to the $(2,2)$ mode,

[TABLE]

The $\hat{A}_{\ell m}(t)$ are extracted over the interval $t\in(-500M,40M)$, where $t=0M$ corresponds to the peak amplitude of the $(2,2)$ mode.111The total mass $M$ of the binary sets the characteristic timescale of the system. In natural units, where $c=1$ and $G=1$, time can be expressed in terms of $M$, with $1M_{\odot}\approx 5\times 10^{-6}$ s. Figure 1 illustrates this interval on the real and imaginary parts of the full strain $h$ [Eq. (1)] for the simulation SXS:BBH:0305, which appears in Fig. 1 of the GW150914 discovery paper [33]. We highlight the situation of this interval near and around the merger.

The primary dynamical parameter for a nonspinning system is its mass ratio $q=\frac{m_{1}}{m_{2}}$. Using the convention $m_{1}\geq m_{2}$, such that $q\geq 1$, we can also define the symmetric mass ratio $\eta$ and the mass asymmetry parameter $\delta$,

[TABLE]

In aligned-spin binaries, where the spin vectors are parallel to the orbital angular momentum, the waveform is determined by $q$ and the spin magnitudes of the two bodies, given in dimensionless form as $\chi_{i}=\frac{S_{i}}{m_{i}^{2}}$, with $i=1,2$ labeling the object in the binary. It is also useful to define the symmetric and antisymmetric spin combinations,

[TABLE]

II.2 Numerical relativity simulations

Several NR groups have performed simulations of BBH mergers using different formalisms to set initial conditions and evolve the spacetime metric [23, 24, 26, 27, 28, 29, 34, 22]. In principle, simulations with identical initial conditions should produce the same gravitational waveform. However, differences in numerical precision, gauge choices, and waveform extrapolation techniques, among other factors, can introduce discrepancies between catalogs [35].

To compare results across catalogs, we use simulations from the SXS [23, 24], RIT [26, 27], and MAYA [28, 29] catalogs. In this section, we briefly review the methods employed by each catalog and outline our simulation selection criteria.

II.2.1 SXS catalog

The SXS catalog utilizes the spectral Einstein code [36] to perform simulations. The initial data are generated by solving the extended conformal thin-sandwich equations, which are discretized on a grid and solved using a spectral elliptic solver. Iterative tuning is applied to achieve the desired BBH properties, and the initial data are evolved for several orbits to reduce eccentricity and produce quasicircular orbits [24].

The evolution of the initial data is carried out using a first-order version of the generalized harmonic formulation [37, 38, 39] of Einstein’s equations, with constraint damping. A multidomain spectral method is applied, where the computational domain is partitioned into subdomains at various stages. A spectral adaptive mesh refinement algorithm determines when grid points should be added or subtracted in a given subdomain, while a control system and the associated maps shift the boundaries of the various subdomains with respect to each other in order to accommodate the requirement that the excision boundary must stay inside the apparent horizons at all times [40, 41]. These techniques rely on a set of specified error tolerances that are determined by the resolution of a simulation [23].

In this work, we use the highest available resolution level for each SXS simulation. The waveforms were read using the sxs python package [42, 43].

II.2.2 RIT catalog

The simulations in the RIT catalog are generated using the LazEv code [44]. The initial data for RIT simulations are derived from the Bowen-York solution and computed using a generalized version of the TwoPunctures code [45], which solves a coupled system of the Hamiltonian and momentum constraints. RIT simulations use PN approximations [46] to produce quasicircular orbits.

The LazEv code implements the Baumgarte-Shapiro-Shibata-Nakamura-Oohara-Kojima [47, 48, 49] formulation of Einstein’s equations. The code utilizes the cactus [50] infrastructure from the EinsteinToolkit [51] and uses the carpet [52] module for mesh refinement. Each simulation is labeled nXYY, where the grid spacing in the wave zone is specified as $M/X.YY$. For example, n120 corresponds to a grid spacing of $M/1.2$ in the wave zone.

Only simulations with a resolution label of n120 or higher are included from the RIT catalog. The simulations were read directly from publicly available strain data files [53].

II.2.3 MAYA catalog

The simulations in the MAYA catalog use the maya code [29]. Similar to LazEv, the maya code uses the cactus/carpet/EinsteinToolkit infrastructure for its simulations. It constructs initial data based on the Bowen-York extrinsic curvature, conformally flat spatial metric, and the TwoPunctures solver. The PN equations of motion are also used to define the initial spins and orbital momentum. The maya code uses the Baumgarte-Shapiro-Shibata-Nakamura [47, 48] formulation to evolve the initial data, employing mesh refinement for higher-resolution simulations.

The publicly available MAYA data files [54] were read using the mayawaves [55] python package.

II.2.4 Selection criteria

The linear shift in the time origin of the coordinates, which arises from the residual motion of the center of mass of the binary system [56], was corrected for all simulations used. Additionally, all simulations are selected to have initial orbital eccentricities of less than 0.002, as reported by their respective catalogs.222Eccentricity is not consistently defined across catalogs. While there are methods [57] for consistently calculating eccentricity using waveform data, reevaluating the eccentricity for all simulations considered here is beyond the scope of this work.

For nonspinning binary systems, simulations with any mass ratio are included, provided that the initial magnitude of each of the six dimensionless spin components is less than $10^{-4}$. Table 2 in Appendix B lists the nonspinning simulations used in this work.

For aligned-spin binary systems, simulations with fixed mass ratios $q\in\{1.0,1.5,2.0\}$ are used. The initial perpendicular components of the dimensionless spins are required to be less than $10^{-4}$. Additionally, $\chi_{i}\geq 0.1$ for at least one of the black holes to enhance the impact of spins on the waveforms. Table 3 lists the aligned-spin simulations used from the SXS, RIT, and MAYA catalogs.

We refer to the selected simulations as “datasets” for each NR catalog. After extracting the mode amplitudes $\hat{A}_{\ell m}$ from each dataset, their dependence on the binary parameters was modeled using fits motivated by the leading-order PN dependence.

II.3 Post-Newtonian expansion

The PN approximation describes the motion of inspiraling compact binaries and their emission of GWs [6]. It is an expansion in the characteristic internal velocity $v$ of the system, making the approximation valid for sources that are slowly moving and weakly stressed. Several works [6, 31, 58, 59] have computed the complex $h_{\ell m}$ modes at various PN expansion orders. The expressions for the relative mode amplitudes $\hat{A}_{\ell m}$ referenced throughout this work are shown in Appendix C and can be obtained from Refs. [6, 31, 59, 60, 58].

For nonspinning binaries, we focus on modes with $\ell\leq 4$ and $0<m\leq\ell$.333Nonprecessing systems exhibit mirror symmetry, meaning that the $-m$ modes are fully determined by the positive-$m$ modes through the relation $h_{\ell,-m}=(-1)^{\ell}h_{\ell m}^{*}$. For this reason, we focus only on modes with $m>0$. Among these, the $(2,2)$, $(2,1)$, $(3,3)$, and $(4,4)$ modes are generally the strongest [61].

For aligned-spin systems, we consider the modes $\{(\ell,m)\}=\{(2,2),(2,1),(3,2),(4,3),(4,1)\}$. We include the dominant $(2,2)$ mode due to its central role in waveform modeling, while the remaining modes are analyzed because they correspond to current multipoles. Modes with $\ell+m=\mathrm{odd}$ are entirely given by the current multipole moment, which is directly related to the system’s total angular momentum [6]. Hence, while spin effects are not present at leading order for any mode, these current multipole modes contain spin-dependent terms at next-to-leading order, making them more suitable for probing spin effects.

II.4 Investigating the leading-order post-Newtonian dependence

We examine how the GW mode amplitudes vary with the mass ratio and spins of the system over a range of times using fit Ansätze that are motivated by PN theory. Near the merger, the PN approximation breaks down. However, Ref. [15] suggested that the leading-order dependence of the mode amplitudes on the intrinsic parameters persists through the merger. The first objective of this study is to derive and reinterpret these results by directly fitting for $v$ in the PN expressions at each time step.

Following Ref. [15], a fit is constructed for each mode based on the leading-order functional form of the PN expressions (see Appendix C). Since this previous work found good agreement between their fits and the NR data at $-100M$, we expect $-500M$ to be a conservative choice for start time. Unlike Ref. [15], we look for correlations across modes by directly replacing powers of $v$ with independent fit coefficients instead of fitting for an overall factor. The fit coefficient that replaces the $v^{n}$ term in $\hat{A}_{\ell m}$ is denoted $a_{\ell m}^{(n)}$. We also define $a_{\ell m}\equiv\sqrt[n]{a_{\ell m}^{(n)}}$.

II.4.1 Leading-order fits for nonspinning simulations

For nonspinning simulations, the resulting leading-order fitting functions are

[TABLE]

To perform the fits, we extract $\hat{A}_{\ell m}$ from the NR simulation at each time step and apply the fitting functions presented above. This procedure yields a set of time-dependent fit coefficients $a_{\ell m}^{(n)}(t)$, which enable us to track their evolution and physical significance.

For nonspinning simulations, we apply the fits to the individual catalogs and to the joint dataset to assess their consistency.

One might expect the $a_{\ell m}(t)$ to agree across modes at early times, when the PN expansion is more valid, because they replace $v$ in the PN expressions, which is common to all modes. However, unless the amplitudes are taken early enough, such that the PN expansion can be truncated at leading order for all modes to good approximation, the fit coefficients need not agree between modes. This is because, on one hand, the velocity itself depends on the mass ratio and spins, and on the other hand, PN corrections beyond leading order can still be relevant at early times. As an illustration, if we write out dependencies on $t$ and $\eta$ explicitly, the $(2,2)$ mode is well approximated in the PN regime by Eq. (7a) if

[TABLE]

Therefore, $a_{22}^{(2)}(t)$ absorbs components of $\left[v(t;\eta)\right]^{2}$ that do not depend on $\eta$ as well as components of higher-order corrections to the amplitude that carry a linear $\eta$ dependence, like the leading-order term. We explore the relationship between $v$ and the $a_{\ell m}$ further by applying the fits at $t=-8000M$ for the subset of nonspinning SXS simulations that begin earlier than this time before merger, well into the early inspiral regime.

II.4.2 Leading-order fits for aligned-spin simulations

For aligned-spin systems, the resulting fits are

[TABLE]

where $\hat{A}_{\ell m}^{\mathrm{ns}}$ is the $(\ell,m)$ mode amplitude of a nonspinning system with the corresponding mass ratio. The amplitudes are extracted at each time step from aligned-spin simulations with fixed mass ratios $q=\{1.0,1.5,2.0\}$. With the time and mass ratio fixed, the amplitudes are treated as functions of $(\chi_{1z},\chi_{2z})$, and the fits above are applied.

In Ref. [15], the authors obtained results for aligned-spin systems separate from those for nonspinning systems. However, suppose the leading-order spin combination [e.g., $\chi_{a}+\delta\chi_{s}$ for the $(2,1)$ mode] truly governs the mode’s spin dependence through the merger. Whenever this effective spin combination vanishes, the amplitude should then closely match that of a nonspinning system with the same mass ratio. Therefore, to fully isolate the spin dependence, we set $\hat{A}_{\ell m}^{\mathrm{ns}}$ in the above fits to the NR amplitudes from SXS:BBH:4434 when $q=1.0$, SXS:BBH:3984 when $q=1.5$, and SXS:BBH:2497 when $q=2.0$. Note that odd-$m$ modes are not excited for equal-mass nonspinning systems, and there is liberty to pick the overall sign of the fit coefficients for systems with $q=1.0$.

II.5 Introducing higher-order dependencies

In previous work, the amplitudes of several modes, most notably the $(3,2)$ mode, but also the $(3,1)$, $(4,3)$, $(4,2)$, and $(4,1)$ modes, exhibited deviations from the leading-order PN fits at late times [15]. The second objective of this study is to develop improved models for the relative mode amplitudes by incorporating higher-order dependencies on the intrinsic parameters.

II.5.1 Higher-order fits for nonspinning simulations

We construct fitting functions for the relative mode amplitudes of nonspinning systems that incorporate higher-order dependencies on the symmetric mass ratio $\eta$. Our primary focus in this section is on the quality of the fit, rather than the physical interpretation of the coefficients. Therefore, we do not explicitly replace powers of $v$ with a fit coefficient. Instead, we allow the coefficients to absorb all overall factors and velocity dependencies, and we determine which modes and times are better described with these higher-order fits.

By inspection of the PN expressions in Appendix C, we observe that for modes with even $m$, higher-order PN terms introduce higher powers of $\eta$. In modes with odd $m$, the amplitudes additionally acquire an overall factor of $\delta$. Thus, we define a family of higher-order fitting functions $\hat{A}_{\ell m}^{(N)}$ for nonspinning systems as follows:

[TABLE]

These fits allow us to determine the minimum degree of higher-order dependencies needed to model deviations from leading-order fits.

II.5.2 Higher-order fits for aligned-spin simulations

Similarly, for aligned-spin systems, we define higher-order fitting functions that incorporate linear and quadratic dependencies on $\chi_{s}$ and $\chi_{a}$,

[TABLE]

As done with the leading-order fits, we fix the mass ratio and use the nonspinning NR reference amplitudes $\hat{A}_{\ell m}^{\mathrm{ns}}$.

II.6 Fitting procedure

The optimal fit coefficients are obtained using SciPy’s differential_evolution algorithm. This global optimization method evolves a population of candidate solutions to minimize a given cost function [62]. We set the cost function to the sum of squared residuals between the fit and the NR data.

Although one might expect the fitted leading-order coefficients, which stand in for $v$, to lie between [math] and $1$, these coefficients do not recover $v$ itself, as previously discussed. Since the leading-order coefficient absorbs higher-order corrections beyond the leading-order behavior, we allow it to vary over the broader range $[-5,5]$. In practice, the coefficients typically take on values much smaller in magnitude than the imposed bounds, indicating that the results are not driven by this prior range.

Because we apply the fit at each time step independently, we take the smooth variation of fit coefficients, as well as consistency across catalogs, as evidence against overfitting, since these indicate that the coefficients capture systematic trends rather than noise.

For higher-order fits, we find that the coefficients may saturate the bounds $[-5,5]$, and we therefore extend the allowed range to $[-10,10]$. As before, the smooth variation of the recovered coefficients suggests that the fits are well behaved. However, allowing such a wide parameter range may increase the risk of overfitting or introduce multiple local minima in the optimization, particularly for higher-order models with a larger number of free parameters. To address these concerns, in Appendix A we outline and perform a full Bayesian regression analysis for the higher-order models near merger. We find that the results generally align with those of the differential evolution algorithm.

The quality of each fit is evaluated using the Pearson correlation coefficient,

[TABLE]

where $\hat{A}_{\ell m}^{(\mathrm{fit})}$ represents the predicted amplitude values from the fitting functions, and $\hat{A}_{\ell m}^{(\mathrm{NR})}$ is the reference amplitude data obtained from NR simulations. Overbars denote the mean values of the datasets.

III Results and Discussion

III.1 Relationship between fit coefficients and PN velocity

To illustrate the relationship between $v$ and the fit coefficients $a_{\ell m}$, Fig. 2 shows the mode amplitudes at $t=-8000M$ for the strongest modes considered. We compare the best-fit curve obtained using Eqs. (7) to the amplitudes obtained by truncating the PN expressions at leading order. As can be seen in the figure, higher-order corrections are small but relevant at this early time, and the leading-order PN term does not fully capture the behavior of the mode amplitudes. Furthermore, the best-fit curves are not equivalent to the leading-order PN amplitudes. As previously mentioned, the fits absorb higher-order corrections, which may differ across modes. Therefore, the recovered fit coefficients, which are shown in Table 1, are close but not equivalent. Note also that fitting for the $a_{\ell m}$ at a given time produces a single number independent of $\eta$, while $v$, in general, depends on $\eta$. Consequently, although the $a_{\ell m}$ stand in for $v$ in the PN expressions and recover similar values at early times, they are not immediately comparable to a single reference value of $v$.

We note that the extrapolated NR waveforms used in this work are not in the same Bondi-Metzner-Sachs (BMS) frame as PN waveforms, and PN parameters need not coincide with quasilocal mass and spin measurements at the apparent horizons that are quoted in NR simulations. Recent work shows one can align the BMS frames when hybridizing NR and PN waveforms by optimizing over the PN parameters to find the best fit to the NR waveform [63, 64]. While this optimization over PN parameters is beyond the scope of this study, we expect that systematic differences between PN parameters and quasilocal NR measurements are implicitly absorbed into the fit coefficients. Discrepancies in the fitted coefficients at early times may partially reflect this parameter mismatch.

When applying the fits near merger, higher-order corrections may cause the optimal $a_{\ell m}^{(n)}$ to be negative. In the following sections, we quote $a_{\ell m}^{(n)}$ rather than $a_{\ell m}$.

III.2 Stability and time evolution of leading-order PN-inspired fits

III.2.1 Nonspinning simulations

Having established the relationship between $v$ and the optimal fit coefficients at early times, Fig. 3 shows the fit coefficients $a_{\ell m}^{(n)}$ obtained near the merger using individual catalogs and the joint dataset. The Pearson correlation coefficient is also shown as a function of time for each mode to quantify the agreement between the data and the fits.

For the all modes, the recovered fit coefficients are approximately consistent across catalogs and vary smoothly in time, suggesting the fits are stable. The Pearson correlation remains above $0.99$ at all times for the $(2,2)$, $(2,1)$, and $(3,3)$ modes across catalogs, except briefly for the MAYA dataset, which contributed the fewest simulations to our dataset.

For the weaker modes, there is significantly more variability across catalogs. The SXS catalog, which contributed the most simulations to our full dataset, shows the most stable behavior. It is followed by the RIT dataset, which shows erratic behavior only in the Pearson correlation for the $(4,2)$ and $(4,1)$ modes, but otherwise performs similarly to the SXS dataset. The MAYA dataset is stable for the $(3,2)$ and $(4,3)$ modes, but unstable for the $(3,1)$, $(4,2)$, and $(4,1)$ modes. These subdominant modes have lower amplitudes and are more susceptible to numerical errors. Additionally, higher-$\ell$ modes have higher frequencies and more complex sourcing, making them more susceptible to resolution-dependent differences [23, 35].

For all datasets, as the time approaches $t=0M$, the coefficients tend to grow as each mode absorbs increasingly significant higher-order corrections. Focusing on the results for the SXS data, we find that the fits maintain high correlations and vary smoothly over time for all modes throughout most of their evolution, up to $t\approx-10M$. However, the fit quality drops sharply around $t=0M$ for all weaker modes. This decline indicates that higher-order corrections introduce markedly different dependencies on $\eta$ that cannot be captured by the leading-order basis alone, and additional $\eta$ dependencies must be introduced.

III.2.2 Aligned-spin simulations

To reduce the impact of cross-catalog differences, we perform fits on the aligned-spin simulations using only the SXS dataset. Nonetheless, Appendix D presents snapshots of the SXS fits with RIT and MAYA data points included for visual comparison across catalogs.

Figure 4 shows the fit parameters over time. Early on, recovered parameters vary little across mass ratios. This is because higher-order corrections depending solely on mass ratio are accounted for by the nonspinning NR amplitudes $\hat{A}_{\ell m}^{(\mathrm{ns})}$ used in the fits. The consistency across mass ratios up to $t\approx-20M$ indicates that higher-order spin corrections introduce little to no additional mass-ratio dependence, beyond that captured at leading order. However, differences become more pronounced after the merger, where higher-order spin effects with distinct mass-ratio dependencies play a more significant role.

Regarding fit quality, the $(2,1)$ mode appears to be well described by its leading-order spin dependence through merger, with $C_{21}>0.99$ for $q=1.0$ and $C_{21}>0.97$ for $q=1.5$ and $q=2.0$ at all times. This is consistent with the results of Ref. [15]. Although the $(2,2)$ and $(3,2)$ modes maintain a high correlation with the fits until very close to merger, they deviate significantly from the PN dependence near merger for all mass ratios. For the $(4,3)$ mode, $C_{43}>0.97$ for $q=1.0$ at all times, but this correlation decreases with increasing mass ratio, especially near merger.

Finally, the $(4,1)$ fit also performs well for $q=1.0$ up to $t\approx 0M$, but drops sharply near merger, where the coefficient suddenly increases in magnitude. Furthermore, for $q=1.5$ and $q=2.0$, the $(4,1)$ fit is unstable throughout the evolution, showing small jumps in $a_{41}^{(4)}$ and fluctuations in $C_{41}$. This suggests the model is not well constrained by the data and is likely prone to overfitting the amplitudes for this mode. This uncertainty in the model is evident in the fit snapshots shown in Figs. 10–12 in Appendix D.

The result in [15] found better agreement for higher mass ratios in the $(4,1)$ mode. However, their analysis allowed a constant offset to vary freely, whereas we fixed this offset by enforcing the nonspinning limit through our $\hat{A}_{\ell m}^{(\mathrm{ns})}$. This suggests that some higher-order spin effects were possibly absorbed into their constant offset and may not be captured in our approach.

III.3 Higher-order fits

We now investigate the extent to which higher-order dependencies are needed to significantly improve fit quality for nonspinning and aligned-spin simulations.

III.3.1 Nonspinning simulations

The impact of higher-order corrections on the quality of fits for nonspinning systems is illustrated in two complementary figures. Figure 5 shows the recovered fit coefficients and Pearson correlations when applying fits of degrees $N=1$ and $N=2$ to each mode, while Fig. 6 provides snapshots of the amplitude data and these fits at selected times. Together, these figures allow for a visual assessment of how the inclusion of higher-order terms improves modeling accuracy across modes and times.

Figure 5 focuses on times $t\in[-100M,40M]$, where higher-order $\eta$ terms are expected to have larger impact. At earlier times ($t\ll 0M$), leading-order fits were already found to perform well across all modes in the results presented in the previous section, as well as previous work in Ref. [15]. Furthermore, up to $t\sim-100M$, the higher-order fits recover coefficients consistent with leading-order fits, as can be seen in the figure.

As the system approaches merger ($t\to 0$), and especially in the postmerger regime ($t>0$), higher-order $\eta$ dependencies become increasingly relevant, leading to noticeable improvements in fit quality for several modes. We find that for all modes except the $(4,4)$ mode, the $N=3$ fit offers minimal improvement over the $N=2$ fit. Furthermore, the recovered parameters for the $N=3$ model vary nonsmoothly in time and increase to large values, indicating possible overfitting. We therefore only show fits for $N=1$ and $N=2$ in Fig. 5, except for the $(4,4)$ mode, which we discuss below. We observe the following trends:

$(2,2)$ mode: As previously found, the leading-order PN-inspired fit remains well behaved through the merger. Adding a constant offset ($N=1$) does not change the fit appreciably, while adding an $\eta^{2}$ dependence slightly improves the correlation for $t\gtrsim-25M$. 2. 2.

$(2,1)$ and $(3,3)$ modes: These modes slightly benefit from the inclusion of a $\delta\eta$ term ($N=1$). Improvements are noticeable after the merger, as seen at $t=15M$ and $t=25M$ in Fig. 6. However, improvements beyond $N=1$ are minimal. 3. 3.

$(4,4)$ mode: Higher-order terms become significantly more important after the merger. In this regime, the $N=2$ and $N=3$ fits yield larger $C_{41}$ values than the lower-order models, and this improvement is clearly visible in Fig. 6. For the $N=3$ model, we find that the fit coefficients $c_{2}$ and $c_{3}$ saturate the bounds $[-10,10]$ near merger, although they vary smoothly and remain well within these bounds at earlier times. To ensure that the fit is not constrained by the imposed bounds, we extend the allowed range to $[-100,100]$ for $t>0M$. While this extension could raise concerns about overfitting, we explicitly test this possibility in Appendix A.1 by performing a full Bayesian linear regression using the $N=3$ model. The Bayesian results are consistent with those obtained using differential evolution, and the preference for an $\eta^{3}$ dependence is supported by both methods. 4. 4.

Subdominant modes with $\ell=3,4$ and $m<\ell$: These modes, previously identified as poorly described by leading-order PN fits near merger, show a strong preference for higher-order fits near and beyond the merger. Nevertheless, the $N=1$ and $N=2$ fits can capture their amplitude evolution with correlations comparable to the fits for the stronger modes. These improvements are apparent at late times in the snapshots in Fig. 6.

Figure 6 also plots data points from the RIT and MAYA catalogs. Although these were not used to generate the fits, they generally align well with the SXS-based fits. However, cross-catalog dispersion increases at later times for subdominant modes, particularly for the $(4,2)$ and $(4,1)$ modes.

III.3.2 Aligned-spin simulations

Figure 7 plots the recovered fit coefficients and Pearson correlations of the higher-order aligned-spin fits $\hat{A}_{lm}^{(i,\text{S})}$ defined in Eqs. (11) for the $(2,2)$, $(2,1)$, and $(3,2)$ modes, with the corresponding values from the leading-order fits included for reference. To further illustrate the behavior of the fits, Figs. 10–12 present snapshots of the amplitude data and their associated fits at selected times for each mass ratio considered.

The $c_{s}$ and $c_{a}$ for the linear and quadratic models closely follow their corresponding values from the leading-order fits up to near merger, as expected given the good performance of the leading-order fits at earlier times. For the $(2,2)$ mode, the $\chi_{s}^{2}$ contributions increase significantly after merger, and the quadratic fit yields a much higher Pearson correlation than both the linear and leading-order fits. For $q=1.0$, the coefficients $c_{a}$, $c_{as}$ and $c_{aa}$ remain small for the $(2,2)$ mode, but increase with increasing mass ratios. Similarly, the spin dependence of the $(3,2)$ mode is better captured by the quadratic fit than the lower-order models. In this case, contributions from $\chi_{s}^{2}$ and $\chi_{aa}$ are significant around merger for equal-mass systems, while contributions from $\chi_{a}$ and $\chi_{a}\chi_{s}$ increase with increasing mass ratio. The improvement of the quadratic fit over the leading-order fit for the $(2,2)$ and $(3,2)$ modes is visually evident in Figs. 10–12. For the $(2,1)$ mode, the leading-order fits already yield high correlations. Nevertheless, the added flexibility of the linear fits increases the correlation for unequal-mass systems near and after merger. Quadratic fits provide further improvement, although, as seen in Figs. 10–12, this improvement is modest.

For the $(4,3)$ and $(4,1)$ modes, we find that the fits are sensitive to the stochastic nature of differential evolution, particularly at higher mass ratios. This sensitivity likely arises from degeneracies among the five fit coefficients combined with relatively sparse coverage of simulations. Rather than using global optimization, we therefore perform Bayesian linear regression for these modes, as outlined in Appendix A. Figure 8 shows the posterior distributions of each parameter as a function of time obtained using the quadratic fit. For equal-mass systems, the parameters of the $(4,3)$ mode are well constrained and consistent with those of the leading-order fit. The inclusion of the $c_{as}$ term slightly improves the Pearson correlation relative to the leading-order fit after the merger. In contrast, for the $(4,1)$ mode in equal-mass systems, in equal-mass systems, the fit exhibits occasional spikes, and improvement over the leading-order fit is not clearly established.

For $q=1.5$ systems, the simulations used do not sufficiently constrain the parameters of either the $(4,3)$ or $(4,1)$ modes. This is likely due to the smaller number of data points at this mass ratio, combined with the increasing importance of higher-order spin dependencies. As illustrated by the $q=2.0$ case, the parameters are better constrained when more simulations are available. However, higher-order spin effects become significant near merger, and the model fails to capture a clean evolution of the coefficients in this regime. These results suggest that higher-order spin dependencies are particularly important around merger for these modes, and that additional data would be necessary to robustly constrain the evolution of their coefficients.

Overall, the inclusion of quadratic spin terms improves the modeling of the $(2,2)$ and $(3,2)$ modes for all mass ratios considered and of the $(4,3)$ mode for mass ratios $q>1$ after the merger. Visual comparison of the RIT and MAYA simulations in Figs. 10–12 confirms cross-catalog consistency with these improved higher-order fits for the $(2,2)$, $(2,1)$, $(3,2)$, and $(4,3)$ modes. For the $(4,1)$ mode, however, there is greater dispersion across catalogs, again likely due to differences in numerical techniques across catalogs.

IV Conclusions

In this study, we analyzed the mode amplitudes of quasicircular, nonprecessing BBH mergers using NR data. Our primary goal was to understand the extent to which the PN dependence on mass ratios and spins is carried over from the inspiral, particularly near and after the merger, where modeling remains challenging.

We obtained fit coefficients for the mode amplitudes relative to the $(2,2)$ mode as a function of time and verified their stability across the waveform, from late inspiral ($t=-500M$) to postmerger ($t=+40M$). The stability and quality of these fits demonstrate that for several modes, the time dependence can be cleanly separated from the dependence on the mass ratio and spin. This separation provides a framework for constructing simplified waveform models in the near- and postmerger regime by directly modeling the time evolution of the fit coefficients.

In the inspiral regime, we connected the fit coefficients to the PN velocity. As the system evolves toward the merger, the coefficients deviate from the PN predictions, but we interpret these deviations as absorptions of higher-order effects that preserve the same overall dependence on mass ratio and spin. This reinterpretation provides intuition behind why these PN-inspired fits can, for some modes, be applied in the strong-field regime.

To explore higher-order parameter dependencies, we tested fits of increasing complexity in both mass ratio and spin. For nonspinning systems, we found that PN-inspired fits with polynomial order $N\leq 3$ can model mode amplitudes accurately even close to the merger. While leading-order fits remain effective through the late inspiral, selectively incorporating higher-order terms substantially improves the modeling of subdominant modes near and after the merger. For aligned-spin simulations, we found that quadratic fits in $\chi_{s}$ and $\chi_{a}$ offer notable improvements for modes that were not well modeled by the leading-order PN fits, such as the $(2,2)$ and $(3,2)$ modes. However, the large number of free parameters in quadratic fits can lead to instabilities. For higher-order modes of unequal-mass systems, for which there are fewer available simulations and higher-order spin dependencies become more important, these models are not well constrained.

We validated these results across the SXS, RIT, and MAYA catalogs. While most fits were performed using only SXS data, we generally found good visual agreement with the other catalogs. Discrepancies in a small number of cases may point to numerical artifacts and can serve as indicators of outlying simulations.

The PN expansion does not capture physical effects that are unique to the strong-field regime. For example, near merger, the deformation of the black holes’ dynamical horizons leads to tidal heating through the flow of energy and angular momentum through the horizons. These processes increase the horizon area and change the masses and spins of the black holes, effects that are not typically incorporated in PN waveforms [65, 66, 67, 68]. We therefore emphasize that the phenomenological fits presented here are merely inspired by the functional forms of leading-order PN mode amplitudes. Their physical interpretation cannot be fully described within a PN framework, as they effectively absorb contributions from strong-field, near-merger physics that lie beyond the formal domain of validity of the PN approximation.

Our conclusions do not address the accuracy of PN-inspired fits for modeling the phase evolution. A complementary analysis of phase, which would require additional methodological considerations, is left for future work. We also aim to extend this analysis to simulations of precessing systems, in which the spin components are time dependent.

Overall, this study advances the understanding of mode amplitude structure in BBH mergers and lays a framework for efficient amplitude-consistent waveform modeling, particularly near the merger.

Acknowledgements.

We wish to acknowledge B. S. Sathyaprakash for his invaluable guidance and insightful comments throughout this work. We are grateful to Vaishak Prasad for his assistance with the SXS catalog and for providing useful comments on this manuscript. We further acknowledge the SXS, RIT, and MAYA Collaborations for their publicly available catalogs of NR simulations, which were used in this study. We also thank Carlos Lousto for guidance in using the RIT catalog, and Deborah Ferguson and Aasim Jan for their help with mayawaves and the MAYA catalog. Finally, we are grateful to Ssohrab Borhanian for answering questions related to his original study and to Koustav Chandra, Rossella Gamba, and Dongze Sun for helpful discussions.

Appendix A Bayesian linear regression

When performing Bayesian regression, we assume the NR amplitudes $\hat{A}_{\ell m}^{(\mathrm{NR})}$ are described by one of the models in Eqs. (9)-(11), which are characterized by a set of fit coefficients $\{c_{i}\}\equiv\vec{c}$. We use a Gaussian likelihood with unknown variance $\mathcal{L}(\hat{A}_{\ell m}^{(\mathrm{NR})}|\,\vec{c}\,)$ and uniform priors $\pi(c_{i})=\mathcal{U}(c_{min},c_{max})$ to compute the posterior distribution for each $c_{i}$ from Bayes’s theorem,

[TABLE]

where $\mathcal{Z}(\hat{A}_{\ell m}^{(\mathrm{NR})})$ is the normalizing evidence,

[TABLE]

which is especially useful for model comparison. In particular, the log Bayes factor between models $1$ and $2$,

[TABLE]

is positive if model $1$ is preferred and negative when model $2$ is preferred.

We compute posteriors using the Dynesty nested-sampling algorithm [69] accessed through the bilby package [70, 71].

A.1 $N=3$ fit for the $(4,4)$ mode

As discussed in Sec. III.3, the coefficients $c_{2}$ and $c_{3}$ of the $N=3$ model rail against the previously imposed bounds on the parameters $[-10,10]$ for the $(4,4)$ mode. To examine whether extending the prior bounds on $c_{2}$ and $c_{3}$ leads to overfitting the data, we perform Bayesian regression using uniform priors $\pi(c_{i})=\mathcal{U}(-100,100)$ for all $\{c_{i}\}$.

Figure 9 shows the posterior distributions obtained at each time step, together with snapshots of the fitted models constructed from posterior samples. All coefficients remain well within the prior range and are consistent with the results presented in Sec. III.3 obtained using the differential evolution algorithm. The coefficients $c_{2}$ and $c_{3}$ exhibit relatively large uncertainties. As shown in the snapshots, this uncertainty mostly affects predictions at small $\eta$ and is likely driven by the limited number of data points in this regime, i.e., the limited number of simulations for systems with high-mass ratio. While the data do not tightly constrain the exact value of $c_{3}$, they favor the presence of the $\eta^{3}$ term over its absence. This is evidenced by the positive value of the log Bayes factor at late times and visually corroborated in the snapshots.

Appendix B Selected simulations

Tables 2 and 3 list the simulations from the SXS, RIT, and MAYA catalogs that were used in this study, following the selection criteria outlined in Sec. II.2. Table 2 lists the nonspinning simulations with their mass ratio, while Table 3 lists the aligned-spin simulations with their mass ratio and dimensionless spin components.

Appendix C Post-Newtonian expressions for relative mode amplitudes

The PN waveform mode amplitudes relative to the $(2,2)$ mode up to $\mathcal{O}(v^{5})$ are given by [6], and we rewrite them as $\hat{A}_{\ell m}$ in Eq. (16),

[TABLE]

where the total mass and luminosity distance have been set to unity to match convention in NR simulations.

The spin contributions to the mode amplitudes can be found in Refs. [59, 60, 58]. We rewrite these relative to the $(2,2)$ mode up to $\mathcal{O}(v^{6})$ in Eq. (17).

[TABLE]

Note that these expressions are given relative to the $(2,2)$ mode, and thus differ from the expressions for individual $(\ell,m)$ mode amplitudes typically reported in the literature.

Appendix D Aligned-Spin Fit Snapshots

Figures 10–12 compare the leading-order PN-inspired fits and the quadratic fits, omitting the linear fits for simplicity. For these figures, the coefficients are inferred using the Bayesian regression procedure outlined in Appendix A, and the fit surfaces are drawn from the resulting posteriors.