The NANOGrav 11 yr Data Set: Limits on Supermassive Black Hole Binaries in Galaxies within 500 Mpc

In the local volume, converting redshift to distance is not trivial because the spectroscopically measured velocities are significantly affected by the motion of galaxies (e.g., due to infall to clusters). Therefore, it is necessary to correct for the distortions in the local velocity field or use direct distance measurements that are available for several local galaxies. For the latter, we used the following catalogs: (A) Cosmicflows-3 (Tully et al. 2016) extracted from the the Extragalactic Distance Database (Tully et al. 2009), and (B) the galaxy groups catalog (Crook et al. 2007).

Cosmicflows-3 is a compilation of distances of nearby galaxies, which were obtained with high-quality methods. These include methods that rely on standard candles, such as (1) the Cepheid period–luminosity relation, in which the intrinsic luminosity is determined from the pulsation period of the star; (2) the tip of the red giant branch, in which the brightest red giants in the Hertzsprung-Russel diagram are used as standard candles; and (3) Type Ia supernovae. Cosmicflows-3 also incorporates distances from additional methods, including (4) the Tully-Fisher relation, an empirical correlation between the rotation velocity of spiral galaxies and their intrinsic luminosity; (5) the fundamental plane of elliptical galaxies, a bivariate correlation between the luminosity, the effective radius, and the dispersion velocity of the galaxy; and (6) the surface brightness fluctuation method, which relies on the variance of the observed light distribution (e.g., distant galaxies have lower variance across pixels and appear smoother). Among the galaxies in the 2MRS sample, 8625 were covered in Cosmicflows-3.

For galaxies not included in Cosmicflows-3, we used the distance estimates from the catalog of galaxy groups. In this catalog, a "friends-of-friends" algorithm was employed to determine galaxy groups by identifying neighboring galaxies with similar distances. The galaxy distances were determined by applying a flow-field model, following Mould et al. (2000). More specifically, the heliocentric velocities were first converted to the Local Group frame. Subsequently, the galaxies in the vicinity of major local clusters (i.e., the Virgo Cluster, the Shapley Supercluster, and the Great Attractor) were assigned the velocity of the respective cluster (see Mould et al. 2000; Crook et al. 2007, for details). Finally, an additional correction was added to account for the gravitational pull toward the aforementioned major galaxy clusters. We used the high-density-contrast (HDC) catalog from Crook et al. (2007), in which a galaxy was considered a group member when the density contrast was 80 or more, and for groups with more than two members, we used the mean distance of the members. The group catalog provided the distances for another 4533 galaxies.

For the remaining galaxies (i.e., galaxies not included in any of the above two catalogs), we followed the prescription from Mould et al. (2000), described above, to correct for velocity distortions in the local volume. We adopted a nominal value of H₀ = 70 km s⁻¹ Mpc⁻¹ and accordingly adjusted the distances from Cosmicflows-3 and galaxy groups, in which the adopted values were 75 and 73 km s⁻¹ Mpc⁻¹, respectively. In our calculations, we did not include the distance uncertainties, primarily because they are negligible compared to the uncertainties of the mass estimates (see below).

For the SMBH masses, we used direct measurements when available and global scaling relationships otherwise. Dynamical mass measurements are observationally expensive and are limited only to a small number of galaxies. We used the catalog compiled by McConnell & Ma (2013), which includes 72 galaxies with high-quality dynamical mass measurements, from observations of stellar, gaseous, or maser kinematics. We added to this catalog recent measurements for the galaxies NGC 4526, M60-UCD1, NGC 1271, NGC 1277, and NGC 1600 (Davis et al. 2013; Seth et al. 2014; Walsh et al. 2015, 2016; Thomas et al. 2016). We also included 29 galaxies from the Active Galactic Nucleus (AGN) Black Hole Mass Database, ⁴⁰ the masses of which were determined with spectroscopic reverberation mapping (Bentz & Katz 2015). However, the above dynamical tracers probe spatial scales that cannot distinguish whether the mass corresponds to one or two SMBHs and thus reflect the total enclosed mass.

For the remaining galaxies, we estimated the SMBH masses from correlations with observable quantities, like the M_BH–σ or the M_BH–M_bulge correlations. First, we extracted measurements of the velocity dispersion σ for 3300 galaxies from the Lyon-Meudon Extragalactic Database (LEDA), using the Extragalactic Distance Database to cross-correlate our sample with LEDA. We used the M_BH–σ relation from McConnell & Ma (2013),

$$\begin{eqnarray}&&{\mathrm{log}}_{10}({M}_{\mathrm{BH}})={a}_{1}+{b}_{1}{\mathrm{log}}_{10}(\sigma /200\,\mathrm{km}\,{{\rm{s}}}^{-1}),\end{eqnarray} \tag{ 3 }$$

with (a₁, b₁) = (8.39, 5.20) for early-type (elliptical and S0) galaxies and (a₁, b₁) = (8.07, 5.06) for late-type (spiral) galaxies. Of the 3300 galaxies, we excluded 810 galaxies because their velocities were beyond the range from which the correlation was established (80–400 km s⁻¹) or their dispersion measurements had significant uncertainty (>10%). Additionally, the morphological types of 284 galaxies were not included in the catalog, or they were classified as irregular. For this subset, we did not use the M_BH–σ relation.

If a measurement of the dispersion velocity was not available (or we excluded it for reasons mentioned above), we used the M_BH–M_bulge relation to estimate the SMBH mass. We first used the K-band magnitude to estimate the stellar mass, which in turn was used to calculate the bulge mass. In particular, we determined the absolute magnitude M_K from the observed m_K as

$$\begin{eqnarray}&&{M}_{K}={m}_{K}-5\times {\mathrm{log}}_{10}(D)-25-0.11\times {A}_{v},\end{eqnarray} \tag{ 4 }$$

where D is the distance in Mpc and A_v is the extragalactic extinction, which we calculated from the reddening relation R_V = A_V /E_B−V = 3.1 (Fitzpatrick 1999), with the color term E_B−V provided by the 2MRS catalog. Next, we obtained the total stellar mass of the galaxy M_* using the correlation (Cappellari 2013)

$$\begin{eqnarray}&&{\mathrm{log}}_{10}({M}_{* })=10.58-0.44\times ({M}_{K}+23).\end{eqnarray} \tag{ 5 }$$

At this step, we excluded 21 quasars, since the K-band luminosity does not correspond to the stellar luminosity but is dominated by dust emission (Barvainis 1987).

Next, we estimated the bulge mass as M_bulge = f_bulge M_*, with f_bulge the fraction of the stars that reside in the bulge. For early-type galaxies, we assigned their total stellar masses to their bulges (f_bulge = 1), whereas for late-type galaxies, we used the bulge-to-total mass ratio, determined from image decomposition (Weinzirl et al. 2009), which depends on the galaxy's subclass. The bulge fraction also depends on the existence of a galactic bar (e.g., Weinzirl et al. 2009, Figure 14), but since the 2MRS catalog typically does not specify whether a galaxy has a bar, we used the average value of barred and unbarred galaxies for each subclass. In Table 1, we summarize the adopted values of bulge-to-total mass ratio (or equivalently the bulge fraction f_bulge), as a function of the galaxy type. ⁴¹ Additionally, the galaxy types are unknown for ∼40% of the galaxies in the catalog. For those, we estimated their bulge fractions assuming two different scenarios for the galaxy type (see Section 2.3).

We estimated the SMBH mass using the correlation from McConnell & Ma (2013),

$$\begin{eqnarray}&&{\mathrm{log}}_{10}({M}_{\mathrm{BH}})={a}_{2}+{b}_{2}{\mathrm{log}}_{10}({M}_{\mathrm{bulge}}/{10}^{11}\,{M}_{\odot }),\end{eqnarray} \tag{ 6 }$$

with (a₂, b₂) = (8.46,1.05). The M_BH–σ and M_BH–M_bulge correlations have intrinsic scatters of = 0.38 and = 0.44, respectively. We used these scatters as a proxy for the uncertainties in the mass estimates, but we did not include additional sources of uncertainty (e.g., due to galaxy classification, bulge fraction, etc).

We calculated the GW strain amplitude from Equation (1), using the estimates of the SMBH masses (regarded as the total masses of the binaries) and the distances obtained above. We considered a galaxy to be within the volume that NANOGrav can probe if the GW strain of the tentative binary in the most optimistic scenario in terms of mass ratio and binary orbit, i.e., equal-mass binary in a circular orbit emitting GWs at NANOGrav's most sensitive frequency, ${h}_{\mathrm{GW}}\left(q=1,\ {f}_{\mathrm{GW}}=8\,\mathrm{nHz}\right)$, was equal to or higher than the upper limit from the 11 yr data set. However, as shown in Figure 1 (top panel), the sensitivity of NANOGrav depends on the spatial distribution of pulsars on the sky and is highly anisotropic. We thus compared the GW strain to the upper limit in the direction of the source, ${h}_{95}\left({\rm{R}}.{\rm{A}}.\ \mathrm{Decl}.\right)$, using a sky map with resolution of 768 pixels (healpix with nside = 8) from Aggarwal et al. (2019). Therefore, a binary could be detectable by NANOGrav if

$$\begin{eqnarray}&&{h}_{\mathrm{GW}}\left({\rm{q}}=1,\ {{\rm{f}}}_{\mathrm{GW}}=8\,\mathrm{nHz}\right)\geqslant {h}_{95}\left(\mathrm{Ra},\ \mathrm{Dec}\right).\end{eqnarray} \tag{ 7 }$$

In Table 2, we present the entire 2MRS galaxy catalog, with the estimated SMBH masses and galaxy distances (including the method with which each quantity was calculated) and the relevant GW upper limits at f_GW = 8 nHz.

Zoom In Zoom Out Reset image size

Among the galaxies with known types, there are 194 galaxies with SMBH mass of 10⁹–10¹⁰ M_⊙ in the volume NANOGrav can probe (up to 500 Mpc). ⁴² The vast majority of those are early-type galaxies, with only four classified as S0/a. Next, we examined the unclassified galaxies and estimated their masses examining two possible scenarios for the galaxy type: (1) We assumed that they are late-type Sa galaxies with f_bulge = 0.31, which corresponds to the most massive SMBHs for the case of spiral galaxies. Under this assumption, none of them are within NANOGrav's reach. (2) We assumed that they are early-type galaxies (f_bulge = 1), in which case another 84 galaxies enter NANOGrav's sensitivity volume. We visually inspected those 84 galaxies, using data from the online database SIMBAD. ⁴³ We found that 22 of those are indeed early-type galaxies, classified as brightest cluster galaxies (BCGs). Six are classified as quasars/AGNs, which we excluded because we cannot reliably estimate the stellar mass from the K-band magnitude (Section 2.2). The remaining 56 galaxies are impossible to classify based on the available images. We kept these galaxies as a separate subsample, the analysis of which is presented in the Appendix.

Figure 1 (bottom panel) shows the 194 galaxies, along with the 22 BCGs among the subset with initially unknown type, that are within NANOGrav's sensitivity volume, color-coded according to their distance, with the marker size denoting the estimated SMBH mass. As expected, there are significantly more galaxies in the most sensitive half of the sky, where most of the systematically monitored pulsars are located. In the least sensitive half of the sky, binaries can be detectable only if they reside in very massive or in relatively nearby galaxies.

Additionally, we calculated the number of galaxies that NANOGrav can probe considering the uncertainty (equal to one standard deviation) in the mass estimates, ${\sigma }_{{M}_{\mathrm{BH}}}$. Assuming the total mass of the binary to be equal to ${M}_{\mathrm{BH}}-{\sigma }_{{M}_{\mathrm{BH}}}$ and ${M}_{\mathrm{BH}}+{\sigma }_{{M}_{\mathrm{BH}}}$, we found 20 and 2313 galaxies with known morphological classification within the sensitivity volume. When we also included the galaxies of unknown type, there are still 20 galaxies in the volume for the lower-mass bound, whereas for the high mass bound the number becomes 2338 and 4545 for the two scenarios described above (i.e., late- vs. early-type galaxies). However, in this case it is impractical to visually examine over 2000 galaxies to verify whether they are indeed elliptical galaxies. The impact of the mass uncertainty is obvious from the large discrepancy between the above numbers. We further explore this issue in Section 5.3.

Identical binaries (in terms of total mass, mass ratio, and distance) may have different probabilities of detection depending on their coordinates (Figure 1). We quantified the probabilities of detection by calculating the S/N for putative binaries in the 216 galaxies that fall within the sensitivity volume of NANOGrav. (The Appendix presents the S/N of the 56 galaxies of unknown type that are detectable by NANOGrav, only under the assumption that they are early-type galaxies.)

We calculated the total S/N of a binary, adding the S/N_j in each individual pulsar, following Rosado et al. (2015). For this purpose, we assumed equal-mass binaries in circular orbits and a fixed GW frequency of 8 nHz. Since the S/N depends on the inclination i, the initial phase Φ₀, and the GW polarization angle ψ, we randomly draw these parameters from uniform distributions, in the range [−1, 1] for $\cos i$, [0, π] for ψ, and [0, 2π] for Φ₀. We estimated the average total S/N from 1000 realizations for each putative binary. In Table 3, we ranked the binaries from the highest to lowest based on their average S/N.

We used the NANOGrav 11 yr data set to derive constraints on hypothetical SMBHBs in the 216 galaxies identified in Section 2. Even though this data set allows us to probe frequencies up to ∼317.8 nHz, we limited our analysis to 10 nHz for computational efficiency. We further explain this choice toward the end of the section. As mentioned above, in order to avoid unnecessary computations, we calculated the GW strain upper limits in 110 pixels (out of the 192 in which we divided the sky).

In Figure 3, we show the 95% upper limits on the mass ratio (which corresponds to the solid line in the bottom panel of Figure 2) as a function of GW frequency for eight galaxies in the sample. The four galaxies on the top and bottom panels are distributed close to the most and least sensitive sky location, respectively. We randomly selected these galaxies, in order to illustrate a wide range of total masses and distances of the hypothetical binaries; these quantities are shown in the brackets next to the abbreviated galaxy names in the legend. For a significant fraction of the galaxies we can place stringent constraints on the mass ratio, for several GW frequencies. Therefore, if there is a circular binary in the center of these galaxies, the mass of the secondary SMBH must be only a few percent that of the primary. We tabulate the constraints on mass ratio in Table 3.

Zoom In Zoom Out Reset image size

It is obvious from Figures 2 and 3 that the mass ratio upper limit is a declining function of the frequency at low frequencies, whereas at high frequencies (beyond a critical frequency) the function becomes flat and the constraints on the mass ratio are independent of frequency. The flattening of the mass ratio upper limit at high frequencies can be derived analytically. Since in our analysis we did not include uncertainty in the distance estimates (i.e., for each galaxy the distance is a fixed value), from Equation (1) we can directly relate the 95% upper limits of the GW strain, h₉₅, with the chirp mass upper limits ${{ \mathcal M }}_{95}$. Solving for chirp mass limits,

$$\begin{eqnarray}&&{{ \mathcal M }}_{95}=\frac{{h}_{95}^{3/5}{\left({c}^{4}D\right)}^{3/5}}{{2}^{3/5}{\left(\pi {f}_{{GW}}\right)}^{2/5}}.\end{eqnarray} \tag{ 10 }$$

For frequencies above 8 nHz (f_GW > 8 nHz), the GW strain upper limit has a power-law dependence on frequency ${h}_{95}({f}_{\mathrm{GW}})\propto {f}_{\mathrm{GW}}^{2/3}$, which comes from the fact that at those frequencies the GW sensitivity is set by the level of white noise in the pulsar residuals (see, e.g., Figures 3 and 4 in Aggarwal et al. 2019). Then, from Equation (10), it follows that

$$\begin{eqnarray}&&{{ \mathcal M }}_{95}\propto \frac{{\left({c}^{4}D\right)}^{3/5}}{{2}^{3/5}{\left(\pi \right)}^{2/5}}.\end{eqnarray} \tag{ 11 }$$

Therefore, the chirp mass upper limit is constant, independent of the frequency for f > 8 nHz. Then, from Equation (9) the mass ratio upper limit is also constant and independent of frequencies for f > 8 nHz. This is numerically demonstrated in Figure 5, where we extended the analysis to higher frequencies. As a result, our conclusions can be extended to higher frequencies without additional computations, which explains our decision to restrict the analysis to 10 nHz.

Since it is not practical to show the mass ratio constraints as a function of frequency for the entire sample, in Figure 4 we present the 95% upper limits for a fixed frequency of 8 nHz for the 216 galaxies in our sample, with the top/bottom panel illustrating galaxies in the most/least sensitive half of the sky. There is an obvious trend from yellow to dark-blue points, which indicates that increasingly unequal-mass binaries are required, when the distances decrease for fixed mass (vertically down) and when the total mass increases for fixed distance (horizontally to the right), as expected. The same general trend is repeated in the least sensitive half of the sky, but the overall curve is shifted down, which means that we are limited to more massive and more nearby galaxies. We emphasize that we selected this particular frequency (f_GW = 8 nHz) for two reasons: (1) this is the most sensitive frequency of NANOGrav—for this reason, the galaxy sample was initially selected assuming tentative binaries emitting GWs at 8 nHz, as shown in Figure 1; and (2) this frequency corresponds to the flat part of the mass ratio upper limit curve, and thus our conclusions can be extended to higher frequencies.

Zoom In Zoom Out Reset image size

For the five top-ranked galaxies based on their S/N, we run the upper limits analysis sampling the likelihood of Equation (8) directly in the mass ratio. This allows us to incorporate the uncertainty in the total binary mass as a prior. In Figure 5, we show the 95% mass ratio upper limits as a function of frequency. For these five galaxies, we run the analysis for the full range of NANOGrav frequencies up to the cutoff frequency, which depends on the binary total mass. ⁴⁴ This figure also confirms that the constraints are roughly constant at high frequencies, except for specific frequencies, for which the sensitivity of NANOGrav is limited (e.g., at ∼30 nHz, which corresponds to a period of 1 yr and is caused by fitting the pulsars' positions).

Zoom In Zoom Out Reset image size

In Figure 5, we also show with red dashed lines the limits derived from the GW strain using the mean SMBH as the total mass of the binary in Section 4.2 for comparison. As expected, including the uncertainty leads to weaker upper limits, which is more pronounced for galaxies with highly uncertain mass measurements, like NGC 4889 and NGC 5062. ⁴⁵ The other three galaxies have dynamically measured SMBH masses and thus relatively small uncertainties. The upper limits here are slightly higher but comparable to the ones obtained above, especially at higher frequencies (flat part of the upper limit curve). Additionally, the constraints we derive for putative binaries in these three galaxies are significant with only very unequal-mass binaries allowed, i.e., the secondary SMBH should be at least 10 times less massive than the primary for the majority of frequencies we examined. We further explore the impact of the uncertainty in the total mass in Section 5.3.

We explored whether the galaxies in our sample could host a binary delivered by a major merger. Major mergers are considered the mergers between galaxies with similar mass (i.e., minimum mass ratio of 1/4) and play an important role in galaxy formation and evolution (Volonteri et al. 2003; Kelley et al. 2017). Since the mass of the central SMBH is tightly correlated with the galaxy mass, we assume that a major merger would result in a binary with a mass ratio q > 1/4. This means that major galaxy mergers are also significant for PTAs, since they produce promising binaries for GW detection.

Based on the constraints on the mass ratio derived in Section 4.2, we excluded binaries delivered by major mergers for galaxies where q₉₅ < 1/4. In Figure 6, we show the range of frequencies for which we were able to exclude the existence of binaries from major mergers for the 50 top-ranked galaxies based on S/N. Each bar corresponds to a distinct galaxy, the abbreviated name of which is shown on the left. The color-coding reflects the distance of the binary, and the size of the marker on the right reflects the total mass of the binary (similar to Figure 1). The gaps in the bars indicate frequencies where NANOGrav's sensitivity is lower and binaries with q > 1/4 cannot be excluded. It is obvious from the figure that the most informative constraints are derived for relatively nearby galaxies (within 100 Mpc).

Zoom In Zoom Out Reset image size

Motivated by the above findings, we estimated the number density of nearby galaxies, for which we can exclude binaries delivered by major mergers (q₉₅ < 1/4) at NANOGrav's most sensitive frequency (f_GW = 8 nHz), as a function of the binary total mass. For this, we first estimated the maximum volume out to which NANOGrav can probe binaries from major mergers; for binaries of fixed chirp mass, i.e., fixed mass ratio (q = 1/4) and total mass (see Equation (1)), the GW upper limit can be converted to a maximum distance. Due to heterogeneous sensitivity, the maximum distance is different for each pixel, and thus the volume NANOGrav can probe is a shell of irregular shape. For each pixel, we calculated the maximum distance and the respective spherical volume that NANOGrav would cover if the sensitivity were homogeneous. Since the pixels have equal area, we approximated the NANOGrav volume with a sphere of an effective volume equal to the average spherical volumes that correspond to each pixel. Then, we calculated the number of galaxies within this volume. Dividing the number of galaxies by the effective volume of the shell, we get the NANOGrav limit on the density of galaxies that could host binaries delivered by major mergers at 8 nHz.

Next, we examined several discrete values for the total binary mass in the range of 10⁹–10¹⁰ M_⊙, where all the galaxies identified in Section 2.3 lie. As the binary total mass increases, the distance NANOGrav can probe becomes larger, which allows more galaxies to enter the sensitivity volume. However, a countereffect is that massive galaxies, with higher SMBH masses (e.g., with ${\mathrm{log}}_{10}({M}_{\mathrm{BH}}/{M}_{\odot })\gt 9.5$), are extremely rare. In Figure 7, we show with a black solid line the NANOGrav 95% upper limit on the density of binaries produced by major mergers. The error bars denote the range of densities, if we consider the 1σ uncertainty in the SMBH mass. For comparison, we also show the overall galaxy density above a certain fixed mass for spherical volumes with radii of 100, 300, and 500 Mpc. We see that the density of galaxies is decreasing when we examine a larger volume (from a radius of 100 Mpc to a radius of 500 Mpc). This may mean that the universe has higher density locally, potentially because we consider a small number of galaxies concentrated in local clusters (Virgo Cluster at 16 Mpc, Coma Cluster at 120 Mpc). Finally, the decreased density is unlikely to be caused by the incompleteness of the 2MRS catalog, since for the galaxies of interest (with M_BH > 10⁹ M_⊙) the catalog is almost complete (∼97.6%) up to 500 Mpc. We note, however, that beyond 300 Mpc the catalog becomes incomplete for lower-mass galaxies.

Zoom In Zoom Out Reset image size

We emphasize that this is the first constraint on major galaxy mergers derived directly by GW data. However, it is worth pointing out that there are several caveats directly connecting a major merger to a binary with q > 1/4. The main caveats are summarized below: (1) The scaling relations between the SMBH and the host galaxy have significant scatter, and thus the binary mass ratio does not necessarily reflect the mass ratio of the initial galaxies. (2) The mass of each SMBH may evolve significantly in the post-merger galaxy, meaning that the binary mass ratio at the final stages may not be directly linked to the initial mass ratio of the two SMBHs. For instance, if the secondary SMBH accretes at a higher rate, as seen in binaries embedded in circumbinary disks, the mass ratio tends to equalize (Farris et al. 2014; Siwek et al. 2020). (3) In the mass ratio upper limits, we did not include the uncertainty of the SMBH mass estimate. Therefore, the q < 1/4 cutoff is used merely as proxy for a binary produced by a major galaxy merger.

SM16 used PPTA and EPTA upper limits and performed a similar analysis on a subset of galaxies. In particular, they focused on galaxies with dynamical mass measurements (77 galaxies) and galaxies from the MASSIVE survey (116 galaxies). The total sample consisted of 185 galaxies, since eight of the MASSIVE galaxies had dynamical mass measurements. In Figure 8, we show with blue diamonds the distribution of these galaxies as a function of distance and SMBH mass. We compare them with the sample of galaxies in this paper, illustrated with gray triangles, and the galaxies in Mingarelli et al. (2017), shown with black squares (see below). We also summarize the galaxy selection and properties, along with key differences between the samples in Table 4. Of the 185 galaxies, 33 were in the sensitivity volume of PPTA and 19 were in the EPTA volume. In Section 2.3, we found that 24 of these galaxies are within the NANOGrav volume at 8 nHz. Since each array has distinct and heterogeneous sensitivity, which depends on the distribution of pulsars each collaboration systematically monitors, it is not surprising that each array can probe a different number of binaries. For instance, PPTA is more sensitive in the southern hemisphere, where most of the PPTA pulsars lie (e.g., Figure 2 in SM16).

Zoom In Zoom Out Reset image size

In Figure 9, we show the mass ratio constraints for four of the top galaxies in SM16 and compare with the 95% upper limits we obtained in Section 4.2. We see that the NANOGrav upper limits are comparable or slightly better (e.g., NCG 4889, NGC 4649) in the flat part of the upper limit curve. On the other hand, at the lowest frequencies, our limits are less constraining. The differences between the results of the two studies are mainly due to a combination of the following reasons: (1) There are slight differences in the assumed properties of the putative binaries in the two studies. The total masses of the binaries in these four galaxies are identical, but the distances in our study are slightly larger. All else being equal, this would naturally lead to slightly weaker upper limits for the binaries considered here. (2) At the most sensitive frequency of 8 nHz, the sky-averaged 95% upper limit of NANOGrav is approximately h₀ < 7.3 × 10⁻¹⁵, which is 1.4 times lower than the EPTA limit (h₀ < 10⁻¹⁴) and 2 times lower than the PPTA limit (h₀ < 1.7 × 10⁻¹⁴) used in SM16. The improvement in the GW limits should lead to more stringent constraints on the mass ratio and partly explains why NANOGrav performs better at higher frequencies. (3) PPTA and EPTA have increased sensitivity at low frequencies, because some of the pulsars have long baselines of timing observations (longer than 17 yr). Therefore, the limits derived with the NANOGrav 11 yr data set are less constraining at these low frequencies. (4) The NANOGrav upper limits vary significantly from the most sensitive to the least sensitive sky location (from h₀ < 2.0 × 10⁻¹⁵ to h₀ < 1.34 × 10⁻¹⁴). This means that the relative position of the galaxy to the pulsars in the array is important, and thus discrepancies between the derived limits from the three arrays are expected. An example of this is NGC 1600, which is located at NANOGrav's least sensitive part of the sky, where EPTA has at least four pulsars in its vicinity.

Zoom In Zoom Out Reset image size

Mingarelli et al. (2017, hereafter M17) modeled the local GW sky, using galaxy merger rates from Illustris and the local massive galaxies as an observational anchor. For this, they extracted a sample of 5110 massive early-type galaxies from 2MRS, ⁴⁶ and they calculated the probability of each galaxy hosting a GW-emitting SMBHB system. These galaxies are illustrated with black squares in Figure 8.

Even though the starting point of both studies is the 2MRS catalog, the goals of the two studies are complementary, which explains why different selection criteria and cuts were employed. For instance, M17 assessed the likelihood of PTAs detecting GWs from individual sources, through multiple realizations of the local universe, and thus the completeness of the galaxy sample is important. On the other hand, we present a purely observational approach, in which we used the current GW upper limits to place constraints on circular binaries. For our purposes, compiling the most comprehensive catalog that includes all the galaxies that NANOGrav can probe is crucial.

Examining the galaxies in the M17 catalog, we found that only 106 lie in the NANOGrav sensitivity volume, a factor of two lower compared to the 216 from the catalog that we use here. This discrepancy can be understood in light of the selection criteria and calculation of the galaxy properties (SMBH masses and distances). Below we summarize the key differences between the catalogs (but see also Table 4):

• 1.  
  We examined the entire 2MRS catalog, whereas M17 imposed cuts on the absolute K-band magnitude (M_K < −25) and the galaxy distance (D < 225 Mpc). The magnitude cut has no effect on the sample of interest (i.e., galaxies with binaries detectable by NANOGrav), as it mostly excludes low-mass galaxies. However, the distance cut is significant, because it excludes massive galaxies at relatively large distances, but within NANOGrav's capabilities (see, e.g., Figure 8).
• 2.  
  We used a multilayered approach for the estimates of SMBH masses. The difference with M17 lies in the inclusion of mass estimates from reverberation mapping for AGNs and the estimates from the M_BH–σ relationship. Both methods provide better estimates compared to the M_BH–M_bulge relationship, because the former provides dynamical mass measurements, whereas the latter relies on fewer assumptions (see also Section 5.3). ⁴⁷
• 3.  
  For the distance estimates, we started from the Cosmicflows-3 catalog, since it compiles the highest-quality distance measurements available in the literature. For galaxies not included in Cosmicflows-3, we used the catalog of galaxy groups (Crook et al. 2007) or corrected for the local velocity field (Mould et al. 2000), as in M17. For the latter subset, the distances in the two catalogs are in excellent agreement, with the exception of galaxies in the vicinity of the Great Attractor, to which we assigned the velocity of the cluster. For the subset of galaxies in Cosmicflows-3, the distances in M17 are overestimated.
• 4.  
  We examined the entire 2MRS catalog, whereas M17 selected early-type (elliptical and S0) galaxies. Additionally, the classification in M17 was based on the extragalactic database Hyperleda, whereas we used the classification from the 2MRS catalog itself. We note that a comparison of galaxy types in the two catalogs indicates that not all 5110 galaxies in M17 are classified as early type in the 2MRS classification; a fraction of them (∼16%) are classified as late-type galaxies, whereas ∼15% are unclassified. With visual inspection of images in SIMBAD, we confirmed that at least some of the galaxies in the M17 sample are indeed late-type (spiral) galaxies. Although it is impossible to determine whether this affects their overall conclusions, if spiral galaxies are misclassified as elliptical, the general trend is that the assumed SMBH mass and in turn the GW strain are overestimated.

This work presented a study on how the NANOGrav data can be used to place constraints on SMBHBs in massive nearby galaxies. For most galaxies, we derived the chirp mass upper limits from the GW strain and converted them to mass ratio upper limits using the mean SMBH mass (measured or estimated from scaling relations) as the total mass of the binary. As explained above, the main reason for this approach is computational efficiency. However, one of the caveats is that it ignores the uncertainty in the total binary mass. As we demonstrated in Section 4.3, incorporating the SMBH mass uncertainty deteriorates the mass ratio upper limits, with larger uncertainties having a stronger impact on these constraints.

Here we further explore the effect of the mass uncertainty on the mass ratio upper limits. For this, we derived the upper limits from the GW strain (as in Section 4.2), but instead of using the mean SMBH mass in Equation (9), we sampled the total mass from the entire mass distribution, which includes the uncertainty. Next, we compared with the upper limits we obtained sampling directly in the mass ratio in Section 4.3, where the mass uncertainty was included as a prior in the analysis. We find that when the uncertainty is small (e.g., NGC 4649, NGC 1600, M87), the upper limits from the different methods are in excellent agreement. However, if the mass uncertainty is large (e.g., NGC 5062, NGC 4889), the mass ratio upper limits tend to be higher when we sample directly in mass ratio, which means that the two methods are not always equivalent. We carefully examined the cause for this discrepancy for the galaxy NGC 4889 at a fixed frequency of 8 nHz, which we illustrate in Figure 10. We see that the posterior distribution of the total mass is skewed toward smaller values compared to the prior, which is shown with a dotted histogram. This results in a tail in the mass ratio distribution, which leads to significantly higher 95% mass ratio upper limits.

Zoom In Zoom Out Reset image size

Currently, this is a significant limitation, since sampling the mass ratio for each galaxy is computationally expensive. This is why we chose to perform this type of analysis only for five tentative binaries. We defer a similar analysis for the entire sample to a future study. Additionally, for the vast majority of galaxies, the SMBH mass (attributed to the total mass of the binary) was estimated from global scaling relationships, the intrinsic scatter of which inevitably results in large uncertainties. The uncertainty is small only for galaxies with precise dynamical mass measurements, which are observationally demanding and limited to a small number of galaxies. Such observations are unlikely to become readily available for a large number of galaxies in the near future. However, in this study we also included a small subset of AGNs with SMBH masses measured from reverberation mapping. This is a promising direction, since planned multiepoch spectroscopic surveys, like SDSS-V (Kollmeier et al. 2017), will soon provide SMBH mass estimates for a large number of AGNs through reverberation mapping. This is a particularly exciting development for future studies, as AGNs are promising candidates for hosting SMBHBs.

The mass ratio upper limits depend on the total mass of the binary (Equation (9)), which was assumed to be equal to the mass of the SMBH and for the vast majority of galaxies was calculated from scaling relationships. In Sections 2.2 and 2.3, we described the assumptions required to determine the SMBH mass, such as the bulge fraction, the galaxy classification, the existence of a bar, etc. We also demonstrated the challenges in determining the SMBH mass, when the galaxy type is unknown. Another effect is related to the choice of the SMBH mass correlations; throughout the analysis, we have adopted the scalings (M_BH–σ and M_BH–M_bulge) from McConnell & Ma (2013). These correlations may provide higher SMBH mass estimates (due to measurements of SMBHs at the high-mass end) compared to other correlations (Graham 2016; Shankar et al. 2016; Sahu et al. 2019). Adopting a shallower correlation would result in lower-mass binaries and weaker mass ratio limits. Even though all of the above affect the determination of the mean SMBH mass (and in turn the mass ratio upper limits), they are likely negligible compared to the effect of the uncertainty due to the intrinsic scatter in the M_BH–σ and M_BH–M_bulge correlations, as shown above.

Another potential source of uncertainty in estimating the mass ratio limits is related to the distance determination. In this work, we did not include the uncertainty in the distance estimates, mainly because this uncertainty is typically smaller than that of the SMBH mass (especially when obtained from scaling relations). Additionally, the dependence of the GW strain on the distance is weaker compared to the dependence on the chirp mass. Overall, even if the impact of the distance uncertainty is likely negligible, it should be incorporated in future work. Future analysis should also explore the uncertainty in the value of the Hubble constant, which is necessary for the distance calculations.

Last but not least, in our analysis we assumed circular binaries, and thus our constraints cannot be generalized to noncircular binaries. In reality, binary orbits may deviate from circular owing to high eccentricity at formation, or because the eccentricity is excited as the binary interacts with the surrounding gaseous environment (Armitage & Natarajan 2005). (We note, however, that GW emission tends to circularize the binary orbits.) In eccentric binaries the GW strain is split among multiple harmonics, resulting in weaker GW signal. In practice, this means that if the galaxies we analyzed hosted eccentric binaries, the limits we would derive would be less constraining than the ones we presented here for circular binaries. In a future study we will derive the upper limits examining eccentric binaries, especially since noncircular binaries may be relatively common.