Numerical Relativity Papers

Understanding when inflationary perturbations become genuinely nonlinear near the horizon crossing requires methods that go beyond both linear perturbation theory and the gradient expansion. In this work, we introduce a nonlinear lattice framework for single-field inflation based on a shear-free, locally Friedmann-Lemaître-Robertson-Walker geometry. This approach captures inhomogeneous local expansion rates, curvature contributions to the local Friedmann equation, and proper-volume weighting at a fraction of the computational cost of full numerical relativity. We construct fully nonlinear $δN$ observables on uniform-density slices, together with other practical time-dependent estimators for the curvature perturbations. After validating the framework in a standard slow-roll regime, we apply it to Starobinsky's linear-potential model featuring an intermittent ultra-slow-roll (USR) phase and a sharp potential transition. During this non-attractor USR regime, the lattice captures the separation of curvature perturbation estimators, the growth and subsequent stabilisation of non-Gaussianity, and a transient weakening of the shear-free approximation when the inflaton velocity becomes very small. Our framework provides a practical intermediate approach between rigid background lattice simulations and full numerical relativity, offering a nonlinear bridge between lattice methods, the $δN$ formalism, and the stochastic inflation formalism.

We present 3D numerical-relativity evolutions of the unstable Ellis-Bronnikov wormhole using GRTeclyn, starting from exact isotropic initial data for the coupled Einstein-phantom-scalar system. With a flat initial lapse (alpha=1) and full phantom support, truncation-level noise eventually drives the rarefactive instability and rapid throat expansion. To force a clean collapse while breaking spherical symmetry, we reduce the phantom stress-energy support to S_support=0.5 and add a quadrupolar scalar-field perturbation (A_phi=+0.02, sigma_phi=0.5). The resulting compressive evolution forms a trapped surface and emits a gravitational-wave signal whose peak propagates between extraction radii at v approx c, distinct from superluminal CCZ4 constraint modes. After horizon formation the swallowed phantom matter triggers a violent rebound ("phantom bounce") that launches an outward curvature shock. For the moderate perturbation amplitude simulated here, an intermediate-mass (10^3 solar mass) wormhole at D=1 Mpc falls slightly below the Advanced LIGO design sensitivity; detection requires either closer sources, larger initial asymmetries, or next-generation detectors.

Accurate inclusion of both spin precession and orbital eccentricity effects in gravitational waveform models represents a key hurdle in our ability to fully characterize the properties of compact binaries. Virtually all efforts to model precession rely on a coprecessing frame transformation, a time-dependent spatial rotation that tracks the dominant emission direction and simplifies the waveform morphology. We assess the utility of the coprecessing frame transformation to separate out the effect of the precession of the orbital plane from the waveform in the presence of non-negligible orbital eccentricity. We rely on 20 numerical relativity simulations, which include the complete physical effects of spin precession and eccentricity in the strong-field, and compare waveforms in both the inertial and coprecessing frames. Comparing against the eccentric, spin-aligned model SEOBNRv5EHM, we find that while the waveform mismatches decrease in the coprecessing frame, they remain above the level required for accurate waveform modeling, $\sim$ 0.01 or higher for large inclinations. Further improvements, e.g., modeling mode asymmetries as already pursued for quasicircular binaries, will likely prove essential. We also find that by removing the dominant amplitude and phase modulations from the waveform, the coprecessing frame facilitates surrogate modeling, achieving lower errors at a fixed number of basis elements compared to the inertial frame. Our results demonstrate both the utility and the limitations of the coprecessing frame as a cornerstone in waveform modeling for eccentric and precessing binaries.

Binary black hole (BBH) mergers detected via gravitational waves are addressing key open questions in astrophysics, cosmology, and fundamental physics. Our scientific conclusions rely on extracting accurate source parameters, for which we require accurate signal modelling. It is well known that current BBH waveform models need to be improved for high-mass-ratio, strongly precessing systems, and in this paper we provide a concrete illustration of this issue, showing that the degradation in model performance is substantially more severe than might have been anticipated. We present numerical relativity (NR) simulations of precessing BBH systems with a mass ratio of 18 and a dimensionless spin of 0.8 on the larger black hole (with the smaller black hole non-spinning), covering five values of spin misalignment. We assess the accuracy of state-of-the-art waveform models in this region of parameter space by computing the standard mismatch between the models and the NR waveforms. We find that all current waveform models often exhibit significant mismatches ($\gtrsim$0.1), indicating poor performance in this regime. We also perform limited parameter estimation using a subset of state-of-the-art waveform models, injecting these NR simulations as signals into the three-detector LIGO-Virgo network. In some cases we find errors in mass measurements of over 100%, dramatically illustrating that substantial improvements are required in existing waveform models. The numerical simulations presented here will be valuable for calibrating future BBH waveform models in this region of parameter space.

The dynamics and observable signatures of neutron star mergers are governed by physics under the most extreme conditions. They are particularly impacted by the high-density equation of state, which for the most sophisticated models is usually available in the form of tables. Numerical relativity codes usually evolve particularly well-behaved numerical ("conservative") variables, but at the price that the physically interesting ("primitive") variables need to be found at every computational element and at every integration sub-step by means of expensive (and not always successful) root-finding algorithms. We have recently developed the Lagrangian numerical relativity code SPHINCS_BSSN which evolves the spacetime on an adaptive mesh with well tested methods, but the fluid is evolved by means of freely moving particles. Since our evolution equations differ from those of conventional numerical relativity, we need to develop new conservative-to-primitive algorithms if we want to use tabulated equations of state. We present here three such algorithms: a 3D and a 2D Newton-Raphson method and a 1D root-finding algorithm based on Ridders' method. We find the 3D method to be very fast and robust with an average failure fraction in a full-blown neutron star merger simulation (with the DD2 equation of state) well below 1%. While we do not find obvious advantages for the 2D method, the 1D Ridders' method is slow, but essentially fail-safe. Therefore, we choose the 3D Newton-Raphson as default and fall back to the 1D Ridders' method as a safe "parachute".

In this work we investigate the weak lensing convergence using an end-to-end nonlinear general relativistic framework. Combining numerical relativity simulations of large-scale structure formation with general relativistic ray-tracing, we compare our nonlinear calculation to the expectation based on perturbation theory for a set of 20 synthetic observers. We focus on large angular scales $\ell < 100$ across a broad range of redshifts with $0.05

The Bondi--van der Burg--Metzner--Sachs (BMS) frame of gravitational waves produced by numerical relativity (NR) simulations is crucial for building accurate waveform models. A proper comparison of NR waveforms with other models requires fixing the arbitrary BMS frame. In this work we improve the center-of-mass (CoM) frame fixing for quasicircular, nonprecessing binary systems. Past work approximated the CoM motion with just a linear fit. We compute a post-Newtonian result of the boosted CoM charge to also capture its physical out-spiraling oscillations. We show that using the analytical results improves the robustness of the fit parameters -- translation and boost vectors -- to the choice of duration and time of the fitting window. Our analysis demonstrates a maximum improvement in robustness when the window is placed at the center of the inspiral. We quantified this improvement by computing the ratio of variances of fit parameters when the fit window size is varied. The largest improvement in robustness of parameters is by a factor of $\sim 25$ for the boost vector and $\sim 20$ for the translation vector. Finally, we incorporate this method into the BMS frame-fixing routine of the python package $\text{scri}$ for waveforms produced with Cauchy-characteristic evolution.

We present BLAZ4R, the first living catalog of confirmed $z>4$ blazars, with a focus on the contribution of eROSITA, on board of the Spectrum Roentgen Gamma (SRG) spacecraft. Blazars at $z>4$ are rare but powerful probes of how active supermassive black holes evolve in connection to relativistic jets, in the first 2 billion years of cosmic history. At these redshifts, X-ray observations are essential for constraining blazars jet power and orientation, enabling effective trace of their parent population. The all-sky surveys provided by eROSITA ensure X-ray detection for BLAZ4R sources and, in some cases, allow spectral and temporal studies of their jetted emission. BLAZ4R includes 54 confirmed blazars, characterized through their X-ray properties, radio spectra and morphology, and multiwavelength profiles. We confirm that jetted sources are significantly more numerous relative to non-jetted counterparts at high-$z$, and that blazars (and therefore the overall jetted population) do not exhibit significantly different features compared to the total active galactic nuclei population in the early Universe. Fast accretion processes that involve relativistic jets are clearly required to justify the existence of fully formed jetted AGN at $z>4$. However, the diverse multiwavelength properties characterizing BLAZ4R do not yet allow us to identify the specific signatures of these processes. We will continue updating BLAZ4R to search for such signatures and ultimately understand the early formation of jetted AGN.

Using state-of-the-art scattering results in post-Minkowskian (PM) gravity, we improve the tidal sector of four different flavors of the effective-one-body (EOB) formalism. We notably explore both adiabatic and post-adiabatic gravitoelectric and gravitomagnetic quadrupolar tidal effects at the next-to-next-to-leading PM-order. When comparing the predictions of the so-constructed Lagrange-PM-tidal version of EOB to recent numerical-relativity data on the scattering of neutron stars, we find improved agreement with respect to existing EOB models and PM expansions. Our work lays the foundation for the development of an accurate tidal sector of the PM EOB models, and points out the need to explore improved resummation schemes in PN EOB for bound and circularized orbits.

We consider a broad class of static, spherically symmetric generalized Schwarzschild-like solutions with multiple non-interacting anisotropic fluid sources and derive the coordinate transformation from Schwarzschild-like (curvature) to isotropic coordinates with conformally flat spatial slices. The isotropic form removes spatial-sector coordinate pathologies at the horizon, clarifies geometric quantities (e.g., ADM mass and curvature invariants), and enables the construction of well-posed initial data on t=const hypersurfaces, suitable for the Hamiltonian and conformal formulations of numerical relativity and for perturbation theory. The backgrounds in isotropic coordinates we develop make it straightforward to separate environmental effects from intrinsic strong-gravity signals and meet the growing interest in non-vacuum black hole phenomenology across scattering, lensing, and waveform modeling.

Testing general relativity in the strong-field and highly dynamical regime is now possible through current gravitational-wave observations, where even a single high-quality detection can place competitive constraints on deviations from Einstein's theory. The parametrized post-Einsteinian framework provides a theory-agnostic approach to search for such deviations, but it typically assumes that systematic uncertainties in the base waveform model, particularly those arising from calibration to numerical relativity, are negligible. In this work, we investigate how calibration errors in the late-inspiral fitting coefficients of the IMRPhenomD waveform model can lead to spurious detections of departures from general relativity in parametrized tests. We use an uncertainty-aware version of IMRPhenomD, recalibrated to a set of numerical relativity surrogate waveforms and equipped with a probabilistic description of its fitting coefficients, to simulate general-relativity-consistent signals. We inject these signals into an O5 ground-based detector network and recover them with the original IMRPhenomD model augmented with a parametrized post-Einsteinian phase deformation. We find that false violations of general relativity using this model arise for network signal-to-noise ratios as low as 60. When the uncertainty-aware model is used instead, the inferred parametrized post-Einsteinian phase deformation remains consistent with zero even for signals with a signal-to-noise ratio up to 330. These results demonstrate the need to account for numerical relativity calibration uncertainty in order to perform reliable inspiral tests of general relativity. They also illustrate that explicitly incorporating numerical relativity calibration uncertainty into the waveform model preserves our ability to robustly test general relativity.

Next-generation gravitational-wave (GW) detectors, such as the Laser Interferometer Space Antenna (LISA), will observe vast numbers of overlapping signals. Disentangling these signals from instrumental noise and from one another constitutes a significant data analysis challenge. We explore a denoising technique based on embedding time series into Hankel matrices: a superposition of $n$ (damped) sinusoids corresponds to a matrix of rank $2n$. Thus, the problem of signal extraction is reduced to a structured low-rank approximation problem. Using synthetic data tailored to GW applications, we benchmark three Hankel-based algorithms: ESPRIT, Cadzow iterations, and iteratively reweighted least squares (IRLS). Our test scenarios include isolated and multi-component monochromatic signals, the resolution of sources with closely spaced frequencies, and the recovery of black hole quasinormal modes (QNM). All three algorithms achieve near-optimal performance consistent with Fisher matrix bounds, evidenced by an inverse-square scaling of the mismatch with the signal-to-noise ratio. Furthermore, a proof-of-concept application to numerical relativity waveforms validates the ability of these algorithms to extract QNM frequencies from ringdown signals. Hankel low-rank approximation therefore offers a transparent, computationally efficient avenue for preprocessing GW time series.

The merger-ringdown radiation from a black hole binary merger is accurately modeled by a sum of linear quasinormal modes (QNMs). Recently, a non-QNM ``direct wave" component of the radiation, associated with prompt emission from a plunging perturber, has been identified. Motivated by the behavior of null geodesics perturbed from the remnant light ring, the Backwards One Body (BOB) approach has been shown to model the full merger-ringdown radiation to high accuracy, while using only a minimal number of parameters. In this work, using the Pöschl--Teller potential, we first show how the BOB amplitude evolution can be recovered from the QNM pole contributions. We then apply rational filters to isolate the non-QNM content in BOB and numerical relativity waveforms. We show that BOB naturally captures the direct wave component of the merger radiation, explaining its accuracy near the waveform peak. Finally, we use BOB to show that the direct wave frequency is largely uncorrelated with the horizon frequency, even for high spin remnants, and instead tracks the News frequency at the time of the peak News amplitude.

Penrose's weak cosmic censorship conjecture asserts that spacetime singularities produced by gravitational collapse are generically hidden behind event horizons, thus preventing them from causally influencing distant observers and preserving the predictability of the exterior region. In this work, we probe this conjecture in a setup that deliberately violates one of its central assumptions - the dominant energy condition - by considering the spherical collapse of a phantom scalar field with negative energy density. In principle, such a field could produce a Schwarzschild geometry with negative mass and therefore no event horizon. Our aim is to assess whether, once the dominant energy condition is abandoned, the fully coupled evolution of matter and geometry can dynamically generate or expose naked singularities, thereby probing the robustness of cosmic censorship. To this end, we perform high-accuracy numerical relativity simulations based on fourth-order finite-difference schemes. Starting from smooth, asymptotically flat initial data representing regular phantom scalar wave packets, we follow their fully nonlinear evolution through collapse or dispersion. While an ordinary (positive-energy) scalar field exhibits the standard Choptuik critical behavior at the threshold of black-hole formation, the phantom field displays qualitatively different dynamics. For all amplitudes considered, we find no evidence for trapped surfaces, naked singularities, or alternative stationary end states. Instead, the phantom scalar field always disperses, suggesting that cosmic censorship remains dynamically preserved even in the presence of negative-energy matter.

In globular clusters, hierarchical mergers are among the most promising pathways to forming massive black holes such as GW231123. A key factor determining whether a merger-remnant black hole will be retained in these environments and thus participate in subsequent hierarchical mergers is the recoil kick velocity. Analytic models for the recoil velocity are currently employed in nearly all population-synthesis frameworks. We instead use a state-of-the-art recoil-kick model gwModel_flow_prec developed from a combination of numerical-relativity and black-hole perturbation-theory data, together with data-driven techniques such as normalizing flows and the post-Newtonian structure of the kick. Employing both back-of-the-envelope estimates and detailed N-body as well as semi-analytical cluster simulations, we show that gwModel_flow_prec leads to a noticeable increase in the retention probability of hierarchical-merger remnants compared to the previously used analytic model and changes the mass and spin distribution of the black holes formed through hierarchical mergers. Additionally, we discuss the implications of our results in the context of massive binaries such as GW231123.

Black hole-neutron star mergers, together with binary neutron star mergers, are key laboratories for neutron star physics. They enable us to probe merger dynamics imprinted in gravitational waves and potential electromagnetic counterparts. These systems link microphysics and macrophysics by placing constraints on the dense matter equations of state, potentially revealing the imprint of hadron-quark phase transitions, clarifying the role of neutrino irradiation in shaping the ejecta, its r-process nucleosynthesis, and kilonova emission, as well as assessing how magnetically driven instabilities affect mass ejection and possible electromagnetic signatures. Despite their importance, black hole-neutron star mergers remain relatively less studied and therefore not yet well understood, largely due to the lack of publicly available numerical relativity setups suitable for such investigations. In this work, we present a fully reproducible black hole-neutron star merger simulation performed exclusively using Einstein Toolkit thorns, targeting the detected event GW230529. The simulations are carried out at three resolutions with finest grid spacings of $162$, $222$ and $310$ meters to assess numerical robustness. The entire setup, from initial data to a parameter file with some of the analysis scripts, is publicly released as a new Einstein Toolkit gallery example and will be distributed as part of the Hypatia release, establishing a reference black hole-neutron star merger configuration within the Einstein Toolkit.

Many HPC applications that solve differential equations rely on the Runge-Kutta family of methods for time integration. Among these methods, the fourth-order accurate RK4 scheme is especially popular. This time integration scheme requires applications to evaluate four intermediate stages to take one time step. Depending on the complexity of the problem being solved, the evaluation of these intermediate stages can be computationally expensive. In this paper we develop explicit fourth-order accurate Multistep Runge-Kutta (MSRK) methods. The advantage of such methods is that they re-use data from previous time steps, thus requiring fewer intermediate stage evaluations and potentially speeding up applications. We outline a procedure to obtain and tune the method's coefficients by adjusting their stability regions in an attempt to maximize the size that a time step can take. We validate and evaluate our new methods in the context of Numerical Relativity applications using the EinsteinToolkit. We believe, however, that these methods and results should generalize to other applications using explicit Runge-Kutta methods.

The last evolutionary stages of massive black hole binaries prior to coalescence is dominated by the emission of gravitational waves, which will be probed by the future Laser Interferometer Space Antenna. If gas is present around the two black holes, however, the associated electromagnetic emission can provide additional information about the binary properties and location before the merger event. For this reason, a proper characterisation of the electromagnetic emission during these phases is of fundamental importance, and requires a detailed description of the gas dynamics close to the event horizon of the two black holes, only achievable via numerical simulations. Within this context, we present the implementation of the Superposed Kerr-Schild dynamic metric in the relativistic scheme in the meshless code GIZMO. Our code can now simulate black hole binaries approaching merger with high computational efficiency and accuracy, taking into account relativistic effects on the gas. To validate our implementation, we perform two tests. First, we explore the case of a relativistic Bondi flow around a binary, finding very good agreement with numerical relativity simulations. Then we explore the case of an inviscid relativistic circumbinary disc, comparing our results with a similar simulation run assuming Newtonian gravity. In this second case, we find moderate differences in the mass accretion rate and in the inflow dynamics, which suggest that the presence of a non-Keplerian potential and of apsidal precession in the orbiting gas trajectories may produce stronger shocks and boost angular momentum transport in the disc. Our work highlights the importance of accounting for relativistic corrections in accretion disc simulations around black hole binaries approaching merger, even at scales much larger than those currently probed by numerical relativity simulations.

Golden dark sirens - exceptionally well-localized gravitational-wave (GW) sources without electromagnetic counterparts - offer a powerful route to precision measurements of the Hubble constant, $H_0$, with next-generation (XG) detectors. The statistical promise of this method, however, places stringent demands on waveform accuracy and detector calibration, as even small systematic errors can dominate over statistical uncertainties at high signal-to-noise ratios. We investigate the impact of waveform-modeling systematics on golden dark siren cosmology using a synthetic population of binary black holes consistent with current GW observations and analyzed in the XG-detector era. By comparing state-of-the-art waveform models against numerical-relativity-based reference signals, we quantify modeling inaccuracies from both modeling and data-analysis perspectives and assess how they propagate into biases in luminosity distance, host-galaxy association, and single-event $H_0$ inference. We find that while current waveform models often allow recovery of statistically consistent $H_0$ posteriors, small waveform-induced biases can significantly affect three-dimensional localization and host galaxy ranking, occasionally leading to incorrect redshift assignments. We further derive order-of-magnitude requirements on detector calibration accuracy needed to ensure that calibration systematics remain subdominant for golden dark sirens observed with XG networks. To realize sub-percent $H_0$ measurements with golden dark sirens will require waveform and calibration accuracies that scale as $\mathcal{O}(ρ^{-2})$ with signal-to-noise ratio, motivating sustained advances in waveform modeling, numerical relativity, and detector calibration for the XG era.

We present the first fully non-linear evolutions of binary neutron star mergers in a moving-punctures approach in Einstein-scalar-Gauss-Bonnet gravity. We study both linear and quadratic-type couplings between the scalar and the Gauss-Bonnet invariant, and uncover new post-merger phenomena. These include an enhancement of the prompt collapse of a long-lived hyper-massive neutron star remnant and cases where the remnant develops a scalar configuration due to different scalarisation instabilities. This study initiates the exploration of beyond-General-Relativistic effects enhanced by the non-linear dynamics of the neutron star's fluid.

A Reissner-Nordström black hole (RNBH) enclosed in a cavity is known to be superradiantly unstable to charged scalar perturbations below a critical frequency. Inspired by the emergence of the QCD axion as a prominent dark matter candidate, we construct a model featuring an axion field coupled to an electromagnetic field that undergoes superradiant growth around an RNBH. Utilizing numerical relativity, we achieve stable, long-term evolution of this system and perform a comparative analysis across various parameter spaces. Our comprehensive investigation reveals the formation of a hairy black hole, whose final state is governed by a diverse set of physical parameters. Notably, the decay constant in the axion potential, representing nonlinear interactions, bifurcates the superradiant instability into two distinct behavioral regimes, leading to more significant dynamical shifts than previously reported. Furthermore, we examine the influence of the scalar field's charge and mass, as well as the mirror's position. We investigate the axionic bosenova process and observe a long-term beating pattern of the axion field induced by nonlinear interactions. By fine-tuning these parameter combinations, we demonstrate that the system can evolve toward a variety of distinct physical endpoints.

Within the framework of scalar-tensor theories of gravity, certain models can evade classical black hole no-hair theorems. A well-known example is Einstein-scalar-Gauss-Bonnet gravity, where black holes carrying a scalar charge can exist. We find that, within this theory, binary black holes initially described by General Relativity can acquire scalar charges once they reach a critical orbital separation ("dynamical scalarization"). We develop a simple semi-analytic model, based on the adiabatic conservation of the total Wald entropy, to estimate the scalar charge evolution during the binary inspiral. We also run fully nonlinear numerical-relativity simulations for different configurations, finding consistent results. The gravitational-wave phase difference between Einstein-scalar-Gauss-Bonnet and General Relativity waveforms, which we use to assess detectability, is also computed. We find that dynamical scalarization might be observable in nearly equal-mass binary black hole mergers with third-generation ground-based gravitational-wave detectors, in a narrow range of the dimensional coupling of the theory.

We study the properties of the principal symbol of the 3+1 equations of motion in Teleparallel Equivalent of General Relativity (TEGR) and assess the conditions for hyperbolicity. We use the Hamiltonian formulation based on the vectorial, antisymmetric, symmetric trace-free, and trace (VAST) decomposition of the canonical variables in the Hamiltonian formalism, and the Hamilton's equations previously presented in the literature. We study the system of differential equations at the linear level, and show that the principal symbol has a sector with imaginary eigenvalues, which renders the system not hyperbolic. This situation persists by taking spatial derivatives in either one or three coordinate directions, and it should be interpreted as a problem of the specific gauge choice instead of a general problem with TEGR. The first practical use of Hamilton's equations in this work can be extended for proving well-posedness in spherical symmetry, and establish numerical relativity setups in TEGR.