{'fraction_non_alphanumeric': 0.03194139924257525, 'fraction_numerical': 0.021177994817620092, 'mean_word_length': 4.642114141130166, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 2, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dependence of Gravitational Wave Transient Rates on Cosmic Star Formation and Metallicity Evolution History

We compare the impacts of uncertainties in both binary population synthesis models and the cosmic star formation history on the predicted rates of Gravitational Wave compact binary merger (GW) events. These uncertainties cause the predicted rates of GW events to vary by up to an order of magnitude. Varying the volume-averaged star formation rate density history of the Universe causes the weakest change to our predictions, while varying the metallicity evolution has the strongest effect. Double neutron-star merger rates are more sensitive to assumed neutron-star kick velocity than the cosmic star formation history. Varying certain parameters affects merger rates in different ways depending on the mass of the merging compact objects; thus some of the degeneracy may be broken by looking at all the event rates rather than restricting ourselves to one class of mergers.


INTRODUCTION
Since the detection of the first confirmed gravitational wave compact binary merger events, astronomers and astrophysicists have been considering how to use this new window on the Universe to place constraints on its contents and our understanding of the underlying physics. (e.g. Abbott et al. 2016a;Abbott et al. 2017c,b;Abbott et al. 2019b). One of the most straightforward observables from GW events is their volumetric rate in the local Universe. A key test of stellar population synthesis codes is to reproduce that observed rate (e.g. Perna et al. 2018;Kruckow et al. 2018;Eldridge et al. 2019). To do this population synthesis codes predict a delay-time distribution: the expected event rate of GW transients versus time for a given amount of star formation. This is then combined with an assumed star formation history, along with its metallicity evolution, to predict a rate at the current epoch (see e.g. Langer & Norman 2006;de Mink & Belczynski 2015;Eldridge et al. 2019). There has been significant study of the how the uncertainties and assumptions in the stellar population models, especially the natalsupernova kick, affect the rate predictions and the observed double neutron star population (e.g Fryer et al. 1998;Wex et al. 2000;Dominik et al. 2013;Lipunov & Pruzhinskaya ntan032@aucklanduni.ac.nz † j.eldridge@auckland.ac.nz 2014; Abbott et al. 2016a;Tauris et al. 2017;Belczynski et al. 2017;Kruckow et al. 2018;Vigna-Gómez et al. 2018;Chruslinska et al. 2018;Mapelli & Giacobbo 2018;Eldridge et al. 2019;Andrews & Mandel 2019). However there are also substantial uncertainties in the star formation history of the Universe and its metallicity evolution which have been largely neglected until recently (Lamberts et al. 2018;Chruslinska et al. 2019;Mapelli et al. 2019;Artale et al. 2019;Neijssel et al. 2019).
In this letter we build upon our initial work in Eldridge & Stanway (2016) and Eldridge et al. (2019) to investigate how varying the early star formation history and metallicity evolution of Universe affects the predicted event rate of GW transients. We compare this to the changes in the event rates from changing some of the parameters within our stellar population synthesis models.   Figure 1. Comparison of our cosmic evolution models to observations. The upper panel shows the volume-averaged star formation rate density evolution, with observed values from (Madau & Dickinson 2014). Red points indicate estimates derived from infrared measurements, cyan points are derived from ultraviolet data. The lower panel shows our model for the median and 16-84th percentile range in metallicity as a function of redshift. These are compared to observed values derived from Damped Lyman-α (DLA) systems: green, blue and orange points from the compilation of Poudel et al. (2019), asterisks from Balestra et al. (2007), small crosses from Rafelski et al. (2012). Models are uniformly offset by -0.5 dex to account for the metallicity deficit of DLAs relative to typical star forming regions.

METHODS, OBSERVATIONS & SIMULATIONS
Population and Spectral Synthesis (BPASS 1 ) v2.2.1 code (Eldridge et al. 2017;Stanway & Eldridge 2018) were combined with a volume-averaged cosmic star formation history (Madau & Dickinson 2014) and a model for the evolution of metallicity in star forming regions (Langer & Norman 2006). We now build upon this method to investigate how varying the star formation history and/or the cosmic metallicity evolution changes the expected rate of compact binary mergers. Specifically we consider neutron star-neutron star (NS-NS), black hole-neutron star (BH-NS) and black hole-black hole (BH-BH) mergers.
We modify our method from Eldridge et al. (2019) in two ways. First, to allow for uncertainties in the high redshift cosmic star formation history we adopt a star formation rate 1 http://bpass.auckland.ac.nz density as a function of redshift as follows, ψ(z) = 0.015 (1 + z) 2.7 1 + ((1 + z)/2.9) u M yr −1 Mpc −3 . (1) This is a modification of the functional form given by Madau & Dickinson (2014), where we have replaced the exponent in the denominator of 5.6 with a variable u. This alters how quickly the density of star formation in the early Universe increases, while having little impact on the density of star formation and the total stellar density observed today. Nonetheless, it could change the GW event rate significantly due to the expected long delay times of GW events. We allow u to take values from 4 to 6, allowing for greater or lower star formation densities respectively at early times. Values of u in the range for 5.5 to 6 give the best agreement with observations as shown in Figure 1. However these data (compiled from multiple sources by Madau & Dickinson 2014) are derived from rest-frame ultraviolet measurements and so are subject to dust correction factors approaching 1 dex. Galaxies at the highest redshifts are currently believed to have very little dust, based on fitting to their spectral energy distributions (Bouwens et al. 2012) but the uncertainties on this are large (Wilkins et al. 2016(Wilkins et al. , 2018. If the intrinsic spectra of these sources are bluer (i.e. younger, lower metallicity or with a higher stellar rotation and/or multiple fraction) than currently estimated, the dust extinction will be underestimated, pushing the function towards lower u. We evaluate a large range of models to allow for this possibility.
In addition we calculate the fraction of star formation at different metallicities using the expression of Langer & Norman (2006): In this expression the value of β determines how quickly the Universe becomes enriched with metals, Z, as a function of redshift z and how broad the metallicity distribution is at each redshift. While in Eldridge et al. (2019) we used β = 2 here we allow this exponent to vary from 1 to 6. A higher β means the Universe was more quickly enriched, with a smaller metallicity scatter. In Figure 1 we compare the model metallicity enrichment for different βs to that determined from observations of Damped Lyman-α systems (DLAs). We offset the models by -0.5 dex to correct for the higher impact parameters (and thus lower measured metallicities) of DLAs relative to measurements in star forming galaxies (Møller & Christensen 2019). The scatter in the data suggests that β = 1 − 2 provides a better match to the observed spread than higher values, although the behaviour at the highest redshifts is largely unconstrained. If star formation occurs preferentially in already-enriched (i.e. more massive, older) dark matter halos in the distant Universe, rather than in sparse regions of the cosmic web, then a narrower range of metallicities might be expected for starbursts than is seen in the vast range of environments probed by DLAs. Again, we consider a broad range of β values to allow for this possibility. Over this cosmic history parameter space we calculate three model grids with different BPASS stellar population synthesis model sets. These are as follows: (i) BPASSv2.1 (Eldridge et al. 2017) models assume every star is in a binary with a flat distribution in mass ratio and the log of initial period. When a supernova occurs, a kick velocity is picked at random from the neutron star kick velocity distribution of Hobbs et al. (2005).
(ii) BPASSv2.2, Hobbs (Stanway & Eldridge 2018) models use the empirical binary population and parameter distributions of Moe & Di Stefano (2017) and are typically more robust for the old stellar populations dominated by low mass stars. The supernova kick is also picked at random from the distribution of Hobbs et al. (2005).
(iii) BPASSv2.2, Bray models use the same v2.2 initial binary parameter distributions as for (ii) but now we use the neutron star kick velocity from the work of Bray & Eldridge (2018): We show the predicted z = 0 GW event rates for these three BPASS model sets assuming the fiducial cosmic history parameters of β = 2 and u = 5.6 in Table 1. Varying the initial binary population, mass ratio and separation distributions has little effect on the NS-NS or NS-BH merger rate, but more than doubles the BH-BH merger rate. This is similar to results found by Mandel & de Mink (2016) and Belczynski et al. (2017). However changing the SN kick model has a more general effect, with the Bray kick (Bray & Eldridge 2018) increasing the NS-NS and NS-BH rates while the BH-BH merger rate is decreased. This indicates that the effects of these different stellar population synthesis assumptions on merger rates in different mass categories are orthogonal and that fitting all event categories at the same time will provide firmer constraints on the underlying physics of stellar models.
The NS-NS and BH-BH merger rates from models (i) and (iii) are consistent with observations from LIGO Observing Run 2 (which are themselves still subject to significant uncertainty, see Table 1), while in model (ii) the BH-BH merger rate is too high. The observed NS-NS merger rate derived from the GW 170817 event is higher than previously predicted (Abbott et al. 2017a). Thus while stellar population synthesis models (i) and (iii) are both consistent with the observed rates within the formal uncertainty, model (iii) is in best agreement with the data, assumping our fiducial star formation and metallicity histories. The next step, of course, is to vary this assumption.

RESULTS
We present the results of varying the early star formation rate parameter u and the metallicity evolution parameter β in Figure 2 and provide quantitative values in Table A1 in the appendix. The trends we observe are relatively simple.
Changing the early amount of star formation by varying u has the weakest impact of our two parameters. Decreasing u leads to greater early star formation and thus increases the number of low metallicity, massive early stars that can produce GW events after a long delay time. As a result the current-epoch rate of BH-BH and BH-NS events at fixed β shows a weak tendency to increase at low values of u.
In comparison the metallicity evolution parameter β has a much stronger effect, changing the merger rate by an order of magnitude for BS-NS and BH-BH mergers. The merger rates are highest at β = 1, indicating that GW transients prefer a low metallicity environment. In all cases a higher β leads to more metal-rich stars that are typically less efficient at creating BH-BH and BH-NS GW events (e.g Belczynski et al. 2010;Eldridge & Stanway 2016). This is because higher metallicity stars have stronger stellar winds so lose more mass creating less massive black holes as well as wider binaries that have longer merger times via gravitational waves.
NS-NS mergers in comparison show much weaker dependence on β. At the highest β values the merger rates reverse their decline and begin to increase instead, especially in models (i) and (iii). This is because as the stellar mass loss rates increase in high metallicity stars they become more likely to produce more neutron star than black holes at core collapse. The NS-NS merger rate increases at the expense of mergers involving black holes. We note that the almost flat dependence of the NS-NS merger rate on the star-formation parameters shows that this rate is more dependent on the stellar evolution model parameters than the star formation history.
This confirms that the relative rates of different classes of events from LIGO/VIRGO Consortium observations will have diagnostic power in distinguishing between physically motivated models for both stellar physics effects and cosmic evolution histories.
On each of the panels we include lines representing the current observational constraints on rates (Abbott et al. 2016b(Abbott et al. , 2019a. For the BH-BH mergers, models with β close to 2 give values in good agreement with observations, suggesting that our fiducial model is a reasonable match to the observed Universe, although a range of β from 1 to 3 are also consistent with the rate uncertainties. The current observationally inferred NS-NS merger rate is significantly higher than most extant population synthesis codes predict Mapelli & Giacobbo 2018;Chruslinska et al. 2019). However, novel kick schemes or a different metallicity evolution history of the Universe have been shown to push predictions close to the observed constraint, with the highest being that from using the kick of Bray & Eldridge (2018) as we see here. Both models (i) and (ii) are at the lower bound of observed rate while model (iii) is at the upper bound. This indicates that the assumed kick distribution may lie between these estimates and that GW event rates may prove an effective constraint on this distribution. Finally all cases with β = 2 yield predictions in good agreement with current upper limits on the NS-BH merger rate.

DISCUSSION & CONCLUSIONS
It is clear from our results that adopting different assumptions for the star formation and metallicity evolution of the Universe can have a significant effect on current epoch compact binary merger rates. These are as important, if not more so for the metallicity evolution, than changing the input physics of the binary interactions. We have found that the highest merger rates are obtained for low values of β with little dependence on the values of u. However this primarily affects the merger rates involving black holes, the NS-NS merger rate depends only very weakly on the star formation history. Thus by using the different merger rates it may be possible to break the degeneracy between binary evolution and cosmic history uncertainties. We note that our results are in line with other recent work in highlighting the importance of the metallicity evolution in predicting the correct merger rates (Chruslinska et al. 2019;Neijssel et al. 2019). For model (i) (BPASS 2.1, Hobbs), varying u and β changes the predicted NS-NS rate from the fiducial value by < 10 per cent in all cases except that of the lowest β which reaches 18 per cent at u = 4. By contrast the NS-BH rate can vary by as much as 67 per cent from the fiducial value and the BH-BH rate can vary by 90 per cent of the fiducial value.
For model (ii) (BPASS 2.2, Hobbs), varying u and β changes the predicted NS-NS rate from the fiducial value by < 7 per cent, the NS-BH rate by up to 77 per cent and the BH-BH rate by up to 80 per cent of the fiducial value.
For model (iii) (BPASS 2.2, Bray), the pattern is similar with the predicted NS-NS rate varying from the fiducial value by < 10 per cent across all parameters other than β = 1 which can reach a 43 per cent variation at low u. The NS-BH and BH-BH rates both vary by a factor of 2 relative to the fiducial value over the parameter range explored.
In all cases, the variation seen (within a factor of two) suggests that the current observation of the NS-BH and BH-BH rates would struggle to discriminate between cosmic histories, given their substantial uncertainties. Measurements with a precision of a few percent will be required to do so, which may also lie beyond the capability of the current observing run (LVC O3).
By contrast, we have demonstrated that discriminating between the underlying assumptions of the stellar population synthesis models may be possible with weaker observational data. Changing the initial binary distribution (from v2.1 to v2.2) modestly reduces NS-NS and NS-BH events by ∼5-10 per cent but causes a dramatic, factor of 2.6, increase in the BH-BH rates due to the larger number of close binaries occurring in a stellar population of a given total initial mass using the updated prescription. More pronounced still is the effect of changing the supernova kick prescription, which boosts the rate of NS-NS mergers by almost an order of magnitude, since more low mass systems survive their two-supernova evolution pathways without being disassociated. Importantly, we notice that the impact of changing cosmic history input parameters and physical model of stellar population synthesis vary by compact binary type. While interpreting event rates of one type (e.g. NS-NS) will lead to degenerate explanations, this degeneracy can be broken by comparing the rates of different types simulataneously. This implies that attempting to match the observed merger rate of these events at the same time will allow us to constrain our understanding of various aspects of stellar physics and also the evolution of the Universe.
A consistent picture from our results is that changing the assumed star formation history or metallicity evolution away from their current fiducial parameters generally leads to a decreased GW event rate. The predicted GW event rate only increases with an (unlikely) increase in the amount of stars formed in the early Universe. This would likely put the cosmic star formation rate density history in tension both with direct observations of the early universe (e.g. Bouwens et al. 2019) and with the observed history of mass assembly (which should follow its integral, see Madau & Dickinson 2014;Wilkins et al. 2019). Alternately a higher rate may be possible with a much slower metallicity evolution that is at odds with the observed early enrichment of the Universe (with relatively high metallicities being measured out to z > 3, see e.g. Sanders et al. 2019b, and Figure 1). The forthcoming James Webb Space Telescope will place much improved direct constraints on the properties of early star formation. In particular, it will improve observations of galaxies in the rest frame optical (with its relative insensitivity to dust extinction) and permit the true star formation rate density and metallicity at high redshift to be determined at unprecedented precision. While this work is already making progress from the ground at z ∼ 2 − 3 (e.g. Steidel et al. 2016;Sanders et al. 2019a), the sensitivity, wavelength coverage and multiplexing capability of NIRSPEC will permit determination of precise metallicities directly in the star forming galaxies which are responsible for the large delaytime tail of the GW event distributions.
We note that predictions from BPASS stellar population synthesis models, when combined with the adopted fiducial SFH and metallicity evolution of u = 5.6 and β = 2 (from Madau & Dickinson (2014) and Langer & Norman (2006) respectively) are in good agreement with the current LIGO/Virgo BH-BH and NS-NS merger rate estimates. The high NS-NS rate derived from GW 170817 favours model (iii), BPASS v2.2 with a Bray & Eldridge (2018) SN kick, in particular. If the NS-NS merger rate derived in the current and future LIGO/VIRGO observing runs is indeed as high as we predict, this will provide further support for adoption of a revised neutron star kick distribution in future work. Combining predicted event rates with the chirp mass distribution will provide further constraints on stellar and cosmic evolution models, although care will be needed to account for the impact of chirp mass on detectability of systems when comparing models to data. Table A1. Gravitational Wave event rates in Gpc −3 yr −1 to 3