Search

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 $\mathcal{O}(10^{-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 $\geq 80\%$ efficiency up to 32768 GPUs and $74\%$ up to 65536 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., Comment: 26 pages, 25 figures; to be submitted to ApJS soon

We present the numerical relativity module within AthenaK, an open source performance-portable astrophysics code designed for exascale computing applications. This module employs the Z4c formulation to solve the Einstein equations. We demonstrate its accuracy through a series of standard numerical relativity tests, including convergence of the gravitational waveform from binary black hole coalescence. Furthermore, we conduct scaling tests on OLCF Frontier and NERSC Perlmutter, where AthenaK exhibits excellent weak scaling efficiency of 80% on up to 65,536 AMD MI250X GPUs on Frontier (relative to 4 GPUs) and strong scaling efficiencies of 84% and 77% on AMD MI250X and NVIDIA A100 GPUs on Frontier and Perlmutter respectively. Additionally, we observe a significant performance boost, with two orders of magnitude speedup ($\gtrsim 200\times$) on a GPU compared to a single CPU core, affirming that AthenaK is well-suited for exascale computing, thereby expanding the potential for breakthroughs in numerical relativity research.

We numerically investigate the imprints of gravitational radiation-reaction driven changes to a black hole's mass and spin on the corresponding ringdown waveform. We do so by comparing the dynamics of a perturbed black hole evolved with the full (nonlinear) versus linearized Einstein equations. As expected, we find that the quasinormal mode amplitudes extracted from nonlinear evolution deviate from their linear counterparts at third order in initial perturbation amplitude. For perturbations leading to a change in the black hole mass and spin of $\sim 5\%$, which is reasonable for a remnant formed in an astrophysical merger, we find that nonlinear distortions to the complex amplitudes of some quasinormal modes can be as large as $\sim 50\%$ at the peak of the waveform. Furthermore, the change in the mass and spin results in a drift in the quasinormal mode frequencies, which for large amplitude perturbations causes the nonlinear waveform to rapidly dephase with respect to its linear counterpart. %These two nonlinear effects together create a large distortion in both the amplitude and phase of the ringdown gravitational waveform. Surprisingly, despite these nonlinear effects creating significant deviations in the nonlinear waveform, we show that a linear quasinormal mode model still performs quite well from close to the peak amplitude onwards.

We investigate quadratic quasinormal mode coupling in black hole spacetime through numerical simulations of single perturbed black holes using both numerical relativity and second-order black hole perturbation theory. Focusing on the dominant $\ell=|m|=2$ quadrupolar modes, we find good agreement (within $\sim10\%$) between these approaches, with discrepancies attributed to truncation error and uncertainties from mode fitting. Our results align with earlier studies extracting the coupling coefficients from select binary black hole merger simulations, showing consistency for the same remnant spins. Notably, the coupling coefficient is insensitive to a diverse range of initial data, including configurations that led to a significant (up to $5\%$) increase in the remnant black hole mass. These findings present opportunities for testing the nonlinear dynamics of general relativity with ground-based gravitational wave observatories. Lastly, we provide evidence of a bifurcation in coupling coefficients between counter-rotating and co-rotating quasinormal modes as black hole spin increases.

The spectroscopic study of black hole quasinormal modes in gravitational-wave ringdown observations is hindered by our ignorance of which modes should dominate astrophysical signals for different binary configurations, limiting tests of general relativity and astrophysics. In this work, we present a description of the quasinormal modes that are excited in the ringdowns of comparable mass, quasi-circular precessing binary black hole coalescences -- a key region of parameter space that has yet to be fully explored within the framework of black hole spectroscopy. We suggest that the remnant perturbation for precessing and non-precessing systems is approximately the same up to a rotation, which implies that the relative amplitudes of the quasinormal modes in both systems are also related by a rotation. We present evidence for this by analyzing an extensive catalog of numerical relativity simulations. Additional structure in the amplitudes is connected to the system's kick velocity and other asymmetries in the orbital dynamics. We find that the ringdowns of precessing systems need not be dominated by the ${(\ell,m)=(2,\pm 2)}$ quasinormal modes, and that instead the $(2,\pm 1)$~or~$(2,0)$ quasinormal modes can dominate. Our results are consistent with a ringdown analysis of the LIGO-Virgo gravitational wave signal GW190521, and may also help in understanding phenomenological inspiral-merger-ringdown waveform model systematics., Comment: Data Release and Analysis Scripts: https://github.com/HengruiPrinceton/precession_ringdown

The intricacies of black hole ringdown analysis are amplified by the absence of a complete set of orthogonal basis functions for quasinormal modes. Although damped sinusoids effectively fit the ringdown signals from binary black hole mergers, the risk of overfitting remains, due to initial transients and nonlinear effects. In light of this challenge, we introduce two methods for extracting quasinormal modes in numerical simulations and qualitatively study how the transient might affect quasinormal mode fitting. In one method, we accurately fit quasinormal modes by using their spatial functional form at constant time hypersurfaces, while in the other method, we exploit both spatial and temporal aspects of the quasinormal modes. Both fitting methods leverage the spatial behavior of quasinormal eigenfunctions to enhance accuracy, outperforming conventional time-only fitting techniques at null infinity. We also show that we can construct an inner product for which the quasinormal eigenfunctions form an orthonormal (but not complete) set. We then conduct numerical experiments involving linearly perturbed Kerr black holes in horizon penetrating, hyperboloidally compactified coordinates, as this setup enables a more precise isolation and examination of the ringdown phenomenon. From solutions to the Teukolsky equation, describing scattering of an ingoing gravitational wave pulse, we find that the contributions from early-time transients can lead to large uncertainties in the fit to the amplitudes of higher overtones ($n\geq 3$). While the methods we discuss here cannot be applied directly to data from merger observations, our findings underscore the persistence of ambiguities in interpreting ringdown signals, even with access to both temporal and spatial information., Comment: 13+6 pages, 12 figures. V2: Minor changes to match the published version

Recent interferometric observations by the Event Horizon Telescope have resolved the horizon-scale emission from sources in the vicinity of nearby supermassive black holes. Future space-based interferometers promise to measure the "photon ring"--a narrow, ring-shaped, lensed feature predicted by general relativity, but not yet observed--and thereby open a new window into strong gravity. Here we present AART: an Adaptive Analytical Ray-Tracing code that exploits the integrability of light propagation in the Kerr spacetime to rapidly compute high-resolution simulated black hole images, together with the corresponding radio visibility accessible on very long space-ground baselines. The code samples images on a nonuniform adaptive grid that is specially tailored to the lensing behavior of the Kerr geometry and is therefore particularly well-suited to studying photon rings. This numerical approach guarantees that interferometric signatures are correctly computed on long baselines, and the modularity of the code allows for detailed studies of equatorial sources with complex emission profiles and time variability. To demonstrate its capabilities, we use AART to simulate a black hole movie of a stochastic, non-stationary, non-axisymmetric equatorial source; by time-averaging the visibility amplitude of each snapshot, we are able to extract the projected diameter of the photon ring and recover the shape predicted by general relativity., Comment: 39 pages, 16 figures. Updated to match published version. V4: corrected two very minor typos (in Eqs. 50 and 70)

In extreme scattering events, the brightness of a compact radio source drops significantly, as light is refracted out of the line of sight by foreground plasma lenses. Despite recent efforts, the nature of these lenses has remained a puzzle, because any roughly round lens would be so highly overpressurized relative to the interstellar medium that it could only exist for about a year. This, combined with a lack of constraints on distances and velocities, has led to a plethora of theoretical models. We present observations of a dramatic double-lensing event in pulsar PSR~B0834+06 and use a novel phase-retrieval technique to show that the data can be reproduced remarkably well with a two-screen model: one screen with many small lenses and another with a single, strong one. We further show that the latter lens is so strong that it would inevitably cause extreme scattering events. Our observations show that the lens moves slowly and is highly elongated on the sky. If similarly elongated along the line of sight, as would arise naturally from a sheet of plasma viewed nearly edge-on, no large over-pressure is required and hence the lens could be long-lived. Our new technique opens up the possibility of probing interstellar plasma structures in detail, leading to understanding crucial for high-precision pulsar timing and the subsequent detection of gravitational waves.

The interstellar medium hosts a population of scattering screens, most of unknown origin. Scintillation studies of pulsars provides a sensitive tool for resolving these scattering screens and a means of measuring their properties. In this paper, we report our analysis of 34 years of Arecibo observations of PSR B1133+16, from which we have obtained high-quality dynamic spectra and their associated scintillation arcs, arising from the scattering screens located along the line of sight to the pulsar. We have identified six individual scattering screens that are responsible for the observed scintillation arcs, which persist for decades. Using the assumption that the scattering screens have not changed significantly in this time, we have modeled the variations in arc curvature throughout the Earth's orbit and extracted information about the placement, orientation, and velocity of five of the six screens, with the highest-precision distance measurement placing a screen at just $5.46^{+0.54}_{-0.59}$ pc from the Earth. We associate the more distant of these screens with an under-dense region of the Local Bubble., Comment: 16 pages, 9 figures, accepted for publication in ApJ

An important result by Agmon implies that an eigenfunction of a Schr\"{o}dinger operator in $\mathbb{R}^n$ with eigenvalue $E$ below the bottom of the essential spectrum decays exponentially if the associated classically allowed region $\{x \in \mathbb{R}^n~:~ V(x) \leq E \}$ is compact. We extend this result to a class of Schr\"{o}dinger operators with eigenvalues, for which the classically allowed region is not necessarily compactly supported: We show that integrability of the characteristic function of the classically allowed region with respect to an increasing weight function of bounded logarithmic derivative leads to $L^2$-decay of the eigenfunction with respect to the same weight. Here, the decay is measured in the Agmon metric, which takes into account anisotropies of the potential. In particular, for a power law (or, respectively, exponential) weight, our main result implies that power law (or, respectively, exponential) decay of "the size of the classically allowed region" allows to conclude power law (or, respectively, exponential) decay, in the Agmon metric, of the eigenfunction.

Radio pulsar signals are significantly perturbed by their propagation through the ionized interstellar medium. In addition to the frequency-dependent pulse times of arrival due to dispersion, pulse shapes are also distorted and shifted, having been scattered by the inhomogeneous interstellar plasma, affecting pulse arrival times. Understanding the degree to which scattering affects pulsar timing is important for gravitational wave detection with pulsar timing arrays (PTAs), which depend on the reliability of pulsars as stable clocks with an uncertainty of ~100ns or less over ~10yr or more. Scattering can be described as a convolution of the intrinsic pulse shape with an impulse response function representing the effects of multipath propagation. In previous studies, the technique of cyclic spectroscopy has been applied to pulsar signals to deconvolve the effects of scattering from the original emitted signals. We present an analysis of simulated data to test the quality of deconvolution using cyclic spectroscopy over a range of parameters characterizing interstellar scattering and pulsar signal-to-noise ratio. We show that cyclic spectroscopy is most effective for high-S/N and/or highly scattered pulsars. We conclude that cyclic spectroscopy could play an important role in scattering correction to distant populations of highly scattered pulsars not currently included in PTAs. For future telescopes and for current instruments such as the Green Bank Telescope upgraded with the ultrawide bandwidth (UWB) receiver, cyclic spectroscopy could potentially double the number of PTA-quality pulsars., Comment: 21 pages, 10 figures. Typo fixed, conclusions unchanged. Accepted for publication in ApJ