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

The delay time distribution (DTD) of binary neutron stars (BNSs) 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 gamma-ray bursts (SGRBs): 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 Γ and minimum delay t_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 ∼ 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 (δz/z ≈ 0.1z) can be overcome through a large number of anticipated detections, ∼10⁵ yr⁻¹. 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 ≈ 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 Section 2 we present the galaxy sample used for this study (the Galaxy and Mass Assembly (GAMA) survey) and its SFH reconstruction. In Section 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 Section 4 we discuss our findings in terms of the sample size needed, and we conclude in Section 5. We adopt the Planck 2015 cosmological parameters (Planck Collaboration et al. 2016): Ω_M = 0.308, Ω_Λ = 0.692, Ω_b = 0.048, and H₀ = 67.8 km s⁻¹ Mpc⁻¹.

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 the third data release (DR3) of the GAMA survey (Driver et al. 2011; Baldry et al. 2018), and includes 21 bands ranging from the far-ultraviolet (GALEX) to the far-infrared (Herschel). The galaxies have spectroscopic redshifts. All galaxies at 0.05 < z < 0.08 with stellar masses M_* > 10⁹ 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-α model within the Prospector inference machine (Johnson & Leja 2017; Leja et al. 2017). This model includes a seven-parameter non-parametric SFH, 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 spectral energy distributions (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 star formation rate (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, log M) 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 log M and continuity priors. These examples highlight both the range of behavior and the associated uncertainties (statistical and systematic) in the reconstructed SFH.

Zoom In Zoom Out Reset image size

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

The expected BNS merger rate at z = 0 for galaxy i with SFH ψ_i(z) is given by

$$\begin{eqnarray}&&{\dot{n}}_{i}={\int }_{{z}_{b}=10}^{{z}_{b}=0}\lambda \displaystyle \frac{{{dP}}_{m}}{{dt}}(t-{t}_{b}-{t}_{\min }){\psi }_{i}({z}_{b})\displaystyle \frac{{dt}}{{dz}}({z}_{b}){{dz}}_{b},\end{eqnarray} \tag{ 1 }$$

where dt/dz = −[(1 + z) E(z) H₀]⁻¹ and $E(z)=\sqrt{{{\rm{\Omega }}}_{m,0}{\left(1+z\right)}^{3}+{{\rm{\Omega }}}_{k,0}{\left(1+z\right)}^{2}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}(z)};$ t_b is the cosmic time corresponding to z_b; λ 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 law³ with index Γ, minimum delay time t_min, and a fixed maximum delay time t_max = 10 Gyr. Our results are not sensitive for a larger value of t_max. This formulation of Equation (1) is identical to that used in Paper I, except that there we took ψ to be a direct function of halo mass, M_h.

To generate the simulated data, we assume a fixed number of ${N}_{\mathrm{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 (Γ and t_min) and SFH (ψ_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 nine representative DTDs, with Γ = [−1.5, −1, −0.5] and t_min = [10, 100, 1000] Myr.

The number of BNS merger events, N_i = {0, 1, ...}, observed for any given galaxy over a given period of time, Δt, follows a Poisson distribution based on the merger rate

$$\begin{eqnarray}&&P({N}_{i}| {\rm{\Gamma }},{t}_{\min },{\psi }_{i}(z))=\mathrm{Poisson}({\dot{n}}_{i}{\rm{\Delta }}t)=\displaystyle \frac{{\left({\dot{n}}_{i}{\rm{\Delta }}t\right)}^{{N}_{i}}{e}^{-{\dot{n}}_{i}{\rm{\Delta }}t}}{{N}_{i}!}.\end{eqnarray} \tag{ 2 }$$

We then can simulate N_i directly by drawing it from the Poisson distribution based on the associated rate

$$\begin{eqnarray}&&{N}_{i}\sim \mathrm{Poisson}({\dot{n}}_{i}{\rm{\Delta }}t).\end{eqnarray} \tag{ 3 }$$

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

$$\begin{eqnarray}&&P(\{{N}_{i}\}| {\rm{\Gamma }},{t}_{\min },\{{\psi }_{i}(z)\})=\displaystyle \prod _{i=1}^{{N}_{\mathrm{gal}}}P({N}_{i}| {\rm{\Gamma }},{t}_{\min },{\psi }_{i}(z)).\end{eqnarray} \tag{ 4 }$$

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

$$\begin{eqnarray}&&P({N}_{i}| {\rm{\Gamma }},{t}_{\min })=\int P({N}_{i}| {\rm{\Gamma }},{t}_{\min },{\psi }_{i}(z))P({\psi }_{i}(z))d[{\psi }_{i}(z)].\end{eqnarray} \tag{ 5 }$$

We can approximate this integral by averaging over N_samp of the samples from the SFH posteriors for each galaxy:

$$\begin{eqnarray}&&P({N}_{i}| {\rm{\Gamma }},{t}_{\min })\approx \displaystyle \frac{1}{{N}_{\mathrm{samp}}}\displaystyle \sum _{j=1}^{{N}_{\mathrm{samp}}}\displaystyle \frac{{\left({\dot{n}}_{i,j}{\rm{\Delta }}t\right)}^{{N}_{i}}{e}^{-{\dot{n}}_{i,j}{\rm{\Delta }}t}}{{N}_{i}!}.\end{eqnarray} \tag{ 6 }$$

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

$$\begin{eqnarray}&&P({\rm{\Gamma }},{t}_{\min }| \{{N}_{i}\})\propto P(\{{N}_{i}\}| {\rm{\Gamma }},{t}_{\min }| )P({\rm{\Gamma }},{t}_{\min })\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\propto \displaystyle \prod _{i=1}^{{N}_{\mathrm{gal}}}\displaystyle \sum _{j=1}^{{N}_{\mathrm{samp}}}\displaystyle \frac{{\left({\dot{n}}_{i,j}{\rm{\Delta }}t\right)}^{{N}_{i}}{e}^{-{\dot{n}}_{i,j}{\rm{\Delta }}t}}{{N}_{i}!}\end{eqnarray} \tag{ 8 }$$

where we have assumed that the prior over Γ and t_min is uniform such that P(Γ, t_min) is a constant.

We are interested in marginalizing over any particular set of events {N_i} associated with the N_gal possible host galaxies to forecast possible constraints on the DTD as a function of the total number of BNS merger events, N_BNS. This gives

$$\begin{eqnarray}&&\begin{array}{l}P({\rm{\Gamma }},{t}_{\min }| {N}_{\mathrm{BNS}})\\ \quad =\int P({\rm{\Gamma }},{t}_{\min }| \{{N}_{i}\},{N}_{\mathrm{BNS}})P(\{{N}_{i}\}| {N}_{\mathrm{BNS}})d(\{{N}_{i}\}).\end{array}\end{eqnarray} \tag{ 9 }$$

We can approximate this integral using N_repeat realizations of the observed BNS merger event counts, conditioned on the total number of events ${\sum }_{i}{N}_{i}$ being equal to N_BNS. Combining this with our previous expression then gives

$$\begin{eqnarray}&&P({\rm{\Gamma }},{t}_{\min }| {N}_{\mathrm{BNS}})\approx \displaystyle \frac{1}{{N}_{\mathrm{repeat}}}\displaystyle \sum _{k=1}^{{N}_{\mathrm{repeat}}}P({\rm{\Gamma }},{t}_{\min }| {\{{N}_{i}\}}_{k})\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\propto \displaystyle \sum _{k=1}^{{N}_{\mathrm{repeat}}}\displaystyle \prod _{i=1}^{{N}_{\mathrm{gal}}}\displaystyle \sum _{j=1}^{{N}_{\mathrm{samp}}}\displaystyle \frac{{\left({\dot{n}}_{i,j}{\rm{\Delta }}t\right)}^{{N}_{i,k}}{e}^{-{\dot{n}}_{i,j}{\rm{\Delta }}t}}{{N}_{i,k}!}\end{eqnarray} \tag{ 11 }$$

where again ${\sum }_{i=1}^{{N}_{\mathrm{gal}}}{N}_{i,k}={N}_{\mathrm{BNS}}$ for each realization.

Because 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_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:

$$\begin{eqnarray}&&P({\rm{\Gamma }},{t}_{\min }| {N}_{\mathrm{BNS}})\sim \displaystyle \sum _{j=1}^{{N}_{\mathrm{repeat}}}\displaystyle \prod _{i=1}^{{N}_{\mathrm{gal}}}\displaystyle \frac{{\left({\dot{n}}_{i,j}{\rm{\Delta }}t\right)}^{{N}_{i,j}}{e}^{-{\dot{n}}_{i,j}{\rm{\Delta }}t}}{{N}_{i,j}!}.\end{eqnarray} \tag{ 12 }$$

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({\rm{\Gamma }},{t}_{\min }| {N}_{\mathrm{BNS}})$ derived above with the biased posterior distribution, $\tilde{P}({\rm{\Gamma }},{t}_{\min }| {N}_{\mathrm{BNS}})$, ignoring the non-detections:

$$\begin{eqnarray}&&\tilde{P}({\rm{\Gamma }},{t}_{\min }| {N}_{\mathrm{BNS}})\sim \displaystyle \sum _{j=1}^{{N}_{\mathrm{repeat}}}\displaystyle \prod _{i=1}^{{N}_{\mathrm{gal}}}{ \mathcal I }({N}_{i,k}\gt 0)\displaystyle \frac{{\left({\dot{n}}_{i,j}{\rm{\Delta }}t\right)}^{{N}_{i,j}}{e}^{-{\dot{n}}_{i,j}{\rm{\Delta }}t}}{{N}_{i,j}!},\end{eqnarray} \tag{ 13 }$$

where ${ \mathcal I }({N}_{i,j}\gt 0)$ is the indicator function that evaluates to 1 if the condition N_i,j > 0 is true and 0 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_gal = 1000 galaxies and N_repeat = 100 realizations for the nine different assumed DTDs, and interpolate between their associated Γ and t_min parameters to obtain the probability for a different pair of Γ and t_min values.

To summarize, our inference procedure is as follows.

• 1.  
  We select a SFH, ψ_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.  
  We then sample the number of BNS merger events, N_i, from the corresponding Poisson distribution based on Equation (2) for a given timescale ${\rm{\Delta }}t$. We repeat this process until the total number of events is ${N}_{\mathrm{BNS}}={\sum }_{i}{N}_{i}$.
• 3.  
  We repeat this procedure N_repeat times to generate many realizations of the BNS merger events, {N_i}, for varying SFHs.
• 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.

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 nine 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: log M). 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 ≲300 Myr, followed by a sharp drop in the SFR (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 log M 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_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 Γ and t_min. Third, the log M 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.

Zoom In Zoom Out Reset image size

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 the Advanced Laser Interferometer Gravitational-Wave Observatory (adLIGO)/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 AdLIGO/Virgo sensitivity, this requires knowledge of the SFHs of ∼10⁶ galaxies, while for A+/Voyager this number increases to ≳10⁷ 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 Γ and t_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 ≳10⁶ galaxies.

Zoom In Zoom Out Reset image size

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 EM 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 ∼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 Γ − t_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 (∼3 × 10⁷ Mpc³ for AdLIGO 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⁻³ yr⁻¹ (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. M.T.S. is thankful to the Center for Astrophysics ∣ Harvard & Smithsonian for hospitality, which made this work possible.