The dynamics and observable signatures of neutron star mergers are governed by physics under the most extreme conditions. They are particularly impacted by the high-density equation of state, which for the most sophisticated models is usually available in the form of tables. Numerical relativity codes usually evolve particularly well-behaved numerical ("conservative") variables, but at the price that the physically interesting ("primitive") variables need to be found at every computational element and at every integration sub-step by means of expensive (and not always successful) root-finding algorithms. We have recently developed the Lagrangian numerical relativity code SPHINCS_BSSN which evolves the spacetime on an adaptive mesh with well tested methods, but the fluid is evolved by means of freely moving particles. Since our evolution equations differ from those of conventional numerical relativity, we need to develop new conservative-to-primitive algorithms if we want to use tabulated equations of state. We present here three such algorithms: a 3D and a 2D Newton-Raphson method and a 1D root-finding algorithm based on Ridders' method. We find the 3D method to be very fast and robust with an average failure fraction in a full-blown neutron star merger simulation (with the DD2 equation of state) well below 1

Neutron star mergers hold the key to several grand challenges of contemporary (astro-)physics. In view of the upcoming next generation of ground-based detectors, it is crucial to keep improving theoretical predictions to harvest the full scientific returns from these investments. We introduce here a substantial update of our Lagrangian numerical relativity code SPHINCS_BSSN. Apart from changing our unit system, we add constraint damping terms to the BSSN spacetime evolution equations. We demonstrate that this measure reduces, without noteworthy computational cost, the Hamiltonian constraint violations by more than an order of magnitude. We further implement contributions to thermal energy and pressure that are based on Fermi liquid theory and contain a parametrization of the Dirac effective mass. These terms can be combined with any cold equation of state, and they enhance the physical realism of our simulations and introduce a physics-based concept of a temperature. In a set of merger simulations, we demonstrate good agreement with other temperature-dependent numerical relativity simulations. We find that different parametrizations of the Dirac effective mass can translate into shifts of ∼ 150 Hz in the dominant post-merger gravitational wave peak frequency.

The ejection of neutron-rich matter is one of the most important consequences of a neutron star merger. While the bulk of the matter is ejected at fast, but non-relativistic velocities (∼0.2c), a small amount of mildly relativistic dynamic ejecta have been seen in a number of numerical simulations. Such ejecta can have far reaching observational consequences ranging from the shock breakout burst of gamma-rays promptly after the merger, to an early (∼ 1 hour post-merger) blue kilonova precursor signal, to synchrotron emission years after the merger ("kilonova afterglow"). These all potentially carry the imprint of the binary system parameters and the equation of state. By analyzing Lagrangian simulations in full General Relativity, performed with the code SPHINCS_BSSN, we identify two ejection mechanisms for fast ejecta: i) about 30% of the ejecta with v> 0.4c are "sprayed out" from the shear interface between the merging stars and escape along the orbital plane and ii) the remaining ∼ 70% of the fast ejecta result from the central object "bouncing back" after strong, general-relativistic compression. This "bounce component" is ejected in a rather isotropic way and reaches larger velocities (by ∼0.1c) so that its faster parts can catch up with and shock slower parts of the spray ejecta. Even for a case that promptly collapses to a black hole, we find fast ejecta with similar properties to the non-collapsing case, while slower matter parts are swallowed by the forming black hole. We discuss observational implications of these fast ejecta, including shock breakout and kilonova afterglow.

We present version 1.0 of our Lagrangian Numerical Relativity code SPHINCS_BSSN. This code evolves the full set of Einstein equations, but contrary to other Numerical Relativity codes, it evolves the matter fluid via Lagrangian particles in the framework of a high-accuracy version of Smooth Particle Hydrodynamics (SPH). The major new elements introduced here are: i) a new method to map the stress–energy tensor (known at the particles) to the spacetime mesh, based on a local regression estimate; ii) additional measures that ensure the robust evolution of a neutron star through its collapse to a black hole; and iii) further refinements in how we place the SPH particles for our initial data. The latter are implemented in our code SPHINCS_ID which now, in addition to LORENE, can also couple to initial data produced by the initial data library FUKA. We discuss several simulations of neutron star mergers performed with SPHINCS_BSSN_v1.0, including irrotational cases with and without prompt collapse and a system where only one of the stars has a large spin (χ = 0.5).

We present the first neutron star merger simulations performed with the newly developed Numerical Relativity code SPHINCS_BSSN. This code evolves the spacetime on a mesh using the BSSN formulation, but matter is evolved via Lagrangian particles according to a high-accuracy version of general-relativistic Smooth Particle Hydrodynamics (SPH). Our code contains a number of new methodological elements compared to other Numerical Relativity codes. The main focus here is on the new elements that were introduced to model neutron star mergers. These include (a) a refinement (fixed in time) of the spacetime-mesh, (b) corresponding changes in the particle–mesh mapping algorithm and (c) a novel way to construct SPH initial data for binary systems via the recently developed “Artificial Pressure Method.” This latter method makes use of the spectral initial data produced by the library LORENE, and is implemented in a new code called SPHINCS_ID. While our main focus is on introducing these new methodological elements and documenting the current status of SPHINCS_BSSN, we also show as a first application a set of neutron star merger simulations employing “soft” ( $$\varGamma =2.00$$ ) and “stiff” ( $$\varGamma =2.75$$ ) polytropic equations of state.

The observation of gravitational waves from compact objects has now become an active part of observational astronomy. For a sound interpretation, one needs to compare such observations against detailed Numerical Relativity simulations, which are essential tools to explore the dynamics and physics of compact binary mergers. To date, essentially all simulation codes that solve the full set of Einstein’s equations are performed in the framework of Eulerian hydrodynamics. The exception is our recently developed Numerical Relativity code SPHINCS_BSSN which solves the commonly used BSSN formulation of the Einstein equations on a structured mesh and the matter equations via Lagrangian particles. We show here, for the first time, SPHINCS_BSSN neutron star merger simulations with piecewise polytropic approximations to four nuclear matter equations of state. In this set of neutron star merger simulations, we focus on perfectly symmetric binary systems that are irrotational and have 1.3 M⊙ masses. We introduce some further methodological refinements (a new way of steering dissipation, an improved particle–mesh mapping), and we explore the impact of the exponent that enters in the calculation of the thermal pressure contribution. We find that it leaves a noticeable imprint on the gravitational wave amplitude (calculated via both quadrupole approximation and the Ψ4 formalism) and has a noticeable impact on the amount of dynamic ejecta. Consistent with earlier findings, we only find a few times 10−3M⊙ as dynamic ejecta in the studied equal mass binary systems, with softer equations of state (which are more prone to shock formation) ejecting larger amounts of matter. In all of the cases, we see a credible high-velocity (∼0.5…0.7c) ejecta component of ∼10−4M⊙ that is launched at contact from the interface between the two neutron stars. Such a high-velocity component has been suggested to produce an early, blue precursor to the main kilonova emission, and it could also potentially cause a kilonova afterglow.

We provide expansions of the Detweiler–Whiting singular field for a particle with a scalar field moving along arbitrary, planar accelerated trajectories in Schwarzschild spacetime. We transcribe these results into mode-sum regularization parameters, computing previously unknown terms that increase the convergence rate of the mode-sum. We test our results by computing the self-force along a variety of accelerated trajectories. For non-uniformly accelerated circular orbits we present results from a new 1+1D discontinuous Galerkin time-domain code which employs an effective source. We also present results for uniformly accelerated circular orbits and accelerated bound eccentric orbits computed within a frequency-domain treatment. Our regularization results will be useful for computing self-consistent self-force inspirals where the particle's worldline is accelerated with respect to the background spacetime.

We introduce a new relativistic astrophysics code, SpECTRE, that combines a discontinuous Galerkin method with a task-based parallelism model. SpECTRE's goal is to achieve more accurate solutions for challenging relativistic astrophysics problems such as core-collapse supernovae and binary neutron star mergers. The robustness of the discontinuous Galerkin method allows for the use of high-resolution shock capturing methods in regions where (relativistic) shocks are found, while exploiting high-order accuracy in smooth regions. A task-based parallelism model allows efficient use of the largest supercomputers for problems with a heterogeneous workload over disparate spatial and temporal scales. We argue that the locality and algorithmic structure of discontinuous Galerkin methods will exhibit good scalability within a task-based parallelism framework. We demonstrate the code on a wide variety of challenging benchmark problems in (non)-relativistic (magneto)-hydrodynamics. We demonstrate the code's scalability including its strong scaling on the NCSA Blue Waters supercomputer up to the machine's full capacity of 22,380 nodes using 671,400 threads.

Due to the numerical complexities of studying evolution in an anisotropic quantum spacetime, in comparison to the isotropic models, the physics of loop quantized anisotropic models has remained largely unexplored. In particular, robustness of bounce and the validity of effective dynamics have so far not been established. Our analysis fills these gaps for the case of vacuum Bianchi-I spacetime. To efficiently solve the quantum Hamiltonian constraint we perform an implementation of the Cactus framework which is conventionally used for applications in numerical relativity. Using high performance computing, numerical simulations for a large number of initial states with a wide variety of fluctuations are performed. Big bang singularity is found to be replaced by anisotropic bounces for all the cases. We find that for initial states which are sharply peaked at the late times in the classical regime and bounce at a mean volume much greater than the Planck volume, effective dynamics is an excellent approximation to the underlying quantum dynamics. Departures of the effective dynamics from the quantum evolution appear for the states probing deep Planck volumes. A detailed analysis of the behavior of this departure reveals a non-monotonic and subtle dependence on fluctuations of the initial states. We find that effective dynamics in almost all of the cases underestimates the volume and hence overestimates the curvature at the bounce, a result in synergy with earlier findings in the isotropic case. The expansion and shear scalars are found to be bounded throughout the evolution.

We report on a set of long-term general-relativistic three-dimensional (3D) multi-group (energy-dependent) neutrino radiation-hydrodynamics simulations of core-collapse supernovae. We employ a full 3D two-moment scheme with the local M1 closure, three neutrino species, and 12 energy groups per species. With this, we follow the post-core-bounce evolution of the core of a nonrotating 27-M-circle dot progenitor in full unconstrained 3D and in octant symmetry for >= 380 ms. We find the development of an asymmetric runaway explosion in our unconstrained simulation. We test the resolution dependence of our results and, in agreement with previous work, find that low resolution artificially aids explosion and leads to an earlier runaway expansion of the shock. At low resolution, the octant and full 3D dynamics are qualitatively very similar, but at high resolution, only the full 3D simulation exhibits the onset of explosion.

Understanding the physics of a Gamma-Ray Burst (GRB) demands truly petascale multiphysics simulations on the largest and most modern High Performance Computing (HPC) systems. To take full advantage of such systems extreme requirements are placed on the computational infrastructure on which such codes are based. In this report we describe the work within the Cactus framework to enable such demanding simulations on Blue Waters.

In recent years, numerical simulations with Gaussian initial states have demonstrated the existence of a quantum bounce in loop quantum cosmology in various models. A key issue pertaining to the robustness of the bounce and the associated physics is to understand the quantum evolution for more general initial states, which may depart significantly from Gaussianity and may have no well defined peakedness properties. The analysis of such states, including squeezed and highly non-Gaussian states, has been computationally challenging until now. In this paper, we overcome these challenges by using the Chimera scheme for the spatially flat, homogeneous and isotropic model sourced with a massless scalar field. We demonstrate that the quantum bounce in this model occurs even for states that are highly squeezed or are non-Gaussian with multiple peaks and with little resemblance to semi-classical states. The existence of the bounce is found to be robust, being independent of the properties of the states. The evolution of squeezed and non-Gaussian states turns out to be qualitatively similar to that of Gaussian states, and satisfies strong constraints on the growth of the relative fluctuations across the bounce. We also compare the results from the effective dynamics and find that, although it captures the qualitative aspects of the evolution for squeezed and highly non-Gaussian states, it always underestimates the bounce volume. We show that various properties of the evolution, such as the energy density at the bounce, are in excellent agreement with the predictions from an exactly solvable loop quantum cosmological model for arbitrary states.

The existence of a quantum bounce in isotropic spacetimes is a key result in loop quantum cosmology (LQC), which has been demonstrated to arise in all the models studied so far. In most of the models, the bounce has been studied using numerical simulations involving states which are sharply peaked and which bounce at volumes much larger than the Planck volume. An important issue is to confirm the existence of the bounce for states which have a wide spread, or which bounce closer to the Planck volume. Numerical simulations with such states demand large computational domains, making them very expensive and practically infeasible with the techniques which have been implemented so far. To overcome these difficulties, we present an efficient hybrid numerical scheme using the property that at the small spacetime curvature, the quantum Hamiltonian constraint in LQC, which is a difference equation with uniform discretization in volume, can be approximated by a Wheeler-DeWitt differential equation. By carefully choosing a hybrid spatial grid allowing the use of partial differential equations at large volumes, and with a simple change of geometrical coordinate, we obtain a surprising reduction in the computational cost. This scheme enables us to explore regimes which were so far unachievable for the isotropic model in LQC. Our approach also promises to significantly reduce the computational cost for numerical simulations in anisotropic LQC using high performance computing.

A key result of isotropic loop quantum cosmology is the existence of a quantum bounce which occurs when the energy density of the matter field approaches a universal maximum close to the Planck density. Though the bounce has been exhibited in various matter models, due to severe computational challenges, some important questions have so far remained unaddressed. These include the demonstration of the bounce for widely spread states, its detailed properties for the states when matter field probes regions close to the Planck volume and the reliability of the continuum effective spacetime description in general. In this manuscript we rigorously answer these questions using the Chimera numerical scheme for the isotropic spatially flat model sourced with a massless scalar field. We show that, as expected from an exactly solvable model, the quantum bounce is a generic feature of states even with a very wide spread, and for those which bounce much closer to the Planck volume. We perform a detailed analysis of the departures from the effective description and find some expected, and some surprising results. At a coarse level of description, the effective dynamics can be regarded as a good approximation to the underlying quantum dynamics unless the states correspond to small scalar field momenta, in which case they bounce closer to the Planck volume or are very widely spread. Quantifying the amount of discrepancy between the quantum and the effective dynamics, we find that the departure between them depends in a subtle and non-monotonic way on the field momentum and different fluctuations. Interestingly, the departures are generically found to be such that the effective dynamics overestimates the spacetime curvature, and underestimates the volume at the bounce.

We revisit the problem of computing the self-force on a scalar charge moving along an eccentric geodesic orbit around a Schwarzschild black hole. This work extends previous scalar self-force calculations for circular orbits, which were based on a regular "effective'' point-particle source and a full 3D evolution code. We find good agreement between our results and previous calculations based on a (1 + 1) time-domain code. Finally, our data visualization is unconventional: we plot the self-force through full radial cycles to create "self-force loops,'' which reveal many interesting features that are less apparent in standard presentations of eccentric-orbit self-force data.