Abstract Two-dimensional axisymmetric simulations of binary neutron star (BNS) merger remnant are a cheap alternative to 3D simulations. To maintain realism for secular timescales, simulations must avoid accumulated errors from drifts in conserved quantities and artificial heating, and they must model turbulent transport in a way that remains plausible throughout the evolution. It is also crucial to avoid numerical artifacts due to the polar coordinate axis singularity. Methods that behave well near the axis often break flux-conservative form of the hydrodynamic equations, resulting in significant drifts in conserved quantities. We present a flux-conservative scheme that maintains smoothness near the axis without sacrificing conservative formulation of the equations or incurring drifts in conserved global quantities. We compare the numerical performance of different treatments of the hydrodynamic equations when evolving a hypermassive neutron star resembling the remnant of a BNS merger. These simulations demonstrate that the new scheme combines the axis smoothness of non-conservative methods with the mass and angular momentum conservation of other conservative methods on ~10 2 ms timescales of viscous and neutrino-driven evolution. Because fluid profiles remain smooth in the remnant interior, it is possible to remove artificial heating by evolving the entropy density. We show how physical heating and cooling terms can be easily calculated from source terms of the conservative evolution variables and demonstrate our implementation. Finally, we discuss and implement improvements to the effective viscosity scheme to better model the effect of magnetohydrodynamic instabilities as the remnant evolves.

The nature of the remnant of a binary neutron star (BNS) merger is uncertain. Although it certainly is a black hole in the cases of the most massive BNSs, X-ray light-curves from short gamma-ray burst afterglows suggest a neutron star (NS) as a viable candidate for the merger remnant and central engine of these transients. When jointly observed with gravitational waves (GWs), X-ray light-curves from BNS merger events could provide critical constraints on the remnant nature. We assess the current and future capabilities for detecting an NS remnant through X-ray observations following GW detections. To this end, we simulated GW signals from BNS mergers and the subsequent X-ray emission from newborn millisecond magnetars. We modeled the GW detectability for the current and next-generation GW interferometers, and we reproduced the X-ray emission using a dedicated numerical code that models magnetar spin-down and ejecta dynamics informed by numerical relativity simulations. In our simulations, 2% - 16% of the BNS mergers form millisecond magnetars. Up to ∼ 70% of these might be detectable, which means up to 1.0 +0.3 _ -0.3 millisecond magnetar detections per year with instruments such as SVOM/MXT during the LIGO Virgo KAGRA LIGO India (LVKI) O5 run. The best detectability occurs about two hours post merger. For next-generation GW interferometers, this rate might increase by up to three orders of magnitude, with the peak detectability three to four hours post merger. We also explored how the magnetar magnetic field strength and observer viewing angle affect detectability, and we discuss optimized observational strategies. Although more likely with upcoming GW interferometers, the detection of the spin-down emission of a millisecond magnetar may already be within reach. This warrants sustained theoretical and observational efforts given the profound implications for mergers, gamma-ray bursts, and NS physics of a single detection.

In abstract visual reasoning, monolithic deep learning models suffer from limited interpretability and generalization, while existing neuro-symbolic approaches fall short in capturing the diversity and systematicity of attribute and relation representations. To address these challenges, we propose a systematic abductive reasoning model with diverse relation representations (Rel-SAR) in vector-symbolic architecture (VSA) to solve Raven's progressive matrices (RPM). To derive attribute representations with symbolic reasoning potential, we introduce not only various types of atomic high-dimensional (HD) encodings that capture numeric, periodic, and logical semantics, but also the structured HD representation (SHDR) for the overall grid component. For systematic reasoning, we further propose novel numerical and logical relation functions and perform rule abduction and execution in a unified framework built upon these relation representations. Experimental results demonstrate that Rel-SAR achieves significant performance improvements on RPM tasks. By synergistically combining HD attribute representations with symbolic reasoning, Rel-SAR enables systematic abductive reasoning with both interpretable and computable semantics.

Abstract Apparent horizon (AH) finders are essential for characterizing black holes and excising their interiors in numerical relativity (NR) simulations. However, open-source AH finders to date are tightly coupled to individual NR codes. We introduce BHaHAHA , the BlackHoles@Home AH Algorithm, the first open-source, infrastructure-agnostic library for AH finding in NR. BHaHAHA implements the first-ever hyperbolic flow-based approach, recasting the elliptic partial differential equation for a marginally outer trapped surface as a damped nonlinear wave equation. To enhance performance, BHaHAHA incorporates a multigrid-inspired refinement strategy, an over-relaxation technique, and OpenMP parallelization. When compared to a naïve hyperbolic relaxation implementation, these enhancements result in 64x speedups for difficult common-horizon finds on a single spacetime slice, enabling BHaHAHA to achieve runtimes within 10% of the widely used (single-core) AHFinderDirect and outperform it on multiple cores. For dynamic horizon tracking with typical core counts on a high-performance-computing cluster, BHaHAHA is approximately 2.1 times faster than AHFinderDirect at accuracies limited by interpolation of metric data from the host NR code. Implemented and tested in both the Einstein Toolkit and BlackHoles@Home , BHaHAHA demonstrates that hyperbolic relaxation can be a robust, versatile, and performant approach for AH finding.

A user-friendly Python interface for the numerical relativity code AMSS-NCKU

Abstract The gravitational collapse of a massless scalar field remains a demanding benchmark for numerical methods in numerical relativity, as it exhibits critical behavior at the boundary between dispersion and black hole formation. In this work, we revisit this problem by relying on Physics-Informed Neural Networks (PINNs) as flexible solvers for partial differential equations, thereby providing a comparative assessment of several recent neural architectures. Building on the Einstein-massless-Klein-Gordon formulation in polar-areal coordinates, we consider four initial-value problems encompassing subcritical, critical, and supercritical regimes and use high-resolution finite-difference simulations as reference solutions. Our study is primarily comparative: we evaluate several state-of-the-art deep learning architectures, including vanilla and high-precision PINNs, sinusoidal-feature and quadratic-residual variants, and Kolmogorov-Arnold Networks, all trained under a common loss design that encodes the field equations, boundary conditions, and causal time-space enforcement, together with a novel adaptive spacetime sampling. Within this framework, we also introduce ModPINN, a modest modification of standard PINNs that augments standard multilayer perceptrons with coordinate embeddings, quadratic layers, and other common ingredients in recent literature. This study shows that deep-learning-based methods can reproduce finite-difference solutions for the scalar field and the spacetime metric with competitive accuracy using significantly fewer collocation points than more traditional methodologies. While no single architecture dominates in all regimes, ModPINN achieves particularly stable and accurate solutions near criticality, indicating that suitably designed embeddings and adaptive sampling can enhance the robustness of PINNs for challenging gravitational-collapse scenarios.

It was recently shown that, in a binary coalescence, the greybody factor of the remnant black hole modulates the postmerger ringdown signal. In this work, we demonstrate that a simple four-parameter model based on the greybody factor accurately reproduces the frequency-domain amplitude of a large set of comparable-mass, aligned-spin numerical relativity waveforms from the SXS catalog, achieving mismatches of order O ( 10 − 5 ) and improving existing models by roughly 2 orders of magnitude. We also identify the optimal initial frequency for applying the model in the frequency domain and provide analytical fits of the model parameters in terms of the progenitor masses and aligned spins. Our results pave the way for new consistency tests of the ringdown phase, complementary to traditional black hole spectroscopy.

Abstract Many HPC applications that solve differential equations rely on the Runge–Kutta (RK) 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 RK 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 RK methods.

Light scalar particles arise naturally in many extensions of the standard model and are compelling dark-matter candidates. Gravitational interactions near black holes can trigger the growth of dense scalar configurations that, if sustained during inspiral, alter binary dynamics and imprint signatures on gravitational-wave signals. Detecting such effects would provide a novel probe of fundamental physics and dark matter. Here, we develop a semianalytic waveform model for binaries in scalar environments, validate it against numerical relativity simulations, and apply it in a Bayesian analysis of the LIGO-Virgo-KAGRA catalog. We obtain physically meaningful upper limits on scalar densities around most compact binaries. For GW190728 and GW190814, vacuum lies outside the 95% credible region. When including superradiance priors, GW190728 shows tentative evidence for a scalar environment with a Bayes factor of ln B_{vac}^{env}≈3.5, consistent with a light scalar of mass ∼10^{-12} eV.

Recent progress in gravitational self-force theory has led to the development of a first postadiabatic (1PA) waveform model for nonspinning, quasicircular compact binaries []. In this paper, we extend that model to allow for a slowly spinning primary black hole and a generic, precessing spin on the secondary object, restricting to the case of small misalignment between the primary spin and the orbital angular momentum. We demonstrate excellent agreement between our waveforms and fully nonlinear numerical relativity simulations for mass ratios q ≳ 5 and primary spins | χ 1 | ≲ 0.1 and arbitrary secondary spin χ 2 ≲ 1 . In particular we present the resummed 1PAT1R waveform model, which significantly improves the accuracy of the original 1PAT1 waveforms for comparable masses and increasing primary spin. Our models are publicly available in the a package.

Abstract Gravitational wave signals from asymmetric binary black hole systems have been shown to exhibit additional chirps beyond the primary merger chirp in the post-merger region of the time-frequency domain. These secondary post-merger chirps correlate to the evolving geometry of the common horizon that forms as the binary merges and were previously studied through numerical relativity simulation in a zero-spin regime. In this work, we investigate the post-merger time-frequency structure in systems with both aligned and precessing spin using widely available waveform models. We find that the inclusion of strong aligned spin (ξ = 0.75) induces further post-merger time-frequency peaks. Additionally we show that even a mild precessing spin (χp = 0.25) strongly affects the distribution of post-merger radiative power across the celestial sky of the final black hole. Our results support the theory of a correlation between the post-merger signal and horizon geometry.

Abstract We model the inspiral and merger dynamics of two co-planar rings in Newtonian mechanics with GR motivated corrections and illustrate their similarity with those of black hole (BH) binary systems on the orbital plane. Our simulation reveals a banana-shape deformation of the ‘BHs’ involved, and a typhoon-like spiral structure in the merger product. Using an eXact One-Body approach, we compute the full gravitational waveform of this process and qualitatively reproduce results consistent with those of numerical relativity (NR). Our simulation offers a transparent link between the feature of gravitational waveforms and the internal structure of BHs, thus a complementary interpretation of physics behind NR.

The last evolutionary stages of massive black hole binaries prior to coalescence are 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, 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; this is 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 a merger with high computational efficiency and accuracy, taking relativistic effects on the gas into account. To validate our implementation, we performed two tests. First, we explored the case of a relativistic Bondi flow around a binary, finding very good agreement with numerical relativity simulations. Then, we explored 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 apsidal precession in the orbiting gas trajectories 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 a merger, even at scales much larger than those currently probed by numerical relativity simulations.

Abstract We describe AthenaK : a new implementation of the Athena++ block-based adaptive mesh refinement framework using the Kokkos programming model. Finite volume methods for Newtonian, special relativistic, and general relativistic (GR) hydrodynamics and magnetohydrodynamics (MHD), and GR-radiation hydrodynamics and MHD, as well as a module for evolving Lagrangian tracer or charged test particles (e.g., cosmic rays) are implemented using the framework. In two companion papers, we describe (1) a new solver for the Einstein equations based on the Z4c formalism, and (2) a GRMHD solver in dynamical spacetimes also implemented using the framework, enabling new applications in numerical relativity. By adopting Kokkos , the code can be run on virtually any hardware, including CPUs, GPUs from multiple vendors, and emerging Advanced RISC Machine processors. AthenaK shows excellent performance and weak scaling, achieving over 1 billion cell updates per second for hydrodynamics in three dimensions on a single NVIDIA Grace Hopper processor. It does this with a typical parallel efficiency of 80% on 65,536 AMD GPUs on the OLCF Frontier system. Such performance portability enables AthenaK to leverage modern exascale computing systems for challenging applications in astrophysical fluid dynamics, numerical relativity, and multimessenger astrophysics.

Interpreting multimessenger signals from neutron stars and black holes requires reliable general relativistic magnetohydrodynamics (GRMHD) simulations across rapidly evolving high-performance computing platforms, yet key algorithms are routinely rewritten within infrastructure-specific numerical relativity codes, hindering verification and reuse. We present the General Relativistic Hydrodynamics Library (GRHayL), a modular, infrastructure-agnostic GR(M)HD library providing conservative-to-primitive recovery, reconstruction, flux/source and induction operators, equations of state, and neutrino leakage through an intuitive interface. GRHayL refactors and extends the mature IllinoisGRMHD code into reusable point- and stencilwise kernels, enabling rapid development and cross-code validation in diverse frameworks, while easing adoption of new microphysics and future accelerators. We implement the same kernels in the Einstein Toolkit (Carpet and Carpetx) and BlackHoles@Home, demonstrating portability with minimal duplication. Validation combines continuous-integration unit tests with cross-infrastructure comparisons of analytic GRMHD Riemann problems, dynamical Tolman-Oppenheimer-Volkoff evolutions, and binary neutron star mergers, showing comparable or improved behavior over legacy IllinoisGRMHD and established Einstein Toolkit codes.

With the upcoming third-generation gravitational-wave detectors comes the need to build complete, faithful, and fast waveform models for asymmetric-mass-ratio compact binaries. Most efforts within the self-force community have focused on modeling these binaries' inspiral regime, but for ground-based detectors the systems' final merger can represent the dominant part of the signal. Recent work by three of us has extended the multiscale self-force framework through the transition-to-plunge and merger-ringdown regimes for nonspinning binaries. In this paper, we generalize the next-to-next-to-leading-order transition-to-plunge waveform model to include the spin of the primary black hole. We also improve the construction of composite inspiral-transition waveform models by performing a change of variables on the binary's mechanical phase space during the transition to plunge. We provide detailed discussions of our numerical implementation and comparisons with numerical relativity simulations.

The Einstein constraint equations form a nonlinear and underdetermined elliptic system whose solution requires specifying appropriate gauge, conformal, and freely chosen geometric data.In this work, we present a hybrid numericalmachine learning framework for solving the Hamiltonian and momentum constraints using a combination of classical solvers, Physics-Informed Neural Networks (PINNs), and Deep Operator Networks (DeepONets).We clarify the mathematical structure of the constraint system, explicitly describe the conformal background and freely chosen components of the initial data, and construct a well-posed elliptic formulation suitable for numerical treatment.To demonstrate feasibility, we implement PINN and DeepONet models for the conformally flat, time-symmetric vacuum constraint and compare the neural solutions with a classical finite-difference reference.The results show that machineassisted PDE solvers can approximate the constraint equations with competitive accuracy while offering mesh-free flexibility.This revised formulation provides a concrete basis for future large-scale simulations in numerical relativity.

Abstract The formation of black holes by the gravitational collapse of stars is known to spontaneously excite particle pairs out of the quantum vacuum. For the canonical vacuum state at past null infinity, the expected number of particles received at future null infinity can be obtained in full closed form at sufficiently late times. However, for intermediate times, or for more complicated astrophysical processes (e.g. binary black hole mergers), the problem is technically challenging and has not yet been resolved. We develop here a numerical approach to study scattering problems of massless quantum fields in asymptotically flat spacetimes, based on the hyperboloidal slice method used in numerical relativity and perturbation theory. This promising approach can reach both past and future null infinities, and therefore it has the potential to address the Hawking scattering problem more rigorously than evolution on the usual Cauchy slices. We test this approach with some dynamical toy models in Minkowski using effective potentials that mimic the effects of gravity, and compute the spectrum of particles received at future null infinity. We finally discuss future prospects for applying this framework in more relevant gravitational scenarios.

We revisit the problem of gravitational-wave extraction in numerical relativity with gauge-invariant metric perturbation theory of spherical spacetimes. Our extraction algorithm allows the computation of even-parity (Zerilli-Moncrief) and odd-parity (Regge-Wheeler) multipoles of the strain from a ( <a:math xmlns:a="http://www.w3.org/1998/Math/MathML" display="inline"> <a:mrow> <a:mn>3</a:mn> <a:mo>+</a:mo> <a:mn>1</a:mn> </a:mrow> </a:math> ) metric without the assumption that the spherical background is in Schwarzschild coordinates. The algorithm is validated with a comprehensive suite of three-dimensional problems including fluid ( <c:math xmlns:c="http://www.w3.org/1998/Math/MathML" display="inline"> <c:mi>f</c:mi> </c:math> -modes) and spacetime ( <e:math xmlns:e="http://www.w3.org/1998/Math/MathML" display="inline"> <e:mi>w</e:mi> </e:math> -modes) perturbations of neutron stars, gravitational collapse of rotating neutron stars, circular binary black holes mergers and black hole dynamical captures, and binary neutron star mergers. We find that metric extraction is robust in all the considered scenarios and delivers waveforms of overall quality similar to curvature (Weyl) extraction. Metric extraction is particularly valuable in identifying waveform systematics for problems in which the reconstruction of the strain from the Weyl multipoles is ambiguous. Direct comparison of different choices for the gauge-invariant master functions show very good agreement in the even-parity sector. Instead, in the odd-parity sector, assuming the background in Schwarzschild coordinates can minimize gauge effects related to the use of the <g:math xmlns:g="http://www.w3.org/1998/Math/MathML" display="inline"> <g:mi mathvariant="normal">Γ</g:mi> </g:math> -driver shift. Moreover, for optimal choices of the extraction radius, a simple extrapolation to null infinity can deliver waveforms compatible to Cauchy-characteristic extrapolated waveforms.

We reinvestigate the stability properties of ultracompact spinning boson stars with a stable light ring using fully nonlinear 3+1 and 2+1 numerical relativity simulations and two different formulations of the Einstein equations. We find no evidence of an instability on timescales of $t μ\sim 10^4$ (in units of the scalar mass), when allowing the star to be perturbed either solely by discretization error or by imposing various types of perturbations to our initial data. We find that the initially imposed perturbations exhibit slow decay, even for magnitudes just below the order where immediate collapse is induced.