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
Mohammadtaher Safarzadeh, Edo Berger, Joel Leja, Joshua S. Speagle
DDraft version May 29, 2019
Preprint typeset using L A TEX style AASTeX6 v. 1.0
MEASURING THE DELAY TIME DISTRIBUTION OF BINARY NEUTRON STARS. III. USING THEINDIVIDUAL STAR FORMATION HISTORIES OF GRAVITATIONAL WAVE EVENT HOST GALAXIES INTHE LOCAL UNIVERSE
Mohammadtaher Safarzadeh , Edo Berger , Joel Leja and Joshua S. Speagle Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA, 02138, USA, [email protected] School of Earth and Space Exploration, Arizona State University, AZ, USA
ABSTRACTIn 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 thestar 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 galaxiescan 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 ∗ > M (cid:12) in the redshift range . < z < . . Weuse the reconstructed SFHs derived from the Prospector code, for two distinct sets of priors (favoringcontinuous and bursty SFHs), and convolve those with power law DTDs characterized by an index Γ and a minimum delay time t min . We find that with this approach O( ) − O( ) host galaxies arerequired to constrain the DTD parameters, with the number depending on the choice of SFH priorand on the parameters of the true DTD. We further show that using only the host galaxies of BNSmergers, 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 nearlyan order of magnitude smaller relative to the approach of using scaling relations, and we expect sucha host galaxy sample to be collected within a decade or two, prior to the advent of third-generationGW detectors. INTRODUCTIONThe delay time distribution (DTD) of binary neutronstars (BNS) is currently poorly constrained, but as wehave recently shown it can be determined using boththe mass distribution of BNS merger host galaxies inthe local universe (Safarzadeh & Berger 2019; hereafter,Paper I) and the redshift distribution of BNS mergersas probed by third-generation gravitational wave (GW)detectors (Safarzadeh et al. 2019; hereafter, Paper II).The former approach takes advantage of galaxy scalingrelations that map halo/stellar mass into star formationhistory (SFH), which when convolved with the DTD leadto a predicted BNS merger host galaxy mass function(this approach was previously proposed and used in thecontext 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 BNSmergers can be identified through the detection of elec-tromagnetic (EM) counterparts, but this is likely onlyachievable within a few hundred Mpc. We found thatfor a power law DTD characterized by index Γ and min-imum delay t min , O( ) host galaxies are required toreasonably constrain the DTD.The latter approach instead relies on a redshift map-ping of the BNS merger rate, which requires GW de-tections to z ∼ few , achievable with the next-generationEinstein Telescope (ET) and Cosmic Explorer (CE). In this approach, it is unlikely that EM counterparts can bedetected, but the individual redshift uncertainties fromthe GW data ( δ z / z ≈ . z ) can be overcome througha large number of anticipated detections, ∼ yr − .We found that with about a year of CE+ET data theDTD parameters, as well as the mass efficiency of BNSproduction, can be determined to about 10%.Here we continue our investigation of the DTD, withan alternative approach to the use of BNS merger hostgalaxies at z ≈ . Namely, unlike in Paper I, which usedscaling relations between mass and SFH, we explore theuse of detailed reconstructed SFHs for the individualhost galaxies. In § § § §
5. We adopt thePlanck 2015 cosmological parameters (Planck Collabo-ration et al. 2016): Ω M = . , Ω Λ = . , Ω b = . ,and H = . km s − Mpc − . GALAXY DATAWe use SFHs inferred from galaxy photometry in Lejaet al. (2019). The photometry is measured with theLAMBDAR code (Wright et al. 2016) from DR3 of the a r X i v : . [ a s t r o - ph . H E ] M a y -5 -4 -3 -2 -1 log(1 + z ) l og ( S F R / ( M fl / y e a r )) ContinuityPriorlogMPrior log( ˙n)[year − ] P D F Γ= − / ,t min =10MyrΓ= − / ,t min =100MyrΓ= − / ,t min =1GyrΓ= − ,t min =10MyrΓ= − ,t min =100Myr Γ= − ,t min =1GyrΓ= − / ,t min =10MyrΓ= − / ,t min =100MyrΓ= − / ,t min =1Gyr log( ˙n)[year − ] P D F Γ= − / ,t min =10MyrΓ= − / ,t min =100MyrΓ= − / ,t min =1GyrΓ= − ,t min =10MyrΓ= − ,t min =100Myr Γ= − ,t min =1GyrΓ= − / ,t min =10MyrΓ= − / ,t min =100MyrΓ= − / ,t min =1Gyr -5 -4 -3 -2 -1 log(1 + z ) l og ( S F R / ( M fl / y e a r )) ContinuityPriorlogMPrior log( ˙n)[year − ] P D F Γ= − / ,t min =10MyrΓ= − / ,t min =100MyrΓ= − / ,t min =1GyrΓ= − ,t min =10MyrΓ= − ,t min =100Myr Γ= − ,t min =1GyrΓ= − / ,t min =10MyrΓ= − / ,t min =100MyrΓ= − / ,t min =1Gyr log( ˙n)[year − ] P D F Γ= − / ,t min =10MyrΓ= − / ,t min =100MyrΓ= − / ,t min =1GyrΓ= − ,t min =10MyrΓ= − ,t min =100Myr Γ= − ,t min =1GyrΓ= − / ,t min =10MyrΓ= − / ,t min =100MyrΓ= − / ,t min =1Gyr Figure 1 . The BNS merger rate for two galaxies from the GAMA survey with distinct SFHs (
Left ) that peak at early (
Top )and later (
Bottom ) cosmic time (based on Equation 1). We show the SFH based on the two types of priors (yellow: logM; grey:continuity), with the solid line indicating the median SFH and the shaded regions marking the range of 16th to 84th percentile.
Middle and
Right are the merger rate PDFs from convolution of the 9 DTDs with the posterior distribution of SFH of the galaxymodeled based on continuity and LogM priors respectively. The vertical bars indicate the median values of the merger ratePDFs for each DTD, shown with the same line style and color.
Galaxy and Mass Assembly (GAMA) survey (Driveret al. 2011; Baldry et al. 2018), and includes 21 bandsranging from the far-UV (
GALEX ) to the far-IR (
Her-schel ). The galaxies have spectroscopic redshifts. Allgalaxies at . < z < . with stellar masses M ∗ > M (cid:12) 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- α modelwithin the Prospector inference machine (Leja et al.2017; Johnson & Leja 2017). This model includes aseven-parameter non-parametric star formation history,as well as flexible dust attenuation model, far-infrareddust re-emission via energy balance, and nebular emis-sion self-consistently powered by the stellar ionizing con-tinuum.A key source of uncertainty in SFH recovery is thechoice of prior (Carnall et al. 2019; Leja et al. 2019).This sensitivity occurs because the SEDs of stellar pop-ulations change slowly as a function of time; specifically,SEDs evolve roughly evenly in each logarithmic timestep (Ocvirk et al. 2006). As a result, the SFR inferredin adjacent time bins is typically highly degenerate. For-tunately, key outputs of galaxy SED-fitting such as the mass-to-light ratio and, to a lesser extent, the recentSFR, are calculated using moments of the SFH, whichare largely insensitive to this degeneracy (Bell et al.2003; Leja et al. 2019). However, the DTD does notinteract with these conserved quantities but instead cou-ples directly to the SFH, and is thus sensitive to this de-generacy. Accordingly, we perform the analysis here us-ing two different priors that assume opposite behaviorsin this degeneracy: one that favors bursty SFHs (here-after, logM) and one that favors smooth SFHs (here-after, continuity) (Leja et al. 2019). Two representativeexamples for galaxies with opposite SFHs are shown inFigure 1, for both the logM and continuity priors. Theseexamples highlight both the range of behavior and theassociated uncertainties (statistical and systematic) inthe reconstructed SFH.The two SFH priors adopted in this work encapsu-late the plausible range of choices (Leja et al. 2019).Star formation histories derived from high-resolutionoptical/near-IR spectroscopy can provide more preciseSFHs than photometry alone, although in most casesthis will be a modest improvement (Pacifici et al. 2012). DTD DETERMINATION METHODThe expected BNS merger rate at z = for galaxy i with star formation history ψ i ( z ) is given by: (cid:219) n i = ∫ z b = z b = λ dP m dt ( t − t b − t min ) ψ i ( z b ) dtdz ( z b ) dz b , (1)where dt / dz = −[( + z ) E ( z ) H ] − and E ( z ) = (cid:112) Ω m , ( + z ) + Ω k , ( + z ) + Ω Λ ( z ) ; t b is the cosmictime corresponding to z b ; λ is the BNS production massefficiency, assumed to be a fixed value of − M − (cid:12) , inde-pendent of redshift or environment; and dP m / dt is theDTD, likewise assumed to be independent of environ-ment. As in Papers I and II, we parameterize the DTDto follow a power law with index Γ , minimum delaytime t min , and a fixed maximum delay time t max = Gyr. Our results are not sensitive for a larger value of t max . This formulation of Equation 1 is identical to thatused in Paper I, except that there we took ψ to be adirect function of halo mass, M h .To generate the simulated data, we assume a fixednumber of N gal = galaxies that can serve as possi-ble hosts for BNS merger events. These are selected toserve as a representative subset of the GAMA samplethat preserves the mass distribution of the full sample.For a given DTD ( Γ and t min ) and SFH ( ψ i ( z ) ) we canthen estimate the mean merger rate, (cid:219) n i , for each galaxyusing Equation 1. This then gives each galaxy a differentprobability to host a BNS merger event for the differentDTDs (see Figure 1). As in Papers I and II, we use aset of 9 representative DTDs, with Γ = [− . , − , − . ] and t min = [ , , ] Myr.The number of BNS merger events, N i = { , , . . . } , observed for any given galaxy over a given period of time, ∆ t , follows a Poisson distribution based on the mergerrate: P ( N i | Γ , t min , ψ i ( z )) = Poisson ( (cid:219) n i ∆ t ) = ( (cid:219) n i ∆ t ) N i e − (cid:219) n i ∆ t N i ! (2)We then can simulate N i directly by drawing it from thePoisson distribution based on the associated rate: N i ∼ Poisson ( (cid:219) n i ∆ t ) (3)Once we have simulated a set of BNS merger events,we determine the constraining power they have on theunderlying DTD. Assuming the BNS merger events ineach galaxy are independent of each other and the SFHsare known, the corresponding likelihood is P ({ N i }| Γ , t min , { ψ i ( z )}) = N gal (cid:214) i = P ( N i | Γ , t min , ψ i ( z )) . (4) In this work (as well as in Papers I and II) we have focusedon a power law DTD, which is well motivated. However, it ispossible that the true DTD might deviate from this power lawform. This could potentially be investigated by comparing theresults of the analysis proposed here (using host galaxies at z ∼ )and the analysis in Paper II (using the redshift distribution ofBNS mergers from third-generation detectors). We further need to marginalize over the uncertainty onthe SFH of each galaxy: P ( N i | Γ , t min ) = ∫ P ( N i | Γ , t min , ψ i ( z )) P ( ψ i ( z )) d [ ψ i ( z )] . (5)We can approximate this integral by averaging over N samp of the samples from the SFH posteriors for eachgalaxy: P ( N i | Γ , t min ) ≈ N samp N samp (cid:213) j = ( (cid:219) n i , j ∆ t ) N i e − (cid:219) n i , j ∆ t N i ! (6)The resulting SFH-marginalized posterior of the DTD istherefore given by: P ( Γ , t min |{ N i }) ∝ P ({ N i }| Γ , t min |) P ( Γ , t min ) (7) ∝ N gal (cid:214) i = N samp (cid:213) j = ( (cid:219) n i , j ∆ t ) N i e − (cid:219) n i , j ∆ t N i ! (8)where we have assumed that the prior over Γ and t min isuniform such that P ( Γ , t min ) is a constant.We are interested in marginalizing over any particularset of events { N i } associated with the N gal possible hostgalaxies to forecast possible constraints on the DTD asa function of the total number of BNS merger events, N BNS . This gives: P ( Γ , t min | N BNS ) = ∫ P ( Γ , t min |{ N i } , N BNS ) P ({ N i }| N BNS ) d ({ N i }) (9)We can approximate this integral using N repeat realiza-tions of the observed BNS merger event counts, condi-tioned on the total number of events (cid:205) i N i being equalto N BNS . Combining this with our previous expressionthen gives: P ( Γ , t min | N BNS ) ≈ N repeat N repeat (cid:213) k = P ( Γ , t min |{ N i } k ) (10) ∝ N repeat (cid:213) k = N gal (cid:214) i = N samp (cid:213) j = ( (cid:219) n i , j ∆ t ) N i , k e − (cid:219) n i , j ∆ t N i , k ! (11)where again (cid:205) N gal i = N i , k = N BNS for each realization.Since the rate per galaxy is generally very small, (cid:219) n i (cid:28) , we expect the uncertainty to be dominatedby variation in the observed counts over the potential N gal host galaxies rather than the uncertainties in eachgalaxy’s SFH. As such, we opt to use the same SFHsused to generate the data to evaluate the posterior ratherthan trying to marginalize over them directly: P ( Γ , t min | N BNS ) ∼ N repeat (cid:213) j = N gal (cid:214) i = ( (cid:219) n i , j ∆ t ) N i , j e − (cid:219) n i , j ∆ t N i , j ! (12)It is critical to note that the posterior above is takenover the observed counts of all potential galaxy hosts,including those with N i = for which no BNS mergerevents have been detected. That is because these log( t min / Myr) Γ N BNS = 30 ContinuityLogM log( t min / Myr) Γ N BNS = 100 ContinuityLogM log( t min / Myr) Γ N BNS = 300 ContinuityLogM log( t min / Myr) Γ N BNS = 30 ContinuityLogM log( t min / Myr) Γ N BNS = 100 ContinuityLogM log( t min / Myr) Γ N BNS = 300 ContinuityLogM log( t min / Myr) Γ N BNS = 30ContinuityLogM log( t min / Myr) Γ N BNS = 100ContinuityLogM log( t min / Myr) Γ N BNS = 300ContinuityLogM
Figure 2 . The constraint achieved on the parameters of the DTD as a function of the number of host galaxies of BNS mergerevents (
Left : 30;
Middle : 100;
Right : 300) and for the two choices of SFH priors (red: logM; blue: continuity). In each row theinput model is marked with a yellow circle.
Top:
An input DTD with Γ = − / and t min = Myr.
Middle:
An input DTD with Γ = − and t min = Myr.
Bottom:
An input DTD with Γ = − / and t min = Myr. The contours show the 90% percentileconfidence. “non-detections” in aggregate contain a non-negligibleamount of information in a regime where (cid:219) n i (cid:28) andmerger events are rare. To illustrate this, we also com-pare the “complete” posterior distribution P ( Γ , t min | N BNS ) derived above with the biased posterior distribution, ˜ P ( Γ , t min | N BNS ) , ignoring the non-detections: ˜ P ( Γ , t min | N BNS ) ∼ N repeat (cid:213) j = N gal (cid:214) i = I( N i , k > ) ( (cid:219) n i , j ∆ t ) N i , j e − (cid:219) n i , j ∆ t N i , j ! , (13) where I( N i , j > ) is the indicator function that evaluatesto if the condition N i , j > is true and otherwise. Ingeneral, we expect that ignoring non-detections will biasthe inferred DTD, which we discuss in the next section.We compute the above posterior using N gal = galaxies and N repeat = realizations for the 9 differentassumed DTDs, and interpolate between their associ-ated Γ and t min parameters to obtain the probability fora different pair of Γ and t min values.To summarize, our inference procedure is as follows:1. We select a star formation history, ψ i ( z ) , for everygalaxy i from its SFH posterior distribution andcompute the corresponding BNS merger rate, (cid:219) n i ,using Equation 1.2. We then sample the number of BNS merger events, N i , from the corresponding Poisson distributionbased on Equation 2 for a given timescale ∆ t . Werepeat this process until the total number of eventsis N BNS = (cid:205) i N i .3. We repeat this procedure N repeat times to generatemany realizations of the BNS merger events, { N i } ,for varying SFHs.4. We then use the simulated BNS merger events andSFHs from these realizations to compute the DTDposteriors including and excluding galaxies thatdid not host detected BNS merger events usingEquations 12 and 13, respectively. RESULTSIn Figure 1 we show the BNS merger rate probabil-ity distribution function (PDF) of two galaxies from theGAMA survey, with SFHs chosen to peak at early andlate cosmic time (right panels), for the 9 different DTDmodels; the vertical lines in each panel indicate the valueof (cid:219) n for the mean SFH. We show the SFHs, and mergerrate PDFs for both sets of priors (middle panel: conti-nuity; right panel: logM). The SFHs from different as-sumed priors often do not overlap, illustrating the chal-lenge of inferring SFHs from photometry. In particular,the galaxy in the lower panel illustrates a specific degen-eracy where the photometry of star forming galaxies canoften be fit equally well with a continuous SFH or with alarge burst at (cid:46) Myr followed by a sharp drop in thestar formation rate (Leja et al. 2019). The figure alsoclearly indicates how the convolution of the DTDs withdifferent SFHs leads to distinct merger rate PDFs andhow the choice of SFH prior affects the resulting mergerrate PDFs.In Figure 2 we show the constraints achieved on theparameters of the DTD model as the number of observedhost galaxies increases from 30 to 100 to 300. We showthe results for three different injected DTD models andfor both the logM and continuity SFH priors. The con-tours mark the 90% confidence region. There are sev-eral key takeaway points from Figure 2. First, in thecase that the true DTD prefers short merger timescales,by having a steep power law slope and short t min (upperrow of Figure 2), the DTD parameters can be reason-ably constrained with fewer than O( ) host galaxies.Second, in other permutations of the DTD, O( ) hostgalaxies may be required to constrain the DTD parame-ters, but with a lingering degeneracy between Γ and t min .Third, the logM prior leads to tighter constraints onthe DTD parameters compared to the continuity priorbecause it allows for bursty SFHs that pick more well-defined timescales when convolved with the DTD; thecontinuity prior smooths the SFH and hence systemati-cally reduces its constraining power. log( t min / Myr) Γ N BNS = 300
DetectedHostGalaxiesOnly
ContinuityLogM
Figure 3 . Same as Figure 2, but for just a single injectedDTD (yellow circle) and a sample size of 300 host galaxies.Here we use the SFHs of only the host galaxies and ignorethe galaxies that did not host BNS mergers (Equation 13).We find that the resulting DTD parameters are biased withrespect to the input model, but not severely.