Recurrent Neutrino Emission from Supermassive Black Hole Mergers

With the first detections of astrophysical high-energy neutrinos (Aartsen et al. 2013, 2014) and gravitational waves (GWs; Abbott et al. 2016), the era of multimessenger astronomy has begun, enabling a new way of exploring the universe. We have learned from GW observations that there exists a significant population of stellar-mass binary black holes with properties we only begin to understand (Abbott et al. 2019a, 2019b). Neutron star mergers can now be investigated using a combination of GWs and the broad photon spectrum from radio up to GeV gamma-rays, in the future possibly adding TeV gamma-rays and neutrinos (Bartos et al. 2013; Abbott et al. 2017; Albert et al. 2017, 2019; Murase & Bartos 2019). We have learned from neutrino observations that the high-energy events are of extragalactic origin given the lack of clustering of events in the Galactic plane (Aartsen et al. 2016, 2019). The detection of a possible correlation of a gamma-ray flare from the blazar TXS 0506+056 with a high-energy neutrino in 2017 September is a first hint that active galaxies may produce a significant fraction of the observed flux of cosmic high-energy neutrinos (Aartsen et al. 2018a; Kun et al. 2020). A dedicated search for further neutrino flares from TXS 0506+056 in the past revealed a 3.5σ evidence for a long-duration (${110}_{-24}^{+35}$ days) flare of TeV neutrinos around 2.8 yr prior to the 2017 measurement (Aartsen et al. 2018b).

What makes the two TXS flares so difficult to explain is that (1) the duration of the 2014/2015 flare is long and the detected high-energy neutrinos have energies around ∼10–100 TeV, while the 2017 flare only consists of one neutrino with much higher energy (close to 300 TeV); (2) while the 2017 flare was detected due to the coincident detection of a GeV gamma-ray flare, the 2014/2015 flare is lacking any sign of correlation between gamma-ray and neutrinos. This makes the modeling of the multimessenger data challenging (Reimer et al. 2019), in particular when considering a one-zone model, although it is not impossible to explain the results (Halzen et al. 2019). Several ideas of how to produce these multimessenger signatures have been presented (see, e.g., Cao et al. 2020; Gao et al. 2019; Kun et al. 2019, 2020; Oikonomou et al. 2019; Rodrigues et al. 2019; Xue et al. 2019).

Transient blazar flares may be produced in the wake of galaxy mergers when two supermassive black holes (SMBHs) inspiral toward each other. The black holes on close orbits reorient their spin in this inspiral phase. Spin precession due to relativistic effects in turn can result in precessing relativistic outflows that periodically change the direction of high-energy radiation (Gergely & Biermann 2009). This scenario leads to the prediction of a population of blazars that are currently in such a state. The first potential hints of such signatures have been identified by Kun et al. (2017). Recent observations of periodicity at radio wavelengths also point toward a precessing jet scenario (Britzen et al. 2019).

In this Letter, we examine the observational consequences of transient blazar flares from jet precession in supermassive binary black holes (SMBBHs) close to merger. In particular, we investigate the possible time structure of high-energy neutrino emission from blazars due to jet precession, focusing on TXS 0506+056. Considering this scenario, we predict the time of the next neutrino flare from TXS 0506+056 and the expected time of the corresponding SMBBH merger. We discuss the possibility of detecting this merger through GWs using the Laser Interferometer Space Antenna (LISA), scheduled to launch in 2034.

Many accreting SMBHs are observed to drive relativistic jets. The rotation angle of the accretion disks around these black holes and consequently the direction of the jets are thought to be aligned with the spin axis of the black holes (Rees 1978). When two SMBHs inspiral toward each other, their spin precession will also change the orientation of the jet periodically.

We model the inspiral dynamics of a SMBH binary close to merger using the post-Newtonian (PN) approximation or the General Theory of Relativity up to 2.5 PN order (for a review see Blanchet 2014). The magnitudes of the spins of the two black holes (i = 1, 2) are defined as

$$\begin{eqnarray}&&{S}_{i}\approx \displaystyle \frac{{{Gm}}_{i}^{2}}{c}{\chi }_{i},\quad i=1,2,\end{eqnarray} \tag{ 1 }$$

with the dimensionless spin-parameter χ_i = V_i/c ∈ [0, 1] and the rotational velocity V_i. We assume that the rotation velocity is maximal, i.e., V_i ≈ c, leading to χ ≈ 1. Comparing S₁ and S₂ leads to

$$\begin{eqnarray}&&\displaystyle \frac{{S}_{2}}{{S}_{1}}\approx {\left(\displaystyle \frac{{m}_{2}}{{m}_{1}}\right)}^{2}={q}^{2}.\end{eqnarray} \tag{ 2 }$$

Here, we can neglect the spin of the lighter black hole as S₂ ≪ S₁ holds for most binaries, in particular concerning the range of mass ratios 0.033 ≲ q ≲ 0.33 that has been identified as most common in Gergely & Biermann (2009) and Caramete & Biermann (2010). This approximation is justified as the spins of the central SMBHs in active galaxies are usually believed to be high (e.g., Daly 2019).

Due to the spin–orbit interaction, which is a 1.5 PN-order effect, the axes of the angular momentum L and the spin S₁ start precessing around the axis of the total angular momentum J (Gergely & Biermann 2009). This is the dominant effect and is just slightly modified by additional interactions, i.e., spin–spin, mass quadrupolar, magnetic dipolar, self-spin, and higher-order spin–orbit effects, which we can neglect (Gergely & Biermann 2009). The precession process begins when the timescale of gravitational radiation becomes comparable to the timescale of dynamical friction (Zier & Biermann 2001, 2002; Bar-Or & Alexander 2014, 2016):

$$\begin{eqnarray}&&{t}_{\mathrm{dyn}}\approx \displaystyle \frac{2\,{c}^{2}}{9\,{G}^{2}}\displaystyle \frac{{r}_{\mathrm{dist}}^{3}}{{M}^{2}}{\varepsilon }^{2/3}\,{\eta }^{-1}{\left(1+q\right)}^{-1}\approx {10}^{8}{\left(\displaystyle \frac{\varepsilon }{{10}^{-3}}\right)}^{3/2}\,\mathrm{yr},\end{eqnarray} \tag{ 3 }$$

where the post-Newtonian parameter ε ≈ v²/c² is introduced. We obtain the numerical value that is close to the one where gravitational radiation becomes important at around ε ∼ 10⁻³ by using typical values for the total mass of the binary M ∼ 3 · 10⁸, and the symmetric mass ratio $\eta ={m}_{1}{m}_{2}/{\left({m}_{1}+{m}_{2}\right)}^{2}\sim 0.25$ and q = 0.33. We further use r_distr = 5 pc for the core radius of the stellar population.

As the black holes' separation gradually decreases, S₁ eventually aligns with J (see Figure 1). On the observable timescale of decades, the signature is a signal with slightly decreasing periods with constant amplitude, only subject to variations by intrinsic nutation changes. At starting and end times of this time epoch of passing over Earth, the signal is expected to be weaker and the period at which the jet points toward Earth shorter, as only a part of the jet passes over Earth. For the signal duration, toward the end of the lifetime of the merger, it could happen that the signal duration is decreased, depending on the exact geometrical setting.

Zoom In Zoom Out Reset image size

The timescale of gravitational radiation, defined here as the remaining time until merger ΔT_GW = t_merger − t, is used as a running variable as it is defined by the angular momentum L and its change rate $\dot{L}$ at a given time t:

$$\begin{eqnarray}&&{\rm{\Delta }}{T}_{\mathrm{GW}}=-\displaystyle \frac{L}{\dot{L}}\approx \displaystyle \frac{5{GM}}{32{c}^{3}}{\varepsilon }^{-4}{\eta }^{-1}.\end{eqnarray} \tag{ 4 }$$

Angular momentum loss induced by changes in the gravitational energy becomes dominant once ΔT_GW becomes comparable to the characteristic time of dynamical friction due to the dense stellar cusp in the galactic center. We can use this to set a lower bound on ε in the precessing regime (Gergely & Biermann 2009):

$$\begin{eqnarray}&&\varepsilon \gt {\varepsilon }^{\star }={\left[\displaystyle \frac{45}{64}(1+q)\right]}^{2/11}{\left(\displaystyle \frac{{GM}}{{c}^{2}{r}_{\mathrm{distr}}}\right)}^{6/11}\approx {10}^{-3}.\end{eqnarray} \tag{ 5 }$$

Our Post-Newtonian approximations of order 2.5 are sufficiently accurate for ε ≲ 0.1, corresponding to ΔT_GW ≳ 0.5 − 2.4 yr before the merger. This regime is sufficient for the purposes of our calculations below, in particular in the case of TXS 0506+056.

Once the system loses enough orbital angular momentum so that L < S₁, the angular velocity of the precession Ω_P can be expressed as (Gergely & Biermann 2009)

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{\rm{P}}}(\varepsilon )\approx \displaystyle \frac{2{c}^{3}}{{GM}}{\varepsilon }^{3}.\end{eqnarray} \tag{ 6 }$$

We use this to obtain the angle ϕ of S₁ as a function of time until merger:

$$\begin{eqnarray}&&\phi ({\rm{\Delta }}{T}_{\mathrm{GW}})=-8{\left(\displaystyle \frac{5c}{32\eta {G}^{1/3}{M}^{1/3}}\right)}^{3/4}{\rm{\Delta }}{T}_{\mathrm{GW}}^{1/4}+{\phi }_{0}.\end{eqnarray} \tag{ 7 }$$

Here, ϕ₀, the integration constant, is fixed for the specific problem.

The blazar TXS 0506+056 produced a detectable neutrino signature both in 2014/2015 (Aartsen et al. 2018b) and then in 2017 September (Aartsen et al. 2018a). We now examine the possibility that TXS 0506+056 is a SMBBH system in which the two observed neutrino emission episodes are due to the precession of the relativistic outflow from the black holes. We can use observations in the context of this scenario to make predictions on upcoming flares, and even the time of the binary merger.

We model TXS 0506+056 as a binary in which the jet of the more massive black hole is precessing around the total angular momentum J. As indicated in Figure 2, we expect a neutrino flux every time the jet points toward Earth. We use the time difference of the two potential neutrino flares to determine the current periodicity of the object as 2.78 ± 0.15 years. To determine this periodicity, we assume that the duration of the two flares is generally approximately the same and we base the length of the episodes on the one of the first flare: this one has a Gaussian width of 55 days, as indicated by the Gaussian analysis of the flare (Aartsen et al. 2018b). For the second flare, we assume that the one detected event could have arrived any time during the duration of the period in which the jet points toward the Earth. We then assume that the duration of the 2017 flare is approximately the same as the 2014/2015 flare. This gives us an uncertainty of the duration of the period of ±110/2 days = ±0.15 yrs, as the neutrino could have arrived earlier or later, with no specific preference. Assuming the emission from TXS 0506+056 to be due to a precessing jet pointing toward Earth is a viable scenario as (1) TXS 0506+056 is a source with boosted emission (Padovani et al. 2019), (2) there are first hints of a binary black hole signature in the radio regime (Britzen et al. 2019), and (3) there are two potential neutrino flares that look quite differently in neutrinos and in the multimessenger data, as discussed in the Section 1. While this could be explained in a one-zone scenario by adjusting only a few parameters (Petropoulou et al. 2020), Luo & Zhang (2020) find that it is questionable if the 2017 neutrino flare is correlated to the gamma-ray flare by jet physics. Even in such a scenario, there are different models that could lead to such a behavior of the source: (i) nutation will produce a deviation from the circular path of the jet projection, which can change the signal strength for each crossing of the jet line of sight with Earth; (ii) the local environment of neutrino production might be a different one for the first and for the second flares. That is a reasonable assumption, as around 3 yr pass in between the two flares. For example, if we assume the propagation of a plasmoid along the jet as done in many of the models, e.g., Gao et al. (2019) and Hoerbe et al. (2020), the blob moves along the jet axis for ∼ 1 pc  ≈ 3 · 10¹⁸ cm, a distance that is much larger than the size of the blob itself, which is  ∼ 10¹⁵ − 10¹⁶ cm. Thus, the two flares can happen in quite different source environments.

Zoom In Zoom Out Reset image size

These observables—i.e., the times and length of the flares—are used to fix the period of the path ϕ(t) and thus link the time variable and ϕ to the model with the two TXS outbursts. This enables us to predict the next neutrino signals and to calculate the merger time of the SMBBH. The other angles α and ϑ as indicated in Figures 1 and 2 are fixed, as the condition α + ϑ =  constant prevails (see Gergely & Biermann 2009 for a discussion). Thus, only ϑ evolves independently and determines for how long a signature can be detected. For all possible parameter settings, we find that the jet periodically points back at Earth for the next ∼15–120 yr.

Parameters that enter the calculation of ϕ(t) are the total mass of the system, which we take to be M = 3 · 10⁸ M_⊙ (Padovani et al. 2019), and the mass ratio of the two black holes, which is varied between 0.033 < q < 0.33 as discussed above.

The analytical description of the system works properly for times at which the angular momentum L is smaller than the dominant spin S₁, t ≳ t (L < S₁). For values q < 0.2, this is granted for all times. For larger values, q ≳ 0.2, we need a description at smaller times and achieve this via extrapolation.

The result of our calculation is shown in Figure 3. The blue shaded areas represent the predictions for the next two expected future periods in which the jet will point toward the Earth and when we would expect a new neutrino flare. The uncertainty band is shown in light blue. As can be seen in the figure, the new time depends on the binary mass ratio and the next flare should have arrived before the end of 2020 October. A dedicated unblinding of neutrino data for the time period 2019–2020 for TXS 0506+056 could thus reveal another flare. The exact timing of such an event would enable us to specify the mass ratio of the two black holes and with high precision predict the flare after the next one. Due to possible intrinsic variations in the jet, however, it is not clear how strong the signal will be. Figure 3 also shows the prediction for the binary merger time (red shaded area), and the expected operation time frame for LISA. For large mass ratios (q ≳ 0.15), if the origin of the two neutrino emission episodes is binary jet precession, LISA will be able to detect GWs from TXS 0506+056 and thus enable the detection of neutrinos and GWs from the same source.

Zoom In Zoom Out Reset image size

In the scenario presented above, the TXS source is a candidate SMBBH merger during the lifetime of LISA. But is this a rare, singular event, or can we expect more sources that are neutrino emitters today and that could be detected as mergers in one or two decades?

To present an order-of-magnitude estimate, we compute the expected merger rate ${R}_{\mathrm{merger}}\sim {N}_{\mathrm{SMBBHs}}/{t}_{\mathrm{merger}}^{\mathrm{av}}$, with N_SMBBH as the number of SMBBHs with their jet pointing toward Earth and ${t}_{\mathrm{merger}}^{\mathrm{av}}$ as the average timescale for the inspiral phase. The latter can be inferred from radio-interferometric data as described in Gergely & Biermann (2009), with an estimate of ${t}_{\mathrm{merger}}^{\mathrm{av}}\sim 5\cdot {10}^{6}$ yr.

Our assumption that most blazars host a SMBBH is justified by observations. Active galaxies are believed to be associated with galaxy mergers (Begelman et al. 1980). This implies the existence of SMBBHs as a standard case for these types of sources and thus also for blazars, with several detections of blazar SMBBHs summarized in Kun et al. (2017). Other galaxies without jets, like Arp 220 and IRAS 14348-1447 (Genzel et al. 2001), as well as SDSS J1010+1413 (Goulding et al. 2019), have been shown to host SMBBHs. A triple nucleus has been identified for the galaxy NGC 6240, being close to merging (Komossa et al. 2003; Kollatschny et al. 2020). For such galaxies, the next merger event could lead to a gain in spin and thus to the new formation of an active galaxy with a jet, with GW and neutrino emission being expected in close temporal correlation as predicted in Yuan et al. (2020). Our model presented in this Letter is very different, as the giant jets already exist and particle acceleration and interaction are ongoing processes, only being modified by precession.

In order to estimate the number of observable SMBBHs, we start with the sample of blazars at GHz frequencies (e.g., Eckart et al. 1986 and Eckart et al. 1987). For blazars as sources with their relativistic jets pointing toward Earth, synchrotron emission is quite intense at these frequencies and thus a good tracer for these sources (see, e.g., Muñoz et al. 2003 and references therein). About 320 sources down to 1 Jy at 5 GHz over all sky are known (Kühr et al. 1981). In the southern sky, the Parkes survey was done at a lower radio frequency, achieving a lower flux density limit (Price & Milne 1965; Otrupcek & Wright 1991).

In order to investigate the full statistics, we need to consider the redshift evolution of the sources (Biermann et al. 2014; Becker 2008). We use the Planck cosmic activity rate (Planck Collaboration 2014) as a redshift distribution of the sources, which indicates that the dominant contribution from blazars comes from the redshift range 1 to 2, with a broad peak up to z = 3 contributing to the ν−GW connection. We apply the SMBH mass function presented in Caramete & Biermann (2010).

Including all sources with a flux down to 10⁻³ Jy at 5 GHz and higher could increase the sample size from about 10² to about 10^6.5±2. If the redshift evolution is weaker than the general one seen in Planck data (Planck Collaboration 2014), this would push the redshift range of the dominant contribution to redshifts z < 1. This way, the number of sources would decrease by an order of magnitude, while the sources would produce stronger signals as they are closer.

Concerning the question on how many of these sources really are neutrino emitters, Halzen (2019) summarizes that a fraction of the population of known GeV gamma-ray blazars is sufficient to explain the neutrino background. If every source would be as strong as TXS, 5% of the sources would be enough. As the other neutrino sources are not identified yet, they are expected to be weaker, and thus the true fraction will lie significantly above 5%. Stacking-searches with IceCube show that only 27% of the GeV blazars can contribute to the IceCube diffuse neutrino flux (Aartsen et al. 2017), so it is expected that the contribution is somewhere between 5% and 27%.

We conclude that the number of active sources that we can expect in the universe to be observable right now as precessing SMBBH jets in electromagnetic waves and neutrinos should be on the order of 10⁵ < N_SMBBH < 10⁸. With the average inspiral time of ${t}_{\mathrm{merger}}^{\mathrm{av}}$, the merger rate can therefore be estimated to 0.03 yr⁻¹ < R_merger < 30 yr⁻¹.

Within these uncertainties, our conclusions are as follows: (1) it is a plausible scenario to expect several merger events with GW signatures that also show electromagnetic and neutrino emission on timescales that suit the LISA time window. (2) Those sources that are close to the merger should now reveal a time periodicity in the electromagnetic and neutrino spectra, already 15–25 yr before the actual merger. To identify such a periodic behavior is only possible if intrinsic variations of the sources are not too strong.

Concerning the periodicity P of an SMBBH with jets, for a typical mass of 10⁸M_⊙ and mass ratio q = 0.1, it lies between 0.3 yr ≤ P ≤ 3000 yr (Gergely & Biermann 2009). The binary period that changes faster for shorter periods is steeply decreasing with time (Haiman et al. 2009); therefore, the systems will spend most of their time at low periodicities that can be identified through observations (e.g., Charisi et al. 2016). For IceCube's current ∼10 yr observing time, this maximum identifiable period period is around a few years, in line with the period we find here indicated for TXS.

Explaining the two neutrino flares of the blazar TXS 0506+056 in the scenario of a jet that is precessing prior to the merger of two SMBHs leads to three key predictions:

• 1.  
  A new neutrino flare is expected any time between late 2019 and the end of 2020 October. Obviously, it is likely that such a flare already would have happened. There was no trigger in the online system of IceCube; thus, no events passed the threshold for such a trigger. As we only expect a small number of events above the atmospheric background, this is probably not a surprise—in the data, there is no significant flare in 2011/2012, when it was expected if our model is correct. For this ongoing time window, the flare can still be found via a dedicated offline analysis of IceCube data—this way, evidence for the 2014/2015 flare was found, while not visible in the online analysis. A blind analysis of the IceCube data or other detectors that search for neutrinos (see, e.g., Aab et al. 2019) with data from the direction of TXS 0506+056 could therefore reveal a flare if the signal is strong enough to be detectable.
• 2.  
  If the mass ratio of the binary system of TXS 0506+056 is q ≳ 0.15, the merger will occur at a time when it can be detected by LISA. For mass ratios q ≲ 0.15, merging will happen at later times, up to 120 yr from now.
• 3.  
  We estimate the merger rate of SMBBHs that drive precessing, relativistic outflows to be ∼1 yr⁻¹, albeit with large uncertainties. Therefore, the explanation of the neutrino detection from TXS as a precursor event for an SMBBH merger is reasonable in the light of these numbers. The example of TXS shows that the periodicity of those sources that will produce GW emission in ∼15–25 yr is around a few years. Thus, in order to identify more neutrino sources that might also be observable through GWs later, such periodicity timescales will be interesting to search for with dedicated analyses using IceCube and future telescopes like Km3NET and IceCube-Gen2.

The above predictions can open a new window to precision multimessenger astrophysics if it is possible to identify sources that show a flaring behavior as discussed here. These observations can have the power to finally provide fundamental parameters of the SMBHs as the total mass, mass ratio, and precession angle. The next years will show if the scenario of a precessing jet is viable for the system TXS 0506+056 and if other sources with jet precession can be detected even in neutrinos.

The authors would like to thank Francis Halzen and Emma Kun for valuable discussions on the topic and Roger Clay and Elisa Lohfink for commenting on the manuscript. I.B. is grateful for the support of the National Science Foundation under award number PHY-1911796 and the Alfred P. Sloan Foundation. J.B.T. acknowledges support from the DFG (grant TJ62/8-1).