Axisymmetric Radiative Transfer Models of Kilonovae

The detailed observations of GW170817 proved for the first time directly that neutron star mergers are a major production site of heavy elements. The observations could be fit by a number of simulations that qualitatively agree, but can quantitatively differ (e.g., in total r-process mass) by an order of magnitude. We categorize kilonova ejecta into several typical morphologies motivated by numerical simulations, and apply a radiative transfer Monte Carlo code to study how the geometric distribution of the ejecta shapes the emitted radiation. We find major impacts on both spectra and light curves. The peak bolometric luminosity can vary by two orders of magnitude and the timing of its peak by a factor of five. These findings provide the crucial implication that the ejecta masses inferred from observations around the peak brightness are uncertain by at least an order of magnitude. Mixed two-component models with lanthanide-rich ejecta are particularly sensitive to geometric distribution. A subset of mixed models shows very strong viewing angle dependence due to lanthanide “curtaining,” which persists even if the relative mass of lanthanide-rich component is small. The angular dependence is weak in the rest of our models, but different geometric combinations of the two components lead to a highly diverse set of light curves. We identify geometry-dependent P Cygni features in late spectra that directly map out strong lines in the simulated opacity of neodymium, which can help to constrain the ejecta geometry and to directly probe the r-process abundances.


INTRODUCTION
The origin of the rapid neutron capture (r-process) elements is one of the longest-standing unsolved problems in Corresponding author: O. Korobkin korobkin@lanl.govnuclear astrophysics with the most popular candidate sites being stellar collapse and compact object mergers.Neutron star mergers have long been argued to be sources of rprocess elements (Lattimer & Schramm 1974;Eichler et al. 1989;Rosswog et al. 1998;Freiburghaus et al. 1999;Rosswog et al. 1999), mostly because their neutron richness effortlessly leads to the production of platinum-peak elements which has been a major challenge for other production sites.
Such rare events (compared to supernovae) that eject large amounts of r-process per occurrence are also supported by geological evidence (Wallner et al. 2015;Hotokezaka et al. 2015).However, without direct observational evidence our understanding of the enrichment of the Universe by r-process elements was primarily based on theoretical models (for a review, see Côté et al. 2017).
The high r-process yields estimated from GW170817 combined with the large merger rate predicted by its detection (Abbott et al. 2017a), indicate that neutron star mergers could be the dominant site for galactic r-process elements (Côté et al. 2018;Rosswog et al. 2018) and, in fact, could actually produce more r-process than is observed.Current chemical evolution models (Côté et al. 2019;Wehmeyer et al. 2019), however, argue that an additional component is needed, at least at early times, to find satisfactory agreement with the observed r-process abundances.
In a supernova context it has been known for a long time that inferring masses from the optical emission is difficult and that a given light curve can be matched by a wide range of ejecta masses (see De La Rosa et al. 2017, for example).Differences in the modeling methods (e.g.opacity implementations, transport schemes) and the applied microphysics (e.g.opacities, shock physics, equilibrium assumptions) both contribute to these errors.The situation is similar for GW170817, where a range of masses has been inferred by different groups (Côté et al. 2018;Ji et al. 2019).To better observationally constrain neutron star merger yields, kilonova emission models need to be refined, with a particular emphasis on exploring what role the ejecta geometry plays in shaping the electromagnetic emission.This is the subject that we address here.
have revealed different ejecta channels and morphologies.A number of models have invoked two (Tanaka et al. 2017;Kawaguchi et al. 2018;Tanvir et al. 2017;Evans et al. 2017;Troja et al. 2017) or even three (Perego et al. 2017;Cowperthwaite et al. 2017;Villar et al. 2017) ejecta components to explain the electromagnetic observation of GW170817.A typical kilonova scenario with multi-component ejecta mechanisms is sketched in Figure 1.The components may include toroidal tidally-expelled highly neutron-rich ejecta, which produces the red (nIR-IR) emission in the kilonova, or a medium neutron-rich wind emerging from the surface of transient hypermassive neutron star and/or an accretion disk (see Shibata & Hotokezaka 2019, for a review on mass ejection in kilonovae).
It appears that the models ignoring the multidimensional character of kilonova ejecta tend to systematically overestimate the amount of ejected mass.On the other hand, the models with true multidimensional treatment of radiative transfer (Troja et al. 2017;Tanvir et al. 2017;Kawaguchi et al. 2018) seem to favor lighter and faster neutron-rich ejecta accompanied by a heavy "wind" with moderately neutron-rich composition (show in Fig. 4 of Ji et al. 2019).This scenario is also supported by numerical simulations (e.g.Perego et al. 2014;Radice et al. 2016;Dietrich et al. 2017;Shibata & Hotokezaka 2019).
The neutron-rich dynamical ejecta, or the neutron-rich outflows from the post-merger accretion disk (Miller et al. 2019a) are expected to contain a large fraction of highopacity lanthanides (Barnes & Kasen 2013;Tanaka & Hotokezaka 2013;Fontes et al. 2015Fontes et al. , 2020;;Gaigalas et al. 2019).The properties of this lanthanide-rich outflow dic-tate how much optical/infrared light will be blocked; in some cases, lanthanide-rich geometric configurations can produce an effective lanthanide "curtain" (Kasen et al. 2015;Wollaeger et al. 2018).In such cases, the kilonova can be highly sensitive to the viewing angle.
Kilonova light curves and spectra are influenced by both ejecta morphology and microphysics, for example complicated, wavelength-dependent opacities (Kasen et al. 2017;Even et al. 2019), nuclear heating (Lippuner & Roberts 2015), and thermalization (Barnes et al. 2016;Rosswog et al. 2017;Hotokezaka & Nakar 2019).Opacities, nuclear heating and thermalization are all sensitive to the density which is specified by geometric distribution of mass.As we shall see, deviations from the spherical shape are responsible for modified light curves not only via modified surface area or "curtaining", but mainly through an indirect effect on the microphysics because of their different density.
Previously, Kasen et al. (2017) studied different ejecta geometries: one with a light r-process composition focused around the axis to represent wind, and one with a heavy rprocess composition and oblate ellipsoidal morphology to represent dynamical ejecta.In these models, the dynamical ejecta, with an ellipsoid axis aspect ratio of 4, the angular variation is on the order of a factor of 2. In the wind ejecta, when focused into ∼ 45 • conical regions about the axis, Kasen et al. (2017) find an angular variation of ∼20% in peak luminosity.The relatively minor angular variation in peak luminosity of the wind is a result of the less extreme variation in the projected surface area of the wind, relative to the dynamical ejecta, so that the high lanthanide opacity of the low electron-fraction ejecta played a diminished role.Recently, Darbha & Kasen (2020) published a detailed Monte Carlo radiative transport study of several non-spherical geometries in gray-opacity approximation.Zhu et al. (2020) investigated several analytically prescribed 3D morphologies using simplified transport models with effective gray opacity.
In this paper, we study the importance of morphology in combination with its effects on microphysics, focusing on the geometrical distributions that depart from spherical symmetry.Ejecta morphologies are generated from the family of Cassini ovals in 2D axisymmetric geometry.For these morphologies, we demonstrate the enhancement to the radiative flux when departing from spherical symmetry.We then proceed to construct two-component models, representing different possible configurations of low electron fraction (Y e ) dynamical ejecta/accretion disk outflows, and high-Y e wind.For each of the components, we employ a respective composition and nuclear heating representative of low-Y e or high-Y e post-nucleosynthetic conditions.Here we demonstrate the effects of radiative interaction between two superimposed components: lanthanide curtaining, photon reprocessing and photon redirection.Our other recent work provides prelimi-nary studies focusing specifically on the effects of composition (Even et al. 2019) and additional energy sources (Wollaeger et al. 2019).
This paper is organized as follows.In Section 2, we discuss the numerical methods and approximations made in our simulations.In Section 2.3, we describe the functional forms and compositions of our model ejecta.In Section 3, we present light curves and spectra of several of our one-and two-component models, and discuss the effect of morphology and superimposing of the components on these observables.We conclude in Section 4 with a summary of the results and discussion of the possible impact of morphology on the interpretation of observations.In Appendix A, we provide absolute AB magnitudes for each model, in r and J bands at days 1, 4, and 8.In Appendix B, we give a detailed analysis of the radiative structure of single-component axisymmetric morphologies and different factors affecting their electromagnetic emission.

METHODS
In this section, we review our approach for modeling kilonovae, which involves computing multi-wavelength opacity tables and using them to simulate radiative transfer through homologously expanding ejecta.This code uses the same method and opacity implementation as described in Wollaeger et al. (2018);Even et al. (2019); Wollaeger et al. (2019); Fontes et al. (2020).Below we briefly review this methodology and then focus primarily on the improvements to the opacity and radiative transfer introduced for this paper.

Radiative transfer modeling
We employ the multidimensional Monte Carlo radiative transfer software SuperNu1 (Wollaeger & van Rossum 2014) to synthesize the spectra from our models.SuperNu uses a Monte Carlo method for thermal radiative transfer that is semi-implicit in time, which introduces an effective scattering term (Fleck Jr & Cummings Jr 1971).A discrete diffusion Monte Carlo (DDMC) optimization is also implemented to replace multiple small scattering steps with larger jumps between spatial cells (Densmore et al. 2012;Abdikamalov et al. 2012;Cleveland & Gentile 2014).DDMC has recently been generalized to include leakage probabilities out of lines for Lyman-α transport (Smith et al. 2018).The geometries available are spherical, axisymmetric, and Cartesian, each in 1, 2, or 3D.In the calculations presented here, we have incorporated an improved treatment of Doppler shift in DDMC (Wollaeger et al 2020, in prep), that is consistent with the opacity-regrouping algorithm described by Wollaeger & van Rossum (2014).This Doppler shift algorithm has permitted more accurate application of DDMC on the binned opacities described by Fontes et al. (2020).
Originally intended for supernova transients (van Rossum et al. 2016;Kozyreva et al. 2017;Wollaeger et al. 2017), SuperNu has been applied to modeling kilonovae in 1D and 2D axisymmetric geometry (Kasliwal et al. 2017b;Troja et al. 2017;Tanvir et al. 2017;Evans et al. 2017;Wollaeger et al. 2018Wollaeger et al. , 2019;;Even et al. 2019).These studies varied the mass and velocity of wind and dynamical ejecta components but always assumed a spherical wind superimposed with ejecta from SPH simulations (Rosswog et al. 2014).The morphologies we explore in this paper more systematically test the multidimensional functions of SuperNu.

Choice of composition and opacities
We use tabulated opacities generated with the LANL suite of atomic physics codes (Fontes et al. 2015).As with Fontes et al. (2017) and Wollaeger et al. (2018), the tables are calculated under the assumption of local thermodynamic equilibium (LTE), where ionization and excitation populations can be obtained from Saha-Boltzmann statistics.For this work, we use an updated approach to the tabulation of opacity, which employs the binning of lines instead of smearing them over multiple wavelength points (Fontes et al. 2020).In conjunction with higher group resolution in SuperNu this approach permits more resolved spectra for lower-velocity outflows.For the radiative transfer, the frequency values for each density and temperature point of these tables are mapped to a logarithmic 1024-group wavelength grid from 10 3 to 1.28 × 10 5 Å.
All models in this study use one of the two compositions, representing a low-and high-electron fraction (Y e ) abundance pattern.The low-Y e abundance pattern follows the Solar rprocess residuals and is the same as in Even et al. (2019), their Fig.1b.It includes all lanthanides, uranium representing the actinides, and lighter r-process elements represented by Se, Br, Te, Pd, Zr and Fe.The high-Y e abundance pattern is the "wind 2" composition studied in Wollaeger et al. (2018) (see their Fig. 3).Throughout this paper, we will refer to these two cases as "low-Y e " and "high-Y e " compositions.Sometimes, for simplicity, we refer to these components as "dynamical ejecta" and "wind", respectively, but in general, low-Y e material does not necessarily come from the dynamical ejection (Fernández & Metzger 2013), while high-Y e material does not have to be generated in the wind (Wanajo et al. 2014).
Nuclear composition influences the amount and type of radioactive nuclear heating produced in the ejecta over time (Lippuner & Roberts 2015;Zhu et al. 2018Zhu et al. , 2019)).We use detailed time-dependent nuclear heating output from a nucleosynthesis network which distinguishes different radiation species (as in Wollaeger et al. 2018).

Ejecta Geometries
We adopt a suite of ejecta geometries motivated by numerical simulations of dynamical ejecta from binary neutron stars and winds from accretion disks and postmerger hypermassive neutron star.We focus on single-and two-component models with uniform composition across individual components.

Single-component models
Although the variety of possible morphologies is limitless, we have picked a few morphologies that are motivated by the astrophysical setting and easy to set up. Figure 2 depicts the five baseline axisymmetric morphologies: two prolate (H "hourglass", P "peanut"), one spherical (S), and two oblate (B "biconcave", T "torus").The density for each model is given in velocity space axisymmetric coordinates {v r , θ} using the level function for a family of Cassini ovals: The value q = 0 gives a "figure-eight" shape, which if rotated around the z-axis produces our "H" and "T" morphology for the prolate and oblate cases, respectively.We use q = 1 for the "P" and q = 1.5 for the "B" morphologies (see Fig. 2).The spherical morphology "S" is specified by the "cubed inverted parabola" density profile, introduced in Wollaeger et al. ( 2018): In equations (1-2), ρ 0 sets the density scale, vr = v r /v 0 , and v 0 sets the velocity scale.These dimensional constants are connected with the ejecta mass m ej and root mean square expansion velocity v ej as follows: where the nondimensional shape factors C v and C ρ are listed in Table 1 for each of the morphologies.Models with the same mass and average expansion velocity can possess very different density profiles.The latter are with R being the cylindrical radius): "hourglass" (H), "peanut" (P), "sphere" (S), "biconcave" (B) and "torus" (T).The contours are equi-spaced by 0.5 dex in log space, down from the density maximum.The spherical and prolate morphologies are motivated by simulations of post-merger outflows (e.g.Perego et al. 2014;Martin et al. 2015), which in general tend to concentrate towards the polar regions.Figure 4 shows a snapshot of material density in velocity space for a wind from a post-merger accretion disk, computed in a general relativistic magnetohydrodynamics simulation (Miller et al. 2019b,a).The snapshot is taken approximately 100 ms after the black hole formation.This wind outflow is mainly driven by magnetic and thermal forces in the disk.Although the wind has a highly elongated shape, we expect it to become more spherical due to heating by r-process, actively operating in the first seconds  (Miller et al. 2019b,a).Outflow from the disk is highly aspherical. of expansion (Rosswog et al. 2014).Another case of prolate morphology is the "squeezed" collisional ejecta in the polar regions (Oechslin et al. 2007;Sekiguchi et al. 2016;Kasen et al. 2017).The oblate morphologies are motivated by simulations of tidally expelled dynamical ejecta, which are naturally concentrated towards the merger plane (Rosswog et al. 2014;Sekiguchi et al. 2016).
Our one-component models assume uniform composition and are constructed by rotating Cassini ovals around vertical axis as specified above.Although we expect multiple components in the true outflows, by simulating these onecomponent models, we can better understand the impact of geometry, such as the effects of the surface area or highdensity regions.This allows for a better assessment of the impact of superimposing morphologies (e.g.lanthanide curtain).In addition, these non-spherical one-component models can be directly compared to other models in the litera-ture (such as those studied in Perego et al. 2017;Darbha & Kasen 2020;Zhu et al. 2020).

Two-component models
We next combine our one-component morphologies to have a suite of two-component (mixed) models, with overlapping high-Y e and low-Y e compositions.In the mixed models, we only pick three basic morphologies: S, P and T (combinations with H instead of P and B instead of T produce very similar results).The selected combinations are outlined in Fig. 5, with the low-Y e and high-Y e compositions shown in red and blue, respectively.In the model names, the first and second letters stand for morphology of the low-Y e and high-Y e components, respectively, and the number designates the first significant digit in the mean velocity of expansion for the corresponding component.More exactly, "1" designates mean expansion velocity of 0.1 c, "2" designates 0.2 c, and "5" designates 0.05 c.For example, T1S2 stands for a model with toroidal low-Y e component with mean expansion velocity 0.1 c and spherical high-Y e component with mean expansion velocity 0.2 c.
Each superposition of the components represents a specific scenario in outflow configurations.For instance, the TS (TP) morphologies can model a case when tidal dynamical ejecta has toroidal shape, while the secondary wind outflow is spherical (axially focused).Indeed, recent ab initio simulations of remnant accretion disks produce axiallyelongated outflows with spatially variable neutron richness, such that high-Y e ejecta expands along the axis, while low-Y e ejecta tends towards the equatorial plane (Fernández & Metzger 2013;Miller et al. 2019a).The bi-spherical models S1S2 (S2S1) represent the simpler scenario of isotropic two-component outflow, in which high-Y e (low-Y e ) overtakes the other component, forming an outer envelope and obscuring it (Fernández & Metzger 2016).Some studies (Wanajo et al. 2014) have suggested that the neutron fraction in a fastmoving interaction component, also known as "squeezed" dynamical ejecta (Goriely et al. 2011;Korobkin et al. 2012), can be reset to higher values by the intense neutrino irradiation from a high-temperature interface (Wanajo et al. 2014).Morphologies T1S2, S1S2 and P1S2 represent this scenario in our models.The models P1S2 and S2P1 represent polar outflow of one type, enveloped by a spherical outflow of the other type.The S2P1 scenario can be realized if initial low-Y e ejecta is fast and isotropic (Sekiguchi et al. 2016), and secondary high-Y e ejecta is slow and prolate.The slow component can be realized by neutrino-driven wind with prolate morphology, either from a transient hypermassive neutron star (Perego et al. 2014), or an accretion disk (Miller et al. 2019a).Overall, these configurations represent a reasonably exhaustive set of outflow scenarios that can happen in neutron star mergers.However, they are still limited to axisymmetry.This may be too constraining for kilonovae from neutron star-black hole mergers, where even more diversity is expected due to strong non-axisymmetric nature of the ejecta in some cases (Kyutoku et al. 2013;Perego et al. 2017;Zhu et al. 2020).

RESULTS
Here we explore the range of possible kilonova signals based on the selected suite of ejecta morphologies and two representative compositions, separately expected to produce the "blue" and "red" kilonova.For the majority of simulations, we take the mass of either one or both components to be 0.01 M , in order to put more focus on the effects of ejecta geometry.Below, we first study single-component models and then proceed to two-component models.

Single-Component Models
Figure 6 shows the evolution of angle-integrated bolometric luminosity for the one-component axisymmetric morphologies.The departure from spherical symmetry can have a large impact on the peak time and luminosity: B and T morphologies turn out to be brighter than the spherical model by a factor of two, and only need half of the time to reach their peak luminosity.The trend in brightness and peak epoch is set by the temperature of the outer radiative layer and, to a lesser degree, by the ejecta surface area.This can be inferred from Figure 7 which shows the ratio of diffusion surface area to the area of the outer boundary, as a function of time.It is clear that the surface area is not the dominant factor in deciding the peak brightness and epoch.Diffusion surface here is defined as enclosing the opaque "core" of the ejecta where photons are escaping slower than local expansion and are therefore trapped (Grossman et al. 2014).As the expansion progresses, the diffusion surface (along with the photosphere) recedes inwards until the entire bulk becomes exposed and transparent.To better understand the uncertainties in what sets the brightness of single-component models, we present a more extended discussion of these aspects in Appendix B.
Our single-component simulations assume the same mass: 0.01 M , and yet produce a wide range of peak times and luminosities, as shown in Table 2. Several groups have developed formulae that relate the peak luminosity and time of peak luminosity to the mass, expansion velocity and opacity of the ejecta (Grossman et al. 2014;Wollaeger et al. 2018).To understand how the variation in the light curve leads to errors in the mass estimates, we have inverted the pair of formulae, Eqs. ( 27)-( 28) from Wollaeger et al. (2018), deriving the mass inferred from the light curves produced in the models: where t peak and L peak are the time and bolometric luminosity at the peak, κ is the average opacity, and the resulting mass is in units of 0.01 M .
Inserting the peak parameters of our model light curves, we see that a mass which could be inferred from an observation using this formula can range over several orders of magnitude (see Table 2).The inferred mass only weakly depends on the effective opacity: for the calculations in Table 2, we assumed κ = 10 cm 2 g −1 .
We conclude that without strong observational constraints on the geometry of the ejecta, one can artificially infer larger masses from non-spherical explosions, which have higher density and correspondingly higher temperature, or higher apparent surface area (see discussion in Appendix B).If we can not constrain the composition (the late-time infra-red will be able to place some constraints), the uncertainty can be even larger.

Two-Component Models
Combining two components with different morphologies enables a much wider class of kilonovae.The range of mass estimates that can be inferred from limited information such as peak brightness and epoch now spans almost three orders of magnitude (see Table 2).The range of inferred masses can vary strongly even within the same model.For example, in the model T2P5 with fast toroidal lanthanide-rich component and slow high-Y e outflow, the mass estimate ranges between 0.002 and 0.23 M , depending on orientation.In the side view, the toroidal component with lanthanide-rich composition has a very high opacity that completely obscures bluer and brighter secondary outflow, which in turn leads to smaller inferred mass.This is the lanthanide "curtaining": high-opacity lanthanide-rich ejecta in front of more luminous outflow acts as like a "curtain", hiding the blue kilonova for a range of viewing angles (Kasen et al. 2015).At the same time, the top view leads to overestimated mass because of the brighter transient, since the secondary outflow is revealed and the toroidal component has large projected area.In Fig. 8, we plot the angle-integrated bolometric luminosity versus time for all mixed models.It shows that the impact of morphology can be a dominant indicator of the brightness, even when summed over viewing angle.In every model where the low-Y e ejecta is more extended than the high-Y e ejecta (S2S1, S2P1 and S2T1), the high lanthanide opacities block all of the optical emission, causing the peak luminosity to be an order of magnitude dimmer than the rest of the models.On the other hand, models in which the high-Y e wind component is spherical and more extended (S1S2, P1S2 and T1S2) are the brightest, but also peak too soon.
The rest of the models lie in between these two extremes.The panels in Figure 8 correspond to three different masses of the lanthanide-rich component: 0.01, 0.005 and 0.002 M , while the mass of the high-Y e component is kept fixed at 0.01 M .In every panel, the models that lie between the two extremes span a wide range of peak times (from 0.2 to 6 days) and two orders in luminosity.For all models, there is a trend of an approximate inverse correlation between peak time and luminosity: which also follows from Eq. ( 4) for fixed mass and opacity.This result shows that the effect of mixing morphologies is degenerate with the expansion velocity or effective opacity.Notice that the relative mass of the lanthanide-rich component has very little effect on the light curves.The only exception is the position of the luminosity peak for models S2P1, S2S1 and S2T1, in which this component is spherical and more extended, hiding the high-Y e outflow.Surprisingly, a lower mass of lanthanide-rich component in these models leads to a brighter and earlier peak, due to the fact that with lower mass it becomes transparent sooner, exposing the brighter high-Y e outflow.This effect is again degenerate with changing the expansion velocity or opacity, which further complicates ejecta mass estimates.
It is unlikely that measuring the light curves in a broad range of bands will be sufficient to lift these degeneracies.Fig. 9 shows the optical g-and nIR K-band light curves for all mixed morphologies and low-Y e component mass of 0.002 M .In the g-band, it is easy to tell apart models with more extended spherical lanthanide-rich ejecta (S2P1, S2S1 and S2T1, black lines): they are very strongly suppressed in this band for all orientations.These models instead produce the brightest red kilonova in the K-band.The rest of the models peak between 0.5 and 2 days in the g-band with magnitudes −15 to −16 mag, and rapidly decay by about 4 days.Models T2P5 and T2S5 present a notable exception: they radiate longest in the g-band for the "top" orientation and are strongly suppressed for the "side" view.This strong angular dependence is due to lanthanide curtaining.
Most of our mixed models do not show pronounced dependence on the viewing angle beyond that of single-component models, where variability of about one mag can be attributed to the projected area of the photosphere (Grossman et al. 2014;Darbha & Kasen 2020).This is because lanthanide curtaining is only partial: T1S1 is a typical representative of this scenario.In Figure 10, we compare its spectra with that of a single-component model T which lacks the high-Y e contribution.At early times, the spectra for model T1S1 are significantly brighter and bluer than for the model T. The impact of the high-Y e outflow also extends to late times, maintaining bluer emission relative to T. The latter behavior and the much weaker angular dependence of the blue wing of T1S1 at late time are the result of optical reprocessing by the wind.At 8 days, the viewing angle-dependent P Cygni features of Nd lines start to form in the mid-IR.However, these features only appear for the "top" orientation, as in this case the differential expansion velocity towards the observer is only about 0.1 c (more about this at the end of this Section).
In some mixed-morphology scenarios, the viewing angle dramatically affects the spectrum and observed kilonova magnitudes.This is illustrated with the model T2S5 (Fig. 10, right column).Here the lanthanide-rich compo-nent is much more extended, curtaining the high-Y e outflow for some orientations.This behavior is unlike that observed for T1S1, where both components are always visible.Moreover, the blue wing of T2S5 for the top orientations is much brighter than that of T1S1, because the lanthaniderich toroidal component additionally redirects the flux towards low-opacity polar regions (also observed in the kilonova models of Kawaguchi et al. 2018).At 8 days, pronounced P Cygni features of Nd lines start to develop at the mid-IR wing of the spectrum, resembling those of model T, but twice as broadened due to higher line-of-sight expansion velocity.This is in contrast with model T1S1 where these features are diluted by reprocessing in the enclosing high-Y e spherical component (bottom middle panel in Fig. 10).
The strong dependence on orientation for some mixed models directly projects onto broadband light curves.As can be seen in Figure 9, models T2P5 and T2S5 (shown in red), are strongly suppressed in optical bands for the "side" orientation.On the other hand, for the "top" orientation, they show the same magnitude in optical bands as the rest of the models.In the right column, the velocity of the high-Ye "wind" vw = 0.05 c is much less than for the low-Ye outflow v d = 0.2 c.The differences in the velocities for two outflows in the right column cause strong variation on the spectrum due to lanthanide curtaining or the high-Ye outflows by the more extended lanthanide-rich toroidal component.This is the case which manifests the strongest angular dependence, while for the rest of the models the angular dependence is very moderate.Black vertical arrows in the bottom row point to the P Cygni features in mid-IR generated by peaks of the wavelength-dependent opacity of Nd.
Moreover, they produce the longest-lasting blue kilonovaedue to the flux redirection mentioned above.
Because of the very high opacity of the lanthanides, a strong angular dependence persists even if the mass of the lanthanide-rich component is very small.In our models, it goes down to 0.002 M .This is shown in Figure 9, where the difference for the g-band is 4 mags between top and side orientations.In nIR bands, such as the K-band shown in the right column of Figure 9, the difference is only 1 mag.
At late times, we observe P Cygni features (Castor & Lamers 1979;Robinson 2007) forming in the mid-IR wing of the spectra of one-component models and in the spectra of models where high-Y e and low-Y e components are well-separated (T2P5 and T2S5).These features form around strong lines in the opacity and can be used to characterize both composition and morphology of the ejecta.Figure 11 overplots the opacity κ λ and the late-time spectra for singlecomponent models T and H, each with velocities of 0.1 c (denoted H1, T1) or 0.2 c (denoted H2, T2).For the toroidal morphology T, the spectrum is shown from a top view and for the hourglass H from the side.The opacity is evaluated at ρ = 10 −16 g cm −3 , T = 1000 K, and scaled by a constant factor to facilitate comparison with the spectra.The low-velocity models H1 and T1 produce narrow features due to less Doppler broadening or blending between emission peaks.The high-velocity models H2 and T2 appear to have These models are moving at an average speed of 0.1 or 0.2c.Relative to the 0.1c simulations, line emission blending is enhanced in the 0.2c simulations.Line features in the spectra are more visible in the axial (edge) view for model T (H), where the velocity towards the observer is slower.A majority of the tall lines above ∼ 4 µm are from Nd, which is setting the IR emission at late time in our models (cf.much less structured spectrum, in particular around 10 µm.These features clearly reflect spikes in the opacity profile, which in turn are produced by Nd at this density and temperature (cf.Fig. 8 in Even et al. 2019).

CONCLUSIONS
We study a set of analytically prescribed axisymmetric and spherically-symmetric morphologies representative of neutron star merger ejecta using the multidimensional, multigroup Monte Carlo code SuperNu (Wollaeger & van Rossum 2014) and detailed opacities from the LANL suite of atomic physics codes (Fontes et al. 2015(Fontes et al. , 2020)).Our five representative basic shapes are obtained from the family of Cassini ovals in axisymmetry (see Fig. 2), rescaled to a given total mass and average expansion velocity.We model the merger ejecta with single-or two-component morphologies with two spatially uniform representative compositions producing "red" and "blue" contributions of the kilonova.Specifically, we focus on two compositions: a low electron fraction (Y e ) lanthanide-rich Solar r-process residuals (Even et al. 2019), and a high-Y e composition with Y e = 0.27 ("wind 2" in Wollaeger et al. 2018).Compositions determine not only specific nuclear heating but also opacities of the material.The former is computed with the WinNET nuclear network (Winteler et al. 2012) with energy partitioning between radiation species and spatially-dependent thermalization (Barnes et al. 2016;Rosswog et al. 2017;Wollaeger et al. 2018).For the latter, we employ a new suite of atomic opacities, which includes a complete set of lanthanides and uranium (Fontes et al. 2020).
The results of our study can be summarized as follows: 1.When detailed, temperature-dependent opacities are used, the morphology affects the light curve more than the mass and velocity of the ejecta.For the same mass and expansion velocity, switching between the different morphologies considered here leads to stronger variation in the peak time and luminosity than if we fixed the morphology and varied mass (by a factor of ten) or velocity (within 0.05 − 0.3 c -see Fig. 6 and 8).
2. Because of this result, the ejecta mass can not be correctly inferred solely from the peak bolometric lumi-nosity and time.Our simulations show that the mass inferred from an effective gray opacity formula can vary by an order of magnitude within one-component morphologies, and by three orders of magnitude for mixed models.Without strong constraints on the geometry of the ejecta, artificially large or small masses can be inferred from non-spherical explosions (see Table 2).
3. Morphological variability affects the peak luminosity in a way similar to varying the ejecta velocity or effective opacity.For a family of morphologies with the same mass and expansion velocity, the peak luminosity is inversely correlated with peak times: L peak ∝ t −5/4 peak (Fig. 5).The effect of mixing morphologies is degenerate with the effect of changing the average expansion velocity or effective gray opacity.
4. Unlike the gray-opacity case, in the models with detailed opacity, the temperature at the surface is more important in defining the peak properties than the projected area (see Appendix B for detailed analysis).
5. Density-dependent thermalization of radioactive heating adds an extra boost to the bolometric luminosity, increasing it by almost a factor of two.Therefore, a proper treatment of thermalization is as important as accurate composition-dependent nuclear heating rates (see also Barnes et al. 2016;Rosswog et al. 2017;Hotokezaka & Nakar 2019).
6.It is difficult to achieve the lanthanide curtaining effect if the two components have similar expansion velocities (Kasen et al. 2015).In this case, only partial lanthanide curtaining is observed.The majority of our models are quite isotropic, except when the high-Y e component is much slower than the lanthanide-rich component (models T2S5 and T2P5 -see Fig. 10).Lanthanide curtaining is observed only in these models at low-latitude angles.
7. For the models with lanthanide curtaining, even a very small mass of 0.002 M of lanthanide-rich material is sufficient to obscure the blue kilonova, because of the exceptionally high opacity of lanthanides (Kasen et al. 2013;Tanaka & Hotokezaka 2013;Fontes et al. 2015).
8. At late epochs ( 8 days) the models in which lanthanide-rich ejecta is more extended exhibit pronounced P Cygni features, allowing the characterization of composition and line-of-sight velocity of the ejecta (Fig. 11).Our results here confirm the findings of Kawaguchi et al. (2018).These P Cygni features are unambiguously unidentified with peaks in opacity of Nd, which has been shown to dominate the opacity of lanthanide-rich material (Even et al. 2019;Fontes et al. 2020).
9. Light reprocessing: in the case when the lanthaniderich component is engulfed in a spherical high-Y e envelope, the blue kilonova becomes almost isotropic and the P Cygni features from Nd disappear (Fig. 10, middle column).
10. Light focusing by toroidal component: in the models with more extended toroidal dynamical ejecta, the light from the blue component is effectively funnelled towards the axis such that the blue kilonova appears brighter (Fig. 10, right column).
Morphological freedom creates sufficient variability to crudely fit the early spectrum of kilonova GW170817 with the limited set of models that we presented in this study.In Figure 12 we demonstrate a "fit" to GW170817 spectra at t = 1 days.The "best-matching" models appear to be T2S2 and T2P2.Clearly, tuning the velocities and masses of individual components, as well as the viewing angle, can significantly improve this result.This concept, yet again, demonstrates that a meaningful interpretation of kilonova light curves must properly account for morphological features which should be informed by numerical simulations.
We summarize the bolometric and broad band magnitudes at 1, 4, and 8 days for all models used in this study in Table 3 below.The complete suite of our models is available from the LANL CTA website. 2 "Gravitational Radiation and Electromagnetic Astrophysical Transients (GREAT)" funded by the Swedish Research council (VR) under Dnr 2016-06012.We gratefully acknowledge support from COST Action CA16104 "Gravitational waves, black holes and fundamental physics" (GWverse) and from COST Action CA16214 "The multi-messenger physics and astrophysics of neutron stars" (PHAROS).Left column: uniform gray opacities (10 and 1 cm 2 g −1 for top and bottom panels, resp.)Middle column: simulations using detailed opacities.Right column: bolometric luminosity estimate based on the temperature and the area of the outer surface (see Eq. B1).All models have mass mej = 0.01 M , median expansion velocity vej = 0.2 c and uniform analytic power-law heating (without thermalization).The labels on the right panel show area of the outer surface for the morphologies relative to the spherical model.

B. RADIATIVE STRUCTURE IN SINGLE-COMPONENT MORPHOLOGIES
One way to understand the features of the kilonova light curves is to break it up into components.In this appendix, we study how the light curve of a single-component morphology is influenced by different factors, such as area of the photosphere or the temperature of the radiative layer.We apply phenomenological analysis of the radiative structure of axisymmetric morphologies and disentangle effects of geometry, density and opacity.For all models in this section, we assume uniform specific heating given by ε(t) = 2 × 10 10 t −1.3 d erg g −1 s −1 , where t d is the time since merger, in days (Metzger et al. 2010;Korobkin et al. 2012).

B.1. Gray-opacity models
In the gray-opacity approximation, models with more mass distributed at low optical depth are expected to be brighter.This is illustrated in Figure 13 (left column).Total luminosity can be estimated then using a simple expression: where ε(t) is the specific volumetric heating rate, and m unc (t) is the "uncovered" mass, or the mass of the layer above the diffusion surface.The latter is defined as a surface at which the diffusion velocity v diff = c/τ equals the velocity differential to the edge of ejecta.Diffusion surface separates the bulk where photons are trapped from the outer envelope from where photons can escape, diffusing faster than the matter is expanding.Here, τ is the optical depth.In every (comoving) point of homologously expanding ejecta with uniform gray opacity, τ (t) ∝ t −2 (Grossman et al. 2014).This model has limited applicability: it does not take into account thermalization in thin optical region above the photosphere and thus overestimates luminosity after the peak is reached.
Figure 14 shows the fraction of uncovered mass as a function of optical depth for axisymmetric morphologies, normalized to the same mass and root mean square (RMS) expansion velocity.This profile m unc (τ ) is only dependent on the density distribution inside each morphology, and is sufficient to estimate the bolometric light curve-using expression (B1) and the fact that τ (t) ∝ t −2 .If more mass is "buried" at high optical depth, the kilonova will peak later and be dimmer.Looking at Figure 14, we can anticipate that model S will peak earliest, followed by models B and T, and the models P and H will peak last.Similar trend is expected for the peak luminosity (in decreasing order).This is indeed observed in simulations with gray opacity in Figure 13 (left column), both for opacity κ = 10 and κ = 1.

B.2. Thin-layer models
When detailed opacities are introduced, their temperature dependence complicates the light curves, shifting peak time and luminosity by a factor of a few (Fig. 13, middle column).To understand these effects, we constructed a simple thin-layer approximation which uses snapshots of radiation temperature recorded by SuperNu during our simulations.Bolometric luminosity is computed as follows: where T out is the recorded temperature at the surface of the morphology (plotted in Figure 15), and A(t) is the uniformly expanding area of this surface: A(t) = A 0 (t/t 0 ) 2 .We use c a T 4 out for the surface flux instead of the usual σ T 4 out (≡ c a T 4 out /4) because in the thin photospheric layer the radiation is free-streaming, such that intensity distribution is strongly peaked towards outward normal rather than being isotropic.In this case, the entire radiative energy in the bulk of photospheric layer is escaping.We also neglect photospheric recession and corresponging decrease in the emitting area.Figure 7 showing the fractional area of diffusion surface (discussed in detail below) proves this to be a reasonable assumption, as for non-spherical morphologies the area of diffusion surface remains at least 0.8 of the area of the outer contour up to about 2 days.For the spherical morphology, this assumption is not very accurate and the "thin-layer" model is expected to overestimate the luminosity.
The resulting evolution of bolometric luminosity L thin (t) is presented in the right column of Figure 13.Comparing it to the middle column on the same plot we see that our model correctly reproduces the features of the bolometric light curves, but overestimates luminosities by about 20 − 40%.It also captures the order of peak times for different morphologies: B peaking first and H peaking the last.The same is only partly true for the order in peak brightness as this approximation overestimates it for morphology H. Overall, our numerical experiment demonstrates that the features in bolometric light curve are primarily dictated by the behaviour of the surface temperature, and less so by the emitting area.In particular, spherical morphology S stands out by being the "coldest" and thus having the least lumonous peak despite its having the largest surface area for the same mass and mean expansion velocity.
Next, we observe that the surface temperature ranking for different morphologies (S < R < T ≈ P < H, Fig. 15) can be traced back to their ranking in density (Fig. 3).Naturally, models with higher density produce more radioactive heat per volume and are expected to be hotter.Moreover, if we compare bolometric luminosity of the low-Y e composition computed with thermalization taken into account (Fig. 8), it boosts the higher-density models even further, as more radioactive energy is thermalized in denser regions.Overall, density-dependent thermalization adds an extra boost by about factor of 2 in luminsity (cf.Fig. 8 and Fig. 13) The thin-layer model qualitatively reproduces the trends of the full RT model because it uses the surface temperatures from the simulations, which proves that the outbound MC flux is consistent with the temperature evolution on the surface.Below we attempt to further validate the "thinness" assumption by explicitly calculating an approximate position of the diffusion surface from simulation data.

B.3. Estimating location of the diffusion surface and photosphere
Figure 16 shows the temperature color maps in our basic morphologies at different times, with overplotted estimated contours of the photosphere (dashed lines) and diffusion surface (solid lines).To locate the photosphere, we used Rosseland mean opacities and integrated the optical depth dτ = κ R (ρ, T )ρd inwards from the outer edge of the expansion, adjusting the path of integration so that it always follows the local density gradient.This results in a family of contours which determine the optical depth globally for every point as the optical depth minimized over all possible escape routes.
The photosphere is then given by the τ = 2/3 contour.It turns out, however, that Rosseland mean significantly underestimates true effective opacity and places the photosphere too deep.Simple numerical integration of the radiative energy density in the layer above the photosphere computed in this manner gives numbers which greatly exceed the observed luminosity output of the kilonova.
A better way to pinpoint location of the radiative layer is given by the diffusion surface.We define the diffusion surface as enclosing the opaque "core" of the ejecta where photons are escaping slower than local expansion and are therefore trapped (Grossman et al. 2014).To simplify the analysis, we make further approximations and compute the diffusion surface as given by an optical depth τ ds (t) such that the integral of radiative energy E rad above this contour (τ < τ ds ) is equal to the total bolometric luminosity: τ <τds(t) E rad dV = L bol (t). (B3) This approximation ignores the anisotropu of radiation flux due to asphericity of our models, but it is nevertheless sufficient to estimate the validity of the thin-layer approximation.Indeed, Figure 16 allows to conclude that for all morphologies except S, up to about 2 days the diffusion surface is not only very close to the outer edge of the expansion but also that the temperature in the radiative layer above it does not change significantly.So taking the surface temperature as a proxy for thin-layer model is justified.On the other hand, contours of optical depth computed using Rosseland mean fail to capture the photosphere: they instead place it inside the diffusion surface, which is unphysical for a typical kilonova scenario.
For any non-gray opacity model, location of the photosphere strongly depends on the wavelength.This point is clearly illustrated on Figure 17 (plots of this type were first introduced in Fontes et al. 2020)).Each plot represents a color map of the wavelength-dependent optical depth experienced by a photon entering the ejecta and moving towards the density maximum on a straight line orthogonal to the longer axis of the morphology.For non-spherical morphologies photons traverse velocity differential of 0.1 c and are blue-shifted in the rest frame of the expansion, therefore experiencing progressively more blue-shifted opacities.This is easy to see on the plots for spherical morphology (middle row) where the velocity differential is 0.4 c and there is a noticeable upward "tilt" of the plot features towards the left edge.
The three-color heat map was selected to show three different regions in opatical depth: low (blue), medium and comparable to the diffusion speed (black) and high (red).White matte regions on the plots are added to highlight the parts of the spectra where most of the photons would be emitted if spectra were blackbody.Transparency of the matte layer is proportional to the Planck function with the local radiation temperature.Solid and dashed vertical lines on each plot show the location of the diffusion surface and photosphere (computed with Rosseland mean).
The common feature on all plots is that the solid line crosses the optical depth map predominantly in the regions where the optical depth is "black", i.e. ∼ 1 − 10.In other words, the surface computed with expression (B3) is located at optical depth of about τ ds ≈ 10, which is where diffusion surface should be for an outflow with expansion velocity of c/10.This provides additional argument in support of using line-binned opacities instead of expansion opacities for non-spherical models (see Fontes et al. 2020, for detailed exposition of the arument).

Figure 2 .
Figure2.Maps of density of the basic building-block morphologies used in this work, shown in the upper polar half-plane (the R − z plane, with R being the cylindrical radius): "hourglass" (H), "peanut" (P), "sphere" (S), "biconcave" (B) and "torus" (T).The contours are equi-spaced by 0.5 dex in log space, down from the density maximum.

Figure 3 .
Figure 3. Density profiles of the basic building-block morphologies along the longer axis.Morphologies are constrained to have the same mass and average (RMS) expansion velocity.In gray: spherically averaged density profiles of the dynamical ejecta from neutron star mergers (Rosswog et al. 2014).

Figure 4 .
Figure 4. Density versus velocity of neutrino-driven accretion disk wind from an ab-initio radiation GRMHD simulation(Miller et al.  2019b,a).Outflow from the disk is highly aspherical.

Figure 5 .Figure 6 .Figure 7 .
Figure5.Schematics of the adopted two-component combinations of morphologies in the meridional plane.The first letter represents abbreviation of the neutron-rich (low-Ye) outflow morphology, while the second letter stands for the higher-Ye outflow shape.The low-Ye (high-Ye) component is represented by the red (blue) contour.The first (second) digit encodes mean expansion velocity of the low-Ye (high-Ye) component: 1 = 0.1 c, 2 = 0.2 c, and "5" designates 0.05 c.This is reflected in the relative sizes of the components.

Figure 8 .
Figure 8. Bolometric luminosity as a function of time for twelve mixed morphologies, and three different masses of the low-Ye (neutron-rich) component.The top plot shows ∝ t −5/4 trend.

Figure 9 .
Figure9.Broadband light curves (g-and K-bands) for the twelve mixed morphologies with the lanthanide-rich (low-Ye) component mass 0.002 M .Top row shows the on-axis, or "top" orientation and bottom panels show the "side" orientation towards the observer.

Figure 10 .
Figure10.Kilonova spectra at 1, 4 and 8 days for a single-component model "T" (left) and two-component "TS" morphologies (middle and right), for the full range of viewing angles.The middle and right components differ by the expansion velocities: in the middle column, both components have the same expansion velocity v d = vw = 0.1 c.In the right column, the velocity of the high-Ye "wind" vw = 0.05 c is much less than for the low-Ye outflow v d = 0.2 c.The differences in the velocities for two outflows in the right column cause strong variation on the spectrum due to lanthanide curtaining or the high-Ye outflows by the more extended lanthanide-rich toroidal component.This is the case which manifests the strongest angular dependence, while for the rest of the models the angular dependence is very moderate.Black vertical arrows in the bottom row point to the P Cygni features in mid-IR generated by peaks of the wavelength-dependent opacity of Nd.

Figure 11 .
Figure11.Axial-view flux for model T and edge-view flux for model H, along with scaled opacity, versus wavelength.These models are moving at an average speed of 0.1 or 0.2c.Relative to the 0.1c simulations, line emission blending is enhanced in the 0.2c simulations.Line features in the spectra are more visible in the axial (edge) view for model T (H), where the velocity towards the observer is slower.A majority of the tall lines above ∼ 4 µm are from Nd, which is setting the IR emission at late time in our models (cf.Fig. 8 in Even et al. 2019).

Fig. 8 in
Figure11.Axial-view flux for model T and edge-view flux for model H, along with scaled opacity, versus wavelength.These models are moving at an average speed of 0.1 or 0.2c.Relative to the 0.1c simulations, line emission blending is enhanced in the 0.2c simulations.Line features in the spectra are more visible in the axial (edge) view for model T (H), where the velocity towards the observer is slower.A majority of the tall lines above ∼ 4 µm are from Nd, which is setting the IR emission at late time in our models (cf.Fig. 8 in Even et al. 2019).

Figure 12 .
Figure 12.Model spectra for various geometries at the epoch t = 1 day in comparison with the spectrum measured for the kilonova AT 2017gtfo accompanying GW170817 (instrument: VLT/Xshooter, Pian et al. 2017).All models have equal mass of the neutron-rich and neutron-poor component, m d = mw = 0.01 M and are shown with on-axis view.Wavelength ranges with poor data quality for the observed spectra are marked in gray (see Pian et al. 2017, for detals).

Figure 13 .
Figure13.Bolometric luminosity as a function of time for the single-component morphologies.The top and bottom rows represent different compositions: low-Ye lanthanide-rich solar r-process residuals (left) vs high-Ye composition (right).Left column: uniform gray opacities (10 and 1 cm 2 g −1 for top and bottom panels, resp.)Middle column: simulations using detailed opacities.Right column: bolometric luminosity estimate based on the temperature and the area of the outer surface (see Eq. B1).All models have mass mej = 0.01 M , median expansion velocity vej = 0.2 c and uniform analytic power-law heating (without thermalization).The labels on the right panel show area of the outer surface for the morphologies relative to the spherical model.

Figure 14 .Figure 15 .
Figure14.Fraction of mass uncovered above certain optical depth for models H, P, S, B, T. Gray opacity κ = 1 cm 2 g −1 was used.The models are constrained to have the same total mass and average (RMS) expansion velocity.

Figure 16 .
Figure 16.Temperature maps for axisymmetric single-component morphologies with solar r-process residuals for four different epochs.Black contours on each model mark estimated location of the diffusion surface (solid line) and the photosphere (dashed line).

Figure 17 .
Figure 17.Wavelength-dependent optical depth for single-component morphologies with solar r-process residuals for four different epochs.Solid lines show location of the diffusion surface.

Table 2 .
Ejecta masses inferred using equation 4, for single-component and mixed morphologies.The second column shows the mass (masses) of the single component (low-Ye component + high-Ye component) in percents of the Solar mass (0.01 M ).