Detection Prospects of Fast-merging Gravitational Wave Sources in M31

Except GW transients, no DCO systems have been reported outside the Milky Way. Here we present estimations on the potential population of LISA sources in M31 by adopting a BPS approach. Only the systems formed via isolated binary evolution are included in our calculations. We employ the BSE code originally developed by Hurley et al. (2002) and significantly modified by Shao & Li (2014) that has been used to predict the detection of Galactic LISA sources (Shao & Li, 2021). In this code, some important modifications are summarized as follows. During the evolution of primordial binaries, it is assumed that the accretion rate of the secondary stars is dependent on their rotational velocities. Under this assumption, small amount of transferred mass is accreted by rapidly rotating secondary stars. Mass accretion onto the secondary stars can cause them to expand and spin up. On the one hand, the significant expansion of the secondary stars due to rapid mass accretion may lead to the formation of contact binaries as the secondary stars fill their Roche lobes (Nelson & Eggleton, 2001). On the other hand, a small amount of accreted mass can spin up the secondary stars to near critical rotation (Packet, 1981). As a consequence, the secondary stars are likely to reach critical rotation at the early stage of mass transfer and then accrete at relatively low rates. In this situation, the maximal mass ratios between the primary stars and the secondary stars for avoiding the formation of contact binaries can reach as high as about 6 (Shao & Li, 2014). For stellar wind mass-loss rates, the prescription of Belczynski et al. (2010b) has been applied to various types of stars, except that for helium stars we use half of the Hamann et al. (1995) rates (Kiel & Hurley, 2006). When dealing with mass-transfer stability in the binaries with a BH accretor and a massive donor, the criteria of Shao & Li (2021) have been adopted to determinate whether such binaries experience CE evolution or stable mass transfer. These criteria are related to the masses of both components and the radii of donor stars, and show that the binaries with heavier BHs are more likely to undergo stable mass transfer. It’s worth noting that these criteria are not strongly dependent on metallicities, except for the systems with donors initially more massive than $\sim 40M_{\odot}$¹¹1The stability of mass transfer is sensitive to the structure of the donor stars (Soberman et al., 1997). At solar metallicity, the donor stars with masses above $\sim 40-60M_{\odot}$ never develop a fully convective envelope which can increase the stability of mass transfer (Pavlovskii et al., 2017; Klencki et al., 2021).. Since the massive donors with increasing metallicities are more likely to experience extensive wind mass loss before mass transfer, the mass transfer tends to be more stable. The use of these criteria allows the formation of tight DCO binaries via stable mass transfer rather than CE evolution (see also e.g., Pavlovskii et al., 2017; van den Heuvel et al., 2017; Marchant et al., 2021).

During a CE phase, we use the standard energy conservation equation (Webbink, 1984) to calculate the orbital evolution of a binary system, i.e.,

where

is the binding energy of the donor’s envelope at the onset of CE evolution. Here $M_{\rm d}$ and $M_{\rm a}$ respectively denote the donor’s and accretor’s masses, $a$ is the orbital separation of the binary system, $M_{\rm d,env}$ is the mass of donor’s envelope that is ejected out of the binary, $R_{\rm d,L}$ is the Roche lobe radius of the donor at the beginning of mass transfer, the subscripts $\mathrm{i}$ and $\mathrm{f}$ correspond to the initial and final stages of a CE phase, respectively. $\lambda$ is the binding energy parameter that is determined by the mass distribution of the envelope (Xu & Li, 2010) and $\alpha_{\mathrm{CE}}$ is the CE ejection efficiency with the orbital energy of the binary unbinding donor’s envelope. When CE evolution is triggered in binaries with Hertzsprung gap donors, we allow them to survive CE phases (see e.g., Dominik et al., 2012).

Equation (1) indicates that all of binary orbital energy is used to eject the envelop of engulfing star when $\alpha_{\mathrm{CE}}=1$. Generally, the larger the $\alpha_{\mathrm{CE}}$, the binaries more likely to survive CE evolution and the longer the separations of post-CE binaries. Recently, Scherbak & Fuller (2023) suggested that $\alpha_{\mathrm{CE}}\sim 1/3$ according to the analyses of Galactic WDWD binaries. This is consistent with the case of $\alpha_{\mathrm{CE}}\sim 0.2-0.3$ for the post-CE binaries with a WD and a main-sequence companion (Zorotovic et al., 2010). Note that $\alpha_{\mathrm{CE}}\geq 3$ seems to be required to account for the formation of some post-CE binaries such as IK Peg with a $\sim 1.2M_{\odot}$ WD and a $\sim 1.7M_{\odot}$ main-sequence star (Davis et al., 2010). Furthermore, investigations on NS binaries indicate higher CE ejection efficiencies varying from $\alpha_{\mathrm{CE}}\sim 1$ (Zuo & Li, 2014) to $\alpha_{\mathrm{CE}}\sim 5$ (Fragos et al., 2019; Yarza et al., 2022). In our calculations, we set three values for CE ejection efficiencies i.e., $\alpha_{\mathrm{CE}}=0.3$, $\alpha_{\mathrm{CE}}=1.0$ and $\alpha_{\mathrm{CE}}=3.0$.

Following Shao & Li (2021), we apply three different mechanisms of supernova explosions. These mechanisms involve the rapid (Fryer et al., 2012), the delayed (Fryer et al., 2012) and the stochastic (Mandel & Müller, 2020) prescriptions. The rapid recipe allows the formation of a $2-5\mathrm{~{}M_{\odot}}$ mass gap among which no compact remnants are formed, while the delayed and stochastic mechanisms do not. So whether existing such a mass gap can provide a clue to constrain the process of supernova explosions (see e.g., Shao, 2022). Both the rapid and the delayed mechanisms treat the remnant masses as a function of CO core masses prior to supernova and fallback-material masses during explosion. In these two mechanisms, it is assumed that the natal kicks of NSs and BHs obey a Maxwellian distribution. We adopt a dispersion velocity of $\sigma=265\mathrm{~{}km}\mathrm{~{}s}^{-1}$ (Hobbs et al., 2005) for NSs formed via core-collapse supernovae. For BHs, we use reduced kick velocities due to the fallback of some materials compared to the case of NS kicks. In the stochastic mechanism, the remnant masses and the natal kicks are assumed to be probabilistic rather than deterministic. It is assumed that the natal kicks of both NSs and BHs follow a Gaussian distribution and a part of BHs form without any kick (Mandel & Müller, 2020). In this case, we further apply the scaling pre-factor calibrated by Kapil et al. (2023) to deal with the velocity distribution of NS kicks. In our calculations, we also include the formation of NSs through electron-capture supernovae. It is assumed that these NSs have the mass of $1.3M_{\odot}$ and their natal kicks follow a Maxwellian distribution with a relatively low dispersion velocity of $\sigma=30\mathrm{~{}km}\mathrm{~{}s}^{-1}$ (see also Shao & Li, 2014, 2018).

For primordial binaries we set a grid of initial parameters: the primary masses $M_{1}$ and the secondary masses $M_{2}$ are both limited in the range of $0.5M_{\odot}-100M_{\odot}$, and the orbital separations $a$ in the range of $3R_{\odot}-10000R_{\odot}$. The fraction of binaries among all stars is taken to be unity and all initial binaries are assumed to have circular orbits. Following Shao & Li (2021), we calculate the birth rate of a specific DCO system $j$ as

where $\rm SFR$ is the star formation rate of the M31 disk, $W_{\mathrm{{b}}}=\Phi(\ln M_{1})\varphi(\ln M_{2})\Psi(\ln a)\delta\ln M_{1}\delta\ln M_{2}\delta\ln a$ is a combination parameter to weight the contribution of the system evolved from a specific primordial binary. Then $\Phi(\ln M_{1})=M_{1}\xi(M_{1})$ represents the distribution of primary masses in logarithmic space and $\xi(M_{1})$ describes the initial mass function (IMF, Kroupa et al., 1993). Here we apply a modified IMF for M31 (Weisz et al., 2015), i.e.,

where $c_{1}=0.2074$ is a normalized parameter. This IMF can slightly increase the fraction of stars with masses above $0.5M_{\odot}$ compared to the Kroupa’s IMF. The secondary masses are assumed to have a uniform distribution (Kobulnicky & Fryer, 2007), then giving $\varphi\left(\ln M_{2}\right)=M_{2}/M_{1}$. It is assumed that the orbital separations follow a flat distribution in logarithmic space with $\Psi(\ln a)=c_{2}$, where $c_{2}=0.12328$ is a normalized parameter. In addition, the size of each interval $\delta\ln\chi$ in logarithmic space is determined as $\delta\ln\chi=\left(\ln\chi_{\rm max}-\ln\chi_{\rm min}\right)/(n_{\chi}-1)$. Here, $n_{\chi}$ represents the number of grid points for a specific parameter $\chi$ and is set to $n^{1/3}$ (where $n$ is the total number of primordial binary samples in the 3-dimensional parameter space mentioned above). Our calculations include 9 models with different supernova mechanisms and CE ejection efficiencies, and run $\sim 10^{7}$ primordial binaries for each evolutionary model.

The evolutionary history of SFR and stellar metallicity $\rm Z$ is an important factor that can make a considerable impact on the population of observable LISA sources in M31. We follow the work of Williams et al. (2017) which gives wide-area maps of the full star-formation history of the M31 disk by fitting resolved stellar photometry coming from the Panchromatic Hubble Andromeda Treasury survey. According to the result of PARSEC model (see Table 3 of Williams et al., 2017), we obtain the total SFR of M31 during different history stages. For the metallicity evolution we use the age-dependent model (Kemp et al., 2022) which linearly interpolates between (age (Gyr), $\rm Z$) coordinates of (0, 0.0002), (2, 0.008), (10, 0.02), (14, 0.022) approximately. The history distribution of $\rm SFR$ and $\rm Z$ of M31 has been mapped by Kemp et al. (2022, see their Figure 4).

Since the metallicity changes over the age, we need to divide all binary stars into different sample spaces according to the metallicity or age. For convenience, we set 12 grids of metallicities as shown in Figure 1. The orange ladder-like curve gives the sample spaces separated by different metallicities. The metallicity of every primordial binary is randomly taken from the 12 individual values. For each metallicity grid, we know the age and the SFR. When calculating $W_{\mathrm{{b}}}$, the corresponding number of the primordial binaries in each grid is then set to be $n/12$. In our calculations, we select the binary samples if the primordial binaries can evolve to be observable LISA sources at the present age of 14 Gyr. For one selected binary sample, its number as LISA sources can be calculated by multiplying the birthrate $R_{j}$ by the timespan spent in the LISA band. As a consequence, we can estimate the total number of detectable LISA DCOs in M31 by accumulating the contribution from all selected binary samples.

For GW signal detected over a mission duration (4 or 10 yr), the characteristic strain $h_{c,n}$ at the $n\mathrm{th}$ harmonic can be expressed as (Barack & Cutler, 2004),

Here $\dot{E}_{n}$ is the power of the $n\mathrm{th}$ harmonic radiation (Peters & Mathews, 1963),

where $\mathcal{M}_{c}$ is the chirp mass of a source and $g(n,e)$ (depending on the harmonic order $n$ and the eccentricity $e$) is the Fourier decomposition of GW signal and gives the relative power of gravitational radiation at the nth harmonic. If binaries have eccentric orbits, they would emit most of their energy at higher harmonics and be detected in a higher frequency band (Nelemans et al., 2001).

The change rate of GW frequency $\dot{f}_{n}$ can be written as

where $F(e)=\left[1+(73/24)e^{2}+(37/96)e^{4}\right]/\left(1-e^{2}\right)^{7/2}$ gives the integrated enhancement factor of gravitational radiation from an eccentric source compared with an equivalent circular source. The luminosity distance $d_{\mathrm{L}}$ of M31 is taken to be $780\mathrm{~{}kpc}$ (Ribas et al., 2005), which reduces the characteristic strain amplitudes of tight DCO systems in M31 by over 30 times compared with those in the Milky Way.

In this work, we use the python package LEGWORK developed by Wagg et al. (2022a) to calculate signal-to-noise ratio (S/N) and identify the DCO systems as LISA sources with the condition of $\mathrm{S/N>5}$ (e.g., Lamberts et al., 2018; Shao & Li, 2021). In the following, we briefly describe the method used in the LEGWORK.

The S/N calculation includes the response of detector to signal and the detector’s noise. For quasi-stationary sources we consider here, they evolve slowly during detector’s observation mission and their sky locations in the detector-based coordinate system are always changing. The averaged spectral power of GW signal is used to determine which sources can be detected and has the expression of

Here $\mathcal{R}(f)$ is sky and polarization averaged signal response function of the instrument. It has been numerically computed by Larson et al. (2000), and then fitted by Robson et al. (2019) as

The second term in right side of Equation (8) gives the spectral powers of plus and cross GW,

where $A(f)$ and $\Psi(f)$ represent the amplitude and phase of GW, respectively, $\iota$ describes the inclination of the orbital angular momentum of a binary with respect to the line of sight.

For sources with known positions, we should consider GW signal modulation to obtain more precise S/N. The orbital motion of LISA induces a power spreading effect, where the power spectrum of the GW signal exhibits sidebands with a bandwidth of $10^{-4}\,\mathrm{mHz}$ around the GW frequency $f$. As a result, the strain amplitude at the true frequency is diminished. This effect is more pronounced at higher frequencies due to the larger bandwidth. Consequently, WDWD binaries are particularly impacted due to their relatively higher detectable frequencies. We make the assumption that all sources detectable by LISA share the same sky location as M31, with the ecliptic coordinates ($\theta=33.35^{\circ},\phi=27.85^{\circ}$). Additionally, we consider their polarization $\psi$ and inclination $\iota$ to be uniformly distributed within the allowed range. In this study, we utilize the orbit-averaged amplitude modulation, originally defined by Cornish & Larson (2003) and subsequently revised by Wagg et al. (2022a), as

where $\mathcal{A}$ is the intrinsic amplitude of the source and $A$ is the amplitude recorded in the detector. The amplitude modulation is given as

where the orbit-averaged detector responses $\langle F_{+}^{2}\rangle$ and $\langle F_{\times}^{2}\rangle$ are the functions of sky location and polarization of sources in the ecliptic coordinates.

For the noise, we apply LISA’s sensitivity curve constructed by Robson et al. (2019). The effective power spectral density of the noise $S_{\mathrm{n}}(f)$ is expressed as

where $P_{n}(f)$ is the power spectral density of the LISA noise, including the single-link optical metrology noise $P_{\mathrm{OMS}}$ and the single test mass acceleration noise $P_{\mathrm{acc}}$,

where $L=2.5\mathrm{~{}Gm}$ is the arm-length of LISA and $f_{*}=c/(2\pi L)=19.09\mathrm{~{}mHz}$ is LISA transfer frequency. Actually, $P_{\mathrm{acc}}$ dominates $P_{n}(f)$ when $f$ is below $\sim\mathrm{5mHz}$ otherwise $P_{\mathrm{OMS}}$ dominates. For convenience, we ignore the term $f^{4}$ in $P_{\mathrm{acc}}$ and the term $f^{-4}$ in $P_{\mathrm{OMS}}$ when calculating $P_{n}(f)$. In addition, $S_{c}(f)$ is the effective noise source resulted from unresolved galactic binaries (Cornish & Robson, 2017; Babak et al., 2021) and can be alleviated by the increase of observation time. In our calculations, we apply the $S_{c}(f)$ calculated by Cornish & Robson (2017).

After the analyses of signal and noise, we can calculate the sky, inclination and polarization averaged S/N (Wagg et al., 2022a)

The equation above for the $\langle\rm S/N\rangle$ is suitable whether the binary is circular or eccentric and stationary or non-stationary.

It is noted that the amplitude modulation derived by Cornish & Larson (2003) is just suitable for quasi-circular systems. For WDWD systems detected by LISA, it is proper to add the amplitude modulation when calculating S/N. While for eccentric binaries with BH and/or NS components, we take the averaged value of $A_{\mathrm{mod}}$. After calculation, we find that the ratio of S/N with amplitude modulation to S/N without amplitude modulation is about $0.6$ for the same circular source.

After analysing our obtained data, we find that the most common evolutionary path for the formation of LISA DCO binaries involves the canonical channel proposed by the reviews of Bhattacharya & van den Heuvel (1991) and Tauris & van den Heuvel (2023). Evolved from a primordial binary, the primary star firstly evolves to fill its Roche lobe and stably transfers its envelope to the secondary star. After mass transfer, the primary star becomes a helium star and then evolves to be the first compact object. When the secondary star also evolves to fill its own Roche lobe, the binary enters the X-ray binary phase. Since the mass ratio between the secondary star and the compact star is relatively large, the mass transfer is dynamically unstable leading to the CE phase. During CE evolution, the secondary star fully loses its hydrogen envelope to become another helium star and finally evolves to be the second compact object. This CE phase leads to significant orbit decay of the binary so that the formed DCO is able to appear as a GW source. In some cases, the formation of tight DCO binaries with BH components does not involve this CE phase and the mass transfer is stable during the X-ray binary stage (see also e.g., Pavlovskii et al., 2017; Shao & Li, 2021). According to our calculations, the formation of LISA BHBH and BHNS binaries can be dominated by the channel of stable Roche lobe overflow if one takes $\alpha_{\rm{CE}}=0.3$ and $\alpha_{\rm{CE}}=1.0$. In the $\alpha_{\rm{CE}}=3.0$ case, all kinds of LISA DCO systems are predominantly formed through the channel of CE evolution.

Figure 2 presents estimated numbers of various types of DCOs in M31 that can be detected by LISA for a mission duration of 4 (top panel) or 10 yr (bottom panel), by assuming different models related to CE ejection efficiencies and supernova-explosion mechanisms. These numbers are also given in the Appendix. We can see that, for a 4 yr observation duration, at most $\sim 3$ DCO systems are likely to be identified by LISA. In this case, the systems with BH components dominate overall GW sources albeit in small numbers. For a 10 yr observation duration, the total number of LISA sources significantly increases to nearly a dozen. The numbers of the LISA binaries with BH components show a slight increase, while the numbers of the systems with WD components (i.e. WDWDs and NSWDs) increase by a factor of $\gtrsim 10$. For the rapid supernova mechanism (indicated by black dots in the diagram), no BHWD binaries can form within the LISA band. As this mechanism is only able to produce BHs with masses of $>5M_{\odot}$ and the masses of WD’s progenitors are less than $\sim 8M_{\odot}$, the mass ratios of the WD’s progenitors to the BHs are so small that mass transfer occurs via stable Roche lobe overflow rather than CE evolution, eventually leading to the formation of long-period BHWD binaries (see also Shao & Li, 2021). Overall, the numbers of different types of LISA sources are affected by adopted CE ejection efficiencies and supernova-explosion mechanisms. There is a common feature that stochastic supernova mechanism tends to produce more LISA sources as relatively small natal kicks are imparted to NSs and BHs (Mandel & Müller, 2020).

Figure 3 shows the $\rm S/N$ distribution of M31 LISA sources as a function of orbital periods and chirp masses for 4 (top panel) or 10 yr (bottom panel) observation duration. The black solid curve gives a boundary of $\rm S/N=5$, indicating that only the DCO systems located above this curve are considered detectable by LISA. The blue and red dashed lines correspond to the chirp masses of typical NSWD and WDWD systems, respectively. In the case with $T_{\rm obs}=4\rm\,yr$, the red dashed line is below the black solid curve, indicating a lack of detectability of WDWD binaries. In addition, the blue dashed line is very close to the black solid curve. This means that the maximal $\rm S/N$ of merging NSWD binaries can just reach LISA’s detection threshold. In the case with $T_{\rm obs}=10\rm\,yr$, we see that the regions where WDWD and NSWD systems can be detected are obviously enlarged. By contrast, the detection of the systems with BH components is not sensitive to the options of instrument’s observation lifetime since these systems have relatively large chirp masses.

Considering that a number of GW sources are expected to be observed by LISA for a 10 yr mission lifetime, we will analyse their parameter distributions in the following.

Figure 4 presents the probability density functions of M31 LISA WDWDs as a function of component masses and orbital periods. Since WDWD systems do not undergo supernova explosions during their formation, we discuss the calculated results with different CE ejection efficiencies. Our results show that almost all of these LISA sources consist of CO WDs and/or ONe WDs. Note that the Galactic WDWD sample identified via electromagnetic observations mostly include He WDs (Korol et al., 2022, and references therein). These binaries, characterized by small chirp masses, would not be detectable as M31 LISA sources. Under the assumption of different $\alpha_{\mathrm{CE}}$, the masses of the secondary WDs are distributed with a peak at $\sim 0.7-0.8M_{\odot}$ and the masses of the primary WDs have a relatively flat distribution varying in the range of $\sim 0.8-1.4M_{\odot}$. Importantly, these WDWD systems possess total masses exceeding the Chandrasekhar limit, indicating they are promising progenitors of type Ia supernovae.

For $\alpha_{\mathrm{CE}}=0.3$, there is a remarkable subgroup of systems consisting of a $\sim 1.0-1.2M_{\odot}$ ONe WD and a $\sim 1.2-1.4M_{\odot}$ CO WD. These systems are evolved from the primordial binaries with both component masses of $\sim 6-8M_{\odot}$ in $\sim 20$ d orbits. When $\alpha_{\mathrm{CE}}=1.0$ or 3.0, the primordial binaries with similar parameters always evolve into wide WDWD systems that do not contribute to the population of LISA sources.

The orbital periods of WDWD binaries are always less than 6 min (see also Figure 3) and distributed with a peak at 2 min. Note that the $\alpha_{\mathrm{CE}}=0.3$ scenario leads to a bulge at $\sim 4-6$ min in the curve of orbital period distribution, corresponding to the formation of very massive WDWD systems as mentioned above. The orbital period distributions of post-CE binaries are expected to be influenced by the CE ejection efficiencies. However, it should be noted that very close systems detectable by LISA have undergone significant orbital shrinkage due to GW radiation. Therefore, different assumptions of $\alpha_{\mathrm{CE}}$ do not significantly alter the shape of the orbital period distributions. For other types of LISA DCO systems, we focus on the role of supernova-explosion mechanisms in their parameter distributions and take $\alpha_{\mathrm{CE}}=1.0$.

In Figure 5, we show the probability density functions of M31 LISA NSWDs as a function of component masses and orbital periods under the assumption of $\alpha_{\mathrm{CE}}=1.0$. It can be clearly seen that different supernova mechanisms we adopted tend to produce similar parameter distributions for LISA NSWD systems. The NS masses are distributed with a sharp peak at $\sim 1.3M_{\odot}$, and the WD masses are mainly distributed in the range of $\sim 0.7-1.2M_{\odot}$ with a peak at $\sim 0.8M_{\odot}$. The orbital periods of these binaries are narrowly distributed within the range of $\lesssim 6$ min with a peak at $\sim 4$ min, which are greatly lower than the minimal orbital period of $\sim 0.1$ d for Galactic WD binaries with a pulsar companion (Manchester et al., 2005).

For these LISA sources, almost all NSs are formed via electron-capture rather than core-collapse supernovae. For core-collapse NSs, the resulting NS binaries with main-sequence companions of intermediate mass ($\sim 2-7M_{\odot}$) have orbital periods of $\lesssim 100$ d and most of them have orbital periods of less than a few days; while for electron-capture NSs, the corresponding systems possess orbital periods of $\lesssim 1000$ d and the majority are relatively wide binaries (see Figure 2 of Shao & Li, 2015). It is possible that the former binaries with relatively close orbits undergo mergers in subsequent CE phases, avoiding the formation of tight NSWD binaries.

Figures 6 and 7 present the probability density functions of LISA BHNS and BHBH systems in M31, respectively, as a function of component masses and orbital periods. Here the CE ejection efficiency is taken to be 1.0. Since the rapid supernova mechanism allows the existence of a $2-5$ mass gap between NSs and BHs, it predicts that the BH masses in both BHNS and BHBH systems are always larger than $5M_{\odot}$ and distributed with a peak at $\sim 7-8M_{\odot}$. On the contrary, the delayed and stochastic mechanisms are likely to form the systems with low-mass ($<5M_{\odot}$) BH components. Compared with the LISA sources with WD components, the BHNS (BHBH) systems have wider orbital period distributions in the range of $\lesssim 20$ min ($\lesssim 40$ min) with peaks at $\sim 10$ min ($\sim 20$ min).