David Radice - Homepage

Neutron star merger simulations are extremely computationally challenging. WhiskyTHC requires 6 hours of calculation on 512 cores to simulate one millisecond of evolution of a binary. This is for our standard resolution of 180 meters. Doubling the resolution increases the computational cost by a factor 16. Unfortunately, we know that, during the merger, there are small-scale magnetohydrodynamical instabilities that generate turbulence with typical driving scales of few meters or less. These are scales that are unaccessible to our simulations, even when running on the World's largest supercomputers. How can we model their effect in simulations?

Our approach is to use turbulence models initially developed for weather forecast and other "terrestrial" applications. Indeed, WhiskyTHC is the first, and so far only, code to include an effective, sub-grid scale, treatment of turbulence. This is based on the general-relativistic (GR) extension of the large-eddy simulations (LES) method, which we recently introduced.

For our first GRLES simulations, we adopted the so-called mixing-length closure of turbulence and we studied the impact of (sub-grid scale) turbulent mixing in the evolution of the merger remnant. We found that turbulence could have a potentially very important role in the evolution of the binary after contact and might alter the gravitational-wave and neutrino signals from the binary. On the other hand, for the most conservative values of the mixing length parameter ℓ_mix, we found that these effects are minor. This is very good news for all existing models of the gravitational-wave emissions from merging neutron star, which did not include any turbulence model. However, the final words will have to wait until we will be able to actually measure ℓ_mix by means of more restricted, but much higher resolution, GRMHD simulations.

Gravitational wave strain and spectra generated during the collision of two neutron stars modeled with two different nuclear equations of state. The DD2 equation of state includes only "normal" nuclear matter, while the BHBΛφ also includes Λ-hyperons.

Can we use coming gravitational-waves observations of merging neutron stars to detect the appearance of exotic particles or phase transitions in nuclear matter at several times nuclear density?

The answer to this question is "yes" according to our most recent simulations. The typical outcome of neutron star mergers is the formation of an object, called hypermassive neutron star, which survives only for short time before collapsing to form a black hole. The densities in this remnant can be several times larger than the density in the neutron stars before they collide. Now, we expect that phase transitions or other changes in the nature of matter in the remnant will be the result of the system re-arranging to find a new energy minimum. Where does the energy difference go? In gravitational waves! As a consequence, gravitational waves can be used to confirm or rule out the appearance of exotic states of matter in these extreme conditions.

We recently performed a series of high-resolution simulations of the late inspiral and merger of neutron star binaries using the high-order methods implemented in WhiskyTHC.

We studied the development and saturation of an hydrodynamical instability that has been recently discovered in eccentric neutron star mergers that included neutron star spin by Paschalidis and collaborators. We found that this instability is not restricted to eccentric or spinning neutron star mergers, but that the development of the one-armed spiral instability is a generic outcome of the merger. Furthermore, if detected in gravitational-waves, this instability would help constrain unknown parameters of the equation of state of matter at super-nuclear densities.

Complete gravitational wave spectrum for an M1 = M2 = 1.35 Msun binary neutron star merger at 10 Mpc with edge-on orientation. We use the MS1b interaction model to describe the equation of state of neutron stars (the "stiff" case above).

For this study, we constructed hybrid waveforms by combining our numerical-relativity waveforms with state-of-the-art effective-one-body waveforms we generated using a publicly available code. Using these waveforms, we studied the detectability of this instability and of other components of the gravitational-wave signal by ground-based gravitational wave observatories. Unfortunately, we find that the detection of this instability by Advanced LIGO is unlikely, but might be possible for third generation detectors such as the Einstein Telescope.

We released the complete hybrid gravitational waveforms on Zenodo.

Rest mass density in the orbital plane for a binary coomposed of a 1.4 and a 1.2 solar mass neutron stars. The top panels show a simulation that does not account for turbulent viscosity, while the bottom panels show a simulation that included a treatment for turbulent viscosity. The black countours enclose matter that is gravitationally unbound. Viscous heating of the tidal stream between the primary and the secondary stars can significantly enhance the mass loss from the binary.

Light curve modeling of the kilonova that followed GW170817 revealed the presence of a few percent of solar mass of material ejected with large velocity of about 0.3 c. Extant neutron star merger simulations find that material is indeed ejected with large velocities during the collision between the stars. However, the amount of the fast-moving ejecta is significantly smaller than that required by the kilonova models. It appears that something is missing in the current simulations.

We have recently performed the first simulations of unequal mass binary neutron star mergers that included a sub-grid treatment of turbulent viscosity. These new simulations show a new mechanism for mass ejection during mergers. We find that mass exchange streams between the neutron star prior to merger become super-heated due to the effects of turbulent viscosity and reach temperatures as high as 20 MeV (230 billion Kelvins!). The thermal pressure then drives powerful outflows with high velocities. Viscous-dynamical ejecta might explain the fast outflows seen for GW170817.

Our simulations also make two predictions that could be tested. First, they predict that some of the ejecta expand so rapidly that some of the neutrons might escape capture by seed nuclei during the decompression. The neutrons would then decay on a timescale of minutes and produce a bump in the UV bands on a time scale of an hour. in the first hour of the merger. Second, the radio fluency from GW170817 might reveal the presence of a substantial amount of fast-moving ejecta in the next few months or years.

We present a comprehensive study of mass ejection, nucleosynthesis, and associated electromagnetic counterparts from binary neutron star mergers. In the course of 2 years, we have performed 59 neutron star merger simulations with a microphysical treatment of neutron star merger. We find that typically a few 10^-3 solar masses of material are ejected dynamically during the mergers. The amount and the properties of these outflow depend on binary parameters and on the NS equation of state (EOS). A small fraction of these ejecta, typically 10^-6 solar masses, is accelerated by shocks formed shortly after merger to velocities larger than 0.6 c and produces bright radio flares on timescales of weeks, months, or years after merger. Their observation could constrain the strength with which the NSs bounce after merger and, consequently, the EOS of matter at extreme densities. The dynamical ejecta robustly produce second and third r-process peak nuclei with relative isotopic abundances close to solar. The production of light r-process elements is instead sensitive to the binary mass ratio and the neutrino radiation treatment. Accretion disks of up to 0.2 solar masses are formed after merger, depending on the lifetime of the remnant. In most cases, neutrino- and viscously-driven winds from these disks dominate the overall outflow. Finally, we generate synthetic kilonova light curves and find that kilonovae depend on the merger outcome and could be used to constrain the NS EOS.

Estimated outcomes for the viscous evolution of a binary neutron star system producing a stable remnant after merger. The grey shaded area shows the set of all rigidly-rotating equilibrium configurations. The solid line is a conservative estimate of the mass ejection and a possible trajectory for the viscous evolution. The blue shaded region denotes the range of all possible outcomes of the viscous evolution, for which we identify two possible regimes. The first corresponds to the "normal" ejection of matter due to nuclear recombination of the disk. The second to the case in which more massive outflows are produced as the remnant settles to a rigidly-rotating equilibrium. We find that the merger remnant has enough angular momentum to unbind more than 0.1 solar masses of material.

The outcome of neutron star mergers depends on the total mass of the system and on the poorly known equation of state of dense nuclear matter. Binaries with mass significantly larger than the maximum mass for a nonrotating neutron star result in prompt black-hole formation. Binaries with lower masses, but above the maximum mass of isolated rigidly rotating neutron stars, result in the formation of so-called hypermassive neutron stars, objects temporarily supported against gravitational collapse by the large differential rotation. Lower mass systems could produce remnants that are stable even after the effective viscosity due to turbulence has removed the differential rotation.

The identification of the outcome of the merger of binary neutron star systems with different masses would yield a precise measurement of the maximum mass of neutron stars. This, in turn, would constraint the equation of state of matter at extreme densities. It is therefore important to identify signatures indicative of the formation of long-lived remnants.

In a recent work, we studied the formation of long-lived merger remnants in our simulations. We showed that these are born with so much angular momentum that they will inevitabily become unstable under the action of turbulent angular momentum transport, which will try to bring them to solid-body rotation. The result is expected to be the launching of massive neutron rich winds. These, in turn, could produce particularly luminous kilonova counterparts that would be smoking gun evidence for the formation of massive or such objects if detected by future UV/optical/infrared follow ups on gravitational waves events or short gamma-ray bursts.

Ejecta and disk mass (top panel) and black-hole formation time (lower panel) as a function of the tidal parameter Λ for binary neutron star systems. Each point condenses the results of a full 3D simulation.

Neutron star mergers generically result in the ejection of a small fraction (0.1% - 1%) of a solar mass of neutron rich material. As these ejecta expands and cools, they may undergo the so-called r-process and produce heavy nuclei, like gold. The radioactive decay of by-products of the r-process can power an UV/optical/infrared transient known as kilonova.

During the merger, material is ejected from the neutron stars because of tidal torques and shocks. Most of this material is bound and forms an accretion disk around the merger remnant, a massive neutron star or a black hole. A small fraction achieves velocities in excess of the escape velocity and is unbound. Additionally, a substantial fraction of the disk (10% - 40%) becomes unbound over a timescale of a second due to magnetic and neutrino processes in the disk.

We constructed a large database of neutron star merger simulations, the largest to date with fully general-relativistic simulations and with a microphysical treatment of the neutron star equation of state. For each simulation we estimated the dynamical ejecta and the disk masses. We used this data to provide a unified interpretation of gravitational-wave and electromagnetic observations of GW170817. Our results showed that the observation of a kilonova implies a constraint on the presently unknown neutron star equation of state, complementary to that derived by LIGO/Virgo on the basis of the gravitational-wave data alone.

With the aim of studying this process, we performed binary neutron star merger simulations using a temperature- and composition-dependent nuclear equation of state and including neutrino cooling and heating. The focus of this study was to understand the role of weak reactions in shaping the morphology and composition of the outflows and the final abundances.

Final abundances. The HY runs included no weak interactions. The LK simulations included neutrino emission and electron/positron captures. Finally, the M0 simulations included both neutrino emission and absorption, as well as electron/positron captures.

We performed simulations with different microphysics modules in WhiskyTHC switched on or off to be able to measure their relative importance. We used the nuclear reaction network code SkyNet to estimate the final yields of the ejecta for each of these runs.

We found the final abundances of nuclei with atomic mass number A > 120 to be insensitive to weak reactions. On the other hand, the abundances of lighter nuclei are affected in a qualitative and quantitative way by weak reactions. The reason for this is that heavy nuclei are mostly synthesized in the cold and catalyzed part of the ejecta due to tidal interactions, while lighter nuclei are also synthesized in shock- and neutrino-driven outflows, which undergoes weak reprocessing within the timescale of our simulations.

Neutrino-driven convection plays an important role in the explosion of core-collapse supernovae. However, since it results in the creation of small scale flow structures that are difficult to capture numerically, it is also a serious obstacle to the development of accurate supernova simulations. Quantifying these aspects has been the aim of my research.

One of the main finding of our work is that the numerical viscosity inherent to any scheme can interfere with the turbulent cascade in such a way as to trap energy at large scale (the so-called bottleneck effect). This explains why low resolution simulations are artificially more prone to explosion, something that had been already observed before. It also suggests that some peculiar properties of neutrino driven convection reported in the literature might be numerical artifacts. One in particular is the rather flat kinetic energy spectrum observed in simulations that is difficult to explain from a turbulence theory point of view.

Neutrino-driven convection. Compensated velocity power-spectra at different resolutions. As the resolution increases, the slope predicted by Kolmogorov's theory is recovered.

We used the high-order methods implemented in WhiskyTHC to attack this problem and carry out simulations of neutrino-driven convection in core-collapse supernovae at unprecedented resolutions. We found that the predictions from Kolmogorov's classical theory of turbulence are indeed recovered. However, this happens only at extremely high resolutions, 15 times or more higher than current high-resolutions global simulations.

Driven relativistic turbulence. Logarithm of the Lorentz factor minus 1.

WhiskyTHC was used to perform the first systematic study of stationary, isotropic, turbulence in the relativistic regime.

For that study, we performed 4 batches of simulations, labelled "A", "B", "C", and "D". For the A runs we simulated sub-relativistic turbulence, while the B, C, and D simulations had increasingly relativistic flows. We studied the way in which turbulence changes as the flow becomes relativistic using data from all of these simulations, which we performed at multiple resolutions.

Driven relativistic turbulence. Compensated velocity power-spectra for different resolutions (64^3, 128^3, 256^3, and 512^3 grid points) and for increasingly relativistic flows (A, B, C, and D).

We found that turbulence becomes increasingly intermittent as the fluid becomes relativistic. However, somewhat surprisingly, we find low order velocity statistics, such as the power-spectrum, to be in good agreement with the prediction of Kolmogorov's classical theory of turbulence.