The Impact of metallicity evolution of the universe on the maximum mass of LIGO binary black holes
DDraft version September 5, 2019
Preprint typeset using L A TEX style AASTeX6 v. 1.0
THE IMPACT OF METALLICITY EVOLUTION OF THE UNIVERSE ON THE MAXIMUM MASS OF LIGOBINARY BLACK HOLES
Mohammadtaher Safarzadeh & Will M. Farr Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge MA 02138, USA [email protected] School of Earth and Space Exploration, Arizona State University Tempe AZ 85287, USA Department of Physics and Astronomy, Stony Brook University, Stony Brook, NY 11794, USA Center for Computational Astronomy, Flatiron Institute, New York, NY 10010, USA
ABSTRACTWe can be biased against observing massive black holes to merge in the local universe as the bounds onthe maximum black hole mass ( M max BH ) depends on the assumptions regarding the metallicity evolutionof the star forming gas across the cosmic time. We investigate the bounds on the metallicity evolution,mass distribution and delay times of the binary black hole sources based on the ten observed eventsby LIGO. We parametrize M max BH to be a function of metallicity which itself is modeled to evolve withredshift in either a modest or rapid fashion. Rapid metallicity evolution models predict a stringentbound of M max BH = + − M (cid:12) , while the bound on M max BH in the models with modest metallicity evolution is M max BH = + − M (cid:12) . Therefore, inferring M max BH from GW data depends on the assumed metal enrichmenthistory of the universe that is not severely constrained at the moment. INTRODUCTIONDetection of binary black holes (BBHs) byLIGO/Virgo has opened a new era in astronomy.Much effort has been focused on characterizing theformation scenario of these systems, whether they areborn in the field or assembled dynamically. There havebeen studies of the properties of the progenitors ofthese systems, largely based on the population synthesismodels which rely on uncertain physics in large parts.One of the key questions is whether there exists an up-per mass limit for black holes formed through stellar evo-lution. The theoretical models anticipate larger blackhole masses to be formed at lower metallicities sincethe line-driven winds would be quenched, and there-fore, a larger mass is available for collapse (Kudritzki& Puls 2000; Vink et al. 2001; Brott et al. 2011; Fryeret al. 2012). On the other hand, it is believed that pair-instability supernovae (PISN) creates a gap in BH massdistribution, with its location set by the pulsational pairinstability supernovae (PPISN) that consequently de-termine an upper limit on the most massive BHs thatcan potentially form at the lowest metallicities (Hegeret al. 2003; Belczynski et al. 2016; Yoshida et al. 2016;Woosley 2017; Marchant et al. 2018; Leung, Nomoto &Blinnikov 2019) due to the mass loss from pulsationspre-supernovae. This leads to the so-called second massgap between ≈
50 and 135 M (cid:12) for BHs formed from stel-lar core collapse. Given that the space-time volume thatLIGO is sensitive to probe scales with the primary massof the BBH as m / , if there is a cut off at around 50 M (cid:12) , the evidence for this should be there in the LIGOdata. There has been claims in the literature that the LIGOdata so far suggest the presence of a strong upper masscut for the black holes (Fishbach & Holz 2017; Talbot& Thrane 2018; Roulet & Zaldarriaga 2019). Fishbach& Holz (2017) conclude M max BH =
40 M (cid:12) , and a power lawindex of α < based on about 6 early BBH systemsdetected by LIGO. Roulet & Zaldarriaga (2019) arriveat M max BH = + − M (cid:12) and α ≈ by analyzing the tenobserved systems. LIGO collaboration analysis of theten events suggests that no more than 1% of black holesare more massive than 45 M (cid:12) (Abbott et al . et al. 2018).Moreover, they constrain the power law index of theprimary black hole to be α = . + . − . (90% credibility).One caveat that has been missing in the literaturewith regards to the M max BH is the influence of the metal-licity evolution of the universe. If black holes closeto the M max BH limit are born at the lowest metallicitiesof log ( Z / Z (cid:12) ) < − , then in order to detect the limit,we need the universe to have gone through such lowmetallicities for enough extended times to provide uswith observables. In other words, if PISN is active at log ( Z / Z (cid:12) ) < − , then if the universe lasted half of itsage at such low metallicities, then we would have ampleevidence for the presence of the upper mass limit. How-ever, if the universe spent only an insignificant lifetimeat such low metallicities, then there would have been notmuch star formation at such low metallicities, and there-fore, our power to detect the evidence for the presenceof such a mechanism would diminish.In this paper, we parametrize the distribution of theBBHs with six different parameters, and investigate theconstraining power in the ten observed events on them. a r X i v : . [ a s t r o - ph . H E ] S e p In our model, we tie the maximum black hole mass to themetallicity of the star forming gas, and we parametrizethe star forming gas metallicity to evolve either rapidlyor slowly with redshift. The M max BH is considered to be themaximum mass born at zero metallicity and thereforehow much time the universe is assumed to have spent atsuch low metallicities will determine the expected birthrate of such massive black holes.The structure of this paper is as follows: In § § § METHOD2.1. calculating the merger rate of the BBHs
The BBH formation rate as a function redshift percomoving volume per source frame time is defined as: dN form dm dm dt f dV c = λ BBH m − α m − β / C ( α, β ) ψ ( z ) , (1)where C ( α, β ) is the normalization constant given by: C ( α, β ) = ∫ m − α m − β dm dm (2) ψ ( z ) is the cosmic star formation rate density adoptedfrom Madau & Dickinson (2014): ψ ( z ) = . ( + z ) . + [( + z )/ . ] . M (cid:12) yr − Mpc − . (3)Here, λ BBH is the currently unknown BBH mass effi-ciency assumed not to evolve with redshift. The corre-sponding merger rate is given by: dN merge dm dm dt m dV c = ∫ ∞ t m ( z m ) P ( t m | t f ) dN form dm dm dt f dV c dt f (4) P ( t m | t f ) is the delay time distribution of the BBHs thatsets the probability of merging after t m of time is pastsince the formation of the binary. We set a minimumdelay time, t min = Myr and impose a maximum delaytime of 10 Gyr. P ( t m | t f ) = t − κ m / C ( κ ) (5)where C ( κ ) is the normalization constant given by: C ( κ ) = ∫ t max t min t − κ m dt . (6)Subsequently the merger rate in the detector frame is: dN merge dm dm dt d dz = dN merge dm dm dt m dV dV c dz + z (7)where the redshift derivative of the comoving volume is dV c / dz = ( π c / H )[ D L /( + z ) E ( z )] , where D L is the luminosity distance to the source, and H is the Hubbleconstant.In this framework, M min < m < m < m max1 , where M min = (cid:12) and M max BH is set by the metallicity as: m max1 = ( M max BH − c ) e − bZ ( z ,γ ) + c (8)with constants b = . , and c = . . This is shown inthe top panel of Figure 1. This parametrization matchesthe maximum mass of a blackhole as a function of metal-licity as derived in Belczynski et al. (2010) when we con-sider the maximum black hole mass to be 80 M (cid:12) . In laterseries of papers, Belczynski et al. (2016) have includedthe impact of pair-instability mass loss on black holebinaries which creates a second mass gap between 50-150 M (cid:12) for black holes. Our approach here is whetherthe presence of PISN could be inferred from the LIGOBBH systems. Z ( z , γ ) defines the metallicity evolution with redshiftthat we parametrize in two different ways: The firstmodel is a metallicity evolution in which the metallicitydrops exponentially with redshift, i.e., Z / Z (cid:12) = e − γ z . Inthe second model the metallicity is modeled as Z / Z (cid:12) = ( + z ) − γ where the metallicity evolution is much moremodest. While the impact of metallicity on the LIGOblack holes has been explored recently in other works(Kovetz et al. 2017; Neijssel et al. 2019), here we ex-plore its impact on the maximum BH mass that LIGOwould infer.The bottom panel of Figure 1 shows the different mod-els for the metallicity evolution of the universe that isadopted in this work. The metallicity in this work refersto the star formation rate weighted metallicity of the gasin the galaxies in which the BBHs are born.From observational perspective, metallicity studiesof damped Lyman alpha (DLA) systems at high red-shifts suggest a modest evolution at redshifts between . − (Pettini et al. 1997; Prochaska & Wolfe 2000;Cen et al. 2002; Kulkarni & Fall 2002; Prochaska et al.2003; Berg et al. 2016). If the star forming gas metallic-ity evolves in a similar manner, then low values of γ inour parametrization would be the closest model to theobserved metallicity evolution.The observed BBH merger rate is: dN obs dm dm dt d dz = dN merge dm dm dt d dz P det ( m , m , z ) , (9)where P det ( m , m , z ) is the detection probability of aBBH with masses of m , m , at redshift z . We note thatin this work we have assumed the mergers come fromthe same formation channel, and as such would followthe same λ BBH parameter.2.2. inference analysis
To perform our inference analysis, we proceed as fol-lows: Our model has 6 parameters that we fit for θ =( λ BBH , α, β, γ, κ, M max BH ). The posterior distribution ofthese parameters given the data is: P ( θ | data ) = P ( data | θ ) P ( θ ) (10) Z / Z fl M m a x B H [ M fl ] Belczynski + 10M maxBH = 80 M fl M maxBH = 40 M fl z -4 -3 -2 -1 Z / Z fl e − γz , γ = 1 . z ) − γ , γ = 1 . e − γz , γ = 2 . z ) − γ , γ = 2 . e − γz , γ = 3 . z ) − γ , γ = 3 . Figure 1 . Top panel: the parametrized maximum black holemass as a function of metallicity. The blue line shows theresults of Belczynski et al. (2010). The green and red linesshow the parametrized M max BH curves when the M max BH is setto 80, and 40 M (cid:12) respectively. Bottom panel: two differentmetallicity evolution models adopted in this work, one thatdrops exponentially with redshift (solid lines), and a power-law model (dashed lines).
The prior on our parameter is such that they arebound between < α, β, γ, κ < , and < M max BH / M (cid:12) < . We approximate P ( data | θ ) which we call for brevity P ( d | θ ) as follow: For each BBH event i , we havethe P ( m i , m i , z i | d i ) from the waveform analysis done byLIGO team. We have P ( m i , m i , z i | d i ) = P ( d i | m i , m i , z i ) P ( m i , m i , z i ) (11) and P ( d i | θ ) = P ( d i | m i , m i , z i ) P ( m i , m i , z i | θ ) (12)where by combining the last two equations we arrive at P ( d i | θ ) ∝ P ( m i , m i , z i | d i ) P ( m i , m i , z i | θ ) (13)Therefore, to compute P ( d i | θ ) , we draw N sample of ( m j , m j , z j ) i pairs from the posterior P ( m i , m i , z i | d i ) andcalculate P ( m j , m j , z j | θ ) . For each event i , we have P ( d i | θ ) = / N sample j = N sample (cid:213) j = dN merge dm dm dt d dz ( m j , m j , z j ) i | θ (14)The posterior distribution from N obs events is: P ( θ | d ) ∝ e − N eff | θ i = N obs (cid:214) i = P ( d i | θ ) , (15)where N eff | θ is the expected number of events given θ defined as : N eff | θ = ∫ M max BH ∫ m ∫ ∞ ∫ t obs dN obs dm dm dt d dz | θ dm dm dzdt (16)Where t obs is the total observing time by LIGO in O1and O2 runs. RESULTSFigure 2 shows the posterior distribution for the sixparameters of our model when the metallicity evolutionis modeled as Z / Z (cid:12) = e − γ z . The median BBH efficiencyis predicted to be ≈ × − / M (cid:12) . This birth rate isrobust and is not affected by our metallicity evolutionparametrization. However, we note that we have as-sumed λ BBH to be constant in this work, but BBH for-mation is intrinsically tied to this parameter through thewind mass loss. So the reader should note that the infer-ences on λ BBH in this work is with the prior assumptionthat λ BBH is non-evolving with redshift itself.While the simulation can not put stringent constraintson α , β , and γ , one can say large values for β , andsmall values of γ are disfavored. The posterior on κ issuggestive of a shallow slope and therefore a preferencefor long delay times for the BBHs. The anti-correlationbetween the birth efficiency λ BBH , and κ is due to thefact that if a model with long delay times is chosen,then it should be balanced out with lower birth rateefficiency since long delays increase the number densityof the BBH mergers in the local universe. Of all theparameters in our model, it is the M max BH that is very wellconstrained to be M max BH = + − in this model.Figure 3 shows the same results but for the modelwith metallicity evolution modeled as Z / Z (cid:12) = ( + z ) γ .It appears that all the parameters expect M max BH have thesame posterior distribution. The bounds on the M max BH that is less constrained and is M max BH = + − . Not onlythe median value is larger, but the upper bound extendsto a much larger value.The impact of the assumptions about the metallicityevolution on the M max BH should be understood as follows:In our model, the maximum black hole mass enters ourcalculation in a non-trivial manner. M max BH sets the maxi-mum mass that a black hole can have at zero metallicity.If in one model, the metallicity evolution is modest, andbarely touches very low metallicities, then to explain theLIGO black holes one needs to push the M max BH to largevalues to open the room for the model to fit the mas-sive LIGO systems such as GW170729, and GW170823.This is the case when the metallicity evolution follows ( + z ) − γ . However, if the universe spends much of itscosmic time at very low metallicities, then one can easilyexplain GW170729 and GW170823 by the star forma-tion at high redshifts. The idea can be best seen in theanti-correlation between γ and M max BH in Figure 3: Largevalues of γ which translate into a faster drop in metal-licity, leads to lower values of M max BH and vice versa.A different perspective on our results is provided inFigure 4. In the left panel the thin black lines show pos-terior draws from the M max BH and γ from the model withexponential metallicity evolution with redshift. Thethick red line shows the median predicted evolution. Foreach of the ten observed BBH systems, we show thebounds on the mass and redshift of the primary (moremassive) black hole. We are not fitting these data points,we are showing them to be compared to the maximumpossible black holes that could be formed above a certainredshift range, such that after a delay time they mergerin the local universe. The right panel shows the samebut for the model with power-law metallicity evolutionwith redshift.The left panel of Figure 5 shows the bounds on themetallicity evolution itself in the two models. The solidlines and the shaded region of the same color show themedian and the 16th-84th percentile range for the eachof the metallicity models. The evolution shown withblue is more consistent with the observations of the DLAsystems at high redshifts which suggest a modest evo-lution of their metallicity with redshift (Pettini et al.1997; Prochaska & Wolfe 2000; Cen et al. 2002; Kulka-rni & Fall 2002; Prochaska et al. 2003; Berg et al. 2016).Right panel of Figure 5 shows the posterior BBHmerger rate as a function of redshift for the modelwith Z / Z (cid:12) ∝ e − γ z (red shaded region showing the 16th-84th percentile range. The blue line and shaded regionshow the same for the model with metallicity evolutionparametrized as ∝ ( + z ) − γ . The dashed black line is the λ BBH ψ ( z ) , which shows what the merger rate wouldbe if there is no delay time for the BBHs. The differentmetallicity evolution models did not have a discernibleimpact on the merger rate of the binaries, and thereforeon the maximum mass would be the best probe of themetallicity evolution in this picture. SUMMARY & DISCUSSIONOur results can be summarized as follows: If maxi-mum black hole mass is set at close to zero metallicity,then in order to infer it from data, it is crucial to have alarge part of the cosmic time to have a metallicity closeto zero to generate BBH systems that can probe themass limit. In other words, if for example, we lived ina universe in which the metallicity never dropped below0.1 Z (cid:12) , then there would have been little hope to con-strain a parameter that requires probing metallicitiesclose to − Z (cid:12) . In our two models, one prescriptionof the metallicity evolves rapidly with redshift and theother evolves rather smoothly. The bounds on the M max BH are much more stringent in the model with a rapid dropof metallicity with redshift (i.e., Z ∝ e γ z ), compared tothe model in which metallicity is modeled as Z ∝ ( + z ) − γ . Similarly, if we lived in a universe in which the veryheavy black holes tend to be born in close binaries,and therefore merger rapidly, then we would be bi-ased against finding them in the local universe. Sucha parametrization is not considered in this work, butit would have resulted in the same conclusions that wehave reached so far.Therefore, any claim as to the presence of an upperlimit on the M max BH should be taken with the caveat thatwe can be easily biased against them, and the bound onthe M max BH depends on our assumptions with regard to (i)how these systems are born (metallicity range) and howdoes the universe on average evolve in metallicity, and(ii) whether the more massive systems tend to clusterin a parameter space in delay times that we would bebiased against them.This work was supported by the National ScienceFoundation under grant AST14-07835 and by NASA un-der theory grant NNX15AK82G as well as a JTF grant.MTS is grateful to the Center for Computational Astro-physics for hospitality during the course of this work.REFERENCES Abbott et al ., B. P., Collaboration, t. V., Abbott, B. P., et al.2018, 1811.12940Belczynski, K., Bulik, T., Fryer, C. L., et al. 2010, TheAstrophysical Journal, 714, 1217Belczynski, K., Heger, A., Gladysz, W., et al. 2016, Astronomy& Astrophysics, A97Berg, T. A. M., Ellison, S. L., S´anchez-Ram´ırez, R., et al. 2016,Monthly Notices of the Royal Astronomical Society, 3021 Brott, I., de Mink, S. E., Cantiello, M., et al. 2011, Astronomy &Astrophysics, A115Cen, R., Ostriker, J. P., Prochaska, J. X., & Wolfe, A. M. 2002,The Astrophysical Journal, 598, 741Fishbach, M., & Holz, D. E. 2017, The Astrophysical JournalLetters, 851, L25Fryer, C. L., Belczynski, K., Wiktorowicz, G., et al. 2012, TheAstrophysical Journal, 749, 91 log( λ ) = − . +0 . − . . . . . α α = . +0 . − . . . . . β β = . +0 . − . . . . . γ γ = . +1 . − . . . . . = . +0 . − . . . . . log( λ ) M m a x B H . . . . α . . . . β . . . . γ . . . . M maxBH M maxBH = . +9 . − . Figure 2 . Results of MCMC simulation on six parameters in our model by fitting 10 LIGO events in O1 and O2 observing runs.In this model the metallicity evolution is modeled as Z / Z (cid:12) = e γ z . Heger, A., Fryer, C. L., Woosley, S. E., Langer, N., & Hartmann,D. H. 2003, The Astrophysical Journal, 591, 288Kovetz E. D., Cholis I., Breysse P. C., Kamionkowski M., 2017,PhRvD, 95, 103010Kudritzki, R.-P., & Puls, J. 2000, Annual Review of Astronomyand Astrophysics, 38, 613Kulkarni, V. P., & Fall, S. M. 2002, The Astrophysical Journal,580, 732Leung S.-C., Nomoto K., Blinnikov S., 2019, arXiv,arXiv:1901.11136Madau, P., & Dickinson, M. 2014, Annual Review of Astronomyand Astrophysics, 52, 415 Marchant P., Renzo M., Farmer R., Pappas K. M. W., TaamR. E., de Mink S., Kalogera V., 2018, arXiv, arXiv:1810.13412Neijssel C. J., et al., 2019, arXiv, arXiv:1906.08136Pettini, M., Smith, L. J., King, D. L., & Hunstead, R. W. 1997,The Astrophysical Journal, 486, 665Prochaska, J. X., Gawiser, E., Wolfe, A. M., Castro, S., &Djorgovski, S. G. 2003, The Astrophysical Journal, 595, L9Prochaska, J. X., & Wolfe, A. M. 2000, The AstrophysicalJournal, 533, L5Roulet, J., & Zaldarriaga, M. 2019, Monthly Notices of the RoyalAstronomical Society, 484, 4216 log( λ ) = − . +0 . − . . . . . α α = . +0 . − . . . . . β β = . +0 . − . . . . . γ γ = . +0 . − . . . . . = . +0 . − . . . . . log( λ ) M m a x B H . . . . α . . . . β . . . . γ . . . . M maxBH M maxBH = . +16 . − . Figure 3 . Same as in Figure 2 but the metallicity evolution is modeled as Z / Z (cid:12) = ( + z ) − γ . Since in this parametrization themetallicity evolves more gradually with redshift, the impact is evident in the detail of the M max BH posterior as larger black holemasses would be allowed in this model. Talbot, C., & Thrane, E. 2018, The Astrophysical Journal, 856,173Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001,Astronomy & Astrophysics, 369, 574 Woosley, S. E. 2017, The Astrophysical Journal, 836, 244Yoshida T., Umeda H., Maeda K., Ishii T., 2016, MNRAS, 457,351 z M m a x B H GW151012GW151226GW170104GW170608GW170729 GW170809GW170814GW170818GW170823GW150914 z M m a x B H GW151012GW151226GW170104GW170608GW170729 GW170809GW170814GW170818GW170823GW150914
Figure 4 . Left Panel: thin black lines show posterior draws from the M max BH and γ from the model with exponential metallicityevolution with redshift. The thick red line shows the median predicted evolution. The heavier black hole mass and redshift ofthe ten observed BBH systems in O1 and O2 observing runs are plotted. The M max BH = + − in this model. Right panel:
Thesame as left panel, but for a metallicity evolution parametrized as ∝ ( + z ) γ . The bounds on the maximum mass in this modelis less constrained and is predicted to be M max BH = + − . z -4 -3 -2 -1 Z / Z fl Z/Z fl ∝ e − γz Z/Z fl ∝ (1 + z ) γ z -1 BB H M e r g e r R a t e ( G p c − y e a r − ) Z/Z fl ∝ e − γz Z/Z fl ∝ (1 + z ) γ Figure 5 . Left Panel: thin black lines show posterior draws of the metallicity evolution in the model with exponential metallicityevolution with redshift. The red line and the shaded region show the median and the 16th-84th percentile range.
Right panel: posterior BBH merger rate as a function of redshift for the model with Z / Z (cid:12) ∝ e − γ z (red shaded region showing the 16th-84thpercentile range. The blue line and shaded region show the same for the model with metallicity evolution parametrized as ∝ ( + z ) − γ . The dashed black line is the λ BBH ψ ( z ))