Dynamical ejecta from binary neutron star mergers: Impact of a small residual eccentricity and of the equation of state implementation

Dynamical ejecta from binary neutron star mergers: Impact of a small

residual eccentricity and of the equation of state implementation

Francois Foucart ,

Matthew D. Duez ,

Lawrence E. Kidder ,

Harald P. Pfeiffer,

and Mark A. Scheel

Department of Physics and Astronomy,

University of New Hampshire

9 Library Way, Durham, New Hampshire 03824, USA

Department of Physics and Astronomy,

Washington State University

, Pullman, Washington 99164, USA

Cornell Center for Astrophysics and Planetary Science,

Cornell University

, Ithaca, New York 14853, USA

Max Planck Institute for Gravitational Physics (Albert Einstein Institute)

Am M

}

uhlenberg 1, 14476 Potsdam, Germany

TAPIR, Walter Burke Institute for Theoretical Physics, MC 350-17,

California Institute of Technology

Pasadena, California 91125, USA

(Received 1 May 2024; accepted 7 June 2024; published 2 July 2024)

Predicting the properties of the matter ejected during and after a neutron star merger is crucial to our

ability to use electromagnetic observations of these mergers to constrain the masses of the neutron stars, the

equation of state of dense matter, and the role of neutron star mergers in the enrichment of the Universe in

heavy elements. Our ability to reliably provide such predictions is however limited by a broad range of

factors, including the finite resolution of numerical simulations, their treatment of magnetic fields,

neutrinos, and neutrino-matter interactions, and the approximate modeling of the equation of state of dense

matter. To disentangle some of these uncertainties, in this manuscript we focus narrowly on the impact that

two small variations in the setup of numerical simulations can have on the matter ejected during the merger.

Specifically, we perform 4 simulations: one at higher residual eccentricity (

∼

), two with reduced

eccentricity (

∼

003

–

004

), and slightly different equation of state implementations, and a fourth at

higher resolution to estimate numerical errors. We perform these simulations for a single physical

configuration, a

⊙

−

⊙

binary neutron star system using the SFHo equation of state, ignoring

magnetic fields and neutrino effects. We find that a residual eccentricity

∼

, as measured

∼

–

orbits

before merger, causes

–

30%

relative changes in the amount of ejected mass. These differences in

mass ejection appear mainly due to changes in the amount of matter ejected as a result of core bounces

during merger. We note that

residual eccentricities have regularly been used in binary neutron star

merger simulations as proxy for circular binaries, potentially creating an additional source of error in

predictions for the mass of the dynamical ejecta. On the other hand, we find that simulations performed

with slightly different implementations of the same equation of state are consistent within estimated

numerical errors.

DOI:

10.1103/PhysRevD.110.024003

I. INTRODUCTION

Over the last few years, neutron star mergers have

become an important source of information for theoretical

nuclear physics. The observation of the first binary neutron

star merger GW170817, through gravitational waves

[1]

provided us with constraints on the equation of state of

neutron-rich matter above nuclear saturation density

[2]

.In

that respect, gravitational wave observations complement

measurements of neutron star radii by NICER

[3,4]

astrophysical constraints on the maximum mass of neutron

stars

[5]

, as well as laboratory measurements of heavy,

neutron-rich nuclei

[6,7]

. The ejection of neutron-rich

matter during and after merger also provides us with

important information about neutron stars, and their

potential role in enriching the Universe in heavy elements.

Neutron-rich ejecta from binary neutron star and black

hole

–

neutron star mergers undergo

-process nucleosyn-

thesis, and radioactive decay of the ashes of the

-process

power a UV/optical/infrared signal, hours to weeks after the

merger

—

kilonova

. Kilonovae carry information about the

properties of the merging compact objects (including

additional information about the equation of state of dense

matter)

[8,9]

, and about the elements synthesized in the

ejecta

[10]

Extracting information from kilonovae is a complex

problem,requiringmergersimulationsaccountingforgeneral

relativity, magnetic fields, and neutrino-matter interactions

(see, e.g.,

[11

–

17]

), as well as radiation-hydrodynamics

simulation of the evolution of the ejected matter as it is

PHYSICAL REVIEW D

110,

024003 (2024)

2470-0010

2024

110(2)

024003(11)

024003-1

continuously heated by nuclear reactions

[18

–

21]

.Our

current ability to model kilonovae remains limited by both

the accuracy of the simulations as well as the manner in which

they incorporate (often approximately) important nuclear

physics: dense-matter equation of state, neutrino transport,

neutrino-matter interactions, neutrino oscillations, out-of-

equilibrium nuclear reactions, photon transport, and opacities

in the ejecta... While these issues are well known, their actual

impact on our ability to extract information from kilonovae

remains uncertain.

In this manuscript, we begin an investigation of the

robustness of simulation results to changes in the simu-

lation setup and the input nuclear physics. To avoid

confusing the many potential sources of uncertainty affect-

ing neutron star mergers, we begin here with general

relativistic hydrodynamical simulations with neither mag-

netic fields nor neutrinos, and explore the sensitivity of the

inspiral merger and early-matter ejection to small changes

in the initial conditions (binary eccentricity), exact imple-

mentation of the dense-matter equation of state, and

numerical resolution. Previous simulations have already

shown that large eccentricities can lead to the ejection of

much larger amount of mass during tidal disruption of a

neutron star (see, e.g.,

[22,23]

). Our results show that small

changes

in the initial eccentricity of the binary can

impact predictions for the ejected matter well above

estimated numerical errors. This is an important consid-

eration when attempting to predict the mass ejected by any

given binary system, as many existing numerical relativity

simulations have

residual eccentricity in their initial

conditions. In addition, the simulations presented here can

serve as a baseline for further investigation of the impact of

nuclear inputs and/or numerical methods on the main

observables extracted from merger simulations.

II. EQUATIONS OF STATE

All simulations presented in this manuscript are per-

formed with the SFHo equation of state

[24]

. We note

however that different implementations of this equation of

state are currently available, and that some postprocessing

is applied to publicly available equations of state in order to

make them usable by numerical simulations. Specifically,

our numerical simulations require subluminal sound speeds

and a one-to-one map between physical variables (density,

temperature, electron fraction) and evolved variables (con-

served mass density, energy density, momentum density).

We will see that different implementations of the same

equation of state lead to

variations in the fluid

properties, e.g., the pressure. One of the objectives of this

paper is to test whether such variations matter for the

predictions made for the dynamical ejecta. We will show

that for the simulations presented here and within our

estimated numerical errors, the two versions of the SFHo

equation of state used here provide consistent results.

To study the impact of these implementation details, we

consider two distinct tabulated versions of the SFHo

equation of state. First, we start from the table available

in the Compose database with a minimal amount of

postprocessing. We read the pressure

, specific internal

energy

, and sound-speed square

from the table. The

values are tabulated as functions of the logarithm of the

baryon number density

, logarithm of the temperature

and electron fraction

. Whenever a point in the table has

sound-speed square

, we correct the value

by simply taking the average value of all

neighboring points in the table (in 3D; thus from 8

neighboring points). We additionally correct the pressure

to guarantee that

dP=dT

≥

everywhere, through linear

interpolation of the logarithm of

from neighboring point

of the table (in 1D only here at constant density and

electron fraction; thus from 2 neighboring points).

Our second version starts from the stellarCollapse

[25]

table. As for the Compose equation of state table, we

correct the pressure to guarantee that

dP=dT

≥

every-

where. The two equations of state are close to each other

(by the standard of the methods typically used to produce

equations of state matching specific nuclear physics mod-

els), but not exactly identical. Table

shows the gravita-

tional (ADM) mass, baryon mass, and radii of the stars

evolved in this manuscript, while Fig.

shows the mass-

radius relation of cold isolated neutron stars for both

equations of state, the pressure for 2 one-dimensional

slices through the table, and the relative differences in that

pressure. We see differences of

or less between the

two implementations. The two tables also differ in their

definition of the baryon mass at the

level, due to

different choices for the reference mass of a single baryon

(see below). We note that while we will refer to these

equations of state as

“

StellarCollapse

”

and

“

Compose

”

TABLE I. Properties of neutron stars evolved in this manu-

script. The table from which we construct the SFHo equation of

state (EoS) is either from StellarCollapse (SC)

[25]

or Compose

(Comp)

[26,27]

, postprocessed for use in SpEC simulations. We

also quote the ADM mass, baryon mass (defined through

, with

the number of baryons and

the neutron

mass), and circular radius of each neutron star, as well as the

simulations in which they were used (using the notation of

Table

). The circular radius is defined so that

is the

circumference of an isolated neutron star at the equator.

EoS table

ADM

circ

Simulations

Units

⊙

1.30

1.44

11.92

SC-MR/SC-HR

1.40

1.56

11.88

SC-MR/SC-HR

Comp

1.30

1.44

11.94

Comp-no-ecc

Comp

1.40

1.56

11.90

Comp-no-ecc

Comp

1.29

1.42

11.95

Comp-ecc

Comp

1.39

1.55

11.90

Comp-ecc

FOUCART, DUEZ, KIDDER, PFEIFFER, and SCHEEL

PHYS. REV. D

110,

024003 (2024)

024003-2

here, the differences between tables used in our simulations

are likely due in part at least to interpolation from tables

with different grids in density, temperature, and electron

fraction

—

i.e., they may tell us more about errors due to the

finite resolution of the tables than about differences

between StellarCollapse and Compose.

III. INITIAL DATA AND ECCENTRICITY

We generate initial conditions using the Spells initial

data solver

[28,29]

. Spells begins by generating quasie-

quilibrium binaries with zero initial radial velocity (

“

qua-

sicircular

”

binaries, with neutron stars distorted according

to equilibrium tides) for user-provided values of the neutron

star baryon masses and binary separation. This typically

results in orbits with residual eccentricity

∼

, which

can be further reduced through an iterative procedure

solving for the angular and radial velocity of a binary

with the same baryon masses and separation but zero

eccentricity

[30,31]

We note that eccentricity is not a

uniquely defined quantity for inspiraling binaries. Here, we

calculate the eccentricity as in

[31]

; i.e., we measure the

orbital angular velocity

from the first 2.5 orbits of the

trajectory of the center of mass of the neutron stars, and fit

the functional form

Þ¼

−

cos

;

with free parameters

;

. By analogy with

the eccentricity of Newtonian orbits, we define

. Alternative fits that impose

and/or

can also be used and provide equivalent results for

this simple equal mass, nonspinning system.

In this manuscript, all simulations are started

∼

orbits

before merger. We begin by applying the eccentricity

reduction procedure to the StellarCollapse equation of state

table, with neutron star baryon masses

⊙

;

⊙

corresponding to gravitational masses

⊙

;

⊙

for neutron stars in isolation. All neutron stars are initially

nonspinning. This results in an eccentricity

∼

003

in the

simulation continued through merger. We then use the same

baryon masses, angular velocity, and radial velocity for the

Compose equation of state table. We note however that the

two tables use different definitions of the baryon mass.

Taking the Compose definition

, with

the

number of baryons and

the mass of the neutron, the

baryon masses of the StellarCollapse stars are actually

⊙

;

⊙

. Due to this and other small differences

in the equation of state implementations, the Compose

simulation has neutron stars with gravitational masses

⊙

;

⊙

. As gravitational wave emission and

orbital evolution are primarily driven by the gravitational

masses, this results in slightly different orbits, with residual

FIG. 1. Top: mass-radius relationship for the two equations of

state used in this manuscript. For the range of masses used in our

simulations, the radii differ by

∼

m or less. Middle: pressure as

a function of energy density for the same equations of state. We

consider pressure at

MeV and for either the initial

equilibrium

or for fixed

. We show results in units

⊙

, for the range of densities accessed by the

simulations. Bottom: relative error between the two equations of

state, for the same quantities as in the middle panel.

We note that the initial angular and radial velocity of the

binary are the actual free parameters that can be set in our

simulations, not the eccentricity or mean anomaly, which we can

only approximately measure after a few orbits of evolution.

DYNAMICAL EJECTA FROM BINARY NEUTRON STAR

...

PHYS. REV. D

110,

024003 (2024)

024003-3

eccentricity

∼

. While residual eccentricities of that

order are commonly used as initial data for binary neutron

star mergers in the numerical relativity community, we will

see that this results in differences in the merger time and

dynamical mass ejection beyond the estimated numerical

errors. To better differentiate the effect of equation of state

implementation and eccentricity, we thus perform an addi-

tional simulation with the Compose equation of state table,

now keeping the gravitational mass of the neutron stars

constant. This brings the eccentricity down to

∼

004

We will see that this latter simulation is in better agreement

with the simulation performed using the StellarCollapse

equation of state table. We note that matching baryon

masses for equations of state using different definitions of

that quantity was obviously not the right choice to make.

The two simulations performed here with stars of the same

gravitational mass in fact also have nearly the same baryon

mass (see Table

) when the same reference baryon mass

is used consistently for both simulations. This error in

our initial choice of matching initial conditions however

ends up allowing us to test the outcome of simulations with

very similar neutron stars but different eccentricities.

Finally, in order to compare the effect of equation of state

implementation with the expected numerical error, we also

perform the simulation using the StellarCollapse equation

of state table at higher resolution. Specifically, our standard

simulation setup is such that our finite-difference grid

initially has

∼

grid cells across the diameter of each

neutron star, while the high-resolution simulation has

∼

113

grid cells across the diameter of each neutron star. These are

fairly standard resolutions for pure hydrodynamics or

radiation-hydrodynamics simulations, though much

coarser than what would be required to perform accurate

magnetohydrodynamical simulations of these systems.

During the evolution, the grid spacing decreases by as

much as 20% as the grid contracts to follow the inspiral of

the binary (i.e., the resolution used during the evolution is

at most times better than the initial resolution quoted here).

The pseudospectral grid used to evolve Einstein

’

s equations

uses adaptive mesh refinement, with the number of basis

functions in each element chosen to reach a target trunca-

tion error in the spectral expansion

[32,33]

. When going

from standard to high resolution, that target error is

multiplied by

, with

the grid

spacing on the finite-difference grid in our standard

–

high-resolution simulation.

To summarize, we perform 4 simulations in this manu-

script (see also Table

(i) The merger of a system with

003

using the

StellarCollapse version of the SFHo equation of

state, at standard resolution (SC-SR)

(ii) The merger of a system with

003

using the

StellarCollapse version of the SFHo equation of

state, at high resolution (SC-HR)

(iii) The merger of a system with

004

using the

Compose version of the SFHo equation of state at

standard resolution (Comp-no-ecc)

(iv) The merger of a system with

010

using the

Compose version of the SFHo equation of state at

standard resolution (Comp-ecc)

all started 6 orbits before merger with nonspinning neutron

stars. The first three simulations have neutron stars of

gravitational mass

⊙

and

⊙

; the fourth has

gravitational masses

⊙

and

⊙

(see Table

for

other properties of these stars).

We note that we choose these initial conditions because

they allow us to study a wide range of physical processes

observed in neutron star merger simulations. The mass

asymmetry is enough to lead to the production of small

amounts of tidal ejecta,

−

⊙

. The system is of low

enough mass that the remnant does not immediately

collapse to a black hole, leading to the ejection of matter

associated with the collision of the two neutron stars, but

high enough mass that it is expected to collapse on tens of

milliseconds timescales. The same system was for example

evolved to black hole collapse in

[34]

, with 2-moment

neutrino transport. In this manuscript, as we lack any

microphysics in the simulations, we focus solely on the

dynamical ejecta and formation of the early neutron star

remnant, leaving the study of black hole formation for

simulations including some means of modeling angular

momentum transport in the remnant.

IV. NUMERICAL METHODS

We evolve the binary neutron star systems presented here

using the Spectral Einstein Code (SpEC)

[35]

. SpEC

evolves Einstein

’

s equations in the generalized harmonics

formalism

[36]

on a pseudospectral grid, and the general

relativistic equations of hydrodynamics in the Valencia

formalism

[37]

on a separate finite-volume grid.

Specifically, the fluid evolution assumes a perfect fluid

with stress-energy tensor

μν

;

TABLE II. List of simulations performed for this manuscript.

The EoS table is defined as in Table

. The initial grid spacing of

the finite-difference grid

;

and initial eccentricity

are also

provided. Note that in the coordinates of the simulation (quasi-

isotropic coordinate), the radius of each neutron star is only

∼

km; additionally, the grid resolution at later times is always

finer than the initial resolution (by as much as

∼

20%

Sim

EoS table

;

[m]

SC-MR

205

0.003

SC-HR

165

0.003

Comp-no-ecc

Comp

205

0.004

Comp-ecc

Comp

205

0.01

FOUCART, DUEZ, KIDDER, PFEIFFER, and SCHEEL

PHYS. REV. D

110,

024003 (2024)

024003-4

with

the fluid baryon density,

its pressure,

its specific enthalpy,

its specific internal

energy,

its 4-velocity, and

μν

the inverse 4-metric. We

evolve the densitized baryon density

ffiffiffi

;

with

the unit normal perpendicular to a constant time

slice and

the determination of the spatial 3-metric

.We

also evolve the projections of the stress-energy tensor:

ffiffiffi

μν

−

;

−

ffiffiffi

;

and the conserved electron lepton number density

with

the electron fraction of the fluid. The evolution

equations are

∂

∂

Þ¼

;

∂

Þþ

∂

Þ¼

;

∂

∂

ffiffiffi

−

Þ¼

−

ffiffiffi

μν

∇

;

∂

∂

ffiffiffi

Þ¼

ffiffiffi

μν

∂

μν

Here,

ffiffiffiffiffiffiffiffi

−

is the lapse scalar, and the 3-velocity

is defined such that

, with

the Lorentz

factor and

. The equation of state closes this

system of equation by providing us with a prescription for

;T;Y

, with

the temperature.

We use a third-order Runge-Kutta method for the time

discretization, and the WENO5

[38

–

40]

reconstruction and

HLL

[41]

approximate Riemann solver to capture shocks in

the fluid evolution. Details of the methods used to evolve

binary neutron star systems in SpEC can be found in

[33,42]

. In these simulations, we ignore magnetic fields and

neutrinos.

The initial finite-volume grid covers a rectangle

−

;

−

;

−

;

km. The initial

location of the neutron star centers is at

km,

on the

axis, and the rectangle rotates with the binary.

During tidal disruption, a second level of refinement is

created, twice as large and with half the grid resolution.

Finally, after merger, we switch to a cubic grid with four

levels of refinement. The finest level has sides of length

40 km, with each coarser level increasing in size by a factor

of 2. Each refinement level uses

256

cells (

312

at high

resolution), for a grid resolution on the finest grid of 156 m

(128 m at high resolution). We switch from the rectangular

to cubical grid as soon as matter outflows at the boundary

rises above

⊙

As we focus on the production of

dynamical ejecta, we follow all systems to 3.5-ms post-

merger only. Here, the merger time is defined (semiarbi-

trarily) as the time at which the maximum baryon density

on the grid rises 3% above its value at

(due to the

collision of the cores).

We note that here and in the rest of the manuscript,

distances are reported in the coordinates of the simulations,

which evolve according to the generalized harmonic

equations. Distances are thus not coordinate-invariant

quantities here, and are mostly useful to guide our under-

standing of the dynamics of mergers.

V. DEFINITION OF UNBOUND MATERIAL

We measure the unbound mass following the methods

described in

[43]

. For ease of comparison with other work,

we report results using the geodesic criterion

−

, the

Bernoulli criterion

−

∞

−

(for our equation of

state table), as well as the most advanced model of

[43]

that

accounts for heating from

-process nucleosynthesis, cool-

ing from neutrino emission during the

-process, and the

finite amount of time available for bound material to gain

energy from the

-process before it reaches apoastron.

The geodesic criterion typically provides a lower bound

on the mass and kinetic energy of the ejecta, while the

Bernoulli criterion overestimates the ejecta. Indeed, the

geodesic criterion ignores

-process heating and adiabatic

expansion of the fluid, while the Bernoulli criterion ignores

cooling and the finite time available for

-process heating.

The more complex criterion is far from perfect, as it still

ignores nonlocal effects in the evolution of the fluid.

[43]

, we showed that for the one comparison point

available with a full 3D time evolution

[44]

, the latter

criterion predicted the mass and kinetic energy of the ejecta

accurately. It did not however accurately predict its geom-

etry or velocity distribution. We infer from this that our

criterion appears to do well at predicting global quantities,

but is not an appropriate replacement to long-term hydro-

dynamics evolution of the ejecta for more detailed proper-

ties of the outflows.

We also calculate the mass-weighted root-mean-square

average velocity

and mass-weighted average electron

fraction

, using the more complex criterion for unbound

mass. The average electron fraction is not particularly

meaningful for these simulations, however, as we do not

include neutrino-matter interactions here and thus ignore

what should be the main source of variation of the electron

fraction. All simulations find

∈

035

;

037

, but the

correct electron fraction would certainly be higher when

We note that this mass loss occurs over a timescale of a

fraction of a millisecond. In all of our simulations, the mass loss

at the boundary during that time is less than

−

⊙

. The mass

loss is computed by integrating the mass fluid

across a

surface 6 grid cells away from the domain

’

s spatial boundary.

DYNAMICAL EJECTA FROM BINARY NEUTRON STAR

...

PHYS. REV. D

110,

024003 (2024)

024003-5

accounting for weak interactions (e.g.,

[34]

finds

∼

for this system).

VI. RESULTS

At the beginning of our simulations, the main difference

observed between our four cases is due to the residual

eccentricity of these systems (or, equivalently, the differ-

ence in initial gravitational masses of the neutron stars).

The system with

∼

merges about half an orbit

(1.2 ms) later than the other simulations. For comparison,

the merger times of the other three simulations, defined in

SpEC as the time at which the maximum density of the grid

reaches

max

merger

Þ¼

max

, are all within

s of each other. We illustrate these differences in

Fig.

, which shows the coordinate separation between the

two neutron stars as a function of the time to merger.

In all systems, the remnant is a neutron star supported by

differential rotation, surrounded by an accretion disk of a

few percent of a solar mass. These remnants are, in most

respects, very similar to each other. We can for example

look at the properties of the disk (defined as the mass with

cc) about 2.5-ms postmerger, and find that all

simulations show disk masses

disk

∈

–

−

⊙

with average temperature

∈

–

MeV. The

differences in these quantities are largest between the

medium and high-resolution simulations, indicating that

any impact of the eccentricity or equation of state imple-

mentation on these quantities is not resolved. These

quantities would nearly certainly be more significantly

impacted by the introduction of realistic magnetic-field

evolution

[12,45,46]

than by any of the effects studied here.

Specifically, we would expect the rapid growth of magnetic

field in the contact regions between the two neutron stars

due to the Kelvin-Helmholtz instability, with potentially

significant impact on the heating of the remnant, angular

momentum transport in the remnant, and collapse time. We

visualize the evolution of the binary through merger in

Fig.

. We clearly see the heating of the contact regions to

tens of MeV, and the formation of hot tidal arms (a few

MeVs) outside the remnant. We note that what we call

“

disk

”

in these systems is far from the conventional

axisymmetric torus around a compact object. In such a

merger, the system is not expected to become symmetric

until after collapse to a black hole.

The most important output from our simulations that

appears impacted by their residual eccentricities is the

dynamical ejecta, depicted in Figs.

and

. In this system,

unbound matter is continuously ejected by the remnant,

first in the tidal tail, then due to core bounce and possibly

tidal arm instabilities

[47]

. The amount of cold material

from the tidal tail marked as unbound stabilizes about 1 ms

postmerger at

tail

∈

–

−

⊙

. The effect of

eccentricity, of numerical resolution, and of the criterion

used to define unbound matter are, at this step, compa-

rable.

We note that our interpretation of the material as

coming from the

“

tidal tail

”

“

core bounce

”

is largely

based on the very different temperature of the two compo-

nents, as demonstrated in Fig.

The amount of matter marked as unbound after the first

core bounce is most cleanly measured around (2

–

2.5) ms

postmerger (see Fig.

). As this matter is hotter, the choice

of unbound criteria now becomes particularly important

—

the geodesic criteria, which ignore the internal energy of

the fluid, predict significantly less ejecta than the other

criteria. A summary of the unbound mass for each

simulation and each criterion at that time is presented in

Table

III

. We see that especially for criteria that do include

the internal energy of the fluid, uncertainties in the

predicted amount of ejecta are now clearly dominated by

residual eccentricity. We see

∼

–

30%

differences between

the two simulations using the Compose table, with the

004

simulation using the Compose table being closer

to the

003

simulation using the StellarCollapse table

than then

simulation using the Compose table.

A clear difference is also observable in the velocity of

FIG. 2. Coordinate separation between the center of mass of the

neutron star as a function of the time to merger. The high-

resolution SpEC simulation would, on this scale, be indistin-

guishable from the other two low-eccentricity evolution.

The mass quoted here is for the more advanced criterion of

[43]

, as will generally be the case here unless otherwise noted. On

the other hand, the figure shows results for the Bernoulli criterion,

which typically overestimates outflow masses. We use that simple

criterion in the figure because the more advanced criterion was

only computed as a postprocessing step on specific simulation

snapshots, and thus we do not have densely sampled data for that

criterion. We only show densely sampled results for the Stellar-

Collapse table for the same reason: due to the different definition

of the baryon mass in the Compose table, the outflow mass in

those simulations was only computed as a postprocessing step on

a few specific snapshots, shown in the figure.

FOUCART, DUEZ, KIDDER, PFEIFFER, and SCHEEL

PHYS. REV. D

110,

024003 (2024)

024003-6

FIG. 3. Visualization of the baryon density (left) and temperature (right) of simulation Comp-ecc at contact (top, 0.75-ms premerger),

at merger (middle), and 2.5-ms postmerger (bottom).

DYNAMICAL EJECTA FROM BINARY NEUTRON STAR

...

PHYS. REV. D

110,

024003 (2024)

024003-7