Adaptive mesh refinement with subcycling in time enables different grid levels to advance using their own time steps, ensuring finer grids employ smaller steps for accuracy while coarser grids take larger steps to improve computational efficiency. We present the development, validation, and performance analysis of a subcycling in time algorithm implemented within the carpetx driver in the einstein toolkit framework. This new approach significantly improves upon the previous subcycling implementation in the carpet driver by achieving higher-order convergence—fourth order in time instead of second order—and enhanced scaling performance. The key innovation lies in optimizing the exchange of ghost points at refinement boundaries, limiting it to the same number as those at interprocess boundaries using dense output from coarser levels, thereby reducing computational and communication overhead compared to the implementation in carpet, which required a larger number of buffer zones. To validate the algorithm, we first demonstrate its fourth-order convergence using a scalar wave test. We then apply the algorithm to binary black hole simulations, confirming its robustness and accuracy in a realistic astrophysical scenario. The results show excellent agreement with the well-established lazev code. Scaling tests on CPU (Frontera) and GPU (Vista) clusters reveal significant performance gains, with the new implementation achieving improved speed and scalability compared to the carpet-based version.

We present a novel systematic framework for engineering energy-momentum tensors in linearized gravity that generate gravitational fields with a prescribed multipolar structure, aimed at modeling black hole mimickers. As a concrete example, we analyze an anisotropic rotating fluid with a Gaussian-like energy-density profile which, at linearized order, satisfies all energy and causality conditions while generating a gravitational field that exactly reproduces the multipolar structure of a Kerr black hole. We then investigate the phenomenology of both the source and the induced metric, comparing our findings with the corresponding black hole limit.

The entangling properties of models of classical gravity interacting with quantum matter (i.e., hybrid models of gravity) are investigated in the context of the experimental proposals to detect gravitationally induced entanglement. We prove that entanglement generation can indeed take place within these models, and characterize it quantitatively. We identify the root of entanglement generation in these models in the presence of some underlying nonlocality of the theories. Nonetheless, by focusing in particular on the Diósi-Penrose model for two gravitationally interacting masses, we show that entanglement-based experiments have the potential to either falsify the model entirely or constrain the free parameter of the model $R_0$ up to values 6 orders of magnitude above the current state of the art.

We evolve for the first time in full general relativity a small, collisional $N$-body black hole (BH) cluster of arbitrary total mass $M$. The bound cluster is initially compact (radius $R/M\sim 10$), stable, and consists of 25 equal-mass, nonspinning BHs. The dynamical interactions of compact objects in $N$-body clusters is of great interest for the formation of BHs in the upper mass gap as well as intermediate and supermassive BHs. These are potential sources of gravitational waves that may be detected by both current and future observatories. Unlike previous $N$-body Newtonian and post-Newtonian simulations, no “subgrid physics” is required to handle collisions and mergers. We can therefore confirm in full general relativity several predictions from these simulations and analytic estimates: the runaway growth of a large BH via repeated mergers; spindown of the central BH with increasing captures; the ejection of a BH with a large asymptotic velocity due to a several-body interaction; and a regime where mergers occur primarily via direct collisions on highly eccentric orbits instead of quasicircular inspirals. We extract the gravitational wave signal and find it has several distinct features associated with the compact cluster regime. Our results suggest the signal is sufficiently loud that next generation observatories would likely be able to detect similar events across most of the observable universe. This work is a preliminary proof-of-principle study that we hope will open up a new arena for numerical relativity and the study of $N$-body compact systems.

The Laser Interferometer Space Antenna (LISA) will detect mHz gravitational waves from many astrophysical sources, including millions of compact binaries in the Galaxy, thousands of which may be individually resolvable. The large number of signals overlapping in the LISA dataset requires a global fit in which an unknown number of sources are modeled simultaneously. This introduces a label-switching ambiguity for sources in the same class, making it challenging to distill a traditional astronomical catalog from global-fit posteriors. We present a method to construct a catalog by post-processing the global-fit posterior, relabeling samples to minimize the statistical divergence between the global fit and a factorized catalog representation. The resulting catalog consists of the source posterior distributions and their probabilities of having an astrophysical origin. We demonstrate our algorithm on two toy models and on a small simulated LISA dataset of Galactic binaries. Our method is implemented in the open-source python package petra_catalogs, and it can be applied in postprocessing to the output of any global-fit sampler.

Recent work has highlighted the potentially detectable gravitational-lensing effect of neutrino halos on cosmic-microwave-background (CMB) fluctuations with upcoming instruments like SO, CMB-S4, and CMB-HD. Accurate modeling of neutrino-halo density profiles are essential for making theory predictions of their cosmological effects. Yet, they are computationally intensive, particularly in the nonlinear regime. In this work, we present an efficient numerical framework for computing neutrino profiles based on $N-1$–body simulations within a flat Friedmann-Robertson-Walker Universe. Our approach enables highly parallelized, rapid calculations of neutrino trajectories near spherically symmetric dark-matter halos, delivering results within seconds. In addition to neutrinos, we demonstrate an application to model the clustering of other noncold relics, such as a freeze-in dark-matter component. The framework is flexible in its definitions of cosmological and dark-matter-halo parameters, which can be particularly valuable for rapid-scanning tasks. It can also seamlessly incorporate new physics, as we demonstrate with examples of a time-varying gravitational constant and nonstandard neutrino phase-space distributions.

Quasiuniversal relations are known to exist among various neutron star observables that do not depend sensitively on the underlying nuclear matter equations of state. For example, some of these relations imply that the tidally induced multipole moments are approximately characterized by the electric-type quadrupolar tidal deformability. Such relations can be used to reduce the number of independent tidal parameters in gravitational-waveform modeling, thereby allowing us to infer extreme nuclear matter properties more accurately and test general relativity in an insensitive manner to uncertainties in nuclear physics. We present a comprehensive theoretical investigation into approximate universality of neutron stars. Our approach employs a semianalytic relativistic stellar interior model, which extends the Tolman VII solution, thereby enabling a refined exploration of the tidal properties of nonrotating stars within a semianalytic framework. The derived power-law relations among various tidal deformabilities—referred to as the universal Love relations—agree well with expressions in previous work found empirically. We elucidate how the equation-of-state dependence is suppressed in a particular combination of macroscopic physical parameters induced by perturbations and demonstrate that the relation between the moment of inertia and electric-type quadrupolar tidal deformability (I-Love relation) rests on the same underlying mechanism. Our findings indicate that the approximate universality of neutron stars can be attributed to low compressibility, consistent with some of the previous studies on the possible origin of the universality.

A novel methodology to obtain global transonic solutions around compact objects is reported here. A unified methodology to obtain accretion as well as wind solutions around these objects has been presented. Flows around compact objects are dissipative, and the conservation equations are therefore stiff. In such conditions, obtaining of sonic point(s) and hence, the transonic solution is not trivial. The conserved equations of motion fail to integrate in the presence of realistic viscosity, thereby making it difficult to obtain a global solution. This inhibits one from getting an actual picture of an astrophysical flow. The current work addresses this longstanding issue of obtaining solutions for both accretion and wind. The methodology developed utilizes the inner boundary conditions and takes recourse to Implicit-Explicit (ImEx) integration schemes, to obtain general global transonic accretion and wind solutions. This is the first time such an attempt has been made. Current work considers the different cooling processes like bremsstrahlung, synchrotron and their inverse-Comptonizations, which are found to affect the thermodynamics of the flow. This methodology could successfully generate all topologies of global solutions, multiple sonic point regime, as well as shocks. A broad parameter space study has been done in this work. In an upcoming part II of the paper, a detailed discussion on the spectra and luminosity of the accretion and wind solutions has been presented.

We explore the Faddeev-Jackiw symplectic Hamiltonian reduction of the $sl(2)$ affine Toda model coupled to matter, which includes new parametrizations for a scalar field and a Grassmannian fermionic field. The structure of constraints and symplectic potentials primarily dictates the strong-weak dual coupling sectors of the theory, ensuring the equivalence between the Noether and topological currents. The analytical calculations encompass the fermion-kink classical solution, the excited fermion bound states localized on the kink, and the scattering states, all of which account for the fermion backreaction on the soliton. The total energy, which includes the classical fermion-soliton interaction energy, the bound-state fermion energy, and the fermion vacuum polarization energy, is determined by the topological charge of the kink. This system satisfies first-order differential equations and a chiral current conservation equation. Our results demonstrate that the excited fermion bound states and scattering states significantly alter the properties of the kink. Notably, they give rise to a pumping mechanism for the topological charge of the in-gap kink due to fermionic backreaction, as well as the appearance of kink states in the continuum.

We study a variant of the three-Higgs-doublet model (3HDM), referred to as the Branco-Grimus-Lavoura–(BGL-)3HDM, incorporating a $U(1{)}_1\times U(1{)}_2$ symmetry, which can distinguish the primary sources of mass for different fermion generations. In the version considered here, the Yukawa matrices in the down-quark and charged lepton sectors are diagonal, thereby eliminating tree-level flavor changing neutral currents (FCNCs) in these sectors. FCNC interactions mediated by neutral nonstandard Higgses are confined to the up-quark sector only. No new Beyond-Standard-Model parameters are introduced by the Yukawa sector of the model, making it as economical as the natural flavor conservation (NFC) versions of 3HDM with a $U(1{)}_1\times U(1{)}_2$ symmetry in terms of the number of free parameters. However, even in the down-quark and in the charged lepton sectors, flavor diagonal but nonuniversal Higgs couplings set this model apart from the NFC versions of the 3HDM.

The recent observation of a diphoton excess around 95 GeV by the CMS and ATLAS Collaborations, along with the $b\overline{b}$ excess reported by the LEP Collaboration in the same mass region, has drawn significant interest in the possibility of new physics beyond the Standard Model (SM). The Georgi-Machacek (GM) model, which extends the Higgs sector of the SM by introducing additional triplet scalars while preserving custodial symmetry at tree level, provides a compelling framework to explain both excesses simultaneously via a light custodial singlet Higgs. In this work, we investigate whether the GM model can still accommodate these excesses when taking into account newly proposed vacuum stability constraints, particularly the positive definiteness conditions. Our numerical analysis not only confirms the existence of a viable parameter space capable of explaining the 95 GeV excesses, but also demonstrates that, compared to traditional tree-level constraints at the electroweak scale, the positive definiteness conditions further expand the allowed parameter space, thereby enhancing the viability of the GM model. Furthermore, we emphasize that future collider experiments will play a crucial role in testing this interpretation by refining Higgs coupling measurements and searching for additional Higgs bosons.

The dynamical diquark model describes multiquark exotic hadrons in terms of diquark components nucleated by heavy quarks and successfully explains multiple features of hidden-charm and hidden-bottom exotics. Here, we apply the model to the marginally heavy case of hidden-strange states to probe whether mesons near 2 GeV with peculiar properties, such as $\varphi (2170)$, $f_2(2340)$, and $X(2370)$, are possible tetraquark candidates. We calculate spin-multiplet average masses using potentials obtained through lattice simulations and quark models, and we also describe the detailed spectra of the expected multiplets as a diagnostic to discern the nature of future hadrons likely to be discovered in this mass region by experiments at facilities such as BESIII, JLab, and the EIC.

Inspired by the BESIII observation of exotic hadron $G(3900)$ in the $e^{+}{e}^{-}\to D\overline{D}$ process [M. Ablikim et al. (BESIII Collaboration), Phys. Rev. Lett. 133, 081901 (2024)], we use the parton and hadron cascade model PACIAE to simulate the productions of charmonium-like exotic hadrons, including $X(3872)$, $Z_c(3900{)}^0$, and $G(3900)$ in the $e^{+}{e}^{-}$ annihilations at $\sqrt{s}=4.95\mathrm{GeV}$. The charmonium-like candidates are recombined by the dynamically constrained phase-space coalescence model using component mesons $D\overline{D}$ [for $G(3900)$] or $D{\overline{D}}^{*}/\overline{D}{D}^{*}$ [for $X(3872)$, $Z_c(3900{)}^0$, and $G(3900)$] in the PACIAE simulated final hadronic state. We then calculate, for the first time, the charmonium-like candidate’s orbital angular momentum quantum number in its rest frame and perform the spectral classification for each of the above candidates. For $G(3900)$, as its $J^{PC}$ is $1^{--}$, it is identified as the $P$-wave $D\overline{D}$ or $D{\overline{D}}^{*}/\overline{D}{D}^{*}$ state. Meanwhile, for $X(3872)$ and $Z_c(3900{)}^0$, as their $J^{PC}\mathrm{s}$ are, respectively, $1^{++}$ and $1^{+-}$, they are identified as the $S$-wave $D{\overline{D}}^{*}/\overline{D}{D}^{*}$ candidates. It is observed that the yield of the $P$-wave $G(3900)$ composed of $D{\overline{D}}^{*}/\overline{D}{D}^{*}$ is significantly larger than that composed of $D\overline{D}$. Moreover, in the $D{\overline{D}}^{*}/\overline{D}{D}^{*}$ channel, the yield of the $P$-wave $G(3900)$ is around three times as large as the yield of the $S$-wave $Z_c(3900{)}^0$, and it is around two orders of magnitude as large as the yield of the $S$-wave $X(3872)$. Finally, significant discrepancies are observed in the transverse momentum spectra and the rapidity distributions among the $X(3872)$, $Z_c(3900{)}^0$, and $G(3900)$ states. These discrepancies are proposed as valuable criteria for identifying the different charmonium-like exotic hadrons from each other.

We construct light cone sum rules (LCSR) for the $B\to K^{*}$ form factors in the large recoil region using vacuum-to-$B$-meson correlation functions, and systematically calculate subleading-power corrections to these form factors at tree level, including next-to-leading power contributions from the hard-collinear propagator, the subleading effective current $\overline{q}\mathrm{\Gamma }[i{\cancel{D}}_{\perp }/(2m_b)]h_v$, and twist-five/six four-particle higher-twist effects. By incorporating the available leading-power results at $\mathcal{O}({\alpha }_s)$ and the corrections to higher-twist $B$-meson light cone distribution amplitudes from our previous work, we improve the precision of theoretical predictions for $B\to K^{*}$ form factors and find that the subleading-power contributions amount to 30% of the corresponding leading-power results. Employing the Bourrely-Caprini-Lellouch parametrization, we determine the numerical results for $B\to K^{*}$ form factors across the full kinematic range through a combined fit of LCSR predictions in the large recoil region and lattice QCD results in the small recoil region. Using the newly obtained $B\to K^{*}$ form factors, we compute the branching fractions for the rare decays $B\to K^{*}{\nu }_{\ell }{\overline{\nu }}_{\ell }$ in the Standard Model, obtaining $\mathcal{B}\mathcal{R}({\overline{B}}^0\to {\overline{K}}^{*0}{\nu }_{\ell }{\overline{\nu }}_{\ell })=8.09(96)\times {10}^{-6}$ and $\mathcal{B}\mathcal{R}({\overline{B}}^{+}\to {\overline{K}}^{*+}{\nu }_{\ell }{\overline{\nu }}_{\ell })=9.95(1.05)\times {10}^{-6}$. Additionally, we predict that the longitudinal $K^{*}$ polarization fraction is $F_L=0.44(4)$.

We measure the time-integrated $CP$ asymmetry in $D^0\to K_{\mathrm{S}}^0{K}_{\mathrm{S}}^0$ decays reconstructed in $e^{+}{e}^{-}\to c\overline{c}$ events collected by the Belle and Belle II experiments. The corresponding data samples have integrated luminosities of 980 and $428{\mathrm{fb}}^{-1}$, respectively. To infer the flavor of the $D^0$ meson, we exploit the correlation between the flavor of the reconstructed decay and the electric charges of particles reconstructed in the rest of the $e^{+}{e}^{-}\to c\overline{c}$ event. This results in a sample which is independent from any other previously used at Belle or Belle II. The result, $A_{CP}(D^0\to K_{\mathrm{S}}^0{K}_{\mathrm{S}}^0)=(1.3\pm 2.0\pm 0.2)\%$, where the first uncertainty is statistical and the second systematic, is consistent with previous determinations and with $CP$ symmetry.

We develop the pole-skipping structure in de Sitter (dS) spacetime and find that their leading frequencies satisfy the relation ${\omega }_{dS}=i2\pi T_{dS}(1-s)$, where $T_{dS}=1/2\pi L$ and $s$ denotes spin. In the two-dimensional dS spacetime, the pole-skipping points near the cosmic horizon $r=L$ for the scalar field of spin 0 and the fermionic field of spin $\frac{1}{2}$ correspond one to one with those in the classical limit as $\lambda \to 0$ in double-scaled Sachdev-Ye-Kitaev model when the temperature is infinite (${\mathrm{DSSYK}}_{\infty }$). This provides a numerical correspondence between quantum gravity in the static patch of two-dimensional dS spacetime and a one-dimensional quantum system, which we consider as a realization of the DS/dS correspondence.

We use Ising-like models to probe the thermal nature of Euclidean spacetime backgrounds. We determine which properties of the background—curvature, the presence of a horizon, or temperature—play a role in phase transitions. The geometries we use are Euclidean Schwarzschild, Rindler, extremal Reissner-Nördstrom (ERN), anti–de Sitter (AdS), and de Sitter. Among these, Rindler is flat, AdS does not have a horizon, and both AdS and ERN have zero temperatures. We find second-order phase transitions as the metric parameter is varied in all cases except for Rindler. Specifically, we find that the transition from order to disorder occurs as the curvature—or Euclidean gravity—increases. This supports our conjecture that Euclidean gravity is an essential ingredient for these phase transitions, as opposed to the presence of a horizon or temperature. Separately, since the selected geometries are position-dependent, the Ising-like models constructed on them are inhomogeneous, whereby they generalize the standard Ising model. We find that a consequence of this is that criticality does not correspond to maximal correlation lengths and scale invariance.

This study investigates polarized images of an equatorial synchrotron-emitting ring surrounding hairy black holes within the Einstein-Maxwell-scalar theory. Our analysis demonstrates qualitative similarities between the polarization patterns of hairy black holes and Schwarzschild black holes. However, due to the nonminimal coupling between the scalar and electromagnetic fields, an increase in black hole charge and coupling constant can substantially amplify polarization intensity and induce deviations in the electric vector position angle. These effects may offer observational signatures to distinguish hairy black holes from Schwarzschild black holes.

The regularity of black-hole solutions, embedded in an expanding universe, is studied in a subclass of Horndeski theories, namely the sum of the simplest quadratic, cubic, and quintic actions. We find that in the presence of a time derivative of the scalar field, driven by the cosmological expansion, this regularity generically imposes large scalar charges for black holes, even when assuming strictly no direct coupling of matter to the scalar field. Such charges cause a significant accretion of the scalar field by the black holes, driving its local time derivative to a small value. This phenomenon, together with the Vainshtein screening typical of these theories, strongly suppresses observable scalar effects. We show that this full class of models is consistent with LIGO/Virgo detections of gravitational waves, but that the LISA mission should be able to constrain the coefficient of the quintic term at the ${10}^{-30}$ level in a self-acceleration scenario, an improvement by 16 orders of magnitude with respect to what is imposed by the speed of gravitational waves.

We investigate the gravitational wave emission from extreme mass ratio inspirals, key targets for the upcoming space-based detector LISA, considering the scenario where the lighter black hole in the binary is endowed by a long-lived, time-dependent scalar field configuration, known as a scalar wig. We develop a formalism to compute scalar perturbations for extreme mass ratio inspirals on circular orbits around Schwarzschild and Kerr black holes, and apply this framework to compute additional fluxes induced by the scalar wig as well as their dependence on the scalar field properties. Our computation provides strong indications that in this scenario the presence of the scalar field does not significantly affect the orbital motion and the gravitational waveform.