Frontiers in Numerical Relativity 2022 (FNR2022)

Both for studies of cosmic censorship and for practical purposes in gravitational wave astronomy, it is desirable to include future null-infinity in the computational domain. Extending formulations of general relativity known to behave well in the strong-field regime out to infinity with compactification is, however, a subtle game. In my presentation I will explain how the competition between decay of fields near infinity and growth of coefficients (due to compactification) plays out in dual-frame generalized harmonic gauge. I will discuss the numerical implementation of the resulting PDEs.

I will review the use of numerical relativity to probe cosmological spacetimes, highlighting the particular challenges faced in such simulations compared to the asymptotically flat case associated with black holes. I will focus on its potential to inform the "inflation versus bounce" debate, as an illustration of these challenges.

Gravitational wave radiation, our window for probing the strong field and dynamical regime of gravity, is unambiguously defined only at future null infinity - the location in spacetime where light rays arrive and thus where signals and global properties of spacetimes can be measured. A convenient way to reach it within numerical relativity simulations is by means of hyperboloidal foliations, smooth spacelike slices that reach future null infinity. Among the approaches to the hyperboloidal problem, here I will focus on conformal compactification, based on an idea by Nobel-laureate Roger Penrose. Our implementation uses the BSSN and Z4 formulations of the Einstein equations and has provided some very promising numerical results of a massless scalar field coupled to gravity in spherical symmetry. I will give an update on the current ongoing work towards a 3D generalization, as the final goal of this work is to provide a far-field numerical framework that includes null infinity for simulations of compact object mergers with accurate gravitational wave extraction.

One of the challenges in numerical relativity is to include future null infinity in the computational domain with a well-posed formulation. Success will not only enable us to evolve any system of astrophysical interest, e.g. binary black holes and extracting the gravitational wave signal at future null infinity, with any desired accuracy, but also help in studying various phenomena of fundamental interest. One proposal is to use hyperboloidal slices. In this talk, I will present our ongoing efforts for obtaining a well-posed formulation of the Einstein Field Equations on hyperboloidal slices, all in spherical symmetry. The natural extension will be to generalize these methods to full 3d.

The Einstein and Maxwell equations are both systems of hyperbolic equations which need to satisfy a set of elliptic constraints throughout evolution. However, while electrodynamics and magnetohydrodynamics have benefited from a large number of evolution schemes that are able to enforce these constraints and are easily applicable to curvilinear coordinates, unstructured meshes, or N-body simulations, many of these techniques cannot be straightforwardly applied to existing formulations of the Einstein equations. We develop a 3+1 a formulation of the Einstein equations which shows a striking resemblance to the equations of relativistic MHD and to EM in material media. The fundamental variables of this formulation are the frame fields, their exterior derivatives, and the Nester-Witten and Sparling forms. These mirror the roles of the electromagnetic 4-potential, the electromagnetic field strengths, the field excitations and the electric current. The role of the lapse function and shift vector, corresponds exactly to that of the scalar electric potential. The formulation is manifestly first order and flux-conservative, which makes it suitable for high-resolution shock capturing schemes and finite-element methods. Being derived as a system of equations in exterior derivatives, it is directly applicable to any coordinate system and to unstructured meshes, and leads to a natural discretisation potentially suitable for the use of machine-precision constraint propagation techniques such as the Yee algorithm and constrained transport. Due to these properties, we expect this new formulation to be beneficial in simulations of many astrophysical systems, such as binary compact objects and core-collapse supernovae as well as cosmological simulations of the early universe.

In this talk I will present the equations of the Einstein-scalar-Gauss-Bonnet theory, a modified theory of gravity, in the 3+1 formulation with a modified gauge that proves to be strongly hyperbolic. Then I'll show some of the numerical results that we have obtained from the implementation of these equations with GRChombo, a numerical general relativity code with fully adaptive mesh refinement.

The study of vacuum critical collapse has acquired an aura of unapproachability in its 30 years of existence. It is often thought that only highly sophisticated numerical methods together with immense supercomputer clusters will suffice for treating this problem. In our recent work we refute this folk lore by using only open-source code running on single machines to discover universal echoing in near-critical axisymmetric vacuum spacetimes. In this talk we discuss the technical aspects of our work, focusing on the implementation of our "quasi-maximal" slicing condition and tuning it to make the most of our limited computing resources.

We use our adapted pseudo-spectral code bamps, with its new hp adaptive mesh refinement, to tune close to the barrier between gravitational collapse and dispersed fields, in order to study the critical phenomena that emerges near that threshold. In spherical symmetry, we observe critical phenomena by evolving massless scalar fields. This has allowed us to assess the adequate choice of gauge source functions and damping parameters for collapse evolutions performed with generalised harmonic gauge. On the other hand in axisymmetry, we evolve six one-parameter families of Brill wave initial data in vacuum: three prolate and three oblate, including two centred and four off-centred. In agreement with previous work we find two disjoint centres of collapse close to the threshold. We also found no evidence for strict DSS nor universal power-law scaling of curvature quantities. Nevertheless, approximately universal strong curvature features and echoes were found within all six families close to the threshold.

We studied the asymptotic behavior of perturbed Kerr initial data by solving the evolutionary form of the vacuum Einstein constraint equations. Unlike in the elliptic formulation, solving the constraints in their evolutionary form, we have direct control over the constrained data only on a single 2-surface. This immediately raises the question of whether it can guarantee asymptotic flatness of initial data. Previous works based on perturbed Schwarzschild initial data have already demonstrated that it is possible to guarantee the asymptotic flatness of solutions. Integrating the equations numerically, we found that results on the modified parabolic-hyperbolic system can be generalized to perturbed Kerr initial data. We have also investigated a new, more robust method for the algebraic-hyperbolic system. Part of this work is funded by NSF CAREER Award PHY–2047382 and NKFIH Grant No. K-115434.

Detection and parameter inference of gravitational-wave signals relies on the comparison of the incoming detector strain data to waveform templates for the gravitational-wave strain $h(t)$ which ultimately rely on the resolution of Einstein's equations via numerical relativity simulations. These, however, commonly output a quantity known as the Newman-Penrose scalar $\psi_4(t)$, related to the strain by $\psi_4(t)=\mathrm{d}^2h(t) / \mathrm{d}t^2$. Therefore, obtaining strain templates involves a double time-integration that introduces artefacts that need to be eased in a rather manual way. By taking second-order finite differences on the detector data and inferring the corresponding noise distribution, we develop a framework to perform gravitational-wave data analysis directly using $\psi_4(t)$ templates. I will first demonstrate this formalism by recovering numerically simulated signals from head-on collisions of Proca stars injected in Advanced LIGO noise. Next, I will show a re-analysis of GW190521 under the Proca-star merger scenario. Our framework removes the need to obtain the strain from numerical relativity simulations therefore avoiding the associated systematic errors.

Compact binary coalescence waveforms generated by numerical relativity codes and associated phenomenological models are a critical part of gravitational-wave searches in detector data. However, the computational cost of the conventional matched-filtering searches is steeply rising as increased sensitivity of detectors expands the waveform parameter space. Lately, machine-learning based search algorithms have gained popularity as a potential efficient complement to the conventional searches. The MLGWSC-1 has been created in order to build a unified testing network for these algorithms, which rely on NR and associated software. This contribution presents the idea, organization and results of the challenge.

The formation of Primordial black holes is naturally enhanced during the QCD phase transition, because of the softening of the equation of state: at a scale between 1 and 3 solar masses, the threshold is reduced of about 10% with a corresponding abundance of primordial black holes increased by more than 100 times. Such black holes could be an interesting source of gravitational waves emitted during black hole mergers, detected by LIGO/VIRGO.

In this talk, I will present the recent progress of numerical modeling of compact binary mergers, which involves at least one neutron star in numerical relativity. Our group develops self-consistent modelings by performing a long-term numerical relativity neutrino-radiation magnetohydrodynamic simulation. I will present the details of the simulation results and perspective.

Numerical-relativity simulations are prerequisites for a reliable interpretation of multi-messenger events such as binary neutron star or black hole neutron star mergers. When using simulation results to interpret observational data, it is of uttermost importance to also ensure a proper discussion of the uncertainties of the simulations. Keeping this in mind, we show some of our most recent results and identify how simulation uncertainties could affect the development of gravitational-wave and electromagnetic models. Finally, we will discuss some of the most recent developments of our numerical-relativity code BAM and review its capability and shortcomings in simulating compact binary mergers.

We discuss the postmerger stage of binary neutron star coalescences. We present an analytic model of postmerger gravitational-wave emission, which achieves an overall good description of the gravitational-wave signal. The physical parameters of the model are useful to understand the dynamics of the system, and we identify new mechanisms shaping the gravitational-wave spectrum. Moreover, we address the deviations of individual data points in quasi-universal relations of postmerger gravitational-wave features encoding the equation of state dependence. We also discuss the threshold for prompt black-hole formation and the impact of strong phase transitions on postmerger properties.

GW190425 was the second gravitational wave (GW) signal compatible with a binary neutron star (BNS) merger detected by the Advanced LIGO and Advanced Virgo detectors. Despite intense follow-up campaigns, no electromagnetic counterpart was identified. Whether the associated kilonova was too dim or the localisation area too broad is still an open question. We simulate 28 BNS mergers with the chirp mass of GW190425 and spanning a mass ratio 1 ≤ q ≤ 1.67, using numerical-relativity simulations with finite temperature, composition dependent nuclear equation of state (EOS) and neutrino radiation. The energy emitted in GWs is less than 0.083 Mo 𝑐^2, where Mo is the solar mass, with a peak luminosity of 1.1 − 2.4 × 10 58 erg s^−1 /(1 + q)^2 . Dynamical ejecta and disc mass are relatively small, the former ranging between 5 × 10^−6 and 10^−3 solar masses and the latter between 10^−5 and 0.1 solar masses . Asymmetric mergers, especially in the case of stiff EOSs, are able to unbind more matter and to form heavier discs compared to equal mass binaries. The angular momentum of the disc is 8 − 10 Mo the disc mass over three orders of magnitude in the disk mass. While the associated nucleosynthesis shows no peculiarity, the simulated kilonovae are relatively dim compared with the GW170817 event. In particular, for distances compatible with GW190425, we find AB magnitudes always dimmer than 20 mag for the B, r and K bands, with brighter kilonovae associated to more asymmetric binaries and stiffer EOSs. We suggest that, even assuming a good coverage of GW190425’s sky location, the kilonova signal could hardly have been detected by present wide-field surveys and no firm constraints on the binary parameters or neutron star (NS) EOS can be argued from the lack of the detection.

The impact of nuclear reactions on observed signals from neutron star mergers is uncertain. We make a first attempt at quantifying that uncertainty by studying two cases that, intuitively, are the most extreme. In one case we assume that the reactions happen instantaneously (on timescales much faster than can be resolved in simulations). In the other we assume reactions do not happen (or occur on timescales much slower than the orbital, or simulation, timescale). We find that the resulting differences in the gravitational wave signal are small, but are potentially detectable by next generation observatories. These differences will impact on parameter recovery efforts driven by matching to templates extracted from numerical relativity simulations. Failure to take the effects discussed here into account may lead to systematic errors of unknown size being introduced into equation of state parameters inferred from gravitation wave signals.

Modelling the transport of neutrinos is crucial for core-collapse supernavoae and neutron star mergers. Evolving the complete Boltzmann equations for radiative transport is extremely expensive, so several more efficient approximation formalisms have been developed to face such difficulty. In this talk, I will focus on the truncated-moment formalism, considering only the first two moments (M1) within the gray approximation, which reduces Boltzmann seven-dimensional equation to a system of 3+1 equations. Just the first two moments (radiation energy density and radiation flux) are evolved, while all the other higher-order moments are provided by an approximate algebraic closure relation as a function of low-order moments. Solving this approximated formalism is still mathematically challenging because it is necessary to model the interaction of radiation and matter in several regimes where the evolution equations become stiff and behaves as an advection-diffusion problem. Here, I will present different high-order time integration schemes based on Implicit-Explicit Runge-Kutta (IMEX) methods designed to overcome the time-step restriction caused by such behaviour and analyze their performance in several numerical tests. Simulations were carried out using Julia, a high performance dynamic programming language that allows the resulting code to be easily executed on GPUs.

Authors: Nigel T. Bishop, Monos Naidoo, and Petrus J. van der Walt. Using linearized perturbations within the Bondi-Sachs formalism, we have considered the problem of a gravitational wave (GW) source surrounded by a spherical dust shell. It was shown that the shell causes the GWs to be modified both in magnitude and phase; and that if the shell is viscous, then the shear induced in the velocity field results in an energy transfer so damping the magnitude of the GWs. Both effects can be significant if the shell radius r is much smaller than the GW wavelength λ, and there are astrophysical scenarios for which the modification to the GW signal is large enough to be measurable. These scenarios include core collapse supernovae (CCSNe), quasinormal mode emission from the remnant of a binary neutron star (BNS) merger, and primordial GWs. In particular, there are feasible values of the CCSNe parameters for which viscous damping of the GW signal would be almost complete. In numerical relativity, an evolution of the full Einstein and matter field equations with GW extraction far from the source will properly include all effects described above. However, in situations such as CCSNe and BNS the complexity of the matter physics may necessitate an approximate treatment of the gravitational field or of the GW extraction. Results obtained in these cases may need to be corrected. The presentation builds on our previous work: 1. N.T. Bishop, M. Naidoo, and P.J. van der Walt (2020), Gen. Rel. Grav. 52:92, https://doi.org/10.1007/s10714-020-02740-9 arXiv:1912.08289 2. M. Naidoo, N.T. Bishop, and P.J. van der Walt (2021), Gen. Rel. Grav. 53:97, https://doi.org/10.1007/s10714-021-02841-z arXiv:2102.00060

In my talk, I will present the latest results of general relativistic hydrodynamic simulations of binary neutron star mergers with the recently introduced finite-temperature extension of the V-QCD equation of state. The V-QCD model is based on the gauge/gravity duality and provides a consistent description of nuclear and quark matter at densities beyond nuclear saturation which are realized in neutron stars. We simulate binaries that are consistent with GW170817 and study the imprint that a hadron-quark phase transition may have in the gravitational-wave emission and in the thermodynamical properties of the post-merger remnant.

The detection of a binary neutron star merger in 2017 through both gravitational waves and electromagnetic emission opened a new era of multimessenger astronomy. During the merger, several mechanisms like the Kelvin-Helmholtz instability, the winding up effect and the MRI, can amplify the initial magnetic field in the remnant to be powerful enoguh for launching a jet, with an associated short GRB. When performing simulations, simplified assumptions arise for the initial magnetic field strength and topology of the merging neutron stars. Here I will show convergent results by using high-resolution, large-eddy simulations of binary neutron star mergers, following the newly formed remnant for up to 30 milliseconds. I will specifically compare simulations with different initial magnetic field strenghts and configurations, going beyond the widespread-used aligned dipole confined within each star. I will show that the magnetic field is always amplified up to ~10^16 G in the bulk region of the remnant, while the initial topology is quickly forgotten in a timescale of few miliseconds after the merger due to the Kelvin-Helmholtz instability.

Binary neutron star mergers provide the unique opportunity to study matter at densities and temperatures unreachable in laboratories on earth. Its properties are encoded in the equation of state whose influence on the graviational wave signal of merging neutron stars can be used to constrain the physics of the strong interaction. But it is not only the equation of state which carries information about the underlying microphysics of neutron stars. Transport effects such as bulk viscosity are hypothesized to influence the dynamics binary neutron star mergers as well as opening another potential window to study and constrain matter in such extrem conditions. In this talk I will present the first steps towards the inclusion of transport effects in fully general-relativistic numerical simulations. Here, I will introduce a conservative 3+1 formulation of the equations of dissipative general-relativistic hydrodynamics which enables the inclusion of dissipative effects in numerical relativity codes. As a first application of this formulation I will show standard tests in special relativity such as the shock-tube problem. To further move on to general relativity, the problem of stationary, spherically symmetric black hole accretion including bulk viscosity is solved. Finally, we present first results of oscillating neutron stars with bulk viscosity in a three-dimensional setup with a dynamically evolved spacetime.

The construction of constraint-satisfying initial data is an essential element for the numerical exploration of the dynamics of compact-object binaries. While several codes have been developed over the years to compute generic quasi-equilibrium configurations of binaries comprising either two black holes, or two neutron stars, or a black hole and a neutron star, these codes are often not publicly available or they provide only a limited capability in terms of mass ratios and spins of the components in the binary. We here present a new open-source collection of spectral elliptic solvers that are capable of exploring the major parameter space of binary black holes (BBHs), binary neutron stars (BNSs), and mixed binaries of black holes and neutron stars (BHNSs). Particularly important is the ability of the spectral-solver library to handle neutron stars that are either irrotational or with an intrinsic spin angular momentum that is parallel to the orbital one. By supporting both analytic and tabulated equations of state at zero or finite temperature, the new infrastructure is particularly geared towards allowing for the construction of BHNS and BNS binaries. For the latter, we show that the new solvers are able to reach the most extreme corners in the physically plausible space of parameters, including extreme mass ratios and spin asymmetries, thus representing the most extreme BNS computed to date. Through a systematic series of examples, we demonstrate that the solvers are able to construct quasi-equilibrium and eccentricity-reduced initial data for BBHs, BNSs, and BHNSs, achieving spectral convergence in all cases. Furthermore, using such initial data, we have carried out evolutions of these systems from the inspiral to after the merger, obtaining evolutions with eccentricities ≲10−4−10−3, and accurate gravitational waveforms.

The solution of elliptic equations is an integral part of computational physics. This talk presents a new elliptic solver which is based on a discontinuous Galerkin scheme applicable to a wide class of elliptic partial differential equations. The solver is employing task-based parallelism and scales to a few thousand of compute cores. We show applications to the construction of initial data for binary black holes and to the study of thermal noise in the mirrors of interferometric gravitational wave detectors.

Binary black hole simulations become increasingly computationally expensive with smaller mass ratios because the Courant-Friedrich-Lewy condition imposes smaller time steps by the need to resolve the small black hole. Here we propose and explore a method for alleviating the scale disparity in simulations with mass ratios in the intermediate astrophysical range, where purely perturbative methods may not be adequate. A region much larger than the smaller black hole is excised from the numerical domain, and replaced with an analytical model approximating a tidally deformed black hole. We apply this idea to a toy model of a scalar charge in a circular geodesic orbit around a Schwarzschild black hole, solving for the massless Klein-Gordon field in a 3+1D framework using SpECTRE, the new discontinuous Galerkin code of the SXS collaboration.

Conservative finite difference methods have proven extremely robust and reliable for magnetohydrodynamics simulations of binary neutron star mergers. However, finite difference methods are generally less accurate and efficient than spectral methods when the solution is smooth, e.g. away from the stellar surfaces. The attractiveness of spectral methods has been demonstrated by thousands of long and highly accurate binary black hole simulations. Discontinuous Galerkin methods seek to provide the accuracy of spectral methods while also robustly capturing shocks in hydrodynamics simulations. I will give an overview of a discontinuous Galerkin-finite-difference hybrid method that inherits the best properties of both spectral and finite difference methods. I will show longterm simulations of magnetized and rotating neutron stars, as well as test problems to demonstrate the reliability and accuracy of the hybrid method.

Numerical relativity has been a frontier problem in computational science for over fifty years. The complexity of this research field spans multiple areas, from the properties of partial differential equations to high performance computing. As gravitational wave detectors continually improve, and new detectors come online, a new challenge will be to compute waveforms with higher fidelity and for a broader array of possible sources. We have developed Dendro-GR, a new computational platform for numerical relativity that uses a highly scalable octree with sparse, wavelet-based refinement. I will present some of initial results obtained with Dendro-GR for merging black hole binaries with different mass ratios. I will also talk about some possible new directions, including a new faster GPU-based extension of Dendro-GR, and some possible approaches to refine solutions in time as well as space. Performance and scaling results for these codes will be presented.

Numerical simulations are the only way to calculate exact gravitational waveforms from binary neutron star mergers and to design templates for gravitational-wave astronomy. Our knowledge about the physical properties of the inspiraling neutron stars, the mechanisms inherent in the inspiral and the massive object produced after the merger depends crucially on the accuracy of these numerical simulations. In this talk, I will present an entropy-based flux-limiting scheme that can be used in high-order, convergent simulations of neutron star spacetimes. First, I will demonstrate how the entropy produced during the inspiral can be used to detect the presence of shocks and will incorporate this idea into the design of a flux-limiter. Next, it will be shown that this entropy based flux-limiting scheme effectively tracks the stellar surface and the physical shocks during the inspiral. Finally, I will present the first neutron star merger simulations with such a method and will demonstrate up to fourth-order convergence in the gravitational waveform phase. The observed convergence rates are two orders higher than the one achieved by state-of-the-art high-order characteristic schemes.

GR-Athena++ is a general-relativistic, high-order, vertex-centered solver that extends the oct-tree, adaptive mesh refinement capabilities of the astrophysical (radiation) magnetohydrodynamics code Athena++. To simulate dynamical spacetimes GR-Athena++ uses the Z4c evolution scheme of numerical relativity coupled to the moving puncture gauge. Stable and accurate binary black hole merger evolutions are demonstrated in convergence testing, cross-code validation, and verification against state-of-the-art effective-one-body waveforms. GR-Athena++ leverages the task-based parallelism paradigm of Athena++ to achieve excellent scalability. Strong scaling efficiencies above 95% for up to 1.2×1e4 CPUs and excellent weak scaling up to 1e5 CPUs in a production binary black hole setup with adaptive mesh refinement are measured. GR-Athena++ thus allows for the robust simulation of compact binary coalescences and and offers a viable path towards numerical relativity at exascale.

We demonstrate the performance of the new code GR-Athena++ in evolving general relativistic hydrodynamics (GRHD) in a dynamically evolving spacetime. GR-Athena++ utilises the task-based parallelism and block based adaptive mesh refinement of the Athena++ code, as well as its approach to solving GRHD problems in stationary spacetimes; combined with new functionality to solve the Einstein equations in the Z4c formulation. We demonstrate the performance of this new code by simulating the evolution of Neutron Stars in GR-Athena++, removing the Cowling approximation assumed in previous work, presenting a fully dynamical spacetime evolution.

Numerical relativity plays a key role in the study of the dynamical and strong-field gravity regime of black hole and neutron star binary systems. The accuracy required nowadays for gravitational waves modeling and to investigate the outcome of compact binaries involving neutron stars and the matter and light emission can be obtained only through high-resolution simulations. Such simulations can be very costly and thus high-performance and scalable codes, capable of efficiently making use of the modern massively-parallel architectures available nowadays, are needed. It has become increasingly important to write codes that are portable, i.e. able to run on several architectures while keeping a high parallel efficicency. A possibility to achieve this goal is offered by the Kokkos C++ library, which allows one to easily write readable codes that can be compiled on virtually any architecture, included GPU nodes. We present the first results of single puncture evolutions with a new version of the oct-tree based code GR-Athena++, rewritten with the Kokkos paradigm. We compare CPUs and GPUs performance and we discuss the feasibility of binary black hole runs on GPUs, which are challenging due to the small memory available on GPUs. We achieve a speed-up of $\sim 50$ comparing the new Kokkos-based code run on a GPU with respect to GR-Athena++ and we argue that it is possible to run a medium/high resolution binary black hole runs using $\sim 4$ GPUs. These results demonstrate that a Kokkos version of GR-Athena++ can in principle be used to perform binary black hole evolution with a significant gain in speed, providing a valid option for exa-scale numerical relativity on the newest architectures.

We will quickly review the last 30 years of open general codes, platforms, and frameworks for Numerical Relativity, with special emphasis on the lessons that we have learnt and the challenges that lie ahead. We will also describe results and features of the last versions of different codes, with special emphasis on Simflowny: an open platform which automatically generates efficient parallel AMR code of scientific dynamical models for different simulation frameworks. It incorporates Finite-difference based methods for both smooth and non-smooth solutions, coupled with particle-based methods. It allows for generic systems of PDEs, and there are already implemented different formulations of the Einstein Equations, the full relativistic MHD equations and both the leakage and the M1 approximation to model the neutrino transport. Finally, we will discuss the changes in development languages, from Fortran and C++ legacy codes to Julia programming and GPUs.

The numerical solution of a system of hyperbolic PDEs all the way to future null infinity requires the knowledge of asymptotics. The Good-Bad-Ugly-F model is known to mimic the asymptotic properties of Einstein equations in generalized harmonic gauge. In this talk I will present the results of numerical evolution of this system, both in spherical symmetry and full 3D, with the scope of using these ideas to incorporate the hyperboloidal approach in numerical relativity simulations.

Since the first general-relativistic simulation of a binary neutron-star (BNS) merger less than a decade ago, the numerical relativity community has embraced increasingly sophisticated and accurate descriptions of the physics of these systems, most notably, the gravitational-wave signal. On the other hand, the electromagnetic signals, which are crucial to understand in the novel context of multi-messenger astrophysics, are significantly more challenging to predict. Indeed, these signals stem from low-density expanding ejecta that must be captured in its dynamics and composition for long timescales. Neutrinos play a prominent role in the cooling phase of the merger remnant and during the expansion of the neutron-rich ejecta, thus affecting both dynamics and composition of the outflow. Therefore, they are a key ingredient to achieve a realistic description of the merger phenomenon and its electromagnetic counterparts. In this talk, I will present a new implementation of a gray two-moment (M1) radiation transport scheme in full general relativity, aiming to capture the neutrino physics in BNS mergers. I will highlight the technical challenges involved in the development of what is arguably a cutting edge tool for its modelling capability in the context of BNS mergers, while also presenting a comprehensive set of tests that allow for cross-comparison and validation of such codes in the research landscape.

We have implemented fully adaptive hp mesh refinement into our pseudospectral NR code 'bamps', using heuristic indicators to drive dynamic refinement both in terms of cell size and cell resolution, while retaining both strong and weak scaling across several thousand processors. Because the mesh refinement infrastructure is independent of the evolved equations, it can be applied to a variety of physics problems. Its first major application has been the study of critical collapse of gravitational waves, where it allowed tuning towards criticality at a much higher efficiency compared to previous studies, letting us examine several new families of initial data that had not previously been studied.

We discuss the construction of quasi-equilibrium configurations of compact binaries, in which each component of the binary is modeled as a mixture of two ideal fluids. For the first fluid we use an ordinary baryonic-matter equation-of-state and for the second fluid, describing dark matter, we trial different dark-matter equation-of-state. We use the obtained quasi-equilibrium configurations as the basis for numerical relativity merger simulations and analyse the impact of the dark matter on the system's dynamical properties and the waveforms of the emitted gravitational radiation.

Fuzzy dark matter (FDM) is an exciting alternative to the standard cold dark matter (CDM) paradigm, reproducing its many successful large scale predictions, but solving most of the existing tension with small scale (galactic) observations. FDM models postulate that dark matter (DM) is constituted by light bosons of mass~$m_\psi\sim 10^{-22}\,\text{eV}$, which behave like CDM on scales larger than the de Broglie wavelength~$\lambda_\text{dB}$ ($\sim 1\, \text{kpc}$, in a typical galaxy), but are subjected to a \emph{wave pressure} on smaller scales. One of the most solid predictions from FDM cosmological numerical evolutions is the condensation of a soliton core at the center of haloes. These solitons -- also known as boson stars (BSs) or oscillatons -- are stationary (or long-lived), self-gravitating scalar field configurations supported by wave pressure. However, if these solitons host a \emph{parasite} supermassive black hole (SMBH) they will be swallowed by their guest, as implied by the no-hair theorem. It is thus crucial to understand in detail the accretion process and the timescales involved. Several works have assumed (based on crude estimates) that the soliton remains largely unaltered by the presence of a central SMBH and used then, e.g., inner galactic rotation curves to put constraints on the mass~$m_\psi$ of the FDM particles. In this work we use numerical relativity to self-consistently solve the problem of accretion of a BS by a central SMBH. We find the different stages of accretion, estimate the BS lifetime.

A promising way to study the properties and behaviour of dark matter on small scales is through the effects it may have on binary mergers. In particular, its presence around black holes may lead to a distinctive dephasing of the signal due to dynamical friction. In this talk I will present how we calculate this force numerically and discuss how our findings extend some already existing analytic results to the relativistic regime, where this effect can become detectable.

The multi-messenger detection of GW170817, GRB170817A, and AT2017gfo originating from the merger of two neutron stars have been a scientific breakthrough. Under the assumption that dark matter accumulates in and around neutron stars, multi-messenger observations of compact binary mergers will provide a new way to search for and constrain the nature of dark matter. In this context, we extended the numerical-relativity code BAM to enable the simulation of neutron stars that contain ordinary baryonic matter, but also a mirror dark matter component. We performed single star tests to verify our code changes and perform the binary neutron star simulations using a two-fluid approach. Despite the drawback that our simulations are based on superimposed initial configurations, we find that the presence of dark matter reduces noticeably the lifetime of the remnant and changes the amount of ejected material. Consequently, electromagnetic signals that arise due to the merger of binary neutron star mergers admixed with dark matter will be different than those from `pure' neutron stars. Given the increasing sensitivity of gravitational-wave interferometers and other observational facilities, our analysis gives a new perspective on how to probe the presence of dark matter.

Scalar fields around compact objects are of interest for scalar-tensor theories of gravity and dark matter models consisting of a massive scalar, e.g. axions. These fields can form long-lived clouds around black holes via superradiance or simple gravitational accretion. Black hole binary mergers occuring inside such clouds, or in dense scalar dark matter environments, may modify the gravitational wave signal compared to the vacuum case. I will discuss our numerical relativity simulations of scalar field accretion, and binary mergers in dense scalar field environments, and the challenges we face in producing realistic templates of such “dirty” binary black hole mergers for the next generation of gravitational wave observatories. In particular I will explore the impact of the scalar field initial data.

Supplementing General Relativity (GR) with an additional scalar degree of freedom yields to a simple and popular class of gravitational theories known as scalar-tensor theories. While scalar-tensor theories can produce self-accelerated cosmic expansion without a cosmological constant, they typically produce also local deviations from GR on small scales. However, some theories possess “screening mechanisms” that locally produce a GR-like phenomenology, potentially passing existing constraints. Here we will study numerically the evolution of neutron star systems on two specific scalar-tensor theories (i.e., one with this screening mechanism and the other without). Both cases show a interesting non-linear dynamics which could help to constrain these alternative theories by using current and near future gravitational wave observations.

In this talk I will discuss recent progress on modelling black hole binaries in higher derivative theories of gravity. I shall consider first cubic Horndeski theories for which the higher derivatives are in the scalar sector only. In these theories in the weakly coupled regime, even though the differences with standard general relativity can be locally small, they accumulate over the lifetime of a binary resulting in waveforms that significantly deviate from those in general relativity. In the second part of the talk, I will discuss binaries in a certain eight-derivative theory of gravity. I will first discuss how such a theory can be rendered well-posed and then will present results on binaries as well as pointing out generic issues that may arise when modelling black hole binaries in such theories.

We present numerical simulations of boson-star binary systems and extract the gravitational waves emitted by these systems. We discuss numerical methods and with particular focus on improved initial data and the spurious effects that may result from inappropriate choices for the construction of initial data.

Recently, a new class of fully nonlinear numerical solutions to the static, spherically symmetric Einstein-Klein-Gordon system for a collection of an arbitrary odd number of complex scalar fields with an internal U(N) symmetry was found. These solutions are parametrized by an angular momentum number l, an excitation number n, and a continuous parameter representing the amplitude of the fields. They are regular at every point, possess a finite total mass, and they generalize the standard spherically symmetric boson stars to arbitrary values of l. In this talk, we discuss the main properties of these l-boson stars, including their stability with respect to small perturbations, their compactness, the large l limit, their pressure anisotropy and some preliminary results regarding their interpretation as a solution of the semi-classical Einstein-Klein-Gordon equations.

Initial data play a vital role in numerical simulations -- the ability to construct constraint satisfying initial data closely representing the physical system determines the quality of the subsequent numerical evolution. It has been shown by Helfer et. al. (arxiv:1802.06733) that plainly superposed initial data can lead to significant constraint violations that affect the physics of the evolution of boson stars in the equal-mass case. In particular, plainly superposed initial data may lead to pre-mature collapse of the stars to a black hole and altered gravitational wave signals. However, the method proposed in Helfer et.al. is limited to the equal-mass case. In this talk we develop a generalisation of their method to unequal-mass binaries that recovers the Helfer construction in the equal-mass limit. Whilst a constraint solver is needed to obtain fully constraint satisfying initial data, we can still significantly reduce constraint violations by using our proposed construction. Furthermore, the significant reductions in the constraint violations around the stars avoid unphysical features in the boson star dynamics and the gravitational wave signals. We will present the methodology of the construction and show its applications to unequal-mass boson star collisions.

We are in the era of gravitational wave (GW) astronomy, with the rapidly growing catalogue of gravitational wave events detected by the LIGO-Virgo collaboration. These detections allow us to study the strong regime of gravity through GW signals produced by the coalescence of compact objects. In this talk, I will present our numerical simulations on the coalescence of binary Exotic Compact Objects, considering two different types of binaries: boson stars and fermion-boson stars (a model for neutron stars that contain a small fraction of dark matter). These binaries lead to different dynamics and gravitational signatures during their coalescence respect to standard compact objects, which might be crucial to distinguish them from other binaries with current/future detections by ground-based interferometers.

New light fundamental fields are natural candidates for all or a fraction of dark matter. Self-gravitating structures of such fields might be common objects in the universe, and could comprise even galactic haloes. These structures would interact gravitationally with black holes, process of the utmost importance, since it dictates their lifetime, the black hole motion and possible gravitational radiation emission. I will present a study on the dynamics of a black hole piercing through a much larger fully relativistic boson star, made of a complex minimally coupled massive scalar without self-interactions. As the black hole pierces through the bosonic structure, it is slowed down by accretion and dynamical friction, giving rise to gravitational wave emission. We find evidence of a "gravitational atom" left behind as a product of the process.

Particle physics models of dark matter, and extensions to the Standard Model, predict the existence of a large abundance of light scalar degrees of freedom in the universe. From a diffuse cloud, these can form into clumps of energy - boson stars. Additionally, due to their high compactness, close to that of black holes, these solutions serve as test beds to study the non-linear dynamics of a large class of ultra-compact objects. We outline a method to obtain constraint satisfying binary Boson star initial data, and discuss the dynamics of merging binaries, focusing specifically on the formation of rotating Boson star remnants.

Einstein-dilaton-Gauss-Bonnet is a theory of modified gravity in which a scalar field, called dilaton, is nonminimally coupled to the metric via an exponential function. Black holes (BHs) in this theory are particularly interesting since they possess a critical configuration with minimum mass and finite Hawking temperature. This means that a critical BH loses mass due to Hawking's radiation, but it is not clear what is its fate after this process, since it cannot reach a final configuration with lower mass. In a recent work we studied this problem by means of fully nonlinear numerical evolutions of spherically symetric BH spacetimes (Corelli, De Amicis, Ikeda, and Pani, to appear on arXiv in the next few days). Specifically, by simulating the collapse of wave packets of a phantom scalar field we have been able to dynamically reduce the BH mass, reproducing the effect of the Hawking's evaporation process. In this talk I will present our results focusing on the case in which the BH mass falls below the critical value. In particular I will show that in the formalism we used an elliptic region with high curvature appears outside the apparent horizon, and I will discuss how this could possibly hint to a violation of the weak cosmic censorship conjecture.

It was recently shown that a broad class of gravity theories that couples a dynamical scalar field to the Gauss-Bonnet invariant can lead to spontaneous scalarization of black holes, allowing these objects to grow "scalar hair" once certain conditions are met and to remain "bald" otherwise. While most works on the topic have focused on isolated black holes, progress has recently been made in understanding this effect in binary black hole systems. I will present an overview of what has been achieved so far in this context. I will discuss new phenomena that takes place in black-hole binaries and explore some of the potential observational consequences of these results.

In recent years, gravitational wave observation of black holes have furnished new opportunities to test our understanding of gravity in the strong field, highly dynamical regime. However, in order to perform model dependent tests of General Relativity with these observations one needs accurate waveforms for modified gravity theories which is still an outstanding theoretical problem. In this talk, we will discuss some recent progress in this regard, in particular the use of the modified harmonic formulation of Horndeski theories of gravity to numerically solve for the dynamics of binary black holes systems in such theories, including the merger phase. We will focus on the particular case of Einstein-Scalar-Gauss-Bonnet gravity as it is not only motivated by effective field theory but also admits scalar hairy black hole solutions.

The black hole merger in scalar-Gauss-Bonnet gravity can lead to dynamical descalarization this is a spontaneous release of the scalar hair of the newly formed black hole. Depending on the exact form of the Gauss-Bonnet coupling function, the stable scalarized solutions can be either continuously connected to the Schwarzschild black hole, or the transitions between the two can happen with a jump. In the present talk we will discuss the black hole head-on collisions in scalar-Gauss-Bonnet gravity prone to dynamical descalization and the gravitational wave manifestation of such a jump. The distinct signature in the gravitational wave signal will share similarities with the effects expected from first order matter phase transitions happening during neutron star binary mergers.

Though General Relativity has been successfully tested so far, concepts such as dark energy and string theory suggest the need of modifying it. Scalar-tensor theory is one of the most popular alternatives discussed. We produce studies of stellar core collapse in spherical symmetry that were performed by adapting the numerical code GR1D to the case of massive scalar-tensor gravity. We systematically explore the parameter space that characterizes the progenitor stars, the equation of state and the scalar-tensor theory of the core collapse events. We identify a remarkably simple and straightforward classification scheme of the resulting collapse events. For any given set of parameters, the collapse leads to one of three end states, a weakly scalarized neutron star, a strongly scalarized neutron star or a black hole, possibly formed in multiple stages. The latter two end states can lead to strong gravitational-wave signals that may be detectable in present continuous-wave searches with ground-based detectors. We identify a very sharp boundary in the parameter space that separates events with strong gravitational-wave emission from those with negligible radiation.

I present the first numerical-relativity simulations for black holes in Quadratic Gravity, i.e., including the leading-order (quadratic) curvature corrections to General Relativity. I review the nonlinear degrees of freedom and discuss a well-posed initial value formulation for Quadratic Gravity, both in spherical symmetry and in (3+1) dimensions. In spherical symmetry, self-convergence tests and decaying constraint violations confirm numerically stable evolution of perturbations on Minkowski as well as Schwarzschild spacetime. For the first time, I will also present full (3+1) simulations of black-hole initial data, which enable ongoing work to simulate black-hole binary systems in Quadratic Gravity.

Scalar-tensor theories with first-derivative self interactions, known as k-essence, may provide interesting phenomenology on cosmological scales. On smaller scales, however, initial value evolutions (which are crucial for predicting the behavior of astrophysical systems, such as binaries of compact objects) may run into instabilities related to the Cauchy problem becoming potentially ill-posed. Moreover, on local scales the dynamics may enter in the nonlinear regime, which may lie beyond the range of validity of the infrared theory. Completions of k-essence in the ultraviolet, when they are known to exist, mitigate these problems, as they both render Cauchy evolutions well-posed at all times, and allow for checking the relation between nonlinearities and the low energy theory's range of validity. In this talk, I will present an exploration of these issues explicitly by considering an ultraviolet completion to k-essence and performing vacuum 1+1 dynamical evolutions within it. I will argue that these results are consistent with the low-energy theory, and will compare with evolutions obtained with the low-energy theory suitably deformed with a phenomenological "fixing the equations" approach.