GW and EM transient rates

A consistent estimate for Gravitational Wave and Electromagnetic transient rates

J. J. Eldridge , E. R. Stanway and Petra N. Tang
Department of Physics, University of Auckland, Private Bag 92019, Auckland, New Zealand
Physics Department, University of Warwick, Gibbet Hill Road, Coventry, CV4 7AL, United Kingdom
j.eldridge@auckland.ac.nz
Accepted XXX. Received YYY; in original form ZZZ
Abstract

Gravitational wave transients, resulting from the merger of two stellar remnants, are now detectable. The properties and rates of these directly relates to the stellar population which gave rise to their progenitors, and thus to other, electromagnetic transients which result from stellar death. We aim to estimate simultaneously the event rates and delay time distribution of gravitational wave-driven compact object mergers together with the rates of core collapse and thermonuclear supernovae within a single consistent stellar population synthesis paradigm. We combine event delay-time distributions at different metallicities from the Binary Population and Spectral Synthesis (BPASS) models with an analytic model of the volume-averaged cosmic star formation rate density and chemical evolution to determine the volume-averaged rates of each event rate at the current time. We estimate rates in excellent agreement with extant observational constraints on core-collapse supernovae, thermonuclear supernovae and long GRBs. We predict rates for gravitational wave mergers based on the same stellar populations, and find rates consistent with current LIGO estimates. We note that tighter constraints on the rates of these events will be required before it is possible to determine their redshift evolution, progenitor metallicity dependence or constrain uncertain aspects of stellar evolution.

keywords:
Methods: numerical – Gamma-ray burst: general – supernovae: general – Gravitational waves
pubyear: 2018pagerange: A consistent estimate for Gravitational Wave and Electromagnetic transient ratesA consistent estimate for Gravitational Wave and Electromagnetic transient rates

1 Introduction

The transient sky is being studied in ever more detail, by a rapidly increasing number of facilities and surveys. These identify astrophysical sources which vary significantly in luminosity on short timescales - typically ranging from minutes to months. While many such variations are associated with recurrent sources, such as cataclysmic variables or accretion onto active galactic nuclei, others are associated with the dramatic energy release events which accompany stellar death.

The dramatic types of stellar death fall into two broad categories, based on of the driving processes that cause them. Either they are the result of evolutionary processes affecting the core of stars, or they result from orbital evolution as gravitational wave emission carries away energy and angular momentum, driving two compact remnants into contact. The first is best exemplified by the cessation of nuclear burning in an iron-group elements (or oxygen-neon) core of a massive star that then undergoes core collapse resulting in the emission of an optically-luminous transient classified as a supernova (SN). Where hydrogen remains in the stellar envelope these are classified as Type II, while hydrogen-free core-collapse supernovae (CCSNe) are classified as Type Ib or Ic. In a subset of stripped-envelope core-collapse events, a relativistic jet may be launched causing beamed, high-energy emission which may be detected as a Long Gamma-Ray burst (LGRB) by observers on the jet emission axis (e.g Heger et al., 2003; Langer, 2012; Smartt, 2015).

Spectacularly, the LIGO and VIRGO facilities have now opened a new window onto the transient sky, detecting the gravitational wave (GW) strain variations associated with the second category of transients: the merger of compact binaries. In these events, either two black holes (BH) or neutron stars (NS) lose angular momentum through GW emission and spiral towards a final merger and associated ‘chirp’ event (Abbott et al., 2016, 2017b). While these events had been theorised, and the mechanism of short GRBs (SGRBs) tentatively assigned to GW-driven NS-NS mergers (Paczynski, 1986), the detection of gravitational wave chirps has finally provided conclusive evidence for this interpretation.

Type Ia SNe may have progenitor pathways that belong in both these categories. White dwarfs, supported by electron degeneracy pressure, form the end state in the life of 95% of stars. However either accretion onto a white dwarf from a companion or the merger of two white dwarfs in a close system, due to gravitational waves, can cause the star to exceed the Chandrasekhar limit and if the conditions are correct ignite, giving rise to a thermonuclear detonation known as a Type Ia SN (e.g Howell, 2011; Maoz et al., 2014).

Current supernova surveys such as PanSTARRS (Chambers et al., 2016), ATLAS (Tonry et al., 2018) and ASAS-SN (Holoien et al., 2017) report large samples of supernovae on a weekly basis, while forthcoming surveys, notably those associated with the Large Synoptic Survey Telescope (LSST, Abell et al., 2009), will produce data on thousands such sources per night. Gamma-Ray bursts are rarer, being detected at a rate of a few per week. GW transients are currently the rarest known, although this is likely a result of the enormous technical difficulty of detecting such events rather than their intrinsic rate.

While the mechanism behind each type of transient, the progenitor stars and the detection methods vary, fundamentally all these events represent a probe of the stellar population which gives rise to them, and our understanding of the stellar evolution processes involved. Their relative rates in a population of known age and metallicity can provide insight into the relative fraction of stars in different mass ranges. As the sample grows over the coming years it will be possible to invert this, by assuming an initial mass function and adopting models for the delay time between stellar birth and transient event, to yield information on the cosmic star formation history and metallicity evolution.

The strongest such constraints will be obtained from a joint analysis of all available transient samples. However simultaneous modelling of the varying progenitor systems is a challenging task, requiring detailed handling of the physics of stellar evolution across a broad mass range. Crucially, it also requires the modelling of stellar binary interactions, which play a key role in the evolution of nearly all transient progenitors (e.g. Podsiadlowski et al., 1992; De Donder & Vanbeveren, 1998; Langer et al., 2000; Belczynski et al., 2008; Mennekens et al., 2010; Ruiter et al., 2011; Eldridge & Stanway, 2016; Mennekens & Vanbeveren, 2016; Zapartas et al., 2017; Wang et al., 2017; Liu et al., 2017). Few stellar population synthesis models exist which are capable of a simultaneous, systematic, uniform analysis of the CCSN, GRB, thermonuclear SN and compact merger transient rates arising from a stellar population.

In this paper we utilise the Binary Population and Spectral Synthesis (BPASS) code to present a simultaneous estimate of the transient event rates expected in both gravitational waves and electromagnetic light, arising from the same underlying fundamental stellar population, which reflects the star formation and chemical evolution histories of the Universe on cosmic scales. The paper is structured as follows: in section 2 we describe the stellar population synthesis models used for this analysis and the delay time distribution arising from single-metallicity populations; in section 3 we produce a metallicity and star formation weighted model for the volume-averaged stellar population as a function of redshift and compare this to observational constraints on the rates of different transient types; in section 4 we discuss the implications of our analysis and prospects for future observatories. Finally in section 5 we present a summary of our conclusions

Where necessary we assume a standard CDM cosmology with =0.3, =0.7 and =100  km s Mpc with =0.7

2 Stellar Transients Population Synthesis

Figure 1: The predicted delay time distribution derived from BPASS for electromagnetic and gravitational wave transients at two representative metallicities. We compare the predicted delay times at both metallicities with the same set of observational constraints. The square points are the observed core-collapse supernova rates reported by Maoz & Badenes (2010). The remaining data points indicate type Ia SN rates from Maoz et al. (2010, crosses) and Totani et al. (2008, triangles).

2.1 Binary Population and Spectral Synthesis

To determine the relative rates of different transient types arising from a stellar population of known age and metallicity we make use of the Binary Population and Spectral Synthesis (BPASS) models111accessible at bpass.auckland.ac.nz. These combine a set of custom, detailed, theoretical stellar evolution tracks at thirteen input metallicities. Crucially a substantial subset of the evolution tracks include effects due to binary interactions, with both stars in a binary system of known initial period, mass and mass ratio evolved through detailed modelling in a 1D spherical approximation. Prescriptions are adopted not only for stellar winds but also for mass transfer through Roche Lobe overflow, common envelope evolution, rejuvenation by mass transfer and, more rarely, resulting chemically homogeneous evolution.

The stellar evolution tracks are combined through weighting by an initial mass function (for which we use a broken power-law based on the Kroupa (2001) IMF) and by an initial period and mass ratio distribution as a function of primary star mass (for which we use the empirical distributions determined by Moe & Di Stefano (2017)). The underlying physical prescriptions and their applications are fully described in Eldridge et al. (2017) and Stanway & Eldridge (2018). The current version is v2.2, and we present results using our default IMF with a Salpeter-like upper mass slope of -2.35 and a maximum mass cut off of 300 M. We use the 13 metallicities that have been computed to date that range from to , and assume .

There are two key points to note about binary interaction models in BPASS. First, for common envelope evolution we use our own prescription as described in Eldridge et al. (2017). It is unique because we use a detailed stellar evolution code so cannot simply remove the hydrogen envelope in one time-step. We therefore remove the envelope as quickly as possible, and also remove its binding energy from the binary orbit. This means we naturally take account of the structure of the star in estimating the outcome of the interaction. Second, during more mild interactions of Roche lobe overflow the companion star is limited in how much it can accrete. For normal stars the limiting factor is the thermal timescale, for white dwarfs and neutron stars we limit the accretion rate by the Eddington luminosity. We do not limit the accretion onto black holes.

BPASS models have previously been used to estimate GW transient properties from simple stellar populations (Eldridge & Stanway, 2016) and to explore supernova rates (Xiao & Eldridge, 2015) and progenitor populations (Xiao et al., 2018). They have also been used to evaluate the properties of GRB host galaxies (e.g. Stanway et al., 2015; Eldridge et al., 2017)

Stellar death is tracked in a number of ways in the models, which identify core-collapse, the formation of remnants and the further evolution of the secondary star beyond the death of the primary. Many remnants are formed with natal kicks as a result of impulse from a supernova. We adopt a prescription for kicks as described in Hobbs et al. (2005). In brief, we use a Maxwell-Boltzmann distribution with a velocity of 265, sampling this distribution through repeated model iterations to understand the full range of possible outcomes from each SNe. We do not currently employ the kick of Bray & Eldridge (2016) or Bray & Eldridge (2018). This kick will only have an impact on our predicted GW event merger rates as shown in Bray & Eldridge (2018), and not on supernova rates.

Key pathways for stellar transients are as follows:

2.1.1 Core Collapse

We identify a core-collapse event if central carbon burning has taken place and the CO core mass exceeds 1.38 M while the total stellar mass is greater than 1.5 M. We classify the event type as described in Eldridge et al. (2011); Eldridge et al. (2013) and Eldridge et al. (2017). In brief, an event is classified as type II if the mass of Hydrogen is  M, and Type Ib or c otherwise.

The remnant mass is determined by calculating the residual mass bound to the star after ejection of material given a typical supernova energy injection of  ergs. This method was first outlined in Eldridge & Tout (2004) which provides more details. Core-collapse events which produce a remnant mass exceeding 3 M, and in which the progenitor has experienced chemically homogeneous evolution due to mass transfer at very low metallicity (), are deemed to explode as long gamma-ray bursts (see Eldridge et al., 2017, for details). While this is a probable LGRB progenitor pathway, there are likely others which are not identified within BPASS. In particular, this mechanism is unable to produce progenitors in stellar populations with a metallicity above 0.2 Z, while such events are known to occur within the observed sample. We intend to explore other pathways in future but improved treatment of stellar rotation, and possibly tidal forces, will be required. The current BPASS LGRB rate estimates should therefore be treated as lower limits.

Finally if the helium core mass at the end of evolution is between 64 and 133 M  a pair-instability supernova (PISN) is assumed to occur which completely disrupts the star and leaves no remnant (Heger & Woosley, 2002).

2.1.2 Compact Binary Coalescence

When the second star in a binary system forms a compact remnant, the orbital evolution of the resultant compact binary is presumed to be dominated by gravitational wave radiation. For each compact binary formed we identify the total delay time, which incorporates both stellar lifetime and the timescale for merger given the remnant masses and the binary period at compact object formation. We adopt the Peters (1964) prescription for evolution of binary orbits due to gravitational wave emission.

As in Eldridge & Stanway (2016) we use the near-circular and highly eccentric approximate formulae from Peters (1964) and interpolate between these for intermediate eccentricities. We note that we may be underestimating the eccentricity of the initial binary because after the first SN in a binary we assume the orbit is circular.

We calculate rates separately for WD-WD, NS-WD, NS-NS, NS-BH and BH-BH compact binaries. All of these are technically sources of gravitational waves, although for WD-WD and NS-WD systems the frequency range and signal strength puts mergers below the sensitivity threshold of current and planned detectors. WD-WD mergers are considered type Ia supernovae as described below. The remaining combinations are considered GW transient sources. NS-NS, and perhaps NS-BH, systems are also likely to be electromagnetic transients, with on-axis jetted emission giving rise to a short GRB and more isotropic emission in the form of a fainter kilonova (Abbott et al., 2017d; Tanvir et al., 2013).

2.1.3 Type Ia (Thermonuclear) events

We identify Type Ia thermonuclear supernovae which occur through two complementary channels, both of which result from binary interactions. In the single-degenerate channel a white dwarf accretes material from a main-sequence, brown dwarf or giant companion and detonates when it reaches the Chandrasekhar limit. Here we assume that if a white dwarf with a mass less than 1.2 accretes material and reaches a mass of 1.4 the star will explode as a type Ia SN.

The double-degenerate channel causes a type Ia supernova when two white dwarfs evolve and remain bound in a binary system. Gravitational wave radiation is then assumed to dominate the dynamical evolution of the binary and the merger timescale is calculated using the same method as for NS-NS, NS-BH and BH-BH binaries but assuming a circular orbit for the binary.

Our rates are dominated by the single-degenerate channel. Other channels are hypothesised for subsets of type Ia events, which result in partial detonations (e.g. the sub-luminous events which give rise to hypervelocity white dwarfs, Raddi et al., 2018). However these represent a small fraction of the observed population.

One key aspect of modelling type Ia SNe from the single-degenerate channel may explain several of these minor pathways: the accretion rate of material onto the white dwarf (e.g Nomoto & Kondo, 1991; Langer et al., 2000; Kato & Hachisu, 2004; Ruiter et al., 2009; Ruiter et al., 2011; Wang et al., 2017). If the rate is too low then material builds up on the surface of the white dwarf until it ignites in a nova. If it is too high then the Eddington luminosity can be exceeded or a hydrogen envelope can be formed around the white dwarf returning it to the giant branch. There is however a narrow accretion rate range where the material can burn at the surface of the white dwarf as it is accreted, thus increasing the mass of the white dwarf. This mass transfer rate is estimated as around a few (e.g Ruiter et al., 2011). The average accretion rate in our models from start of mass transfer to explosion cover a range of mass transfer rates from to a few . However the models that contribute most of the weight to our Ia rates lie within the narrower range of to a few .

While we still stress that our predicted type Ia rates should be used with caution, the fact they lie within the expected accretion rate window for white dwarf growth, as well as matching observational constraints as shown in this article, suggest that they represent a fair estimate. Further analysis of the progenitor systems is beyond the scope of this paper. The improvement in our low mass stellar grid and prescriptions in the current BPASS v2.2 (described in Stanway et al. 2018) included the calculation of a significantly larger grid of stellar models which produce WD binaries, while the revised binary parameter distributions also adopted in v2.2 have further improved our estimates for type Ia rates, which we now consider a robust output of the BPASS code.

2.2 Delay Time distribution

A key prediction of any stellar population synthesis model is the distribution in delay time, i.e. the total time required from the onset of star formation to allow stellar evolution to occur and a transient to be triggered. This distribution will differ by transient type and also with metallicity and as a function of binary interaction parameters.

In figure 1 we show the delay time distributions predicted by BPASS for electromagnetic and gravitational wave transients at two representative metallicities (Z and 0.1 Z). We compare the predicted delay times at both metallicities with the same set of local Universe observational constraints (for which detailed metallicity information on progenitors is unavailable). The square points are the observed core-collapse supernova rates reported by Maoz & Badenes (2010) as a function of host stellar population age. The remaining data points indicate type Ia SN rates from Maoz et al. (2010, crosses) and Totani et al. (2008, squares).

Core collapse supernovae occur in stellar population ages as high as  Myr in a binary evolution model, as opposed to the few 10 Myr limit in single star evolution models. This is a consistent prediction between binary population synthesis models as shown by Zapartas et al. (2017). Furthermore as would be expected, the sources requiring more massive progenitors (LGRBs and PISNe) are restricted to younger ages and lower metallicities.

Thermonuclear supernovae begin to appear in a population at ages of  Myr, but extend up to (and indeed beyond) the age of the Universe. For both core collapse and type Ia events, the distribution of observed rates in the local Universe requires the presence of binary pathways in the transient progenitors, and is in good agreement with BPASS predictions. The high rate of observed type Ia events at short timescales of a few 100 Myr also marginally favours a sub-Solar mean metallicity for their progenitors.

By contrast the delay time distribution in gravitational wave transients is considerably broader, with the most rapid, prompt events occurring within a few Myr of the formation of their progenitor stars, and the distribution extending to  Gyr. We note that while event rate per year is low at late times for a given initial mass of star formation, the large width of these time bins means that they actually dominate the total number of events from the population.

3 Cosmic Evolution of Transient Rates

3.1 Input Cosmic History

The star formation history and chemical evolution of individual galaxies are sensitive to their merger trees and environment. A detailed estimate of transient rates in a given galaxy, such as the Milky Way, must account for this detailed evolutionary history. However on sufficiently large scales, such that both dense and sparse regions of the large scale structure are well sampled, these variations can be averaged out to determine a smooth evolution. Madau & Dickinson (2014) considered a compilation of galaxy survey data to determine analytic forms for the volume-averaged evolution in both star formation rate density (SFRD) and metallicity mass fraction (Z) over cosmic time.

Here we use the analytic expressions for SFRD, , from Madau & Dickinson (2014) to follow star-formation through cosmic history. This is given by,

(1)

We combine this with the analytic expressions for metallicity evolution used by Langer & Norman (2006). Here the fractional mass density of star formation at and below metallicity mass fraction of is given by,

(2)

where and are the incomplete and complete Gamma functions and we have used the fiducial values for the model from Langer & Norman (2006). We use this parameterization because it provides a distribution of metallicity at each redshift, rather than other methods requiring a spread around the mean of a metallicity.

To calculate the the rate of transients at each redshift step, , we first calculate the star-formation density at each redshift bin. We then split this star formation between each of our input metallicities. We also calculate the time interval represented by each redshift bin, as well as the age of the Universe in each redshift bin. We then convolve this with the look back time distributions of each transient type at each metallicity, starting at each redshift bin and working back through time to work out the contribution to the current redshift bin from all previous star-formation.

From the resultant, mixed metallicity, mixed age population, we calculate the volume-averaged mean rate of each transient type, , as a function of redshift:

(3)

where represents the metallicity distribution at and R() is the final event rate per year per () Gpc.

We show the resulting star-formation rates in each of our metallicities versus look-back time in Figure 2. As we can see today the higher metallicities dominate current star formation. At the peak of the star-formation rate history in the Universe, the metallicity is dominated by the approximately half-solar metallicities. Only at the earliest times do the lowest metallicities dominate the star-formation.

Figure 2: The input cosmic star-formation rate and metallicity evolution used to calculate our rates through cosmic history. The dashed line represents the total star-formation rate of all metallicities combined. Each line is colour coded to a specific metallicity.

3.2 Rate evolution

Given this mixed stellar population, we present the rate evolution of each transient type against redshift and lookback time in figure 3.

Core collapse SNe, which arise from massive stars with short delay timescales and which show a relatively weak metallicity dependence, closely track the volume-averaged mean star formation rate. Their volume-averaged population peaks just below , an epoch known as Cosmic Noon due to the intense galaxy merger, AGN and star formation activity at this redshift.

Type Ia SNe are the second most frequent event type at most redshifts. The long delay time associated with these events (see figure 1) shifts the peak in event rate to lower redshifts, , below which the rate declines. At , the volume-averaged type Ia rate is 21% of that for CCSNe. At the current time the type Ia rate is dominated by the single-degenerate channel. The rate of the much slower double-degenerate channel is still rising and it will dominate the rate in the future.

By contrast, the strong metallicity dependence of long-GRB progenitors in our model shifts the peak in GRB emission rate to higher redshifts, with a near-constant volume-averaged rate at . The inclusion of higher metallicity GRB pathways may modify this redshift evolution in future, but will have little effect on the shape of the evolution at low redshift.

Pair-instability supernovae are very rare events in our population synthesis, and are also heavily biased towards very low metallicities and only a few stellar models in our grid, resulting in a predicted rate that is highly stochastic with redshift. We show the predicted rates smoothed over three time bins. The population peaks at in our model, but at a volume-averaged rate more than 1 dex below that of LGRBs.

The delay times for gravitational wave transients involving black holes and neutron stars are shorter than those for typical white dwarf mergers, but still substantial. As a result, the number density in predicted rate peaks between and for all three flavours of GW event considered here. Low metallicities (i.e. higher formation redshifts) favour the retention of mass during evolution of a stellar progenitor and thus the formation of a black hole, and this manifests as a slightly more rapid fall off in NS-BH events than NS-NS events towards .

Figure 3: The predicted rate density for electromagnetic and gravitational wave transients, as a function of source redshift and lookback time, given the volume-averaged history of cosmic star formation and chemical enrichment. We show both the total SN Ia rate and that due to the double degenerate pathway (DD) alone.

3.2.1 Observational Constraints

In figure 4 we compare the predictions with a compilation of observational constraints from a range of surveys and data sources. For observed SN Ia rates we use the data compilation of Melinder et al. (2012) and supplement this with rates more recently reported by Okumura et al. (2014), Cappellaro et al. (2015), Rodney et al. (2014) and Graur et al. (2014). For CCSNe we again use the compilation of Melinder et al. (2012) and additionally show the rates reported by Taylor et al. (2014), Strolger et al. (2015) and Petrushevska et al. (2016).

Broadly speaking the agreement between observational constraints and BPASS predictions as a function of redshift is good. For type Ia SNe, the predicted rates consistently tend towards the upper end of the range of observed values, although it agrees with a reasonable fraction of the observed data points, and with the observed redshift trend. A slight overestimate may arise because we have too loose a definition of what leads to a type Ia SN and some fraction of our progenitor models may in fact give rise to some other kind of low luminosity transient as discussed earlier. Type Ia rates are known to differ in cluster environments relative to the field, so a certain amount of scatter is expected. It is also, of course, possible that there are completeness corrections or other systematic offsets in some of the observational data, but overall the predictions, based purely on stellar population synthesis and the volume-averaged cosmic history are excellent.

The predictions also agree well with observed rates for core collapse supernovae as a function of redshift. The high redshift core-collapse supernova rates reported by Strolger et al (2015, red diamonds on figure 4) are systematically lower than other estimates, and than the BPASS predictions. These are intrinsically challenging deep field observations and the analysis from the CLASH and CANDELS fields necessarily surveys a smaller volume than most other supernova surveys. As a result, these rates may be subject to uncertainties due to cosmic variance. Nonetheless, if these lower CCSN rates are supported by future observations, the implication is that a smaller fraction of core collapse events are resulting in a luminous transient than expected - perhaps supporting the hypothesis that some black-hole forming events are dark (Adams et al., 2017).

The volumetric rate of long-GRBs is extremely hard to derive from existing data sets, given the sensitivity to the energy bands and trigger levels in the observed events, together with the extremely non-uniform follow-up and redshift determination for these sources. Further complication is added by the fact that these events are relativistically-beamed and any comparison to theoretical rates requires an estimate of the jet opening angle. In figure 4 we show the volumetric rates of events with an isotropic-equivalent energy E ergs reported by Perley et al. (2016) as a function of redshift, having corrected these to an event rate assuming a typical jet opening angle of 8 degrees (triangle) or 12 degrees (points).

Interestingly, despite considering only one plausible pathway for the formation of Long GRBs, and likely missing others, the BPASS rate predictions are higher than those determined by observations by a factor of even assuming a narrow jet. We note however that the Perley et al. (2016) study only estimated the rate for the most luminous GRBs. The energy range of observed bursts extends well below the cut-off used by the Perley et al study, and is parameterised by a broken power-law with slope -1.2 below a break luminosity of  ergs and -1.92 above this (Pescalli et al., 2016). Given this luminosity function, we would expect only a few percent of the total population of Long GRBs to be included in the Perley et al (2016) rates. To illustrate this, we also show the Perley et al rates adjusted to account for the luminosity function integrated over the range  ergs, but with a somewhat broader jet opening angle of 20 degrees (crosses). These exceed the rate estimates from BPASS by a factor of , confirming that there is scope for additional evolution pathways to be identified.

Figure 4: The predicted evolution in volume-averaged event rates predicted by BPASS (solid lines) and derived from a compilation of observational data for CCSNe, Type Ia SNe and LGRBs (see text for references). LGRB rates have been corrected to isotropic rates assuming jet opening angles of 8° (triangles) or 12° (points), and integrated over the luminosity function for a broader 20° opening angle (crosses).

For reference, we provide the inferred event rates at and , accounting for the global volume-averaged star formation and chemical enrichment history, in table 1. Our fiducial model is also available for download from the BPASS websites.

Transient Type log R() log R()
CCSNe 5.118 0.001 5.600 0.001
Ia 4.456 0.003 4.878 0.002
LGRB 2.47 0.02 3.07 0.01
PISN 1.01 0.12 1.63 0.06
NSWD 1.44 0.08 1.60 0.06
NSNS 2.61 0.02 2.89 0.02
NSBH 2.35 0.03 2.63 0.01
BHBH 1.80 0.05 1.96 0.04

Rates are given per type in units of events yr Gpc

Table 1: Local Universe transient event rates, given the volume-averaged history of cosmic star formation and chemical enrichment. Quoted uncertainties give Poisson uncertainties on observed predicted event rate based on our model grid, but do not account for statistical offsets or missing pathways.

Observational constraints on the local rate of NS-NS and BH-BH mergers are shown in figure 5. The majority of constraints on NS-NS mergers arise from the short gamma-ray burst population. As was the case for LGRBs, the conversion from an observed rate to an intrinsic rate density involves an uncertainty on the opening angle assumed for the relativistically-beamed emission. We show the volumetric rate density and its uncertainty assuming a jet opening angle of 10° and also the range of possible densities for angles in the range 3 to 26° (after Paul, 2018). This also encompasses the opening angle estimate of Fong et al. (2015, °) We also show the astrophysical rates inferred from the detection of gravitational wave transient GW 170817 Abbott et al. (2017b) and the handful of BH-BH mergers identified to date Abbott et al. (2017a, c, e).

As figures 4 and 5 demonstrate, the same stellar population synthesis model which successfully reproduces the rate of core collapse and type Ia supernovae also provides a good fit to the observed rate of binary NS and BH mergers.

No known constraints on the gravitational-wave driven merger of NS-WD binaries exists. These events have been proposed as one possible progenitor for the Ca-rich class of transients (Metzger, 2012). The merger rates predicted for NS-WD events by BPASS are substantially lower than current best estimates of the local volumetric rate of Ca-rich SNe (Frohmaier et al., 2018) suggesting that either this progenitor pathway is a negligible fraction of observed events, or that some alternate mechanism is required for hardening NS-WD binaries before the GW-dominated inspiral phase.

We finally note that our predicted rates are dependent on the assumed star-formation history and metallicity evolution of the Universe. For the star-formation history, the total star formation that occurred at the highest redshifts, beyond , is the most uncertain. This only has an impact on the rate of those events with the longest delay times, i.e. the GW events and double-degenerate type Ia SNe. The uncertainties in the cosmic metallicity evolution has the strongest effect on the events with the strongest dependence on metallicity, i.e. LGRBs and PISNe. Here their rate evolution depends strongly on how quickly the Universe became metal-rich. However, considering the good qualitative agreement between our current predictions and observational constraints there is no evidence that the fiducial model for this metallicity and star-formation history should be changed.

Figure 5: The predicted BPASS merger rates of NSNS and BHBH binaries at (horizontal solid line), and (dotted line), compared to local rate estimates of events in the literature. Points indicate constraints from short GRBs, while boxes indicate constraints from gravitational wave transient detections. For short GRBs, the extended grey bar indicates the range of number densities assuming jet opening angles ° while bold points indicate the number density and its uncertainty inferred at °, except in the case of ref. 12 where the uncertainty in opening angles is already accounted for in the error bars shown. Data sources: GW170817 - Abbott et al. (2017b), BHBH rate estimates - Abbott et al. (2017a), SGRBs: 1 - Ando (2004), 2 - Guetta & Piran (2006), 3 - Nakar et al. (2006), 4 - Guetta & Stella (2009), 5 - Dietz (2011), 6 - Coward et al. (2012), 7 - Siellez et al. (2014), 8 - Wanderman & Piran (2015), 9 - Yonetoku et al. (2004), 10 - Ghirlanda et al. (2016), 11 - Paul (2018), 12 - Fong et al. (2015).

4 Discussion

4.1 Implications for Gravitational Wave Observatories

For NS-NS merger transients, the mass range is relatively narrow. Nonetheless, the rate will be sensitive to the assumed metallicity history in a similar way; if more compact objects retain sufficient mass to become black holes, there may be a reduction in the relative rate of NS-NS mergers. Future observations of increased samples of short GRBs may constrain this through the occurrence rates of electromagnetic transients. However further complication in predicting gravitational wave transient rates is offered by the relatively small horizon of LIGO to these sources, which may be insufficient to average over the cosmic variance in the environs of the Local Group and nearby clusters. A more detailed analysis of large scale structure effects and the most probable host galaxies for such events will require a more detailed simulation informed by the patchy, merger driven star formation and metal enrichment histories and spreads on galactic scales. Combining BPASS event rates with an N-body based, semi-analytic galaxy evolution model is likely the most promising route to such simulations.

The precise rates of observed GW events, and the mass distribution of their progenitors, will remain subject to small number uncertainties for some time due to the low rate of detectable events. The current upgrades to both LIGO and VIRGO, and potentially the addition of KAGRA and LIGO-India, will allow fainter triggers to be confidently identified, and so there will likely be improvement in the number statistics and in the sensitivity horizon (and thus volume probed). The rates of GW events are expected to vary by less than a factor of two between and , and so the detailed redshift evolution (and hence metallicity dependence) will remain uncertain for some time to come, unless events can be localised to galaxies with well-constrained metallicity.

4.2 Future constraints

There are three key inputs to a physical analysis of transient rates: the intertwined cosmic star formation and metallicity histories, the stellar evolution models, and the observational selection effects.

The cosmic history of star formation and metal enrichment is being constantly refined, and independently determined using different samples and indicators. The rise of large area infrared surveys, including those obtained using the Herschel telescope, has allowed the hitherto hidden dust-obscured components of star formation to be added to those observed in visible and ultraviolet light. The history of galaxy formation is also being modelled through ever more precise N-body and hydrodynamical simulations. Uncertainties do remain. There is likely to be considerable variation on galactic scales, and in cluster environments. Nonetheless, the volume-averaged results are relatively robust and will be further refined by future observations with the James Webb Space Telescope (JWST) and other large scale facilities.

The corrections applied to determine volumetric rate estimates of transients from observed number counts are also moderately well understood. Current uncertainties on observational event rates are dominated by corrections for the electromagnetic transient luminosity function and, in the case of beamed sources, the jet opening angle. Survey sensitivity and completeness corrections can introduce systematic offsets, which suggests that a uniform survey is desirable, but these are now being modelled using sophisticated simulation and the scatter on rate estimates from different surveys has decreased over time.

The Large Synoptic Survey Telescope (LSST) is expected to begin observations in 2022. The 8 m class telescope is expected to identify and provide redshifts for of order one hundred thousand of type Ia supernova to through their photometric colours and light curves, with the primary goal of constraining dark energy. The number of CCSNe detected will likely be several times higher (LSST Science Collaboration et al., 2009, 2017; Zhan & Tyson, 2018). These will simultaneously provide a uniform and precise measurement of the volumetric density of these events which will allow the redshift evolution of the supernova population to be measured with a precision hitherto impossible. Between and , we expect the core collapse supernova rate to evolve by a factor of six, while the type Ia rate declines by a much smaller factor of two. Thus reducing their relative volumetric rate uncertainties to (requiring samples of 1000 events in each class per redshift bin) would allow testing of this theoretical scenario and any alternatives as necessary.

In samples this size, statistical corrections for observational selection will dominate rate uncertainties and here a single, uniform survey has obvious benefits. However the precision of LSST-derived rates will suffer from the limited plans for rapid spectroscopic follow up of these relatively faint transients. If a significant fraction of events is followed up spectroscopically, and classified, it will also be possible to contrast the volumetric rates in Type II (II-L, II-P, II-n), Type Ib and Type Ic supernovae - a sensitive probe of the mass loss and stripping of massive stars in binary systems. LSST’s deep, high-cadence supernova survey may also facilitate identification of rare transient sources such as PISNe, and may provide an alternate route to identify some compact mergers through kilonova emission, although this is strongly dependent on the final survey cadence selected.

Rate estimates for GRBs (both long and short) will likely remain uncertain for some time. Most of the current constraints arise from the Neil Gehrels Swift Observatory, due to its sensitivity and ability to localise a large fraction gamma-ray transients. The Fermi Gamma-Ray Space Telescope and INTErnational Gamma-Ray Astrophysics Laboratory (INTEGRAL) also continue to provide burst alerts, with different energy sensitivity ranges. While this makes them potentially more sensitive to short (hard) GRBs than Swift, both suffer from poor burst localisation and limited follow-up, such that very few bursts have redshifts. It is possible that the new generation of ground-based, wide field images designed for GW transient follow-up (such as the Gravitational Optical Transient Observatory, GOTO, Steeghs, 2017) will be able to localise a larger fraction of future events, but it is not clear whether these will be priorised.

The Space-based multi-band astronomical Variable Objects Monitor mission (SVOM, due to launch in 2021, Wei et al., 2016) will produce a GRB sample smaller than that accumulated by current Swift observations, but incorporates systematic rapid follow up, and aims to substantially improve number counts for distant () bursts, which are strongly sensitive to the metallicity enrichment history within the first Gyr of cosmic history. In this sense, it will be complementary to improved measurements of star formation rate variation and metal enrichment at these epochs with the JWST telescope. SVOM also aims to have increased sensitivity to the hard spectra of short bursts, their thermalised prompt emission and their time evolution, when compared to Swift. Hence, it is expected to provide constraints on NS-NS and potentially NS-BH mergers.

The more ambitious, and expensive, Theseus mission (currently being considered for ESA’s M5 mission to launch in 2032, Amati et al., 2018) will address the question in more detail. Theseus, while a longer term prospect, expects to detect GRBs at an event rate an order of magnitude higher than the existing Swift telescope, meaning that it will equal the precision of existing GRB rate estimates within one year, and far exceed it over the mission lifetime. It will also greatly improve constraints on the GRB luminosity function. while existing instruments such as Fermi and Integral are already probing events with energies outside of the peak sensitivity range of Swift.

So, given these sources of information, is there potential for transient rates to constrain the third input: the physics of stellar evolution models, particularly at the low metallicities and in the binary environments required to create the most massive transients?

Barrett et al. (2018) suggested that future GW event detections will be able to highly constrain the input physics of binary population synthesis codes. However they ignored the fact that we also must understand the metallicity and star-formation evolution of the Universe. Here by using the same model codes to predict the rate of SNe we have a second complimentary dataset that is mostly sensitive to the star-formation history of the Universe. The CCSN rate effectively directly measures the star-formation rate while the type Ia SN provides a measurement of older star-formation history if we can reproduce the delay-time distribution of these events. The fact here we can reproduce the observed rates and delay-time distributions for both SNe types means that if our predicted GW event rates do not match constraints from future observations then it must be our stellar models that are incorrect.

The relevant uncertain aspects of the stellar evolution models are the mass-loss rates, common envelope evolution and the neutron-star/black-hole kicks. These all have a minor effect on the overall supernova rate and primarily will change the SNe to be of different types observationally. However they all have an effect on the rate of GW events as they determine how close two compact remnants are when they explode, and thus how long they will take to merge via gravitational radiation. Only by such a multi-event type approach as we use here is it possible to start to test aspects of the stellar evolution physics.

BPASS is a detailed stellar evolution code in which each individual stellar model takes several minutes to run and models feed into the population synthesis analysis required to produce figure 3. The computing resources required to examine a grid of such models with varying stellar astrophysics is only today becoming available. However simulations based on the rapid population synthesis codes, e.g. BSE and COMPAS (which uses semi-analytic formalisms to estimate the evolution of stars with physical properties between fixed grid points), have already suggest that reasonable changes to the input stellar physics can cause up to an order of magnitude variation in compact binary merger rates and that stellar physics constraints are obtainable (e.g. Taylor & Gerosa, 2018; Giacobbo et al., 2018). Given the ability of binary stellar population synthesis codes, such as BPASS, to model multiple flavours of transients, the increased precision in rate density measurements which should be obtained from both LSST and LIGO may be sufficient to start to constrain different physical assumptions. This will be explored in future work using BPASS models.

5 Conclusions

Our main conclusions can be summarised as follows:

  1. We combine event rates predicted by BPASS population synthesis models as a function of age and metallicity with analytic prescriptions for the cosmic volume-averaged star formation rate and chemical evolution history.

  2. We determine the redshift dependence of transient event rates for both electromagnetic and gravitational wave transients from the same stellar population synthesis model.

  3. These are in excellent agreement with current observational constraints, and with the rate estimates determined from gravitational wave sources, albeit with uncertainties on the observed source luminosity functions and the opening angle of relativistically-beamed emission.

We aim to continue using BPASS to constrain the joint evolution of different astrophysical transients as data statistics improve. We also plan to refine our models. In particular, we aim to analyse the effect of variations in the analytic star formation and metallicity histories in this study, and also to replace volume averaged rates with metallicity spread and star formation environment as determined from cosmological models. Ultimately, we will use this to look at transient rates in cluster versus field environments and as a function of galaxy/cluster mass, as well as exploring the effects of stellar physics and multiplicity statistics.

acknowledgements

JJE and ERS acknowledge travel funding and support from the University of Auckland. ERS also acknowledges support from the Australian Astronomical Observatory through the Shaw Distinguished Visitor scheme. We thank Joe Lyman for useful conversations. BPASS is enabled by the resources of the NeSI Pan Cluster. New Zealand’s national facilities are provided by the NZ eScience Infrastructure and funded jointly by NeSI’s collaborator institutions and through the Ministry of Business, Innovation & Employment’s Research Infrastructure programme. URL: https://www.nesi.org.nz.

References

  • Abbott et al. (2016) Abbott B. P., et al., 2016, Physical Review Letters, 116, 061102
  • Abbott et al. (2017a) Abbott B. P., et al., 2017a, Physical Review Letters, 118, 221101
  • Abbott et al. (2017b) Abbott B. P., et al., 2017b, Phys. Rev. Lett., 119, 161101
  • Abbott et al. (2017c) Abbott B. P., et al., 2017c, Physical Review Letters, 119, 141101
  • Abbott et al. (2017d) Abbott B. P., et al., 2017d, ApJ, 848, L12
  • Abbott et al. (2017e) Abbott B. P., et al., 2017e, ApJ, 851, L35
  • Abell et al. (2009) Abell P. A., et al., 2009
  • Adams et al. (2017) Adams S. M., Kochanek C. S., Gerke J. R., Stanek K. Z., Dai X., 2017, MNRAS, 468, 4968
  • Amati et al. (2018) Amati L., et al., 2018, Advances in Space Research, 62, 191
  • Ando (2004) Ando S., 2004, J. Cosmology Astropart. Phys., 6, 007
  • Barrett et al. (2018) Barrett J. W., Gaebel S. M., Neijssel C. J., Vigna-Gómez A., Stevenson S., Berry C. P. L., Farr W. M., Mandel I., 2018, MNRAS, 477, 4685
  • Belczynski et al. (2008) Belczynski K., Kalogera V., Rasio F. A., Taam R. E., Zezas A., Bulik T., Maccarone T. J., Ivanova N., 2008, ApJS, 174, 223
  • Bray & Eldridge (2016) Bray J. C., Eldridge J. J., 2016, MNRAS, 461, 3747
  • Bray & Eldridge (2018) Bray J. C., Eldridge J. J., 2018, preprint, (arXiv:1804.04414)
  • Cappellaro et al. (2015) Cappellaro E., et al., 2015, A&A, 584, A62
  • Chambers et al. (2016) Chambers K. C., et al., 2016, preprint, (arXiv:1612.05560)
  • Coward et al. (2012) Coward D. M., et al., 2012, MNRAS, 425, 2668
  • De Donder & Vanbeveren (1998) De Donder E., Vanbeveren D., 1998, A&A, 333, 557
  • Dietz (2011) Dietz A., 2011, A&A, 529, A97
  • Eldridge & Stanway (2016) Eldridge J. J., Stanway E. R., 2016, MNRAS, 462, 3302
  • Eldridge & Tout (2004) Eldridge J. J., Tout C. A., 2004, MNRAS, 353, 87
  • Eldridge et al. (2011) Eldridge J. J., Langer N., Tout C. A., 2011, MNRAS, 414, 3501
  • Eldridge et al. (2013) Eldridge J. J., Fraser M., Smartt S. J., Maund J. R., Crockett R. M., 2013, MNRAS, 436, 774
  • Eldridge et al. (2017) Eldridge J. J., Stanway E. R., Xiao L., McClelland L. A. S., Taylor G., Ng M., Greis S. M. L., Bray J. C., 2017, Publ. Astron. Soc. Australia, 34, e058
  • Fong et al. (2015) Fong W., Berger E., Margutti R., Zauderer B. A., 2015, ApJ, 815, 102
  • Frohmaier et al. (2018) Frohmaier C., Sullivan M., Maguire K., Nugent P., 2018, ApJ, 858, 50
  • Ghirlanda et al. (2016) Ghirlanda G., et al., 2016, A&A, 594, A84
  • Giacobbo et al. (2018) Giacobbo N., Mapelli M., Spera M., 2018, MNRAS, 474, 2959
  • Graur et al. (2014) Graur O., et al., 2014, ApJ, 783, 28
  • Guetta & Piran (2006) Guetta D., Piran T., 2006, A&A, 453, 823
  • Guetta & Stella (2009) Guetta D., Stella L., 2009, A&A, 498, 329
  • Heger & Woosley (2002) Heger A., Woosley S. E., 2002, ApJ, 567, 532
  • Heger et al. (2003) Heger A., Fryer C. L., Woosley S. E., Langer N., Hartmann D. H., 2003, ApJ, 591, 288
  • Hobbs et al. (2005) Hobbs G., Lorimer D. R., Lyne A. G., Kramer M., 2005, MNRAS, 360, 974
  • Holoien et al. (2017) Holoien T. W.-S., et al., 2017, MNRAS, 464, 2672
  • Howell (2011) Howell D. A., 2011, Nature Communications, 2, 350
  • Kato & Hachisu (2004) Kato M., Hachisu I., 2004, ApJ, 613, L129
  • Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
  • LSST Science Collaboration et al. (2009) LSST Science Collaboration et al., 2009, preprint, (arXiv:0912.0201)
  • LSST Science Collaboration et al. (2017) LSST Science Collaboration et al., 2017, preprint, (arXiv:1708.04058)
  • Langer (2012) Langer N., 2012, ARA&A, 50, 107
  • Langer & Norman (2006) Langer N., Norman C. A., 2006, ApJ, 638, L63
  • Langer et al. (2000) Langer N., Deutschmann A., Wellstein S., Höflich P., 2000, A&A, 362, 1046
  • Liu et al. (2017) Liu D., Wang B., Wu C., Han Z., 2017, A&A, 606, A136
  • Madau & Dickinson (2014) Madau P., Dickinson M., 2014, ARA&A, 52, 415
  • Maoz & Badenes (2010) Maoz D., Badenes C., 2010, MNRAS, 407, 1314
  • Maoz et al. (2010) Maoz D., Sharon K., Gal-Yam A., 2010, ApJ, 722, 1879
  • Maoz et al. (2014) Maoz D., Mannucci F., Nelemans G., 2014, ARA&A, 52, 107
  • Melinder et al. (2012) Melinder J., et al., 2012, A&A, 545, A96
  • Mennekens & Vanbeveren (2016) Mennekens N., Vanbeveren D., 2016, A&A, 589, A64
  • Mennekens et al. (2010) Mennekens N., Vanbeveren D., De Greve J. P., De Donder E., 2010, A&A, 515, A89
  • Metzger (2012) Metzger B. D., 2012, MNRAS, 419, 827
  • Moe & Di Stefano (2017) Moe M., Di Stefano R., 2017, ApJS, 230, 15
  • Nakar et al. (2006) Nakar E., Gal-Yam A., Fox D. B., 2006, ApJ, 650, 281
  • Nomoto & Kondo (1991) Nomoto K., Kondo Y., 1991, ApJ, 367, L19
  • Okumura et al. (2014) Okumura J. E., et al., 2014, PASJ, 66, 49
  • Paczynski (1986) Paczynski B., 1986, ApJ, 308, L43
  • Paul (2018) Paul D., 2018, MNRAS, 477, 4275
  • Perley et al. (2016) Perley D. A., et al., 2016, ApJ, 817, 7
  • Pescalli et al. (2016) Pescalli A., et al., 2016, A&A, 587, A40
  • Peters (1964) Peters P. C., 1964, Physical Review, 136, 1224
  • Petrushevska et al. (2016) Petrushevska T., et al., 2016, A&A, 594, A54
  • Podsiadlowski et al. (1992) Podsiadlowski P., Joss P. C., Hsu J. J. L., 1992, ApJ, 391, 246
  • Raddi et al. (2018) Raddi R., Hollands M. A., Koester D., Gänsicke B. T., Gentile Fusillo N. P., Hermes J. J., Townsley D. M., 2018, ApJ, 858, 3
  • Rodney et al. (2014) Rodney S. A., et al., 2014, AJ, 148, 13
  • Ruiter et al. (2009) Ruiter A. J., Belczynski K., Fryer C., 2009, ApJ, 699, 2026
  • Ruiter et al. (2011) Ruiter A. J., Belczynski K., Sim S. A., Hillebrandt W., Fryer C. L., Fink M., Kromer M., 2011, MNRAS, 417, 408
  • Siellez et al. (2014) Siellez K., Boër M., Gendre B., 2014, MNRAS, 437, 649
  • Smartt (2015) Smartt S. J., 2015, Publ. Astron. Soc. Australia, 32, e016
  • Stanway & Eldridge (2018) Stanway E. R., Eldridge J. J., 2018, MNRAS,
  • Stanway et al. (2015) Stanway E. R., Levan A. J., Tanvir N., Wiersema K., van der Horst A., Mundell C. G., Guidorzi C., 2015, MNRAS, 446, 3911
  • Steeghs (2017) Steeghs D., 2017, Nature Astronomy, 1, 741
  • Strolger et al. (2015) Strolger L.-G., et al., 2015, ApJ, 813, 93
  • Tanvir et al. (2013) Tanvir N. R., Levan A. J., Fruchter A. S., Hjorth J., Hounsell R. A., Wiersema K., Tunnicliffe R. L., 2013, Nature, 500, 547
  • Taylor & Gerosa (2018) Taylor S. R., Gerosa D., 2018, preprint, (arXiv:1806.08365)
  • Taylor et al. (2014) Taylor M., et al., 2014, ApJ, 792, 135
  • Tonry et al. (2018) Tonry J. L., et al., 2018, PASP, 130, 064505
  • Totani et al. (2008) Totani T., Morokuma T., Oda T., Doi M., Yasuda N., 2008, PASJ, 60, 1327
  • Wanderman & Piran (2015) Wanderman D., Piran T., 2015, MNRAS, 448, 3026
  • Wang et al. (2017) Wang B., Podsiadlowski P., Han Z., 2017, MNRAS, 472, 1593
  • Wei et al. (2016) Wei J., et al., 2016, preprint, (arXiv:1610.06892)
  • Xiao & Eldridge (2015) Xiao L., Eldridge J. J., 2015, MNRAS, 452, 2597
  • Xiao et al. (2018) Xiao L., Stanway E., Eldridge J. J., 2018, preprint, (arXiv:1801.07068)
  • Yonetoku et al. (2004) Yonetoku D., Murakami T., Nakamura T., Yamazaki R., Inoue A. K., Ioka K., 2004, ApJ, 609, 935
  • Zapartas et al. (2017) Zapartas E., et al., 2017, A&A, 601, A29
  • Zhan & Tyson (2018) Zhan H., Tyson J. A., 2018, Reports on Progress in Physics, 81, 066901
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
290358
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description