High-energy Neutrino Emission Associated with GWs from Binary Black Hole Mergers in AGN Accretion Disks

The search for multimessenger signals of binary black hole (BBH) mergers is crucial to understanding the merger process of BBH and the relative astrophysical environment. Considering BBH mergers occurring in the active galactic nuclei (AGN) accretion disks, we focus on the accompanying high-energy neutrino production from the interaction between the jet launched by the postmerger remnant BH and disk materials. Particles can be accelerated by the shocks generated from the jet–disk interaction and subsequently interact with the disk gas and radiations to produce high-energy neutrinos through hadronic processes. We demonstrate that the identification of the high-energy neutrino signal from BBH merger in AGN disks is feasible. In addition, the joint BBH gravitational wave (GW) and neutrino detection rate is derived, which can be used to constrain the BBH merger rate and the accretion rate of the remnant BH based on the future associated detections of GWs and neutrinos. To date, an upper limit of BBH merger rate density in AGN disks of R 0 ≲ 3 Gpc−3 yr−1 is derived for the fiducial parameter values based on the current null association of GWs and neutrinos.

Especially, AGN accretion disk can be a promising factory of compact-stellar binaries (CNBs) including BBH systems since it generally contains various stars and compact objects, such as white dwarfs (WDs), neutron stars (NSs), and stellar mass black holes (BHs).
Due to the high density of AGN accretion disk, many stars and compact objects may be captured from nuclear star clusters or migrate from the outer self-gravitating region via gravitational instability (Syer et al. 1991;Artymowicz et al. 1993;Kolykhalov & Syunyaev 1980;Goodman & Tan 2004).Due to the rich number of stars and compact objects in the accretion disk, there can be a high probability of forming binary systems, including binary neutron stars (BNSs), neutron star-black holes (NSBHs), and BBHs (Bartos et al. 2017;Leigh et al. 2018).
The simultaneous detections of GW and electromagnetic (EM) signals have been long expected.Unlike BNS and NSBH mergers, the BBH mergers are usually believed not to generate EM radiations due to a lack of accretion materials to power a jet.However, BBH mergers embedded in AGN disks can potentially power EM counterparts by accreting a significant amount of disk gas and then interacting with the disk medium (Bartos et al. 2017;McKernan et al. 2019;Wang et al. 2021;Kimura et al. 2021;Tagawa et al. 2022).In the observational aspect, EM follow-up observations have been implemented for LIGO/Virgo BBH merge events and seven potential AGN flares statistically associated with one or more of nine LIGO/Virgo events have been sug-gested (Graham et al. 2023).In particular, the optical flare ZTF19abanrhr with a luminosity of ∼ 10 45 erg/s detected by the Zwicky Transient Facility (ZTF) in the AGN J124942.3+344929 was reported to be spatially coincidence with GW190521 (Graham et al. 2020), which is the heaviest BBH merger with a total mass of M BBH ∼ 100M ⊙ detected so far (Abbott et al. 2020).
Although the joint observation of EM and GW signals is exciting and will help to better study the BBH merger process, it is still under debate whether the accompanying EM flare can be identified from the bright AGN disk emission.For instance, if the shock is too weak, the subsequent EM emission is too dim to be identified against bright AGNs and can be disclosed only for lowerluminosity AGNs (McKernan et al. 2019).Besides, premerger outflows may also create a cavity around the BBH merger remnant to make the accretion rate insufficient to produce high luminosity (Kimura et al. 2021).
In this letter, we propose that the high-energy neutrino can be an alternative probe for the joint multimessenger study for BBH merger embedded in AGN disks.Due to the negligible absorptions by AGN disk materials, the high-energy neutrinos can easily escape from the complicated AGN disk environment.In addition, unlike the bright EM emission background, the AGN disk generally does not produce high-energy neutrino emission unless the shocks are generated by the perturbation of disk materials induced by some stellar activities therein (Zhu et al. 2021;Zhu et al. 2021a;Zhou et al. 2023).However, the study of the highenergy neutrino production mechanism for BBH mergers insider AGN disks is still lacking although some searches of associated high-energy neutrino events have been explored (Abbasi et al. 2023).In our scenario, the merged BH remnant is kicked out from the cavity and into the intact AGN disk environment and then accretes surrounding gas at the Bondi-Hoyle-Lyttleton rate (Kimura et al. 2021;Comerford et al. 2019;Wang et al. 2021;Graham et al. 2023;Tagawa et al. 2023).The interaction between the post-merger Blandford-Znajek jet and the AGN disk medium can create various shocks, which in turn accelerates high-energy cosmic rays and produce neutrinos through hadronic processes.It is also expected to produce a detectable EM signal when shock breaks out from the disk (Wang et al. 2021;Tagawa et al. 2023).Therefore, the BBH merger inside AGN disks can be the potential triplet messenger (GW, EM, and neutrino) source.

JET STRUCTURE AND SHOCK ACCELERATION
BBH mergers are expected to occur at about r ∼ 10 3 R g in the migration traps (Bellovary et al. 2016), where R g = GM SMBH /c 2 is the Schwarzschild radius of central supermassive BH (SMBH), G is the gravitational constant, M SMBH is the SMBH mass, and c is the speed of light.For a gas-pressure-dominated disk, the disk density is ρ (h) = ρ 0 exp (−h/H), where ρ 0 is the density of mid-plane, h is the vertical distance, and H is the typical disk height.For a SMBH with mass 10 8 M ⊙ and the disk with aspect ratio r/H ∼ 0.01, one can get H = 1.5 × 10 14 cm.The mid-plane density near the migration traps is approximately 10 −10 g cm −3 .
A cocoon usually accompanies a jet in a dense environment.We consider a Blandford-Znajek jet launched from rapidly accreting and spinning remnant BH in an AGN disk.The Bondi-Hoyle-Lyttleton rate is adopted as (Comerford et al. 2019;Wang et al. 2021 (1) Through simulation using a hierarchical population analysis framework, (Gayathri et al. 2023) considered binary formation in AGN disks along with phenomenological models and found that the high-mass and highmass-ratio binaries appear more likely to have an AGN origin.By referring to the GW event parameters (Gayathri et al. 2021), we adopt a total mass of M BBH ∼ 100M ⊙ for the BBH and assume that the relative velocity is given by v rel = v k + c s .Here we use a gas sound speed of c s ∼ 50km s −1 and a kick velocity of v k ∼ 200 km s −1 (Graham et al. 2020).The jet kinetic luminosity L j can be expressed as L j = η j ṁc 2 , where η j = 0.5 is the adopted jet conversion efficiency and the accretion rate can be parameterized as ṁ = f acc ṁBHL .We will use f = [0.1, 1, 10] respectively in our calculations.Then the jet kinetic luminosity can be written as (2) where M ⊙ is the solar mass and ρ 0,−10 = (ρ 0 /10 −10 )g cm −3 .
The jet is initially uncollimated because the pressure of the cocoon on the jet is not high enough, it will become collimated as it travels.As suggested by (Bromberg et al. 2011), the jet will be collimated if , where θ 0 ≈ 0.2 is the jet opening angle and ρ  We study the high-energy neutrino emission from four possible sites, i.e., internal shocks, collimation shocks, forward shocks, and reverse shocks, formed during the relativistic jet of post-merger BH interacting with the disk atmosphere (Yuan et al. 2020).Fig. 1 schematically describes the evolution of the structure of the jet-cocoon system as well as the shocks inside the jet.Particle acceleration is driven by these shocks, and high-energy neutrinos are produced by pp and pγ interaction processes.
If the ambient medium is non-relativistic and the reverse shock is strong, the velocity of the jet head is At the initial time t 0 , one can get the initial height of the jet head h 0 ≈ ct 0 , the initial average density of the disk ρ0 = ρ (h 0 ), the initial L0 = Lj ρ0t 2 0 θ 2 0 c 5 and the initial jet head velocity β h,0 = βj 1+ L0 −1/2 .Therefore, the dynamic evolution process of the jet can be calculated with time.Noted that L = . Through this method, we obtained the evolution of jet head height R h = h, L, jet head velocity β h over the jet eruption time.The pressure of the cocoon in the collimated regime can be written as P c ≃ t −4/5 L 2/5 j ρ3/5 θ 2/5 0 , then the height of the collimation shock can be calculated by We calculate the structural evolution of the jet and find the jet becomes collimated at 100 s after the jet eruption, which is a quite short time compared with the duration of a BZ jet.At this moment, the height of the jet head now is R h ∼ 10 12 cm ≪ H for the fiducial f acc = 1.We assume the Lorentz factor of the unshocked material in the pre-collimation region to be Γ j , In this particular region, internal shocks emerge due to fluctu-ations in velocity within the outflow, resulting in the creation of gas shells exhibiting differential speeds.We may approximate the height of the internal shocks to be R is = min R cs , 2Γ 2 j ct var (Yuan et al. 2020), where t var ≃ 1 s is the adopted variability timescale.Shock acceleration can be efficient only if the shock is collisionless.Therefore, we can obtain a radiation constraint on the upstream of the shock for efficient Fermi acceleration, which is described as (Murase & Ioka 2013) where τ u is the upstream optical depth, n u is the comoving number density of upstream material, σ T is the Thomson cross section, l u is the length scale of the upstream fluid, Γ rel represents the relative Lorentz factor between the shock downstream and upstream.It means that an efficient particle acceleration can only occur when the shock has a sufficiently strong jump between the upstream and downstream.
To better represent the velocity relationship between different regions of the jet, we consider an idealized jet diagram in Fig. 1 where the upstreams of the collimation shock and the reverse shock are downstream of the internal shock and the collimation shock respectively.The jet head is downstream of the forward shock and the reverse shock.Roughly considering that the position of the collimation shock is always at the position where the jet begins to collimate, we can write the comoving number density of the upstream as n cs,u = L iso /(4πΓ 2 rel R 2 cs m p c 3 ), where L iso ≈ 2L j /θ 2 0 is the isotropic equivalent one-side jet luminosity and m p is the mass of proton.Here Γ j ∼ 20 we used is the Lorentz factor of the unshocked material.Upstream optical depth can be calculated as (5) The relative Lorentz factor of internal shock between the shock downstream and upstream is Γ rel,is ≈ Γ r /2Γ j ≈ 5, and therefore upstream optical depth internal shock then can be derived by Similarly, one can get the upstream optical depth for the reverse shock as where Γ rel,cs ≈ Γ j /2Γ j1 ≈ Γ j θ 0 /2 ≃ 2. For the forward shock, the upstream optical depth is Comparing with Equation 2, we find that the forward shock is almost always inefficient in accelerating particles, resulting in inefficient neutrino production at the forward shock site, while the other three shocks always efficiently accelerate particles.Therefore, the neutrino production from forward shock is neglected in the next calculations.

NEUTRINO PRODUCTION
By the shock jump conditions (Piran et al. 1995), we can describe the internal energy and density evolution of the upstream and downstream shock by e sh,d /n sh,d m p c 2 = Γ rel − 1 and n sh,d /n sh,u ≈ 4Γ rel , where n sh,u and n sh,d are the comoving number density in the upstream and downstream shock respectively, e sh is the comoving internal energy density in the shock.So the energy density of downstream of the shock is e sh,d = 4Γ rel (Γ rel − 1)n sh,u m p c 2 .The photon temperature of the downstream of the shock is k B T = (15 3 c 3 ε e e sh,d /π 2 ) 1/4 , where ε e ≈ 0.1 is the electron energy fraction (Zhu et al. 2021a).During the next calculations, these thermal photons are treated as the background photon field for inverse Compton scattering (ICS), photomeson production pγ, and Bethe-Heitler processes.
To calculate the neutrino emission efficiency, we need to estimate the cooling and acceleration timescales of the protons.The acceleration timescale is given by t p,acc = ǫ p /eBc, where B = 8πε B e sh,d is the downstream magnetic field intensity and the magnetic field energy fraction is adopted as ε B = 0.1.High energy protons cooling mainly include synchrotron radiation and ICS as radiative processes, pp, pγ, and Bethe-Heitler processes as hadronic interaction processes, and finally, the adiabatic process.
For radiative processes, the cooling timescale of synchrotron radiation is and ICS has the cooling timescale where ε γ = 2.7k B T is the average thermal photon energy downstream and n γ = ε e e sh,d /ε γ is the average thermal photon density downstream.For hadronic cooling mechanisms, the pp scattering timescale is given by t pp = 1/cσ pp n p κ pp , where σ pp (Kelner et al. 2006) is the cross section and κ pp ≃ 0.5 is the inelasticity.The cooling efficiency of pγ process can be calculated by (Murase 2007) dǫǫ −2 dn dǫ (11) where γ p = ǫ p /mc 2 , ǭ is the photon energy in the rest frame of the proton and dn dǫ is the photon number density.ǭth,pγ = 145MeV is the threshold energy for pγ process, σ pγ and κ pγ represent the cross section (Kelner & Aharonian 2008) and inelasticity (Stecker 1968), respectively.By replacing the cross section, inelasticity, and threshold energy in the Equation 11 with those of the Bethe-Heitler process, one can get the cooling efficiency of the Bethe-Heitler process (Chodorowski et al. 1992).Finally, the adiabatic cooling timescale is t p,ad = R sh /cΓ sh,d .
Considering that only pp and pγ processes can produce high-energy neutrinos, the other processes suppress the production of neutrinos.We can write the proton suppression factor by involving various cooling processes as (12) Besides, neutrinos are produced by the decay of pions and kaons created through pp and pγ processes, which will be suppressed by other cooling processes.The suppression factor of these mesons can be calculated by (Zhu et al. 2021a) where i represents the meson produced by the pp or pγ processes.Neutrino fluence for a single event can be obtained by calculating the summation of each neutrino channel by ln (ǫ p,max /ǫ p,min ) dt (14) where N i is the energy fraction that the jet energy converts to the neutrinos, and i represents different neutrino production channels.We assume t end = 5 × 10 6 s is the duration of the jet corresponding to the observational ZTF19abanrhr flare of GW190521, which lasts around tens of days (Graham et al. 2020).ǫ p,max is the maximum neutrino energy calculated by t acc = t cool and ǫ p,min ≈ Γ sh,d m p c 2 is the minimum neutrino energy.A collimation shock at the base of the jet is evident in the simulations whenever the jet is collimated by the surrounding cocoon.
The produced neutrino fluence of a single BBH merger event occurred at D L = 100 Mpc is shown in Figure 2 including the reverse shock, collimation shock, and internal shock.We can see that the neutrino production of the collimation shock is comparable with that of the reverse shock above around few×10 6 GeV, while below this energy the reverse shock contributes more neutrino production.For three shock sites, at the lower energy part ( 10 5 GeV), all neutrino production suffers the significant suppression of the adiabatic cooling process, showing a lower fluence than the high-energy part.Note that here we consider the same jet duration t end for different accretion rates, however, it may last a shorter time for higher accretion rate (Wang et al. 2021).
Figure 3. Cumulative detected neutrino number of a single source by IceCube with different accretion rates at two different locations, i.e., 100 Mpc and 200 Mpc.The green line indicates three detected neutrinos so that the neutrino detection by IceCube can be identified at ∼ 3σ significance level.

NEUTRINO AND JOINT GW+NEUTRINO DETECTION
The all-flavor neutrino detection number can be caculated by where A eff ǫ νµ is the effective area (100 GeV-100 PeV) of IceCube for a point source (Aartsen et al. 2020).The accumulative neutrino number with time is presented in Fig. 3, it can be seen that the Ice-Cube can receive at least three neutrinos within 5 × 10 6 s within a distance of 200 Mpc.In a relatively optimistic situation (f acc = 10 and D L = 100Mpc), we can expect Ice-Cube to receive three neutrinos within 6 hours which has reached a relatively high level of confidence.In the case of a high accretion rate (f acc = 10), the detection distance of the IceCube can reach 1 Gpc.However, as shown in Fig. 3, for the high accretion rate, the particle acceleration tends to be forbidden at the early stage so that the neutrino production is low at the beginning.
We calculate the joint detection rate of neutrinos and GW by generating 3 × 10 4 random BBH merger events.f acc = 1 is used in the calculation.GW signal is treated as isotropic and high-energy neutrinos from the jet are beamed with a beaming correction factor f b = θ 2 0 /2.As for the mass distribution of BBH, we adopt the results of (Gayathri et al. 2023), which presented a one-parameter Figure 4. Expected all-flavor diffuse neutrino fluence contributed from BBH mergers in AGN disks.Diverse local event rates and accretion rates are considered.Red points and upper limits are the observed diffuse neutrino fluence by IceCube (Aartsen et al. 2015(Aartsen et al. , 2021a)).The gray lines are 90% upper limit of the cosmogenic neutrino for diverse instruments.
model for BBH formation and merger within an AGN disk and parameterized by the maximum mass m max of the natal BH distribution.The mass distribution for m max = 75 is adopted for calculations.We use Python's module PYCBC to generate the waveform of the GW signal and calculate the signal-to-noise ratio (SNR) of each event by (S/R) Sn(f ) df (Zhu et al. 2021b).The threshold SNR=8 is involved to confirm the detection of GW signals, and the neutrino detection number N ν = 1 is roughly employed as the threshold for confirming the detection of IceCube or IceCube-Gen2.The effective area of IceCube-Gen2 is taken as 6 times larger than IceCube's (Aartsen et al. 2021b).The redshift evolution factor we adopted is f (z) = (1 + z) 5.7η + 1+z 0.36 1.3η + 1+z 3.3 −9.5η + 1+z 3.3 −24.5η 1/η (Sun et al. 2015).For multiple gravitational wave detectors, we calculated the joint detection rates separately.The detector sensitivities can be expressed as amplitude spectral densities (ASD, https://dcc.ligo.org/LIGO-T1500293/public).The result is shown in Table 1, where R 0 the local BBH mergers event rate.It can be seen that although different GW detectors have different detection capabilities, the calculated joint detection rates are largely identical with only minor differences for the same neutrino detector.This indicates that the detection ability of the neutrino detector is the main factor affecting the joint detection rate, and indeed, the joint detection rate significantly increases for IceCube-Gen2 with a larger effective area.
Several studies have given the BBH merger event rates.The estimated GW rate density associated with BBH mergers lies in the range R ∼ (0.002 − 18) Gpc −3 yr −1 is given by (Gröbner et al. 2020).Recently, some literature gives the event rate to be 2.5 Gpc −3 yr −1 or dozens Gpc −3 yr −1 (Gayathri et al. 2023(Gayathri et al. , 2021)).Here three typical values of [0.002,0.2,20]as the possible local event rates are used.We sum the neutrino fluences of each type of shock and use the mass distribution of BBH mentioned above.Our result shows that BBH mergers in AGN disk contribute to the neutrino background at the relatively high energy above 10 6 GeV, while they provide a relatively little contribution to the neutrino background at the low energy.

DISCUSSIONS AND CONCLUSIONS
Neutrinos from remnant BH jet-induced shock acceleration make it possible to predict a BBH merger without the detection of an EM signal if the jet breakout brightness is overshadowed by the AGN accretion disk.In this letter, we consider a persistent BZ jet and investigate its structural evolution.Collimation shocks and reverse shocks contribute more neutrino production compared with internal shocks.The accumulative neutrino number over time is derived as well, and we find that it is possible to receive enough neutrinos (exceeding three neutrinos) by IceCube within tens of days if the BBH mergers in AGN disk take place within a distance of a few hundred of Mpc and the remnant BH has an accretion rate of f acc 0.1.Note that the jet eruption duration will affect the neutrino detection significantly.Although a comparable observational duration of ZTF19abanrhr flare, i.e., tens of days, is involved for evaluation, the real jet eruption duration or accretion timescale may be shorter or even episodic (Wang et al. 2021).For a shorter jet duration, e.g., 10 5 s, the neutrino detection of a single BBH merger even by IceCube is feasible only for a higher accretion rate or closer distance.
Based on 83 LIGO/Virgo BBH and lower-mass-gap merger alters (Graham et al. 2023), the number of GW events within 500 Mpc is ∼ 7. Considering the beaming correction, the source number with the neutrino detection by IceCube for these GW events can be estimated by ∼ 7f b ≃ 0.14(θ 0 /0.2) 2 if all BBH mergers occur in AGN disks and the neutrino emission from theses source can be identified, which is consistent with the current null association of GW and neutrino events recently reported in Abbasi et al. (2023).
In addition, we calculate the joint GW + Neutrino detection rate by combining the diverse GW detectors and IceCube (or IceCube-Gen2).The uncertainty of joint detection is still quite large due to the event rate of BBH merger in AGN disks is still quite unclear.However, the detection or non-detection of joint BBH GW + neutrino association in the future can be used to constrain the BBH merger rate within AGN disks.For instance, the BBH merger rate in AGN disks should be lower if no association has been observed.Based on our calculations of diffuse neutrino background from BBH mergers in AGN disks (Fig. 4), high BBH merger rate and high accretion rate of remnant BH tend to be excluded due to the conflict with the neutrino observations.BBH mergers in AGN disks are the ideal targets for multi-messenger observations and joint observations of EM, neutrino, and GW signals from them have caused more and more attention and have also been frequently explored recently.
In the aspect of high-energy neutrinos, they can play an important role in the multi-messenger study of BBH mergers.The next-generation neutrino telescopes, e.g., IceCube-Gen2, Huge Underwater high-energy Neutrino Telescope (HUNT) (Huang et al. 2023), The tRopIcal DEep-sea Neutrino Telescope (TRIDENT) (Ye et al. 2022), and the radio-Cherenkov neutrino detector AR-IANNA (Anker et al. 2020) and ARA (Allison et al. 2019) combining with the next more advanced GW detectors can help to understand the nature of BBH mergers in AGN disks.
We thank Jin-Ping Zhu for helpful discussions.Note that when the paper is ready to be submitted we noticed Jin-Ping Zhu also submitted a work on a similar topic, but both works are carried out independently.We acknowledge support from the National Natural Science Foundation of China under grant No.12003007 and the Fundamental Research Funds for the Central Universities (No. 2020kfyXJJS039).
dh is the average density of the medium.

Figure 1 .
Figure 1.The three pictures on the left are schematic pictures of jet propagation in the AGN disk.There are four types of shocks in the jet, including forward shock (FS), reverse shock (RS), collimation shock (CS) and internal shock (IS).The picture on the far right is the schematic picture of an idealized jet used in our calculations.Each color corresponds to a same Lorentz factor in that region.The shocks are located at the junction of each shock area.