Towards inferring the geometry of kilonovae

Recent analysis of the kilonova, AT2017gfo, has indicated that this event was highly spherical. This may challenge hydrodynamics simulations of binary neutron star mergers, which usually predict a range of asymmetries


INTRODUCTION
The kilonova AT2017gfo, which was coincident with the gravitational wave signal GW170817 (Abbott et al. 2017), has provided us with a wealth of observations (e.g., Smartt et al. 2017;Pian et al. 2017;Villar et al. 2017;Coulter et al. 2017).Recently, Sneppen et al. (2023b) presented evidence that AT2017gfo was highly spherical, which is surprising given that binary neutron star merger simulations show strong asymmetries (Bauswein et al. 2013;Sekiguchi et al. 2015;Bovard et al. 2017;Radice et al. 2018;Foucart et al. 2023;Combi & Siegel 2023).Sneppen et al. (2023b) inferred the expansion velocity of the ejecta by analysing the most prominent feature in the spectra of AT2017gfo, which has been suggested to be a Sr ii P-Cygni feature (Watson et al. 2019;Gillanders et al. 2022;Domoto et al. 2021Domoto et al. , 2022, although see Tarumi et al. 2023 for discussion of He i as an alternative explanation for this feature).Such line profile analysis is predominantly sensitive to the line of sight velocity component.Sneppen et al. (2023b) also inferred a photospheric radius using the expanding photosphere ★ E-mail: c.collins@gsi.demethod (EPM; Baade 1926;Kirshner & Kwan 1974;Eastman et al. 1996) and, assuming homologous expansion of the ejecta, converted this to an expansion velocity of the photosphere.This method is primarily sensitive to the expansion velocity perpendicular to the line of sight.They found remarkable consistency between the velocities obtained via these two methods across two observed epochs, suggesting consistency between line of sight and perpendicular expansion velocities, which is indicative of a near-spherical explosion.To quantify the inferred sphericity, Sneppen et al. (2023b) defined a zero-centred asymmetry index, Υ =  ⊥ − ∥  ⊥ + ∥ , where  ∥ is the expansion velocity of the photosphere primarily along the line of sight (obtained from the P-Cygni profile analysis) and  ⊥ is the expansion velocity of the photosphere primarily in the direction perpendicular to the line of sight (obtained from the EPM), and found Υ = 0.00 ± 0.02 1 .In addition to this analysis, they identify that the shape of the P-Cygni feature is best matched by assuming a spherical photosphere across multiple epochs (see Sneppen et al. 2023b).Shingles et al. (2023) presented a three-dimensional radiative transfer calculation using ejecta from a binary neutron star merger simulation coupled to a nucleosynthesis network.In the polar directions, the spectra resembled the observations of AT2017gfo remarkably well, considering the model was not tuned in any way to match this event, although at earlier times than those observed (see Shingles et al. 2023 for discussion).The synthetic observables showed variations with both polar and azimuthal viewing angles, and therefore are not isotropic.Here, we analyse this simulation to determine whether we can infer information about the underlying symmetry of the ejecta from the synthetic observables.
The aim of this paper is to determine whether the geometry and overall degree of sphericity of the ejecta is encoded in the observables as suggested by Sneppen et al. (2023b).To this end, we consider the synthetic observables predicted by the 3D radiative transfer calculation by Shingles et al. (2023) from an asymmetric merger simulation to analyse the velocities that would be derived from these observables following the method laid out by Sneppen et al. (2023b), and investigate whether the apparent consistency of  ∥ and  ⊥ is a fundamental challenge to the explosion model.Determining whether the geometry of simulations is compatible with observations will further our understanding of binary neutron star mergers, including the high density Equation of State, the dynamics of matter ejection, including the role of neutrinos, and the underlying rapid neutron capture (r-process) nucleosynthesis.
The EPM is typically used to measure luminosity distances, D L , to supernovae, and can be used to obtain distance measurements independently of the cosmological distance ladder, which can lead to independent measurements of the Hubble constant (de Jaeger et al. 2017;Gall et al. 2018;Sneppen et al. 2023a,b).Sneppen et al. (2023b) applied the EPM to measure the distance to AT2017gfo and found D L to be consistent with previous distance measurements.In addition to inferring the geometry of our simulations, we aim to test how accurately the luminosity distance to our synthetic spectra can be measured using the EPM.

Merger ejecta symmetry
We consider the 3D binary neutron star merger simulation of equal mass 1.35 M ⊙ neutron stars used by Shingles et al. (2023) and Collins et al. (2023).This merger simulation was carried out with a 3D general relativistic smooth particle hydrodynamics (SPH) code (Oechslin et al. 2002;Bauswein et al. 2013), which adopts the conformal flatness condition on the spatial metric (Isenberg & Nester 1980;Wilson et al. 1996) and includes an advanced neutrino leakage treatment (Ardevol-Pulpillo et al. 2019).We employ the SFHo equation of state (Steiner et al. 2013).As discussed by Collins et al. (2023), the ejecta from this merger simulation show asymmetries both in the density structure and in the distribution of   (see figures 1 and 2 in Collins et al. 2023).A higher   is found near the poles, primarily due to the inclusion of a neutrino treatment, and a lower   is found near the equator.Nuclear network calculations were carried out for each SPH particle trajectory (see Collins et al. 2023 for details of the calculation) to model the spatial distribution and time dependence of the early heating and elemental composition.The evolution of the composition and energy released due to radioactive decays is then followed by artis (see Shingles et al. 2023 for details).The mass ejected into equal solid-angle bins, spaced by polar angle, for this model is shown in Figure 1.The simulation predicts more mass per solid angle near the equator compared to near the poles.There is a mild equatorial (i.e., north vs. south) asymmetry, visible in Figure 1, which to some degree may be physical (as a result of stochastic fluctuations and hydrodynamic instabilities in the matter flow), but may also be amplified by numerical effects (e.g.particle noise).We discuss these aspects in more detail in Appendix A. We note that the merger simulation provides only dynamical ejecta (evolved until 20 ms after the merger).The merger simulation shows some variation in asymmetry with time, and it terminates before the matter ejection ceases and the ejecta configuration reaches its final state.We thus do not expect an overall match with AT2017gfo.The impact of the continued mass ejection (secular ejecta) on our simulation results should be investigated in future.However, given that this dynamical ejecta model is able to produce spectra comparable to those observed for AT2017gfo (see Shingles et al. 2023 and discussion in Section 3.1) we use this ejecta model as the basis for our analysis.
Also shown is the mass of material ejected in the velocity range 0.15c <  < 0.45c, M  , which is approximately in the line forming region during the first day of the kilonova (Section 3.3).Note that 44 per cent of the ejecta mass is below 0.15c.Figure 2 shows how the ejecta mass is distributed with polar angle and radial velocity (see figure 1 in Collins et al. 2023 for a 3D rendering of the ejecta structure).The ejecta exhibit asymmetry both in the total mass ejected into different directions, and in the velocity at which the bulk of the mass is ejected.
To quantify the asymmetry of the ejecta in the simulation, we Table 1.Mass ejected within equal solid-angle bins, defined by polar angle, at the poles and around the equator.We define  ej to be the range in polar-angle of each bin in the ejecta, which each have a width of cos(  ) = 0.2.Included is the total mass, M, within the solid-angle and the mass, M  , within the velocity range (assuming homologous expansion) 0.15c <  < 0.45c, which is approximately within the spectral line-forming region during the first day of the kilonova.This is also shown for the mass of Sr, M Sr , and the mass of Sr within the velocity range, M Sr  .define a zero-centred asymmetry parameter2 for the ejected mass, where M eq is the mass within a solid-angle at the equator and M pole is the mass within an equal solid-angle at the pole.The masses within equal solid-angles (plotted in Figures 1 and 2) are listed in Table 1 for solid-angles near the poles and equator.The values of Υ M are listed in Table 2.The values of Υ M indicate a moderate level of asymmetry in the mass per solid-angle ejected near the poles compared to the equator, to quantify the sphericity in the inferred velocities of the photosphere.Similarly, we use zero-centred asymmetry parameters to quantify sphericity, however the measurements on the sphericity of the mass and luminosity are not directly comparable to the quantity inferred by Sneppen et al. (2023b).
however, Υ M, (considering only the mass ejected within the velocity range 0.15c <  < 0.45c) indicates that the mass approximately within the line forming region within the first day of the kilonova is more symmetric, which may lead to a higher level of symmetry in the synthetic observables than indicated by Υ M .
The abundances of synthesised elements in the ejecta also show asymmetries.We show the distribution of Sr synthesised in the ejecta in Figure 2, which is predominantly responsible for the strongest feature in our simulated spectra (Shingles et al. 2023).Note that the feature is also composed of contributions from Y and Zr, which show a similar distribution in the ejecta to Sr.We show the total mass of these and other representative elements ejected by polar angle in Figure B1 in Appendix B. Higher masses of Sr, Y and Zr are synthesised near the poles than near the equator (due to the higher   near the poles, see figure 2 of Collins et al. 2023).Sneppen et al. (2023b) suggest that the line shape of the spectral feature in AT2017gfo indicates a near spherical distribution of Sr, which is not shown by the distribution of Sr synthesised in our model.Using the asymmetry parameter, Υ Sr M , listed in Table 2, the mass distribution of Sr shows a moderate level of asymmetry.Within the velocity range 0.15c <  < 0.45c, the distribution of Sr with polar angle is very asymmetric (e.g., Υ Sr M, = −0.71).

Radiative transfer simulation
Given the asymmetry of the ejecta model, we aim to test whether this level of asymmetry can be consistent with the observational constraints inferred by Sneppen et al. (2023b) for AT2017gfo, which appear to favour near spherical symmetry.For this, we analyse the 3D radiative transfer kilonova simulation, 3D AD2, carried out by Shingles et al. (2023) to identify whether, if subjected to a similar analysis as data from actual observations, it can match (or be ruled out by) constraints on its apparent sphericity, such as the Υ parameter used by Sneppen et al. (2023b).This simulation was carried out using the multi-dimensional, time-dependent Monte Carlo radiative transfer code artis (Sim 2007;Kromer et al. 2010;Shingles et al. 2020Shingles et al. , 2023, based on the methods of Lucy 2002Lucy , 2003Lucy , 2005)).We note that in artis there is no photosphere defined in the simulation.We can therefore infer the apparent photosphere from the synthetic observables (using the same methods as for observations) without the constraint of a photospheric boundary being imposed in the simulation.To reduce Monte Carlo noise in the orientation-dependent synthetic observables, escaping packets of radiation are binned into ten equal solid-angle bins, defined by polar angle.

Synthetic observables
The synthetic light curves and spectra for this simulation are presented by Shingles et al. (2023).The simulated light curves were fainter than observed for AT2017gfo, likely due to the lower ejecta mass in the model (0.005 M ⊙ ) than that inferred for AT2017gfo (0.04 ± 0.01 M ⊙ was estimated by Smartt et al. 2017).The spectral evolution in the polar direction was similar to that observed in AT2017gfo, showing a prominent feature primarily shaped by Sr ii, however, the simulated spectra show a more rapid evolution than observed in AT2017gfo.In this section, we focus on the observer orientation dependence and the level of isotropy shown by the synthetic observables.

Observer orientation variation in bolometric luminosity
As a first quantification of the degree to which the simulation predicts observer orientation dependencies we compare the bolometric light curves in different directions.As presented by Shingles et al. (2023), the light curves in the lines of sight at the poles are brighter than those at the equator, due to lower densities and lower   near the poles (see Collins et al. 2023).In Figure 3 we plot the angle-dependent bolometric light curves as isotropic-equivalent luminosities (i.e., from the simulation we record the energy emitted per second into each solid angle bin to obtain light curves in erg s −1 sr −1 for each orientation; we then scale these to an equivalent isotropic luminosity by multiplication by the full-sphere solid angle of 4•sr).Note that the solid-angle bins around the poles encompass the inferred observer angle of AT2017gfo (between 19 • and 25 • ; Mooley et al. 2022).
We define a zero-centred asymmetry index (similar to Sneppen et al. 2023b, but not directly comparable), where L pole is the luminosity at either pole and L eq is the luminosity near the equator.The luminosity emitted into the solid-angle bins above (0 < cos() < 0.2) and below (−0.2 < cos() < 0) the equator is almost identical -see Figure 3.The maximum deviation of Υ bol between the poles and the equator is Υ bol ≈ −0.18 (Figure 3).This clearly demonstrates that the synthetic observables for this ejecta model are not isotropic, however, according to this simple asymmetry metric, the synthetic observables are not as asymmetric as the mass distribution of the ejecta.We discuss this further in Section 3.1.2.

Observer orientation variation in band-limited light curves
We show the  light curves for this simulation in Figure 4. Similarly to the bolometric light curves, the band-limited light curves are not isotropic, but also do not show a very strong observer orientation dependence.In all directions, the light curves initially peak in the bluer bands, and then become brighter in the redder bands, showing similar behaviour to the observed blue to red colour evolution of AT2017gfo.This was also found for this ejecta model in an approximate sense by Collins et al. (2023).Despite having a composition with a higher lanthanide fraction at the equator than at the poles, we do not predict a significantly redder colour at the equator than the poles (see Figure B2 in Appendix B).Although, we note that the atomic data considered in this simulation does not include actinides (see Shingles et al. 2023 for details of the atomic data), and therefore the opacity may be underestimated.This suggests that in a 3D simulation, a high lanthanide fraction at the equator does not necessarily lead to a significantly redder spectral energy distribution (SED) viewing from the equator than viewing from the poles.
The reason for this behaviour can be seen in Figure 5, where we show the location in the ejecta (in velocity space, assuming homolo-  gous expansion) where radiation last interacted before being emitted towards an observer viewing from a polar or an equatorial direction.
The radiation viewed from a given direction has been emitted from a broad range of ejecta, both parallel and perpendicular to the line of sight.Viewed from an equatorial direction, radiation is emitted from high opacity, lanthanide rich regions of the ejecta near the equator, but in addition to this, radiation is also emitted towards an observer at the equator from lower opacity regions (with lower lanthanide fractions) of the ejecta near the poles.We therefore find that the asymmetry of the ejecta (in mass distribution and the variation in   ) does not strongly influence the anisotropy of the light curves for varying observer orientations, since an observer does not view radiation emitted from only one region of the ejecta.

Observer orientation variation in spectra
As discussed by Shingles et al. (2023), the model predicts that the spectra would not appear the same to an observer viewing towards the poles as to an observer near the equator.We show the viewingangle dependent spectra in Figure 6.The model spectra show phases that resemble the observed spectra of AT2017gfo in the direction of the poles (see Shingles et al. 2023), however, at earlier times than those observed.The more rapid evolution is likely due to the lower mass of our model compared to the inferred mass of AT2017gfo (see discussion by Shingles et al. 2023).The spectra in the directions of the poles show a feature resembling a P-Cygni profile, which in the spectra of AT2017gfo has been suggested to be Sr ii.As indicated by the band-limited light curves, the SEDs at the poles and equator peak at similar wavelengths, however, at the equator the spectra are relatively featureless (overlapping emission and absorption as well as Doppler broadening leads to no clear signatures of individual species, see Shingles et al. 2023) and fainter than the directions near the poles (Figure 6).Therefore, our model spectra are significantly dependent on observer orientation.We now examine whether determination of the asymmetry parameter, Υ, for the inferred photospheric velocities from these synthetic spectra can yield results consistent with observational constraints.

Inferring photospheric velocities
In this section we apply the same method as Sneppen et al. (2023b) to our simulated spectra to determine the level of symmetry that would be inferred.We use our model spectra at 0.4 days and 0.6 days, since these resemble the spectra of AT2017gfo at 1.43 days and 2.43 days, which were analysed by Sneppen et al. (2023b).We refer to the model spectrum at 0.4 days as epoch 1 and at 0.6 days as epoch 2.

Photospheric velocity from P-Cygni feature
The most prominent feature in the spectra of AT2017gfo has been suggested to be a Sr ii P-Cygni feature.Our simulated spectra show a similar feature (Figure 6), as discussed by Shingles et al. (2023).We infer  ∥ from the simulated feature, assuming it can be modelled as a simple P-Cygni feature dominated by Sr ii, however, we note that in the simulation the feature is actually a blend of features, predominantly Sr ii, Y ii and Zr ii. Figure 7 shows the P-Cygni profiles used to infer  ∥ .The P-Cygni profiles were generated using a line profile calculator 3 based on the Elementary Supernova Model of Jeffery & Branch (1990).We assume the Sr ii triplet lines to be of similar strength, and that one line (which we choose to be the midwavelength line at 10327.31Å) is representative of the triplet.The P-Cygni profile is characterised by the rest wavelength, the optical depth and the velocity of the ejecta.The velocities measured from the P-Cygni feature are listed in Table 3.Since the spectra in the equatorial lines of sight do not show any clear features, a velocity cannot be obtained in the same way from these synthetic spectra.

Photospheric velocity from EPM
Following Sneppen et al. (2023b), we infer  ⊥ using the expanding photosphere method (EPM).Assuming the emitting region to be a sphere of radius  ph , emitting as a blackbody (,  ′ ), at wavelength, , where  ′ is the inferred temperature in the co-moving frame of the ejecta, the inferred luminosity,  BB  , is given by where the blackbody flux emitted in the co-moving frame of the ejecta has been transformed into the rest frame of the observer (see Sneppen et al. 2023b;Sneppen 2023).The EPM assumes that  BB  can be equated to the pseudo-bolometric luminosity, where   is the flux and  L is the luminosity distance.The EPM assumes homologous expansion, where  ph is the velocity of the photosphere and  is the time since explosion.
Using the EPM to infer  ph , we obtain  ⊥ from the inferred blackbody temperature (by matching blackbody distributions to the synthetic spectra -see Figure 7) and the luminosity obtained from the simulation.The velocities inferred are listed in Table 3.We focus on fitting the blackbody distributions to the UV and optical, rather than to the IR since the simulated spectra appear to be well described by a blackbody at bluer wavelengths.However, compared to this blackbody fitting the spectra show a flux deficit in the IR (unlike AT2017gfo).If a blackbody is fit to the IR flux, the peak of the blackbody is much too blue compared to the synthetic spectra.
At epoch 1, the values of  ⊥ obtained from the EPM are similar to the velocities,  ∥ , inferred from the simulated spectral feature (e.g., 0.34c compared to 0.33c at the positive pole).At epoch 2, however, the velocities,  ⊥ , inferred using the EPM are much higher than those inferred from the spectral feature,  ∥ , (e.g., 0.38c compared to 0.26c at the positive pole).Additionally, the inferred photospheric velocity from the EPM increases from epoch 1 to epoch 2, in contradiction to the velocities from the spectral feature, the ejecta velocity radiation was emitted from (see Section 3.3), and the expectation that the photosphere would most likely recede with time.At our simulated epoch 2, the spectrum is not well matched by a blackbody distribution at wavelengths redder than ∼ 10000 Å (unlike AT2017gfo).This likely explains why the EPM does not produce reasonable values of the photospheric velocity for our simulation at epoch 2. To an observer, it would be apparent that the spectra are not well represented by a blackbody at this time, and that the photospheric velocity inferred from the EPM is not realistic.Dessart & Hillier (2005) noted that the EPM is best used at early times when the spectrum is closest to a blackbody.It is possible that the synthetic spectra show poorer agreement with a blackbody compared to AT2017gfo because the atomic data is not complete, for example, we do not include actinides.Therefore opacity could be missing from the simulation.Additionally, we only consider dynamical ejecta, which may also be responsible for the larger deviation from a blackbody.

Geometry implied by 𝑣 ∥ and 𝑣 ⊥
Following Sneppen et al. (2023b), we use the asymmetry index, Table 3. Quantities for inferring the geometry of our simulation for an observer viewing the simulation from a direction,  obs , where we give the range in polar-angle of  obs .The parameters listed include the following.Temperature in the co-moving frame of the ejecta, T ′ , inferred from fitting a blackbody (, T ′ ), transformed into the rest frame of the observer, to the simulated spectra (see Figure 7).Inferred photospheric velocity,  ∥ , from the P-Cygni feature, and the luminosity that would be inferred for these photospheric velocities using the expanding photosphere method ( BB  ), which can be compared to the simulated luminosity  bol .Perpendicular velocity,  ⊥ , inferred from the simulated luminosity ( bol ) using the EPM.Asymmetry index, Υ ,ph , inferred from  ∥ and  ⊥ .Note that by 0.6 days the EPM does not predict reliable velocities for our simulation, which is the reason for the higher level of inferred asymmetry at this time.observer directions ( obs , listed in each panel).Dashed lines show the blackbody distribution (in the co-moving frame of the ejecta transformed into the rest frame of the observer, scaled to match the brightness of the synthetic spectra) best matching the spectra, where the ejecta temperature, T ′ , and the parallel photospheric velocity,  ∥ are listed in each panel.The blackbody distribution is modified to include a P-Cygni profile for Sr ii, considering spherical, prolate or oblate ejecta.For the spectra at the equator, indicative photospheric velocities are chosen for the relativistic correction to the blackbody distribution, since photospheric velocity cannot be inferred from a spectral feature at the equator.

Epoch Time
Table 4. Luminosity distance, D L , estimated using the EPM in the direction  obs .We set the distance to our simulated spectra as 1 Mpc.The photospheric velocity considered is that inferred from the P-Cygni profile ( to quantify the degree of sphericity implied by our synthetic spectra, which we show in Table 3.At epoch 1, a high level of symmetry is inferred from both polar directions.Indeed, the inferred symmetry at the positive pole is within the uncertainty of the sphericity inferred by Sneppen et al. (2023b) for AT2017gfo (Υ = 0.00 ± 0.02).This demonstrates that aspherical ejecta can lead to directions that can appear near-symmetrical to an observer, as quantified by the Υ ,ph measurement.
At epoch 2, however, the inferred values of Υ ,ph are much higher.As discussed in Section 3.2.2, this is likely because the simulated spectra are no longer well matched by a blackbody and the inferred perpendicular velocities are much higher than expected.This was not the case for AT2017gfo, during the epochs analysed by Sneppen et al. (2023b), where the spectra continued to resemble a blackbody until later times than in our simulation.

Distance estimate
The EPM can be used to infer the luminosity D L , by equating Equations 3 and 4. We test how accurately the distance can be inferred for our simulated spectra.We set the distance to the simulated spectra as 1 Mpc, and the inferred distances using the EPM are listed in Table 4.For this calculation, we assume a spherical photosphere ( ∥ =  ⊥ ) and use the velocity,  ∥ , inferred from the P-Cygni profile (Table 3) to determine R ph .We use the temperatures inferred from the blackbody distributions matching the simulated spectra (Table 3), and the flux from the simulated spectra.
At epoch 1, the distance can be inferred to a good degree of accuracy (within 4-7 per cent), however, at epoch 2 when the spectra are no longer well matched by a blackbody the distance estimate is much more uncertain (>25 per cent error) and underestimated.For each epoch, the direction with the more spherical value of Υ ,ph (Table 3) gives a closer estimate of distance to the actual distance, possibly indicating that testing the sphericity implied by  ∥ and  ⊥ could be a test for how good a distance estimate can be obtained, but this should be investigated for more models in future.
Both the P-Cygni profile analysis and the EPM assume that the photosphere is a sharp boundary at a single velocity.We discuss in Section 3.3 that in our simulation radiation escapes from a broad range of ejecta velocities, which is likely also the case for AT2017gfo, indicating that the photosphere is not a sharp boundary.However, even with this simple assumption an accurate distance estimate can be obtained while the spectra resemble a blackbody.

Last-interaction velocities
To give an indication of how well the asymmetry parameter, Υ ,ph , from the inferred photospheric velocities represents the symmetry of radiation leaving the simulation, we extract from the artis calculation the ejecta radial velocity ( i ) at which each Monte Carlo packet last interacted with the ejecta before escaping towards an observer.artis does not impose any photospheric boundary condition on the simulation, but the distribution of  i provides an indication of where radiation-matter interactions are occurring and thus the location of the spectrum forming region (see Figure 5).The mean, energy weighted radial velocity at which radiation packets last interacted with the ejecta ( vi ) is listed in Table 5 for radiation that has escaped the ejecta, travelling towards an observer at a polar or equatorial orientation (i.e., packets of radiation that have escaped into the given solid-angle bin,  obs ) at epochs 1 and 2. The range of  i of packets that have escaped in a given direction can be seen in Figures 5 and 8 for epoch 1, giving an indication of the extent of the spectrum forming region.In particular, at the equator the radiation escapes from an extremely broad range of velocities, as can be seen by the more flat topped distribution of  i in Figure 8 at the equator, compared to at the pole which is more sharply peaked around 0.3c.
We also show in Figures 8 and 9 the range of ejecta velocities where the last absorption process underwent by a packet was with the Sr ii triplet ( Sr i ), with wavelengths of 10036.65,10327.31and 10914.87Å, immediately before being re-emitted towards an observer in a given direction.The mean, energy weighted radial velocities where radiation was last absorbed by the Sr ii triplet, vSr i , are listed in Table 5.
As discussed by Shingles et al. (2023), Sr ii is not the only species responsible for shaping the predicted feature.The Sr ii triplet absorption is, however, at the velocities required to match the feature.The values of vSr i are higher than the values of vi (e.g., 0.42c compared to 0.35c, see Table 5).This is broadly consistent with the general principles of the simple P-Cygni model adopted to fit this feature, i.e., photon interactions in the Sr line are generally occurring in a spatially extended line-forming region that extends above the region in which the pseudo-continuum forms.
As expected, the values of vi and vSr i decrease from epoch 1 to 2, indicating that the spectrum-forming region recedes from epoch 1 to 2 in the simulation, verifying that the velocities inferred at epoch 2 using the EPM are too high.

Symmetry of last-interaction velocities
To quantify the level of symmetry shown by the mean last-interaction velocities, we compare vi of radiation that escaped towards an observer at a polar orientation ( vi,pole ) to vi of radiation escaping towards an observer at an equatorial orientation ( vi,eq ), using which is listed in Table 6.At both epochs, Υ v,i indicates a high degree of symmetry, with relatively low values for Υ v,i .However, Υ v,i indicates a slightly lower level of symmetry than determined in Section 3.2.3 from the inferred photospheric velocities, Υ ,ph , (at epoch 1), where Υ v,i > Υ ,ph .Υ v,i indicates that the level of symmetry increases from epoch 1 to 2 in this model.
In agreement with Υ ,ph inferred from spectra at either pole at  7) for the mean last-interaction velocities ( vi ), listed in Table 5, of packets escaping into the solid-angle  obs , in a polar orientation ( vi,pole ) and an equatorial orientation ( vi,eq ).This is shown for packets of radiation escaping at all wavelengths, and for radiation that was last absorbed by the Sr ii triplet.
Epoch Time  obs (pole) ) is found at the positive pole (0.8 < cos() < 1), than at the negative pole (−1 < cos() < −0.8).The level of symmetry inferred from radiation escaping at all wavelengths Υ v,i , is similar to that shown by only radiation that was last absorbed by the Sr ii triplet, Υ Sr v,i .Even though the distribution of Sr in the ejecta is highly asymmetrical (as quantified by Υ Sr M and Υ Sr M, ), Υ Sr v,i is highly spherical.The high level of symmetry shown by the mean last-interaction velocities, as quantified by Υ v,i , supports the high level of symmetry inferred from  ∥ and  ⊥ at epoch 1, measured from the P-Cygni feature and EPM.Υ v,i remains highly spherical at epoch 2, showing that the simulation still exhibits a high level of symmetry beyond epoch 1.This implies that if the EPM was still valid at epoch 2 to infer  ⊥ then a high degree of sphericity could be found for this simulation..This demonstrates that even when the mass distribution and elemental abundances in the ejecta are asymmetric, the line forming region can be relatively close to spherical.

Inferring geometry from P-Cygni profile
Following Sneppen et al. (2023b), we analyse the shape of the P-Cygni profile in our synthetic spectra to explore how well it can be used to constrain the geometry of our simulation.Using the simple P-Cygni model, we show the comparison of P-Cygni features with spherical and elliptical photospheres 4 to our simulated spectra in 4 See Sneppen et al. 2023b for details of the modifications to the line profile calculator to represent an elliptical photosphere.Note that we did not apply Figure 7.We focus on matching the wavelengths of the absorption component of the P-Cygni profile and the blue side of the emission component, given that redward of the emission feature the spectra are not well matched by the blackbody distribution (which was not the case for AT2017gfo).We also note that the red side of the emission component is strongly blended with Ce iii in our simulation (see Shingles et al. 2023, figure 4).Across both epochs, the feature in the simulated spectra at the negative pole is best matched by a P-Cygni feature with a spherical photosphere, and the spectral feature at the positive pole is best matched by a P-Cygni feature with a prolate photosphere.This is surprising since the positive pole was inferred to show a higher degree of symmetry in Section 3.2.3(Table 3).Since opposite poles can not be fit by assuming the same geometry for the photosphere, this demonstrates that assuming a single photospheric geometry for the entire model is likely too simple an approximation.The range of  i (Figure 8) shows that radiation is emitted from a broad distribution of ejecta, which may not be well captured by assuming the photosphere can be modelled by a single velocity.We note, however, that it was possible for AT2017gfo to obtain a better fit with a P-Cygni profile than for our synthetic spectra (see Sneppen et al. 2023b) and this analysis should be tested for models more closely resembling a blackbody in future.
The distributions of Sr, Y and Zr, which are predominantly responsible for shaping the simulated spectral feature, do not show a spherical distribution since significantly lower masses of these elements are ejected at the equator (Figure B1).Since there is a direction for our asymmetric ejecta model at which the spectral feature is best the MCMC fitting used by Sneppen et al. (2023b) to fit the P-Cygni profile, but rather varied parameters manually to match the P-Cygni feature.i , weighted by the energy represented by the escaped packets of radiation, arriving at an observer at the pole ( obs : 0.8 < cos(  ) < 1.0) or equator ( obs : 0 < cos(  ) < 0.2) at 0.4 days.
matched by assuming a spherical photosphere, this demonstrates that the underlying ejecta does not necessarily have to be symmetrical for the P-Cygni profile shape to appear consistent with a spherical model.

P-Cygni profile shape for spherically symmetric ejecta
Having shown that an aspherical model can produce line profiles that appear consistent with spherical ejecta, we now explore whether the line shape predicted from a model with spherical ejecta is also well matched if fit with a simple P-Cygni model.Shingles et al. (2023) present a model where they enforce spherically symmetric ejecta by constructing a 1D, spherically averaged version of the 3D model we consider in this paper (1D AD2 from Shingles et al. 2023).Both the mass distribution and elemental abundances are spherically averaged.Imposing spherical symmetry on this model results in spectra and light curves that do not resemble those predicted for any observer direction from the full 3D simulation.The 3D structures in this model are important in shaping the spectra.Only when the 3D structure is included is the model able to produce synthetic observables that resemble AT2017gfo for this model.We note, however, that 1D simulations can produce spectra resembling AT2017gfo (e.g., Watson et al. 2019;Gillanders et al. 2022).
As with our 3D model, we fit the spectral feature produced by the 1D spherically symmetric simulation at ∼ 6000 -8000 Å with P-Cygni profiles assuming a spherical, prolate or oblate photosphere (Figure 10).Note that in this 1D simulation, the 'photosphere' (and spectrum-forming region) must be spherical.The best match is with the P-Cygni profile assuming a prolate photosphere.The emission component (on the blue side) increases too steeply to be matched by a spherical or oblate photosphere.However, given that the continuum is not well described by a blackbody, the exact contribution of the P-Cygni profile is more uncertain and less constrained than in the 3D case.Since we know that this simulation is spherical, this indicates that the geometry of the photosphere assumed to calculate a simple P-Cygni profile may not be a good test for spherically symmetric ejecta.
In both the 1D spherical simulation and the 3D asymmetric simulation, the spectral feature is actually a blend of features (primarily Sr ii, Y ii, Zr ii and Ce iii -see Shingles et al. 2023), and the simple Sr ii P-Cygni calculation here is likely unable to capture the complex interplay of these spectral features.

CONCLUSIONS
We have analysed a 3D radiative transfer simulation of a kilonova carried out by Shingles et al. (2023) to determine whether the simulation is compatible with the inferred symmetry constraints suggested for AT2017gfo by Sneppen et al. (2023b).
We have shown that although the ejecta from the neutron star merger model have moderate asymmetries in the mass ejected per solid-angle (e.g., Υ M =0.33), the synthetic light curves produced show a lower level of asymmetry (e.g.Υ bol = −0.12 at 0.4 days) than the level of asymmetry in the mass-distribution of the ejecta.The mass ejected within the velocity range that is approximately within the line forming region for the first day of the kilonova shows a high degree of symmetry (e.g., Υ M, = 0.037 at the negative pole), which may lead to a higher level of symmetry in the synthetic observables than suggested by Υ M , however, the distribution of elements synthesised in the ejecta in this velocity region is very asymmetric (e.g., Υ Sr M, = −0.71 at the positive pole).Radiation observed in a given line of sight is emitted from a broad range of regions in the ejecta, which has the effect of decreasing the observed anisotropy.For example, the Sr ii triplet absorption is highly spherical (e.g., Υ Sr ,i = 0.023 at the positive pole) despite the asymmetric distribution of Sr in the ejecta.However, the spectra are not predicted to appear the same at the poles as at the equator.The equatorial spectra are relatively featureless in comparison to the spectra produced near the poles (as discussed by Shingles et al. 2023).
Following Sneppen et al. (2023b), we quantify the level of symmetry that would be inferred from our simulation from the photospheric velocities,  ∥ and  ⊥ , obtained via simple fitting of the synthetic spectra (i.e., adopting similar methods as can be readily applied to real observations).At the first epoch considered (0.4 days), the values inferred for  ⊥ are similar to those obtained for  ∥ , indicating a high degree of sphericity, as quantified by Υ ,ph .From one pole, Υ ,ph is within the uncertainty inferred for AT2017gfo by Sneppen et al. (2023b).This demonstrates that the synthetic observables can appear .Spectrum at 0.4 days from spherically symmetric ejecta, compared to P-Cygni profiles (assuming the feature can be fit by a simple Sr ii P-Cygni feature) with a spherical, oblate or prolate photosphere.The dashed line shows a blackbody distribution.To plot the P-Cygni profiles, we modify the blackbody with the calculated line profile.However, we note that in the 1D spherical simulation, the spectra are not well represented by a blackbody.
consistent with spherical ejecta when viewed from certain directions, even when the ejecta are asymmetric.At the second epoch (0.6 days) the synthetic spectra are no longer well represented by a blackbody, and the inferred values of  ⊥ are too high.However, this is not what was found by Sneppen et al. (2023b) for the observations of AT2017gfo, where the spectra continue to resemble a blackbody until later times and the inferred photospheric velocities from AT2017gfo using the EPM were similar to those obtained from the Sr ii P-Cygni feature over the two epochs considered by Sneppen et al. (2023b).
The high degree of symmetry determined from the inferred photospheric velocities, quantified by Υ ,ph , at epoch 1 is supported by the level of symmetry in the mean last-interaction velocities extracted from the simulation, quantified by Υ v,i .The level of symmetry indicated by Υ v,i increases from epoch 1 to 2.
At epoch 1, while the spectra are well represented by a blackbody, the EPM can be used to infer the distance to our synthetic spectra to a good degree of accuracy (4-7 per cent).At epoch 2, however, the distance is underestimated (>25 per cent).In our simulation, where a higher degree of sphericity, Υ ,ph , was inferred from the synthetic spectra, a more accurate estimate of the distance D L was obtained.This may suggest that Υ ,ph could be used as a test for how accurate a distance estimate can be obtained from the EPM, however, this should be investigated with future simulations.
We compare simple P-Cygni profile models, assuming the photosphere to be spherical or elliptical, to our simulated spectral feature.The shape of the spectral feature produced at one pole is best matched by a P-Cygni profile assuming a spherical photosphere.However, from the opposite pole the shape was best matched by a P-Cygni model assuming a prolate photosphere.Additionally, we found that for a spectrum from a simulation with spherically symmetric ejecta, the feature was best matched by a P-Cygni profile assuming a prolate photosphere.These simulations suggest that fitting the profile shape alone may not be a robust test of spherical symmetry, although we note that the P-Cygni fits to our simulations are less certain than for AT2017gfo, due to the poorer match to a blackbody redward of the P-Cygni profile.
In our 3D simulation, there are lines of sight for which it could be inferred that the kilonova is highly spherical using the methods suggested by Sneppen et al. (2023b).The level of symmetry of the last-interaction velocities (Υ ,i ) of radiation escaping from our simulation is highly spherical, despite the asymmetric ejecta.This indicates that the combination of  ∥ and  ⊥ can indicate sphericity of the escaping radiation, however, this should not be interpreted to mean that the ejecta are necessarily spherically symmetric.Not all directions in the simulation would appear as spherical as inferred for AT2017gfo.If our simulation is representative of a real kilonova event, then we would expect that a high level of symmetry would not be inferred for all future observations (e.g., at the equator, where we do not predict a spectral feature from which to infer  ∥ ).However, if all future observations appear as spherical as AT2017gfo then this could suggest that this simulations is too anisotropic.More observations and simulations are required to understand the geometry of kilonovae.

APPENDIX A: DEVIATIONS FROM REFLECTION SYMMETRY
The merger simulation used in our artis calculation features a sizeable asymmetry with respect to the orbital plane.For instance, the ejecta masses in the Northern and Southern polar bins differ by a few 10% (see Figure A1).It is not straightforward to assess the extent to which this deviation from reflection symmetry is purely numerical or could be expected to occur in nature.
To better understand the asymmetries in our hydrodynamics tool, we consider a larger set of simulations.We analyse simulations where we vary the total binary mass in small steps for two different equations of state (SFHo and SFHx).The SFHX configurations cover a relatively large total binary mass range with a fine sampling rate.In all these simulations neutrinos are not included and we employ the Wendland SPH kernel (Schaback & Wendland 2006;Rosswog 2015).The ejecta masses for a given model differ from the simulations used in the main text, which include neutrinos and use the cubic spline kernel (Monaghan & Lattanzio 1985).In Figure A1 we show the total ejecta mass and the ejecta mass in the polar bins (as defined in the main text).Additionally, we consider only those fluid elements in the polar directions with a radial velocity in the range 0.15 ≤  ≤ 0.45.All simulations are analysed 12 ms after merging.Considering that we vary the total binary mass only within a relatively small range, we can expect that bulk properties of the merger would only change slightly.In Figure A1 we find that the different near-by models lead to stochastic increase/decrease of the ejecta in the Northern and Southern polar bin.The difference of the mass of polar ejecta are typically about ∼10% in these models implying that the asymmetry of our model used in the main paper may be somewhat on the extreme side of what is a typical symmetry breaking of the simulation code.It is reassuring that there is apparently no systematic trend favouring the ejection in a certain direction.
The magnitude of symmetry breaking visible in Figure A1 may be a result of the specific numerical scheme and is likely not informative about the amount of asymmetry in a real merger.Other numerical codes may lead to a different strength of symmetry breaking.Most grid-based merger simulations (presumably) impose reflection symmetry and by surveying the literature we could not easily compare our data to other calculations.The histograms of the ejecta mass at different polar angles in Palenzuela et al. (2015) and Lehner et al. (2016) also show a mild asymmetry of the mass ejection with respect to the orbital plane.It is not obvious to determine if a numerical scheme would overestimate asymmetries or possibly even damp the emergence of deviations from reflection symmetry.At any rate our simulations show that the exact ejecta distribution is apparently very sensitive to small perturbations.
In addition, we examine the strength of symmetry breaking in simulations with the moving-mesh code arepo (Springel 2010), which was recently extended to describe general relativistic systems (Lioutas et al. 2024).We consider equal-mass mergers for MPA1 (Müther et al. 1987;Read et al. 2009) ( tot = 2.5, 2.7, 2.9 M ⊙ ), DD2 (Hempel & Schaffner-Bielich 2010;Typel et al. 2010) ( tot = 2.7 M ⊙ ), SFHO (Steiner et al. 2013) ( tot = 2.7 M ⊙ ) and FSU2H (Kochankovski et al. 2022) ( tot = 2.8 M ⊙ ).In all cases, the microphysical EOS is modeled as a zero-temperature barotrope supplemented by a thermal ideal-gas component with Γ th = 1.75.We employ a resolution of  cell,0 = 1.5 × 10 −6 M ⊙ (see section 5.1 in Lioutas et al. 2024, for more details).For the system described by MPA1 with  tot = 2.7 M ⊙ , we perform additional simulations with lower ( cell,0 = 3.75×10 −6 M ⊙ ) and higher ( cell,0 = 6×10 −7 M ⊙ ) resolution as well.We find that these simulations also feature asym-metries of up to some 10% in the ejection along the polar direction.For the setups considered, we do not identify a favoured direction for the polar ejection nor any systematic trend with respect to resolution.Overall, this provides additional support by an independent code that there may be asymmetries with respect to the poles and that they might be of the order 10% in mass.
It is possible that in nature, as in our simulations, relatively small seed perturbations could be enhanced by the stochasticity of the violent and dynamical merger process and could lead to a considerable asymmetry of the ejecta, implying that a kilonova would appear differently from opposite poles.Initial perturbations could for instance stem from inhomogenities of the pre-merger magnetic field or from the intrinsic rotation of the neutron stars with one or both spin axes being misaligned with the orbital axis.The latter may introduce a significant mechanism for symmetry breaking even for relatively long rotation periods.If the asymmetries found in our models were representative of actual mergers, the kilonova calculations presented here provide an estimate of the variations that could be expected in nature.Generally, it may not seem unreasonable to expect at least some degree of asymmetry from kilonovae in nature.

APPENDIX B: SUPPLEMENTAL FIGURES
The mass distribution of select elements is shown in Figure B1, demonstrating that heavier, higher opacity elements are predominantly produced near the equator.This higher abundance of heavy elements does not lead to a significantly redder colour in the light curves at the equator, as demonstrated by Figure B2.In all directions, the light curves initially exhibit blue colours and evolve to redder colours over time.

Figure 1 .
Figure 1.Mass ejected by polar angle, where each bin has an equal solidangle.The height of each bar indicates the mass ejected within the solid-angle bin.The black lines indicate the mass ejected into the velocity range which is approximately in the line forming region of the kilonova within the first day (assuming homologous expansion).

Figure 2 .
Figure 2. Total mass (a) and mass of Sr (b) ejected into polar angle bins, where each angle bin has an equal solid-angle with a width of cos() = 0.2.The colour scale indicates the mass lying within each radial zone.Note a lower cut has been placed on mass to highlight the model structure.

Figure 3 .
Figure3.Bolometric light curves averaged over azimuthal angle, showing the variation with polar angle.Also shown is the zero-centred asymmetry index, Υ bol comparing the luminosity at the poles to the luminosity at the equator.In the lower plot, the blue line shows the negative pole compared to the equator (−1 < cos() < −0.8 and −0.2 < cos() < 0) and orange shows the positive pole compared to the equator (0.8 < cos() < 1 and 0 < cos() < 0.2).

Figure 6 .
Figure 6.Simulated spectra in polar and equatorial directions ( obs is listed in the figure) at 0.4 days (upper) and 0.6 days (lower) compared to the spectra of AT2017gfo at 1.43 days and 2.42 days, respectively.

Figure 7 .
Figure7.Model spectra at epoch 1 (0.4 days; left panels) and epoch 2 (0.6 days; right panels) for polar (upper and middle panels) and equatorial (lower panel) observer directions ( obs , listed in each panel).Dashed lines show the blackbody distribution (in the co-moving frame of the ejecta transformed into the rest frame of the observer, scaled to match the brightness of the synthetic spectra) best matching the spectra, where the ejecta temperature, T ′ , and the parallel photospheric velocity,  ∥ are listed in each panel.The blackbody distribution is modified to include a P-Cygni profile for Sr ii, considering spherical, prolate or oblate ejecta.For the spectra at the equator, indicative photospheric velocities are chosen for the relativistic correction to the blackbody distribution, since photospheric velocity cannot be inferred from a spectral feature at the equator.

Figure 9 .
Figure 9. Radial velocity and polar angle of the location in the ejecta ( ej ) where radiation was last absorbed by the Sr ii triplet ( Sr i ) before being re-emitted towards an observer in the direction,  obs , (arriving at the observer at 0.4 days) viewing towards the poles (a and b) and towards the equator (c).The arrow indicates the direction,  obs , of the observer.

Figure A1 .
Figure A1.Mass of ejecta from neutron star merger simulations with varying total binary mass assuming the SFHo or SFHx equation of state.The left panels show the total mass ejected by each simulation.The middle panels show the mass ejected into solid-angles at the Northern or Southern pole (with a solid-angle width of cos(  ) = 0.2).The right panels show the mass within the Northern and Southern solid-angle bins ejected with radial velocities ( r ) in the range 0.15c <  r < 0.45c.

Table 2 .
Asymmetry parameter comparing mass ejected within equal solid-angle bins near the equator and at the poles using the values in Table1.The range in polar-angle of each bin,  ej , (with a width of cos(  ) = 0.2) is listed.

Table 5 .
Mean ejecta radial velocities of the last interaction underwent by a Monte Carlo packet of radiation ( vi ) before escaping towards an observer in the direction,  obs .vi gives an indication of the spectrum-forming region in the ejecta.We show this for packets of radiation escaping at all wavelengths and for packets of radiation last absorbed by the Sr ii triplet ( vSr i ) before being re-emitted and escaping towards an observer.Also listed is the standard deviation  in  i and  Sr i .