The recently reported binary black hole merger, GW231123, has unusual properties that make it hard to explain astrophysically. Parameter estimation studies are consistent with maximally spinning black holes and the dimensionless spin of the more massive component is constrained to be $\chi_1\gtrsim 0.8$. Analysis of data also revealed potential systematics that could not be fully replicated with simulated studies. We explore the possibility that these measurements are biased due to unmodeled non-Gaussian noise in the detectors, and that the actual black hole spins are more modest. We present evidence for a population of \textit{microglitches} in LIGO gravitational-wave strain data that can lead to biases in the parameter estimation of short-duration signals such as GW231123. Using simulated data of a massive event like GW231123, we demonstrate how microglitches can bias our measurements of black hole spins toward $\chi\approx1$ with negligible posterior support for the true value of $\chi\approx0.7$. We develop a noise model to account for microglitches and show that this model successfully reduces biases in the recovery of signal parameters. We characterize the microglitch population in real interferometer data surrounding GW231123 and find a single detector glitch duty cycle of $0.57_{-0.19}^{+0.21}$, which implies nearly a $100\%$ probability that at least one event through the fourth gravitational wave transient catalog coincides with microglitches in two detectors. We argue that further investigations are required before we can have a confident picture of the astrophysical properties of GW231123.

This is an introduction to Bayesian inference with a focus on hierarchical models and hyper-parameters. We write primarily for an audience of Bayesian novices, but we hope to provide useful insights for seasoned veterans as well. Examples are drawn from gravitational-wave astronomy, though we endeavor for the presentation to be understandable to a broader audience. We begin with a review of the fundamentals: likelihoods, priors, and posteriors. Next, we discuss Bayesian evidence, Bayes factors, odds ratios, and model selection. From there, we describe how posteriors are estimated using samplers such as Markov Chain Monte Carlo algorithms and nested sampling. Finally, we generalize the formalism to discuss hyper-parameters and hierarchical models. We include extensive appendices discussing the creation of credible intervals, Gaussian noise, explicit marginalization, posterior predictive distributions, and selection effects.

Stellar theory predicts a forbidden range of black-hole masses between ${\sim}50$--$130\,M_\odot$ due to pair-instability supernovae, but evidence for such a gap in the mass distribution from gravitational-wave astronomy has proved elusive. Early hints of a cutoff in black-hole masses at ${\sim} 45\,M_\odot$ disappeared with the subsequent discovery of more massive binary black holes. Here, we report evidence of the pair-instability gap in LIGO--Virgo--KAGRA's fourth gravitational wave transient catalog (GWTC-4), with a lower boundary of $45_{-4}^{+5} M_\odot$ (90\% credibility). While the gap is not present in the distribution of \textit{primary} masses $m_1$ (the bigger of the two black holes in a binary system), it appears unambiguously in the distribution of \textit{secondary} masses $m_2$, where $m_2 \leq m_1$. The location of the gap lines up well with a previously identified transition in the binary black-hole spin distribution; binaries with primary components in the gap tend to spin more rapidly than those below the gap. We interpret these findings as evidence for a subpopulation of hierarchical mergers: binaries where the primary component is the product of a previous black-hole merger and thus populates the gap. Our measurement of the location of the pair-instability gap constrains the $S$-factor for $^{12}\rm{C}(\alpha,\gamma)^{16}\rm{O}$ at 300keV to $256_{-104}^{+197}$ keV barns.

With the release of the fourth LIGO--Virgo--KAGRA gravitational-wave catalog (GWTC-4), we are starting to gain a detailed view of the population of merging binary black holes. The formation channels of these black holes is not clearly understood, but different formation mechanisms may lead to subpopulations with different properties visible in gravitational-wave data. Adopting a phenomenological approach, we find GWTC-4 data supports the presence of at least three subpopulations, each associated with a different range of black hole mass and with sharp transition boundaries between them. Each subpopulation is characterized by different distributions for either the mass ratios, the black-hole spin magnitudes or both. Subpopulation A with primary mass $m_1 \leq 27.7^{+4.1}_{-3.4} M_{\odot}$ ($90 \%$ credibility), is characterized by a nearly flat mass ratio distribution $q=m_2/m_1$, and by small spin magnitudes ($\chi \leq 0.5^{+0.1}_{-0.1}$). Subpopulation B with $27.7^{+4.1}_{-3.4} M_{\odot} \leq m_1 \leq 40.2^{+4.7}_{-3.2} M_{\odot}$, has a much sharper preference for mass ratio $q \approx 1$. Subpopulation C, with $m_1 \geq 40.2^{+4.7}_{-3.2} M_{\odot}$, has support for large spin magnitudes, and tentative support for mass ratios $q\approx0.5$. We interpret these transitions as evidence for multiple subpopulations, each potentially associated with a different formation pathways. We suggest potential formation scenarios for each subpopulations, and suggest that Subpopulation B may be associated with chemically homogeneous evolution or population III stars. Our findings for Subpopulation C are largely consistent with recent claims of hierarchical mergers, but with some curious differences in properties.

Time series analysis is ubiquitous in many fields of science including gravitational-wave astronomy, where strain time series are analyzed to infer the nature of gravitational-wave sources, e.g., black holes and neutron stars. It is common in gravitational-wave transient studies to apply a tapered window function to reduce the effects of spectral artifacts from the sharp edges of data segments. We show that the conventional analysis of tapered data fails to take into account covariance between frequency bins, which arises for all finite time series -- no matter the choice of window function. We discuss the origin of this covariance and show that as the number of gravitational-wave detections grows, and as we gain access to more high signal-to-noise ratio events, this covariance will become a non-negligible source of systematic error. We derive a framework that models the correlation induced by the window function and demonstrate this solution using both data from the first LIGO--Virgo transient catalog and simulated Gaussian noise.

Gravitational-wave observations of massive, rapidly spinning binary black holes mergers provide increasing evidence for the dynamical origin of some mergers. Previous studies have interpreted the mergers with primary mass $\gtrsim45\,M_\odot$ as being dominated by hierarchical, second-generation mergers, with rapidly spinning primaries being the products of previous black hole mergers assembled in dense stellar clusters. In this work, we reveal confident evidence of another subpopulation with rapid and isotropic spins at low mass containing the two exceptional events GW241011 and GW241110, consistent with a hierarchical merger hypothesis. Our result suggests the mass distribution of the second-generation black holes is peaked at low primary masses of $\sim16\,M_\odot$ rather than $\gtrsim45\,M_\odot$ in the pair-instability gap. Such low-mass second-generation black holes must be formed from the merger of even lighter first-generation black holes, implying that dense, metal-rich stellar environments contribute to the binary black hole population. By separating the contamination of higher-generation black holes, our result reveals the primary mass distribution of first-generation black holes formed from stellar collapse, which shows a significant dip between $\sim12\,M_\odot$ to $\sim20\,M_\odot$. This may indicate a dearth of black holes due to variation in the core compactness of the progenitor.

The ringdown portion of a binary black hole merger consists of a sum of modes, each containing an infinite number of tones that are exponentially damped sinusoids. In principle, these can be measured as gravitational-waves with observatories like LIGO/Virgo/KAGRA, however in practice it is unclear how many tones can be meaningfully resolved. We investigate the consistency and resolvability of the overtones of the quadrupolar $\ell = m = 2$ mode by starting at late times when the gravitational waveform is expected to be well-approximated by the $\ell m n = 220$ tone alone. We present a Bayesian inference framework to measure the tones in numerical relativity data. We measure tones at different start times, checking for consistency: we classify a tone as stably recovered if and only if the 95\% credible intervals for amplitude and phase at time $t$ overlap with the credible intervals at all subsequent times. We test a set of tones including the first four overtones of the fundamental mode and the 320 tone and find that the 220 and 221 tones can be measured consistently with the inclusion of additional overtones. The 222 tone measurements can be stabilised when we include the 223 tone, but only in a narrow time window, after which it is too weak to measure. The 223 tone recovery appears to be unstable, and does not become stable with the introduction of the 224 tone. We find that $N=3$ tones can be stably recovered simultaneously. However, when analysing $N \geq 4$ tones, the amplitude of one tone is consistent with zero. Thus, within our framework, one can identify only $N=3$ tones with non-zero amplitude that are simultaneously stable.

Current population models of binary black hole distributions are difficult to interpret because standard population inferences hinge on modeling choices, which can mask or mimic real structure. The maximum population likelihood ``$\pistroke$ formalism'' provides a means to investigate and interpret features in the distribution of binary black holes using only data -- without specifying a population model. It tells us if features inferred from current population models are truly present in the data or if they arise from model misspecification. It also provides guidance for developing new models by highlighting previously unnoticed features. In this study, we utilize the $\pistroke$ formalism to examine the binary black hole population in the LIGO--Virgo--KAGRA (LVK) fourth Gravitational-Wave Transient Catalog (GWTC-4). Our analysis supports the existence of a gap around $45\,M_\odot$ in the secondary black hole mass distribution and identifies a widening in the distribution of the effective inspiral spin parameter $\chi_\text{eff}$ near this mass as recently reported by Tong et al. (2025). Similar to earlier studies, we find support for an anti-correlation between $\chi_\text{eff}$ and mass ratio. However, we argue that this may be a spurious correlation arising from misspecification of the joint distribution of black hole masses. Furthermore, we identify support for dimensionless black hole spin magnitudes at approximately $\chi \approx 0.2$ and $\chi\approx0.7$. The data support the existence of a correlation between the spin magnitudes $\chi_1$ and $\chi_2$, though subsequent study is required to determine if this feature is statistically significant. The accompanying data release includes $\pistroke$ samples, which can be used to compare theoretical predictions to LVK data and to assess assumptions in parameterised models.

The detection of gravitational waves from merging black holes with masses $\sim\,80-150\,\mathrm{M_\odot}$ suggests that some proportion of black hole binary systems form hierarchically in dense astrophysical environments, as most stellar evolution models cannot explain the origin of these massive black holes through isolated binary evolution. A significant fraction of such mergers could occur in Active Galactic Nuclei disks (AGN), however connecting individual black hole mergers to host galaxies is a challenging endeavor due to large localization uncertainties. We assess the feasibility of determining the fraction of hierarchically merging black hole binaries by computing the angular cross-correlation between gravitational wave localization posteriors and galaxy catalog skymaps. We forecast when the clustering of gravitational wave sky localizations can be measured accurately enough to distinguish the AGN origin scenario from hierarchical mergers in galaxies that do not host AGN. We find that if the observed merging population is dominated by binaries formed dynamically in AGN, then this could be determined with $\mathcal{O}(5000)$ mergers detected at the sensitivity that is projected for the upcoming A\# gravitational wave detectors.

SN\,2025kg, linked to EP250108a, is among the brightest broad-lined Type Ic supernova (SN Ic-BL) known, showing unique helium absorptions, a late-time broad H$\alpha$, and an early bump. In this {\em{Letter}}, we propose a jet-cocoon origin to explain EP250108a as off-axis cooling emission from a mildly relativistic inner cocoon viewed at $\sim45^\circ$ and the early bump of SN\,2025kg as the outer cocoon cooling emission, both constraining an energy of $\sim(1-2)\times10^{52}{\rm{erg}}$ and a progenitor radius of $\sim5\,R_\odot$. To explain SN\,2025kg's exceptionally luminous peak, potential energy injection into the $\sim2.5\,M_\odot$ ejecta from a magnetar with initial period $\sim1.7\,{\rm{ms}}$ and magnetic field $\sim2\times10^{15}{\rm{G}}$ may be required, implying a rapidly rotating $\sim4\,M_\odot$ progenitor. Thus, the progenitor may be a low-mass helium star with an extended helium envelope, supported by helium absorption lines and an inferred weak pre-SN wind. Hydrogen-rich material may reside in the inner ejecta layers, as suggested by the late-time broad H$\alpha$, possibly originating from main-sequence companion material evaporated by the magnetar wind. Since the observed near-solar metallicity challenges the popular quasi-chemically homogeneous evolution channel, the rapidly rotating helium-star progenitor of EP250108a/SN\,2025kg might attain angular momentum by being tidally spun up by a main-sequence companion in a close binary formed through isolated binary evolution.

It has become increasingly useful to answer questions in gravitational-wave astronomy using transdimensional models where the number of free parameters can be varied depending on the complexity required to fit the data. Given the growing interest in transdimensional inference, we introduce a new package for the Bayesian inference Library (Bilby) called tBilby. The tBilby package allows users to set up transdimensional inference calculations using the existing Bilby architecture with off-the-shelf nested samplers and/or Markov Chain Monte Carlo algorithms. Transdimensional models are particularly helpful when we seek to test theoretically uncertain predictions described by phenomenological models. For example, bursts of gravitational waves can be modelled using a superposition of N wavelets where N is itself a free parameter. Short pulses are modelled with small values of N whereas longer, more complicated signals are represented with a large number of wavelets stitched together. Other transdimensional models have found use describing instrumental noise and the population properties of gravitational-wave sources. We provide a few demonstrations of tBilby, including fitting the gravitational-wave signal GW150914 with a superposition of N sine-Gaussian wavelets. We outline our plans to further develop the tbilby code suite for a broader range of transdimensional problems.

Primordial gravitational waves are expected to create a stochastic background encoding information about the early Universe that may not be accessible by other means. However, the primordial background is obscured by an astrophysical foreground consisting of gravitational waves from compact binaries. We demonstrate a Bayesian method for estimating the primordial background in the presence of an astrophysical foreground. Since the background and foreground signal parameters are estimated simultaneously, there is no subtraction step, and therefore we avoid astrophysical contamination of the primordial measurement, sometimes referred to as "residuals". Additionally, since we include the non-Gaussianity of the astrophysical foreground in our model, this method represents the statistically optimal approach to the simultaneous detection of a multi-component stochastic background.

Bayesian parameter estimation is fast becoming the language of gravitational-wave astronomy. It is the method by which gravitational-wave data is used to infer the sources' astrophysical properties. We introduce a user-friendly Bayesian inference library for gravitational-wave astronomy, Bilby. This python code provides expert-level parameter estimation infrastructure with straightforward syntax and tools that facilitate use by beginners. It allows users to perform accurate and reliable gravitational-wave parameter estimation on both real, freely-available data from LIGO/Virgo, and simulated data. We provide a suite of examples for the analysis of compact binary mergers and other types of signal model including supernovae and the remnants of binary neutron star mergers. These examples illustrate how to change the signal model, how to implement new likelihood functions, and how to add new detectors. Bilby has additional functionality to do population studies using hierarchical Bayesian modelling. We provide an example in which we infer the shape of the black hole mass distribution from an ensemble of observations of binary black hole mergers.

We carry out astrophysical inference for compact binary merger events in LIGO-Virgo's first gravitational-wave transient catalog (GWTC-1) using a physically motivated calibration model. We demonstrate that importance sampling can be used to reduce the cost of what would otherwise be a computationally challenging analysis. We show that including the physical estimate for the calibration error distribution has negligible impact on the inference of parameters for the events in GWTC-1. Studying a simulated signal with matched filter signal-to-noise ratio $\text{SNR}=200$, we project that a calibration error estimate typical of GWTC-1 is likely to be negligible for the current generation of gravitational-wave detectors. We argue that other sources of systematic error---from waveforms, prior distributions, and noise modelling---are likely to be more important. Finally, using the events in GWTC-1 as standard sirens, we infer an astrophysically-informed improvement on the estimate of the calibration error in the LIGO interferometers.

The next generation of gravitational-wave observatories will achieve unprecedented strain sensitivities with an expanded observing band. They will detect ${\cal O}(10^5)$ binary neutron star (BNS) mergers every year, the loudest of which will be in the band for $\approx 90$ minutes with signal-to-noise ratios $\approx 1500$. Current techniques will not be able to determine the astrophysical parameters of the loudest of next-gen BNS signals. We show that subtleties arising from the rotation of the Earth and the free-spectral range of gravitational-wave interferometers dramatically increases the complexity of next-gen BNS signals compared to the one-minute signals seen by LIGO--Virgo. Various compression methods currently relied upon to speed up the most expensive BNS calculations -- reduced-order quadrature, multi-banding, and relative binning -- will no longer be effective. We carry out reduced-order inference on a simulated next-gen BNS signal taking into account the Earth's rotation and the observatories' free-spectral range. We show that standard data compression techniques become impractical, and the full problem becomes computationally infeasible, when we include data below $\approx 16$Hz -- a part of the observing band that is critical for precise sky localisation. We discuss potential paths towards solving this complex problem.

We present 3D radiation hydrodynamics simulations of common-envelope (CE) evolution involving a 12 solar mass red supergiant donor and a 3 solar mass companion. Existing 3D simulations are predominantly adiabatic, focusing strongly on low-mass donors on the red giant and asymptotic giant branches. However, the adiabatic assumption breaks down once the perturbed CE material becomes optically thin or when entering a longer-timescale evolutionary phase after the dynamical plunge-in. This is especially important for high-mass red supergiant donors, which have short thermal timescales, adding significant uncertainty to our understanding of how massive binary stars evolve into gravitational-wave sources, X-ray binaries, stripped-envelope supernovae, and more. We compare our radiation hydrodynamics simulations with an adiabatic simulation from Paper I that is otherwise identical, finding that radiative diffusion strongly inhibits CE ejection. The fraction of ejected mass is roughly half that of the adiabatic case without accounting for recombination energy release. Almost no material is ejected during the dynamical plunge-in, and longer-timescale ejection during the slow spiral-in is suppressed. However, the orbital separation reached at the end of the dynamical plunge-in does not differ significantly. The large amount of remaining bound mass tentatively supports the emerging view that the dynamical plunge-in is followed by a non-adiabatic phase, during which a substantial fraction of the envelope is ejected and the binary orbit may continue to evolve.

Gravitational-wave astronomy provides a promising avenue for the discovery of new physics beyond general relativity as it probes extreme curvature and ultra-relativistic dynamics. However, in the absence of a compelling alternative to general relativity, it is difficult to carry out an analysis that allows for a wide range of deviations. To that end, we introduce a Gaussian process framework to search for deviations from general relativity in gravitational-wave signals from binary black hole mergers with minimal assumptions. We employ a kernel that enforces our prior beliefs that - if gravitational waveforms deviate from the predictions of general relativity - the deviation is likely to be localised in time near the merger with some characteristic frequency. We demonstrate this formalism with simulated data and apply it to events from Gravitational-Wave Transient Catalog 3. We find no evidence for a deviation from general relativity. We limit the fractional deviation in gravitational-wave strain to as low as 7% (90% credibility) of the strain of GW190701_203306.

Intermediate mass ratio inspiral (IMRI) binaries -- containing stellar-mass black holes coalescing into intermediate-mass black holes ($M>100M_{\odot}$) -- are a highly anticipated source of gravitational waves (GWs) for Advanced LIGO/Virgo. Their detection and source characterization would provide a unique probe of strong-field gravity and stellar evolution. Due to the asymmetric component masses and the large primary, these systems generically excite subdominant modes while reducing the importance of the dominant quadrupole mode. Including higher order harmonics can also result in a $10\%-25\%$ increase in signal-to-noise ratio for IMRIs, which may help to detect these systems. We show that by including subdominant GW modes into the analysis we can achieve a precise characterization of IMRI source properties. For example, we find that the source properties for IMRIs can be measured to within $2\%-15\%$ accuracy at a fiducial signal-to-noise ratio of 25 if subdominant modes are included. When subdominant modes are neglected, the accuracy degrades to $9\%-44\%$ and significant biases are seen in chirp mass, mass ratio, primary spin and luminosity distances. We further demonstrate that including subdominant modes in the waveform model can enable an informative measurement of both individual spin components and improve the source localization by a factor of $\sim$10. We discuss some important astrophysical implications of high-precision source characterization enabled by subdominant modes such as constraining the mass gap and probing formation channels.