Bound Debris Expulsion from Neutron Star Merger Remnants

Many studies have found that neutron star mergers leave a fraction of the stars’ mass in bound orbits surrounding the resulting massive neutron star or black hole. This mass is a site of r-process nucleosynthesis and can generate a wind that contributes to a kilonova. However, comparatively little is known about the dynamics determining its mass or initial structure. Here we begin to investigate these questions, starting with the origin of the disk mass. Using tracer particle as well as discretized fluid data from numerical simulations, we identify where in the neutron stars the debris came from, the paths it takes in order to escape from the neutron stars’ interiors, and the times and locations at which its orbital properties diverge from those of neighboring fluid elements that end up remaining in the merged neutron star.

1. INTRODUCTION On August 17, 2017, the LIGO/Virgo collaboration (LVC) detected the first gravitational wave (GW) signal from a binary neutron star (BNS) merger, the GW170817 event (Abbott et al. 2017a;Villar et al. 2017).This was also the first multi-messenger detection of such an event, as it was seen across the electromagnetic (EM) spectrum (Fermi, Chandra, UV, optical, IR, radio).
From GW170817 and similar BNS mergers, we can learn about many fundamental questions in general relativity, cosmology, and astrophysics (Coulter et al. 2017;Guidorzi et al. 2017;Haggard et al. 2017;Margutti et al. 2017).For example, the distance to this event inferred from GWs was consistent with entirely independent measurements of the Hubble constant (Abbott et al. 2017b).Analysis of the GW170817 GW form led to important constraints on the nuclear equation of state (EOS) (Damour et al. 2012;Radice et al. 2018b;Bernuzzi 2020).Its optical/IR emission (optical/IR analysis papers) gave support to the idea that such mergers could be a main channel of rapid neutron-capture (r-process) nucleosynthesis (e.g, Lattimer & Schramm 1974;Eichler et al. 1989;Arcones & Thielemann 2013;Cowan et al. 2021;Pian 2023).In addition, it also demonstrated that, as predicted many years before, neutron star (NS) mergers can create short γ-ray bursts (see, Eichler et al. 1989;Rosswog & Ramirez-Ruiz 2002;Nakar 2007;Berger 2014;Goldstein et al. 2017) There is a strong theoretical consensus that a portion of the neutron stars' original mass is not held in the remnant, whether it is a massive neutron star or a black hole.Mass can be separated during numerous stages of a BNS merger: the late inspiral phase, the merger stage, and the subsequent evolution of the initially bound debris.Multiple mechanisms can contribute, including tidal gravity, shock heating, neutrino heating, nuclear heating, magneto-centrifugal winds (from both a debris disk and highly-magnetized massive neutron star merger remnant), and magnetohydrodynamic (MHD) turbulent dissipation (see, Metzger et al. 2010;Hotokezaka et al. 2013;Fernández et al. 2015;Radice et al. 2020;Margutti & Chornock 2021;Sarin & Lasky 2021;Most & Quataert 2023;Combi & Siegel 2023;Shibata et al. 2023).Some of the matter that escapes from the remnant is ejected from the system altogether during the late inspiral and the merger itself.This part is called "dynamical ejecta", and it is generally ascribed to a combination of tidal forces and shock heating when the two stars come into direct contact (Hotokezaka et al. 2013).Its mass is rather small, ∼10 −3 M (Radice et al. 2016;Siegel & Metzger 2017;Fujibayashi et al. 2018;Shibata et al. 2023).The dynamical ejecta have received much attention for their possible role in producing a kilonova (KN) event, particularly those features indicating especially high outflow velocities (Metzger et al. 2010; Barnes & Kasen 2013;Radice et al. 2018a;Domoto et al. 2021;Curtis et al. 2023).
However, the same studies predict that the rest of the matter left outside the remnant is, at least temporarily, bound; its mass is rather larger, ∼0.01-0.1M .Here we call this portion the "bound debris".It, too, likely plays a major role in neutron star merger phenomenology.The magnetic flux and confinement pressure needed by relativistic jets may be supplied by this debris.In addition, on rather longer timescales than dynamical ejection, the heat released inside the bound debris by nucleosynthesis and dissipation of MHD turbulence may drive a wind carrying away neutron-enriched isotopes.The photosphere of such a wind could also contribute to kilonova radiation (e.g, Rosswog & Ramirez-Ruiz 2002;Fernández et al. 2015;Just et al. 2015;Korobkin et al. 2020).
Despite its likely importance to the phenomenology of neutron star mergers, there has been little effort to understand the mechanisms governing how much mass is in bound orbits or the nature of its dynamical state, i.e., its distribution of angular momentum and energy: Most et al. (2021) studied the correlations between the angular momentum distribution of debris and other properties, but only for NS-BH mergers.The distribution of mass with angular momentum has special interest because it determines the initial surface density profile of the bound debris, which is almost certainly far from inflow equilibrium.The initial angular momentum distribution, therefore, governs the bound debris' subsequent evolution.
In this paper, we will begin the study of these issues, starting with the most basic questions: From which part(s) of the NSs is the bound debris drawn?How, and at what point in the merger, is the bound debris separated from the majority of the merged NS mass, the part that remains in the newlyformed NS or is carried into the black hole (BH) if one is created?

CALCULATION 2.1. Overview
To place our calculation in context, we first review the stages of a BNS merger (see Fig. 1).Once formed, a NS binary can evolve only by emission of gravitational radiation.During very nearly the entire span of this period, the stars' separation is so large that there is no direct interaction between them, and even gravitational tidal forces are negligible.When, however, their separation shrinks to only a few times their radii, tidal forces begin to distort them; this is the period we call "late inspiral".Eventually, they touch, beginning the "merger proper", which creates a very massive NS (MNS).Depending upon the mass and rotational state of this newly-formed star, it may collapse to a BH (almost immediately or after some delay) or live for an indefinitely long time.In the course of these events, a small fraction of the initial total mass is expelled.
When the mass distribution is substantially asymmetric (the inspiral, merger, and possible collapse phases), there is no simple description of the spacetime.
Consequently, it is necessary to employ numerical relativity methods in order to follow gravitational dynamics.On the other hand, once the mass distribution relaxes to axisymmetry, either as a stable NS or a Kerr BH, there is no need for a further solution of the Einstein field equations.For this reason, we divide our simulation of such events into two pieces.The first part is performed with the numerical relativity plus MHD code IllinoisGRMHD (Etienne et al. 2015;Etienne et al. 2020;Werneck et al. 2023); the second part uses HARM3D+NUC (Murguia-Berthier et al. 2021).A "hand-off" procedure assures consistency between the two in terms of the spacetime and the fluid properties (Lopez Armengol et al. 2022).In the present paper, we make use of data taken only from the first part; future work in this project will employ data from the second part for studies of the bound debris' evolution, including jet-launching, nucleosynthesis, and wind-driving.

Physics treated and equations solved
IllinoisGRMHD solves the equations of MHD, corresponding to the conservation of the baryon number, conservation of Lepton number, conservation of energymomentum, and the two homogeneous Maxwell's equations, respectively.In coordination with the solution of these equations, the BSSN evolution thorn Baikal solves the Einstein Field Equations.In these equations, n b (n e ) is the baryon (lepton) number density and u µ is the fluid four-velocity.The net rate of change in the lepton number R is the difference between the creation rates of electron anti-neutrinos and neutrinos.The cooling rate Q is the rate at which neutrinos (and antineutrinos) carry energy out of the gas.It is calculated from a local emission model and a leakage rate formalism (see Murguia-Berthier et al. 2021;Werneck et al. 2023, for details about the calculation of R and Q).
The energy-momentum tensor is assumed to be that of a perfect fluid plus an EM contribution, where ρ b = m b n b is the baryon density, m b is the mean baryon mass, * F µν = (1/2)˜ µνρσ F ρσ is the dual of the Faraday tensor F µν , and ˜ µνρσ is the Levi-Civita tensor.The parameter h = 1 + + P/ρ b is the specific enthalpy, is the specific internal energy, and P is the fluid pressure, which is given through a tabulated form of the SFHo EOS for hot degenerate matter (O'Connor & Ott 2010).1 Further, b µ = (4π) −1/2 B µ (u) is the rescaled magnetic 4vector in the fluid frame, where Here B i is the magnetic field in the frame normal to the hypersurface, b 2 ≡ b µ b µ is the magnetic energy density, g µν is the spacetime metric, and α is the lapse function.The lapse, the shift vector, and the metric are all defined on spatial hypersurfaces of constant coordinate time t.Ideal MHD (u µ F µ = 0) is assumed throughout.
In addition to solving these equations, we also follow the locations and velocities of a large number of tracer particles.The initial positions of the tracer particles are selected by randomly choosing points p = (x, y, z) within a radius R seed = 75 km from the system center of mass and accepting the particle if where ζ ∈ [0, 1] is a random number.The process is repeated iteratively until all N particles = 50, 000 have been seeded.The fact that the NSs are extremely compact and have central densities many orders of magnitude larger than that of the surrounding gas ensures that the vast majority of the tracer particles are seeded within each NS, as shown in Fig. ??.
Once the positions p n = (x n , y n , z n ) of the particles are known, where n = 0, 1, . . ., N particles − 1, we update them in time using where v i (p n ) is the fluid three-velocity v i = u i /u 0 interpolated to the particle's position p n .The tracer particle positions are updated every 16 local time steps, and the threevelocities are obtained using fourth-order Lagrange interpolation.

Numerical setup
We perform BNS simulations using the Einstein Toolkit (Loffler et al. 2012;Haas et al. 2022) 2 applied to a Cartesian grid with adaptive mesh refinement (AMR) provided by Carpet (Schnetter et al. 2004).Initially, the grid contains eight refinement levels differing by factors of two, such that the resolution at the finest refinement level is ≈185 m.Once the minimum value of the lapse function drops below 0.1, BH formation is assumed to be imminent and two additional refinement levels are added, such that the resolution of the finest refinement level becomes ≈46 m.
At the outer boundary, located at ≈5670 km, we apply radiation boundary conditions to the metric quantities using Newrad (Alcubierre et al. 2003;Loffler et al. 2012) and a simple copy boundary condition for the MHD fields.
The spacetime is evolved using Baikal (Haas et al. 2022), while the MHD fields are evolved using a newly developed version of IllinoisGRMHD (Werneck et al. 2023) that supports finite-temperature, microphysical EOSs, and neutrino physics via a leakage scheme.In addition to the MHD fields, we evolve the entropy by assuming that it is conserved (see e.g., Noble et al. 2009, for a similar strategy), an approximation that, while poor at shocks, allows us to use the entropy as a backup variable during primitive recovery (see Werneck et al. 2023, for more details).
Initial data for an equal-mass binary is constructed using LORENE (Gourgoulhon et al. 2001;Feo et al. 2017;Gourgoulhon et al. Accessed in 2022), with each NS having a baryonic (gravitational) mass of 1.550M (1.348M ; see Table 1).The initial separation of the stars is 45 km and the stellar radii are ≈9.3 km.Each is given a strong poloidal magnetic field (see e.g., Appendix A of Etienne et al. 2015) such that max Where polar coordinates are more easily interpretable, we define such a system through the simple coordinate transformation r 2 = (x 2 + y 2 + z 2 ), cos θ = z/r, tan φ = y/x, and the xy-plane is identical to the initial binary orbital plane.In these coordinates, the polar component of angular momentum is u φ ≡ xu y − yu x becomes a conserved quantity when the spacetime is axisymmetric.

Original location of matter escaping the merger remnant
The most basic distinction between the tracer particles is their post-merger fate: unbound, bound, or confined within either a stable NS or a BH.We define unbound as having hu t < −1, where h is the specific enthalpy of the fluid surrounding the particle and u t is the particle's conserved energy.A bound particle is one with hu t > −1, but located (in this case) outside the innermost stable circular orbit (ISCO) of the BH that forms (for this case, 3.7r g = 14.8 km from the center of the BH).We call the latter two categories "survivors".
The simplest description of the particles' post-merger position is their radial distance from the system center-of-mass, R f ; here "post-merger" means at the end of the simulation, ∼18 ms after its start.One might expect that particles initially placed near the surface of a NS would be more likely to "survive" than particles whose starting position was deeper inside one of the stars.The relation between R f and the distance R i between each survivor particle and the center of its home NS is shown in Figure 2.Only 950 of the particles survive, 1.9% of the original number.Roughly half are found outside 15r g 60km, and ∼20% outside 100km.Almost all those outside 200 km ( 8% of the surviving particles) are unbound.Most importantly for our present purposes, there is no correlation at all between R i and R f except for a slight (and statistically weak) tendency for the unbound particles to originate at larger R i than the bound particles.Note that both R i and R f are coordinate distances defined by , where (x c , y c , z c ) are either the coordinates of the center-of-mass (x c = y c = z c = 0) or the coordinates of the nearest NS center.However, there is a sense in which particles whose initial home was at the surface of a NS are more likely to survive from the remnant.As shown by the upper panel of Figure 3, the probability that a tracer particle escapes grow for R i > 8 km, i.e., within ∼1 km of the surface.Compared to particles with smaller R i this probability is greater by ∼3-4.The reason why Figure 2 shows no correlation is revealed by the bottom panel of Figure 3: the probability density for the initial radius of a surviving particle is considerably larger over the range R i ≈ 2-8 km than for R i 8 km.Even if a large fraction of all the surface particles escapes, while only a small fraction of interior particles do, there is so little mass near the surface compared to the mass in the interior that the great majority of surviving particles originate in the deep interior of a NS.

Escape of debris from the NSs' interiors
Although answering the first question we raised, this result prompts a second: how does the bulk of the debris find its way from the deep interior of a NS to an orbit well outside the merged NS or BH? Put another way, how is it that ∼2% of the matter in the interior of a NS is somehow separated from its immediate neighbors and given enough energy and angular momentum to be released from the merged NS?
In order to approach this question, we begin by identifying the time from the beginning of the simulation at which this separation occurs.We do so by choosing a small sample of surviving particles and then pairing each with the closest non-surviving particle in the initial state of the system.These were selected by finding survivors for which R * − 2 km < R i < R * in a single neutron star and were all located in different quadrants.We also constructed comparison samples varying these choices, but the results were quite similar.
In Figure 4, we display how the positions of these particles projected onto the x-y (orbital) plane evolve from t = 13.0 ms, 0.2 ms before the two stars first touch, to 17.62 ms, 1 ms after the black hole forms.The oscillations reflect the rotation of the NSs, first as they orbit one another, and later as they rotate as an asymmetric merged object.As this figure shows, four of the six pairs begin to separate at 14.5 ms, while the other two do not substantially part ways until 15.5 ms.This small sample suggests the critical time for distinguishing surviving particles from captured particles is 14.5-15.5 ms; as we will see shortly, this suggestion is confirmed when we examine the data more closely.
To see what is happening during this crucial millisecond, we present images of the density and Y e maps at that time (Fig. 5).The first of the density images (t = 14.66 ms) shows the system 0.5 ms after the star's first touch.Over the course of the following ∼2 ms, the configuration of the NSs changes from one in which the two stars remain partially distinguishable but are in contact over a sizable surface to one in which most of the stars' mass has taken on a much rounder shape, but a noticeable amount of gas has been expelled into spiral arms.
Complementary images of the specific internal energy density are shown in Figure 6.For these, we begin at a slightly earlier time (t = 14.50 ms) and follow it with the same three times shown in Fig. 5.The internal energy images make shock locations easily visible.In particular, at the two earlier times, they show clearly the shock at the contact surface, while at the two later times they who how it weakens.They further point out shocks in surrounding gas driven by both polar outflow and the rotating asymmetric merged star (Hotokezaka et al. 2013).
Thus, the moment at which the paths of escaping and captured particles diverge is the moment when the two stars are in direct contact, but have not yet combined their masses into a single (in this case, short-lived) NS.During the first part of this stage, the mass distribution is not only highly asymmetric; it also rotates around the center of mass, so spacetime is both asymmetric and time-dependent.However, during the latter part, the mass distribution becomes significantly more axisymmetric, weakening the azimuthal component of gravitational acceleration.Throughout this period, as shown by the numerous narrow regions of elevated internal energy per unit mass in Figure 6, there are many shocks, where kinetic energy is dissipated into heat, and deflection can exchange angular momentum between different fluid elements.
To understand more deeply why the surviving particles diverge from their neighboring captured particles, we now take a closer look at the six surviving particles in our sample.Figure 7 displays their x-y positions at t = 14.50 ms in relation to the fluid density in the orbital plane.Four of the six particles' locations trace a curve in the x-y plane running along the outer edge of the merged star; the other two are found along the "seam-line" between the two merging stars.Their z-coordinates inform us that these particles are typically ∼ (0.1 − 0.3)R * out of the orbital plane, where R * is the radius of the original NSs.
To interpret this 3D geometry, we first plot the x-y locations of all the tracer particles at this same time (blue dots in the top panel of Fig. 8).The four particles in our six-member sample found along the edge of the merged star are all found in the lower hook of the "S"-curve traced out by the entire ensemble of particles; the densest parts of both hooks coincide with the plumes of gas density visible in Figure 7.The other two are in the central bar of the "S"-curve, which coincides closely with the seam-line.It, therefore, appears that particles escape when they move outward through the seam region and then reach the surface.The curved portions of the "S"-curve result from the time sequence at which the particles crossed through the surface of the merged NS.Note that the orbital plane velocities of the survivor particles are primarily azimuthal, and the radial gradient in their speeds indicates that they acquire greater angular momentum as they move outward.
Further information about vertical structure can be derived from examining the locations of the particles within the bar at t = 14.50 ms.Their orbital-plane coordinates all lie along the linear feature in the "S", so that their 3D positions all lie in a vertical plane defined by the line feature and the z−axis (the simple geometry of this structure is likely due to our choice of an equal-mass merger).The bottom panel of Figure 8 shows their distribution in this plane.All of the particles lie a short distance inside the outline of the merged star, a figure that is flattened in the vertical direction (see the images in the third column of Fig. 5).Moreover, they are weakly concentrated toward the orbital plane and tend to avoid the polar axis.The evolution of this structure in time is illustrated by the difference between the figures outlined by the blue and gold dots in the upper panel of Figure 8.At t = 14.50 ms, numerous particles can be found in the plane of the inter-star seam, while the "hooks", taken together, extend over ∼ 3π/2 in az-imuthal angle.Only 0.2 ms later, hardly any particles remain in the seam, while the hooks now cover all 2π in azimuth.Thus, the time period we have singled out as critical in terms of separating surviving particles from initially neighboring captured particles coincides with the time when the expulsion of matter from the seam region is completed.
The variety of paths these particles follow in order is illustrated in Figure 9, which shows three particle trajectories during the ∼ 2 ms spanning the time of first contact between the stars.Just before contact, the tracer particles follow their star's orbit, but moving in and out relative to the system center-of-mass as the star rotates.Once the two stars come into contact (at t = 13.2 ms), the character of their motion changes: now they are part of the merged NS, partaking in its rotation around the system center of mass.However, the axial fluid motions within the contact surface permit particles  Their z-coordinates are z = 2.558, −2.225, −3.097, 2.182, 1.738, 3.504 km for the particles shown as a plus-sign, a square, a circle, a dot, an 'X', and a caret, respectively.
to travel both outward and inward; the red path in Figure 9 is an example of the latter option.Large amplitude motions within this zone are expected because the strong shear at the boundary between the two stars excites turbulence through the Kelvin-Helmholtz instability (Rasio & Shapiro 1999;Kiuchi et al. 2014;Aguilera-Miret et al. 2022).By t = 14.7 ms, only a few particles remain in the seam-plane.Instead, the particles destined to survive that had been inside the merger remnant at 14.5 ms are almost all found in a shell surrounding the remnant.The survivor particles that at 14.5 ms had been in the hooks of the "S" are now in either of two places: roughly half are in the shell surrounding the remnant, many of them having moved inward to do so; the other half are outside the boundary of the figure.

Separation of surviving particles from neighboring particles
The surviving particles are 2% of the NSs' mass; as can be seen in Figure 8, 90% have reached positions just outside the star by ∼ 1 ms after contact.From the data illustrated in Figure 4, it is clear that the captured particles that began close to the survivors are still close to them at this time.It immediately follows that the distinction between survivors and the captured matter is accomplished in two stages.First, a few percent of the total mass of NSs is conveyed to locations outside the stars by motion through the contact surface.Second, when this matter reaches locations immediately outside the merged star, local effects divide this mass into two subsets, one that escapes, the other remaining with the merged star.How this distinction is made is portrayed in Figure 10.For the first ∼ 2 ms shown, the angular momenta of the surviving particles are very similar to those of their matching captured particles; both oscillate as a result of the NSs' rotation.At around 14 − 15 ms, the angular momentum of surviving particles begins to rise (as shown by the arrows in Fig. 8), while the angular momentum of the neighboring particles that remain with the star for the long run do not change appreciably.This growth in survivor particle angular momentum finishes shortly before collapse to a BH (t = 16.6 ms).By then, the surviving particles' angular momentum is ∼ 2 − 10× the angular momentum of the captured particles that started out from locations very near them.
That these results are characteristic of all the surviving particles are shown in Figure 11.A small number of particles have large excursions in u φ during the inspiral, but the over- whelming majority of particles follow the pattern of the small sample we have chosen for a particular study.
The torques creating the change in angular momentum could be the result of non-axisymmetric gravity (Hotokezaka et al. 2013;Radice et al. 2018b;Fujibayashi et al. 2018;Shibata et al. 2023), smooth pressure gradients, or shocks (Hotokezaka et al. 2011(Hotokezaka et al. , 2013;;Sekiguchi et al. 2016).If the magnetic field is as strong as Kiuchi et al. (2014) found (but not seen here), magnetic forces might also contribute.Shocks can be recognized by the sharp steps in angular momentum evident in Figure 10 as well as the cusps seen in the trajectories of Figure 9, which might be associated with shock features visible in the later images of Figure 6.Both gravitational acceleration and smooth pressure gradients could lead to smooth portions of trajectories.The short timescale smaller fluctuations in u φ shown in Figure 10 are more likely associated with pressure gradients than with gravity because there are many short length-scale density gradients, whereas the gravitational field is the result of the global stress-energy distribution.The general upward trend could be associated with non-axisymmetric gravity, and its end at t ≈ 16 ms roughly coincides with the axisymmetrizing of the merged NS (see Fig. 5).On the other hand, it is unclear how a gravitational effect would distinguish the particles that ultimately survive from the nearby particles that are ultimately captured; non-axisymmetric gravity could become important only after other dynamics place fluid elements in the regions where gravitational torques exist.In this context, it is also worth noting that the captured particles show much less short timescale fluctuation in u φ , suggesting that the regions through which they travel have much smoother density and pressure distributions.This would, of course, be a natural consequence of staying closer to the merged star.The electron fraction Y e , the ratio of the number of protons n p to the total number of nucleons, is one of the key quantities for determining the elements synthesized over longer timescales in the bound debris.Through this influence, it is also the key to determining the opacity of the optical/UV photosphere where the KN emission is radiated (see, Barnes & Kasen 2013;Hotokezaka et al. 2013;Siegel & Metzger 2017;Radice et al. 2018a).
If there were special locations in the NSs that contributed disproportionately to the bound debris, the distribution of mass with electron fraction dM/dY e in its initial state could be quite different from dM/dY e in the NSs as a whole.However, we have shown that, in rough terms, the probability that a fluid element inside one of the NSs is expelled from the remnant is very nearly constant for all radii between 1 and 8 km, a region containing the great majority of mass in the star.There are only two exceptions to this constant probability.One is a spherical shell 1 km thick at the outside of the star, where the probability of escape rises to a few times greater than elsewhere.The other is a sphere of radius ∼ 1 km at the center of the star, where the escape probability is 0.1× the value in the star's bulk.Consequently, dM/dY e in the pre-merger NSs should be a very good predictor of dM/dY e in the bound debris immediately after it leaves the merged NS.
On the other hand, it should also be borne in mind that Y e can be affected by events happening during the merger and while the bound debris travel from the merged NS to its place in the surrounding disk.Shock heating in the merger and during the subsequent evolution of the merger remnant can raise the matter temperature above 10MeV.In such a hightemperature environment, electron-positron pair creation can substantially raise the rate of the reaction e + + n → p + + νe , thereby increasing Y e .Neutrino irradiation, within the merged star, as the gas travels to the disk, and in the disk can also raise Y e via reactions like ν e +n → p + +e − .The degree to which any of these mechanisms may or may not alter Y e requires further investigation (Shibata et al. 2023).

Parameter-dependence
Quantitative details in the single calculation on which we base our results could change for different values of the principal system parameters, the mass ratio q and the total mass M .In particular, the bilateral symmetry so evident in the images of Figures 5, 7, and 8 will be broken when q = 1.It has also long been known (Oechslin & Janka 2006;Shibata et al. 2023) that the total mass removed from the merger remnant increases for smaller values of q, so the total amount of bound debris mass in this simulation is on the lower end of what might be expected for a merger of this total mass.
The results should also depend on M .If it is so large that the merger results in immediate collapse to a BH, the amount of expelled mass, both in dynamical ejecta and bound debris, would be much smaller.Moreover, because the tidal gravity contribution to the expelled mass would be a larger fraction of the total, the bound fraction would likely be considerably smaller.On the other hand, if M were small enough for the  merged NS to survive for a longer time, the expelled mass could be rather larger.
The combination of different M and different q can alter the results quantitatively in several ways.Smaller values of both imply particularly small mass for the lessmassive star, giving it both lower mean density and a less centrally-concentrated internal density profile (see, Lattimer & Prakash 2007;Bauswein et al. 2012;Shibata et al. 2021;Raithel et al. 2022).Relatively large tidal mass-loss would likely follow.Larger M and smaller q could enhance the contribution of the contact shock if the larger mass star has 2M because for such large masses, the radius of the NS diminishes, whereas it hardly changes (for most EOSs) for masses between ∼ 0.7 and ∼ 2M (Kashyap et al. 2022).The sense in which these changes tip the balance between bound and unbound debris requires a detailed calculation to determine.

Uncertainty in the equation of state
Although the gravitational waveform measured in the event GW170817 has helped prune the long list of proposed nuclear EOSs (see, Hotokezaka et al. 2011;Bauswein et al. 2013;Radice et al. 2018b), there still remain many candidates.Our quantitative results could be altered in several ways if a different EOS applies.First, the internal density profile of an equilibrium NS is sensitive to the EOS, and that could alter the probability of escape as a function of radius.Second, when the two stars come into contact, the pressure of the shocked region between the two depends on the EOS.Stiffer EOSs support larger NSs that strike each other with lower velocity, diminishing the post-shock temperature (Hotokezaka et al. 2011); on the other hand, a larger ratio between pressure and internal energy density would make the shocked region wider and enlarge the pressure contrast between the shocked region and the surface of the merged star.Third, the EOS governs both the critical mass dividing quick collapse from a longer-lived merger NS and the time to collapse if rapid collapse is, indeed, the result (Bauswein et al. 2013;Radice et al. 2018b;Margalit et al. 2022).A longer lifetime for the merged NS likely leads to more escaping mass, and the ratio of bound to unbound mass could change over this span of time.Thus, different EOSs could lead to changes in both the total mass removed from the remnant and how it is divided into bound and unbound portions.

CONCLUSIONS
In much of the literature on NS mergers, when matter that escapes being confined within the merger remnant is dis-cussed, the focus is solely on the unbound fraction.This fraction is then further divided into two categories, "dynamical" and "secular" (sometimes called "post-merger") (Metzger et al. 2010;Fernández et al. 2015;Siegel & Metzger 2017;Fujibayashi et al. 2018;Radice et al. 2018a)."Dynamical" refers to the material expelled during the merger itself, whether by tidal forces, shocks heating gas, or neutrinos that heat or push the gas."Secular" refers to matter driven off the disk of bound debris many milliseconds after the merger.
However, the majority of the matter outside the merger remnant remains bound to the system, forming the disk that is the source of the secular ejecta and may also support any jet that forms.Hitherto, the mechanisms by which this material leaves the merging neutron stars and is deposited in the debris disk have gone largely unexamined.Here we have determined the origin of the bound debris within the NSs and examined how its fate is distinguished from that of matter remaining inside the merged star.
In particular, by making use of the complementary information provided by Eulerian fluid data and Lagrangian tracer particle data, we have uncovered several surprising facts about the process by which matter is expelled from merging neutron stars.
• The probability that a given parcel of mass leaves the merger remnant is almost, but not quite, independent of its location inside the pre-merger NSs.Matter placed near the surface of a NS is considerably more likely than most matter to be expelled, but there is so much more matter in the bulk of the star that, despite its lower escape probability, nearly all the escaping matter, whether bound or unbound, has its origin deep in a star's interior.
• Turbulent motions within the "seam" surface formed when the two stars touch can mix matter from a wide range of initial depths within the stars.These same motions transport some of this matter through the seam plane to the surface of the merged star.
• A significant amount of matter is brought out to the surface of the merged NS but does not escape.Which matter is captured and which escapes depends on details of the shock structure in the region where the matter that has been squeezed outward by the shock between the two merging stars reaches the surface.
• The surviving matter acquires significant additional angular momentum as it leaves the remnant.It is likely that torques due to the non-axisymmetric gravity of the merged NS, smooth pressure gradients, and shocks all contribute.
• All of these processes occur within a few milliseconds after the stars first touch.
Thus, we have shown that the process by which matter escapes the remnant created by a neutron star merger is considerably more complex and subtle than previously envisioned.Because the disk of bound matter circulating around the merged NS or BH is far from inflow equilibrium, its subsequent evolution depends on what sort of non-equilibrium structure it is given.Understanding the mechanisms regulating this structure is, therefore, an essential part of understanding many of the dramatic observational properties of NS mergers-jets, kNe, and r−process nucleosynthesis.

Figure 1 .
Figure 1.Schematic presentation of the different stages and outcomes of BNS mergers: from left to right, inspiral, merger, post-merger remnant plus debris disk; and, possibly, a stable NS.Timescales for each stage are shown below the corresponding sketches.The bottom line lists the GW frequency and the principal physical mechanisms of each stage.

Figure 2 .
Figure 2. The final distance R f (from the system center-of-mass) of the surviving tracer particles plotted against their initial distance Ri from the center of their home neutron star.Red circles distinguish unbound particles.

Figure 3 .
Figure 3. Upper panel: The probability that a particle with initial radius Ri survives black hole collapse.Lower panel: The probability density (in arbitrary units) that a surviving particle had initial radius Ri.

Figure 4 .
Figure4.The evolution of the separation between six surviving particles and their initial close neighbors.The top panel shows how their y-coordinates evolve, the bottom panel shows the evolution of their x-coordinates.Each pair's starting location (X, Y ) is shifted by an arbitrary amount relative to the other pairs so that they can be visually distinguished.In some cases (e.g., the bottom two), the members of the pair start out so close together that their curves are superposed.

Figure 5 .
Figure 5. Evolution of the electron fraction, xy-plane mass density, and xz-plane mass density.Ye is shown on a linear color scale (xy-plane),the mass densities on a logarithmic scale.Each quantity is presented at times 14.66 ms, 15.29 ms, 16.31 ms, and 17.25 ms; the last is shortly after the collapse to a BH.Note that the colorscale for Ye at t = 14.66 ms is different from the colorscale at the other times.

Figure 6 .
Figure 6.Evolution of the specific internal energy int in the xyand xz-planes, shown on a logarithmic scale at times 14.50 ms, 14.66 ms, 15.37 ms; 16.39 ms.

Figure 7 .
Figure 7.The x-y locations of the same six surviving particles shown in Fig. 4 superposed on the t = 14.50ms density map.Their z-coordinates are z = 2.558, −2.225, −3.097, 2.182, 1.738, 3.504 km for the particles shown as a plus-sign, a square, a circle, a dot, an 'X', and a caret, respectively.

Figure 8 .
Figure 8. (Top): Position in the xy-plane of the surviving particles at t = 14.506 msec (blue dots) and at t = 14.705 msec (yellow crosses).The line labeled "lx" is fit to the line-like feature in the middle of the "S"-shape at t = 14.50 msec.The red arrows illustrate the xy-velocity of a random sample of surviving particles at t = 14.50 ms; the arrows' lengths are proportional to the speeds of the particles' orbital plane motion.(Bottom): Position in the lxz-plane of the surviving particles in the line-like feature at t = 14.50 ms.

Figure 9 .
Figure 9. Motion in (x, y) of three surviving particles (blue, green, and red dotted paths) from t = 12 ms to t = 16.6 ms, when the star collapses to a BH.The black circle shows the outline of their parent star at t = 12 ms.Cusps are the signature of shock deflection.

Figure 10 .
Figure10.Time-dependence of the angular momentum u φ for our sample of six pairs of particles.Survivors are shown with solid lines, and captured particles with dashed lines; matching colors identify the pairs.The time of collapse to a BH, 16.6 ms is marked with a gray dashed vertical line; for captured particles after that time, the physical significance of u φ is unclear.

Figure 11 .
Figure11.Evolution of u φ , the conserved angular momentum component in an axisymmetric spacetime, for a large sample of the survivor tracer particles: 220 particles selected at random.The small number of surviving particles with higher u φ and noisy u φ histories (blue and green lines) were initially located very close to the surface of one or the other of the NSs.

Table 1 .
Simulation parameters for an equal-mass, magnetized BNS merger performed with IllinoisGRMHD using a microphysical, finite-temperature EOS and a neutrino leakage scheme.