GRMHD simulations of neutron-star mergers with weak interactions: r-process nucleosynthesis and electromagnetic signatures of dynamical ejecta

Fast material ejected dynamically over $<10$ ms during the merger of a binary neutron-star (BNS) system can give rise to distinctive electromagnetic counterparts to the system's gravitational-wave emission that can serve as a"smoking gun"to distinguish between a BNS and a NS-black-hole merger. We present novel ab-initio modeling of the associated kilonova precursor and kilonova afterglow based on three-dimensional general-relativistic magneto-hydrodynamic simulations of BNS mergers with tabulated, composition-dependent, finite-temperature equations of state (EOSs), weak interactions, and approximate neutrino transport. We analyze dynamical mass ejection from 1.35-1.35Msun binaries, typical of Galactic double-NS systems and consistent with properties of the first observed BNS merger GW170817, using three nuclear EOSs that span the range of allowed compactness. Nuclear reaction network calculations yield a robust 2nd-to-3rd-peak r-process. We find few 1e-6Msun of fast ($v>0.6$c) ejecta that give rise to broad-band synchrotron emission on ~yr timescales, consistent with recent tentative evidence for excess X-ray/radio emission following GW170817. We find 2e-5Msun of free neutrons that power a kilonova precursor on<h timescale. A boost in early UV/optical brightness by a factor of a few due to previously neglected relativistic effects, with appreciable enhancements up to 10h post-merger, provides promising prospects for future detection with UV/optical telescopes such as Swift or ULTRASAT out to 250Mpc. We find that a recently predicted opacity boost due to highly ionized lanthanides at ~70000K is unlikely to affect the early kilonova lightcurve based on the obtained ejecta structures. Azimuthal inhomogeneities in dynamical ejecta composition for soft EOSs found here ("lanthanide/actinide pockets") may have observable consequences for both early kilonova and late-time nebular emission.


INTRODUCTION
The first detection of gravitational waves from the merger of a binary neutron star (BNS) system, called GW170817 (Abbott et al. 2017a), was associated with fireworks of electromagnetic (EM) counterparts across the electromagnetic spectrum. The largest EM followup campaign ever conducted (Abbott et al. 2017b) revealed an associated short gamma-ray burst (GRB) GRB170817, followed by broad-band GRB afterglow emission, as well as the quasi-thermal counterpart AT2017gfo-the first unambiguous detection of a kilonova (Li & Paczynski 1998;Metzger et al. 2010; Metzger 2020), powered by radioac- * CITA National Fellow tive heating from the production of heavy elements via the rapid neutron-capture process (r-process; Burbidge et al. 1957;Cameron 1957) in dense, neutron-rich plasma ejected from the merger site.
The GRB afterglow of GW170817 showed several interesting and unusual properties and has been followed up for several years until today. The faintness of the prompt GRB gamma-ray emission given the proximity of only 40 Mpc, the unusual rising X-ray flux appearing after nine days (Troja et al. 2017), peaking at ≈ 160 days after the merger, followed by a steep decay (Haggard et al. 2017;Hallinan et al. 2017;D'Avanzo et al. 2018;Lyman et al. 2018;Margutti et al. 2017;Mooley et al. 2018;Troja et al. 2018;Ruan et al. 2018), as well as the observation of superluminal motion of the radio centroid, are interpreted as a structured jet ex-panding into the interstellar medium (ISM) viewed at an angle of θ obs ∼ 20 o − 30 o from the jet core (Lamb & Kobayashi 2017;Lamb et al. 2018;Alexander et al. 2017;Hotokezaka et al. 2018;Wu & MacFadyen 2019;Fong et al. 2017;Ghirlanda et al. 2019;Hajela et al. 2019;Lamb et al. 2019;Troja et al. 2020). After several years, the broad-band synchrotron afterglow originating in the decelerating external shock driven into the ISM has faded to levels at which other possible emission components may be revealed, with first observational indications (Hajela et al. 2022;Troja et al. 2022). One robust additional emission component expected on timescales of years after the event is the broad-band synchrotron emission from a transrelativistic shock of fast dynamical ejecta expanding into the ISM (e.g., Nakar & Piran 2011;Hotokezaka et al. 2018;Nedora et al. 2021a;Balasubramanian et al. 2022), which we refer to here as the 'kilonova afterglow'. Late-time fall-back accretion onto a remnant black hole offers a possible alternative explanation to a rebrightening ). The existence of a rebrightening in GW170817 is still a matter of debate Balasubramanian et al. 2022).
The astrophysical sites that give rise to the synthesis of heavy elements via the r-process remain a topic of active debate (see Horowitz et al. 2019;Cowan et al. 2021;Siegel 2022 for recent reviews). Several lines of evidence, including measurements of radioactive isotopes in deposits on the deep seafloor (e.g., Wallner et al. 2015;, abundances of metal-poor stars formed in the smallest dwarf galaxies (e.g., Ji et al. 2016;Tsujimoto et al. 2017), in the halo of the Milky Way (Macias & Ramirez-Ruiz 2018), as well as in globular clusters (Zevin et al. 2019;Kirby et al. 2020), point to a low-rate, high-yield source, much rarer than ordinary core-collapse supernovae both in early and recent Galactic history (see Fig. 3 of Siegel 2019 for a compilation of observational constraints). From the most promising candidates, which include neutron-star mergers (Lattimer & Schramm 1974;Symbalisty & Schramm 1982), rare types of core-collapse supernovae producing rapidly spinning, strongly magnetized proto-neutron stars (magnetorotational supernovae; Winteler et al. 2012;Symbalisty et al. 1985;Nishimura et al. 2006), and hyper-accreting black holes ("collapsars"; Pruet et al. 2003;Surman et al. 2006;Siegel et al. 2019Siegel et al. , 2021; see also Grichener & Soker 2019), definitive evidence for the production of r-process elements has only been obtained in the case of neutron-star mergers through the GW170817 kilonova so far.
Several mechanisms can give rise to the ejection of neutron-rich material conducive to undergoing r-process nucleosynthesis and powering kilonovae in BNS mergers 1 . In the dynamical merger phase (typically lasting ∼ < 10 ms), unbound matter is ejected from the system by tidal forces (sometimes forming 'tidal tails'; Rosswog et al. 1999;Rosswog 2013;Korobkin et al. 2012;Goriely et al. 2011), by shocks generated at the collision interface of the two stars (Hotokezaka et al. 2013a;Oechslin et al. 2007;Sekiguchi et al. 2016), and by a sequence of quasi-radial oscillations ('bounces'; Bauswein et al. 2013b;Radice et al. 2018) of a double-core structure of the remnant that forms immediately upon merger (see Fig. 1; 'dynamical ejecta').
After merger, a remnant NS may be formed, temporarily supported from collapsing into a black hole by solid-body or differential rotation (Duez et al. 2006a;Kaplan et al. 2014;Siegel et al. 2014), which can unbind material in a combination of neutrino-driven (Dessart et al. 2009;Perego et al. 2014;Desai et al. 2022) and magnetically driven winds (Siegel et al. 2014;Ciolfi & Kalinani 2020;Metzger et al. 2018;Mösta et al. 2020;Curtis et al. 2021) over the lifetime and Kelvin-Helmholtz cooling time of the remnant (∼ 10 ms − 1 s). Bound material surrounding the remnant circularizes and forms a neutrino-cooled accretion disk (Fernández & Metzger 2013;Just et al. 2015;Siegel & Metzger 2017;Lippuner et al. 2017;Fujibayashi et al. 2018;Fernández et al. 2019;Miller et al. 2019a;Nedora et al. 2021b) with typical masses of 0.01-0.3 M . A combination of heating-cooling imbalance in the disk corona in the presence of MHD turbulence (Siegel & Metzger 2017 at early times and viscously-driven outflows at later evolution stages (Fernández et al. 2019), together with nuclear binding energy release as seed particles for the rprocess form, can eject up to ∼ 30 − 40% of the initial disk mass.
The GW170817 kilonova was likely generated by a combination of ≈ 0.01 M of dynamical ejecta and winds from the remnant NS for the early blue emission, and by ≈ 0.05 M post-merger disk ejecta for the red emission (see, e.g., Siegel 2019; Metzger 2020; Radice et al. 2020;Margutti & Chornock 2021;Pian 2021 for reviews overseeing the interpretation of the event). The likely dominance of disk outflows over other ejecta channels may extend to the population of merging BNS systems Figure 1. Volume renderings of rest-mass density layers from our APR4 BNS simulation illustrating the different merger stages and associated ejecta components. From left to right: (i) The two neutron stars reach high angular speeds before the collision and tidally deform; (ii) At merger, the neutron star cores bounce off of each other, generating shock-heated material that is ejected at high speeds, represented here as a red density layer; (iii) Due to tidal forces during merger, 'tidal tails' develop that eject cold and highly neutron-rich material, represented here by a blue density layer; (iv) After merger much of the debris material circularizes around the remnant neutron star to form an accretion disk. Winds launched from this disk eject a significant fraction of the total unbound material of this BNS merger. The hot remnant neutron star additionally launches a neutrino and magnetically driven wind in the polar direction, highlighted here by a green density layer.
in general if these systems follow a distribution similar to the Galactic double-neutron star systems observable as radio pulsars (Siegel 2022). The population of merging BNS systems is a matter of active debate (e.g., Landry et al. 2020) and future combined gravitational wave and EM observations will be necessary to infer the observed population. In a population sense, the contribution of neutron-star-black-hole (NS-BH) systems to Galactic rprocess nucleosynthesis is likely by far subdominant relative to BNS mergers (Chen et al. 2021).
Although likely subdominant with respect to other ejecta channels in many scenarios, dynamical ejecta can give rise to distinctive EM signals associated with its most energetic component. The outermost ∼ 10 −8 M of dynamical ejecta might reach ultra-relativistic velocities and serve as a breakout medium for a relativistic jet and give rise to prompt gamma-ray emission similar to GRB 170817A as envisioned by Beloborodov et al. (2020). The transrelativistic tail of ejecta with velocity ∼ > (0.5 − 0.6)c drives a shock into the ISM, giving rise to the aforementioned broad-band 'kilonova afterglow'. Furthermore, this transrelativistic component typically has a sufficiently rapid expansion timescale (< 5 ms) to produce free neutrons, preventing rapid neutron capture to occur in these fastest layers. Free-neutron decay on a timescale of ≈ 15 min acts with a specific heating rate roughly an order of magnitude larger than typical r-process heating at a similar epoch in the outermost layers of the ejecta, giving rise to a short-lived thermal UV/optical kilonova "precursor" transient on a timescale of ∼ < 1 hr after merger (Kulkarni 2005;Metzger et al. 2015;Ishii et al. 2018). Such early blue emission might overlap (and be confused) with additional heating of the outermost ejecta layers by a cocoon of hot shock-heated material surrounding a jet propagating into the ejecta envelope (e.g., Kasliwal et al. 2017;Gottlieb et al. 2018). The aforementioned signatures, if detected early on after merger in future gravitationalwave events, will help to distinguish between a BNS and an NS-BH merger, a distinction that is presently not possible based on gravitational-wave data alone (except when exclusions can be made based on plausible ranges for inferred component masses). This motivates a detailed investigation of the physical mechanisms related to dynamical ejecta, a characterization of its properties, and self-consistent modeling of its unique EM signatures based on first-principle numerical simulations, as attempted here.
In this paper, we present a set of general-relativistic magnetohydrodynamic (GRMHD) simulations of equalmass BNSs typical of the Galactic double-neutron star systems and consistent with the inferred source parameters of GW170817. These simulations self-consistently combine several physical ingredients from inspiral to post-merger, including magnetic fields, tabulated (finite temperature, composition-dependent) equations of state (EOS), weak interactions, and approximate neutrino transport via a ray-by-ray scheme. We employ these simulations to explore in detail the physical ejection mechanisms of dynamical ejecta, the nucleosynthesis it gives rise to, and its distinctive EM signatures including neutron precursor emission and the kilonova afterglow. We focus here on building novel, self-consistent models for the EM signatures that are directly based on ab-initio numerical simulations of BNS mergers.
The paper is organized as follows. In Sec. 2, we describe our simulation setup, including numerical methods, initial data, and simulation diagnostics. In Sec. 3, we discuss simulation results on dynamical ejecta, provide a comprehensive discussion of ejection mechanisms, and conclude with a brief description of the post-merger phase. Section 4 reflects on nucleosynthesis results of the dynamical ejecta. In Sec. 5, we present models of the kilonova emission including the contribution of freeneutron heating (Secs. 5.1 and 5.2), of the non-thermal kilonova afterglow (Sec. 5.5), and discuss their application in the context of GW170817 (Secs. 5.2.2 and 5.5.3). We also discuss the role of exceptionally high opacities associated with high ionization states of lanthanides recently calculated by Banerjee et al. (2022) with respect to early kilonova and precursor emission (Sec. 5.3). We also compute high-energy gamma-ray emission from inverse-Compton and synchrotron self-Compton processes associated with the kilonova afterglow (Sec. 5.5.4). Conclusions are presented in Sec. 6. Finally, two appendices elaborate on numerical details regarding convergence, ejecta analysis, and the use of passive tracer particles.

SIMULATION DETAILS
In the following subsections, we briefly discuss the analytic foundations (Secs. 2.1 and 2.2), numerical methods (Sec. 2.3), initial data (Sec. 2.4) as well as simulation dignostics (Sec. 2.5) employed here for simulating BNS mergers. Our methods closely follow Siegel & Metzger (2018) for the GRMHD part with weak interactions and Radice et al. (2018) for modeling neutrino transport via a one-moment closure scheme ('M0 scheme').

GRMHD with weak interactions
We model BNS mergers within ideal GRMHD coupled to weak interactions. The equations of ideal GRMHD comprise energy-momentum conservation, baryon and lepton number conservation, Maxwell's equations in the ideal MHD regime, and Einstein's equations: and Here, G µν is the Einstein tensor, u µ denotes the fourvelocity of the fluid, n b , and n e are the baryon and electron number density, F * µν is the dual of the Faraday electromagnetic tensor, and the quantities Q and R represent source terms due to weak interaction and neutrino transport (neutrino absorption and emission). The energy-momentum tensor is given by where ρ = n b m b is the rest-mass density and m b is the baryon mass, p is the pressure, h = 1 + + p/ρ denotes the specific enthalpy, is the specific internal energy, b µ ≡ (4π) −1/2 F * µν u ν is the magnetic field vector in the frame comoving with the fluid, b 2 /2 ≡ b µ b µ /2 is the energy density of the magnetic field, and g µν is the spacetime metric. To close the equations, we assume a threeparameter, finite-temperature, composition-dependent EOS, i.e., we assume that every dependent thermodynamic variable P is a function of mass density, temperature, and electron fraction, P ≡ P(ρ, T, Y e ). The electron or proton fraction is defined by Y e = n p /(n p + n n ), where n p and n n denote the proton and neutron number density, respectively. In accordance with the tabulated EOS employed here (Sec. 2.4), we set the baryon mass m b to the free neutron mass for APR and LS220, and to the atomic mass unit in the case of the SFHo EOS (Sec. 2.4). Using a 3+1 decomposition of spacetime into nonintersecting space-like hypersurfaces of constant coordinate time t and timelike unit normal n µ (Lichnerowicz 1944;Arnowitt et al. 2008) , we can write Eqs. (1)-(5) in conservative form as where f (i) (p, q) and s(p) represent the fluxes and the source terms, respectively (see Eqs. (18) and (19) in Siegel & Metzger 2018), α denotes the lapse function, γ is the determinant of the spatial metric γ ij on the hypersurfaces, denotes the vector of conserved (evolved) variables, and is the vector of primitive (physical) variables. The conserved vector is composed of the conserved mass density D = ρW , the conserved momenta S i = −n µ T µ i , the conserved energy τ = n µ n ν T µν − D, the three-vector components of the magnetic field B µ = (4π) −1/2 F * µν n ν , as well as the conserved electron fraction DY e -as measured by the Eulerian observer, defined as the normal observer moving with four velocity n µ perpendicular to the spatial hypersurfaces. This observer measures a Lorentz factor of the fluid of W = −u µ n µ , which moves with three-velocity v i = u i /W on the hypersurfaces.

Neutrino leakage scheme and neutrino transport
The composition of the fluid as traced by Y e , as well as its internal energy and momentum, can change due to weak interactions (i.e., due to the emission and absorption of neutrinos). In the evolution equations, the effect of these interactions is represented by the source terms Q and R of the right-hand side of Eqs. (1) and (3), which correspond to the net energy loss/deposition and net lepton number emission/absorption rate per unit volume, respectively. We use an energy-averaged leakage scheme to treat neutrino emission at finite optical depth as described in Siegel & Metzger (2018), closely following the methods of Galeazzi et al. (2013) and Radice et al. (2016b), which are, in turn, based on Ruffert et al. (1996).
In more detail, we track the reactions involving electron neutrinos, ν e , electron anti-neutrinos,ν e , and the heavy-lepton neutrinos ν µ and ν τ summarized as ν x . In the rest frame of the fluid, the net neutrino heating/cooling rate per unit volume and the net lepton emission/absorption rate per unit volume are given as a local balance of absorption and emission of freestreaming neutrinos: Here, ν i = {ν e ,ν e , ν x } and κ νi , n νi , and E νi denote the corresponding absorption opacities, number densities, and mean energies of the free-streaming neutrinos in the rest frame of the fluid, respectively. The effective emission rates R eff νi and Q eff νi are calculated from intrinsic (free) emission rates R νi and Q νi by taking into account the finite optical depth due to the surrounding medium: (13) For a given neutrino species ν i , we define t diff,νi = D diff τ 2 νi k −1 νi as the local diffusion timescale, where τ νi is the optical depth, and D diff ≡ 6 is the diffusion normalization factor (O'Connor & Ott 2010). We define the local neutrino number and energy emission timescales as where n eq νi is the neutrino number density in chemical equilibrium, and e νi denotes the corresponding neutrino energy densities. The intrinsic emission rates take into account charged-current β−processes, electron-position pair annihilation, and plasmon decay (see Eqs. (25)-(32) in Siegel & Metzger 2018). Neutrino opacities κ νi include contributions from absorption of electron and anti-electron neutrinos and from scattering on heavy nuclei A and on free nucleons (see Eqs. (33)-(39) in Siegel & Metzger 2018 and Eq. (8) of Radice et al. 2018). We compute optical depths τ νi using the effective local approach of Neilsen et al. (2014).
Finally, using a zeroth-moment (M0) scheme to approximate the general-relativistic Boltzmann transport equation (Radice et al. 2016b), we evolve the number density n νi and average energy E νi of free-streaming neutrinos. In this scheme, only one moment of the neutrino distribution functions is used, the neutrino number currents J µ νi . Neglecting scattering, the first moment of the Boltzmann equation yields an equation for lepton number conservation, which is closed by assuming that free-streaming neutrinos propagate on null radial rays, J µ νi = n νi k µ , with k µ a null four-vector proportional to the radial coordinate of the grid, and normalized such that k µ u µ = −1. In a stationary spacetime with time-like Killing vector ∂ t and in the absence of interactions with the fluid, the mean neutrino energy E νi ≡ −p µ ∂ µ t as seen by an observer with four-velocity ∂ µ t is conserved along the neutrino worldline (parametrized by the affine parameter s), dE νi /ds = d(−p µ ∂ µ t )/ds = 0. Here, p µ = E νi k µ is the average neutrino four-momentum. In the presence of heating and cooling, this becomes where Q νi,h = −k µ ∂ µ t (Q eff νi /n νi ) is the local effective average energy deposition rate per neutrino due to neutrino absorptions and Q νi,c = E νi (R eff νi /n νi ) is the average energy loss rate per neutrino due to neutrino emission. The assumption of a stationary spacetime is approximately satisfied once the neutron stars merge and the spacetime becomes approximately stationary. We switch on this transport scheme only as the initially cold neutron stars start to merge and matter is heated up to at least several MeV to copiously produce neutrinos. The neutrino number densities n νi and mean energies E νi are obtained by evolving Eqs. (15) and (16) on radial rays (radial worldlines) using a separate spherical grid that requires interpolation at each timestep onto the Cartesian grid on which the GRMHD equations are evolved.
This combined neutrino leakage and M0 transport approach include general-relativistic and Doppler kinematic effects and it is computationally far cheaper than more complex schemes such as M1 schemes in GRMHD (Foucart 2018;Li & Siegel 2021;Radice et al. 2022) or Monte-Carlo based radiation transport schemes (Miller et al. 2019a;Foucart et al. 2020) and thus allow for longer evolution timescales and wider parameter space studies. While quantitative differences exist between these methods, it is reassuring that quantitative comparisons between the methods Foucart et al. 2020) reveal remarkable agreement for most physical observables, in particular with respect to ejecta mass, velocity, and nucleosynthesis. Deviations in individual physical quantities are typically on the ∼ 10% level (except for ν x neutrinos, which, however, are not the focus here), within the error budget of other assumptions and approximations.

Numerical set-up
We evolve Einstein's equations coupled to the generalrelativistic MHD equations in conservative form using the open-source framework of the EinsteinToolkit 2 (Goodale et al. 2003;Schnetter et al. 2004;Thornburg 2004;Loffler et al. 2012;Babiuc-Hamilton et al. 2019). Regarding the matter part, we use an enhanced version of the public GRMHD code GRHydro (Baiotti et al. 2003;Mösta et al. 2013) presented in Siegel & Metzger (2018) and . We implement a finite-volume scheme using piecewise parabolic reconstruction (Colella & Woodward 1984) and the approximate HLLE Rieman solver (Harten et al. 1987). The magnetic field is evolved using the FluxCT method (Tóth 2000) to maintain the solenoidal constraint during evolution. The recovery of primitive variables is carried out using the framework presented in , which provides support for general finite-temperature, composition-dependent, three-parameter EOS. We set a static atmosphere at a rest-mass density of 5% above the minimum density of the tabulated EOS. For most of our models, this is ρ atmo 6×10 2 g cm −3 , which is sufficiently small to capture the most tenuous mass outflows from the system. Indeed, at a typical extraction radius of r ≈ 450 km, the total mass of the atmosphere enclosed at that radius is (4/3)πr 3 ρ atmo ≈ 1 × 10 −7 M , which is an order of magnitude smaller than the total amount of material in the fast ejecta tails (v > 0.6c; see below).
2 http://einsteintoolkit.org Our neutrino leakage and M0 transport scheme is based on the implementation in WhiskyTHC (Radice & Rezzolla 2012;Radice et al. 2016b). For the M0 transport, we use a uniform spherical grid of 2048 radial rays extending up to 765 km, and a radial resolution of ∆r = 240 m. We activate neutrino heating when the center of mass of the stars is at 12 km, i.e. just after the NSs touch.
Spacetime is evolved using the BSSNOK formulation (Nakamura et al. 1987;Shibata & Nakamura 1995;Baumgarte et al. 1999) as implemented in the McLachlan code (Brown et al. 2009). We use the 1 + log gauge and the Γ-driver shift condition with the damping parameter set to η = 0.75. We use a fourth-order finitedifference scheme with fifth-order Kreiss-Oliger dissipation and dissipation parameter = 0.1. At the time of black-hole formation, we start to track the apparent horizon with the methods of the AHFinderDirect code (Thornburg 2004) and we set the hydrodynamic variables to atmosphere floor values well inside the horizon. We also excise the M0 transport grid within a radius close to the apparent horizon and we keep the origin of the spherical transport grid centered on the black hole and move it with the black hole if the latter receives a 'kick' after merger due to non-axisymmetric effects.
Time integration of the coupled Einstein and GRMHD equations is performed with the method of lines, using a third-order Runge-Kutta scheme (Gottlieb et al. 2009). The Courant-Friedrichs-Lewy (CFL) condition is set using a CFL factor of 0.25.
In all simulations, we use a fixed Cartesian grid hierarchy composed of six nested Berger-Oliger refinement boxes, doubling the resolution on each level, as provided by the Carpet infrastructure (Schnetter et al. 2004). The finest mesh grid covers a radius of 38 km during the entire evolution and initially contains both NSs. We use fixed refinement boxes and overlap zones between levels to keep violations of the magnetic solenoidal constraints under control (Mösta et al. 2014) and confined to the refinement boundaries without impacting the dynamics. The outer boundary is placed at a radius of 1528 km. For each EOS (see below), we perform two sets of simulations: (i) a set of high-resolution simulations, setting the finest spatial resolution to ∆x = 180 m and imposing reflection symmetry across the equatorial (orbital) plane for computational efficiency, and (ii) a set of mid-resolution simulations with finest spatial resolution of ∆x = 220 m without imposing any symmetries. This allows us to test how ejecta properties change with resolution and to analyze the influence of reflection symmetry with regard to the properties of the magnetic field.

Models: EOS and initial data
The properties of BNS mergers can be very sensitive to the choice of the EOS (e.g., Bauswein et al. 2013b;Radice et al. 2016bRadice et al. , 2018Hotokezaka et al. 2013b;Nedora et al. 2021b). The compactness of the star, as well as finite-temperature and composition-dependent effects, impact the lifetime of a merger remnant, the properties of matter outflows, the emitted gravitational waves, etc. We use three different nuclear EOSs in our set of simulations as tabulated by O'Connor & Ott (2010) and Schneider et al. (2017), available at stellarcollapse.org. These tables tabulate all necessary thermodynamic quantities as a function of the independent variables (ρ, T, Y e ): • LS220: We use the liquid-drop model with Skyrme interactions as implemented in Schneider et al. (2017), using a parametrization corresponding to Lattimer & Swesty (1991). At low densities, the EOS transitions from a single-nucleus approximation to nuclear statistical equilibrium (NSE) using 3335 nuclides (Schneider et al. 2017).
• SFHo: This EOS parametrizes nucleonic matter with a relativistic mean field model (Steiner et al. 2013). The EOS describes nuclei and nonuniform nuclear matter using the statistical model of Hempel & Schaffner-Bielich (2010).
• APR: This EOS is based on the potential model of Akmal et al. (1998) built to fit results from nucleon-nucleon scattering and properties of light nuclei. We use the new 3D, finite-temperature implementation presented in Schneider et al. (2019).
All three EOS transition to a nuclear statistical equilibrium including heavy nuclei in a thermodynamically consistent way, thus including the release of binding energy as alpha-particles and heavier nuclei form. This recombination energy of ∼ 7 MeV per baryon can be an important source of energy to accelerate outflows. This effect is important, in particular, in the context of post-merger accretion disks (Fernández & Metzger 2013;Siegel & Metzger 2017. The full 3D tabulated LS220 and SFHo EOSs have been used in general-relativistic hydrodynamic (GRHD) simulations of BNS mergers (Sekiguchi et al. 2015;Nedora et al. 2021b;Radice et al. 2018) and GRMHD simulations of the post-merger phase (Miller et al. 2019a;Li & Siegel 2021;Mösta et al. 2020), including weak interactions. APR-type EOSs have been used in both GRHD and GRMHD simulations but mostly using an ideal-gas approximation for the finite-temperature part (Ciolfi et al. 2017) or without weak interactions (Hammond et al. 2021). In this work, we present the first simulations of BNS mergers and their post-merger evolution with the APR EOS using the full 3D tabulated EOS and weak interactions.
The choice of EOSs roughly spans the current range of allowed neutron-star radii of ≈ 11 − 13 km (De et al. 2018;Miller et al. 2019b;Riley et al. 2019;Capano et al. 2020;Landry et al. 2020;Dietrich et al. 2020) for a 1.4M neutron star (Fig. 2), thus covering the range of different allowed compactness scenarios at given mass. The radii of a non-rotating, cold 1.4M neutron star in β-equilibrium are 11.6, km, 11.9 km, and 12.7 km for APR, SFHo, and LS220, respectively; the maximum masses of non-rotating neutron stars for these EOS are 2.2 M , 2.06 M and 2.04 M , respectively. With differences in neutron-star compactness and maximum mass, these three EOSs give rise to a range of merger and post-merger phenomena. First, the 'softer' APR and SFHo EOSs are expected to produce more violent collisions, owing to the faster relative velocities of the stars at merger (Bauswein et al. 2013b), and thus faster ejecta emanating from the collision interface than in the case of the 'stiffer' LS220 EOS. Furthermore, the APR EOS leads to long-lived remnants for a typical 1.35M − 1.35M binary (Ciolfi et al. 2019;Ciolfi & Kalinani 2020), owing to the relatively high maximum TOV mass, while, for similar parameters, the SFHo and LS220 EOSs lead to short-lived remnants with light and heavy post-merger disks, respectively (Nedora et al. 2021b;Mösta et al. 2020).
We simulate a set of equal-mass binary neutron star systems on quasi-circular orbits starting at initial separations of 40 km.
At this separation, the stars inspiral for ≈3 orbits before they merge. The relatively small separation introduces a small eccentricity to the quasi-circular orbit, affecting the GW waveforms; this, however, has a small impact on the hydrodynamics of the merger, compared to other properties, such as initial compactness of the stars (which is the focus of the present paper). Individual stars have ADM masses of M ∞ = 1.35M (at infinite separation), typical of galactic merging double neutron-star systems (Özel & Freire 2016). Our selection of binary models are compatible with the inferred source parameters of GW170817, with chirp mass of M c = 1.188 +0.004 −0.002 M and total gravitational mass of 2.74 +0.04 −0.01 (low-spin prior) and mass ratio of 0.7-1.0 at 90% confidence (Abbott et al. 2017a). Moreover, our simulations do not result in prompt collapse upon merger, which is strongly disfavored in the case of GW170817 by electromagnetic observations of J0348+0432 J0740+6620 APR SFHo LS220 GW+PSRs+X-rays Figure 2. Gravitational mass as a function of stellar radius of non-rotating neutron stars in β-equilibrium at zero temperature for different EOSs. Dots mark models with 1.4 M and corresponding radii 11.6,km, 11.9 km, and 12.7 km for APR, SFHo, and LS220, respectively. The grey shaded area represents mass-radius constraints (90% credible intervals) from current non-parametric inferences using heavy pulsar measurements, gravitational-wave observations of GW170817 and GW190425, and NICER data on J0030+0451 (Landry et al. 2020). The blue shaded area indicates the 68% credible mass estimates from the heaviest pulsars measured so far (Cromartie et al. 2020;Demorest et al. 2010) the GW170817 kilonova (e.g., Margalit & Metzger 2017;Bauswein et al. 2017;Siegel 2019). The initial data is built using the open-source code LORENE (Gourgoulhon et al. 2001) 3 assuming that the stars are at zero temperature, β-equilibrated, and nonrotating. For this construction, we slice the 3D EOS table assuming the aforementioned conditions, subtracting the pressure contribution of photons (Radice et al. 2016b), and generate a high-resolution table as a function of rest-mass density (i.e. an effectively barotropic table). The initial data is imported into the evolution code adapting the methods in WhiskyTHC 4 .
After setting the hydrodynamical variables, we initialize the magnetic field as a (dynamically) weak poloidal seed field buried inside the stars. For this purpose, we specify the vector potential at each star according to where r cyl is the cylindrical radius relative to the star's center, ρ max is the initial maximum rest-mass density, and {A b , n p , n s , p cut } are free parameters (Liu et al. 3 We use the secant fix suggested in https://ccrgpages.rit.edu/ ∼ jfaber/BNSID/ 4 https://bitbucket.org/FreeTHC 2008). We choose n p = 0 to set the maximum strength of the field at the center of the star, n s = 3, and p cut is set to p cut = 0.4p max in order to keep the field within the stellar interior and well off the stellar surface. Here, p max is the initial maximum pressure within the star. Finally, we set A b such that the initial maximum field strength is B max = 5 × 10 15 G for all BNS models. This high value of the magnetic field strength is unlikely present in typical BNS, in which surface values of the poloidal component are expected to be closer to B pole ∼ 10 12 G, as observed in radio pulsars (Tauris et al. 2017), while in our case, the strength of the magnetic field near the surface (where p = 0.4p max ) is ≈ 3 × 10 14 G.
Our initial magnetic fields anticipate field strengths similar to those expected promptly after merger without the need to resolve the computationally extremely costly amplification process. Magnetic fields are amplified during the merger process to an equipartition level of ∼ > 10 15 − 10 16 G via the Kevin-Helmholtz instability (e.g., Price & Rosswog 2006;Kiuchi et al. 2014;Anderson et al. 2008;Zrake & MacFadyen 2013 and the magnetorotational instability (e.g., Balbus & Hawley 1991;Duez et al. 2006b;Siegel et al. 2013) on ms timescales. Specifying an initial poloidal seed field is an acceptable assumption motivated by the fact that the topology is 'reset' during merger due to global dynamical effects (see below) as well as by the aforementioned small-scale amplification effects (Palenzuela et al. 2021;Aguilera-Miret et al. 2020). Although large at absolute value, the initial magnetic fields are sufficiently weak to be dynamically insignificant and they do not alter the inspiral dynamics of the stars.
Finally, our initial magnetic field configuration neglects the external magnetosphere and large-scale fields likely present in the inspiral phase of a BNS system. Simulating regions of extremely high magnetic pressure in ideal MHD is challenging; various strategies have been pursued, such as imposing high-density atmospheres (Paschalidis et al. 2015;Ruiz et al. 2016) to reduce the plasma β parameter and to maintain stability, improving the conservative-to-primitive scheme , or using non-ideal (resistitive) MHD approximations (Andersson et al. 2022). Our conservativeto-primitive scheme has already been significantly optimized to handle large plasma-β. However, in order to handle magnetospheric field strengths of ∼ 10 14 G using the first approach, one would still need to increase the atmospheric density to levels at which it absorbs much of the high-velocity ejecta. This would thus preclude a detailed analysis of dynamical ejecta, which is one of the foci of the present paper. While the large-scale structure of the pre-merger field could have an important role in launching a short GRB after merger (Mösta et al. 2020;Ruiz et al. 2016), the properties of the magnetosphere in BNS are still largely unknown. For the purposes of analyzing the dynamical ejecta, we expect magnetospheric effects to be at least of second-order importance, compared to other hydrodynamical effects.

Simulation diagnostics and analysis
In the following, we briefly describe some key simulation diagnostic tools.

Ejecta properties
This work focuses on dynamical ejecta from the merger process itself and the electromagnetic radiation it can give rise to. There are various proposed criteria to determine whether a fluid element becomes unbound from the merging system. Assuming a fluid element moves along a geodesic in a stationary, asymptotically flat spacetime, neglecting fluid pressure and self-gravity, a fluid element can reach spatial infinity if the specific kinetic energy at spatial infinity is E ∞ := −u t + 1 > 0, where −u t is conserved along the geodesic. The fluid then has an asymptotic escape Lorentz factor given by and the escape velocity is v ∞ = 1 − 1/W 2 ∞ . This is known as the geodesic criterion. Because conversion of thermal energy into kinetic energy by pressure gradients can occur in outflows, one guaranteed source of internal energy being the recombination of nucleons into alpha-particles and heavier nuclei, this criterion imposes a lower limit on the total unbound mass.
Assuming a stationary fluid, a stationary spacetime, and an asymptotic value of the specific enthalpy of h ∞ , an unbound fluid element fulfills the Bernoulli criterion, if the specific energy at infinity satisfies −hu t > h ∞ . Typically, the asymptotic enthalpy is h ∞ = 1, if the EOS is independent of composition (e.g., for a polytropic EOS or ideal gas). In general, however, the asymptotic enthalpy depends on Y e . The asymptotic value of Y e varies from fluid element to fluid element if one takes into account neutrino interactions and r-process nucleosynthesis. The former necessarily plays an important role for most ejecta elements due to the high (∼ 10 MeV) temperatures reached during the merger process, and the latter necessarily sets in as the ejected fluid element decompresses and/or moves out of nuclear statistical equilibrium (NSE) as it is being ejected from the merger site. We calculate h ∞ taking into account the contribution of r-process heating as h ∞ = 1 + ∞ , where ∞ is the average binding energy of the nuclei formed by the r-process (see also Foucart et al. 2021;Fujibayashi et al. 2020). This number depends on how the binding energy is defined in the EOS table, and it is approximately independent of Y e if we ignore neutrino cooling during the r-process. For our SFHo table, where the reference mass is the atomic mass unit, we have h ∞ 1, while LS220 and APR use the free neutron mass, which corresponds to h ∞ 0.992. Although this difference is small, it has a non-negligible effect on the measured kinetic energy of the ejecta. Using this criterion, the asymptotic Lorentz factor is: We note that the Bernoulli condition might overestimate the total unbound mass if the system is not stationary (Kastaun & Galeazzi 2015). We analyze the physical properties of the ejecta using two methods: discretizing the outflow on spherical detector surfaces and computing surface integrals, and by injecting passive tracer particles into the simulation domain that record plasma properties along their Lagrangian trajectories. We place different spherical surface detectors at radii R D ≈ {300, 450, 600, 750, 900, 1050, 1350} km, and extract various quantities of the fluid on spherical surface grids with resolution (n θ , n φ ) = (56, 96). The cumulative ejected mass across a detector (coordinate) sphere with radius R D is given by where dA = R 2 D sin θ dθdφ, and W(θ, φ) is a function that is one (zero) if the fluid at a given grid point is unbound (bound). We also calculate mass-averaged properties of the fluid (such as Y e and specific entropy) through the detectors as where dM i,j is the amount of mass passing through the detector during a time interval dt i for the fluid component with property P ∈ [P(j) − δ, P(j) + δ], where δ refers to an appropriately chosen bin width; M i is the total mass crossing the detector surface in dt i , and T is the total time interval considered. When using surface integration to calculate outflow properties, one ideally uses large surface radii for extraction to capture the entire unbound outflow component and to ensure that the fluid is approximately stationary. This requires, in particular, to evolve the system for sufficiently long timescales for outflows to cross the relevant detector spheres. Here, we evolve the systems for more than 30 ms post-merger, which is sufficient to yield approximate convergence of the unbound component at large radii using the geodesic criterion (see Appendix A).

Tracers
Ejected material from BNS mergers undergoes rapid neutron-capture nucleosynthesis and produces heavy elements. We compute nucleosynthesis abundances arising from the simulation ejecta using a nuclear reaction network. Specifically, we sample thermodynamic properties such as ρ, T, Y e , and W using several families of passive tracer particles injected into the simulation domain and input the recorded Lagrangian fluid trajectory histories into the nuclear reaction network in a postprocessing step (Sec. 2.5.3).
We place a total of ∼ 10 5 tracer particles into the crust of each star using a probability function proportional to rest-mass density. At the time of merger, but before ejection of material commences, we allocate 1/4 of all tracers into a second family of tracers and resample regions of high specific entropy (s > 30k b ) using the same probability function, to ensure sufficient sampling of the shock-heated ejecta component that originates in the collision interface. Finally, after merger, when the system has settled into a quasi-stationary state (∼ 20 ms after merger), we allocate 1/4 of all tracers into a third family to adequately track the outflows from the post-merger accretion disks. These tracers are placed within a radius of r ∈ (20, 100) km and a polar angle of θ ∈ (40, 140) using the same probability function as before. This three-stage resampling method allows us to accurately identify, track, characterize, and distinguish between different components of outflows, which are ejected by different physical mechanisms.
In order to compute mass distributions and massweighted averages of ejecta properties based on tracers, we (re)assign mass to each tracer once it becomes part of an outflow. For a set of tracers crossing a sphere of fixed radius (R D = 440 km) at a given time t, we define the mass of those tracers integrating the mass-flux using the rest-mass density and velocity recorded by the where ∆t is the output frequency of the tracer data (see also Bovard & Rezzolla 2017). We find that this method of assigning mass very well reproduces the mass distribution of outflows as extracted by the grid-based spherical surface outflow detectors. In Appendix B, we present a comparison between the Y e and v ∞ distributions as obtained by the grid-based and tracer methods, respectively.

Nucleosynthesis setup
Using the full set of unbound tracers, we perform nuclear reaction network calculations with the open-source reaction network SkyNet 5 (Lippuner & Roberts 2017) in post-processing. SkyNet includes 7852 isotopes from free neutrons to 336 Cn, including 140000 nuclear reactions. Nuclear masses (experimental data where available, data from the finite range droplet macroscopic model (FRDM; Möller et al. 2016) otherwise), partition functions, and forward strong reaction rates are taken from the JINA REACLIB database (Cyburt et al. 2010). Inverse strong reaction rates are derived from detailed balance. Where available, weak reaction rates are taken from Fuller et al. (1982), Oda et al. (1994), Langanke & Martınez-Pinedo (2000), and from REACLIB otherwise. Spontaneous and neutron-induced fission rates are obtained from Frankel & Metropolis (1947), Mamdouh et al. (2001), Panov et al. (2010), and Wahl (2002).
The Lagrangian trajectories of unbound tracers are followed through the grid until the end of each simulation or until they exit the grid. At that point, we smoothly transition and extrapolate their trajectories assuming homologous expansion (ρ ∝ t −3 ). We start evolving the composition with SkyNet at t = t 7GK , once the temperature drops below 7 GK. At this point, a fluid element represented by a tracer particle is still in NSE, and SkyNet initially evolves the composition in NSE accordingly until the temperature drops further to 5 GK. At t = t 5GK NSE starts to break down and SkyNet automatically switches over to a full network evolution. As a result the Y e of a fluid element then decouples from the Y e evolution on the simulation grid. Reaction network calculations are performed up to 10 9 s, when a 'stable' final abundance pattern has emerged. Using the full distribution of tracer particles instead of, e.g., feeding the outflow history recorded by spherical detectors into the nuclear reaction network (e.g., Radice et al. 2018;Nedora et al. 2021b) allows us to distinguish between different ejecta components and to determine the initial conditions for nucleosynthesis self-consistently for each ejecta component and fluid element of the flow.

Overview of simulations
In each BNS simulation, the NSs inspiral for approximately three orbits, emitting GWs, before merging into a massive remnant neutron star. All remnants created in these runs are gravitationally unstable and eventually collapse into a black hole upon removal or redistri- . Key dynamical ejecta properties as measured at a radius of 440 km according to the geodesic criterion: histograms of estimated asymptotic ejecta velocity (v∞), electron fraction (Ye), and specific entropy (s), for each high-resolution simulation. We choose the geodesic criterion here to largely exclude secular ejecta from the remnant NS (not of interest here) and to focus on dynamical ejecta only. The high-velocity tails of the ejecta distributions that give rise to free-neutron decay and associated kilonova precursor emission (Secs. 4.2 and 5.2) are indicated as color-shaded areas. Table 1. Properties of our suite of high-resolution BNS simulations. Here, M1,2 denotes the gravitational mass of the neutron stars, RNS their radius, d12 is the initial binary separation, and Bmax(0) is the initial maximum magnetic field strength; tBH is the time of black-hole formation, Mej is the total dynamical ejecta mass, Mej(v > 0.6) is the amount of fast moving dynamical ejecta, and Mej,n is the amount of ejecta mass in free neutrons. In brackets, we show several mass and time-averaged ejecta properties, such as the electron fraction (Ye), the specific entropy (s), the asymptotic velocity (v∞), the angular direction of the ejecta measured from the pole (θ), and the electron fraction extracted at T = 5 GK by tracer particles (Ye(T = 5 GK)). bution of sufficient angular momentum due to magnetic stresses, spanning a wide range of lifetimes (cf. Tab. 1).
The remnant in the APR case is a long-lived star with mass in the supramassive 6 regime and a lifetime of likely more than a few hundred milliseconds (Ciolfi et al. 2019). In contrast, the LS220 and SFHo binaries considered here lead to stars in the hypermassive 7 regime with short lifetimes of ≈ 15 ms and ≈ 21 ms, respectively. Our simulations self-consistently incorporate weak interactions and approximate neutrino transport, which is pivotal for accurately modeling ejecta properties, the compositional distribution represented by Y e , in particular. Furthermore, our simulations include magnetic fields, which play a key role in generating angular moment transport and outflows in the post-merger phase. In our setup, magnetic fields are initialized well inside the stars (cf. Sec. 2.4) and only 'leak' out of the stars during the inspiral in an insignificant way. At merger, β −1 ≡ b 2 /p remains small, and the 'buried' fields do not influence the ejection of dynamical ejecta. In this early stage of the merger process, our results resemble closely purely hydrodynamic simulations that include weak interactions and approximate neutrino absorption, but neglect magnetic fields (e.g., Sekiguchi et al. 2016;Radice et al. 2018).
In this paper, we focus on ejection mechanisms, ejecta properties, and observables of material ejected during the dynamical phase of the merger itself. We consider dynamical ejecta only, defined as material that is unbound by global dynamical processes. Table 1 provides an overview of the mass-averaged properties of the dynamical ejecta. Corresponding distributions of ejecta mass relevant for observables according to composition (Y e ), asymptotic escape speed (v ∞ ), and specific entropy (s ∞ ) are summarized in Fig. 3. We extract physical quantities at a radius of R = 300 M 440 km, where M is the total binary mass, using the geodesic criterion. We mainly focus on the geodesic criterion here, since at close separations of 440 km it is somewhat insensitive to secular outflows such as neutrino-driven winds from the merger remnant (not of interest for the present study) and it thus acts as a filter for dynamical ejecta. The merger process during which dynamical ejecta is generated according to the geodesic criterion lasts approximately 10 ms in all our simulations (Sec. 3.2 and Fig. 6). We turn to a discussion of the details of mass ejection in the following subsections.

Ejecta dynamics and fast outflow
Two types of ejecta can be distinguished at merger: tidal and shock-heated ejecta. Tidal torques extract material from the surface of the stars during the final inspiral and merger process, creating spiral arms that expand into the orbital plane as they transport angular momentum outwards, expelling cold, neutron-rich material (Y e ∼ < 0.1) into the interstellar medium (see Fig. 4, first panel). Because neutron stars are more compact in general relativity compared to Newtonian gravity, these tidal tails are not as prominent here as in Newtonian simulations (e.g., Rosswog et al. 1999;Korobkin et al. 2012;Rosswog 2013). Furthermore, for equal-mass binaries one expects a minimum of tidal ejecta: For a given EOS, tuning the binary mass ratio away from unity generally enhances the tidal torque on the lighter companion. This leads to increased tidal ejecta while reducing the shock-heated component originating in the collision interface (Hotokezaka et al. 2013a;Bauswein et al. 2013b;Dietrich et al. 2015;Lehner et al. 2016a;Sekiguchi et al. 2016). This is because the less massive companion becomes tidally elongated and seeks to 'avoid' a (radial) collision. Finally, for a given binary mass ratio, changing the EOS from stiff (large NS radii) to soft (small NS radii) one expects the shock-heated component to be enhanced while reducing the tidal component (Hotokezaka et al. 2013a;Bauswein et al. 2013b;Dietrich et al. 2015;Lehner et al. 2016a;Sekiguchi et al. 2016;Palenzuela et al. 2015). This is because tidal forces are smaller for less extended objects, and NSs with smaller radii approach closer before merger, reaching higher orbital velocities at the collision, thus enhancing the shock power and associated ejecta mass.
With our NSs spanning the compactness range of currently allowed EOSs for typical galactic double neutron star masses, we find our runs span a range of dynamical mass ejection phenomena. A detailed analysis shows (see below) that for all systems considered here, by far most of the ejecta is expelled by shock waves produced in quasi-radial bounces of an oscillating double-core remnant structure that forms after the onset of the merger, with only a negligible amount of material being ejected by tidal tails (see Fig. 5 and below; Sec. 3.3). This ejecta material is, in general, faster and more protonrich than tidal ejecta. A large fraction of the material released in such waves has been heated considerably due to hydrodynamical shocks at the collision interface during the merger process and is further heated as it shocks into slower surrounding merger debris. Associated neutrino emission in such a neutron-rich environment favors electron captures onto neutrons over positron captures on protons, raising the electron fraction and causing the electron fraction to widen from the initially cold and neutron-rich conditions of the individual NSs (Y e ∼ < 0.1) into broad distributions with prominent high-Y e tails (Fig. 3). The stronger the heating, as indicated by the specific entropy distributions (Fig. 3), the higher the resulting final Y e of the ejecta. Neutrino reabsorption, which also contributes to increasing Y e in a neutron-rich environment, only plays a minor role in setting Y e at this early dynamical stage. Overall, the more compact NS binaries (APR and SFHo) eject a factor of 2 − 3 more ejecta than the less compact LS220 binary (cf. Tab. 1), since tidal ejecta is subdominant and post-merger bounces are more energetic and lead to more violent shocks unbinding more material for more compact NSs. Figure 6 summarizes a more detailed analysis of the matter ejection mechanism post-merger, correlating the generation of unbound material with radial oscillations of the remnant. Unbound mass, as sampled by passive tracers according to the Bernoulli criterion, initially increases in 'steps' when the maximum density reaches a minimum during oscillations.
At this point of maximum decompression and least compactness, some amount of material becomes unbound and is ejected, the remnant undergoes a bounce, and the maximum density starts to increase again. In the softest EOS simulations, nearly 75% of the outflow becomes unbound after the first three bounces, followed by a secular growth that lasts up to 10 ms after  merger, reflecting the integrated action of several weaker bounces. In the stiffer LS220 case, only two individual mass ejection episodes can be discerned by this tracer method ( 30% of the total ejected mass). The Bernoulli criterion used here is valid for stationary fluid flows only and might overestimate the total ejecta mass near shocks (see Kastaun & Galeazzi 2015 for a discussion). However, confidence in the present interpretation comes from the fact that all ejecta curves are monotonically increasing over time, together with the fact that the amount of unbound material for both the geodesic and the Bernoulli criteria at large distances agree up to saturation of the dynamical ejecta as measured by the geodesic criterion (cf. Fig. 22). The importance of post-merger bounces for mass ejection has been noted previously in smoothed-particle hydrodynamics (SPH) simulations (Bauswein et al. 2013b) and more recently in grid-based GRHD simulations (Radice et al. 2018;Nedora et al. 2021a).
The three (two) individual mass ejection episodes discussed above for APR and SFHo (LS220) drive shock waves into the surrounding medium (see Fig. 4) and manifest themselves in step-wise increases of high specific entropy mass flux (s > 20k B baryon −1 ) at large distances (see Fig. 5). Figure 5 also indicates that the first high-entropy mass flux events for LS220 and APR are associated with predominantly equatorial outflows (within 45 • from the plane), while the first event for SFHo also contains a significant fraction of high-entropy material in the polar direction (within 45 • from the polar axis). In general, mass ejection in equatorial directions dominates over po-lar directions, as shown by Figs. 5 and 8. However, the relative ratio of equatorial to polar ejecta is larger for the softer EOSs (SFHo and APR), owing to stronger bounces and associated shock waves, cf. Figs. 5 and 8). Figure 5 also illustrates that the first two (one) mass ejection episodes discussed above are (is) associated with the ejection of fast material with asymptotic speed v ∞ > 0.6c. We find a fast-moving ejecta component with v ∞ > 0.6c of mass ≈ (4 − 5) × 10 −6 M for the APR and SFHo runs (cf. Tab. 1), which is more than an order of magnitude larger than the corresponding sweptup atmosphere mass and is thus well captured; for the stiff LS220 case, we also find a fast component of mass ≈ 6×10 −7 M , which is still larger than the swept-up atmosphere mass by a factor ∼ > 7. Measuring the amount of fast ejecta at a closer distance of R D = 300 km, we find that the total mass of this fast tail is larger by about 5% for both APR and SFHo and by approximately 30% for LS220. The extracted values for the maximum velocity remain essentially insensitive to the detector radius (see also Fig. 22 for a convergence test).
This fast outflow component is ejected by the first two bounces after merger, consistent with other grid-based GRHD simulations (Nedora et al. 2021a). For SFHo, we observe in Fig. 5 that these fast outflows are ejected in two impulsive trains of similar strength, while for the even softer APR EOS, the first ejection event is much stronger. For the stiffer LS220 EOS, there is only a single fast ejection event because the bounces are generally weaker.
The angular distribution of cumulative kinetic energy E K = (Γ − 1)M ej of the ejecta as a function of βΓ We distinguish between unbound ejecta with high specific entropy (s > 20 kB baryon −1 ), fast velocities (asymptotic speed v∞ > 0.6c), as well as between polar (within 45 • from the polar axis) and equatorial (within 45 • from the equatorial plane) ejecta. We choose the geodesic criterion here to largely exclude secular wind ejecta from the remnant NS (not of interest here) and focus on dynamical ejecta only. The total amount of dynamical ejecta converges around 14 ms at this distance (see also Fig. 22) and the mass flux drops.
is shown in Fig. 7. Here, β denotes the ejecta speed in units of the speed of light and Γ the corresponding Lorentz factor. We observe that the most energetic outflows for SFHo are ejected at latitudes close to the pole (polar angles ≈ 0 • − 35 • ), consistent with previous findings (Radice et al. 2018). For the softest EOS APR, the energetics of the ejecta are dominated by near-equatorial outflows (polar angles ≈ 60 • − 90 • ). For the stiffer EOS LS220, we obtain milder energetics (mostly Γβ < 1), and a fairly homogenous angular distribution, with more energetic outflows near the equatorial plane. The total amount of fast ejecta measured in simulations is sensitive to different factors such as (a) grid resolution or sampling of high-velocity ejecta tails with a small number of particles in SPH simulations, (b) the grid topology, (c) atmosphere levels and magnetosphere treatment, and (d) definition of asymptotic velocity. In our present study (a) is fixed by constraints arising from computational cost. High-resolution Newtonian convergence studies (Dean et al. 2021) suggest that at our resolution of 180 m convergence of the high-velocity ejecta (v > 0.6c) may, in principle, be achieved within measurement uncertainties of a factor of ∼ 2. Aspect (b) is fixed by the Cartesian grid paradigm our code is based on. We elaborate on aspects (c) and (d) below.
In this work, we use lower floor/atmosphere levels than previously reported BNS simulations (e.g., smaller by a factor of 10 compared to Radice et al. 2018) to minimize the effects of an artificial atmosphere on the properties of the fast ejecta (at relevant extraction radii of ≈ 400 km). 8 The asymptotic velocity v ∞ of unbound outflows are computed using the Lorentz factor at infinity W ∞ = −u t , where W ∞ = (1 − v 2 ∞ ) −1/2 . Using the Bernoulli criterion instead, according to which W ∞ = −hu t , we find a total amount of fast ejecta 15% larger, which is expected as thermal effects can still accelerate the ejecta at this detector distance. There are other techniques used in the literature to extract outflows, e.g. volume integrations Ciolfi et al. 2017;Hotokezaka et al. 2013a), and other definitions of asymptotic velocity that explicitly include the gravitational potential (Shibata & Hotokezaka 2019). Nedora et al. (2021a,b) compute asymptotic velocities using the Newtonian expression v N,∞ = √ 2E ∞ , where E ∞ is the specific energy at infinity. We find this expression overestimates the speed and amount of fast ejecta: using the latter expression, we find a larger amount (by almost an order of magnitude) of fast ejecta in our simulations, comparable to that of Nedora et al. (2021b).
Previous analyses of fast-ejecta with GRHD simulations using neutrino transport (Radice et al. 2018;Nedora et al. 2021b) and Newtonian simulations at ultra-high resolution in axisymmetry (Dean et al. 2021) report fast-ejecta masses within factors of 5 − 10 higher than found here. SPH simulations find even higher masses than grid-based simulations, of the order of 10 −4 − 10 −5 M (Metzger et al. 2015;Kullmann et al. 2022;Rosswog et al. 2022). Our default mass estimates are similar to the ones obtained in the grid-based ultrahigh resolution simulations (purely hydrodynamic, without weak interactions) of Kiuchi et al. (2017), specifically, their HB EOS case, which leads to a similar NS compactness to APR (see Figure 7 in Hotokezaka et al. (2018)). Differences to other grid-based merger simulations (Nedora et al. 2021b) can be understood in terms of the definition of asymptotic velocity (see above), while the higher level of ejecta in grid-based convergence studies (Dean et al. 2021) may be attributable to differences in the overall setup (such as, e.g., Newtonian vs. general-relativistic dynamics and geometric effects due to assumed symmetries), since ejecta masses are expected to be roughly converged according to these studies at our present resolution (at least up to a factor of a few). Fast ejecta tails are resolved by ∼ < tens of SPH particles in present SPH simulations Rosswog et al. 2022), and details about the dependence of fast ejecta masses with resolution (number of particles) have been reported only recently (Rosswog et al. 2022). Rosswog et al. (2022) show that the mass of fast ejecta decreases by a factor of ≈ 10 when the resolution is increased from 10 6 to 5 × 10 6 SPH particles used in the simulation. The total fast ejecta mass in the latter paper-estimated using the relativistic asymptotic velocity similar to our work-is still a factor of ≈ 2 higher than the value we find for a BNS merger with a similar set-up and using the APR3 EOS, but much closer than previous SPH simulations.
Introducing a dedicated tracer family to track the shock-heated fast ejecta (cf. Sec. 2.5.2) allows us to sample the fast ejecta tail with velocities v > 0.6c by typically ∼ > 500 tracer particles. Appendix B provides more details on tracer sampling of the fast tail of the ejecta distribution in our simulations.

Electron fraction distribution and shock reprocessing
The Y e distribution of the ejecta extracted at 440 km spans a wide range of values from 0.05 to 0.5 and is centered around ≈ 0.3 for all simulations, with a more pronounced tail to higher Y e in the case of the softer APR and SFHo EOS (Fig. 3; Tab. 1). The angular dependence of the mass-weighted Y e is shown in Fig. 8. Near the pole, we observe high values > 0.3 of the electron fraction for all EOSs. This can be attributed mostly to the domi-  Figure 8. Sky map (Moillede) projection of the massweighted electron fraction (north) and of the ejected mass (south) for the total time-integrated dynamical BNS ejecta measured at 440 km and the three different EOS considered here. Equatorial outflows with low-to-medium electron fraction dominate the ejecta. Shock-reprocessing of tidal ejecta leads to azimuthal asymmetries in Ye for soft EOS ('lanthanide/actinide pockets' ).
nance of shock-heated ejecta in polar regions as well as to neutrino irradiation of the polar ejected material by the newly formed remnant. For APR and SFHo, the NS remnant reaches higher temperatures due to larger compressions and oscillations after collision (owing to higher compactness of the individual stars), whereas the remnant in LS220 is cooler. As a result, since the heating rate of material due to neutrino irradiation is roughly the same in all remnants, the entropy of the polar outflows in LS220 is slightly higher.
Near the equatorial plane, shocks reprocess part of the neutron-rich matter to values of Y e ≥ 0.1. For the soft EOSs simulations, we observe a clear azimuthal asymmetry in the Y e angular distribution, while in LS220 case the azimuthal distribution is nearly homogeneous (Fig. 8). This azimuthal asymmetry for the APR and SFHo runs occurs because the first violent bounces of the remnant produce shock waves in a preferential direction dictated by the phase of orbital rotation.
When the two NSs merge, a rotating remnant with a bar-like (m = 2 mode) shape is formed, which undergoes quasi-radial oscillations of its two-core structure (cf. Sec. 3.2, Fig. 6). The first quasi-radial bounce generates a fast shock-wave expelling material that propagates outwards without encountering significant resistance by the ambient medium. The system then develops two spiral arms due to enhanced tidal torques at maximum decompression, polluting the environment with cold neutron-rich material (see Figure 4, first panel). When the remnant bounces again, it generates another violent shock wave directed along the axis formed by the double-core structure (the x-axis in Figure 4, second panel). This fast shock wave reprocesses the previously ejected tidal tail material to high Y e through shock-heating and associated neutrino emission in a relatively wide azimuthal wedge. The resulting angular distribution shows pockets of low-Y e material around the equatorial plane (in directions orthogonal to the propagation of the second shock wave; Fig. 8). These 'lanthanide/actinide pockets' possess much higher opacities than the surrounding shock-reprocessed material (owing to the lower electron fraction, which is conducive to forming high-opacity lanthanide and actinide elements). Observable imprints of these pockets in the kilonova emission can be explored with multidimensional kilonova radiation transport simulations, which are, however, beyond the scope of the present paper.
Since bounces are much less violent in the LS220 run, azimuthal inhomogeneities in Y e as described above are much less prominent. This is because a series of weaker shock waves run into previously (tidally) ejected material, which leads to an overall more homogeneous shock heating of the ejecta across all azimuthal directions and thus to a more homogeneous azimuthal Y e profile. After the first bounces, when most of the dynamical ejecta is expelled, neutrino-driven outflows and tidal torques help launch additional material from the remnant, some of which join previously launched, circularizing debris material in forming a disk with an initial Y e ≈ 0.15 − 0.25 and specific entropy 8-10 k B per baryon (see cf. Fig. 4, third panel). These entropies that arise from a complicated process of shock reprocessing as described above are almost identical to the initial entropies of 8 k B per baryon assumed in previous post-merger disk simulations (e.g., Fernández & Metzger 2013;Siegel & Metzger 2017 Figure 9. Snapshots of the electron fraction (upper part) and magnetic field (lower part) in the meridional plane at different times after merger for the APR EOS run. White contours indicate rest-mass densities of ρ = {10 9 , 10 10 , 10 10.5 , 10 11 , 10 14 } g cm −3 .
In this section, we provide a brief overview of the postmerger evolution for each simulation. We defer a more detailed discussion of the post-merger stage to a forthcoming paper.
After merger, the deformed remnant keeps shedding material to the environment through spiral arms that contribute to forming a disk. The remnant NS cores of the LS220 and SFHo binaries keep oscillating until they merge into a single core and eventually collapse to a BH after 15 ms and 25 ms post-merger, respectively. In contrast, oscillations in the APR case damp on a much shorter (few ms) timescale (Fig. 6). The time of collapse of the remnant to a black hole depends critically on the angular momentum transport mechanisms present in the system. In this set of simulations, we incorporate magnetic fields and neutrino transport self-consistently, which play an important role in the evolution of the angular momentum of the remnant and its accretion flow.
Angular momentum in the disk is transported outwards by spiral waves, powered by the m = 2 azimuthal mode of the remnant, which then develops an m = 1 mode through hydrodynamical instabilities (see also Lehner et al. 2016b;Radice et al. 2016a). The disk increases in mass until ∼ 20 ms post-merger. At this point, the average magnetic field has increased by two orders of magnitude from its initial value and MRIdriven turbulence starts contributing to angular momentum transport in the disk, opposing outward angular momentum transport by the hydrodynamical spiral waves in the disk. As the MRI fully develops, the bulk of the disk becomes more turbulent, as can be seen qualitatively in Fig. 9 (third panel). As a result of the MRI-mediated angular-momentum transport, the accretion rate increases well above the ignition threshold for neutrino cooling in the disk (Chen & Beloborodov 2007;Metzger et al. 2008a,b;De & Siegel 2021). Thus cool-ing becomes energetically significant, the disk height decreases, and the disk neutronizes due to electrons becoming degenerate, keeping the disk midplane at mild electron degeneracy and Y e ∼ 0.1 − 0.15 in a self-regulated process (Chen & Beloborodov 2007;Siegel & Metzger 2017. This can be seen in Fig. 9 (second and third panel), where the disk Y e evolves from ∼ 0.25 to ∼ 0.15, the magnetic field strength increases throughout the disk, and the disk becomes less 'puffy'.
Material from the disk winds and merger debris initially pollutes the polar regions right after merger, suppressing outflows driven by neutrino absorption (Perego et al. 2014;Desai et al. 2022). After ∼ 15 − 20 ms postmerger, a neutrino-driven wind emerges, ejecting only mildly neutron-rich material (Fig. 9, second and third panel). Such a wind does not develop in the LS220 case, in which the remnant collapses to a black hole too early for a wind to develop. In the SFHo simulation, a coherent neutrino-driven wind is launched, but the remnant collapses shortly thereafter. The long-lived remnant formed in the APR case does develop a steady neutrinodriven wind, which is additionally aided by magnetic fields. We defer a more detailed discussion to a forthcoming paper.

r-process abundance pattern
Detailed nucleosynthesis analyses of the ejecta using a nuclear reaction network are performed in a postprocessing step (Sec. 2.5.3) based on passive tracer particles that record various thermodynamic quantities of the flow along their respective Lagrangian trajectories (Sec. 2.5.2). In each simulation, > 10% of the initially placed tracers are unbound during the dynamical phase of the merger, which constitutes approximately 5 × 10 3 − 8 × 10 3 tracers in total.
We find that the injected tracer particles sample the outflow properties very well compared to those of the mass flux recorded by detector spheres on the Cartesian grid of the Eulerian observer. Figure 24 shows the ejecta mass distribution in Y e for the SFHo run as recorded by a spherical outflow detector and as recorded by tracers passing through the same detector, as well as the total number of tracers per Y e bin. The mass of each tracer is assigned as described in Sec. 2.5.2. Similar good agreement is obtained in all other simulations and across other quantities, such as entropy, velocity, etc.; this is of pivotal importance to ensure accurate nucleosynthesis analyses. Merely the mass distribution in velocity at high velocities (v ∼ > 0.6) is typically slightly oversampled (by factors of 2-3) with respect to grid detectors (see Appendix B). This is a result of placing a large number of tracers in initially highly tenuous plasma (Sec. 2.5.2). Since only a relatively tiny fraction of mass resides at these high velocities, such oversampling does not influence the final abundance patterns (cf. Appendix B); it is beneficial for resolving the properties of and nucleosynthetic processes within the fast outflows.
The final mass-averaged abundances for the dynamical ejecta of each simulation are shown in the top panel of Fig. 10, with solar r-process abundances added for comparison as black dots. The second (mass number A ∼ 130) and third (A ∼ 195) r-process peak are wellreproduced in all simulations-a robust 2nd-to-3rd peak r-process is obtained irrespective of the EOS. This is because, as far as dynamical ejecta is concerned, weak interactions involving both emission and absorption of neutrinos (chiefly via e + +n ↔ p+ν e and e − +p ↔ n+ν e ) only have a finite ( ∼ < few ms; Fig. 6) amount of time to increase the electron fraction of merger material from the cold, highly neutron-rich conditions Y e ∼ 0.05 − 0.15 of the colliding NSs via direct shock heating at the collision interface, via reheating of neutron-rich tidal material by shock waves from the oscillating merger remnant (Sec. 3.2), and via absorption of strong neutrino radiation from the hot remnant NS that is being formed. Provided a dominant fraction of the ejecta remains at Y e ∼ < 0.25 around 5 GK when NSE breaks down, as satisfied here (cf. Fig. 10, bottom panel), a sufficiently high neutron-to-seed ratio can be achieved, such that a pile-up of material in the fission region at freeze-out occurs. The subsequent decay via fission then guarantees, depending on the fission model, a robust r-process pattern in the region A ≈ 120 − 180 (Mendoza-Temis et al. 2015). The robustness of the 2nd-to-3rd peak rprocess abundance pattern largely independent of the EOS is in agreement with other recent BNS simulations including weak interactions and approximate neutrino Top: Mass-weighted final nuclear abundances (arbitrary units) as a function of mass number A for each simulation, based on all tracers sampling the dynamical ejecta of the respective run. A robust 2nd-to-3rd peak r-process independent of the EOS is obtained. For comparison, solar r-process abundances from Sneden et al. (2008) are shown as black dots. Bottom: Mass distribution as a function of Ye extracted at T = 5GK using tracer particles, showing the collective composition of dynamical ejecta at the onset of the r-process.
transport (Radice et al. 2018;Kullmann et al. 2022). The robustness of the pattern also holds when extended to non-equal mass mergers, which suppress the amount of high-Y e material and increase the neutron-rich tidal ejecta component (and thus the nucleon-to-seed ratio). For all runs, we find actinide abundances at the level of uranium similar to solar abundances.
Some deviations from the solar abundance pattern visible in Fig. 10 are likely the result of nuclear input data for the nucleosynthesis calculations. Since we employ the FRDM mass model for nucleosynthesis, the third r-process peak is systematically shifted to the right (slightly higher mass numbers) for all simulations, which has been attributed to neutron captures after freeze-out combined with relatively slower β-decays of thirdpeak nuclei in the FRDM model (Eichler et al. 2015;Mendoza-Temis et al. 2015;Caballero et al. 2014). The trough in abundances between A ∼ 140 − 170 relative to solar is likely due to the fission fragment distribution employed here, as pointed out by r-process sensitivity studies (e.g., Eichler et al. 2015;Mendoza-Temis et al. 2015).
The first r-process peak is under-produced in all models, which is due to only partially reprocessed ejecta material with respect to the original cold, neutron-rich matter (Y e ∼ 0.05 − 0.1) of the individual NSs (see above). Reproducing the first solar r-process abundances requires a Y e -distribution that extends well above 0.25 for a significant fraction of the ejecta (e.g., Lippuner & Roberts 2015). Given that the mean of the Y e distribution is slightly higher for SFHo and APR, these models have a larger fraction of first peak material, which is however still under-produced with respect to solar values. Light r-process elements in the first-to-second rprocess peak region are preferentially synthesized in the post-merger phase in winds launched from a remnant neutron star and from the accretion disk around the remnant, which can give rise to broad Y e -distributions (e.g., Perego et al. 2014;Lippuner et al. 2017;Siegel & Metzger 2018;De & Siegel 2021). We defer a more detailed discussion on nucleosynthesis including postmerger ejecta to a forthcoming paper.

Fast ejecta and free neutrons
The possibility of tracing individual fluid elements allows us to evaluate the radioactive heating rate within the ejecta in more detail. In particular, we are interested in the fast portion of the ejecta that generates free neutrons, which then β-decay with a half-life of ≈ 10 min and provide additional heating of the material at timescales of up to hours relative to what would be expected from pure r-process heating. Such early excess heating can power bright UV emission known as a kilonova precursor (Metzger et al. 2015). The amount of free neutrons generated by the ejecta sensitively depends on the expansion timescale (i.e., on the ejecta velocity) and the proton fraction Y e . Figure 11 shows the specific heating rate as recorded by each unbound tracer that samples the dynamical ejecta. Irrespective of the stiffness/softness of the EOS, we find a fast v ∞ ∼ > 0.5c and neutron-rich Y e ∼ < 0.2 component of tracers that generate excess heating on a ∼ 10 min timescale (as expected for free-neutron decay with half-life ≈ 10 min), relative to the standarḋ Q ∝ t −1.3 heating rate of a pure r-process. This excess self-consistently obtained from nuclear reaction network Figure 11. Evolution of the specific heating rate calculated with SkyNet for all unbound tracers of the dynamical ejecta and APR (upper panel; color-coded by asymptotic velocity) as well as LS220 (lower panel; color-coded by electron fraction). Thick lines represent mass averages over all tracers, which closely follow the expected ∝ t −1.3 power-law for r-process heating. Black dashed-lines correspond to an analytic approximation to the heating rate of free-neutron decay (Kulkarni 2005). Heating due to free neutron decay is present over a wide range of EOS softness/stiffness and originates in high-velocity (v∞ ∼ > 0.5c) and low-to-moderate Ye ∼ < 0.2 outflows.
calculations is in good agreement with analytic predictions for free-neutron decay (dashed lines in Fig. 11; Kulkarni 2005). Although excess heating from this light, fast ejecta component does not significantly alter the mass-averaged heating rate (thick solid lines in Fig. 11) with respect to its power-law behavior, this excess heating occurs in the outermost layers of the ejecta where even a small number of free neutrons can give rise to observable emission on timescales of hours (Sec. 5.2). Figure 12 shows the detailed mass distributions of free neutrons according to the asymptotic speed of ejecta. We find that almost all material with velocities faster than ≈ 0.6c produces free neutrons. Slower parts of the ejecta, however, also contribute. These slower layers reside deeper within the ejecta and thus contribute to the kilonova precursor at a higher initial optical depth. Their contribution to the total luminosity is thus dim- mer and peaks at slightly later times relative to injection in the outermost (fastest) layers (Sec. 5.2). The total mass of free neutrons is ≈ 2 × 10 −5 M in all runs (Tab. 1). We find that the stiffer EOS simulation LS220 provides a slightly higher heating rate at t ≈ 10 −2 days than the softer EOSs employed here (Fig. 13), which we mostly attribute to free neutron heating. Although the fraction of fast ejecta producing free neutrons in the former model is smaller (and slower) than those of the latter (Tab. 1, Fig. 12), it is more neutron-rich owing to less shock heating (Y e ∼ 0.2 for the fast ejecta in LS220, compared to Y e ∼ 0.3 for SFHo and APR, see Fig. 11). Overall, this results in more free neutrons nevertheless. This conclusion is expected from parametric explorations of r-process nucleosynthesis (Lippuner et al. 2017) (see their Fig. 3) and also supported by recent parametric studies in merger ejecta using SPH simulations   Fig. 12).

ELECTROMAGNETIC SIGNATURES
Although dynamical ejecta are expected to be subdominant with respect to post-merger ejecta for typical BNS systems such as those considered here, but also more generally looking at the expected BNS population as a whole (Siegel 2019(Siegel , 2022, they give rise to a number of distinctive electromagnetic transients that we intend to model here in a self-consistent way based on ab-initio GRMHD simulations. Based on characteristics of dynamical ejecta extracted from our GRMHD simulations (Sec. 3) and detailed nucleosynthesis analysis (Sec. 4), we present a model to calculate the corre- The inset at t = 10 −2 days = 900 s shows the differences of the mass-weighted heating rates due to free neutron decay.
sponding early kilonova lightcurves taking into account the contribution of free-neutron decay (the kilonova precursor), as well as relativistic effects associated with the fast-moving outflow (Sec. 5.1). We also present a model to calculate the non-thermal afterglow emission related to the mildly relativistic shock wave that the fast ejecta gives rise to as it expands into the interstellar medium (Sec. 5.5). We build a multi-angle model that solves a 1D hydrodynamic shock-wave problem in multiple directions (Sec. 5.5.1) and computes the associated radiation of a non-thermal distribution of electrons accelerated within the shock, solving the time-dependent Fokker-Planck equation (Sec. 5.5.2). Results of both models based on our simulations runs and in the context of GW170817 are discussed in Secs. 5.2 and 5.5.3, respectively. We also investigate the potential impact of recently calculated extremely large opacities of highly ionized lanthanides on the early kilonova emission (Sec. 5.3) as well as signatures of high-energy gamma-ray emission due to inverse Compton and synchrotron self-Compton processes in the kilonova afterglow shock (Sec. 5.5.4).

Kilonova Semi-Analytical Model
We calculate kilonova light curves generated by the dynamical ejecta of our simulations using a new semianalytical model. We formulate this model based on the approach of Hotokezaka & Nakar (2020), which, in turn, closely follows the work of Waxman et al. (2019) and Kasen & Barnes (2019). This model allows us to compute the heating rates of r-process elements using experimentally measured values for the injection energies of nuclear decay chains. Our code builds on the public version of Hotokezaka & Nakar (2020) 9 , which we optimized in terms of computational efficiency. We further developed the framework to additionally include (a) arbitrary mass distributions from numerical simulations, (b) heating rates due to free neutrons, (c) arbitrary opacities, (d) flux calculation in different wavelength bands, and (c) relativistic effects. Furthermore, we embed this 1D framework in a generalized infrastructure for computing a multi-angle kilonova signal following Perego et al. (2017). For this work, however, we restrict ourselves to a spherically symmetric (angular-averaged) approximation and defer a multi-angle approach including postmerger ejecta components to forthcoming work. In the following, we briefly outline the main features of the model.
Given an angular-averaged mass distribution M (v) extracted from a BNS merger simulation, we assume homologous expansion of the outflow and we discretize it into velocity shells v i . Each velocity shell is characterized by a mass M (v i ) and grey opacity κ(v i ). We typically use at least 30 velocity bins, which we find is sufficient to obtain converged light curves.
The energy of each ejecta shell is determined by losses due to pdV work associated with the adiabatic expansion of the fluid (E i /t), radiative loses due to photons escaping the flow (L rad,i ), and heating from nuclear decay products that thermalize in the ejecta (Q i ). From energy conservation, the equation for each shell in the comoving frame reads which we integrate using a fourth-order Runge-Kutta algorithm. This requires the computation of the heating termQ i and of the radiative luminosity L rad,i The ejected outflow contains a variety of heavy neutron-rich isotopes synthesized by the r-process as well as free neutrons (Sec. 4). Unstable lighter nuclei will mainly β-decay, while translead nuclei will also generate α-decays and undergo spontaneous fission. The products of these processes, mainly fast-moving electrons, α-particles, and gamma-rays, will in part be absorbed by the outflow. The deposition of energy from charged particles in radioactive decay occurs through collisional ionization and Coulomb collisions. It is determined by the energy-loss rate of fast particles βK st , where K st is the stopping cross-section per unit mass in units of MeV cm 2 g −1 and β is the velocity of the particle in untis of the speed of light c. As in Hotokezaka & Nakar (2020), we take the experimentally measured values of the stopping cross-sections from ESTAR and ASTAR of the NIST database 10 .
If the timescales of energy loss of the charged particles created in radioactive decays are much shorter than the dynamical time of the flow, their energy is entirely deposited in the ejecta. At later times, thermalization is not as efficient anymore, and we need to follow the energy loss of the particles. Let us consider a nuclear element j that injects an electron with energy e 0 at time t 0 . Taking into account adiabatic losses, the energy e j of the electron follows (Kasen & Barnes 2019) where γ(p) is the adiabatic index of charged particles (Nakar 2020) and ρ is the density of a mass shell. Once we have solved for e j (t; e 0 , t 0 ), the deposited energy into the mass shell is obtained by integrating the ionization losses as a function of time, βK st (t; e 0 , t 0 ). In order to calculate the total heating rate of charged decay products at time t, one also needs to take into account particles injected at earlier times as they may not have previously lost all of their energy. Summing over all nuclear species, the total heating rate per unit mass iṡ where N j (t ) is the number of a radioactive element j per unit mass, τ j the mean lifetime of that element, and t0 j the injection time of the oldest non-thermal particle surviving at t. To compute the initial injection energies and mean lifetimes we use the Evaluated Nuclear Data File library (ENDF/B-VII.1; Chadwick et al. 2011). We use abundance distributions of nuclei with A ≥ 70 directly extracted from our nucleosynthesis calculations (Sec. 4; Fig. 10). We compute Eq. (24) separately for β-decay, α-decay, and fission, using the corresponding stopping cross-sections for each process. Thermalization of gamma-rays produced in radioactive decays occurs through Compton scattering, absorption, and pair creation. The fraction of energy that is deposited by gamma-rays is calculated as f γ = 1 − exp(−τ eff ), where the effective optical depth is defined as τ eff = κ eff Σ(t), with κ eff the absorptive opacity, and Σ(t) the mass weighted column density of the mass shell (see Sec. 3.2 in Hotokezaka & Nakar 2020 for more details). The specific heating rate is thenq th,γ = f γ (t)q γ , where, again, we use data from ENDF/B-VII.1 to calculateq γ .
Summing the contributions to the heating rates from α and β-decays as well and from fission by both charged particles and gamma-rays, we obtainq r , the total net heating rate per unit mass of the synthesized r-process elements. Upon adding the contribution of free neutrons that undergo β-decay, the total heating rate of each shell Q i in Eq. (22) is obtained aṡ where X fn,i is the mass fraction of free neutrons in each velocity shell as extracted from the nuclear reaction network calculations based on the simulation runs; for the specific heating rateq fn (t) of free neutrons we set (Kulkarni 2005) q fn (t) = 3.2 × 10 14 erg g −1 s −1 exp (−t/τ N ), which provides a good approximation to the freeneutron decay heating rate in our reaction network calculations (Sec. 4, Fig. 11). Here, τ N 880 s is the mean lifetime of a free neutron.
Having calculated the effective heating rate for all mass shells, we need to estimate how photons escape from each shell. In general, the frequency-integrated (grey) opacity of r-process ejecta sensitively depends on the mass fraction of lanthanides and actinides (owing to their valence f-shells, which result in an increase in line expansion opacities by orders of magnitude with respect to light r-process elements) as well as on the temperature of the outflow. Since we focus on the early light curve evolution from the outermost (dynamical) ejecta, the photospheric temperatures satisfy T > 4000 − 5000 K, at which opacities are largely insensitive to temperature (e.g., Kasen et al. 2013;Tanaka et al. 2020;Banerjee et al. 2020). Recent opacity calculations by Banerjee et al. (2022) find a strong enhancement in the opacity of lanthanide elements when highly ionized at temperatures ≈ 8 × 10 4 K. We discuss in Sec. 5.3 that this temperature regime is not of concern for the ejecta structure of our binary simulations presented here, as the photosphere never reaches this temperature regime. The opacity is largely controlled by the lanthanide content, which, in our scenario, mainly depends on the electron fraction of the ejecta. We thus ignore the temperature dependence and use a 1D-parametrization of the opacity as a function of Y e proposed by Wu et al. (2022), κ(Y e ) = 9 1 + (4Y e ) 12 cm 2 g −1 , which models the expected range of grey opacity, with a characteristic rapid drop around Y e ≈ 0.25 (above which lanthanide production is suppressed; Lippuner & Roberts 2015). For each velocity shell v i , we calculate a mass-averaged electron fraction Y e to obtain κ(v i ) according to Eq. (27). Although there is some (residual) mixing of material (and thus of different Y e components) within and in-between neighboring velocity shells, we find this mapping is robust when applied at different extraction radii. The radiative luminosity of each shell, L rad,i , is calculated by taking into account, approximately, the trapping, diffusive, and free-streaming radiative regimes. Generalizing the single-zone model of Piro & Nakar (2013), each velocity shell is assigned an energy escape fraction as where f Esc,i = erfc( t diff,i /t) is the energy escape fraction from each shell. Here, erfc denotes the complementary error function, t diff,i = τ i v i t/c is the diffusion time, and the escape time t Esc,i of the shell is defined by t Esc,i = min(t diff , t) + v i t/c (see also Eqs. (30)-(32) in Hotokezaka & Nakar 2020 for more details). The optical depth is determined by τ i (t) = ∞ vit κ(r)ρ(r)dr, where ρ(r) is the density of the shell. The radiative luminosity is assigned a blackbody spectrum L rad,i,ν , defined by the effective temperature that is determined by the shell's radius, T eff,i = [L rad,i /(4πσ SB v 2 i t 2 )] 1/4 , where σ SB is the Stefan-Boltzmann constant. The total luminosity in the comoving frame is obtained by adding the radiative luminosity of all shells, L bol = i dν L rad,i,ν . For reference, we also define the location of the photosphere as the shell i with radius R ph where τ i = 1. The corresponding photospheric temperature is then defined as T ph ≡ [L bol /(4πσ SB R 2 ph ))] 1/4 . Finally, since we are dealing with mildly relativistic outflows, we need to incorporate relativistic effects. Since the observed luminosity strongly scales with the Doppler factor, even mildly relativistic ejecta shell velocities can have an appreciable impact on the observed lightcurves. Luminosity and fluxes in the observer frame are computed taking into account all relevant relativistic effects (the relativistic Doppler effect, the time-offlight effect, and relativistic beaming) following the reconstruction method based on an 'energy package' approach presented by Siegel & Ciolfi (2016). The photosphere at a given time t is discretized in angular bins and the emergent radiation from each patch of the surface is further discretized into frequency bins. The resulting energy packages ∆E(θ k , ν l ) are emitted from the instantaneous photosphere in the comoving frame and are received by the observer at their respective arrival times. The observer light curve is reconstructed by binning these packages into observer time and frequency bins accounting for relativistic effects (see Sec. 5.7 in Siegel & Ciolfi 2016 for details). Figure 14 presents the bolometric luminosity of the kilonova emission and neutron precursor emission resulting from the dynamical ejecta of all our runs as computed with the model discussed in Sec. 5.1. Shown are separate calculations that distinguish the effects of heating due to free neutron decay and relativistic effects from pure r-process heating only.

Bolometric light curves
Heating from free-neutron decay enhances the total luminosity at early times (∼ 0.5 h after merger) by a factor of 8-10. Differences to Metzger et al. (2015), who find an enhancement of the order of 15-20, are likely due to their mass of free neutron material being an order of magnitude larger and distributed at larger velocities (v ∼ 0.8c) than what we find in our present simulations (see the discussion in Sec. 3.2). We expect relativistic effects (chiefly Doppler boosting) to enhance the observed luminosity roughly by the Doppler factor to the third power (cf. Siegel & Ciolfi 2016). 11 We crudely estimate the angular averaged enhancement factor to where we used a typical (solid angle averaged) value of cos θ = π/4. This estimate is consistent with Fig. 14 which shows an additional enhancement by a factor of 3-4. In total, at early times ( ∼ < 1 h), we find the bolometric luminosity of this kilonova with an amount of ≈ 2 × 10 −5 M of free-neutron producing material and typical velocity of v ∼ 0.6c (cf. Sec. 4.2; Tab. 1) is a factor ∼ < 10 − 20 higher than a non-relativistic kilonova powered by freshly synthesized r-process elements only. Relativistic effects subside as the photosphere recedes to lower expansion speeds further within the ejecta; we find noticeable enhancements in observed luminosity up to ∼ < 10 h. The neutron precursor signal is expected to peak in the source frame when the free-neutron producing ejecta becomes transparent to photons, t diff ≈ t, which, for typical values in our simulations, can be roughly estimated to t p ≈ 3M ej,n κ 4παvc in good agreement with the peak timescales of the dashed-dotted curves in Fig. 14. Here, we have assumed a steep outermost ejecta profile with M (v) ∝ v −α , where α ≈ 10 (cf. Figs. 19 and 12). Relativistic effects then decrease that timescale roughly to t p,obs = 1 − β 2 t p ≈ 0.8t p ≈ 0.3 h for v = βc = 0.6c, consistent with the shift evident from Fig. 19. The observed peak timescale is thus closer to the timescale of free neutron heating of about 0.25 h ≈ 900 s (cf. Eq. (26)).
Softer EOSs show a brighter transient at peak light with a slightly earlier peak time, while LS220 shows a somewhat broader peak and slightly delayed peak time. While this behavior is expected from the fact that softer EOSs lead to larger velocities of the outermost layers of the dynamical ejecta (cf. the discussion in Sec. 3.2, Figs. 3, 7), the electron fraction of the outflows complicates the picture. Part of the fastest ejecta is shockheated and thus turns less neutron-rich, with softer EOS giving rise to more shock heating. As free neutrons become increasingly buried within slower ejecta shells, their heating acts at higher initial optical depths, which implies longer peak timescales, hence a broader peak, and a dimmer light curve at peak light (cf. LS220 in Fig. 14). We note that despite the smaller ejecta velocities, heating due to free neutron decay can still be significant for stiffer EOS due to a higher neutron-richness in the fast outflows.

Application to GW170817
Given the high photospheric temperatures, the kilonova precursor emission is most prominent in the U-V/U band. The binary configurations explored here are compatible with the binary parameters inferred for GW170817 (Sec. 2). Figure 15 shows the predictions of our kilonova model (Sec. 5.1) for the UV/U/B bands based on the SFHo and LS220 runs when applied at the distance of GW170817 of ≈ 41 Mpc, together with kilonova observations in these bands starting at 11 h after merger as compiled in Villar et al. (2017). Our results suggest that if earlier observations could have been conducted, the neutron precursor signal would likely have been detected. Based on our modeling, we find that the ≈ 2 × 10 −5 M neutron precursor generated by the binary systems typical of Galactic double neutron-star systems considered here can be detected with planned missions such as the wide-field ULTRASAT satellite (Sagiv et al. 2014) 12 out to ∼ < 250 Mpc, assuming a 5σ limiting magnitude of ≈ 23 for a 1 h target-of-opportunity integration at a wavelength of 220-280 nm.
The double-peaked structure in the UV/U/B bands (λ = 268, 365, 445 nm, respectively) evident here is due to the neutron precursor and conventional r-process heating, respectively. Figure 15 illustrates that the relative strength of the first and second peak (i.e. neutron precursor to dynamical r-process ejecta) encodes information about the softness of the EOS, which might prove useful to place additional EOS constraints in future merger events. As also evident from Fig. 15 (cf. the short peak timescales), the dynamical ejecta in our runs is not massive enough to explain the blue kilonova data in GW170817 . This is consistent with previous conclusions that the blue kilonova emission in GW170817 likely requires a substantial contribution from post-merger ejecta (Siegel 2019;Metzger 2020), such as a magnetar-accelerated neutrino-driven wind Mösta et al. 2020;Curtis et al. 2021).

Impact of highly ionized lanthanides on early kilonova emission
Recent opacity calculations by Banerjee et al. (2022) point out that high ionization states of lanthanide elements reached at temperatures ≈ 80000 K can boost the opacity by a factor of 10-100, up to ≈ 10 3 cm 2 g −1 . If 12 https://www.weizmann.ac.il/ultrasat/  indicating that the local gas temperature stays below the temperature regime of ≈ 80000 K for an opacity boost due to highly ionized states of lanthanides (Banerjee et al. 2022).
this temperature regime is realized in the early evolution of merger ejecta ∼ 1 h that contains only a small fraction of lanthanides, Banerjee et al. (2022) show that the boosted opacity can significantly dim the early ∼ 1 h kilonova lightcurve. This affecting, in particular, the early UV brightness is a potential concern for detection prospects of the neutron precursor emission as discussed in the previous sections. Figure 16 shows the gas temperature of our ejecta structure in the SFHo case as a function of velocity coordinate as well as time. Similar behavior applies to our other merger simulation runs. We reconstruct the gas temperature in each ejecta velocity shell i with mass dM i and volume V i from the kilonova model (Sec. 5.1) assuming contributions of radiation and ions to the internal energy E i , where a is the radiation constant, A is the average mass number as obtained from the nucleosynthesis calculations (Sec. 4), and m u is the atomic mass unit (see also Just et al. (2022)). In agreement with the model calculations of Banerjee et al. (2022), we find a 70000 K temperature regime at low velocity coordinates v ∼ 0.2 c and early times ∼ < 1 h, in which high ionization states of lanthanides can be realized. However, the photosphere in our merger ejecta resides at much lower gas temperatures at all times, leading us to conclude that the opacity boost from highly ionized lanthanides is unlikely to significantly affect the early kilonova emission in practice.

Impact of "lanthanides pockets" on the kilonova and its nebular phase
At late times, the kilonova reaches a nebular phase, in which the ejecta material becomes optically thin and the electromagnetic emission driven by radioactivity is dominated by emission lines (Hotokezaka et al. 2021). The spectral signatures of this late-time phase can contain unique information about atomic species synthesized by neutron capture.
Our simulations show that the azimuthal distribution of the electron fraction for soft EOSs is non-uniform because of the preferred direction of core bounces (see Sec. 3.3 and Fig. 8). The ejected material contains "pockets" of neutron-rich material that synthesizes lanthanides and actinides. The opacity of these pockets is thus significantly higher than the surrounding ejecta, which has undergone stronger heating due to shocks during the dynamical merger phase. The resulting "opacity pockets" are illustrated in the upper panel of Figure 17, extracted at a distance of ≈440 km from the merger site. For the stiff EOS, these lanthanide pockets are absent and the opacity distribution is more uniform, see Figure  18.
If the structure of these pockets is preserved to larger distances (late times), this will have an impact on the observable properties of the nebular emission since the dominant nebular lines are determined by elements that dominate the cooling of the plasma. We estimate angular spreading of the ejecta based on the mass-weighted radial ( v r ) and azimuthal ( v φ ) velocities extracted at ≈440 km, as shown in the lower panel of Fig. 17 in  Same as Figure 17, but for the stiff LS220 EOS simulation. The opacity distribution is more uniform in the azimuthal direction and shows the expected latitudinal dependence.
the case of SFHo. In general, we observe that radial velocities in these pockets dominate by at least a factor of ≈3−5 over azimuthal velocities, even though this low-Y e component arises mostly from the tidal tails. For surrounding material this ratio is larger, typically ∼ >6 − 10. As a result of the ratio v r / v φ , material from within the pockets and from the surroundings will azimuthally spread and tend to disperse over time. A rough estimate for the angular dispersion of opacity structures in Fig. 17 is given by where d φ (t) and r(t) denote azimuthal and radial distance traveled by the ejecta, respectively. It is therefore plausible that such structures may persist at timescales of interest for both early photospheric (∼days) and latetime (∼weeks) nebular emission.

Kilonova afterglows
As the mildly relativistic dynamical ejecta expands into the interstellar medium (ISM), it sweeps up ISM material and generates a long-lived blast-wave (Nakar & Piran 2011;Piran et al. 2013;Margalit & Piran 2015;Hajela et al. 2022;Hotokezaka et al. 2018). Within the shock, randomly oriented magnetic fields are generated, amplified, and particles, mainly electrons, are accelerated to non-thermal distributions and generate synchrotron emission (Sari et al. 1998). In this section, we present a model to calculate the dynamics of the shock propagation and the generation of non-thermal radiation directly based on the results of numerical relativity simulations.

Shock dynamics
To describe the hydrodynamical propagation of the blast wave the ejecta runs into the ISM, we assume the ejecta has entered homologous, quasi-spherical expansion with a mass profile M (v) as, e.g., in Fig. 3. We assume this shock sweeps up ISM material, which remains concentrated in a thin slab close to the shock front, where most of the electrons are accelerated. If we assume that the shock is adiabatic, and the EOS of the fluid is trans-relativistic (Mignone & McKinney 2007;Nava et al. 2013;van Eerten 2013), the evolution of the outermost ejecta with velocity β (in units of c) and associated Lorentz factor Γ is determined by energy conservation (see, e.g., van Eerten 2013; Ryan et al. 2020 for more detailed discussion), and the velocityṘ of the shock at radial position R from the merger site is determined by the Rankine-Hugoniot jump conditions, where u = Γβ is the four-velocity. In a merger event, the outflow has a radial velocity structure (Fig. 3). As the fastest part of the outflow starts to decelerate because of its interaction with the ISM, velocity shells deeper within the ejecta inject kinetic energy into the shock region, leading to a so-called refreshed shock (Panaitescu et al. 1998;. In this case, conservation of energy is expressed by where E K (> u) is the kinetic energy of the flow extracted from simulations (see Fig. 19), M 0 = M (u max ) is the mass of the outer-most part of the fluid (we impose a mass cut-off of M 0 = 10 −8 M ), and the second term is the kinetic plus thermal energy of the shocked ISM material that has been swept up by the blast wave. Here, ρ ISM is the constant density of the ISM. The fluid velocity as a function of shock position R can be obtained from the non-linear algebraic Eq. (34), and then the shock Lorentz factor Γ sh (R) can be found using Eq. (33). We note that in the non-relativistic limit, Eq. (34) reduces to the solution used by Piran et al. (2013) and Chevalier (1982).
Employing an ISM density of n = ρ ISM /m p = 0.001 cm −3 , where m p is the proton mass, we compare the evolution of the shock Lorentz factor Γ sh = 1 − (Ṙ/c) 2 −1/2 between the different simulation runs as a function of R in Fig. 20. The (forward) shock propagates at mildly relativistic velocities at early times, with differences of tens of percent between different EOSs and larger Lorentz factors for softer EOS. The shock waves become Newtonian (Γ sh ≈ 1) once reaching their respective Sedov radius, defined as the radius R d within which the rest mass energy of the material within the blast wave equals the energy E K of the explosion, ≈ 2.5 × 10 18 cm E K 10 50 erg 1/3 n 10 −3 g cm −3 −1/3 , in good agreement with the behavior shown in Fig. 20. Here, we have evaluated the expression for typical isotropic equivalent energies of the dynamical ejecta of our simulations (Fig. 19). The shock waves decelerate more gradually than expected from a single-shell scenario (Γ sh ∝ R −3/2 ) due to the fact that the shock is continuously being refreshed by kinetic energy of shells deeper within the ejecta that catch up with the decelerating shock front.  Figure 20. Evolution of the Lorentz factor Γ sh of the shock driven into the ISM (with density n = 0.001 cm −3 ) by the dynamical ejecta as a function of radial distance to the merger site and for the three different BNS mergers simulated here.
The shock waves decelerate more slowly than a single-shell model (dashed line; Γ sh ∝ R −3/2 ) due to the radially structured ejecta refreshing the shock wave at late times.
We assume that a fraction e of the thermal energy e = 4(Γ − 1)Γnm p c 2 generated in the forward shock at time t accelerates electrons into a power-law distribution in Lorentz factor γ e in the frame comoving with the shock (Sari et al. 1998), Here, Q 0 is a time-dependent normalization factor, p is the power-law index, and Θ denotes the Heaviside function that restricts the injection to an interval between the minimum energy given by normalization constraints, and a maximum energy given by balancing acceleration and synchrotron losses (De Jager et al. 1996) γ max Here, e is the electron charge, σ T is the Thomson crosssection, B is the co-moving magnetic field strength, and acc is the ratio of acceleration rate to the maximum possible particle energy-gain rate that we fix to 0.35 following Zhang et al. (2020). Furthermore, we assume a fraction B of the instantaneously generated postshock thermal energy is transformed into magnetic energy (Sari et al. 1998), which results in a co-moving field field given by The evolution of the electron distribution n e (γ e , t) in the shock front is determined by the continuity equation for particle number. We numerically solve the full resulting Fokker-Planck equation for the distribution n e (γ e , t) in the comoving frame of the shock, ∂n e (γ e , t) ∂t = − ∂ ∂γ e [γ e (γ e , t)n e (γ e , t)] +Q(γ e , t) − n e (γ e , t) t esc .
Here, t esc is the escape time of particles from the shock layer andγ e represents the energy gains/losses determined by various radiative processes. The comoving time t here is related to the evolution of the shock as seen in the inertial frame of the merger site with time To solve this equation numerically and to compute the radiation fluxes we use the code Paramo (Rueda-Becerril 2021), taking into account synchrotron and inverse-Compton cooling, together with synchrotron self-Compton/selfabsorption, as well as all relativistic effects required for computing fluxes in the observer frame.
Most studies of non-thermal afterglows from BNS mergers (e.g., Piran et al. 2013;Hotokezaka et al. 2018) employ analytic approximations for the synchrotron spectrum, e.g., a simply connected power law as a function of frequency (Sari et al. 1998). Here, we instead numerically obtain the SED solving for the particle distribution, taking into account self-consistently radiative cooling due to inverse Compton scattering, and the full synchrotron emissivities (for details on the numerical methods, see Rueda-Becerril et al. 2021;Rueda-Becerril 2021;Davis et al. 2022). While providing more accurate predictions for non-thermal emission in general, this approach, in particular, replaces the need to rely on order-of-magnitude estimates for the cooling frequencies, which critically shape the synchrotron light curves .
In order to take the angular dependence of the ejecta structure and viewing angle effects into account, we follow the approach of Margalit & Piran (2015). We assume an axisymmetric outflow and discretize the outflow into angular bins according to polar angle Θ, equidistantly spaced in cos Θ , each bin i having its kinetic energy distribution E K,i (u, Θ i ) (cf. Fig. 7). We solve the shock and Fokker-Planck equations for each angular bin separately, assuming no lateral spreading, and subsequently add the emergent radiation from each annulus in solid angle in the observer frame. Angular variation in the ejecta profiles can shift the peak time of the light curves, although in our equal-mass models the ejecta is sufficiently isotropic for this effect not to be significant (Margalit & Piran 2015). Representative kilonova non-thermal afterglow light curves generated by the dynamical ejecta of our BNS simulation runs, labeled by the corresponding EOSs employed. Also shown are X-ray (blue dots) and radio (orange dot) observations reported by Hajela et al. (2022), who present tentative evidence for X-ray and radio rebrightening following the GW170817 event. Our models are consistent with the observations for typical microphysical parameters (see the text for details).
The free parameters of our model are {p, e , B , n, D L , z, θ obs , u rad }, where D L is the luminosity distance to the merger site, z the corresponding cosmological redshift, θ obs the observer angle relative to the orbital plane, and u rad the external photon field density relevant for inverse Compton scattering. As a fiducial photon field, we use the thermal photons of the Cosmic Microwave Background, with u rad = 0.26 eV cm −3 . As a baseline test of our afterglow code, we consider a spherically symmetric shock powered by ejecta following a power-law distribution in velocity focusing on synchrotron emission only. Comparing to the public code afterglowpy (Ryan et al. 2020) 13 we find very good agreement regarding the synchrotron emission.

Kilonova Afterglow Light Curves and GW170817
We apply our non-thermal afterglow model described above to the observations of Hajela et al. (2022), who present evidence for a late-time (≈ 1000 days postmerger) excess of X-ray and radio emission above the gamma-ray burst broad-band afterglow of GW170817.
13 https://github.com/geoffryan/afterglowpy The existence of this late-time rebrightening is still a matter of debate (Balasubramanian et al. 2022;) and, if present, could have other possible origins as well (e.g., . We fix e to a standard value of e = 0.1 and, as usual, we vary the microphysical parameters (p, n, B ). We find that our models directly based on numerical relativity simulation data are compatible with the observational data obtained so far over a range of typical parameter values; Fig. 21 shows representative models for each EOS model. The stiffer LS220 simulation generates less energetic outflows and thus tends to prefer a somewhat higher density of the ISM as well as higher magnetic energy compared to the softer EOSs in order to reach the same observed fluxes at ≈ 1230 days. The properties of peak flux (timescale and flux level) sensitively depend on the ISM density and the velocity of the fastest (outermost) ejecta shell. We observe that the APR run with the fastest ejecta profile tends to peak earliest, although all simulations considered here produce overall similar results that are consistent with the data.
Due to the mildly relativistic velocity structure of the outflow, the light curves show long peak timescales of thousands of days and associated broad peaks. The Xray fluxes almost plateau around the peak timescale between ∼ 200 to ∼ 2000 days, before decreasing again. As also evident from Fig. 21, the X-ray light curves tend to have a faster rise and decay timescale than the radio signal, which indicates that the cooling frequency has crossed the X-ray band by ∼ 2000 days. The radio emission plateaus for several years.
Compared to Nedora et al. (2021a) who use a large set of numerical simulations and a semi-analytical afterglow model, we find that our model is compatible with similar ISM densities but tends to require higher magnetic energies. This is because the extraction of asymptotic ejecta velocities in our simulations leads to overall slower speeds (see Sec. 3.2). We also find that our model favors a low spectral index, p ≤ 2.1, in agreement with Hajela et al. (2022). Finally, we note that at this early stage of what may represent excess emission, the gamma-ray burst afterglow may still be contributing in part to the observed light curve (not included in the emission model here), and thus the derived microphysical parameters would likely need further adjustment if the gamma-ray burst afterglow component was included self-consistently in this analysis.

High-energy γ-ray emission
Given the mildly relativistic nature of the shock that the ejecta drives into the ISM with Lorentz factor Γ sh ≈ 2 (cf. Fig. 20), and the low-density environment of photons and baryonic matter, gamma-ray emission from inverse Compton and synchrotron self-Compton processes in these types of afterglows is not expected to be detectable unless the merger event occurs sufficiently nearby and/or the efficiency e of converting kinetic energy of the shock into accelerating electrons is very high (Takami et al. 2014).
However, if the merger occurs within a medium of high radiation density of soft photons u rad ∼ 10 − 10 3 eV cm −3 , such as in globular clusters (Venter et al. 2009), for example, there may be noticeable emission of gamma-rays generated by inverse Compton scattering. For u rad ∼ 10 3 eV cm −3 , we find fluxes of ∼ TeV gammarays of νF ν < 10 −14 erg s −1 cm −2 for D L > 40 Mpc. This still is too faint to be detected with the Cherenkov Telescope Array (CTA) or 10 years of Fermi Large Area Telescope (LAT) data (Atwood et al. 2009). However, even in the absence of direct detection of gamma-rays, inverse Compton processes may shape the observed nonthermal afterglow emission by helping cool electrons at later times. We find that the resulting decrease in X-ray emission is significant if u rad ∼ > 10 2 eV cm −3 at timescales of t ∼ > 1 year.
Furthermore, at early times (< 10 − 15 d post-merger) the ejecta shock is accompanied by thermal radiation from the kilonova with radiation density ∼ > 8.7 × 10 5 eV cm 3 L bol 6 × 10 40 erg s −1 t 10 d Here, we have scaled to the bolometric luminosity of the GW170817 kilonova of ≈ 10 42 erg s −1 at 0.5 d (Drout et al. 2017), rescaled according to ∝ t −0.85 up to 5.5 d and ∝ t −1.33 thereafter (Drout et al. 2017), and have used a typical shock Lorentz factor of Γ sh ∼ < 2 (cf. Fig. 20). In a similar manner as in the case of the gamma-ray burst afterglow shock (Linial & Sari 2019), electrons in the kilonova ejecta shock interact with the low-energy kilonova photon field and cool by inverse Compton scattering. We find that the X-ray radiation of the ejecta shock may decrease by tens of percent at very early times (< 10 − 15 d), provided strong kilonova afterglow emission due to sufficiently fast ejecta.

CONCLUSIONS
We have presented a comprehensive model of electromagnetic signatures arising from the dynamical ejecta in BNS mergers. For this purpose, we performed GRMHD simulations that incorporated, for the first time, tabulated three-parameter nuclear EOS, magnetic fields, and approximate neutrino transport via a ray-by-ray transport scheme. We provide a detailed analysis of the characteristics of dynamical ejecta and its ejection mechanisms. We employ the dynamical ejecta profiles extracted from these simulations, characteristics of the dynamical ejecta as recorded by several families of tracer particles, and results of detailed nucleosynthesis calculations with the nuclear network SkyNet to compute the neutron-precursor thermal transient, kilonova emission, as well as the broad-band non-thermal kilonova afterglow. The modeling of electromagnetic counterparts presented here is directly based on ab-initio numerical simulations of the BNS merger process and provides a comprehensive picture of observable phenomena related to dynamical ejecta in a self-consistent framework.
Our main results and conclusions can be summarized as follows: • We performed a detailed analysis of the dynamical ejecta properties and ejection mechanisms of 1.35-1.35 M binaries typical of the distribution of Galactic BNS systems and consistent with the inferred system parameters of GW170817 using three nuclear EOSs that span the range of NS radii as allowed by current observational constraints ( Fig. 2; Tab. 1). We showed that matter is unbound during the decompression phase in quasiradial bounces of the remnant NS double-core structure. We find that for the soft EOS cases, ∼ 75% of dynamical ejecta is unbound in the first three bounces, while for the stiffer LS220 case, material is unbound more secularly in a combination of torques and shock waves sourced by somewhat weaker bounces.
• All binaries generate a fast tail of material (v > 0.6c) that can be traced back to ejection during the first two bounces of the double-core structure (Fig. 5). While the amount of such fast material (few 10 −6 M , Tab. 1) is in good agreement with some grid-based GRHD simulations of ultrahigh resolution (Kiuchi et al. 2017;Hotokezaka et al. 2018), differences to other grid-based GRHD simulations (Nedora et al. 2021a) can be understood in terms of the criteria to extract unbounded mass and to compute the asymptotic velocities (see Sec. 3.2). We conjecture that based on these results and convergence studies (Dean et al. 2021), the amount of fast ejecta should be roughly converged at our grid resolution and atmosphere level. Remaining differences to Newtonian setups and SPH simulations (Dean et al. 2021;Kullmann et al. 2022;Rosswog et al. 2022) remain an open question, but may likely be attributable to differences in the overall setup, such as Newtonian dynamics, velocity extraction, geometric effects due to assumed symmetries, and not sufficiently converged ejecta tails in SPH simulations with ∼ < tens of particles (see the discussion in Sec. 3.2).
Although our present GRMHD setup includes neutrino absorption and employs a lower atmosphere floor than previous grid-based GRHD simulations with tabulated EOS, we are still not sensitive to a ∼ 10 −8 − 10 −7 M ultra-relativistic ejecta envelope as suggested by Beloborodov et al. (2020), who employ this envelope as a breakout medium for a relativistic jet to explain the prompt gammaray emission of GRB 170817A associated with the BNS merger event GW170817.
• We showed that for soft EOSs, the second quasiradial bounce after merger produces a shock wave that reprocesses the cold, previously ejected neutron-rich tidal tail material, resulting in an asymmetric, azimuthal pattern in Y e , giving rise to 'pockets' of neutron-rich material while other angular parts are shock-heated to higher Y e values (Sec. 3.3, Figs. 4 and 8). For stiffer EOS, however, the action of multiple weaker shock waves generates a more homogeneous Y e pattern as a function of azimuthal angle. Once r-process nucleosynthesis proceeds the low-Y e regions in the soft EOS cases translate into high-opacity 'lanthanide/actinide' pockets. This is illustrated in Fig. 17.
• Future work will need to explore the observational consequences of the azimuthal opacity effect ( Fig. 17) with radiative transfer calculations of the early kilonova emission. If preserved to large radii, these "pockets" could also significantly impact the nebular phase of the kilonova once the ejecta becomes transparent at late times. Indeed, the radial velocities of these pockets are larger than their azimuthal velocities by a factor of at least ∼3 − 5, translating into angular dispersion of ∼ < (10 − 20) • (see Eq. (32)). This indicates that these structures could be preserved to large distances and thus both to the early photospheric (∼days) and late-time nebular (∼weeks) emission phases. . The dominant nebular lines are determined by the elements that dominate the cooling (Hotokezaka et al. 2021), i.e. lanthanide/actiniderich pockets may produce different nebular emission relative to other regions with equal mass but less lanthanide content. Such nebular composi-tion diagnostics may be performed with the James Webb Space Telescope. This may provide an additional avenue to infer the softness/stiffness of the EOS .
• Our setup of several families of passive tracer particles injected into the simulation (Sec. 2.5.2) allows us to track and characterize different ejecta components and to perform detailed nucleosynthesis calculations of the dynamical ejecta with the nuclear reaction network SkyNet (Sec. 4). We obtain a mass-averaged abundance pattern in good agreement with solar abundances for A > 80, including, in particular, a robust pattern between and including the 2nd and 3rd r-process peak irrespective of the EOS (cf. Fig. 10). An underproduction of first-peak elements relative to solar for A < 80 evident here is expected due to limited possibilities to increase Y e in dynamical ejecta via neutrino absorption and emission. Softer EOSs lead to slightly higher abundance in the A = 50−80 region due to somewhat stronger shock heating and the resulting increase in Y e .
Directly post-processing with the nuclear reaction network, we find that in all simulations a typical amount of ≈ 2 × 10 −5 M of the fast-ejecta tail leads to free neutrons (cf. Fig. 12; Tab. 1) whose decay provides early heating on the free-neutron decay timescale of ≈ 900 s in excess to that of the rprocess (cf. Fig. 11). We find that although stiffer EOS tend to generate less fast dynamical ejecta, a similar amount of free neutrons can be generated due to relatively more neutron-rich outflows.
• We present a novel kilonova model that incorporates a detailed method for thermalization of radioactive decay products according to Hotokezaka & Nakar (2020), heating due to free-neutron decay, and relativistic effects to reconstruct the observer lightcurve (Sec. 5.1).
• We find that the combination of free-neutron heating and relativistic effects results in an enhancement of the kilonova emission by a factor 10−20 at ∼ < 1 h for all EOS cases (cf. Fig. 14), with relativistic effects leading to noticeable enhancement up to ∼ < 10 h. A moderate shift in the peak timescale of the precursor emission due to relativistic effects to values closer to the neutron decay timescale of ≈ 0.3 h complicates confusion with possible early cocoon heating (Gottlieb et al. 2018;Kasliwal et al. 2017). We find that this modeling of early kilonova and neutron-precursor emission based on ab-initio numerical simulations is compatible with observations of GW170817 in that i) the main blue emission component of the GW170817 kilonova cannot be explained with dynamical ejecta only as found here and likely requires substantially more ejecta from the post-merger phase and ii) earlier observations in UV and optical bands would likely have detected the precursor emission (cf. Fig. 15). The boost in UV/optical brightness by a factor of a few due to relativistic effects found here provides promising prospects for detection of the neutron precursor following future merger events with UV telescopes such as Swift or with planned missions such as the wide-field ULTRASAT satellite (Sagiv et al. 2014). Based on our precursor modeling, we estimate a detection horizon of ∼ < 250 Mpc for ULTRASAT.
• The ejecta structure in our simulations suggests that in typical merger systems such as those investigated here an opacity boost from highly ionized lanthanide nuclei arising at temperatures ∼ > 70000 K as recently reported by Banerjee et al. (2022) is unlikely to affect the early UV/optical emission on timescales ≈ 1 h (Sec. 5.3, Fig. 16).
• Finally, we present a novel model to compute the non-thermal kilonova afterglow emission generated by a shock wave that the fast dynamical ejecta runs into the interstellar medium (Sec. 5.5). We numerically solve the 1D hydrodynamics of a refreshed shock propagating into the interstellar medium and compute the synchrotron and inverse-Compton radiation by numerically solving the corresponding Fokker-Planck equation. We find that for typical shock microphysics parameters and typical densities of the interstellar medium inferred for GW170817, our model based on ab-initio numerical simulations is compatible with the recent data by Hajela et al. (2022), who report tentative evidence for a rebrightening of the broad-band GW170817 GRB afterglow emission emerging at a timescale of 10 3 days post-merger (Fig. 21). The presence of this rebrightening, however, is still being debated (Balasubramanian et al. 2021;).
• High-energy gamma-ray emission arising from inverse-Compton and synchrotron self-Compton processes associated with the kilonova afterglow is unlikely directly observable in future merger events. However, we find that in the presence of intense photon fields (e.g., the kilonova at early times and/or in environments such as globular clusters) indirect consequences such as additional cooling of electrons leave observable imprints in the X-ray radiation (leading to a decrease by ∼ tens % relative to the radio band). A. CONVERGENCE OF EJECTA In BNS merger simulations it is challenging to obtain a satisfactory convergence level of the total amount of dynamically ejected mass at finite (typical) numerical resolutions employed. This is in part due to non-linear effects involved in the ejection process (see, e.g., Radice et al. 2018) and the requirement for and treatment of an artificial atmosphere floor in grid-based simulations; another role is played by the methods of measuring the amount of unbound ejecta in a given simulation due to finite numerical resolution available for the numerical extraction techniques and ambiguities in unboundedness criteria. In this appendix, we explore how the properties of the dynamical ejecta change with resolution and other parameters such as the extraction radius of outflow detectors. We focus on the SFHo case and we compare results from simulations with a grid spacing of ∆x = 180, 220, 250 m, which we label as high, medium, and low resolution, respectively.
The left panel of Fig. 22 shows a convergence test regarding the total amount of dynamically ejected mass using the geodesic criterion at a detector distance of 300 km. We observe convergence at the 10% level at the highest resolution employed. A comparison of the corresponding detailed ejecta mass distributions of the high and medium resolution runs according to asymptotic velocity, electron fraction, and specific entropy are shown in Fig. 23. Except for moderate differences in the Y e distribution around Y e ≈ 0.25, which we attribute to better or worse resolved shock heating, we find good agreement in the detailed ejecta distributions. The value Y e ≈ 0.25 being a threshold for lanthanide and actinide production (e.g., Lippuner & Roberts 2015), these results emphasize the importance of well-resolved shock dynamics during the merger process.
Finally, we compare the amount of unbound mass according to the Bernoulli and geodesic criteria for different detector radii based on the high-resolution run in the right panel of Fig. 22. As discussed in Sec. 3.2, the unbound mass according to the geodesic criterion saturates at earlier times, because it only captures the relatively faster fluid elements that are dynamically unbound during the merger phase, while the Bernoulli criterion additionally captures slower material of high specific enthalpy unbound through winds during the post-merger phase, giving rise to a steadily growing unbound ejecta mass in the early post-merger phase. Early on, both criteria lead to identical ejecta masses, highlighting the dynamical nature of the underlying ejection process. We find almost indistinguishable variations in the total ejecta mass for detectors at radii larger than ≈ 450 km using both the Bernoulli and the geodesic criterion, once accounting for a shift in arrival times of the outflow at the various detector locations.
B. TRACER SAMPLING In Fig. 24, we compare the distributions of dynamical ejecta mass extracted from our SFHo simulation using tracer particles (Sec. 2.5.2) and computing the unbound mass flux through a spherical detector surface (Sec. 2.5.1). Overall, we find excellent agreement, specifically in thermodynamic quantities such as electron fraction and specific entropy, which is key for accurate nucleosynthesis calculations. We find a moderate oversampling of high-velocity ejecta using the tracer approach, which is due to our method of placing a large number of tracers in the high-entropy, under-dense regions of shock heated ejecta during the merger (Sec. 2.5.2). We thus discard tracers below a certain mass-threshold as 'noise' and exclude them from further analyses. In doing so, we improve the velocity distributions shown in Fig. 24 and obtain an overall good agreement and sampling of ejecta according to asymptotic velocity, with still more than  Figure 24. Mass distribution of electron fraction for dynamical ejecta material in our SFHo run, measured separately at 440 km by analyzing the mass flux on the Cartesian grid using a spherical detector surface ('Outflow'; Sec. 2.5.1) and by analyzing the properties recorded by passive tracers advected with the flow ('Tracers'; Sec. 2.5.2). The total number of tracer particles with scale on the right is shown as dot-dashed black lines. There is moderate oversampling of high-velocity ejecta (first panel), which we fix in subsequent analyses by excluding low-mass tracers that we consider sampling 'noise'. Other distributions that are key to nucleosynthesis analyses, such as Ye and specific entropy, show excellent agreement and are not affected by such oversampling, regardless of whether the full (third panel) or reduced (second panel) set of tracers is used. See the text for details.