On the formation of GW190521-like binary black hole merger systems Free

ABSTRACT

GW190521 is the most massive merging binary black hole (BBH) system detected so far. At least one of the component BHs was measured to lie within the pair-instability supernova (PISN) mass gap (∼50–135 M_⊙), making its formation a mystery. However, the transient observed signal allows alternative posterior distributions. There was suggestion that GW190521 could be an intermediate-mass ratio inspiral, with the component masses m₁ ∼ 170 M_⊙ and m₂ ∼ 16 M_⊙, happening to straddle the PISN mass gap. Under this framework, we perform binary population synthesis to explore the formation of GW190521-like systems via isolated binary evolution. We numerically calculate the binding energy parameter for massive stars at different metallicities, and employ them in our calculation for common envelope evolution. Our results prefer that the progenitor binaries formed in metal-poor environment with |$\rm Z\le 0.0016$|⁠. The predicted merger rate density within redshift z = 1.1 is |${\sim} 4\times 10^{-5}{\text {--}}5\times 10^{-2} \, \rm Gpc^{-3}yr^{-1}$|⁠. We expect that such events are potentially observable by upcoming both space and ground-based gravitational wave detectors.

1 INTRODUCTION

Detection of gravitational waves (GWs) serves us an alternative way to observe the Universe. Since the first GW event GW150914 was discovered by the ground-based detectors advanced LIGO (aLIGO) and later joined advanced Virgo, the number of binary black hole (BBH) merger events has increased to ∼100 (Abbott et al. 2016a, b, 2019, 2021; The LIGO Scientific Collaboration 2021). The observed GW signals are classified as coalescing BBH, binary neutron star, and neutron star-BH systems. GW170817 is the only GW source with electromagnetic counterpart definitely observed (Abbott et al. 2017a).

GW190521, observed on 2019 May 21 at 03:02:29 utc, is the most massive merging BBH system detected so far (Abbott et al. 2020a, b). The association of GW190521 with the candidate counterpart ZTF19abanrhr reported by Zwicky transient facility (ZTF; Graham et al. 2020) is still inconclusive (Ashton et al. 2021; Nitz & Capano 2021; Palmese et al. 2021). Under the assumption that GW190521 is a quasi-circular BBH coalescence, the estimated individual component masses are |$m_1=85_{{-14}}^{+21}{\rm M}_{\odot }$|⁠, |$m_2=66_{{-18}}^{+17}{\rm M}_{\odot }$|⁠, and the total mass |$150_{{-17}}^{+29}$| M_⊙ within 90 per cent credible region, providing direct evidence of intermediate mass BHs (IMBHs; Abbott et al. 2020a, b). Gamba et al. (2023) drew similar results but under hyperbolic orbit hypothesis. Romero-Shaw et al. (2020) claimed that GW190521 may be an eccentric binary merger with aligned spins. Gayathri et al. (2022) interpreted this signal under the combination of both eccentricity and spin precession configuration. Barrera & Bartos (2022) estimated the ancestral mass of GW190521 and also favoured the heaviest parental BH mass in the pair-instability supernova (PISN) mass gap (between ∼50 and 135 M_⊙; Yusof et al. 2013; Belczynski et al. 2016a).

Since the BH masses in GW190521-like events challenge the standard stellar evolutionary theory, there have been various models proposed to interpret their formation, including dynamical binary formation in dense stellar clusters (Rodriguez et al. 2019; Fragione, Loeb & Rasio 2020; Romero-Shaw et al. 2020; Arca-Sedda et al. 2021; Anagnostou, Trenti & Melatos 2022; Rizzuto et al. 2022; Gamba et al. 2023), additional gas accretion, and hierarchical mergers in active galactic nuclei (Tagawa, Haiman & Kocsis 2020; Tagawa et al. 2021; and references therein), and the primordial BH scenarios (De Luca et al. 2021). Alternatively, Palmese & Conselice (2021) suggested that GW190521 may be the merger of central BHs from two ultradwarf galaxies. However, the origin of this event as an isolated binary still cannot be excluded (Belczynski et al. 2020a; Farrell et al. 2021; Kinugawa, Nakamura & Nakano 2021; Tanikawa et al. 2021a). In addition, the exact boundaries of the PISN mass gap are in dispute, due to the uncertainties in stellar evolution and SN simulation, which may entail a reassessment (Woosley 2017; Farmer et al. 2019; Marchant et al. 2019; Mapelli et al. 2020; Vink et al. 2021)

Nevertheless, GW190521 is qualitatively different from previous GW sources, not only because it was the most massive GW source observed to date, but also this transient signal was found with only a short duration of approximately 0.1 s, and only around four cycles in the frequency band 30–80 Hz, so multimodal posterior distributions would be consequently ineluctable (Fishbach & Holz 2020; Bustillo et al. 2021; Nitz & Capano 2021; Estellés et al. 2022). Among them, Nitz & Capano (2021) suggested that GW190521 may be an intermediate-mass-ratio inspiral (IMRI), with the component masses of |$m_1\sim 170\, {\rm M}_{\odot }$| and |$m_2\sim 16\, {\rm M}_{\odot }$|⁠, straddling the PISN mass gap. Comparison of the parameters derived by Abbott et al. (2020a) and Nitz & Capano (2021) are shown in Table 1.

Inspired by the results of Nitz & Capano (2021), here we attempt to interpret the formation of GW190521 assuming that it was an IMRI through isolated binary evolution, and investigate the properties of their progenitor binaries as well as the possible distributions of natal kicks on the two component BHs, which had promoted their coalescence within the Hubble time τ_H. The information on the BH kicks is crucial in understanding the formation of massive BHs.

The paper is structured as follows. In Section 2, we describe the main features of our binary population synthesis (BPS) models. The calculated results of BPS are presented in Section 3. We then discuss our results in Section 4, and summarize our main conclusions in Section 5.

2 MODEL

All calculations are carried out by using the BPS code bse, originally developed by Hurley, Pols & Tout (2000, 2002) and its updated version bseemp¹ (Tanikawa et al. 2020, 2021b), with the extension to very massive (up to |$M\sim 1300\, {\rm M}_{\odot }$|⁠) and extremely metal-poor stars (down to Z = 10⁻⁸Z_⊙), based on the stellar models computed by the hoshi code. The bseemp code also takes advantage of new stellar-wind and remnant-formation prescriptions, as well as the implementation of pair-instability and pulsational-pair-instability supernova (PISN/PPISN), which all play a role in stellar/binary evolution.

2.1 Stellar wind mass-loss

The masses of BHs are predominantly set by their pre-SN masses, which are mainly affected by stellar wind mass-loss and binary interactions. Mass-loss via stellar winds can significantly influence the fate of massive stars (Fryer et al. 2002). Here, we use the semi-empirical stellar wind prescription (referred as Vink et al. winds) in Belczynski et al. (2010), and consider metallicity dependence for luminous blue variables (LBVs; Tanikawa et al. 2021b). It has been demonstrated that this wind prescription results in more massive pre-SN objects and heavier BHs compared with the traditional one (Belczynski et al. 2010). We ignore the influence of stellar rotation on wind loss.

For massive O and B stars, the wind mass-loss rate |$\dot{M}_{\rm W}$| in units of |${\rm M}_{\odot }\, {\rm yr}^{-1}$| is (Vink, de Koter & Lamers 2001)

with |$12\,500\, {\rm K} \le T \le 25\,000$| K. Here, L and M are the luminosity and the mass in solar units, respectively, Z is metallicity, T is the effective temperature of the star, and V = v_∞/v_esc = 1.3 is the ratio of the wind velocity at infinity to the escape velocity from the star.

For hotter stars with |$25\,000\, {\rm K} \le T \le 50\,000$| K,

with V = 2.6.

For LBVs beyond the HD limit (Humphreys R. M., Davidson K.1994) (L > 6 × 10⁵ and 10⁻⁵RL^0.5 > 1.0, where R is the stellar radius in solar units),

where f_LBV = 1.5 is a calibration factor. Obviously, lower f_lbv results in weaker LBV wind, leaving a heavier remnant (Belczynski et al. 2010).

The reduced Wolf–Rayet star mass-loss with small H-envelope mass takes the form of metallicity-dependent power law,

with

m = 0.86 describing the dependence of wind mass-loss on metallicity, and M_He the He core mass of the star (Vink & de Koter 2005).

For other stars, we use the wind prescriptions described in Hurley et al. (2000).

2.2 Black hole formation

Stars are powered by burning their core fuels to heavier elements step by step. For massive stars, this process continues until an iron core is built up in the stellar centre. As the fusion of iron does not produce further energy, burning halts. Later, stars contract on their own weights, leading to accelerating processes of electron capture and core element dissociation. These processes dramatically reduce the pressure that should have resisted their self-gravity, triggering a runaway core-collapse process. Collapse halted by nuclear forces and neutron degeneracy pressure and form a proto-NS. Explosion launches after the ‘bounce’ of the core, and part or all of the expelled stellar envelope will fall back and accrete on to the proto-NS, which may eventually collapse into a BH (Fryer et al. 2012).

There are still many uncertainties associated with the physics of the SN mechanism. Here, we use the delayed SN prescription of Fryer et al. (2012; hereafter F12-delayed), where the explosion did not launch until over ∼250 ms after the collapse. For a massive star with pre-SN mass of M_SN and CO core mass of M_CO, the expected BH remnant mass M_BH is estimated as follows:

where M_BH,bar is the baryonic mass,²M_proto is the proto-NS mass after core collapse,

and M_fb is the amount of material falls back to the proto-NS,

where |$f_1=0.133 - {0.093 \over M-M_{\rm proto}}$|⁠, and f₂ = −11f₁ + 1. We define f_fb  = M_fb/(M_SN − M_proto) as the fallback fraction during the BH formation, which is important in determining the BH’s natal kick in some kick prescriptions.

Stars with He core mass M_He in the range of |${\sim} 35{\!-\!}60\, {\rm M}_{\odot }$| are subjected to PPISNe (Heger & Woosley 2002; Yusof et al. 2013; Belczynski et al. 2016a; Leung, Nomoto & Blinnikov 2019; Marchant et al. 2019; Stevenson et al. 2019), with most of the mass above the core stripped by a set of pulsations, leaving behind the BHs with mass significantly smaller than they would be if only accounting for the core-collapse SNe. We adopt the prescription of PPISNe in Marchant et al. (2019), who computed an array of H-free metal-poor (0.1 Z_⊙) single-star models based on the standard ¹²C(α, γ)O¹⁶ reaction rate to evaluate the PPISN mass-loss, and the BH masses after PPISNe can be estimated as:

where ζ_i are the polynomial fitting coefficients of Marchant et al. (2019)’s PPISN prescription given by Stevenson et al. (2019; as listed in Table 2). Note that the remnant mass is a non-monotonic function of the initial stellar mass.

More massive stars with |$60\, {\rm M}_{\odot }\le M_{\rm He} \le 135\, {\rm M}_{\odot }$| are subjected to PISNe, the entire star is completely disrupted with no remnant left.³ Stars with |$M_{\rm He}\gt 135\, {\rm M}_{\odot }$| are assumed to directly collapse to BHs.

2.3 Supernova kicks

As BH natal kicks suffer from a lack of stringent constraints from both observation and theory (Willems et al. 2005; Fragos et al. 2009; Repetto, Davies & Sigurdsson 2012; Repetto & Nelemans 2015; Belczynski et al. 2016c; Mandel 2016; Repetto, Igoshev & Nelemans 2017), we adopt three different natal kick prescriptions (kick_F) as follows (Banerjee et al. 2020):

1. Standard fallback-controlled kick (hereafter k1)

   The BH natal kick velocities v_kick,BH are scaled linearly with the NS natal kick velocities v_kick,NS by a factor (1 − f_fb) (Fryer et al. 2012; Giacobbo, Mapelli & Spera 2018a),

   $$\begin{eqnarray} v_{\rm kick,BH} = v_{\rm kick,NS}(1-f_{\rm fb}). \end{eqnarray}$$

   (9)

2. Convection-asymmetry-driven natal kick (hereafter k2)

   The BH natal kicks are produced by the convection asymmetries of the collapsing SN core (Scheck et al. 2004; Fryer & Kusenko 2006), so

   $$\begin{eqnarray} v_{\rm kick,BH} = \left\lbrace \begin{array}{ll}v_{\rm kick,NS}\frac{\lt M_{\rm NS}\gt }{M_{\rm BH}}(1-f_{\rm fb}) & {\rm if}\ M_{\rm CO} \le 3.5\, {\rm M}_{\odot },\\ k_{\rm conv}v_{\rm kick,NS}\frac{\lt M_{\rm NS}\gt }{M_{\rm BH}}(1-f_{\rm fb}) & {\rm if}\ M_{\rm CO}\ \gt 3.5\, {\rm M}_{\odot }.\\ \end{array} \right. \end{eqnarray}$$

   (10)

   In this equation, k_conv is an efficiency factor (somewhere between 2 and 10, and we set k_conv  = 5 here), and <M_NS > is a typical NS mass, taken to be |$1.4\, {\rm M}_{\odot }$|⁠.

3. Neutrino-driven natal kick (hereafter k3)

   The BH natal kicks are produced through asymmetric neutrino emission (Fuller et al. 2003; Fryer & Kusenko 2006),

   $$\begin{eqnarray} v_{\rm kick,BH} = v_{\rm kick,NS}\frac{ \min (M_{\rm eff},M_{\rm BH})}{M_{\rm BH}}, \end{eqnarray}$$

   (11)

where M_eff (usually between |$5$| and |$10\, {\rm M}_{\odot }$|⁠) is the effective remnant mass, and we let |$M_{\rm eff} = 7\, {\rm M}_{\odot }$| (Banerjee et al. 2020).

To constrain the allowed velocity range in different kick prescriptions, we take a flat distribution of v_kick,NS in the range of |$0{\!-\!}1000\rm \, kms^{-1}$|⁠, and assume that the SN kicks are isotropically distributed and the mass is instantaneously lost at the moment of SN. Then v_kick,BH can be obtained from the equations mentioned above for different kick prescriptions. When we calculate the merger rate density we adopt a more realistic, predetermined Maxwellian distribution for the NS kick velocity.

Note also that in both k1 and k2 prescriptions, there is no natal kick for BHs formed through direct core collapse (f_fb = 1.0).

2.4 Natal BH spins

The BH spins are modelled following Tanikawa et al. (2021b). We assume zero natal spin of the zero-age main sequence (ZAMS) star, which then evolves due to stellar evolution, stellar winds, and binary interactions. Newborn BHs inherit their progenitor’s spin angular momenta, except for the PPISN events where zero BH spin parameters are assumed. If the spin angular momenta of the BH progenitors are larger than those of extreme Kerr BHs, the BH spin parameters are forced to be unity. Thus, the BH spin parameter can be expressed as:

where |$\overrightarrow{S}$| and M are the spin angular momentum and the mass of the BH progenitor just before its collapse, respectively, c the speed of light, G the gravitational constant, and |$\overrightarrow{L}$| the binary orbital angular momentum.

The BH natal kicks would tilt |$\overrightarrow{\chi }$| from |$\overrightarrow{L}$|⁠. We choose the coordinate in which the z-axis is parallel to the orbital angular momentum vector just before the second BH formed, that is, the normalized orbital angular momentum vector is (0, 0, 1). Then, the normalized spin vectors of the first and second formed BH |$\overrightarrow{\chi _1}$| and |$\overrightarrow{\chi _2}$| can be written as:

where |$\theta _1^{\prime }$| is the angle between the first BH spin vector and binary orbital angular momentum vector just before the second BH forms, |$\phi _1^{\prime }$| is randomly chosen between 0 and 2π. Finally, the angles θ₁ (θ₂) between the first (second) formed BH spin and the final BBH orbital angular momentum vector just after the second BH formation can be expressed as:

and we do not consider possible BH spin alignment with orbital angular momentum due to tides or mass transfer. The effective spin parameter χ_eff of merging BBHs which reflects the spin-orbit alignment is defined as:

where m₁ and m₂ are the merging BBH masses.

2.5 Common-envelope evolution

For semidetached binaries, a common-envelope (CE) phase occurs when the mass transfer becomes dynamically unstable or when the two stars (with at least one of them being a giant-like star) collide at orbital periastron before Roche lobe overflow (RLOF; Hurley et al. 2002). The CE phase plays a fundamental role in the formation of GW190521-like systems. Due to the extreme initial mass ratio q ∼ 0.1, mass transfer from the massive primary star to the secondary star is usually dynamically unstable and leads to CE evolution. In addition, systems with large orbital eccentricities are likely to collide at periastron. Binaries can survive the CE phase if the accretor’s orbital energy is large enough to unbind the stellar envelope.⁴ Here, we adopt the α_CEλ formalism (de Kool 1990), which can be expressed as:

where the envelop’s binding energy,

Here, α_CE is the efficiency of converting the released orbital energy to eject CE, λ the binding energy parameter depending on envelope’s structure (see Ivanova et al. 2013, for details), M₁, M_c, and M_env  = (M₁ − M_c) the masses of the donor, donor’s core, and envelope, respectively, M₂ the mass of the accretor, a_i,CE and a_f,CE the binary separation before and after the CE phase, respectively, R_RL the donor’s RL radius at the onset of the CE phase, U the specific internal energy (including both thermal and recombination energies) of envelope, and α_th the efficiency with which thermal energy can be used to eject the envelope. In this work, α_th = 1 is assumed.

Following Xu & Li (2010) and Wang, Jia & Li (2016), we calculate the binding energy parameter λ for stars more massive than |$60\, {\rm M}_{\odot }$| with metallicity Z = 0.02, 0.001, and 0.0001 using the stellar evolution code mesa (version 11701; Paxton et al. 2011, 2015, 2018, 2019). We then calculate λ_b and λ_g with E_bind including and excluding the internal energy term U, respectively. Detailed models and fitting results of λ are presented in Appendix  A. We only use λ_b in our following population synthesis calculations.

In order to explore the dependence of our results on α_CE, we have run a set of simulations with α_CE = 0.5, 1.0, and 3.0 (α_CE > 1.0 would occur when additional energy or angular momentum depositing into the giant’s envelope is considered, see e.g. Soker 2004).

2.6 Population synthesis

We perform BPS simulations of binary stars with the bseemp code. We assume that all the stars are in binaries, and exclude binary systems with at least one of the two components fills its RL at the beginning of evolution. The initial masses M_i,1 of the primary stars are distributed in the range of |$[300:900]\, {\rm M}_{\odot }$| following the Kroupa (2001) law. For the masses of the secondary stars M_i,2 = qM_i,1, we adopt a flat distribution of the mass ratio q (Sana et al. 2012), and limit |$M_{\rm i,2}= [10:60]\, {\rm M}_{\odot }$|⁠. We consider various metallicities with Z = 0.0002, 0.0004, 0.0008, 0.0016, 0.0032, 0.0063, 0.0126, and 0.02. The initial orbital semimajor axis a_i is assumed to be distributed uniformly in logspace and restricted to |$[3:10^{7}]\, {\rm R}_{\odot }$|⁠. We consider both initially circular (hereafter ‘e_i = 0’ model) and eccentric (hereafter ‘e_i = 0 ∼ 1’ model) orbit configurations in each case. In the latter model, the eccentricity follows a uniform distribution between 0 and 1.

We focus only on GW190521-like systems with the primary and secondary BH masses |$m_1 = [150:180]\, {\rm M}_{\odot }$| and |$m_2 = [10:20]\, {\rm M}_{\odot }$|⁠, so it is reasonable to only simulate binaries with a limited initial parameter ranges, and make sure that the simulated population parameters are complete to form GW190521-like systems.

Incorporating the three kick prescriptions (kick_F  = k1, k2, and k3) and three values of α_CE = (0.5, 1.0, and 3.0), we perform 18 sets of BPS simulations of 10⁷ primordial binaries at each metallicity. We list the initial parameters of models with Z ≤ 0.0016 in Table 3. For higher metallicity, we found that there is no GW190521-like system formed in our simulation.

We assume a non-conservative mass transfer prescription with the accretion efficiency being 0.5. We follow the evolution of the primordial binaries until the formation of BBHs. For a newborn BBH system with component masses m₁ and m₂, orbital semimajor axis a₀ and eccentricity e₀, the inspiral time delay t_inspiral (a₀, e₀), namely the time elapsed between the birth and the merger of the BBH, can be calculated (Peters 1964):

where

and

2.7 Merger rate density

We then follow the procedure of Giacobbo & Mapelli (2018b) to estimate the cumulative merger rate density |$\mathcal {R}(z\le z_{\rm det})$| of GW190521-like BBH systems within a given redshift z_det,

where t_lb(z) is the look-back time for binaries formed at redshift z, |$M_{*} \simeq 0.55\, {\rm M}_{\odot }$| the mean mass of a stellar system in population with the Kroupa (2001) IMF, f_bin = 0.7 the fraction of stars in binaries (Sana et al. 2012), W_b the contribution of specific binaries we are interested in (Hurley et al. 2002), |$\mathcal {SFR}(z)$| the cosmic star formation rate density as a function of z, usually peaked at z ∼ 1.9 and declined exponentially at later time. We assume that in our models the star formation commenced at z = 15.

The BBH progenitor binaries formed at redshift z_f would merge as GW sources at z_m (z_m < z_f) after a delay time t_delay, which is defined as the interval between the formation of the progenitor binary and the coalescence of the BBH, i.e. t_delay = t_inspiral + T, where |$T\lt 10 \,\,\rm Myr$| is the lifetime of progenitor system, if the orbital angular momentum loss is efficient enough. So, we can get the look-back times at their formation

where E(z) = [Ω_m(1 + z)³ + Ω_λ]^1/2, and at their merger

In our calculation, we employ the flat |$\rm \Lambda CDM$| model with |$H_0=67.8\, \rm kms^{-1}Mpc^{-1}$|⁠, Ω_m = 0.3, and Ω_λ = 0.7, where |$\tau _{\rm H} = 1/H_0=14.4\, \rm Gyr$| is the Hubble time (Planck Collaboration 2016). We adopt the cosmic |$\mathcal {SFR}(z)$| density in Madau & Dickinson (2014):

and the metallicity as a function of redshift z in Belczynski et al. (2016b). For the portions of distributions extending beyond the metallicity range [0.0002, 0.02], we use the recorded information of the systems at Z = 0.0001 or 0.02. We exclude the mergers in the near future (t_merg < 0).

2.8 Character strain

The characteristic strain of the GW signals at the nth harmonic can be calculated following Kremer et al. (2019),

where |$M_{c,z} = M_c(1+z) = \frac{(m_1m_2)^{3/5}}{(m_1+m_2)^{1/5}}(1+z)$| is the observed chirp mass at redshift z, and D_L is the luminosity distance to the source calculated by

|$f_{n,z}=\frac{f_n}{1+z}=\frac{nf_{\rm orb}}{1+z}$| is the observed frequency of the nth harmonic (f_n is the frequency of the nth harmonic in the source frame and f_orb is the source-fame orbital frequency), g(n, e) is the function of eccentricity, and F(e) is the eccentricity correction factor defined to be (Peters & Mathews 1963):

As the GW power is sharply peaked at the peak frequency f_peak (Peters & Mathews 1963), we calculate the characteristic strain of our modelled GW sources at the peak frequency f_peak for simplicity (Hamers 2021),

Thus, n = n_peak = f_peak/f_orb (Wang, Tanikawa & Fujii 2022).

During the inspiral, the eccentricity changes due to gravitational radiation (Peters 1964),

and the orbital separation evolves with eccentricity

where c₀ is determined by the initial condition a(e₀) = a₀ (see equation 19).

3 RESULTS

For each model, we regard the merging BBHs with component masses |$m_1=[150:180]\, {\rm M}_{\odot }$| and |$m_2=[10:20]\, {\rm M}_{\odot }$| as GW190521-like systems (as shown in Fig. 1). According to the analysis of Nitz & Capano (2021), the primary BH mass in GW190521 is |${\sim} 170\, {\rm M}_{\odot }$|⁠. We simulate binary evolution with the primary mass |${\le} 900\, {\rm M}_{\odot }$| at different metallicities and find that such massive BHs can only form at relatively low metallicities (Z ≤ 0.0016), with the pre-collapse core-helium masses heavier than |$135\, {\rm M}_{\odot }$|⁠, while stars at higher metallicity (Z ≥ 0.0032) will undergo PISNe, triggering complete disruption of the star and leaving no compact remnants. This is because stars formed in lower metallicity environments can reach higher central temperatures, which results in larger core masses than their counterparts at higher metallicity. So, our following discussions only refer to the data set of the models with Z ≤ 0.0016.

Table 4 lists the predicted numbers of GW190521-like systems and kick velocity distributions under different conditions. Their main features can be summarized as follows:

• Most of GW190521-like systems form through the ‘MT + CE’ channel. Here, ‘MT’ means that stars in a binary interact via stable RLOF or wind accretion, and the ‘CE’ phase is triggered by eccentric collision of both stars at periastron instead of dynamically unstable mass transfer caused by the expansion of the donor star. Binary stars in eccentric orbits may collide at periastron before either one fills its RL, and such collisions lead to CE evolution if at least one of the stars is a giant-like star (Hurley et al. 2002). Thus, only the k3 kick prescription that produces non-zero BH1 natal kick velocities works in this channel. Besides, the number of BBH mergers decreases as α_CE increases, because larger α_CE leads to wider orbits after the CE phase, making it more difficult for the BBHs to merge. To reproduce GW190521-like systems in this channel requires |$v_{\rm kick,1} \simeq 0{\!-\!}50 \, \rm kms^{-1}$| and |$v_{\rm kick,2} \simeq 0{\!-\!}700 \, \rm kms^{-1}$|⁠.

• Only about 0.1 per cent of GW190521-like systems form through ‘MT + MT’ channel. The systems with moderate low eccentricities can avoid collision at periastron until the BBH formation. So, this channel is independent of the value of α_CE. The predicted natal kick velocities are v_kick,1 = 0 and |$v_{\rm kick,2} \simeq 30{\!-\!}90 \, \rm kms^{-1}$| (kick_F = k1, k2), and |$v_{\rm kick,1} \simeq 0{\!-\!}50 \, \rm kms^{-1}$| and |$v_{\rm kick,2} \simeq 16{\!-\!}400 \, \rm kms^{-1}$| (kick_F = k3).

• For very low mass ratio binaries (q < 0.1), the secondary star usually does not have enough energy to drive off the CE if it is triggered by dynamically unstable MT, so there is no system formed via the ‘CE + MT’ channel in the ‘e_i = 0’ model. In the ‘e_i = 0 ∼ 1’ model, the non-zero eccentricity makes CE evolution possible just like in the ‘MT + CE’ channel. The number of surviving systems also decreases with increasing α_CE. The predicted natal kick velocities are v_kick,1 = 0 and |$v_{\rm kick,2}\simeq 100{\!-\!}400 \, \rm kms^{-1}$| (kick_F = k1, k2), and |$v_{\rm kick,1}\simeq 0{\!-\!}45 \, \rm kms^{-1}$| and |$v_{\rm kick,2} \simeq 60{\!-\!}420 \, \rm kms^{-1}$| (kick_F = k3).

• No GW190521-like system form though the ‘CE + CE’ channel.

Table 5 presents the inferred parameters of GW190521-like systems that will merge within z = 1.1 and their progenitors. As most of them are formed with kick_F = k3, we only show the results with the k3 prescription. The number and its subscripts and superscripts represent the 50th, 16th, and 84th percentiles of each parameter. It is seen that the results do not show significant differences in the ‘e_i = 0’ and ‘e_i = 0 ∼ 1’ models.

The analysis of Abbott et al. (2020a, b) suggested that GW190521 merged at the redshift of |$0.82_{-0.34}^{+0.28}$|⁠, while Nitz & Capano (2021) predicted a luminosity distance of |$1.06_{-0.28}^{+1.4}\rm Gpc$| (⁠|$z\simeq 0.21_{-0.05}^{+0.23}$|⁠). According to their restrictions on the redshift, we display the calculated merger rate density |$\mathcal {R}$| of GW190521-like systems at z ≤ 0.48 and z ≤ 1.1 in Table 6, which lie in the range of |$4\times 10^{-5}{\!-\!}5\times 10^{-2} \, \rm Gpc^{-3}yr^{-1}$|⁠. As mentioned above, the merger rate density with the k3 kick prescription is much more than with the k1 or k2 prescription.

Fig. 2 shows the distribution of e₀ and a₀ at the birth of the BBHs (with Z = 0.0002 and kick_F = k3). The upper and lower panels correspond to the ‘e_i = 0’ and ‘e_i = 0 ∼ 1’ models, and the left-, middle-, and right-hand panels correspond to α_CE = 0.5, 1.0, and 3.0, respectively. The distribution of e₀ tends to be wider with increasing α_CE in both ‘e_i = 0’ and ‘e_i = 0 ∼ 1’ models. For α_CE = 3.0, there are many BBHs with e₀ > 0.8, while for α_CE  = 0.5 and 1.0, few BBHs form with extremely eccentric and wide orbits (a₀ ∼ 10⁴R_⊙). Most of the BBHs have e₀ ≤ 0.4 and a₀ concentrated within ∼10–100 R_⊙. BBHs with shorter orbits and moderate eccentricities usually merge earlier than others, as shown by the colour bars.

Figs 3–5 show the results of Z = 0.0004, 0.0008, and 0.0016, respectively. They reflect the similar tendency as in Fig. 2. A comparison of Figs 2–5 shows that, as the metallicity increases, there are more BBHs with large eccentric and wide orbits (e₀ > 0.8 and a₀ > 100 R_⊙). This is because the progenitors of GW190521-like system at higher Z are more massive, and thus have larger size than those at low Z.

Current and upcoming missions such as the ground-based aLIGO, Cosmic Explorer (Reitze et al. 2019), Einstein telescope (ET; Punturo et al. 2010), and the space-borne DECi-hertz Interferometer Gravitational Wave Observatory (DECIGO; Seto, Kawamura & Nakamura 2001), and Laser Interferometer Space Antenna (LISA; Amaro-Seoane et al. 2017) would detect thousands of merger events of BBHs per year (Evans et al. 2021).We explore whether GW190521-like systems could be detected by these instruments. In Fig. 6, the left three panels present the evolution of the eccentricity of GW190521-like systems during the inspiral prior to merge as a function of the peak frequency. The GW emission dominates the evolution of the binary semimajor axis and eccentricity, leading to efficient circularization of the BBH systems before the merger. The right three panels present the GW signal characteristic strain at the n_peakth harmonic along with the orbital evolution as a function of the peak frequency, which is the key ingredient determining whether such mergers can be seen with the GW detectors, in the ‘e_i = 0’ model. The evolutionary tracks (calculated with equations 29 and 30) are gradually overlapped by the sensitivity curves of LISA, DECIGO, ET, Cosmic Explorer, A + LIGO, and aLIGO during the orbital shrinking and circularizing stage. The systems become largely circularized before entering the sensitivity band of ET, Cosmic Explorer, A + LIGO, and aLIGO, and any residual eccentricity is expected to have a negligible effect on their detectability (Mandel et al. 2008). The most distant detectable GW190521-like mergers are at the redshfit ∼4.4. Fig. 7 shows the results in the ‘e_i = 0 ∼ 1’ model, and there are no significant differences compared with Fig. 6.

Another important characteristic of merging BBHs is their spins, which get imprinted in the GW signal (Cutler, Flanagan & É. 1994). The BH progenitors gain and lose their spin angular momenta through stellar evolution, mass transfer, and tidal interactions (Hurley et al. 2002; Belczynski et al. 2020b; Tanikawa et al. 2021b). The spin angular momenta are generally parallel to the orbital angular momentum, until the kicks to the BHs cause their spin axes tilted. According to traditional tidal theory (Zahn 1977; Hut 1981), the torque depends on the ratio of the stellar radius R to the separation a of both stars, that is, ∝(R/a)⁶. Because the progenitor of BH1 is very massive and the initial orbit is very wide, spinning up of the BH1’s progenitor is ineffective, so the merging BBHs generally have small |χ_eff|. In the ‘e_i = 0’ model, χ_eff = 0 ∼ 0.1 (kick_F = k3) and −0.09 ∼ 0.1 (kick_F = k1, k2); in the ‘e_i = 0 ∼ 1’ model, χ_eff = −0.08 ∼ 0.1 (kick_F = k3) and −0.1 ∼ 0.1 (kick_F = k1, k2). They are in contradiction with Nitz & Capano (2021)’s prediction that the BH spin is anti-aligned with the orbital angular momentum and |$\chi _{\rm eff}=-0.51_{-0.11}^{+0.24}$|⁠. However, the mechanism of tidal interactions is not well understood. If we adopt the Geneva model (Eggenberger et al. 2008) in which angular momentum is mainly transported by meridional currents (see also Belczynski et al. 2020a), |$|\overrightarrow{\chi _1}|$| can increase from ∼0 to ∼0.25. As the spin of the merger product is dominated by the contribution of the more massive BH1, the estimated χ_eff changes to be −0.3 ∼ 0.32. We also note that, by using the Tayler–Spruit magnetic dynamo angular transport, Belczynski et al. (2020a) inferred the natal spins of BBHs (m₁ = 84.9 M_⊙, m₂ = 64.6 M_⊙) that would merge within Hubble time to be |$|\overrightarrow{\chi _1}| = 0.052$| and |$|\overrightarrow{\chi _2}| = 0.523$|⁠.

4 DISCUSSION

The natal kick plays a key role in the life of a compact star binary, as it affects not only the orbital parameters and systemic velocity, but also the binary evolutionary path (Brandt & Podsiadlowski 1995). There is a general consensus that NSs are usually born with large kick velocities |$v_{\rm kick,NS}\sim 200{\!-\!}500\, \rm kms^{-1}$| (Lyne & Lorimer 1994). However, the origin of the SN kick is in debate. One possible mechanism is the asymmetric material ejection during the SN explosion, triggered by the large-scale hydrodynamic perturbation or convection instabilities in the SN core (Burrows & Hayes 1996; Goldreich, Lai & Sahrling 1997; Scheck et al. 2004; Nordhaus et al. 2012; Gessner & Janka 2018). Other investigations suggest that it may be related to the anisotropic neutrino emission from the proto-NS induced by strong magnetic field (Kusenko & Segrè 1996; Lai & Qian 1998; Maruyama et al. 2011). Besides, the topological current may be responsible for the natal kick (Charbonneau & Zhitnitsky 2010).

The kick velocity distribution is also in active study. Arzoumanian, Chernoff & Cordes (2002) studied the velocity distribution of radio pulsars based on large-scale 0.4 GHz pulsar surveys, and found a two-component velocity distribution with characteristic velocities of 90 and 500 kms⁻¹. Hobbs et al. (2005) analysed a catalogue of 233 pulsars with proper motion measurements, and suggested the NS natal kick distribution with a Maxwellian one-dimensional dispersion |$\sigma _{\rm NS}=265\, \rm kms^{-1}$|⁠, which is widely used in later studies. From the analysis of the proper motions of 28 pulsars using very long baseline array interferometry data, Verbunt, Igoshev & Cator (2017) showed that a distribution with two Maxwellians improves significantly on a single Maxwellian for the young pulsar velocities.

Whether stellar-mass BHs receive such large kicks is also a matter of debate. A growing number of studies have been devoted to investigate the natal kicks of BHs relying on a variety of methods and data sets, such as the study of massive runaway and walkaway stars (Blaauw 1961; De Donder, Vanbeveren & van Bever 1997; Renzo et al. 2019; Aghakhanloo et al. 2022), BH X-ray binaries (Mirabel et al. 2001; Jonker & Nelemans 2004; Repetto et al. 2012; Wong et al. 2012, 2014; Atri et al. 2019; Kimball et al. 2022), astrometric microlensing (Andrews & Kalogera 2022), and merging BBH GW events (Abbott et al. 2021; The LIGO Scientific Collaboration 2021).

In light of the observational constraints on the NS/BH natal kick velocities, several phenomenological and analytical kick prescriptions are proposed, mainly depending on the SN ejecta mass and remnant mass (e.g. Bray & Eldridge 2018; Giacobbo & Mapelli 2020; Mandel & Müller 2020; Richards et al. 2023). Because the kick-induced orbital eccentricity determines the time-scale over which BBHs are expected to merger via GW radiation, the merging history of BBHs provide a probe to the natal kick received by BHs. Based on the premise that GW190521 is an IMRI with component masses of |${\sim} 170$| and |${\sim} 16\, {\rm M}_{\odot }$| (Nitz & Capano 2021), we examine the isolated binary evolution channel with three kick prescriptions. In the k1 and k2 prescriptions, the BH natal kick is determined by the fallback fraction f_fb so massive BH experienced totally fallback (f_fb = 1.0) would receive no kick, while in k3 prescription BHs always receive a kick produced through asymmetric neutrino emission.

Our calculations indicate that, to produce the merger event, the less massive BH should receive a natal kick with velocity of a few hundred |$\rm kms^{-1}$|⁠, thus preferring the k3 prescription. This is of particular interest since in most cases both BHs formed through totally fallback, and the conclusion is not sensitive to the choice of the CE efficiency α_CE.

We predict the merger rate density of GW190521-like systems |$\mathcal {R}(z\le 1.1)\sim 4\times 10^{-5}{\!-\!}5\times 10^{-2} \, \rm Gpc^{-3}yr^{-1}$| if the BH natal kick is weighted to follow a Maxwell distribution of the NS kick with σ_NS = 265 kms⁻¹. Under the interpretation that GW190521 is an almost equal mass ratio system, LIGO/Virgo collaboration reported the merger rate density of GW190521-like systems to be |$0.13_{-0.11}^{+0.30}\, \rm Gpc^{-3} yr^{-1}$| with the effective spin parameter |$\chi _{\rm eff}=0.08_{-0.36}^{+0.27}$| (Abbott et al. 2020a, b). By employing a new estimate of the PPISN mass-loss, Belczynski et al. (2020a) obtained a merger rate density of |${\sim} 0.04\, \rm Gpc^{-3} yr^{-1}$| for such events via isolated binary evolution. Tanikawa et al. (2021b) estimated the merger rate density of Population III BBHs (with total mass |${0\sim} 130{\!-\!}260\, {\rm M}_{\odot }$| and composing at least one |$130{\!-\!}200\, {\rm M}_{\odot }$| IMBH) about 0.01 |$\rm Gpc^{-3}yr^{-1}$|⁠. Hijikawa et al. (2022) performed a BPS calculation for very massive Population III stars and derived the property of the BBH mergers, adopting constant values for α_CEλ in their CE evolution. In their ‘low mass + high mass’ model, the resultant compact binaries consist of a stellar mass BH (below the PISN mass gap) and an IMBH (above the PISN mass gap) with mass ratio ranging from 0.15 to 0.35. The predicted merger rate density peaks at z ∼ 10 with a value of |$(1-10)\, \rm Gpc^{-3}yr^{-1}$|⁠, and declines to nearly zero at z ≤ 3 because of the very short delay time (less than 10 |$\rm Myr$|⁠). In our Population II evolution channel, the merger rate peaks at z ∼ 2 with the delay time ranging from |${\sim} 1.4{\!-\!}12.1\, \rm Gyr$|⁠.

5 SUMMARY

The third observing run operated by aLIGO and advanced Virgo discovered a massive BBH merger event GW190521, with a remnant total mass of |$150_{{-17}}^{+29}\, {\rm M}_{\odot }$|⁠, falling in the IMBH regime (Abbott et al. 2020a), and the component masses were estimated to be |$(m1, m2) = (85_{{-14}}^{+21}$|⁠, |$66_{{-18}}^{+17})\, {\rm M}_{\odot }$| within 90 per cent credible region (see also Barrera & Bartos 2022; Gamba et al. 2023). Nitz & Capano (2021), however showed that GW190521 may be alternatively an IMRI, with the component masses of m₁ ∼ 170 M_⊙ and m₂ ∼ 16 M_⊙, which happen to straddle the PISN mass gap. In the most recent analysis, Gamba et al. (2023) revealed the BH masses to be |$81_{{-25}}^{+62}\, {\rm M}_{\odot }$| and |$52_{{-32}}^{+32}\, {\rm M}_{\odot }$| under the hypothesis that it was generated by the merger of two non-spinning BHs on hyperbolic orbits. So, the nature of GW190521 is still uncertain.

Assuming the configuration of Nitz & Capano (2021) for GW190521 (or similar systems to be discovered in the future), we perform BPS simulation to interpret the formation of GW190521-like systems via isolated binary evolution channel. Our analyses prefer that this merger event had evolved from primordial binary systems in metal-poor environment with Z ≤ 0.0016. The majority of them are formed via an initial phase of stable RLOF before the formation of the BH1, followed by a CE phase triggered by collision at periastron when BH1’s companion is a giant-like star in a close eccentric orbit. The initial ZAMS progenitor masses are expected to be |$M_{\rm i,1}\sim 300{\!-\!}800\, \rm {\rm M}_{\odot }$| and |$M_{\rm i,2}\sim 20{\!-\!}60\, {\rm M}_{\odot }$|⁠, respectively, which are metallicity-dependent. By using the fallback independent kick prescription, the merger event requires the primary and secondary BHs to receive natal kicks with velocities |$v_{\rm kick,1}\ \lt\ 50\, \rm kms^{-1}$| and |$v_{\rm kick,2} \ \lt\ 700\, \rm kms^{-1}$|⁠. Our results support the hypothesis that BHs formed by direct core collapse receive considerable large natal kick. The predicted merger rate density for GW190521-like systems is |$4\times 10^{-5}{\!-\!}5\times 10^{-2} \, \rm Gpc^{-3}yr^{-1}$| at z_m ≤ 1.1. We also find that, using the traditional treatment of tidal interaction results in very small effective spin parameter, but if using the Geneva model instead, χ_eff ranges from −0.3 to 0.32, roughly located within the interval |$-0.51^{-0.11}_{+0.24}$| estimated by Nitz & Capano (2021).

ACKNOWLEDGEMENTS

We thank the anonymous referee for useful comments, which helped improve the manuscript. We are also grateful to Shi-Jie Gao for help with the calculation of λ. This work was supported by the Natural Science Foundation of China under grant no. 12041301and12121003, and the National Key Research and Development Program of China (grant no. 2021YFA0718500). We acknowledge use of the following python packages: astropy (Astropy Collaboration 2013), matplotlib (Hunter 2007), numpy (van der Walt, Colbert & Varoquaux 2011), and scipy (Virtanen et al. 2020).

DATA AVAILABILITY

All data underlying this article will be shared on reasonable request to the corresponding authors.

Footnotes

References

APPENDIX A: THE CE PARAMETER λ FOR VERY MASSIVE STARS

In our stellar models, we neglect any stellar rotation. For stellar winds, we use the Vink et al. wind prescription as mentioned in Section 2.1, but the result in Nugis & Lamers (2000) for Wolf–Rayet stars. The convective regions are determined by the Schwarzschild criterion using α_MLT = 1.73 (MacLeod et al. 2017), and we assume an exponential overshoot (Herwig 2000) by 1 per cent of the scale height.

Figs A1 and A2 show the evolution of λ_b(the thick solid lines) and λ_g(the thin solid lines) with stellar radius for different masses at Z = 0.02, 0.001, and 0.0001. We make polynomial fitting of the results with

where R is the stellar radius. Tables A1 and A2 provide the coefficients n_{i(i = 0, 10)} for stars with selected masses. The goodness of fit χ² ≃ 0.998.