Magnetized Outflows from Short-lived Neutron Star Merger Remnants Can Produce a Blue Kilonova

We present a 3D general-relativistic magnetohydrodynamic simulation of a short-lived neutron star remnant formed in the aftermath of a binary neutron star merger. The simulation uses an M1 neutrino transport scheme to track neutrino–matter interactions and is well suited to studying the resulting nucleosynthesis and kilonova emission. A magnetized wind is driven from the remnant and ejects neutron-rich material at a quasi-steady-state rate of 0.8 × 10−1 M ⊙s−1. We find that the ejecta in our simulations underproduce r-process abundances beyond the second r-process peak. For sufficiently long-lived remnants, these outflows alone can produce blue kilonovae, including the blue kilonova component observed for AT2017gfo.


Introduction
The observation of AT2017gfo (Arcavi et al. 2017;Coulter et al. 2017;Soares-Santos et al. 2017)-a kilonova arising from merging neutron stars (NSs)-has provided strong evidence that such mergers are a site where heavy elements are produced via the rapid neutron-capture process, or r-process (Kasen et al. 2017;Pian et al. 2017).The optical and infrared spectra of this transient show an early "blue" peak at ultraviolet/optical frequencies and a late "red" peak in the near-infrared (Chornock et al. 2017;Evans et al. 2017;McCully et al. 2017;Nicholl et al. 2017;Tanvir et al. 2017).This behavior is thought to arise from the presence of two or more distinct outflow components, which differ with respect to their total masses as well as velocities and opacities (Cowperthwaite et al. 2017;Drout et al. 2017;Perego et al. 2017;Villar et al. 2017).
The neutron-richness of ejecta is an important quantity for nucleosynthesis, indicated by its electron fraction Y e = N P /N B , where N P is the total number of protons and N B is the total number of baryons.Relatively neutron-rich ejecta (Y e  0.25) typically synthesize a significant mass fraction of heavy elements, including the lanthanides (140 A < 176; Lippuner & Roberts 2015).Lanthanides produce strong optical line blanketing and shift emission toward infrared bands, giving rise to a "red" kilonova (Barnes & Kasen 2013;Kasen et al. 2013;Tanaka & Hotokezaka 2013;Fontes et al. 2015;Even et al. 2020).Less neutron-rich ejecta (Y e  0.25) are lanthanide-poor with a correspondingly lower opacity and thus produce a "blue" kilonova (Metzger & Fernández 2014;Perego et al. 2014;Wanajo et al. 2014).
While the multicomponent fits to AT2017gfo provide reasonable estimates of the ejecta properties, linking them to various outflow components found in merger simulations has been challenging (Shibata & Hotokezaka 2019;Nedora et al. 2021).The red component is thought to arise from the tidal dynamical ejecta and/or accretion disk outflows (Margalit & Metzger 2017;Kawaguchi et al. 2018;Metzger et al. 2018;Siegel 2019), while the exact origin of the blue component remains debatable.The blue component will be the focus of this Letter.
The ejecta properties derived for the blue component include a total mass of ∼2 × 10 −2 M e , high velocities ∼0.27c, and a low opacity compatible with a composition characterized by Y e ∼ 0.25-0.35(Villar et al. 2017).Shock-heated polar dynamical ejecta may have the requisite velocities and composition but are insufficiently massive (Oechslin et al. 2007;Sekiguchi et al. 2016).Accretion disk outflows may produce massive ejecta with the appropriate composition but with much lower velocities (Fahlman & Fernández 2018;Miller et al. 2019).The production of the blue component thus requires the presence of another mechanism.Several possibilities have been proposed to simultaneously fit these estimates, including spiral-wave winds (Nedora et al. 2019), disk winds (Combi & Siegel 2023), and neutrino-heated, magnetically accelerated outflows from a strongly magnetized NS merger remnant (Metzger et al. 2018;Mösta et al. 2020).
For the multimessenger signatures of neutron star mergers, both magnetic fields and neutrino-matter interactions are important.While substantial work has been done to model NS mergers, simulations typically do not employ high-enough resolution to capture the turbulent magnetic field amplification/ evolution during the merger and in the merger remnant (Duez et al. 2006;Price & Rosswog 2006;Anderson et al. 2008;Giacomazzo et al. 2011;Rezzolla et al. 2011;Dionysopoulou et al. 2013;Neilsen et al. 2014;Palenzuela et al. 2015;Ruiz et al. 2016;Ciolfi et al. 2017Ciolfi et al. , 2019;;Ruiz et al. 2019;Ciolfi 2020).Kiuchi et al. (2015Kiuchi et al. ( , 2018) ) performed highresolution GRMHD simulations but the simulations did not include a realistic nuclear equation of state (EOS) or neutrinos.Mösta et al. (2020) modeled an NS remnant with an added initial poloidal magnetic field using 3D GRMHD simulations (including neutrinos via a leakage scheme) and found that turbulence induced by the magnetorotational instability (MRI) in the remnant amplifies the magnetic field to beyond magnetar strengths (10 15 G).They also found an outflow consistent with a magnetized wind, with properties broadly in agreement with those estimated for the blue component of AT2017gfo.An investigation of the emission produced by these wind ejecta in Curtis et al. (2023) found that these outflows produce a red kilonova that peaks on the timescale of 1 day.Recently, Most & Quataert (2023) showed that NS merger remnant outflows can produce flares, jets, and quasiperiodic outflows, and Combi & Siegel (2023) argued that disk winds dominate the kilonova emission, with the NS remnant only being responsible for a blue precursor signal.All of these simulations treated the neutrino-matter interactions using a leakage scheme or onemoment schemes, which do not accurately capture the Y e evolution.Ciolfi & Kalinani (2020) did not include neutrino radiation in their simulations but followed over 250 ms of postmerger evolution, concluding that magnetically driven winds can produce an ejecta component as massive and fast as required by the blue kilonova.
Here, we provide high-resolution dynamical-spacetime GRMHD simulations of short-lived NS remnants expected to form in the aftermath of a merger while also including important microphysics, such as a nuclear EOS and neutrino effects via an advanced M1 neutrino transport scheme.Our primary focus here is on the neutrino microphysics since neutrinos are key to setting the electron fraction of the ejecta, and in turn, the heavy element abundances and kilonova properties (Lippuner & Roberts 2015;Foucart et al. 2016;Radice et al. 2016;Miller et al. 2019;Radice et al. 2022).We perform a 3D GRMHD simulation of a magnetized short-lived NS merger remnant using a two-moment M1 scheme for neutrino transport and find an outflow that reaches quasisteady-state operation in agreement with Mösta et al. (2020).We calculate r-process abundances in the ejecta and predict kilonova light curves.We find that given a sufficiently longlived NS remnant, such outflows alone can produce blue kilonovae, including the blue kilonova component observed for AT2017gfo.
The Letter is organized as follows.In Section 2, we describe our input models as well as numerical methods and codes used to carry out the simulations and to compute abundances and light curves.We present the outflow properties in Section 3.1, the ejecta composition in Section 3.2, and kilonova light curves in Section 3.3.We discuss the implications of our results and future directions in Section 4.

Simulation Setup
We employ the same simulation setup as in Mösta et al. (2020) with the key difference of using a recently developed M1 neutrino transport scheme (Radice et al. 2022).This choice is geared toward investigating the difference of using an improved neutrino transport method and its impact on the outflow properties and observational signatures.The shortlived NS remnant we evolve formed in the merger of an equalmass binary with individual NS masses of 1.35M e at infinity, originally simulated in GRHD by Radice et al. (2018) using the WhiskyTHC code.We map the NS remnant as initial data to our simulation at 17 ms postmerger, adding a poloidal magnetic field of strength B 0 =10 15 G.In the early postmerger evolution, the core of the short-lived remnant NS bounces as it settles down from the merger.This leads to oscillations in the central density that makes mapping the simulation data difficult.We therefore wait until t − t merger = 17 ms to map the simulation.More details regarding the mapping procedure and the motivations behind it can be found in Mösta et al. (2020).
The early postmerger evolution is described in detail in Section 3.1 of Radice et al. (2018), but we summarize the key points here.In the early postmerger evolution the core of the newly formed short-lived remnant NS centrifugally bounces as it settles down from the merger.During that period some material is ejected as part of the dynamical ejecta.The details of the early postmerger evolution are described in Figures 2 and 3 of Radice et al. (2018).Material is ejected from the remnant only in the first few bounces of the remnant's core, and the mass ejection rate is significantly lower than the mass ejection rate measured from the magnetically launched outflows described in this Letter.Additionally, these early outflows are found predominantly in the equatorial region such that we do not expect these dynamics to alter our qualitative conclusions.We choose the magnetic field strength to mimic the field resulting from the MRI and potential dynamo action that typically saturates on the order of 10 15 G due to secondary instabilities destroying channel flows from the MRI.While the field is also amplified in the merger itself by the Kelvin-Helmholtz instability with a saturation of 10 16 G, this amplification will not be able to deliver a global field as assumed here.While we expect the choice of magnetic field to moderately impact the collapse dynamics of the remnant, we have investigated a range of magnetic field strengths and configurations in de Haas et al. (2023) and have found changes in collapse time of 2-3 ms such that we do not expect this to qualitatively change the global remnant evolution.
Our simulation employs ideal GRMHD using the Einstein Toolkit (Löffler et al. 2012;Werneck et al. 2023) and includes the K 0 = 220 MeV variant of the EOS of Lattimer & Swesty (1991).Neutrino-matter interactions are tracked using a recently developed M1 neutrino transport scheme (Radice et al. 2022).We initialize neutrino quantities by assuming equilibrium conditions, where we compute the population of trapped neutrinos assuming thermal equilibrium at the temperature and proton fraction of the simulation.Under these assumptions, the neutrino number and energy densities are given analytically in terms of Fermi integrals and depend only on temperature and chemical potentials.In Mösta et al. (2020), the same NS remnant was evolved using a similar setup, but using a leakage scheme to treat neutrinos.

Tracer Particles and Nuclear Reaction Network
We use tracer particles to extract the thermodynamic conditions of the ejecta.The tracers are uniformly spaced to represent regions of constant volume.At the start of the simulation, each tracer particle is assigned a mass that accounts for the density at its location and the volume the particle covers.We place 96,000 tracer particles to ensure that a sufficient number of tracer particles are present in the outflow.The tracers are advected with the fluid and collected at a surface defined by a radius of r = 150 M e = 222 km, and tracer quantities are frozen once the particle crosses this surface.
We compute the ejecta composition by postprocessing the tracers with the SkyNet nuclear reaction network (Lippuner & Roberts 2017).We include 7852 isotopes up to 337 Cn.The reaction rates, nuclear masses, and partition functions are the same as those employed in Curtis et al. (2023).The network is started in nuclear statistical equilibrium when the particle temperature drops below 20 GK.The network calculates the source terms due to individual nuclear reactions and neutrino interactions and evolves the temperature.The dynamical simulation, and hence the tracer trajectory, ends around 12 ms after the mapping due to the collapse of the remnant.The network continues the nucleosynthesis calculation up to a desired end time by extrapolating the tracer data beyond the end of the trajectory assuming homologous expansion.Additional details about the extrapolation are provided in Appendix A. Our calculations are carried out to 10 9 s, sufficient to generate a stable abundance pattern as a function of mass number.

Radiation Transport
We compute kilonova light curves using SNEC (Morozova et al. 2015;Wu et al. 2022), a spherically symmetric Lagrangian radiation hydrodynamics code capable of simulating the hydrodynamical evolution of merger ejecta and the corresponding bolometric and broadband light curves.As input, SNEC requires the radius, temperature, density, velocity, initial Y e , initial entropy, and expansion timescale of the outflow as a function of mass coordinate.The Y e , entropy, and expansion timescale are used to compute the heating rates, as described in Wu et al. (2022).The effective wavelengthindependent opacity κ is taken to be a function of the prenucleosynthesis Y e : The minimum opacity is 1 cm 2 g −1 , the maximum is 10 cm 2 g −1 , and the above functional form makes the opacity drop steeply near Y e = 0.25.This is motivated by the work of Tanaka et al. (2018) and produces a range of opacity values in keeping with the lanthanide fractions expected for low-, intermediate-, and high-Y e material.
The outflow properties are recorded by measuring the flux of the relevant quantities through a spherical surface at radius r = 100M e = 148 km, using the same approach as employed in Curtis et al. (2023).

Outflow Properties
We present 2D snapshots of the simulation at t = 10.1 ms in Figure 1, along with the corresponding snapshots from the simulation presented in Mösta et al. (2020).As described earlier, the two simulations differ in their treatment of neutrino physics.The ejecta produced by the magnetized wind are concentrated in the polar region, as seen in the panels depicting the Y e of the unbound material.The M1 simulation leads to an earlier remnant NS collapse compared to our previous work (12 ms versus 22 ms).The collapse time of the remnant depends sensitively on the numerical prescription and physics included and is beyond the scope of this Letter to fully investigate.We refer the reader to discussion in Radice et al. (2022) for details of how the collapse time of the remnant can change when including different neutrino transport prescriptions and to Mösta et al. (2020) for a discussion of collapse time with respect to magnetic field evolution and resolution.The overall outflow geometry and mass ejection rate between the two simulations presented here are largely unaffected by the We only plot the Y e for the unbound material.We show snapshots at t − t map = 10.1 ms, where t map = 17 ms after coalescence of the NS binary.The top row shows the evolution of the remnant using an M1 scheme while the bottom row shows a similar evolution using a leakage scheme.The dashed lines are isodensity contours of ρ = 10 7 , 10 8 , 10 9 , and 10 11 g cm −3 .different collapse times and are very similar.Roughly ∼3 × 10 −3 M e of material becomes unbound, and the mass ejection rate for the M1 simulation when the outflow has reached a quasi-steady-state is ;0.08 M e s −1 , roughly consistent with our previous results.
In Figure 2, we present histograms of the radial velocity of the unbound material and its Y e at various times during the two simulations.The unbound material is determined via the Bernoulli criterion −hu t > 1, where h is the relativistic enthalpy of the magnetized fluid.The ejecta show a broad and overall very similar velocity distribution.A significant amount of material has velocities between 0.2c < v r < 0.3c at all times, consistent with the range of velocities estimated for the blue kilonova component.As the system evolves, a small fraction of the ejecta also attains velocities in the range 0.4c < v r < 0.5c.
Examining the Y e distributions, we note the effect of neutrino-matter interactions in the M1 simulations, which drive the ejecta Y e toward higher values over time, as discussed in some detail in Radice et al. (2022) and Zappa et al. (2023).In the leakage simulation, most of the ejected material has Y e values between 0.2 and 0.3, with a peak around Y e ∼ 0.25, but the M1 simulation predicts meaningfully higher ejecta Y e -barely any of the material ejected at later times has Y e  0.3.The leakage scheme used in Mösta et al. (2020) included an approximate energy deposition rate due to neutrino absorption in optically thin conditions but did not directly include the effect of neutrino absorption on the composition.The inclusion of this effect in the M1 simulation results in the increase in Y e observed.
While the two snapshots in time do not correspond to exactly the same stage in the two simulations (due to the different remnant lifetimes) the difference in outflow composition persists across the entire simulated time.The systematically higher Y e values found using the M1 scheme as compared to leakage agree with the results of Radice et al. (2022).This range of Y e values in the outflow is expected to inhibit the synthesis of a significant mass fraction of lanthanides.

r-process Nucleosynthesis
The mean outflow entropy is roughly 20 k B /baryon with typical expansion timescales of the order ∼10 − 20 ms.The lanthanide turn-off point for these ejecta is thus expected to be around around Y e ∼ 0.24 (Lippuner & Roberts 2015).In the top panel of Figure 3, we show the Y e distribution for all ejected tracers when the temperature of the particles is last above 5 GK (the relevant temperature for r-process nucleosynthesis), as computed within SkyNet.This is distinct from the ejecta Y e shown in Figure 2 since the Y e is evolved within the network based on the neutrino luminosity employed during postprocessing.
We plot the Y e distribution for four calculations, each assuming a different constant value of the neutrino luminosity, ranging from 0 to 10 53 erg s −1 .This range of constant neutrino luminosities is used to bracket the possible uncertainties in composition arising from the approximate neutrino transport.For these calculations, we assume L L .The luminosities observed in the simulation are a few 10 52 erg s −1 .These values are in agreement with Cusinato et al. (2022).We will present the detailed neutrino properties in the simulation in a future study.
In general, higher constant neutrino luminosities shift the Y e distribution toward higher values, approaching the equilibrium Y e .The exact value of Y e , and whether weak equilibrium is actually attained, is decided by the competition between the weak interaction timescale and the dynamical timescale.Here, most of the ejected material has Y e  0.3 for all calculations, and for the extreme case where we assume L ν = 10 53 erg s −1 , the Y e peaks around ∼0.44.
In the bottom panel of Figure 3, we show the corresponding abundance patterns.The r-process abundances produced do not match the solar pattern for any of the scenarios, and the abundances of heavy nuclei are increasingly suppressed relative to the solar pattern as we increase the constant neutrino luminosity employed during postprocessing.However, in combination with the dynamical ejecta, the solar pattern can be reproduced.For longer-lived remnants, we expect the contribution of the remnant outflows to dominate the total abundance pattern.
The lanthanide fractions for the four constant luminosity runs, in order of increasing L ν , are 3.1 ×10 −3 , 2.6 ×10 −3 , 1.4 ×10 −3 , and 2.13 ×10 −10 .Typically, for X La  3 × 10 −3 , the kilonova peaks in the near-infrared J band (Even et al. 2020).However, since the Y e of our outflows, and hence the ejecta composition, varies significantly with angle off of the midplane, the nature of the observed kilonova will depend heavily on the viewing angle, in agreement with previous studies based on numerical relativity simulations, radiative transfer simulations, and Bayesian analyses of observational data (e.g., Perego et al. 2017;Kawaguchi et al. 2018;Breschi et al. 2021).

Blue Kilonova
We calculate the kilonova emission using the SNEC code, which is spherically symmetric.However, the ejecta Y e , which directly sets the opacity in our radiation transport treatment, shows a significant dependence on latitude, with relatively low values closer to the equator and values approaching 0.5 closer to the poles (as seen in Figure 1).This behavior clearly cannot be adequately captured by the average outflow profile.We therefore compute the averaged profile of the ejecta in the polar region in addition to the averaged profile of the total ejecta.
Figure 4 presents the velocity, temperature, Y e , and opacity profiles of the outflow, used as input for SNEC.The two rows correspond to the averaged profiles of the total ejecta and the polar ejecta (contained within a 42°polar angle).Examining the total ejecta profile, we find average temperatures around 12 GK and average velocities between 0.16 and 0.2c.The Y e increases as we move inward in mass coordinate, reflecting the increase in the Y e of the ejected material over time.Most of the ejecta have Y e  0.3, but a small amount of material present in the outermost layers of the total ejecta profile has Y e  0.2.Most of this low-Y e material was ejected at early times at lower latitudes.Our opacity prescription assigns opacity values ∼10 cm 2 g −1 to the outermost, low-Y e layer and opacities of ∼1-2 cm 2 g −1 to the bulk of the total ejecta with Y e  0.3.The outermost high-opacity ejecta may serve as a lanthanide curtain (Kasen et al. 2015;Wollaeger et al. 2018;Nativi et al. 2021), absorbing blue light and shifting the kilonova peak to redder bands.Examining the polar ejecta only, we find temperatures around 11 GK and relatively higher velocities between 0.20-0.24c.The corresponding profile shows typical average Y e  0.35, distinctly higher than that of the total ejecta profile.Most importantly, the lanthanide curtain is nearly absent given most of the low-Y e material lies close to the equator.See Appendix B for additional discussion of the Y e structure and evolution as well as the angular cuts.
In Figure 5, we present the AB magnitudes in optical (ugriz) and near-infrared (JHK s ) filters, computed under the assumption of blackbody emission at the photosphere and for the layers above the photosphere.The distance between the observer and the kilonova is taken to be 40 Mpc, same as the approximate distance to AT2017gfo.We present light curves for the total ejecta and the polar ejecta for the simulation presented here, as well as those computed under the assumption of continued mass ejection by a longer-lived NS remnant.
For the total ejecta profile, the brightest emission is seen in the i and z bands.Due to the presence of the lanthanide curtain, by the time these ejecta expand enough to radiate efficiently, they have also cooled down substantially and the majority of the observed emission happens at longer wavelengths.For the polar ejecta, where the outermost layer of lanthanides is nearly nonexistent, the light curve peaks earlier and the brightest emission is seen in the g and r, i.e., "blue" optical bands.If we consider ejecta within a polar angle of 54°, then the light curve peaks in the r and i bands instead.
Owing to the short remnant lifetime, the total ejecta mass in our simulations is roughly 1 order of magnitude lower than that estimated for the blue component of AT2017gfo.For longerlived remnants, we expect the mass ejection to continue at the quasi-steady-state rate until the remnant collapses.We explore this possibility using extrapolated ejecta profiles, assuming appropriate mass ejection rates and mean values for the relevant physical quantities.In a longer-term simulation, the outflow may become faster, less massive, and less neutron-rich over long timescales, as baryons are evacuated from the neutron star atmosphere by neutrino winds (Thompson et al. 2004).However, we do not expect major departures from the average quasi-steady-state properties for a remnant with a lifetime of a few hundred milliseconds.
Assuming the quasi-steady-state outflow persists, we find a range of outcomes with respect to the kilonova peak magnitudes, timescales, and color corresponding to a range of remnant lifetimes.For the case where we obtain ∼2 × 10 −2 M e of ejecta, matching the total mass inferred for the blue component of AT2017gfo, the blue kilonova produced is broadly consistent with the early blue emission observed for AT2017gfo.This is shown in the bottom row of Figure 5.Note that the late time behavior is thought to arise from other ejecta components.

Summary and Discussion
We present a 3D GRMHD simulation of a short-lived NS remnant formed in a binary neutron star merger.We employ a microphysical finite-temperature EOS and an M1 scheme for neutrino transport.A magnetized wind is launched from the remnant and ejects neutron-rich material with a quasi-steadystate rate 0.8 × 10 −1 M e s −1 .
The ejecta Y e distribution is broad, similar to other studies in the literature that employ M1 transport (Radice et al. 2022;Zappa et al. 2023).It peaks at relatively high Y e  0.3 at early times and the peak shifts to Y e  0.4 at later times.Such high Y e values, in combination with the high velocity and mass ejection rate, make magnetically driven NS remnant ejecta a distinct and potentially dominant component of merger ejecta, and a promising source of blue kilonovae (see Mösta et al. 2020 for the importance of MHD for these outflows).
We compute nucleosynthesis yields by postprocessing the ejecta with SkyNet.The production of r-process elements in the ejecta is suppressed relative to the solar pattern.Depending on the neutrino luminosity employed during postprocessing, the ejecta lanthanide fractions vary between 3.1 ×10 −3 to a negligible amount.However, combining the ejecta produced during the course of our simulation with the dynamical ejecta will result in an abundance pattern close to solar.
We map the outflow to SNEC to predict the resulting kilonova.The kilonova color depends on the viewing angle due to the changing ejecta composition as a function of latitude (Perego et al. 2017;Kawaguchi et al. 2018;Breschi et al. 2021).For the polar ejecta, we find peak magnitudes in the "blue" g and r bands.For longer-lived NS remnants, significantly more mass can be ejected by the magnetized wind.Here, we find that a remnant with a lifetime of ∼240 ms can produce a blue kilonova compatible with AT2017gfo.Note that how long the remnant must survive in order to produce sufficiently massive ejecta depends on the mass ejection rate, and is therefore sensitive to the binary properties, the NS EOS, and the physics implemented in the simulation.This is the first demonstration of the production of blue kilonovae compatible with AT2017gfo from NS remnant outflows based on high physical-fidelity GRMHD simulations and including M1 neutrino transport as well as finite-temperature EOS effects.
Since the dynamical merger ejecta synthesize substantial amounts of lanthanides (Radice et al. 2016;Kullmann et al. 2022;Radice et al. 2022;Just et al. 2023), in combination with the NS remnant ejecta, a kilonova with both blue and red components can be produced.While the dynamical ejecta dominate the early emission in the equatorial region, they are not expected to obscure the kilonova produced by the NS remnant outflows found here.These outflows are also much faster than postmerger accretion disk winds.Thus, the production of an early blue kilonova should not be affected by the disk-wind component launched after the NS remnant collapses to a black hole.
In this work, we have produced kilonova light curves for magnetized NS remnant outflows under the assumption of spherical symmetry.In a follow-up work, we will model a longer-lived NS remnant and perform multidimensional kilonova calculations, tracking ejecta from the moment of launch to the kilonova phase.Our work complements ongoing efforts to model populations of kilonovae and determine observing strategies for upcoming surveys such as LSST (Ekanger et al. 2023;Setzer et al. 2023).Such efforts need to take NS remnant outflows into consideration in order to capture the full spectrum of kilonova transients.
Figure 5. AB magnitudes in ugrizJHK s bands.Top two rows: band light curves for the total ejecta and the polar ejecta.From left to right, the band light curves correspond to the ejecta produced during the course of the simulation, and those produced assuming a total remnant NS lifetime of ∼240 ms and ∼1 s, respectively.The black dotted line gives the effective temperature computed at the photosphere.Bottom row: comparison with broadband data for AT2017gfo, for the case where we obtain ∼0.02M e of ejecta.The aggregated broadband data were collected by multiple teams and compiled in Villar et al. (2017).
simulations have been carried out on Snellius at SURF and SuperMUC-NG at GCS@LRZ, Germany.We acknowledge PRACE for awarding us access to SuperMUC-NG at GCS@LRZ, Germany.

Appendix A Nucleosynthesis
Although the tracer trajectories from our simulation end around 12 ms, the nucleosynthesis is not complete at this point.In a longer-term simulation, these ejecta would experience further dynamical and physical evolution, eventually reaching a phase where they are homologously expanding, typically on a timescale of the order of a few seconds.
Here, we take the standard approach in the community and continue the nuclear reaction network calculation up to a desired end time by smoothly extrapolating the particle data beyond the end of the trajectory under the assumption of homologous expansion.The network expands the particle using ρ ∝ t −3 until a minimum temperature is reached.The evolution of the tracer temperature and density within SkyNet is shown for two tracer particle in Figure 6.This approximation may introduce an error into the resulting abundances but is commonly employed in the absence of long-term simulations that follow the ejecta from the moment of launch to the homologous expansion phase.

Appendix B Electron-fraction Evolution
Figure 7 shows how the Y e of the unbound material evolves over time.At early times, there is low-Y e material present at all latitudes, but the lowest values are attained closest to the equator.As the outflow evolves, this low-latitude low-Y e material is no longer seen in the simulation and the magnetized steady-state outflow, concentrated in the polar region, consists almost entirely of higher-Y e material.
When we calculate the kilonova emission from these ejecta in section 3.3, we want to capture not just the average behavior of the ejecta but also that of the (polar) magnetized wind component.Therefore, we also compute kilonovae using the averaged profiles of the material ejected within various polar angular cuts.When we excise material outside a certain polar angle, we omit some combination of low-Y e material ejected early and close to the equator, and potentially some (small) amount of material residing in the wings of the magnetically driven outflow.We retain, however, any low-Y e material that was ejected at any time during the course of the simulation within the chosen polar angle.This is why the lanthanide curtain is not removed completely when we average ejecta contained within a 42°p olar angle in Section 3.3-about 10 −4 M e of high-opacity material remains, affecting the peak timescale of the kilonova.The exact amount of material present in the lanthanide curtain carries some uncertainty as this behavior may depend on the particular system under study and multidimensional mixing also comes into play; however, that investigation is the subject of our future work.

Figure 1 .
Figure1.Meridional slices, i.e., xz-plane with z being the vertical axis, of the density ρ, the Lorentz factor at infinity Γ, the velocity component aligned with the rotation axis v z , and the electron fraction Y e .We only plot the Y e for the unbound material.We show snapshots at t − t map = 10.1 ms, where t map = 17 ms after coalescence of the NS binary.The top row shows the evolution of the remnant using an M1 scheme while the bottom row shows a similar evolution using a leakage scheme.The dashed lines are isodensity contours of ρ = 10 7 , 10 8 , 10 9 , and 10 11 g cm −3 .

Figure 2 .
Figure2.Histograms of the radial velocity v r of the unbound material (top panels), where r is the radius in spherical coordinates, and the Y e of the unbound material (bottom panels), shown at different times during the dynamical simulation (4.2, 8.5, and 12.7ms), for the leakage and M1 simulations.Each column corresponds to a different time in the evolution.

Figure 3 .
Figure 3. Mass-weighted abundance as a function of mass number for the NS remnant outflow.The different colors correspond to different constant neutrino luminosities employed during postprocessing.The solar abundance pattern has been scaled to match the second r-process peak at A ∼ 130 for the zero luminosity setting.The dashed lines show the total abundance pattern produced in combination with the dynamical ejecta from Radice et al. (2018).

Figure 4 .
Figure 4. Input profiles for SNEC showing averaged quantities for the NS remnant outflow as a function of mass coordinate.The top two panels show the radial velocity (solid blue line) and temperature (dashed black line) of the total ejecta on the left, and the corresponding Y e (solid blue line) and opacity (dashed black line) on the right.The bottom two panels show the same quantities plotted for the polar ejecta.

Figure 6 .
Figure6.The evolution of density (black lines) and temperature (blue lines) for two sample tracers.The solid lines show the original tracer particle data, while the dashed lines show the data extrapolated by SkyNet under the assumption of homologous expansion.This evolution corresponds to the case where zero neutrino luminosity is employed during postprocessing.

Figure 7 .
Figure 7. Meridional slices, i.e., xz-plane with z being the vertical axis, showing the time evolution of the Y e for the unbound material.We show six snapshots at t − t map = 1.02, 2.04, 4.08, 8.05, 10.1, and 11.91 ms, where t map = 17 ms after coalescence of the NS binary.