Measuring the delay time distribution of binary neutron stars. III. Using the individual star formation histories of gravitational wave event host galaxies in the local universe

TL;DR

This study improves the measurement of binary neutron star delay time distributions by reconstructing individual host galaxy star formation histories, reducing the number of required host galaxies for accurate constraints, and preparing for future gravitational wave observations.

Contribution

It demonstrates that detailed SFH reconstruction from individual host galaxies enhances DTD constraints, requiring fewer hosts than previous scaling relation methods.

Findings

Reconstructed SFHs improve DTD parameter constraints.

Sample size of 100-300 hosts suffices for accurate DTD estimation.

Method anticipates data from upcoming GW detectors within a decade.

Abstract

In Paper I (Safarzadeh & Berger 2019) we studied the determination of the delay time distribution (DTD) of binary neutron stars (BNS) through scaling relations between halo/stellar mass and the star formation history (SFH) of galaxies hosting gravitational wave (GW) events in the local universe. Here we explore how a detailed reconstruction of the individual SFHs of BNS merger host galaxies can improve on the use of the scaling relations. We use galaxies from the Galaxy and Mass Assembly (GAMA) survey, which is mass complete at $M_*>10^9$ M$_\odot$ in the redshift range $0.05<z<0.08$. We use the reconstructed SFHs derived from the Prospector code, for two distinct sets of priors (favoring continuous and bursty SFHs), and convolve those with power law DTDs characterized by an index $\Gamma$ and a minimum delay time $t_{\rm min}$. We find that with this approach $\mathcal{O}(100)-\mathcal{O}(300)$ host galaxies are required to constrain the DTD parameters, with the number depending on the choice of SFH prior and on the parameters of the true DTD. We further show that using only the host galaxies of BNS mergers, as opposed to the full population of potential host galaxies in the relevant cosmic volume, leads to a minor bias in the recovered DTD parameters. The required host galaxy sample size is nearly an order of magnitude smaller relative to the approach of using scaling relations, and we expect such a host galaxy sample to be collected within a decade or two, prior to the advent of third-generation GW detectors.

Figures

Equations

$${\dot{n}}_i=$$

$$P(N_i|\Gamma ,t_{\min },{\psi }_i(z))=\mathrm{Poisson}({\dot{n}}_i\Delta t)=\frac{({\dot{n}}_i\Delta t{)}^{N_i}{e}^{-{\dot{n}}_i\Delta t}}{N_i!}$$

$$N_i\sim \mathrm{Poisson}({\dot{n}}_i\Delta t)$$

$$P(\{N_i\}|\Gamma ,t_{\min },\{{\psi }_i(z)\})=\prod _{i=1}^{N_{\mathrm{gal}}}P(N_i|\Gamma ,t_{\min },{\psi }_i(z)).$$

$$P(N_i|\Gamma ,t_{\min })=\int P(N_i|\Gamma ,t_{\min },{\psi }_i(z))P({\psi }_i(z))d[{\psi }_i(z)].$$

$$P(N_i|\Gamma ,t_{\min })\approx \frac{1}{N_{\mathrm{samp}}}\sum _{j=1}^{N_{\mathrm{samp}}}\frac{({\dot{n}}_{i,j}\Delta t{)}^{N_i}{e}^{-{\dot{n}}_{i,j}\Delta t}}{N_i!}$$

$$P(\Gamma ,t_{\min }|\{N_i\})$$

$$\propto \prod _{i=1}^{N_{\mathrm{gal}}}\sum _{j=1}^{N_{\mathrm{samp}}}\frac{({\dot{n}}_{i,j}\Delta t{)}^{N_i}{e}^{-{\dot{n}}_{i,j}\Delta t}}{N_i!}$$

$$P(\Gamma ,t_{\min }|N_{\mathrm{BNS}})$$

$$=\int P(\Gamma ,t_{\min }|\{N_i\},N_{\mathrm{BNS}})P(\{N_i\}|N_{\mathrm{BNS}})d(\{N_i\})$$

$$P(\Gamma ,t_{\min }|N_{\mathrm{BNS}})$$

$$\propto \sum _{k=1}^{N_{\mathrm{repeat}}}\prod _{i=1}^{N_{\mathrm{gal}}}\sum _{j=1}^{N_{\mathrm{samp}}}\frac{({\dot{n}}_{i,j}\Delta t{)}^{N_{i,k}}{e}^{-{\dot{n}}_{i,j}\Delta t}}{N_{i,k}!}$$

$$P(\Gamma ,t_{\min }|N_{\mathrm{BNS}})\sim \sum _{j=1}^{N_{\mathrm{repeat}}}\prod _{i=1}^{N_{\mathrm{gal}}}\frac{({\dot{n}}_{i,j}\Delta t{)}^{N_{i,j}}{e}^{-{\dot{n}}_{i,j}\Delta t}}{N_{i,j}!}$$

$$\tilde{P}(\Gamma ,$$

Full text

1 Introduction

The delay time distribution (DTD) of binary neutron stars (BNS) is currently poorly constrained, but as we have recently shown it can be determined using both the mass distribution of BNS merger host galaxies in the local universe (Safarzadeh & Berger 2019; hereafter, Paper I) and the redshift distribution of BNS mergers as probed by third-generation gravitational wave (GW) detectors (Safarzadeh et al. 2019; hereafter, Paper II). The former approach takes advantage of galaxy scaling relations that map halo/stellar mass into star formation history (SFH), which when convolved with the DTD lead to a predicted BNS merger host galaxy mass function (this approach was previously proposed and used in the context of short GRBs: Zheng & Ramirez-Ruiz 2007; Kelley et al. 2010; Leibler & Berger 2010; Fong et al. 2013; Behroozi et al. 2014). The host galaxies of BNS mergers can be identified through the detection of electromagnetic (EM) counterparts, but this is likely only achievable within a few hundred Mpc. We found that for a power law DTD characterized by index $\Gamma$ and minimum delay $t_{\rm min}$, $\mathcal{O}(10^{3})$ host galaxies are required to reasonably constrain the DTD.

The latter approach instead relies on a redshift mapping of the BNS merger rate, which requires GW detections to $z\sim{\rm few}$, achievable with the next-generation Einstein Telescope (ET) and Cosmic Explorer (CE). In this approach, it is unlikely that EM counterparts can be detected, but the individual redshift uncertainties from the GW data ($\delta z/z\approx 0.1z$) can be overcome through a large number of anticipated detections, $\sim 10^{5}$ yr*-1*. We found that with about a year of CE+ET data the DTD parameters, as well as the mass efficiency of BNS production, can be determined to about 10%.

Here we continue our investigation of the DTD, with an alternative approach to the use of BNS merger host galaxies at $z\approx 0$. Namely, unlike in Paper I, which used scaling relations between mass and SFH, we explore the use of detailed reconstructed SFHs for the individual host galaxies. In §2 we present the galaxy sample used for this study (the Galaxy and Mass Assembly survey) and its SFH reconstruction. In §3 we present the method for extracting the DTD from the galaxy SFHs, as well as our approach to evaluating the number of host galaxies required. In §4 we discuss our findings in terms of the sample size needed, and we conclude in §5. We adopt the Planck 2015 cosmological parameters (Planck Collaboration et al., 2016): $\Omega_{M}=0.308$, $\Omega_{\Lambda}=0.692$, $\Omega_{b}=0.048$, and $H_{0}=67.8$ km s*-1* Mpc*-1*.

2 Galaxy Data

We use SFHs inferred from galaxy photometry in Leja et al. (2019). The photometry is measured with the LAMBDAR code (Wright et al., 2016) from DR3 of the Galaxy and Mass Assembly (GAMA) survey (Driver et al., 2011; Baldry et al., 2018), and includes 21 bands ranging from the far-UV (GALEX) to the far-IR (Herschel). The galaxies have spectroscopic redshifts. All galaxies at $0.05<z<0.08$ with stellar masses $M_{*}>10^{9}$ M*⊙* as determined in Taylor et al. (2011) are modeled, resulting in a mass-complete sample of 6134 galaxies.

The photometry was fit using the Prospector-$\alpha$ model within the Prospector inference machine (Leja et al., 2017; Johnson & Leja, 2017). This model includes a seven-parameter non-parametric star formation history, as well as flexible dust attenuation model, far-infrared dust re-emission via energy balance, and nebular emission self-consistently powered by the stellar ionizing continuum.

A key source of uncertainty in SFH recovery is the choice of prior (Carnall et al., 2019; Leja et al., 2019). This sensitivity occurs because the SEDs of stellar populations change slowly as a function of time; specifically, SEDs evolve roughly evenly in each logarithmic time step (Ocvirk et al., 2006). As a result, the SFR inferred in adjacent time bins is typically highly degenerate. Fortunately, key outputs of galaxy SED-fitting such as the mass-to-light ratio and, to a lesser extent, the recent SFR, are calculated using moments of the SFH, which are largely insensitive to this degeneracy (Bell et al., 2003; Leja et al., 2019). However, the DTD does not interact with these conserved quantities but instead couples directly to the SFH, and is thus sensitive to this degeneracy. Accordingly, we perform the analysis here using two different priors that assume opposite behaviors in this degeneracy: one that favors bursty SFHs (hereafter, logM) and one that favors smooth SFHs (hereafter, continuity) (Leja et al., 2019). Two representative examples for galaxies with opposite SFHs are shown in Figure 1, for both the logM and continuity priors. These examples highlight both the range of behavior and the associated uncertainties (statistical and systematic) in the reconstructed SFH.

The two SFH priors adopted in this work encapsulate the plausible range of choices (Leja et al., 2019). Star formation histories derived from high-resolution optical/near-IR spectroscopy can provide more precise SFHs than photometry alone, although in most cases this will be a modest improvement (Pacifici et al., 2012).

3 DTD Determination Method

The expected BNS merger rate at $z=0$ for galaxy $i$ with star formation history $\psi_{i}(z)$ is given by:

[TABLE]

where $dt/dz=-[(1+z)E(z)H_{0}]^{-1}$ and $E(z)=\sqrt{{\Omega}_{m,0}(1+z)^{3}+{\Omega}_{k,0}(1+z)^{2}+{\Omega}_{\Lambda}(z)}$; $t_{b}$ is the cosmic time corresponding to $z_{b}$; $\lambda$ is the BNS production mass efficiency, assumed to be a fixed value of $10^{-5}$ M${}_{\odot}^{-1}$, independent of redshift or environment; and $dP_{m}/dt$ is the DTD, likewise assumed to be independent of environment. As in Papers I and II, we parameterize the DTD to follow a power law111In this work (as well as in Papers I and II) we have focused on a power law DTD, which is well motivated. However, it is possible that the true DTD might deviate from this power law form. This could potentially be investigated by comparing the results of the analysis proposed here (using host galaxies at $z\sim 0$) and the analysis in Paper II (using the redshift distribution of BNS mergers from third-generation detectors). with index $\Gamma$, minimum delay time $t_{\rm min}$, and a fixed maximum delay time $t_{\rm max}=10$ Gyr. Our results are not sensitive for a larger value of $t_{\rm max}$. This formulation of Equation 1 is identical to that used in Paper I, except that there we took $\psi$ to be a direct function of halo mass, $M_{h}$.

To generate the simulated data, we assume a fixed number of $N_{\rm gal}=1000$ galaxies that can serve as possible hosts for BNS merger events. These are selected to serve as a representative subset of the GAMA sample that preserves the mass distribution of the full sample. For a given DTD ($\Gamma$ and $t_{\rm min}$) and SFH ($\psi_{i}(z)$) we can then estimate the mean merger rate, $\dot{n}_{i}$, for each galaxy using Equation 1. This then gives each galaxy a different probability to host a BNS merger event for the different DTDs (see Figure 1). As in Papers I and II, we use a set of 9 representative DTDs, with $\Gamma=[-1.5,-1,-0.5]$ and $t_{\rm min}=[10,100,1000]$ Myr.

The number of BNS merger events, $N_{i}=\{0,1,\dots\}$, observed for any given galaxy over a given period of time, $\Delta t$, follows a Poisson distribution based on the merger rate:

[TABLE]

We then can simulate $N_{i}$ directly by drawing it from the Poisson distribution based on the associated rate:

[TABLE]

Once we have simulated a set of BNS merger events, we determine the constraining power they have on the underlying DTD. Assuming the BNS merger events in each galaxy are independent of each other and the SFHs are known, the corresponding likelihood is

[TABLE]

We further need to marginalize over the uncertainty on the SFH of each galaxy:

[TABLE]

We can approximate this integral by averaging over $N_{\rm samp}$ of the samples from the SFH posteriors for each galaxy:

[TABLE]

The resulting SFH-marginalized posterior of the DTD is therefore given by:

[TABLE]

where we have assumed that the prior over $\Gamma$ and $t_{\rm min}$ is uniform such that $P(\Gamma,t_{\rm min})$ is a constant.

We are interested in marginalizing over any particular set of events $\{N_{i}\}$ associated with the $N_{\rm gal}$ possible host galaxies to forecast possible constraints on the DTD as a function of the total number of BNS merger events, $N_{\rm BNS}$. This gives:

[TABLE]

We can approximate this integral using $N_{\rm repeat}$ realizations of the observed BNS merger event counts, conditioned on the total number of events $\sum_{i}N_{i}$ being equal to $N_{\rm BNS}$. Combining this with our previous expression then gives:

[TABLE]

where again $\sum_{i=1}^{N_{\rm gal}}N_{i,k}=N_{\rm BNS}$ for each realization.

Since the rate per galaxy is generally very small, $\dot{n}_{i}\ll 1$, we expect the uncertainty to be dominated by variation in the observed counts over the potential $N_{\rm gal}$ host galaxies rather than the uncertainties in each galaxy’s SFH. As such, we opt to use the same SFHs used to generate the data to evaluate the posterior rather than trying to marginalize over them directly:

[TABLE]

It is critical to note that the posterior above is taken over the observed counts of all potential galaxy hosts, including those with $N_{i}=0$ for which no BNS merger events have been detected. That is because these “non-detections” in aggregate contain a non-negligible amount of information in a regime where $\dot{n}_{i}\ll 1$ and merger events are rare. To illustrate this, we also compare the “complete” posterior distribution $P(\Gamma,t_{\rm min}|N_{\rm BNS})$ derived above with the biased posterior distribution, $\tilde{P}(\Gamma,t_{\rm min}|N_{\rm BNS})$, ignoring the non-detections:

[TABLE]

where $\mathcal{I}(N_{i,j}>0)$ is the indicator function that evaluates to $1$ if the condition $N_{i,j}>0$ is true and [math] otherwise. In general, we expect that ignoring non-detections will bias the inferred DTD, which we discuss in the next section.

We compute the above posterior using $N_{\rm gal}=1000$ galaxies and $N_{\rm repeat}=100$ realizations for the 9 different assumed DTDs, and interpolate between their associated $\Gamma$ and $t_{\rm min}$ parameters to obtain the probability for a different pair of $\Gamma$ and $t_{\rm min}$ values.

To summarize, our inference procedure is as follows:

We select a star formation history, $\psi_{i}(z)$, for every galaxy $i$ from its SFH posterior distribution and compute the corresponding BNS merger rate, $\dot{n}_{i}$, using Equation 1. 2. 2.

We then sample the number of BNS merger events, $N_{i}$, from the corresponding Poisson distribution based on Equation 2 for a given timescale $\Delta t$. We repeat this process until the total number of events is $N_{\rm BNS}=\sum_{i}N_{i}$. 3. 3.

We repeat this procedure $N_{\rm repeat}$ times to generate many realizations of the BNS merger events, $\{N_{i}\}$, for varying SFHs. 4. 4.

We then use the simulated BNS merger events and SFHs from these realizations to compute the DTD posteriors including and excluding galaxies that did not host detected BNS merger events using Equations 12 and 13, respectively.

4 Results

In Figure 1 we show the BNS merger rate probability distribution function (PDF) of two galaxies from the GAMA survey, with SFHs chosen to peak at early and late cosmic time (right panels), for the 9 different DTD models; the vertical lines in each panel indicate the value of $\dot{n}$ for the mean SFH. We show the SFHs, and merger rate PDFs for both sets of priors (middle panel: continuity; right panel: logM). The SFHs from different assumed priors often do not overlap, illustrating the challenge of inferring SFHs from photometry. In particular, the galaxy in the lower panel illustrates a specific degeneracy where the photometry of star forming galaxies can often be fit equally well with a continuous SFH or with a large burst at $\lesssim 300$ Myr followed by a sharp drop in the star formation rate (Leja et al., 2019). The figure also clearly indicates how the convolution of the DTDs with different SFHs leads to distinct merger rate PDFs and how the choice of SFH prior affects the resulting merger rate PDFs.

In Figure 2 we show the constraints achieved on the parameters of the DTD model as the number of observed host galaxies increases from 30 to 100 to 300. We show the results for three different injected DTD models and for both the logM and continuity SFH priors. The contours mark the 90% confidence region. There are several key takeaway points from Figure 2. First, in the case that the true DTD prefers short merger timescales, by having a steep power law slope and short $t_{\rm min}$ (upper row of Figure 2), the DTD parameters can be reasonably constrained with fewer than $\mathcal{O}(100)$ host galaxies. Second, in other permutations of the DTD, $\mathcal{O}(300)$ host galaxies may be required to constrain the DTD parameters, but with a lingering degeneracy between $\Gamma$ and $t_{\rm min}$. Third, the logM prior leads to tighter constraints on the DTD parameters compared to the continuity prior because it allows for bursty SFHs that pick more well-defined timescales when convolved with the DTD; the continuity prior smooths the SFH and hence systematically reduces its constraining power.

We note that the results in Figure 2 assume knowledge of the SFHs of all potential host galaxies in the cosmic volume that contains the BNS merger events. In the case of Advanced LIGO/Virgo at design sensitivity, the luminosity distance range for BNS merger detection is about 200 Mpc, while for the planned A+ and Voyager upgrades this distance is expected to be at least twice as large. Thus, even at Advanced LIGO/Virgo sensitivity, this requires knowledge of the SFHs of $\sim 10^{6}$ galaxies, while for A+/Voyager this number increases to $\gtrsim 10^{7}$ galaxies. In Figure 3 we demonstrate the effect of neglecting the galaxies that did not host BNS mergers, using instead the SFHs of only the actual host galaxies. As expected, the resulting reconstructed DTD parameter distribution is biased with respect to the input model. However, this bias is not severe, and the resulting degenerate range of $\Gamma$ and $t_{\rm min}$ contains the “true” answer. This bias may be acceptable given the need to model the SFHs of only a few hundred galaxies as opposed to $\gtrsim 10^{6}$ galaxies.

5 Summary and Conclusions

As we have argued in Papers I and II, the DTD of BNS systems can be constrained in two primary ways using GW events: (i) using the properties of BNS merger host galaxies in the local universe, identified via an associated electromagnetic counterpart (Paper I and here); and (ii) using the BNS merger rate as a function of redshift, which requires third-generation GW detectors (Paper II). Here we expand on the method of Paper I, in which we used galaxy scaling relations to relate the mass function of BNS merger host galaxies to the parameters of the DTD. In particular, we explore the use of the actual SFHs of individual BNS merger host galaxies.

We find that the SFH reconstruction method improves on the use of scaling relations, reducing the required sample size by a factor of $\sim 3-10$, to $\mathcal{O}(100)-\mathcal{O}(300)$. The exact level of improvement depends on the choice of SFH prior, as well as on the location of the true DTD in the $\Gamma-t_{\rm min}$ parameter space. We further note that accurate reconstruction of the DTD requires knowledge of the SFHs of not only the actual host galaxies but also of the general galaxy population within the relevant cosmic volume ($\sim 3\times 10^{7}$ Mpc3 for Advanced LIGO and an order of magnitude larger for A+). However, while using the SFHs of only the host galaxies results in a bias, we find that this bias is not severe.

With the currently allowed range of the BNS merger rate, $110-3840$ Gpc*-3* yr*-1* (The LIGO Scientific Collaboration et al., 2018), it may take a decade or longer to collect a sample of $\mathcal{O}(100)-\mathcal{O}(300)$ host galaxies. However, there is a reasonable chance that such a sample will be available prior to the advent of third-generation GW detectors. This will allow for independent determinations of the DTD from the properties of host galaxies in the local universe (Papers I and III) and from the full merger rate redshift distribution (Paper II).

We are thankful to Enrico Ramirez-Ruiz, Evan Scannapieco, and Doug Finkbeiner for helpful discussions. This work was supported by the National Science Foundation under grant AST14-07835 and by NASA under theory grant NNX15AK82G. The Berger Time-Domain Group at Harvard is supported in part by NSF under grant AST-1714498 and by NASA under grant NNX15AE50G. J.L. is supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-1701487. MTS is thankful to the Center for Astrophysics | Harvard & Smithsonian for hospitality, which made this work possible.