Numerical Relativity Research Articles

Performance-portable Binary Neutron Star Mergers with AthenaK

Abstract We introduce an extension to the AthenaK code for general-relativistic magnetohydrodynamics (GRMHD) in dynamical spacetimes using a 3+1 conservative Eulerian formulation. Like the fixed-spacetime GRMHD solver, we use standard finite-volume methods to evolve the fluid and a constrained-transport scheme to preserve the divergence-free constraint for the magnetic field. We also utilize a first-order flux correction (FOFC) scheme to reduce the need for an artificial atmosphere and optionally enforce a maximum principle to improve robustness. We demonstrate the accuracy of AthenaK using a set of standard tests in flat and curved spacetimes. Using a SANE accretion disk around a Kerr black hole, we compare the new solver to the existing solver for stationary spacetimes using the so-called “HARM-like” formulation. We find that both formulations converge to similar results. We also include the first published binary neutron star (BNS) mergers performed on graphical processing units (GPUs). Thanks to the FOFC scheme, our BNS mergers maintain a relative error of O ( 1 0 − 11 ) or better in baryon mass conservation up to collapse. Finally, we perform scaling tests of AthenaK on OLCF Frontier, where we show excellent weak scaling of ≥80% efficiency up to 32,768 GPUs and 74% up to 65,536 GPUs for a GRMHD problem in dynamical spacetimes with six levels of mesh refinement. AthenaK achieves an order-of-magnitude speedup using GPUs compared to CPUs, demonstrating that it is suitable for performing numerical relativity problems on modern exascale resources.

Formulating the complete initial boundary value problem in numerical relativity to model black hole echoes

Abstract In an attempt to simulate black hole echoes (generated by potential quantum-gravitational structure) in numerical relativity, we recently described how to implement a reflecting boundary outside of the horizon of a black hole in spherical symmetry. Here, we generalize this approach to spacetimes with no symmetries and implement it numerically using the generalized harmonic formulation. We cast the evolution equations and the numerical implementation into a Summation By Parts scheme, which seats our method closer to a class of provably numerically stable systems. We implement an embedded boundary numerical framework that allows for arbitrarily shaped domains on a rectangular grid and even boundaries that evolve and move across the grid. As a demonstration of this framework, we study the evolution of gravitational wave scattering off a boundary either inside, or just outside, the horizon of a black hole. This marks a big leap toward the goal of a generic framework to obtain gravitational waveforms for behaviors motivated by quantum gravity near the horizons of merging black holes.

First investigation of void statistics in numerical relativity simulations

First investigation of void statistics in numerical relativity simulations

On the Origins, Remnant, and Multimessenger Prospects of the Compact Binary Merger GW230529

This study investigates the origins of GW230529, delving into its formation from massive stars within isolated binary systems. Utilizing population-synthesis models, we present compelling evidence that the neutron star component forms second. However, the event’s low signal-to-noise ratio introduces complexities in identifying the underlying physical mechanisms driving its formation. Augmenting our analysis with insights from numerical relativity, we estimate the final black hole mass and spin to be approximately 5.3 M ⊙ and 0.53, respectively. Furthermore, we employ the obtained posterior samples to calculate the ejecta mass and kilonova light curves resulting from r-process nucleosynthesis. We find the ejecta mass to be within 0–0.06 M ⊙, contingent on the neutron star equation of state. The peak brightness of the kilonova light curves indicates that targeted follow-up observations with a Rubin-like observatory may have detected this emission.

Resolving elliptic boundary value problem and integral equation in relational partial metric space

In the present study, a novel version of the contraction theory on a relational partial metric space associated with a binary relation is exhibited. In this process, we observe that numerous relations can be used in many well-known fixed point theorems earlier introduced by multiple authors. We solve an elliptic boundary value problem and an integral equation to make the significance of the main conclusion evident. Exploring a nonlinear system, governed by nonlinear equations, underscores the importance of comprehending fixed points - pivotal states where the dynamics of a system remain unaltered.

Revisiting the apparent horizon finding problem with multigrid methods

Abstract Apparent horizon plays an important role in numerical relativity as it provides a tool to characterize the existence and properties of black holes on three-dimensional spatial slices in 3+1 numerical spacetimes. Apparent horizon finders based on different techniques have been developed. In this paper, we revisit the apparent horizon finding problem in numerical relativity using multigrid-based algorithms. We formulate the nonlinear elliptic apparent horizon equation as a linear Poisson-type equation with a nonlinear source, and solve it using a multigrid algorithm with Gauss-Seidel line relaxation. A fourth order compact finite difference scheme in spherical coordinates is derived and employed to reduce the complexity of the line relaxation operator to a tri-diagonal matrix inversion. The multigrid-based apparent horizon finder developed in this work is capable of locating apparent horizons in generic spatial hypersurfaces without any symmetries. The finder is tested with both analytic data, such as Brill-Lindquist multiple black hole data, and numerical data, including off-centered Kerr-Schild data and dynamical inspiraling binary black hole data. The obtained results are compared with those generated by the current fastest finder AHFinderDirect (Thornburg, Class. Quantum Grav. 21, 743, 2003), which is the default finder in the open source code Einstein Toolkit. Our finder performs comparatively in terms of accuracy, and starts to outperform AHFinderDirect at high angular resolutions (∼1°) in terms of speed. Our finder is also more flexible to initial guess, as opposed to the Newton’s method used in AHFinderDirect. This suggests that the multigrid algorithm provides an alternative option for studying apparent horizons, especially when high resolutions are needed.

Post-Minkowskian Theory Meets the Spinning Effective-One-Body Approach for Bound-Orbit Waveforms.

Driven by advances in scattering amplitudes and worldline-based methods, recent years have seen significant progress in our ability to calculate gravitational two-body scattering observables. These observables effectively encapsulate the gravitational two-body problem in the weak-field and high-velocity regime [post-Minkowskian (PM)], with applications to the bound two-body problem and gravitational-wave modeling. We leverage PM data to construct a complete inspiral-merger-ringdown waveform model for nonprecessing spinning black holes within the effective-one-body (EOB) formalism SEOBNR-PM. This model is closely based on the highly successful SEOBNRv5 model, used by the LIGO-Virgo-KAGRA Collaboration, with its key new feature being an EOB Hamiltonian derived by matching the two-body scattering angle in a perturbative PM expansion. The model performs remarkably well, showing a median mismatch against 441 numerical-relativity (NR) simulations that is somewhat lower than a similarly calibrated version of SEOBNRv5. Comparisons of the binding energy with NR also demonstrate better agreement than SEOBNRv5, despite the latter containing additional calibration to NR simulations.

ExaGRyPE: Numerical general relativity solvers based upon the hyperbolic PDEs solver engine ExaHyPE

ExaGRyPE describes a suite of solvers and solver ingredients for numerical relativity that are based upon ExaHyPE 2, the second generation of our Exascale Hyperbolic PDE Engine. Numerical relativity simulations are crucial in resolving astrophysical phenomena in strong gravitational fields and are fundamental in analyzing and understanding gravitational wave emissions. The presented generation of ExaGRyPE solves the Einstein field equations in the standard CCZ4 formulation under a 3+1 foliation and focuses on black hole space-times. It employs a block-structured Cartesian grid carrying a higher-order Finite Difference scheme with full support of adaptive mesh refinement (AMR), it facilitates massive parallelism combining message passing, domain decomposition and task parallelism, and it supports the injection of particles into the grid as static data probes or as moving tracers. We introduce the ExaGRyPE-specific building blocks within ExaHyPE 2, and discuss its software architecture and compute-n-feel.For this, we formalize the creation of any specific astrophysical simulation with ExaGRyPE as a sequence of lowering operations, where a few abstract logical tasks are successively broken down into finer and finer tasks until we obtain an abstraction level which can directly be mapped onto a C++ executable. The overall program logic is fully specified via a domain-specific Python interface, we automatically map this logic onto a more detailed set of numerical tasks, subsequently lower this representation onto technical tasks that the underlying ExaHyPE engine uses to parallelize the application, before eventually the technical tasks in turn are mapped onto small task graphs including the actual astrophysical PDE term evaluations, initial conditions, boundary conditions, and so forth. These can be injected manually by the user, or users might instruct the solver on the most abstract user interface level to use out-of-the-box ExaGRyPE implementations. We end up with a rigorous separation of concerns which shields ExaGRyPE users from technical details and hence simplifies the development of novel physical models. We present the simulations and data for the gauge wave, static single black holes and rotating binary black hole systems, demonstrating that the code base is mature and usable. However, we also uncover domain-specific numerical challenges that need further study by the community in future work.

Binary neutron star mergers in massive scalar-tensor theory: Properties of postmerger remnants

We investigate the properties of postmerger remnants of binary neutron star mergers in the framework of Damour–Esposito-Farese-type scalar-tensor theory of gravity with a massive scalar field by numerical relativity simulation. It is found that the threshold mass for prompt collapse is raised in the presence of the excited scalar field. Our simulation results also suggest the existence of a long-lived ϕ mode in hypermassive neutron stars due to the presence of the massive scalar field that enhances the quasiradial oscillation in the remnant. We investigate the descalarization condition in hypermassive neutron stars and discover a distinctive signature in postmerger gravitational waves. Published by the American Physical Society 2024

Phenomenological model of gravitational self-force enhanced tides in inspiraling binary neutron stars

Gravitational waves from inspiraling binary neutron stars provide unique access to ultradense nuclear matter and offer the ability to constrain the currently unknown neutron star equation-of-state through tidal measurements. This, however, requires the availability of accurate and efficient tidal waveform models. In this paper we present henom, a new phenomenological tidal phase model for the inspiral of neutron stars in the frequency domain, which captures the gravitational self-force informed tidal contributions of the time-domain effective-one-body model esum. henom is highly faithful and computationally efficient, and by choosing a modular approach, it can be used in conjunction with frequency-domain binary black hole waveform model to generate the complete phase for a binary neutron star inspiral. henom is valid for neutron star binaries with unequal masses and mass ratios between 1 and 3, and dimensionless tidal deformabilities up to 5000. Furthermore, henom does not assume universal relations or parametrized equations of state, hence allowing for exotic matter analyses and beyond standard model physics investigations. We demonstrate the efficacy and accuracy of our model through comparisons against esum, numerical relativity waveforms and full Bayesian inference, including a reanalysis of the binary neutron star observation GW170817. Published by the American Physical Society 2024

Mapping eccentricity evolutions between numerical relativity and effective-one-body gravitational waveforms

Orbital eccentricity in compact binaries is considered to be a key tracer of their astrophysical origin, and can be inferred from gravitational-wave observations due to its imprint on the emitted signal. For a robust measurement, accurate waveform models are needed. However, ambiguities in the definition of eccentricity can obfuscate the physical meaning and result in seemingly discrepant measurements. In this work we present a suite of 28 new numerical relativity simulations of eccentric, aligned-spin binary black holes with mass ratios between 1 and 6 and initial post-Newtonian eccentricities between 0.05 and 0.3. We then develop a robust pipeline for measuring the eccentricity evolution as a function of frequency from gravitational-wave observables that is applicable even to signals that span at least ≳7 orbits. We assess the reliability of our procedure and quantify its robustness under different assumptions on the data. Using the eccentricity measured at the first apastron, we initialize effective-one-body waveforms and quantify how the precision in the eccentricity measurement, and therefore the choice of the initial conditions, impacts the agreement with the numerical data. We find that even small deviations in the initial eccentricity can lead to non-negligible differences in the phase and amplitude of the waveforms. However, we demonstrate that we can reliably map the eccentricities between the simulation data and analytic models, which is crucial for robustly building eccentric hybrid waveforms, and to improve the accuracy of eccentric waveform models in the strong-field regime. Published by the American Physical Society 2024

Well-balanced High-order Finite Difference Weighted Essentially Nonoscillatory Schemes for a First-order Z4 Formulation of the Einstein Field Equations

We develop a new class of high-order accurate well-balanced finite difference (FD) weighted essentially nonoscillatory (WENO) methods for numerical general relativity (GR), which can be applied to any first-order reduction of the Einstein field equations, even if nonconservative terms are present. We choose the first-order nonconservative Z4 formulation of the Einstein equations, which has a built-in cleaning procedure that accounts for the Einstein constraints and that has already shown its ability in keeping stationary solutions stable over long timescales. By introducing auxiliary variables, the vacuum Einstein equations in first-order form constitute a partial differential equation system of 54 equations that is naturally nonconservative. We show how to design FD-WENO schemes that can handle nonconservative products. Different variants of FD WENO are discussed, with an eye to their suitability for higher-order accurate formulations for numerical GR. We successfully solve a set of fundamental tests of numerical GR with up to ninth-order spatial accuracy. Due to their intrinsic robustness, flexibility, and ease of implementation, FD-WENO schemes can effectively replace traditional central finite differencing in any first-order formulation of the Einstein field equations, without any artificial viscosity. When used in combination with well-balancing, the new numerical schemes preserve stationary equilibrium solutions of the Einstein equations exactly. This is particularly relevant in view of the numerical study of the quasi-normal modes of oscillations of relevant astrophysical sources. In conclusion, general relativistic high-energy astrophysics could benefit from this new class of numerical schemes and the ecosystem of desirable capabilities built around them.

Height-function-based 4D reference metrics for hyperboloidal evolution

Hyperboloidal slices are spacelike slices that reach future null infinity. Their asymptotic behaviour is different from Cauchy slices, which are traditionally used in numerical relativity simulations. This work uses free evolution of the formally-singular conformally compactified Einstein equations in spherical symmetry. One way to construct gauge conditions suitable for this approach relies on building the gauge source functions from a time-independent background spacetime metric. This background reference metric is set using the height function approach to provide the correct asymptotics of hyperboloidal slices of Minkowski spacetime. The present objective is to study the effect of different choices of height function on hyperboloidal evolutions via the reference metrics used in the gauge conditions. A total of 10 reference metrics for Minkowski are explored, identifying some of their desired features. They include 3 hyperboloidal layer constructions, evolved with the non-linear Einstein equations for the first time. Focus is put on long-term numerical stability of the evolutions, including small initial gauge perturbations. The results will be relevant for future (puncture-type) hyperboloidal evolutions, 3D simulations and the development of coinciding Cauchy and hyperboloidal data, among other applications.

Reconstructing the Genealogy of LIGO-Virgo Black Holes

We propose a Bayesian inference framework to predict the merger history of LIGO-Virgo binary black holes (BHs), whose binary components may have undergone hierarchical mergers in the past. The framework relies on numerical relativity predictions for the mass, spin, and kick velocity of the remnant BHs. This proposed framework computes the masses, spins, and kicks imparted to the remnant of the parent binaries, given the initial masses and spin magnitudes of the binary constituents. We validate our approach by performing an “injection study” based on a constructed sequence of hierarchically formed binaries. Noise is added to the final binary in the sequence, and the parameters of the “parent” and “grandparent” binaries in the merger chain are then reconstructed. This method is then applied to three GWTC-3 events: GW190521, GW200220_061928, and GW190426_190642. These events were selected because at least one of the binary companions lies in the putative pair-instability supernova mass gap, in which stellar processes alone cannot produce BHs. Hierarchical mergers offer a natural explanation for the formation of BHs in the pair-instability mass gap. We use the backward evolution framework to predict the parameters of the parents of the primary companion of these three binaries. For instance, the parent binary of GW190521 has masses 72−22+32M⊙ and 31−23+24M⊙ within the 90% credible interval. Astrophysical environments with escape speeds ≥100 km s−1 are preferred sites to host these events. Our approach can be readily applied to future high-mass gravitational wave events to predict their formation history under the hierarchical merger assumption.

Acoustic flow resonances in installed rectangular jets

Zonal Large Eddy Simulations (LES) of complex installed rectangular jet flows are performed for subsonic jet conditions of the Central Institute for Aviation Motors (CIAM) experiment, where a strong tone was reported. For far-field noise modelling, the zonal LES solutions are coupled with the permeable formulation of the Ffowcs Williams and Hawkings (FW-H) method. The obtained noise spectra predictions are compared with the acoustic measurements in the CIAM facility. To reduce statistical noise of the relatively short LES time series, contributions due to several dominant low frequency tones and the remaining broadband signal are computed separately, using a tailored Fourier transform technique for each signal type. For cross-validation, noise spectra measurements of an equivalent isolated round jet in the CIAM facility are compared with the empirical sJet model calibrated on the NASA Small Hot Jet Acoustic Rig (SHJAR) jet noise database. While the two datasets agree reasonably well for mid to high frequencies, larger discrepancies are observed for low frequencies, which are attributed to acoustic reflections in the non-anechoic CIAM facility. The differences in noise levels between the CIAM dataset and the sJet model representing the NASA jet noise data are used to determine the experimental uncertainty of the CIAM noise measurements for each frequency and observer angle. For the installed rectangular jet, it is shown that far-field noise spectra predictions of the zonal LES-FW-H method for both the fundamental tone, its sixth harmonic corresponding to the jet Strouhal number of around 1, and the broadband component are broadly within the experimental error bar for most observer angles. Some discrepancies between predictions of the zonal LES method and the CIAM measurements for the main low frequency tone for the 30° and 90°, where the experimental uncertainty at low frequencies is largest are also noted. After validating the acoustic solutions, spectral and conditional averaging analysis methods are applied to elucidate mechanisms of the acoustic-flow resonance in the installed rectangular jet flow. The numerical dispersion relation is obtained to analyse phase velocities of the upstream and downstream travelling pressure waves in the stream-wise direction. The conditional analysis of pressure fluctuations inside the jet reveals the emergence of internal reflection points corresponding to a rapid change of the acoustic wave phase as well as the saddle points interpreted as localised zones of vorticity shed from the side-wall edges. To independently investigate the effect of trailing edges and corner points of the jet embodiment geometry on closing the acoustic-flow resonance cycle, several modifications of the baseline rectangular nozzle are considered, such as introducing chevron-like lobes and cuts on the trailing edges of the side walls as well as smoothing the sharp corner points of the wall edges. Following the series of additional zonal LES runs with the modified installed nozzle geometries, an optimized tone-less modification of the original installed nozzle was obtained, in broad agreement with the conditional analysis results. The obtained tone-less modification of the installed rectangular jet is further analyzed to provide insights into its similarity and differences from the aeroacoustics of the reference isolated round jet in the same CIAM facility.

A review of gravitational memory and BMS frame fixing in numerical relativity

AbstractGravitational memory effects and the BMS freedoms exhibited at future null infinity have recently been resolved and utilized in numerical relativity simulations. With this, gravitational wave models and our understanding of the fundamental nature of general relativity have been vastly improved. In this paper, we review the history and intuition behind memory effects and BMS symmetries, how they manifest in gravitational waves, and how controlling the infinite number of BMS freedoms of numerical relativity simulations can crucially improve the waveform models that are used by gravitational wave detectors. We reiterate the fact that, with memory effects and BMS symmetries, not only can these next-generation numerical waveforms be used to observe never-before-seen physics, but they can also be used to test GR and learn new astrophysical information about our Universe.

Backreaction in Numerical Relativity: Averaging on Newtonian gauge-like hypersurfaces in Einstein Toolkit cosmological simulations

We introduce a spatial averaging scheme and use it to study the evolution of spatial averages in large-scale simulations of cosmological structure formation performed with the Einstein Toolkit. The averages are performed on the spatial hypersurfaces of the simulation setup which, at least initially, represent the hypersurfaces of statistical homogeneity and isotropy. We find only negligible cosmic backreaction on these hypersurfaces even on very small scales, but find significant curvature fluctuations of up to 10% in ΩR for sub-volumes with radius ∼200 Mpc and even larger fluctuations in smaller sub-volumes. In addition, we quantify fluid flow in and out of these sub-volumes. We find this to be significant, up to a 5% change in the density between redshift z=1 and z=0 of a single sphere of radius ∼200 Mpc (and larger for smaller spheres). We suggest this may be important for studies basing averages on volumes co-moving with the simulation hypersurfaces.

Worldtube excision method for intermediate-mass-ratio inspirals: Self-consistent evolution in a scalar-charge model

This is a third installment in a program to develop a method for alleviating the scale disparity in binary black hole simulations with mass ratios in the intermediate astrophysical range, where simulation cost is prohibitive while purely perturbative methods may not be adequate. The method is based on excising a “worldtube” around the smaller object, much larger than the object itself, replacing it with an analytical model that approximates a tidally deformed black hole. Previously [N. A. Wittek , Worldtube excision method for intermediate-mass-ratio inspirals: Scalar-field model in 3+1 dimensions, ] we have tested the idea in a toy model of a scalar charge in a fixed circular geodesic orbit around a Schwarzschild black hole, solving for the massless Klein-Gordon field in 3+1 dimensions on the p platform. Here we take the significant further step of allowing the orbit to evolve radiatively, in a self-consistent manner, under the effect of back-reaction from the scalar field. We compute the inspiral orbit and the emitted scalar-field waveform, showing a good agreement with perturbative calculations in the adiabatic approximation. We also demonstrate how our simulations accurately resolve postadiabatic effects (for which we do not have perturbative results). In this work we focus on quasicircular inspirals. Our implementation is publicly accessible in the p numerical relativity code. Published by the American Physical Society 2024

Self-interacting scalar dark matter around binary black holes

Gravitational waves can provide crucial insights about the environments in which black holes live. In this work, we use numerical relativity simulations to study the behavior of self-interacting scalar (wavelike) dark matter clouds accreting onto isolated and binary black holes. We find that repulsive self-interactions smoothen the “spike” of an isolated black hole and saturate the density. Attractive self-interactions enhance the growth and result in more cuspy profiles, but can become unstable and undergo explosions akin to the superradiant bosenova that reduce the local cloud density. We quantify the impact of self-interactions on an equal-mass black hole merger by computing the dephasing of the gravitational-wave signal for a range of couplings. We find that repulsive self-interactions saturate the density of the cloud, thereby reducing the dephasing. For attractive self-interactions, the dephasing may be larger, but if these interactions dominate prior to the merger, the dark matter can undergo bosenova during the inspiral phase, disrupting the cloud and subsequently reducing the dephasing. Published by the American Physical Society 2024

Element Formation in Radiation-hydrodynamics Simulations of Kilonovae

Understanding the details of r-process nucleosynthesis in binary neutron star merger (BNSM) ejecta is key to interpreting kilonova observations and identifying the role of BNSMs in the origin of heavy elements. We present a self-consistent, two-dimensional, ray-by-ray radiation-hydrodynamic evolution of BNSM ejecta with an online nuclear network (NN) up to a timescale of days. For the first time, an initial numerical relativity ejecta profile composed of the dynamical component and spiral-wave and disk winds is evolved including detailed r-process reactions and nuclear heating effects. A simple model for the jet energy deposition is also included. Our simulation highlights that the common approach of relating in postprocessing the final nucleosynthesis yields to the initial thermodynamic profile of the ejecta can lead to inaccurate predictions. Moreover, we find that neglecting the details of the radiation-hydrodynamic evolution of the ejecta in nuclear calculations can introduce deviations of up to 1 order of magnitude in the final abundances of several elements, including very light and second r-process peak elements. The presence of a jet affects element production only in the innermost part of the polar ejecta, and it does not alter the global nucleosynthesis results. Overall, our analysis shows that employing an online NN improves the reliability of nucleosynthesis and kilonova light-curve predictions.