Can jets make the radioactively powered emission from neutron star mergers bluer?

Neutron star mergers eject neutron-rich matter in which heavy elements are synthesised. The decay of these freshly synthesised elements powers electromagnetic transients ("macronovae"or"kilonovae") whose luminosity and colour strongly depend on their nuclear composition. If the ejecta are very neutron-rich (electron fraction $Y_\mathrm{e}<0.25$), they contain fair amounts of lanthanides and actinides which have large opacities and therefore efficiently trap the radiation inside the ejecta so that the emission peaks in the red part of the spectrum. Even small amounts of this high-opacity material can obscure emission from lower lying material and therefore act as a"lanthanide curtain". Here, we investigate how a relativistic jet that punches through the ejecta can potentially push away a significant fraction of the high opacity material before the macronova begins to shine. We use the results of detailed neutrino-driven wind studies as initial conditions and explore with 3D special relativistic hydrodynamic simulations how jets are propagating through these winds. Subsequently, we perform Monte Carlo radiative transfer calculations to explore the resulting macronova emission. We find that the hole punched by the jet makes the macronova brighter and bluer for on-axis observers during the first few days of emission, and that more powerful jets have larger impacts on the macronova.


INTRODUCTION
The first joint detection of gravitational and electromagnetic waves from a binary neutron star merger on 17 August 2017 marked the beginning of a new era of astrophysics (Abbott et al. 2017a,b,c). About two seconds after the peak of the gravitational wave (GW) signal, a short gamma-ray burst (sGRB) was detected (Goldstein et al. 2017;Savchenko et al. 2017) and the remarkable event was followed up during the subsequent days and weeks all across the electromagnetic spectrum, starting with early UV, optical and IR signals (e.g. Arcavi et al. 2017;Cowperthwaite et al. 2017;Evans et al. 2017;Drout et al. 2017;Kasliwal et al. 2017;Pian et al. 2017;Smartt et al. 2017;Soares-Santos et al. 2017;Utsumi et al. 2017) later followed by X-rays (Troja et al. 2017;D'Avanzo et al. 2018;Margutti et al. 2018) and radio emission (Alexander et al. 2017;Hallinan et al. 2017;Troja et al. 2019). The UVOIR signal, observed from the first day up to two weeks, was broadly consistent with the emission expected from a macronova (or kilonova, hereafter MKN) (Kasen et al. 2017;Perego et al. 2017;Tanaka et al. 2017;Rosswog et al. 2018), a thermal transient powered by the radioactive decay of the freshly synthesized r-process elements (Li & Paczyński 1998;Kulkarni 2005;Rosswog 2005;Metzger et al. 2010;Metzger 2017). The emission in X-rays and the radio band was interpreted as being produced by a relativistic jet, observed slightly off-axis, interacting with previously ejected matter (Alexander et al. 2017;Margutti et al. 2017;Hallinan et al. 2017;Kasliwal et al. 2017;Lazzati et al. 2018;Margutti et al. 2018;Mooley et al. 2018;Lamb & Kobayashi 2018;Kathirgamaraju et al. 2018;Ghirlanda et al. 2019;D'Avanzo et al. 2018). The origin of the sGRB is so far unknown, but could plausibly have been produced inside the jet, or in a shock breakout event when the jet emerges from the ejecta (Nakar & Piran 2017;Lamb & Kobayashi 2017;Gottlieb et al. 2018a,b;Beloborodov et al. 2018).
Soon after the discovery of the first neutron star binary (Hulse & Taylor 1975) it was realized that such binaries would be driven through GW emission towards a violent collision which could potentially eject neutron-rich matter (Lattimer & Schramm 1974). It remained, however, an open question for more than two decades whether such extremely gravitationally bound systems can actually eject any mass at all. The first hydrodynamic-plus-nucleosynthesis calculations (Rosswog et al. 1998;Rosswog et al. 1999;Freiburghaus et al. 1999) showed that ∼ 0.01 M of neutron-rich matter is dynamically ejected during a merger and the nuclear network calculations demonstrated that the extremely neutron-rich ejecta effortlessly reproduce the elements up to and beyond the 3rd r-process/"platinum" peak (A=195) (Freiburghaus et al. 1999). This had been a major challenge for all other suggested r-process production sites. These results immediately triggered the discussion on how a neutron star merger would appear electromagnetically (Li & Paczyński 1998).
While these early studies demonstrated the viability of neutron star mergers as a major r-process site, they identified only one ejection channel: "dynamical ejecta" which are tidally flung out by gravitational torques. Since they are never substantially heated, these ejecta carry their original β−equilibrium electron fraction from the original neu-tron star, Ye ≈ 0.05, and this enormous neutron-richness allows them to undergo a "fission cycling" process (Goriely et al. 2011;Korobkin et al. 2012a) which produces a very robust r-process abundance distribution close to the solar pattern for A ≥ 130, but hardly any lighter r-process elements. Oechslin et al. (2007) pointed out that there is a second channel of mass ejection that also happens on a dynamical time scale: shock-heated matter from the interface where the stars come into contact. As of today, many more mass ejection channels have been discussed: matter that becomes unbound on secular time scales (∼ 1 s) from the post-merger accretion torus (Metzger et al. 2008;Beloborodov 2008;Fernandez & Metzger 2013;Fernandez et al. 2015;Just et al. 2015;Siegel & Metzger 2017Miller et al. 2019a;Fernandez et al. 2019), as MHD-driven winds (Siegel & Ciolfi 2015) and by viscous effects (Shibata et al. 2017;Radice et al. 2018a;Shibata & Hotokezaka 2019) from a long-lived neutron star merger remnant. Similar to the case of proto-neutron stars, the enormous neutrino luminosities (> 10 53 erg s −1 ) after a neutron star merger can also drive substantial matter outflows (Ruffert et al. 1997;Rosswog & Ramirez-Ruiz 2002;Dessart et al. 2009;Perego et al. 2014;Martin et al. 2015;Radice et al. 2018b). The secular torus ejecta contain approximately 40% of the initial torus mass and, although the latter may vary substantially from case to case, they likely contribute the lion's share to the total ejecta mass. Due to their different thermal histories and exposure times to neutrinos, the ejecta channels can have different electron fractions Ye and therefore different nucleosynthesis yields 1 . For electron fractions below a critical value Y crit e ≈ 0.25 (Korobkin et al. 2012a;Lippuner & Roberts 2015) lanthanides and actinides are efficiently produced which, due their open f-shells, have particularly high bound-bound opacities Tanaka & Hotokezaka 2013;Tanaka et al. 2020) and therefore lead to red transients that peak days after the merger. Ejecta with electron fractions above Y crit e , in contrast, only produce "lighter" elements with lower opacities and thus result in bluer transients that peak after about one day. Opaque, low-Ye ejecta blocking the view on high-Ye ejecta can lead to a "lanthanide curtaining" effect Wollaeger et al. 2018) which will efficiently block blue light. Therefore it is important to understand the layering, dynamics, interaction and potential mixing of different ejecta channels.
The multi-messenger detection of GW170817 provided evidence that neutron star mergers can produce short GRBs 2 . Given the expected complexity of the matter distribution engulfing the remnant, it is interesting to understand under which conditions a relativistic jet can successfully drill through the ejecta cloud (Murguia-Berthier et al. 2014;Beniamini et al. 2020) and whether/how it affects the mixing/interaction between the different components. This could have substantial consequences for the layering and in-teraction of the ejecta and it can have potentially large effects on the "lanthanide curtaining".
Here we explore the hydrodynamic interaction of a relativistic jet with previously launched neutrino-driven winds from a long-lived neutron star merger remnant. Contrary to earlier studies (Zhang et al. 2003;Mizuta & Aloy 2009;Mizuta & Ioka 2013;Nagakura et al. 2014;Murguia-Berthier et al. 2014;Duffell et al. 2018;Harrison et al. 2018) we use actual simulation results (Perego et al. 2014) as initial conditions for the surrounding wind structure and dynamics. We are particularly interested in the question how jets of different power impact the observable MKN broadband light curves. The light curves are obtained by running 3D Monte-Carlo radiative transfer simulations on a homologously expanding matter background. For this matter background we use the results of our wind-plus-jet simulations and we add an additional component that is meant to represent the likely important secular disk ejecta which cannot be modelled self-consistently together with our large-scale jet simulation. The jet interaction with the ejecta produces an additional contribution, the cocoon emission. While propagating the jet inflates a pressurised cocoon that leads to an additional electromagnetic signal on similar time scales as MKN, and over a relatively wide range of viewing angles (Gottlieb et al. 2018a). This contribution is not considered in the present work.
Our paper is structured as follows. We begin in Section 2 with an overview of our numerical methods and briefly describe the relativistic adaptive mesh-refinement hydrodynamics code amun and the radiative transfer code possis. Our simulation setup is explained in Section 3 and we present and discuss our results in Section 4. A concise summary is offered in Section 5.

Special-relativistic AMR hydrodynamics
The evolution of an ideal, relativistic fluid is governed by the conservation of baryon number and four-momentum. The corresponding equations form a set of hyperbolic equations which can be written as: These equations are solved for a set of six conserved variables u = (D, S i , E, Xa) T , which are related to the physical variables proper rest-mass density, velocity, pressure and passive scalar q = (ρ, v i , p, xa) T by the relations where Γ = (1 − v i vi) −1/2 is the local Lorentz factor, h the specific enthalpy and we have used units in which the speed of light c = 1. The system is closed by assuming an adiabatic equation of state (EoS) p ∝ ρ γ for which we use an ideal gas law so that the specific enthalpy reads In the current set of simulations the specific heat ratio is assumed to take the constant value for a radiation dominated gas γ = 4/3.
We perform this study with amun (https://gitlab. com/gkowal/amun-code), a parallel, special-relativistic, Eulerian (magneto-)hydrodynamics code. The evolution scheme follows a Godunov approach based on an oct-tree hierarchical Cartesian structured grid with adaptive mesh refinement (Quirk 1991;DeZeeuw & Powell 1993). We further use a 3rd order Strong Stability Preserving Runge-Kutta time integration algorithm (Gottlieb et al. 2011) with the Courant-Friedrichs-Lewy (CFL) parameter Press et al. (1992) set to 0.5, a second order TVD reconstruction with the MinMod limiter (Toro 2009), and an HLL Riemann solver (Harten et al. 1983) to compute the fluxes between adjacent cells.

Radiative transfer with the POSSIS code
Broad-band light curves for the models investigated in this study are calculated with the 3-D time-dependent Monte Carlo radiative transfer code possis (Bulla 2019). Assuming homologous expansion for the ejected material, possis simulates Monte Carlo photons propagating throughout the expanding ejecta and interacting with matter via either electron scattering or bound-bound interactions (bound-free and free-free processes are subdominant at the wavelengths investigated in this study, Tanaka et al. 2018). Synthetic observables including spectral energy distributions (SEDs) and light curves can be extracted for different viewing angles defined by their polar (θ obs ) and azimuthal (φ obs ) angles (where z is the jet direction and xy is the orbital plane). We use analytic functions based on state-of-the-art calculations (Tanaka et al. 2018) for the wavelength-and timedependence of opacities (Bulla 2019). Compared to simulations in Bulla (2019), we adopt an improved version of possis where the temperature is no longer parameterized, but rather estimated from the mean intensity of the radiation field at each time and in each zone. We refer the reader to Bulla (2019) for more details about the possis code, more details about our simulations follow in Section 3.2.

Relativistic hydrodynamics
As initial conditions we use a simulated neutrino-driven wind model from Perego et al. (2014). At ≈ 105 ms after the first contact between the two neutron stars the amount of mass ablated by the wind is ≈ 2×10 −3 M . At this stage, the winds have not yet reached a complete steady state and they could still be evolved further (Martin et al. 2015). Here we use a simulation methodology that is different from Perego et al. (2014): they used Newtonian hydrodynamics with selfgravity, a spectral neutrino leakage scheme and a tabulated, nuclear EoS, while our simulations are special relativistic with a point-mass source of gravity, no neutrino transport and an adiabatic EoS. Both simulation methodologies therefore find slightly different equilibria, with the result that our wind model is initially slightly out of equilibrium. One of the major deliverables of Perego et al. (2014)'s wind simulation is the spatial distribution of the electron fraction Ye. We start with exactly this Ye distribution and advect it with the flow as a passive scalar. Both the dynamics and the electron fraction distribution for this initial configuration are shown in Fig. 1. The pressure normalisation is given by the original wind simulation. The whole system is evolved in the gravitational field of the hyper-massive neutron star, which we approximate as a Newtonian point mass of MHMNS = 2.7 M located at the origin.

Computational hydrodynamics grid
The computational grid is roughly shaped as a box (−10 5 ≤ x ≤ 10 5 , −10 5 ≤ y ≤ 10 5 , 0 ≤ z ≤ 1.6 × 10 5 , measured in km) that is located above the equatorial plane (at z = 0). Since the merger can be assumed to be symmetric about the equatorial plane, we use a reflective boundary condition at the bottom of the computational domain and outflow boundaries elsewhere. Around the wind (r 2 × 10 3 km) and the rotational axis (r 5 × 10 2 km) we fix the numerical resolution to the highest refinement level (resolution length ≈ 6 km), and we use lower resolutions elsewhere. Before outflowing matter can reach a boundary, we re-map the matter configuration into a larger domain to allow further, unhindered evolution. The simulations are stopped once the wind material is expanding roughly homologously (at t 1 s). Once this stage is reached, we can scale the matter homologously to the larger distances so that it can be straight-forwardly used in the radiative transfer simulations (see Section 3.2).

The ambient medium
The evolution of the original wind (Perego et al. 2014) has been performed with a background medium of ρ amb = 5 × 10 3 g cm −3 , a value set by the bottom value of the tabulated nuclear EoS that is used in those simulations. Since we are using a simple polytropic equation of state this bound does not apply here and we embed the wind initial data in a steeply decreasing background density with ρ(r) = ρ0 (R0/r) 4 , where ρ0 = 10 −6 g cm −3 and R0 is close to the upper wind boundary at 2000 km. We choose this environment with the purpose of reaching quickly and without discontinuities very low matter densities. This profile is steeper than profiles of stationary winds (∝ r −2 ), and it has no impact on the results as long as its energetic contribution in the system is negligible compared to the one from the wind and the jet. We set the electron fraction inside the background material to an artificially high value of Y amb e = 1 so that it is easily identifiable throughout the entire simulation.

The jet
We assume here that a relativistic jet has already formed and has reached a height of 40 km above the remnant, and from this point onward propagates into the wind (Gottlieb et al. 2018a;Harrison et al. 2018;Mizuta & Aloy 2009;Mizuta & Ioka 2013).
We model the jet as an unmagnetized conical outflow with an opening angle of θ0 = 5 • . We inject it through inflow boundary conditions close to the grid origin. The jet is parameterised by its (total) luminosity Lj, its initial, Γ0, and asymptotic Lorentz factor Γ∞ = h0Γ0, where h0 is the Table 1. Jet parameters for the current simulations: initial opening angle θ 0 , height of the injection region z 0 , effective jet duration (from launching) ∆t inj , number of cells covering the injection region in the x direction N inj,x , initial luminosity L j and initial Γ 0 and asymptotic Lorentz factor Γ∞.
initial specific enthalpy. The three components of the speed are obtained from the jet geometry and Γ0, while the density is obtained from where Σj is the cross-sectional area of the jet at the top of the injection region z0. The pressure in the inlet region is set by the previous parameters together with the EoS as: As for the ambient medium, we set the electron fraction within the jet to Ye,j = 1 to keep it easily recognisable at later times. The jet injection is kept at full power from the beginning of the simulation. The jet is injected from the beginning of the simulation. Because of the very specific initial conditions around the launching region a way to recognise if the jet manages to propagate is required. To do so we choose to set a minimum Lorentz factor Γ = 2 and a height above the origin. Once the head has reached that height we recognise the jet as "launched", and the jet is pushed further for a time of ∆tinj = 100 ms. After that time the luminosity decays exponentially.
We run three simulations: two with jets of different luminosities (Lj = 10 49 erg s −1 , Jet49, and Lj = 10 51 erg s −1 , Jet51) and, as a reference case, we evolve in one simulation (Wind) only the wind without injecting an additional jet. Our chosen values for Lj, Γ0 and Γ∞ are representative for low-and high-luminosity jets in GRBs (Fong et al. 2015). All our jet parameters are listed in Table 1.

Radiative transfer setup
Radiative transfer simulations are performed for the models Wind, Jet49 and Jet51 introduced in Section 3.1.3. In particular, the grid domain is restricted for the two jet models to be the same as in the Wind model, with a maximum velocity of ∼ 0.35 c (maximum spatial coordinate of ∼ 10 5 km at 1 s after the merger) 3 . Each model grid is symmetrised about the orbital plane and downgraded to a uniform Cartesian grid with 128 3 cells. As mentioned in the introduction, the potentially dominant contributions to the ejecta come from rather slow matter parts that are unbound on a secular time scale from the accretion torus. Their simulation is computationally extremely expensive and to date the properties of such outflows are still not entirely settled. There is agreement that about 40% of the initial torus mass becomes unbound (Just et al. 2015;Siegel & Metzger 2017Miller et al. 2019b;Fernandez et al. 2019), but there is no consensus yet, whether the outflows are lanthanide-rich (Siegel & Metzger 2018) or lanthanide-poor (Miller et al. 2019a). We therefore model this ejecta component as a spherical mass distribution with a density profile as in equation (10) of Wollaeger et al. (2018) and total mass of 0.072 M (corresponding to 40% of our initial torus mass). We release this mass 1 s after merger with velocities distributed between 0.03 to 0.1 c. Since the composition is currently unsettled, we perform each time two simulations once assuming a lanthanidepoor and once a lanthanide-rich composition. We note that the disk ejecta contribute to ∼ 99.5, 99.5 and 98.8 % of the total mass for the Wind, Jet49 and Jet51 model, respectively. The final grid is expanded to ti = 0.1 d and densities are scaled as ∝ t −3 according to homologous expansion. The simulations with lanthanide-poor (-rich) compositions are carried out with N ph = 2 × 10 7 (1.5 × 10 7 ) photons. Viewing-angle dependent SEDs are computed from ti = 0.1 to t f = 20 d after the merger and used to create synthetic ugriz light curves from different orientations. All the models investigated are sufficiently close to axial symmetry about the jet axis z and symmetric about the orbital plane (xy) by construction. Therefore, we restrict our analysis to orientations in the xz plane (φ obs = 0 • ) and extract light curves for N obs = 11 observers from the jet axis (θ obs = 0 • ) to the orbital plane (θ obs = 90 • ) equally spaced in cosine (∆ cos θ obs = 0.1). Different opacities are assumed (Bulla 2019) depending on whether the composition of the ejecta is lanthanide-poor (Ye ≥ 0.25) or lanthanide-rich (Ye < 0.25). The nuclear heating rates are taken from equation (4) Fig. 2 (right). The jet undergoes a strong first collimation shock and stays collimated after breaking out from the ejecta. (The leading shock is an artifact from our chosen density and pressure gradients in the ambient medium, but carries essentially no mass and energy and therefore has no impact on the simulation.)

Hydrodynamic evolution
The hydrodynamic evolution is broadly consistent with what is expected from the literature (Zhang et al. 2003;Mizuta & Aloy 2009;Bromberg et al. 2011;Mizuta & Ioka 2013;Nagakura et al. 2014;Murguia-Berthier et al. 2014;Gottlieb et al. 2018a;Duffell et al. 2018) and the aforementioned theoretical framework provides all that is needed to understand the main features of jet propagation.
Our initial conditions as obtained from the neutrinodriven wind simulations of Perego et al. (2014) have a peculiar Ye distribution: their leading edge is made of very low Ye material that could potentially block emission as a "lanthanide curtain" Wollaeger et al. 2018). It is interesting to see what happens to this effective "high opacity skin" during the further hydrodynamic evolution and what its impact is on the electromagnetic emission.
At the beginning of our simulations, the matter contains expanding wings that flow outwards at angles of ∼ 30 degrees from the rotation axis (see Fig. 1). Furthermore, matter without centrifugal support is falling towards the central neutron star along the rotational axis. The downfall compresses the gas in the injection region, which initially hampers a jet launch. A successful launch requires the total jet momentum flux at the jet head to overcome the gas pressure and ram pressures, Pj,ram > max(Pgas, Pgas,ram). In nature, this condition may be reached once enough energy has been deposited to push the gas aside or the central neutron star collapses into a black hole so that the region around the rotation axis is cleaned from polluting baryons (MacFadyen & Woosley 1999). In our simulations, the jets only manage to overcome surrounding pressure after times of tinj ≈ 5 ms in the Jet51 model and tinj ≈ 50 ms in the Jet49 model, approximatively 100 and 150 ms from the first contact respectively. These numbers are broadly consistent with the estimates of Beniamini et al. (2020) for jet formation time scales.
As shown in Fig. 2, both jet models fill pressurized cocoons, wider and more energetic in the high-luminosity case. Despite their (two orders of magnitude) different jet powers, there are successful breakouts in both cases. We follow the approach of Duffell et al. (2018) to estimate which minimum energy is needed for a successful jet break out. Using our simulation results 4 , we find a threshold value of Ecrit ≈ 5×10 45 erg, i.e. even for our low luminosity case with Ejet = 10 48 erg Ecrit, a successful breakout is expected. When breaking out, the jets push aside the high opacity skin and make the inner, lower opacity regions more accessible to potential observations. This will be discussed in more detail in the next section.
In passing we also want to briefly mention jet collimation by the surrounding high-pressure cocoon. Both jets start with a constant opening angle, and their initial conical shape gets roughly cylindrical during propagation, after the collimation from the surroundings. To estimate whether collimation is expected or not, one can resort to the dimensionless jet luminosity parameter (Bromberg et al. 2011): 4 We use the coefficient calibrated on simulations in Duffell et al. (2018), there named κ ≈ 0.05.L ≡ ρjhjΓ 2 j ρa .
For the Jet51 model we findL 5 θ −4/3 0 and therefore at least one collimation shock is therefore expected, consistent with the Lorentz factor map in Fig. 3. The jet is significantly decelerated as a result of a very strong first collimation shock. After the breakout, the jet remains at first collimated, but starts later to spread sideways, consistent with the results by Mizuta & Ioka (2013) and Nagakura et al. (2014). In model Jet49, collimation is even stronger and it happens earlier and closer to the injection region (Matzner 2003;Bromberg et al. 2011;Mizuta & Ioka 2013;Harrison et al. 2018). It changes quickly from conical to cylindrical and stays very narrow until breakout, experiencing multiple recollimation shocks. In this model, a shocked jet is observed spreading immediately after the breakout and returning to a quasi-conical configuration. Further studies of these features are left for more dedicated future investigations.

Radiative signatures
In this Section, we present broad-band light curves extracted with the radiative transfer code possis for the Wind, Jet49 and Jet51 models. In Section 4.2.1 we present the simulations with a lanthanide-poor disk while those with a lanthanide-rich disk are discussed in Section 4.2.2. Fig. 4 shows ugriz light curves predicted for the wind and jet models assuming a lanthanide-poor composition for the disk ejecta. In general, we do see variations in both the brightness and viewing-angle dependence between the Wind model (left panels) and the two models with jets (middle and right panels). As we describe in the following, these differences can be understood by a close inspection of the opacity maps, which are given in the top panels of Fig. 4 for each model. In the Wind model, material at high latitudes (close to the z axis) is characterised by lower Ye and thus higher opacities from lanthanides compared to material in the orbital plane. As a consequence, radiation can escape more easily in the orbital plane than through this low-Ye "curtain" close to the rotational axis. This effect leads to a clear viewingangle dependence in the light curves, with orientations close to the orbital plane (cos θ obs = 0, dark red light curves) associated to brighter MKNe compared to those along the z axis (cos θ obs = 1, dark blue light curves). Since bound-bound line opacities become larger when moving to shorter wavelengths (Bulla 2019), the viewing-angle effect is stronger in bluer filters. For instance, the peak brightness in the u band is ∼ 1 mag brighter for an observer in the orbital plane compared to one on the z axis. This difference is, instead, small (∆m 0.3 mag) in the redder iz filters.

Lanthanide-poor disk
In the Jet49 and Jet51 models, the low-Ye "curtain" close to the rotation axis in the Wind model is "punched" away and decreased in density by the jet (see "arm"-like structures in the top panels of Fig. 4). This effectively reduces the opacities of the material surrounding the disk  ejecta, leading to a clear imprint of the jet on the final observables. Since the "punch-away" mechanism is restricted to regions at high latitudes, an increase in brightness is found for orientations close to the rotation axis (blue lines). As highlighted in the left panels of Fig. 6, this increase is re-stricted to the first ∼ 3 days after the merger while negligible at later epochs when ejecta become optically thin outside the spherical viscous ejecta. While the jet makes MKN brighter at high latitudes, the increase in brightness is nearly absent for orientations in the orbital plane (red lines), therefore Predictions are shown as a function of time since merger, for models with lanthanide-poor (left) and lanthanide-rich (right) disks. Top panels refer to the Jet49 models, while bottom panels to the Jet51 models. The increase in brightness due to the jet punching is stronger for more powerful jets, bluer filters and lanthanide-rich disks.
decreasing the viewing-angle dependence seen in the Wind model. It is worth noting that this would correspond to a jet introducing a viewing-angle dependence for an initial spherical distribution of the winds. In addition, we find that the increase in brightness and decrease in viewing-angle dependence is stronger moving from redder to bluer filters due to line opacities being larger at shorter compared to longer wavelengths. Hence, MKNe are made bluer by the jet in the first ∼ 3 days after merger, with ∆(g − r) ∼ 0.1 mag and ∆(g − i) ∼ 0.3 mag (see Fig. 6). Fig. 6 also highlights how the "punch-away" mechanism and the corresponding impact on the observables are stronger for the Jet51 model characterised by a more powerful jet. For an observer along the jet axis (polar view), the increase in the gr filters in the first day after the merger is of ∼ 1 − 1.5 mag for the Jet51 model while of ∼ 0.5 − 1 mag for the Jet49 model.
The bottom panels of Fig. 4 include ugriz light curves for the MKN associated with GW 170817, AT 2017gfo (Andreoni et al. 2017;Arcavi et al. 2017;Cowperthwaite et al. 2017;Kasliwal et al. 2017;Pian et al. 2017;Smartt et al. 2017;Tanvir et al. 2017;Valenti et al. 2017;Evans et al. 2017). While the overall peak brightness predicted by all the models is comparable to that observed in AT 2017gfo, we note a clear difference in the light curve evolution. In particular, models are systematically fainter than the data in the first day after the merger (especially in gri optical filters) and systematically brighter afterwards. Fig. 5 shows ugriz light curves predicted for the wind and jet models assuming a lanthanide-rich composition for the viscous ejecta. The behaviours identified in Section 4.2.1 for the lanthanide-poor case are visible in the lanthanide-rich case as well. Namely, we find a modest viewing-angle dependence in the bluer filters of the Wind model, while an increase in brightness and decrease in viewing-angle dependence when a jet is launched in the system. As shown in Fig. 6, the increase in gri magnitudes for a polar viewing angle are relatively similar between jet models with a lanthanide-poor and lanthanide-rich disk. The increase in brightness is limited to the first ∼ 3 days and of the order of ∼ 0.5 mag in the Jet49 and ∼ 0.5 − 1.5 mag in the Jet51 models in the first day.

Lanthanide-rich disk
Although the increase in brightness due to the presence of a jet is similar regardless of the disk composition, we do see a clear difference between the two sets of models. In particular, a lanthanide-rich composition of the disk ejecta is associated with higher opacities in the inner regions of the ejecta. The predicted MKN light curves are therefore much fainter than in the lanthanide-poor case, in all filters and especially after peak. As a consequence, the MKN light curves predicted by the Wind, Jet49 and Jet51 are fainter than those observed in AT 2017gfo, especially at early epochs. Although our models are not tailored to provide good fits to AT 2017gfo, we do note that the Jet51 model with a lanthanide-rich disk produces the closest agreement to the observed MKN. In particular, the predicted light curves peak at ∼ 0.5 − 1 mag lower magnitudes in all filters and have very similar decays to those observed (cf with the lanthanide-poor case, Fig. 4).

DISCUSSION AND CONCLUSIONS
We have performed a set of 3D special-relativistic simulations where we study the interaction of jets with the neutrino-driven winds from a neutron star merger remnant. Our particular focus is on the question how this interaction impacts on the resulting MKN light curves. The hydrodynamic simulations have been performed with the special relativistic AMR hydrodynamics code amun and the radiative signatures are subsequently extracted with the Monte-Carlo radiative transfer code possis.
The electron fraction Ye plays a crucial role for neutron star merger ejecta: it determines the nucleosynthesis, the nuclear heating rates, the opacities within the ejected matter and therefore the properties of the resulting electromagnetic transients. Despite of this, most previous studies started from highly idealised and mostly guessed initial conditions. In our study we use a realistic post-merger environment based on the neutrino-driven wind simulations of Perego et al. (2014) as initial conditions. These winds provide us with information about the dynamics and in particular with a peculiar and non-trivial Ye distribution. There are downflows along the original binary rotation axis towards the central neutron star which create a high-pressure environment that initially hampers jet formation so that we need to trigger the jets in our numerical models for a while before they are finally launched.
We explore two different types of jet models with Lj = 10 49 erg s −1 (Jet49) and Lj = 10 51 erg s −1 (Jet51) representative of typical low and high luminosity GRBs. As a reference we have also run a model where only the ejecta are evolved and no jet is injected (model Wind). The dynamical evolution of our jet models is generally in good agreement with the expectations from the literature (Zhang et al. 2003;Mizuta & Aloy 2009;Bromberg et al. 2011;Mizuta & Ioka 2013;Nagakura et al. 2014;Murguia-Berthier et al. 2014;Gottlieb et al. 2018a;Duffell et al. 2018;Murguia-Berthier et al. 2020). Both our models show recollimation shocks, but especially the features in the Jet49 model suggest a strong collimation which might potentially leave an imprint in the early gamma-ray signal.
We use the radiative transfer code possis (Bulla 2019) to predict viewing-angle dependent MKN light curves for the Wind, Jet49 and Jet51 models. The radiative transfer simulations are performed on a matter background for which we use the results from our relativistic hydrodynamics simulations enhanced by a matter component that models the secular ejecta from a central accretion torus. In summary, we find the following: • The models with no jet (Wind models) show some viewingangle dependence in the predicted light curves, being fainter for orientations close to the rotation axis. This effect, caused by low-Ye/high opacity material ("curtain") near the jet axis, is stronger when moving to shorter wavelengths (e.g. ∆m ∼ 1 mag in the g-band while ∆m 0.3 mag in the i−band).
• In the Jet49 and Jet51 models, the low-Ye/high opacity "curtain" close to the axis is "punched-away" and decreased in density by the jet. Compared to the Wind model, radiation can therefore escape more easily along the jet axis, an effect that is stronger in the near-UV and at short optical wavelengths where line opacities are typically higher; • As a consequence, the presence of a jet makes MKNe brighter and bluer in the first ∼ 3 days after the merger for observer orientations close to the jet axis. The presence of a strong jet seems to erase the viewing angle dependence seen in the Wind model, hence making the emission appear more isotropic.
• The increase in brightness is stronger for a more powerful jet, with the Jet49 and Jet51 models being ∼ 0.7 and 1.5 mag brighter than the Wind model in the first day after the merger in g and r filters; • Since to date there is no consensus about the composition of secular ejecta from the inner accretion torus, we have each time performed a simulation with a lanthanidefree and a lanthanide-rich composition. Although we had no intention to specifically model the first detected neutron star merger GW 170817/AT 2017gfo, it is worth stating that the models with a lanthanide-rich inner disk are fainter and in better agreement with the observations than models with a lanthanide-poor disk.
The jet, with its quasi-radial velocity profile and high pressure, tends to transport material away from the axis, which decreases the on-axis opacity. The jet also gives energy to the material that it interacts with, so that the material reaches a higher ballistic speed, which translates to a lower density at the time of the MKN emission, and therefore a lower opacity. However, in this work we have not considered the details of jet formation. There exists a possibility that the process of jet formation could drag with it large amounts of high opacity material that would otherwise stay bound, effectively forming an additional merger mass-loss channel. The net effect of such mass-loss on the MKN emission is not known. This process could be studied through jet formation simulations, and is beyond the scope of this paper.
Clearly, these findings depend on the presence of a thin low-Ye and high-opacity "skin layer" (∼ 10 −5 M ) in our initial data. 3D neutrino-hydrodynamic simulations of neutron star mergers, of their remnants and emerging winds are major computational challenges of contemporary astrophysics and today's results, where comparable, do not (yet) agree on all aspects. Therefore, it is justified to ask how real this lanthanide curtain is.
The original neutrino-driven wind simulations of Perego et al. (2014) were obtained with the FISH code (Kaeppeli et al. 2009) enhanced by an advanced leakage scheme (Perego et al. 2016) that accounts for neutrino absorption in optically thin conditions. As initial condition, served merger remnant configurations (Price & Rosswog 2006) obtained with the SPH code MAGMA (Rosswog & Price 2007). These latter simulations accounted for weak interactions and in particular for Ye-changes due to electron-/positron captures, but they did not include neutrino absorption (Rosswog & Liebendörfer 2003). Nevertheless we believe that these simulations provide a fair representation of what happens in nature to Ye inside the remnant in the early post-merger evolution. The initial torus is formed from matter coming from the outer core of the merging neutron stars. The disc formation timescale is expected to be smaller than the timescale over which weak interactions operate inside the disc. Thus, weak interactions hardly have time to substantially change Ye, which stays close to its initial value, i.e. Ye 0.1 (Perego et al. 2019).
As the accretion process onto the massive neutron star continues over ∼ 100 ms, neutrinos emitted from the central remnant and inner torus regions get continuously absorbed by matter lying further out. They change the electron fraction and deposit their energy there, thus driving a wind with increased Ye. Inside the wind, a spatial gradient in the electron fraction is expected: material that has first reached large distances from the remnant, such as the outer torus layers, has smaller chances to absorb neutrinos and is therefore likely closer to its original, low Ye value. This low-Ye matter is pushed outwards by deeper layers that have captured neutrinos and therefore forms the leading edge of the wind.
While we think that there are good physical reasons for its presence, we cannot safely exclude a numerical origin of this thin lanthanide curtain or its presence only for a subset of cases. We are not aware of other studies that have discussed this layer, but this is not surprising given that a) only few studies include the relevant neutrino physics and reach comparable time scales (∼ 100 ms), b) the layer is at the leading edge of the ejecta which first exits the outer boundaries of the computational domain in a Eulerian simulation and may therefore easily go unnoticed and c) it may only contain a small amount of mass (∼ few 10 −5 M in our case), hard to track through Lagrangian tracers or fix mesh refinements. But even such a small amount can, as we have demonstrated here, have substantial observational consequences.
We close by noting that we did not attempt to include the very neutron rich (Ye ∼ 0.05), first ejected tidal dynamical ejecta Bauswein et al. 2013;Hotokezaka et al. 2013;Rosswog 2013;Radice et al. 2018a) which cover mostly the binary orbital plane. This exploration is left for future work. While most of the presented results are of qualitative nature, they clearly illustrate that seemingly "small details" can have substantial impacts on observable signatures. As a corollary, this implies that we have to expect a large variety of electromagnetic transients after the merger of two neutron stars.