The Population III origin of GW190521

For simplicity, we here estimate the efficiency of producing BHB mergers like GW190521 from Pop III stars via classical binary stellar evolution without detailed modelling of the mass transfer and CE evolution, as well as SN kicks. Instead, we identify all binaries with Roche lobe overflow of either one of the two stars that is massive enough ($M_{\rm ZAMS}>50~{}\rm M_{\odot}$, Kinugawa et al. 2020a) to have a convective envelope. We assume that these binaries will undergo CE evolution due to unstable mass transfer, and merge within a Hubble time, $t_{\rm H}\simeq 13.7$ Gyr, following a power-law delay time distribution $\mathcal{P}(t_{\rm GW})$ with a slope of -1 in the range of 3 to $10^{4}$ Myr, based on previous studies (e.g., Belczynski et al. 2017)¹¹1The results are not sensitive to the detailed shape of $\mathcal{P}(t_{\rm GW})$ as long as it is dominated by mergers with short delay times (a few Myr), which is a common outcome in classical binary stellar evolution models (e.g., Belczynski et al. 2017; Kinugawa et al. 2020a; Tanikawa et al. 2020b)..

Because close binaries with $a\lesssim 10$ AU in the FD model are rare, the efficiency of Pop III BHB mergers²²2The efficiency of BHB mergers for a stellar population is defined as the average number of BHBs that can merge within a Hubble time per unit stellar mass. similar to GW190521 via CE evolution with PPI is low, $f_{\rm BHB}=3\times 10^{-6}\ \rm M_{\odot}^{-1}$. Considering other models in Liu et al. (2020) with much (a factor of $\gtrsim 15$) smaller cluster sizes and more close binaries, we obtain an extreme upper limit of $f_{\rm BHB}\sim 10^{-4}\ \rm M_{\odot}^{-1}$ in the optimistic case. Note that we do not include SN kicks for the primary, which can enhance the chance of CE evolution by shrinking the binary orbit, especially important for BHBs in the PPISN mass gap (see fig. 33 of Tanikawa et al. 2020b). Therefore, we may have underestimated the efficiency of BHB mergers driven by CE evolution.

We next consider how BHB evolution depends on environment, which can drive initially wide binaries to merge without a CE phase that is only available for initially close binaries. Actually, the dynamics of BHB coalescence via environmental effects is very complex, involving various astrophysical aspects such as a clumpy interstellar medium, dark matter distribution, nuclear star clusters (NSCs), supermassive BHs and AGN discs (e.g., Roškar et al. 2015; Antonini & Rasio 2016; Tamfal et al. 2018; Leigh et al. 2018; Choksi et al. 2019; Ogiya et al. 2019; Zhang et al. 2019; Secunda et al. 2019). For simplicity, we focus on one mechanism particularly important for massive Pop III BHBs: the dynamical hardening (DH) by 3-body interactions with surrounding (low-mass) stars in dense star clusters (Antonini & Rasio, 2016; Leigh et al., 2018). Below we estimate the delay time for a BHB merger with DH.

Following Liu & Bromm (2020b), the evolution of the semi-major axis $a$ for a BHB driven by 3-body interactions with surrounding stars and gravitational wave (GW) emission, can be written as (Sesana & Khan, 2015)

The first term represents 3-body interactions, the second the energy and angular momentum loss via GWs, and

Here, $m_{1}$ and $m_{2}$ are the masses of the primary and secondary, $m\equiv m_{1}+m_{2}$, $\sigma_{\star}$ and $\rho_{\star}$ are the velocity dispersion and stellar density around the BHB, $F(e)=(1-e^{2})^{-7/2}[1+(73/24)e^{2}+(37/96)e^{4}]$ given the eccentricity $e$, and $H\sim 15-20$ is a dimensionless parameter, which we set to $17.5$ for simplicity.

The binary system first undergoes a phase dominated by DH with a characteristic semimajor axis, estimated by imposing $(da/dt)|_{\mathrm{3B}}=(da/dt)|_{\mathrm{GW}}$: $a_{\star/\mathrm{GW}}=\left(B/A\right)^{1/5}$. We assume a constant eccentricity in the DH-dominated phase for simplicity. The duration of this stage for $a$ to reach $a_{\star/\mathrm{GW}}$, the DH timescale, can be estimated with

where $a_{0}$ is the initial semi-major axis. We set $t_{\rm DH}=0$, $a_{\star/\rm GW}=a_{0}$ for $a_{0}<(B/A)^{1/5}$.

After the DH-dominated phase ($a\lesssim a_{\star/\rm GW}$), the evolution is dominated by GW emission. The time spent in GW-driven inspiral before the final coalescence can be estimated with

The total time of inspiral in the dense star cluster is then given by $t_{\rm SC}=t_{\rm DH}+t_{\rm col}$. Other than the binary parameters ($m_{1}$, $m_{2}$, $e$ and $a_{0}$), the key parameters that determine $t_{\rm SC}$ are $\rho_{\star}$ and $\sigma_{\star}$. Besides, only hard binaries will be further hardened (instead of destroyed) by 3-body interactions, which must satisfy $Gm_{1}m_{2}(1-e)/[a(1+e)]>m_{\star}\sigma_{\star}^{2}$ in order to survive disruptions close to the apocenter, where $m_{\star}=1\ \rm M_{\odot}$ is the typical mass of surrounding stars.

To apply the DH machinery to Pop III BHBs, we need to know the properties of host clusters ($\rho_{\star}$ and $\sigma_{\star}$), initial binary parameters and dynamics of Pop III BHs/BHBs falling into dense star clusters. Here we consider Pop III BHBs in nuclear star clusters (NSCs), which occupy the innermost regions of most galaxies (Neumayer et al. 2020), and are also common in high-$z$ atomic-cooling halos with $M_{\rm h}\gtrsim 10^{8}\ \rm M_{\odot}$ (e.g., Devecchi & Volonteri 2009; Devecchi et al. 2010, 2012).

For NSC properties, we explore three cases with $\sigma_{\star}\simeq 10,\ 33\text{ and }100\ \rm km\ s^{-1}$, $\rho_{\star}=10^{4},\ 7\times 10^{5}\text{ and }5\times 10^{7}\ \rm M_{\odot}\ pc^{-3}$, corresponding to NSCs in halos with masses $M_{\rm h}\sim 10^{8},\ 10^{10}\text{ and }10^{12}\ \rm M_{\odot}$, respectively, which dominate the halo populations at $z_{\rm crit}\sim 9,\ 5.5$ and 4 for $2\sigma$ peaks. These three models are meant to cover most ($\gtrsim 90\%$) of the region occupied by observed NSCs with masses $M_{\rm NSC}\lesssim 10^{7}\ \rm M_{\odot}$ for $M_{\rm h}\lesssim 10^{12}\ \rm M_{\odot}$ (see fig. 7 of Neumayer et al. 2020) in the parameter space. In reality, the host NSC properties exhibit scatters and also evolve with redshift. The realistic situation will be a mixture of our models. To populate the $\rho_{\star}$-$\sigma_{\star}$ space with observed NSCs (see Fig. 4 in Appendix A), we estimate $\sigma_{\star}$ with $\sqrt{(1/2)GM_{\rm NSC}/r_{\rm eff}}$, given the half-mass(light) radius $r_{\rm eff}$. We set $\rho_{\star}$ to the stellar density within a characteristic radius $r_{\rm BH}$ where the enclosed stellar mass is twice the mass of the infalling BHB ($m\sim 140\ \rm M_{\odot}$), assuming that the stars follow the Dehnen profile (Dehnen, 1993) with an inner slope of $\gamma=1$ (see Sec. 4.1.2 in Liu & Bromm 2020b)³³3$\gamma=1$ is a conservative choice, which is also consistent with the trend $r_{\rm eff}\propto M_{\rm NSC}^{1/2}$ (Neumayer et al., 2020).. Here $r_{\rm BH}$ describes how deep the BHB can sink into the cluster by mass segregation.

We substitute the candidate BHB progenitors of GW190521 from the FD model into the DH scheme (Equ. 1-4), in which only hard binaries without CE evolution are considered. In this way, we derive the NSC inspiral time distributions $p(t_{\rm SC})$ for the three NSC models. The results with PPI can be well fitted by $p(t_{\rm SC})\propto t_{\rm SC}^{\alpha}$ with $\alpha=2$, 4 and 8, in the ranges of $2.5-13$ Gyr, $0.44-1.1$ Gyr and $69-90$ Myr, for $M_{\rm h}\sim 10^{8},\ 10^{10}\text{ and }10^{12}\ \rm M_{\odot}$, respectively (see Fig. 5 in Appendix A). The corresponding efficiencies of mergers are $f_{\rm BHB,0}=1.3\times 10^{-4}$, $7.6\times 10^{-5}$ and $7.8\times 10^{-6}\ \rm M_{\odot}$. Turning off PPI will enhance $f_{\rm BHB,0}$ by a factor of $\sim 5$, with $p(t_{\rm SC})$ almost unchanged. The final efficiency is given by $f_{\rm BHB}=f_{\rm NSC}f_{\rm BHB,0}$, where $f_{\rm NSC}\sim 0.2-0.8$ (Neumayer et al., 2020) is the fraction of galaxies that host NSCs.

Now, taking into account the BH/BHB dynamics in host halos during cosmic structure formation, the final delay time $t_{\rm GW}$ can be written as $t_{\rm GW}=t_{\rm SC}+t_{\rm IF}+t_{\rm DF}$, where $t_{\rm IF}$ and $t_{\rm DF}$ are the timescales for the BHB to fall into the final host halo and sink into the NSC via dynamical friction of gas. In general, Pop III star formation peaks at $z\sim 10$ in low-mass halos ($M_{\rm h}\sim 10^{7}-10^{8}\ \rm M_{\odot}$) at $2\sigma$ peaks, and the host halos of Pop III BHs will then grow or merge into more massive halos, in which the Pop III BHs can further sink into NSCs. The halo growth/merger process is captured by the timescale $t_{\rm IF}$. Given the mass $M_{\rm h}$ of the final host halo, $t_{\rm IF}$ can be estimated with $t(z=z_{\rm crit})-t(z=10)$, where $t$ is the age of the universe as a function of redshift, and $z_{\rm crit}$ is the redshift at which $M_{\rm h}$ matches the critical mass at $2\sigma$ peaks. We have $t_{\rm IF}\simeq 0.07$, 0.56 and 1.06 Gyr, for $M_{\rm h}\sim 10^{8},\ 10^{10}\text{ and }10^{12}\ \rm M_{\odot}$, respectively.

The (gas) dynamical friction timescale $t_{\rm DF}$ depends on the (initial) distance $r$ of the BHB to the galaxy center and the gas density profile. We adopt the $t_{\rm DF}$ formula in Boco et al. (2020, see their equ. 18 and table 1) with fixed Coulomb logarithm and assuming circular orbits for the BHBs. This formula can be applied to different density profiles, parameterized by the gas mass $M_{\rm gas}$ and half-mass radius $R_{e}$. For the $M_{\rm h}=10^{8}\ \rm M_{\odot}$ model at $z=9$, which is a typical high-$z$ atomic cooling halo, the singular isothermal sphere (SIS) is a good approximation to the density profiles found in simulations at $r\gtrsim 1$ pc (see equ. 2 in Safarzadeh & Haiman 2020). We consider the entire halo⁴⁴4The $t_{\rm DF}$ formula in Boco et al. 2020 does not include the effect of dark matter, which is valid for the inner regions ($\lesssim 1$ kpc) of massive halos ($M_{\rm h}\gtrsim 10^{10}\ \rm M_{\odot}$). However, in our case with $M_{\rm h}=10^{8}\ \rm M_{\odot}$ at $z=9$, one has to consider dark matter even for $r\lesssim 1$ kpc, as the virial radius is $R_{\rm vir}\simeq 1.44$ kpc. Therefore, we boost the initial specific angular momentum of the BHB by a factor of $(\Omega_{\rm m}/\Omega_{\rm b})^{1/2}$ to include the gravity of dark matter. The dynamical friction timescale obtained in this way is consistent with that in Safarzadeh & Haiman (2020) within a factor of 2. with $M_{\rm gas}\sim M_{\rm h}\Omega_{\rm b}/\Omega_{\rm m}=1.55\times 10^{7}\ \rm M_{\odot}$ and $R_{e}=0.5R_{\rm vir}=0.72$ kpc, which gives $t_{\rm DF}\sim 140\ \mathrm{Myr}\times(r/10\ \mathrm{pc})^{2}$. For the $M_{\rm h}=10^{12}\ \rm M_{\odot}$ model at $z=4$, as a typical progenitor for local early-type galaxies, we follow Boco et al. (2020) to focus on the gas-dominated inner region with $R_{e}=1$ kpc and $M_{\rm gas}\sim 2(R_{e}/R_{\rm vir})^{0.3}M_{\rm h}\Omega_{\rm b}/\Omega_{\rm m}=9\times 10^{10}\ \rm M_{\odot}$, under a Sersic profile with index $n=1.5$, such that $t_{\rm DF}\sim 150\ \mathrm{Myr}\times(r/10\ \mathrm{pc})^{2.7}$. Finally, for the intermediate case with $M_{\rm h}=10^{10}\ \rm M_{\odot}$ at $z=5.5$, we adopt the Hernquist profile, which is between the aforementioned SIS and Sersic models in terms of compactness. Again, we consider the gas-dominated inner region with $R_{e}=1$ kpc and $M_{\rm gas}\sim 2(R_{e}/R_{\rm vir})^{0.5}M_{\rm h}\Omega_{\rm b}/\Omega_{\rm m}=9.7\times 10^{8}\ \rm M_{\odot}$, leading to $t_{\rm DF}\sim 40\ \mathrm{Myr}\times(r/10\ \mathrm{pc})^{2.5}$. Although our modelling of gas distributions is highly idealized, it captures the trend that gas is more concentrated due to stronger cooling (with higher metallicities) in more massive halos at lower redshifts.

Finally, we consider the distribution of Pop III BHs/BHBs in high-$z$ galaxies. We derive the fraction of Pop III BHs $f_{\rm BH}(<r)$ enclosed within $r$ around galaxy centers in atomic-cooling halos ($M_{\rm h}\gtrsim 10^{8}\ \rm M_{\odot}$), from the cosmological simulation in Liu & Bromm (2020b, see Fig. 6 in Appendix A), where the dynamical friction of gas is unresolved. Then the (normalized) distribution of $t_{\rm GW}$ takes the form $\mathcal{P}(t_{\rm GW}|M_{\rm h})=\int_{r_{\min}}^{r_{\max}}dr\left[p(t_{\rm GW}-t_{\rm IF}(M_{\rm h})-t_{\rm DF}(r|M_{\rm h}))df_{\rm BH}(<r)/dr\right]$, where we adopt $r_{\min}=10$ pc and $r_{\max}=10^{3}$ pc, since $f_{\rm BH}(<r)$ at $r<10$ pc is negligible ($\lesssim 10^{-3}$), and $t_{\rm DF}(r|M_{\rm h})>t_{\rm H}$ for $r>10^{3}$ pc. The resulting distributions are shown in Fig. 2.