{'fraction_non_alphanumeric': 0.03166465621230398, 'fraction_numerical': 0.02065741857659831, 'mean_word_length': 4.665432801822323, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Kilonova rates from spherical and axisymmetrical models

Detecting the thermal emission from double neutron star merger events is a challenging task because of the quick fading of the observed flux. In order to create an efficient observing strategy for their observing method it is crucial to know their intrinsic rate. Unfortunately, the numerous models existing today predict this rate on a vary wide range. Hence, our goal in this paper is to investigate the effect of different level of approximations on the relative rate predictions. Also, we study the effect of distinct ejecta mass lay-outs on the light curve. We find that the ratio of the expected kilonova detections of the spherical to axisymmetrical models is 6:1 (or 2:1, depending on the input parameter set applied in our work). Nevertheless, the light curve shape is only slightly affected by the various ejecta alignments. This means that different ejecta layouts can produce light curves with similar shapes making it a challenging task to infer the structure of the matter outflow. Thus, we conclude that the uncertainty in the rate predictions arising from the various ejecta mass distribution models is negligible compared to the errors present in other input parameters (e.g. binary neutron star merger rate). In addition, we show that up to moderate redshifts (z<0.2) the redshift distribution type (observed or uniform in volume) does not affect the expected relative rate estimations.


INTRODUCTION
The discovery of the GW170717 gravitational wave source (Abbott et al. 2017a) in association with the GRB 170817A (Abbott et al. 2017c) and its kilonova (hereafter KN) emission (Abbott et al. 2017b) was a significant step for multi-messenger astronomy. The observation of the multiwaveband radiation enabled detailed modelling of the underlying physics in order to constrain the physical parameters shaping the lightcurves.
Having an efficient observing strategy for a sky survey program (e.g. Large Synoptic Survey Telescope (hereafter LSST)) can greatly enhance the probability to detect KNe. There are numerous works which give an estimation for KN rates for the future sky surveys (e.g. Wollaeger et al. (2018); Tan et al. (2018) . When calculating these rates one has to take into account the short gammaray burst (sGRB) rate determined from observations (e.g. Sun et al. (2015); Wanderman & Piran (2015); Paul (2018); Ruffini et al. (2018) Nakar et al. (2006)). In addition, the particular model (e.g. composition, structure) employed in the study can significantly affect the results. Hence, it is not surprising that the predictions for KN rates cover a wide range of values depending what is chosen to be the source of a KN (all binary neutron star systems or only those producing short gamma-ray bursts (sGRB)).
Using an arbitrary detetection limit of mAB = 24.4 (the rband single-visit 5-σ limiting depth of the LSST) our aim in this work is to determine the relative KN rates calculated with spherical and axisymmetrical models. Moreover, the effect of models with different ejecta alignments on the relative rate is investigated.

LIGHTCURVE SIMULATIONS
The light curve simulations were carried out with the MOS-FiT software package. 1 This is an open-source Python-based code developed by Guillochon et al. (2018), which is able to generate synthetic light curves of various types of transient phenomena.

Spherical (V17) model
The spherical model, implemented by Villar et al. (2017) and first outlined by Metzger (2017), was motivated by multiband observations and is composed of two or three ejecta component: • "blue" ejecta: most likely it can be identified as a relatively proton-rich (high electron fraction, Ye) polar dynamical ejecta created by the shock from the collision between the merging neutron stars, with opacity of κ ≈ 0.5 cm 2 g −1 due to the lanthanide-poor material, • "purple" ejecta: probably a delayed neutron-rich (Ye < 0.25 − 0.3) outflow from the accretion disk formed in the merger, with opacity of κ ≈ 3 cm 2 g −1 , • "red" ejecta: essentially the same as the purple component, but composed of lanthanide-rich material, hence, the opacity is larger, κ ≈ 10 cm 2 g −1 .
The emission of the kilonova is powered by the radioactive decay of formerly generated r-process nuclei. However, only a fraction of this luminosity is converted into the observed flux: the thermal efficiency, th < 1, can be approximated with analytical functions derived from numerical simulations (Barnes et al. 2016). The detected spectrum is the sum of the blackbody radiation emerging from each ejecta component. Fitting the KN lightcurve of the GW170817 event V17 came to the conclusion that the three component model gives a slightly better fit to the light curve, thus in our work we use this three ejecta model.

Axisymmetrical (P17) model
Since the ejecta in the V17 model have spherical arrangement, as a second-order approximation, we implemented the axisymmetric ejecta structure into this model created by P17. The main characteristic of this approach is that the mass and opacity of the ejecta have a specific angular distribution inferred from general-relativistic hydrodynamical simulations. As in the spherical model above P17 uses three ejecta components: • dynamical (blue), the mass distribution can be approximated by F (θ) = sin 2 θ, Figure 1. The figure shows the structural alignment of the axisymmetrical model created by P17. The outermost layer is called the dynamic one, the middle layer is labeled as wind and the innermost layer is the secular one (see Section 2.2). In this model not only the ejecta have their specific mass distribution (see Eq. 1), but the opacity is changing along the polar axis. The blue color correspond to a low opacity matter, the orange indicates mild opacities, while the red marks the high opacity matter. The progenitor object may be a hypermassive neutron star or a black hole, around which an accretion disk can be formed. The sketch is based on Figure 2 from P17.
• wind (red), most likely polar emission with uniform distribution in mass, F (θ) ≈ const for θ θw ≈ π/3, • secular (purple), equatorial-dominated flow, F (θ) = sin 2 θ, where F (θ) characterizes the mass distribution according to where A is the normalization constant. The radiation mechanism only slightly differs from the one adopted in the spherical V17 model: the radiation emerging from the two innermost ejecta is reprocessed and emitted again by the outermost envelope. Apart from the structured ejecta and the reprocessed radiation of the innermost ejecta the model is identical to the spherical V17 one. A sketch illustrating the ejecta structure in the axisymmetrical P17 model can be seen in Figure 1.

Ejecta alignments
Since the assigning of a particular mass component in the spherical V17 model to a mass component in the axisymmetrical P17 model is not entirely straightforward, we generated KN samples with various mass distribution layouts (see Table 1). For example, in the P17 model the dynamic blue ejecta has an equatorial dominated mass distri- . The M disk indicates the mass of the disk produced by the massive neutron star before collapsing into a black hole. It has a range of 10 −2 M < M disk < 10 −1 M . The mass is in unit of 10 −2 M , the speed is in c and the opacity is in cm −2 g −1 . All quantities are sampled uniformly from their allowed intervals.

parset P17 parset C19
Parameter Range Range 0.5 -1 10 bution and the red wind component is enclosed into a cone around the rotational axis, in the literature it is a common practice to consider the inverted arrangement as well: the blue/dynamical ejecta is the polar component (beside the tidal tail part), while the red ejecta has the equatorial dominated mass distribution often identified as the wind component. In addition, some works suggest that both of the blue and red radiation component can have the same dynamical ejecta source (e.g. Fernández et al. 2015;Wanajo et al. 2014). Also, some authors propose (Arcavi 2018;Roberts et al. 2011) that because of the uncertainties in the photometric observations it is impossible to distinguish between the different models.

Parameter ranges
In Table 2 we show the two distinct parameter sets we used in the simulations. The first one, denoted as parset P17, was used by P17 to infer the physical parameters of the event AT 2017gfo. The second one, indicated as parset C19, was used by C19 for making predictions for the LSST and the WFD survey. In the case of the spherical V17 model we applied only the parset C19 set, however, in the case of the axisymmetrical P17 model we generated light curve samples using both of the parameters sets (parset C19, parset P17). All quantities are sampled uniformly from their allowed intervals.

Redshift distribution
When simulating the observed redshift distribution of short GRBs from the (a) intrinsic redshift distribution of their progenitor compact merger system one has to take into account the (b) time delay between the formation of the system and the inspiral. The former one, (a), more or less traces the star formation history (e.g. Wanderman & Piran 2015), while the latter one, (b), is usually approximated with empirical functions. Such an approximation was done by Sun et al. (2015), who considered the star formation history from Yüksel et al. (2008), and convolved it with three different time delay distribution models: the power-law, the Gaussian and the lognormal delay model. However, the power-law delay model has been shown to be inconsistent with the observations (e.g. Tan et al. 2018). We used the formula for the lognormal time delay model they give to simulate the redshift distribution for our kilonova samples (Eq. 21 in Sun et al. 2015). The maximum distance of events in this work is Dmax ≈ 1 Gpc (z ≈ 0.2).
Since the maximum brightness of the KNe simulated in various works makes them detectable only up to z 0.1 it is a common practice to neglect the shape of the redshift distribution function and assume a uniform distribution in comoving volume in the Universe with a density (e.g. Wollaeger et al. 2018;Scolnic et al. 2018;Saleem et al. 2018). Thus, we also generated light curve samples for both the V17 and P17scs models where the redshift is distributed uniformly in comoving volume between 100 and 740 Mpc.

Binary neutron star merger rate
There are basically two different approaches when considering double neutron star merger rates: determining the rate from short GRB observations or calculating it by employing binary compact object population synthesis methods. The former method is less reliable, since it is believed to suffer from the uncertainties of various parameters, like the beaming factor, redshift, minimum luminosity and the time-delay distribution. If we assume that all of the short GRBs produce KNe then the rate lies between 0.2 (Ghirlanda et al. 2016) and 40 Gpc −3 yr −1 per f −1 b = 1 − cos θj (Nakar et al. 2006), where f b is the beaming factor and θj is the jet halfopening angle. Much higher KN rate can be derived if we assume that all of the binary neutron star merger events result in KN, namely as high as 1540 Gpc −3 yr −1 (Abbott et al. 2017c). The lowest value in this case is 316 Gpc −3 yr −1 , which is taken from Della Valle et al. (2018) and is based on the work of Belczynski et al. (2008). The minimum and maximum short GRB and binary neutron star merger rates are listed in Table 3.

RESULTS
Having the simulated KN sample we can now determine the expected KN rate and compare these rates in the light of the different models. Table 3. The lowest and highest short GRB rates (Gpc −3 yr −1 per f −1 b = 1) and binary neutron star merger rates (Gpc −3 yr −1 ) in the literature cover a wide range of values. We take these numbers as lower and upper limits when calculating expected KN rates.

Spherical vs. axisymmetrical model
Although, as pointed out in V17 their spherical model captures the main characteristics of a KN, there are two main drawbacks. First, the spherical structure of the ejecta does not reflect the angular alignment of the progenitor stars rotational axis. The result of this is that when the observer faces the polar or equatorial region of the progenitor the radiation coming from an otherwise obscured ejecta component is thus overestimated, or similarly, underestimated. This can lead to false predictions for the underlying physical parameters, e.g. ejecta mass. The difference is clearly visible on the r-band light curve (see Figure 2). Second, allowing only the mass to vary (see Table 2) strongly narrows down the potential progenitor system's physical parameter space and hence the width of the resulting maximum brightness distribution. This can be seen in Figure 3 and in Figure  4: the former one shows the maximum brightness against Figure 3. The plot shows the maximum r-band brightness against the redshift for the distinct samples. It is observable that allowing only the ejecta mass to vary in the spherical V17 model (red dots) results in a very narrow maximum brightness distribution in the r-band. In the case of the axisymmetrical P17 model, the parset C19 sample denoted with green triangles, produces brighter events contrary to the parset P17 sample (blue dots) where the input parameter values are lower compared to the former case (see Section 2.2.2).
the redshift of the V17 and P17scs samples, while the latter one displays the maximum brightness histograms. As it can be observed in Figure 3 the two different parameter sets (Table 2) used for the axisymmetrical P17 model produces slightly distinct maximum light curve brightness distributions. While the parameters from P17 (parset P17) can generate events (marked with blue in Figure 3) as bright as the spherical V17 model only up to z ∼ 0.1, simulations with the higher input values adopted by C19 (parset C19) can clearly reproduce the bright events (indicated with green) observed also with the spherical V17 model (red in Figure 3) up to z 0.2. The lower brightness of events with higher redshift is a consequence of three factors which are related to each other: first, the shifting of the spectrum out of the observed band, second, the different contributions of the blue, purple and red ejecta, and third, the observing angle.

Mass distribution lay-outs
As we mentioned it earlier in Section 2, the KN models commonly used in the literature differ not only in the structural composition (spherical vs. axisymmetrical), but in the different ejecta mass distribution types.  of the ejecta components, changing the opacity accordingly. Our simulations suggest that the various ejecta configurations do not effect significantly the light curve shapes. This result is observable in Figure 5, where the light curves with different colors and linetypes correspond to different models: the thick grey line shows the spherical V17 case, the dashed red line denotes the P17ssc model, the green dashed-dotted line marks the P17scs case and the blue two-dashed line indicates the P17css. Nevertheless, examining the maximum brightness of the light curves it can be seen that the difference in the maximum brightness can reach up to ∼ 1 magnitude: the P17css model produces dimmer events than the P17ssc model. The steeper rising and decreasing of the observed flux along with the earlier peak time might be the consquence of fixing the opacity which should be time and wavelength dependent. However, since we are interested only in the maximum brightness of the events (the ratio of the detectable events simulated with the spherical V17 and axisymmetrical P17 models depends only on the maximum brightness), we believe that the former shortcomings do not affect our results. Improving our model with time and wavelength dependent opacity is the scope of a future work.

Expected kilonova rate ratios
The main problem when calculating the expected KN rate for a specific sky survey program, in our opinion, is that the uncertainties of the input parameters (e.g. binary neutron star merger rate, mass component distribution) are relatively high making the estimation unreliable even in the case when the observing cadence is simulated properly. Because of this we think it is more appropriate to calculate the relative predicted rate of KNe simulated with different approximations. However, we still need an arbitrary detection limit in order to be able to make an estimation. For this purpose we use the detection threshold of mAB = 24.4, which is the r-band single-visit 5-σ limiting depth of the LSST (Ivezić et al. 2008).
Now we can compare the expected rates from different approximations based on the following aspects: a) spherical V17 vs. axisymmetrical (P17) model -we find that the ratio of the rates with the spherical to the axisymmetrical model is 6 : 1. However, if the higher input parameter values are adopted from C19 the ratio becomes smaller, 2 : 1.
b) observed vs. uniform in comoving volume redshift distribution -considering the spherical V17 model, the ratio of detectable kilonovae simulated with the uniform redshift distribution in volume up to z 0.15 to kilonovae simulated with the observed redshift distribution is ∼ 1. This is also true for the axisymmetrical (P17) model. c) different mass distribution lay-outs -the relative rates of events generated with different mass distribution lay-outs cover a relatively wide range: rate(P17scs)/rate(P17ssc) ≈ 0.6, rate(P17scs)/rate(P17css) ≈ 2, rate(P17ssc)/rate(P17css) ≈ 3. These ratios can be explained by the fact that the discrepancy of the maximum brightness of the events simulated with the different models can be as large as ∼1 magnitudes. However, the light curves have similar shapes for the distinct models. This can be observed in Figure 5.

CONCLUSIONS
In this paper we simulated KN events following double neutron star merger events (along with short gamma-ray bursts). In order to investigate the effect of the structural compositon of the underlying physical model on the relative kilonova rate predictions we created samples with spherical and axisymmetrical models. In addition, different ejecta composition lay-outs were explored to analyze their impact on the light curve shape. Since the kilonova rate predictions suffer from large errors because of the uncertainties of the input parameters we calculated only the relative rates of the different models. Our results show that • the ratio of the rates of the events calculated with the spherical and axisymmetrical ejecta structure can be as small as 2:1. This means that the structural alignment of the kilonova ejecta is not a significant factor when determining the expected kilonova rates, • changing the distribution type of the distinct mass components does not alter significantly the light curve shape. Nevertheless, the difference in maximum brightness of the light curves produced with different lay-outs can reach up to ∼1 magnitudes, • our simulations support the common assumption that a uniform in comoving volume redshift distribution can be used when simulating KN events because of their proximity: the number of observable KNe (up to an arbitrary detetction limit) generated with the spherical V17 and axisymmetrical P17 model is the same with uniform and observed redshift distribution.
Based on our work it can be seen that when calculating the expected kilonova rates the chosen physical model (spherical, axisymmetrical) does not affect significantly the results. In contrast, the uncertainties in the model input parameters, such as the binary neutron star merger rate (where the predictions cover 2 ∼ 3 orders of magnitude), play a much more important role in determining the expected number of events. Thus, in the light of the above results we think that the ultimate test for the kilonova rates will be made by the upcoming sky surveys.