Astroparticles from X-ray Binary Coronae

K. Fang Department of Physics, Wisconsin IceCube Particle Astrophysics Center, University of Wisconsin, Madison, WI, 53706 Francis Halzen Department of Physics, Wisconsin IceCube Particle Astrophysics Center, University of Wisconsin, Madison, WI, 53706 Sebastian Heinz Department of Astronomy, University of Wisconsin, Madison, WI, 53706 John S. Gallagher Department of Astronomy, University of Wisconsin, Madison, WI, 53706 Department of Physics and Astronomy, Macalester College, St. Paul., MN, 55105

Abstract

The recent observation of high-energy neutrinos from the Galactic plane implies an abundant population of hadronic cosmic-ray sources in the Milky Way. We explore the role of the coronae of accreting stellar-mass black holes as such astroparticle emitters. We show that the particle acceleration and interaction timescales in the coronal region are tied to the compactness of the X-ray source. Thus, neutrino emission processes may similarly happen in the cores of active galactic nuclei and black hole X-ray binaries (XRB), despite of their drastically different masses and physical sizes. We apply the model to the well-measured XRB Cygnus X-1 and find that the cascaded gamma rays accompanying the neutrino emission naturally explain the GeV emission that only presents during the source’s hard state, while the state-averaged gamma-ray emission explains the LHAASO observation above 20 TeV. We show that XRB coronae could contribute significantly to the Galactic cosmic-ray and Galactic plane neutrino fluxes. Our model predicts variable high-energy neutrino emission from bright Galactic XRBs that may be observed by IceCube and future neutrino observatories.

1 introduction

Neutrino emission from the Galactic plane was recently identified at the $4.5\sigma$ level of significance using ten years of data from the IceCube Observatory (IceCube Collaboration, 2023). The neutrino flux is consistent with diffuse neutrino emission by the cosmic-ray sea of the Milky Way but could also arise from a population of unresolved sources (Fang & Halzen, 2024). When converting the neutrino flux to a pionic gamma-ray flux, the IceCube Galactic plane emission is consistent with the Galactic diffuse emission (GDE) measured by Fermi-LAT at around 1 TeV (Ackermann et al., 2012) and the Tibet AS $\gamma$ Observatory above 100 TeV (Tibet AS ${\gamma}$ Collaboration et al., 2021). Above $\sim 30$ TeV, it is comparable to, or higher than, the flux sum of the GDE and individual sources that are not associated with pulsars measured by the Large High Altitude Air Shower Observatory (LHAASO) (Fang & Murase, 2023). This implies two possibilities: 1) the GDE and a significant fraction of the LHAASO sources originate from hadronic interactions; or 2) there exists a population of $\gamma$ -ray-obscured neutrino sources in the Milky Way.

Gamma-ray-obscured sources are evidently present in the extragalactic neutrino sky. The energy flux of the diffuse cosmic neutrino emission is higher than that of the isotropic $\gamma$ -ray emission (Murase et al., 2016; Fang et al., 2022). This leads to a tension between the Fermi-LAT observation (Ackermann et al., 2015) above $\sim 10$ GeV and the pionic $\gamma$ rays that leave the neutrino sources without attenuation and cascade down to GeV-TeV energies during their propagation in the extragalactic background light. Individual neutrino source observations also hint a source environment that is optically thick to GeV-TeV $\gamma$ rays. The neutrino flux of the active galaxy NGC 1068 is more than an order of magnitude higher than its GeV-TeV $\gamma$ -ray flux (Abbasi et al., 2022a), suggesting that the neutrino emission site must be highly $\gamma$ -ray-obscured.

The coronal region in the proximity of the central black hole provides an appealing site for both neutrino emission and $\gamma$ -ray attenuation. Models of neutrino emission from the cores of active galactic nuclei (AGN) have long been proposed (Berezinskii & Ginzburg, 1981; Stecker et al., 1991, 1992) and recently been developed to explain the NGC 1068 observations (Inoue et al., 2020; Murase et al., 2020; Inoue et al., 2021; Murase, 2022).

Despite of the vast difference in scales, stellar-mass black hole binaries such as Cygnus X-1 are close analogues to luminous AGNs as a result of their similar compactness (Merloni et al., 2003; Falcke et al., 2004; Fabian et al., 2015). The coronal compactness is defined as (Guilbert et al., 1983):

\ell=\frac{L_{X}}{R}\frac{\sigma_{T}}{m_{e}c^{3}}=4\pi\frac{m_{p}}{m_{e}}\frac{R_{g}}{R}\frac{L_{X}}{L_{E}}

(1)

where $L_{X}$ and $R$ are the luminosity and radius of the X-ray source, $R_{g}\equiv GM_{\rm BH}/c^{2}$ is the gravitational radius, $L_{E}=4\pi GM_{\rm BH}m_{p}c/\sigma_{T}$ is the Eddington luminosity, $m_{p}$ and $m_{e}$ are the mass of proton and electron, respectively, and $\sigma_{T}$ is the Thomson cross section. For comparison, at $R=10\,R_{g}$ , the compactness of NGC 1068 is $\ell\approx 17$ with an intrinsic X-ray luminosity $4.6\times 10^{43}\,\rm erg\,s^{-1}$ (Bauer et al., 2015) and a supermassive black hole mass of $\sim 5\times 10^{7}\,M_{\odot}$ (Minezaki & Matsushita, 2015), very similar to that of the Galactic X-ray binary (XRB) Cygnus X-1, which has $\ell\approx 19$ with an X-ray luminosity $2.2\times 10^{37}\,\rm erg\,s^{-1}$ during the hard state (McConnell et al., 2002) and a stellar mass black hole of $21.2\,M_{\odot}$ (Miller-Jones et al., 2021).

The shorter timescales of accretion on to black hole XRBs provide opportunities to study different accretion states and the associated high-energy radiation. Gamma rays have been observed from several XRBs, including LS 5039 (Aharonian et al., 2005, 2006), Cygnus X-3 (Fermi-LAT Collaboration, 2009)), and Cygnus X-1 (Zdziarski et al., 2017; Abdollahi et al., 2020). Orbital modulation and absorption by the radiation of the companion star is observed, suggesting that the $\gamma$ -ray emission must originate from a region smaller than the binary separation (Dubus, 2013). Neutrino emission from XRBs have previously been proposed in the context of jets interacting with matter of the wind or companion star (Stecker et al., 1985; Romero et al., 2003; Bednarek, 2005). Below we focus on the canonical and well-measured Cygnus X-1 and a scenario where its high-energy emission arises from the accretion flow.

2 The Variable Cygnus X-1

Cygnus X-1 is a persistent black hole XRB at a distance of 2.22 kpc (Miller-Jones et al., 2021). Its X-ray emission varies on timescales of several weeks between two discrete levels: the “hard state”, during which the source spends $\sim 90\%$ of its time, characterized by a relatively low flux of soft X-rays ( $\sim 2-10$ keV) and a relatively high flux of hard X-rays ( $\sim 20-100$ keV); and the “soft state”, characterized by a relatively high soft X-ray flux and a relatively low hard X-ray flux (see Figure 2 and 3 for the X-ray spectra during the two states).

The transition of the states is believed to relate to different configurations of the accreting matter (McConnell et al., 2002; Done et al., 2007). A general theoretical picture of the accretion flow structures in Cygnus X-1 includes an inner hot, optically thin, geometrically thick accretion flow (the corona) and an outer cool, optically thick, geometrically thin disk. During the hard state, the source has low accretion rates. The accretion flow extends out to $\sim 100\,R_{g}$ and significantly contributes to the hard X-ray emission. Conversely, when the mass accretion rate increases, the flow becomes optically thick and collapses into a Shakura–Sunyaev disk. This explains why the blackbody radiation of the disk dominates the X-ray spectrum during the soft state.

While the soft and hard X-ray emission may be well reproduced by a multicolour disk model and hard thermal Comptonized spectrum, respectively, Cygnus X-1 also presents a high-energy tail above the thermal-Compton spectrum starting at $\sim 0.1-1$ MeV. This component is explained as the scattering of nonthermal electrons accelerated by the accretion flow (McConnell et al., 2002; Romero et al., 2014) or a compact jet (Cangemi et al., 2021). At even higher energies, $\gamma$ -ray emission at $0.1-10$ GeV is detected from Cygnus X-1 during its hard state (Zdziarski et al., 2017; Abdollahi et al., 2020). Analyses of the Fermi-LAT data during its soft state only found upper limits above 100 MeV, but obtained a tentative detection of a soft spectrum at the 40–80 MeV range (Zdziarski et al., 2017).

Various models have been implemented to explain the $\gamma$ -ray emission by Cygnus X-1. All of them invoke relativistic electrons or protons produced in the source’s steady or transient radio jets. In a leptonic scenario, electrons accelerated in the jets produce GeV emission by up-scattering the blackbody photons from the companion star (Zdziarski et al., 2014), the thermal-Compton emission by the accretion flow, or their own synchrotron radiation in the jets (Malyshev et al., 2013; Zdziarski et al., 2017). These leptonic models are subjected to various constraints. For example, as the GeV flux of Cygnus X-1 does not present a significant modulation (Zdziarski et al., 2017), the inverse Compton emission of stellar photons cannot dominantly contribute to the $\gamma$ -ray emission. A synchrotron self-Compton model requires a highly clumped jet, while an external Compton model on disk photons needs a jet emission region very close to the accretion disk such that the synchrotron flux is comparable to the $\gamma$ -ray flux (Zdziarski et al., 2017). Lepto-hadronic and hadronic scenarios have also been proposed, where non-thermal protons are accelerated in the jets in addition to electrons, producing GeV or TeV $\gamma$ rays by interacting with thermal protons or photons in the jet and in the stellar wind of the companion star (Pepe et al., 2015; Papavasileiou et al., 2023; Kantzas et al., 2023).

In this letter, we study an alternative scenario inspired by the modeling of extragalactic neutrino sources like NGC1068. We propose that relativistic protons are produced in the vicinity of the black hole, interact with the gas and photons in the coronal region, and produce neutrinos and $\gamma$ rays. In the hard state, 0.1-100 GeV $\gamma$ rays may partially survive as a result of the lower density of the disk photons. In the soft state, TeV $\gamma$ rays pair produce with the disk and coronal photons and cascade down to 1–100 MeV energies. Our model naturally connects the variation of the MeV-GeV $\gamma$ -ray flux to the transition of the accretion states. It also suggests that XRBs may be the hadronic accelerators that dominate the cosmic ray and high-energy neutrino emission of the Galaxy (Sudoh & Beacom, 2023a, b).

3 Coronal Neutrino and $\gamma$ -ray production

3.1 Proton Acceleration and Escape Timescales

Proton acceleration in advection-dominated accretion outflows has long been suggested (for example see the review of Yuan & Narayan, 2014) and more recently studied in the context of NGC 1068 (Murase et al., 2020; Mbarek et al., 2023). Processes like turbulent dissipation may accelerate a fraction of electrons and protons in the two-temperature plasma of the corona into a non-thermal power-law distribution. In addition, a reconnection layer could form in the vicinity of a black hole, power the coronal X-ray emission by rapidly converting the magnetic energy into particle energy, and drive the hard/soft state transitions (e.g., Dexter et al., 2014). Below we explore a simple scenario where particles get accelerated and interact inside the corona. The calculation would also be valid for a scenario where a population of protons get accelerated elsewhere and injected into the coronal region.

We model the corona as a spherical region with radius $R\equiv{\cal{R}}_{\rm rel}R_{g}$ , where $R_{g}=3.1\times 10^{6}\,(M_{\rm BH}/21.2\,M_{\odot})\,{\rm cm}$ , and assume that the coronal emission is subject to irradiation with soft photons from the outer accretion disk, which is assume to subtend a solid angle $\sim 2\pi$ . The magnetic energy density in the corona may be estimated as $u_{B}=B^{2}/8\pi=\xi_{B}\,u_{X}$ , where $u_{X}=L_{d}/2\pi R^{2}c+L_{c}/4\pi R^{2}c$ is the X-ray energy density radiated by the disk and corona. The scaling parameter $\xi_{B}\lesssim 1$ in regions dominated by magnetic turbulence (Mbarek et al., 2023) and $\xi_{B}\sim 1/\beta_{\rm rec}$ in the case of magnetic reconnection, where $\beta_{\rm rec}$ is the plasma inflow speed in the reconnection layer in the unit of the speed of light, and $\beta_{\rm rec}\sim 0.1$ in the collionless relativistic regime (Sironi et al., 2015). The field strength in the corona of Cygnus X-1 can then be estimated as $B\approx 1.2\times 10^{6}\,\xi_{B}^{1/2}L_{X,37.3}^{1/2}{\cal R}_{\rm rel,1}^{-1}\,{\rm G}$ .

Comparisons of thermal Comptonization models with observations also suggest $B\lesssim 10^{6}$ G in the soft state and $B\sim 10^{5}-10^{7}$ G in the hard state (Del Santo et al., 2013).

Refer to caption — Figure 1: Timescales of proton acceleration (solid curves) and escape (dash-dotted curves) in the units of light crossing time $t_{\rm cross}\equiv R/c$ as a function of the proton energy. In the hard state (blue), stochastic acceleration (equation 3 with $\eta=10$ ) and particle diffusion (equation 4) in a turbulent magnetic field are assumed to occur in a region of size $R=100\,R_{g}$ . The magnetic field strength is computed using the X-ray luminosity in Supplementary Table 1 for the assumed coronal size. In the soft state (red), magnetic reconnection (equation 5 and 6 with $\beta_{\rm rec}=0.1$ is assumed. The ratios of the timescales to the light crossing time in this scenario do not depend on the coronal size.

The coronal electron number density is $n_{e}\approx\tau_{T}/\sigma_{T}R=5\times 10^{16}\,\tau_{T}{\cal R}_{\rm rel,1}^{-1}\,\rm cm^{-3}$ , where $\tau_{T}$ is the Thomson optical depth, and $\tau_{T}\sim 2$ in hard state and $\tau_{T}\sim 0.1$ in soft state (McConnell et al., 2002). The pair magnetization may be parameterized as a function of the compactness,

\displaystyle\sigma_{\pm}=\frac{B^{2}}{4\pi n_{e}m_{e}c^{2}}=\frac{\xi_{B}}{2\pi\,\tau_{T}}\ell.

(2)

The proton magnetization is then $\sigma_{p}={B^{2}}/{4\pi n_{p}m_{p}c^{2}}=\sigma_{\pm}n_{e}m_{e}/n_{p}m_{p}$ . The proton number density is less well known. In a plasma where ions do not appreciably contribute to the mass density, $n_{p}m_{p}\ll n_{e}m_{e}$ . In a mildly relativistic scenario $\langle\gamma_{p}\rangle\sim\sigma_{p}\sim 1$ , $n_{p}m_{p}\sim n_{e}m_{e}$ .

Stochastic proton acceleration may occur in the turbulent accretion outflow over a timescale $t_{\rm acc}^{\rm turb}\approx\eta(R/c)\,(c/v_{A})^{2}\,(r_{L}/R)^{2-q}$ (Murase et al., 2020), where $\eta$ is a parameter related to the turbulence level $\delta B/B$ , $r_{L}=E_{p}/eB$ is the Larmor radius, $v_{A}=c(\sigma/(1+\sigma))^{1/2}$ is the Alfvén velocity, $\sigma=B^{2}/(4\pi(n_{p}m_{p}c^{2}+n_{e}m_{e}c^{2}))$ is the magnetization parameter, and $\sigma\sim\sigma_{\pm}$ in a plasma where ions do not dominate the mass density. Taking $q=5/3$ for Kolmogorov turbulence and comparing to a light crossing time $t_{\rm cross}\equiv R/c=2.0\times 10^{-3}\,{\cal R}_{\rm rel,1}\,\rm s$ , we obtain

\displaystyle\frac{t_{\rm acc}^{\rm turb}}{t_{\rm cross}}=\eta\frac{1+\sigma}{\sigma}\left(\frac{E_{p}}{E_{m}}\right)^{1/3}

(3)

where $E_{m}\equiv eBR\approx 9.4\,B_{6}{\cal R}_{\rm rel,1}\,\rm PeV$ .

Accelerated particles diffuse away from the turbulent regions. The duration that protons are confined in the turbulent magnetic field can be estimated as $t_{\rm esc}^{\rm turb}\sim(R/c)(r_{L}/eBR)^{-1/3}$ (Effenberger & Petrosian, 2018; Murase et al., 2020; Mbarek et al., 2023), therefore,

\frac{t_{\rm esc}^{\rm turb}}{t_{\rm cross}}\approx\left(\frac{E_{p}}{E_{m}}\right)^{-1/3}.

(4)

In the absence of cooling processes, ions accelerated in the turbulent magnetic field may reach up to $t_{\rm acc,turb}<t_{\rm esc,turb}$ , or $E_{p,\rm max}^{\rm turb}=300\,B_{6}{\cal R}_{\rm rel,1}\,\eta_{1}^{-3/2}(\sigma/1+\sigma)^{3/2}\,\rm TeV$ .

In the case that $\sigma\gg 1$ , magnetic reconnection processes may plausibly occur (Fiorillo et al., 2024). Protons may be accelerated in the magnetic reconnection layer over a timescale $t_{\rm acc}^{\rm rec}\approx E_{p}/\beta_{\rm rec}eBc$ , or

\frac{t_{\rm acc}^{\rm rec}}{t_{\rm cross}}=\frac{r_{L}}{R}\frac{1}{\beta_{\rm rec}}.

(5)

Particles escape the magnetic reconnection region by advecting out of the reconnection layer after $t_{\rm ad}\approx L/c$ , where $L\sim\beta_{\rm rec}R$ is the length of the reconnection layer (Beloborodov, 2017). This yields

\frac{t_{\rm esc}^{\rm rec}}{t_{\rm cross}}=\beta_{\rm rec}.

(6)

Comparing the acceleration and escape time leads to a maximum proton energy $E_{p,\rm max}^{\rm rec}=E_{m}\beta_{\rm rec}^{2}=100\,B_{6}{\cal R}_{\rm rel,1}\,\beta_{\rm rec,-1}^{2}\,\rm TeV$ .

In both turbulence and reconnection scenarios, protons may also escape the acceleration region by being advected to inside the innermost stable orbit. The fallback time is $t_{\rm fall}=R/\alpha v_{K}$ , where $\alpha\sim 0.1$ is the viscous parameter and $v_{K}=(GM/R)^{1/2}$ is the disk velocity in circular Keplerian orbits (Shakura & Sunyaev, 1973), $t_{\rm fall}/t_{\rm cross}\approx 32\,{\cal R}_{\rm rel,1}^{1/2}$ . This timescale is usually longer, and thus less important, than the diffusion or advection time. A relativistic proton may be confined in the acceleration region by

t_{\rm conf}=\min(t_{\rm esc},t_{\rm fall})

(7)

with $t_{\rm esc}$ defined in equations 4 and 6 for turbulent and reconnection scenarios, respectively. Figure 1 presents the acceleration and escape time for a Cygnus X-1-like corona.

3.2 Interaction Timescales

Relativistic protons cool by interacting with the ambient gas and radiation fields. The optical depths for protons due to photopion interaction ( $p\gamma_{b}\rightarrow N\pi$ ) can be estimated as

\tau_{p\gamma}=\frac{t_{\rm conf}}{t_{\rm cross}}\frac{t_{\rm cross}}{t_{p\gamma}}\approx\frac{t_{\rm conf}}{t_{\rm cross}}\frac{\tau_{T}+1}{\Omega}\frac{m_{e}c^{2}}{\epsilon_{X}}\frac{\sigma_{p\gamma}}{\sigma_{T}}\ell

(8)

where $t_{p\gamma}\approx(n_{X}\sigma_{p\gamma}c)^{-1}$ is the $p\gamma$ interaction time, $\sigma_{p\gamma}$ is the inelastic photonpion interaction cross section, and the photon number density is estimated as $n_{X}\sim(\tau_{T}+1)L_{X}/\Omega R^{2}c\epsilon_{X}$ , where $\Omega\sim 4\pi$ if $L_{X}$ is dominated by the corona and $\Omega\sim 2\pi$ is $L_{X}$ is dominated by the disk. The same scaling relation applies to the Bethe-Heitler ( $p\gamma_{b}\rightarrow pe^{+}e^{-}$ ) and pair production ( $\gamma\gamma_{b}\rightarrow e^{+}e^{-}$ ) processes after replacing $\sigma_{p\gamma}$ by the inelastic Bethe-Heitler cross section $\sigma_{\rm BH}$ and pair production cross section $\sigma_{\gamma\gamma}$ .

The proton-proton interaction ( $pp\rightarrow N\pi$ ) time may be written as

\tau_{pp}=\frac{t_{\rm conf}}{t_{\rm cross}}\frac{t_{\rm cross}}{t_{pp}}=\frac{t_{\rm conf}}{t_{\rm cross}}\frac{n_{p}}{n_{e}}\frac{\sigma_{pp}}{\sigma_{T}}\tau_{T}

(9)

where $\sigma_{pp}$ is the inelastic proton-proton cross section.

Equations 7, 8, and 9 show that the optical depth for proton interaction in the corona region only depends on the dimensionless compactness parameters and the Thomson optical depth, thus the conditions for proton interaction and secondary production must be similar in the coronae of XRBs and AGNs, despite of their drastic differences in sizes and timescales.

Figure 4 in Appendix B presents the optical depth to relativistic protons and high-energy $\gamma$ -rays computed using the radiation background of Cygnus X-1.

3.3 Gamma-ray and neutrino spectra

We compute the neutrino and $\gamma$ -ray spectra by iteratively solving a set of transport equations (as listed in Appendix C) that account for the particle escape, energy loss and injection due to proton-proton, photopion, Bethe-Heitler, synchrotron, synchrotron self-absorption, inverse Compton, and photon-photon pair production processes.

For the hard state, we assume $R=100\,R_{g}$ , yielding $\sigma_{\pm}\sim 0.1$ . Stochastic acceleration of protons from a thermal pool or a pre-accelerated population is possible in such a regime. We assume a proton injection spectrum that follows $Q_{p}^{\rm inj}=dN_{p}/dE_{p}\propto E_{p}^{-s}$ from the proton rest mass to $E_{p,\rm max}^{\rm turb}=10$ TeV. We adopt a spectral index $s=2$ though the resulting secondary spectra are similar within the range $s\sim 2-2.3$ . The diffusion of protons takes longer than the light crossing time and efficiently confines the particles for interactions. As a result, neutrinos may be produced in both $pp$ and $p\gamma$ processes.

For the soft state, we assume $R=30\,R_{g}$ . The magnetization parameter $\sigma_{\pm}$ reaches $\sim 60$ , making magnetic reconnection a favored scenario. The proton spectrum may follow a broken power-law with $dN_{p}/dE_{p}\propto E_{p}^{-1}$ below a break energy $E_{p,\rm br}\sim\sigma_{p}m_{p}c^{2}$ and $dN/dE_{p}\propto E_{p}^{-s}$ with $s>2$ (Fiorillo et al., 2024). We thus inject a proton spectrum $Q_{p}^{\rm inj}\propto E_{p}^{-3}$ between $E_{p,\rm min}=60\,\rm GeV$ and $E_{p,\rm max}^{\rm rec}=300$ TeV. Protons leave the reconnection layer rather fast and only $p\gamma$ interaction is efficient in this case.

In both states we set the proton number density to $n_{p}m_{p}=n_{e}m_{e}$ and use the radiation background computed from the observed spectrum as described in Appendix A. The energy density of nonthermal protons is normalized to the X-ray energy density, $u_{p}=u_{X}$ .

Figure 2 compares the resulting $\gamma$ -ray and neutrino spectra to the broad-band spectral energy distribution of Cygnux X-1. In the hard state, GeV and TeV $\gamma$ rays and electrons are produced by $p\gamma$ and $pp$ interactions. Most of them cascade down to 10-100 MeV energies, explaining the MeV $\gamma$ -ray tail observed by COMPTEL. A fraction of the attenuated photons show up in the 1-100 GeV energy range as observed by LAT. In the soft state, fewer interaction happens as a result of the poorer confinement of charged particles in the reconnection layer. Meanwhile, the high compactness makes the accretion flow opaque to $\gamma$ rays above $0.1$ GeV, and photons and pairs from $10-100$ TeV energies cascade down all the way to 10-100 MeV, which explains the observations of COMPTEL and LAT. The $\gamma$ -ray spectrum, comparing with the observations by MAGIC, HAWC, and LHAASO, is presented in Figure 5.

Our model predicts a high-energy neutrino spectrum that cuts off at $\sim 1$ TeV during the hard state and peaks at $\sim 10$ TeV during the soft state. The difference in these cutoff energies is mostly due to the difference in the maximum proton energies of the two states. In the hard state, neutrinos below and above $\sim 50$ GeV mostly come from $pp$ and $p\gamma$ processes, respectively, as suggested by the proton optical depth in Figure 4. In the soft state, neutrinos dominantly come from $p\gamma$ interactions. The neutrino spectra are consistent with the non-detection of IceCube for a point source at Cygnus X-1’s declination but could be observed in the future with combined data samples from different event selections and across experiments.

4 conclusions and discussion

Motivated by the $\gamma$ -ray-opaque neutrino sources in the extragalactic sky, we propose that astroparticles may be produced in the accretion flow of Galactic black hole XRBs. We show that the acceleration and interaction timescales of relativistic protons are tied to the compactness parameter of an X-ray source. Therefore, the same cosmic-ray acceleration and neutrino production processes that happen in an AGN core may also occur in our local accreting black holes.

We have assumed a spherical geometry for the coronal region. While the actual structure of a corona is highly unknown, its shape has been modeled as a slab, wedge, or cone in addition to a sphere. Recent polarimetric observations of Cygnus X-1 start to reveal the geometric structure of the corona and the base of the jet (Krawczynski et al., 2022; Chattopadhyay et al., 2024). Moreover, we have assumed that the magnetic field strength, densities of the radiation field and the plasma are all constant within the coronal region. A more realistic picture would include a density gradient and a radially dependent magnetic field, though it would also introduce more free parameters to the model of particle acceleration and interaction. For example, we have verified that a similar solution like the one in Figure 2 may be found for a slab-like coronal region. We stress that the key results of this work, which only rely on the similarities in the compactness of AGN and XRB coronae, as well as the imprints of the accretion state transitions in the astroparticle production, remain robust under these uncertainties.

Coronal protons may significantly contribute to the Galactic cosmic-ray population. The total X-ray luminosity of Galactic black hole XRBs is at the level of $L_{X}\sim(1-11)\times 10^{40}\,\rm erg\,s^{-1}$ (Grimm et al., 2002; Heinz & Grimm, 2005; Cooper et al., 2020). The coronae of these XRBs provide a total cosmic-ray power $L_{\rm CR,XRB}\lesssim L_{X}$ . This may explain the Galactic cosmic-ray spectrum, which requires $L_{\rm CR,obs}=u_{\rm CR}V_{\rm gal}/\tau_{\rm CR}\sim 5\times 10^{40}\,\rm erg\,s^{-1}$ when assuming cosmic-ray energy density $u_{\rm CR}\approx 1\,\rm eV\,cm^{-3}$ above $1\,\rm GeV$ , escape time $\tau_{\rm CR}\sim 3\,\rm Myr$ (Blasi, 2013) and volume of the Galactic disk $V_{\rm gal}\approx\pi R_{\rm gal}^{2}d_{\rm gal}$ with $R_{\rm gal}\sim 15\,\rm kpc$ and $d_{\rm gal}\sim 0.2\,\rm kpc$ . Unlike Cygnus X-1 which spends most of its time in the hard state due to a high accretion rate, black hole XRBs are not commonly found to stay in the hard state. Considering some average of the soft and hard state, a cosmic-ray injection spectral index $s\sim 2.2-2.5$ could be expected, which would explain the observed index $s_{\rm obs}\sim 2.7$ after the cosmic rays diffuse in the Galactic magnetic field.

The optical depth for the relativistic protons is at the level of $\tau_{p\gamma}\sim 10^{-2}$ around $E_{p}\sim 100$ TeV, and depends on the size of the corona and confinement time at lower proton energy. Taking $\tau_{p\gamma}\sim 10^{-2}$ , the total neutrino luminosity emitted by the XRB coronae is $L_{\nu,\rm XRB}\sim 3/8L_{CR,\rm XRB}\tau_{p\gamma}\sim 6\times 10^{36}\,\rm erg\,s^{-1}$ above $1$ TeV. This corresponds to an all-flavor neutrino flux for an observer in the solar neighborhood (Fang & Halzen, 2024) at the level of $E_{\nu}^{2}F_{\nu}\approx 3L_{\nu,\rm XRB}/4\pi r_{\odot}^{2}=1.3\times 10^{-6}\,\rm GeV\,cm^{-2}s^{-1}$ at 1 TeV, which is comparable to the Galactic plane flux measured by IceCube (IceCube Collaboration, 2023).

The work of K.F and F.H is supported by the Office of the Vice Chancellor for Research at the University of Wisconsin–Madison with funding from the Wisconsin Alumni Research Foundation. K.F. acknowledges support from the National Science Foundation (PHY-2110821, PHY-2238916) and NASA (NMH211ZDA001N-Fermi). This work was supported by a grant from the Simons Foundation (00001470, KF). K.F acknowledges the support of the Sloan Research Fellowship. The research of F.H was also supported in part by the U.S. National Science Foundation under grants PHY-2209445 and OPP-2042807. J.S.G. thanks the University of Wisconsin College of Letters and Science for partial support of his IceCube-related research. SH acknowledges funding from NASA Astrophysics Theory Grant NNX17AJ98G.

References

Aartsen et al. (2020) Aartsen, M. G., et al. 2020, Phys. Rev. Lett., 124, 051103, doi: 10.1103/PhysRevLett.124.051103
Abbasi et al. (2022a) Abbasi, R., et al. 2022a, Science, 378, 538, doi: 10.1126/science.abg3395
Abbasi et al. (2022b) —. 2022b, Astrophys. J. Lett., 930, L24, doi: 10.3847/2041-8213/ac67d8
Abdollahi et al. (2020) Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33, doi: 10.3847/1538-4365/ab6bcb
Ackermann et al. (2015) Ackermann, M., Ajello, M., Albert, A., et al. 2015, ApJ, 799, 86, doi: 10.1088/0004-637X/799/1/86
Ackermann et al. (2012) Ackermann, M., et al. 2012, Astrophys. J., 750, 3, doi: 10.1088/0004-637X/750/1/3
Aharonian et al. (2005) Aharonian, F., et al. 2005, Science, 309, 746, doi: 10.1126/science.1113764
Aharonian et al. (2006) —. 2006, Astron. Astrophys., 460, 743, doi: 10.1051/0004-6361:20065940
Albert et al. (2021) Albert, A., et al. 2021, Astrophys. J. Lett., 912, L4, doi: 10.3847/2041-8213/abf35a
Albert et al. (2021) Albert, A., Alfaro, R., Alvarez, C., et al. 2021, ApJ, 912, L4, doi: 10.3847/2041-8213/abf35a
Bauer et al. (2015) Bauer, F. E., Arévalo, P., Walton, D. J., et al. 2015, ApJ, 812, 116, doi: 10.1088/0004-637X/812/2/116
Becker et al. (2006) Becker, P. A., Le, T., & Dermer, C. D. 2006, Astrophys. J., 647, 539, doi: 10.1086/505319
Bednarek (2005) Bednarek, W. 2005, Astrophys. J., 631, 466, doi: 10.1086/432411
Beloborodov (2017) Beloborodov, A. M. 2017, Astrophys. J., 850, 141, doi: 10.3847/1538-4357/aa8f4f
Berezinskii & Ginzburg (1981) Berezinskii, V. S., & Ginzburg, V. L. 1981, Mon. Not. Roy. Astron. Soc., 194, 3, doi: 10.1093/mnras/194.1.3
Blasi (2013) Blasi, P. 2013, Astron. Astrophys. Rev., 21, 70, doi: 10.1007/s00159-013-0070-7
BLUMENTHAL & GOULD (1970) BLUMENTHAL, G. R., & GOULD, R. J. 1970, Rev. Mod. Phys., 42, 237, doi: 10.1103/RevModPhys.42.237
Boettcher et al. (2012) Boettcher, M., Harris, D. E., & Krawczynski, H. 2012, Relativistic Jets from Active Galactic Nuclei
Bottcher & Schlickeiser (1997) Bottcher, M., & Schlickeiser, R. 1997, Astron. Astrophys., 325, 866. https://arxiv.org/abs/astro-ph/9703069
Cangemi et al. (2021) Cangemi, F., et al. 2021, Astron. Astrophys., 650, A93, doi: 10.1051/0004-6361/202038604
Cao et al. (2024) Cao, Z., et al. 2024. https://arxiv.org/abs/2410.08988
Chattopadhyay et al. (2024) Chattopadhyay, T., Kumar, A., Rao, A. R., et al. 2024, ApJ, 960, L2, doi: 10.3847/2041-8213/ad118d
Chodorowski et al. (1992) Chodorowski, M. J., Zdziarski, A. A., & Sikora, M. 1992, ApJ, 400, 181, doi: 10.1086/171984
Cooper et al. (2020) Cooper, A. J., Gaggero, D., Markoff, S., & Zhang, S. 2020, MNRAS, 493, 3212, doi: 10.1093/mnras/staa373
Coppi & Blandford (1990) Coppi, P. S., & Blandford, R. D. 1990, MNRAS, 245, 453, doi: 10.1093/mnras/245.3.453
Del Santo et al. (2013) Del Santo, M., Malzac, J., Belmont, R., Bouchet, L., & De Cesare, G. 2013, MNRAS, 430, 209, doi: 10.1093/mnras/sts574
Dermer & Menon (2009) Dermer, C. D., & Menon, G. 2009, High Energy Radiation from Black Holes: Gamma Rays, Cosmic Rays, and Neutrinos
Dexter et al. (2014) Dexter, J., McKinney, J. C., Markoff, S., & Tchekhovskoy, A. 2014, MNRAS, 440, 2185, doi: 10.1093/mnras/stu581
Di Salvo et al. (2001) Di Salvo, T., Done, C., Zycki, P. T., Burderi, L., & Robba, N. R. 2001, Astrophys. J., 547, 1024, doi: 10.1086/318396
Done et al. (2007) Done, C., Gierlinski, M., & Kubota, A. 2007, Astron. Astrophys. Rev., 15, 1, doi: 10.1007/s00159-007-0006-1
Dubus (2013) Dubus, G. 2013, A&A Rev., 21, 64, doi: 10.1007/s00159-013-0064-5
Effenberger & Petrosian (2018) Effenberger, F., & Petrosian, V. 2018, ApJ, 868, L28, doi: 10.3847/2041-8213/aaedb3
Fabian et al. (2015) Fabian, A. C., Lohfink, A., Kara, E., et al. 2015, MNRAS, 451, 4375, doi: 10.1093/mnras/stv1218
Falcke et al. (2004) Falcke, H., Körding, E., & Markoff, S. 2004, A&A, 414, 895, doi: 10.1051/0004-6361:20031683
Fang et al. (2022) Fang, K., Gallagher, J. S., & Halzen, F. 2022, Astrophys. J., 933, 190, doi: 10.3847/1538-4357/ac7649
Fang & Halzen (2024) Fang, K., & Halzen, F. 2024, arXiv e-prints, arXiv:2404.15944, doi: 10.48550/arXiv.2404.15944
Fang & Murase (2023) Fang, K., & Murase, K. 2023, Astrophys. J. Lett., 957, L6, doi: 10.3847/2041-8213/ad012f
Fender et al. (2006) Fender, R. P., Stirling, A. M., Spencer, R. E., et al. 2006, MNRAS, 369, 603, doi: 10.1111/j.1365-2966.2006.10193.x
Fermi-LAT Collaboration (2009) Fermi-LAT Collaboration. 2009, Science, 326, 1512. http://www.jstor.org/stable/27736637
Fernández-Barral et al. (2018) Fernández-Barral, A., et al. 2018, PoS, ICRC2017, 734, doi: 10.22323/1.301.0734
Fiorillo et al. (2024) Fiorillo, D. F. G., Petropoulou, M., Comisso, L., Peretti, E., & Sironi, L. 2024, Astrophys. J., 961, L14, doi: 10.3847/2041-8213/ad192b
Frontera et al. (2001) Frontera, F., et al. 2001, Astrophys. J., 546, 1027, doi: 10.1086/318304
Grimm et al. (2002) Grimm, H. J., Gilfanov, M., & Sunyaev, R. 2002, A&A, 391, 923, doi: 10.1051/0004-6361:20020826
Grinberg et al. (2013) Grinberg, V., et al. 2013, Astron. Astrophys., 554, A88, doi: 10.1051/0004-6361/201321128
Guilbert et al. (1983) Guilbert, P. W., Fabian, A. C., & Rees, M. J. 1983, MNRAS, 205, 593, doi: 10.1093/mnras/205.3.593
Heinz & Grimm (2005) Heinz, S., & Grimm, H. J. 2005, Astrophys. J., 633, 384, doi: 10.1086/452624
IceCube Collaboration (2023) IceCube Collaboration. 2023, Science, 380, 1338, doi: 10.1126/science.adc9818
Inoue et al. (2020) Inoue, Y., Khangulyan, D., & Doi, A. 2020, ApJ, 891, L33, doi: 10.3847/2041-8213/ab7661
Inoue et al. (2021) —. 2021, Galaxies, 9, 36, doi: 10.3390/galaxies9020036
Kantzas et al. (2023) Kantzas, D., Markoff, S., Cooper, A. J., et al. 2023, MNRAS, 524, 1326, doi: 10.1093/mnras/stad1909
Kelner & Aharonian (2008) Kelner, S. R., & Aharonian, F. A. 2008, Phys. Rev. D, 78, 034013, doi: 10.1103/PhysRevD.82.099901
Koldobskiy et al. (2021) Koldobskiy, S., Kachelrieß, M., Lskavyan, A., et al. 2021, Phys. Rev. D, 104, 123027, doi: 10.1103/PhysRevD.104.123027
Krawczynski et al. (2022) Krawczynski, H., et al. 2022, Science, 378, add5399, doi: 10.1126/science.add5399
Malyshev et al. (2013) Malyshev, D., Zdziarski, A. A., & Chernyakova, M. 2013, MNRAS, 434, 2380, doi: 10.1093/mnras/stt1184
Mbarek et al. (2023) Mbarek, R., Philippov, A., Chernoglazov, A., Levinson, A., & Mushotzky, R. 2023, arXiv e-prints, arXiv:2310.15222, doi: 10.48550/arXiv.2310.15222
McConnell et al. (2002) McConnell, M. L., et al. 2002, Astrophys. J., 572, 984, doi: 10.1086/340436
Merloni et al. (2003) Merloni, A., Heinz, S., & di Matteo, T. 2003, MNRAS, 345, 1057, doi: 10.1046/j.1365-2966.2003.07017.x
Miller-Jones et al. (2021) Miller-Jones, J. C. A., Bahramian, A., Orosz, J. A., et al. 2021, Science, 371, 1046, doi: 10.1126/science.abb3363
Minezaki & Matsushita (2015) Minezaki, T., & Matsushita, K. 2015, ApJ, 802, 98, doi: 10.1088/0004-637X/802/2/98
Mirabel et al. (1996) Mirabel, I. F., Claret, A., Cesarsky, C. J., Boulade, O., & Cesarsky, D. A. 1996, A&A, 315, L113
Murase (2022) Murase, K. 2022, ApJ, 941, L17, doi: 10.3847/2041-8213/aca53c
Murase et al. (2016) Murase, K., Guetta, D., & Ahlers, M. 2016, Phys. Rev. Lett., 116, 071101, doi: 10.1103/PhysRevLett.116.071101
Murase et al. (2020) Murase, K., Kimura, S. S., & Mészáros, P. 2020, Phys. Rev. Lett., 125, 011101, doi: 10.1103/PhysRevLett.125.011101
Papavasileiou et al. (2023) Papavasileiou, T. V., Kosmas, O. T., & Sinatkas, I. 2023, A&A, 673, A162, doi: 10.1051/0004-6361/202345869
Pepe et al. (2015) Pepe, C., Vila, G. S., & Romero, G. E. 2015, Astron. Astrophys., 584, A95, doi: 10.1051/0004-6361/201527156
Persi et al. (1980) Persi, P., Ferrari-Toniolo, M., Grasdalen, G. L., & Spada, G. 1980, A&A, 92, 238
Poutanen & Vurm (2009) Poutanen, J., & Vurm, I. 2009, Astrophys. J. Lett., 690, L97, doi: 10.1088/0004-637X/690/2/L97
Romero et al. (2003) Romero, G. E., Torres, D. F., Bernado, M. M. K., & Mirabel, I. F. 2003, Astron. Astrophys., 410, L1, doi: 10.1051/0004-6361:20031314-1
Romero et al. (2014) Romero, G. E., Vieyro, F. L., & Chaty, S. 2014, Astron. Astrophys., 562, L7, doi: 10.1051/0004-6361/201323316
Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
Sironi et al. (2015) Sironi, L., Petropoulou, M., & Giannios, D. 2015, Mon. Not. Roy. Astron. Soc., 450, 183, doi: 10.1093/mnras/stv641
Stathopoulos et al. (2024) Stathopoulos, S. I., Petropoulou, M., Vasilopoulos, G., & Mastichiadis, A. 2024, Astron. Astrophys., 683, A225, doi: 10.1051/0004-6361/202347277
Stawarz & Petrosian (2008) Stawarz, L., & Petrosian, V. 2008, Astrophys. J., 681, 1725, doi: 10.1086/588813
Stecker et al. (1991) Stecker, F. W., Done, C., Salamon, M. H., & Sommers, P. 1991, Phys. Rev. Lett., 66, 2697, doi: 10.1103/PhysRevLett.66.2697
Stecker et al. (1992) —. 1992, Phys. Rev. Lett., 69, 2738, doi: 10.1103/PhysRevLett.69.2738
Stecker et al. (1985) Stecker, F. W., Harding, A. K., & Barnard, J. J. 1985, Nature, 316, 418, doi: 10.1038/316418a0
Stirling et al. (2001) Stirling, A. M., Spencer, R. E., de la Force, C. J., et al. 2001, MNRAS, 327, 1273, doi: 10.1046/j.1365-8711.2001.04821.x
Sudoh & Beacom (2023a) Sudoh, T., & Beacom, J. F. 2023a, Phys. Rev. D, 107, 043002, doi: 10.1103/PhysRevD.107.043002
Sudoh & Beacom (2023b) —. 2023b, Phys. Rev. D, 108, 043016, doi: 10.1103/PhysRevD.108.043016
Tibet AS ${\gamma}$ Collaboration et al. (2021) Tibet AS ${\gamma}$ Collaboration, Amenomori, M., Bao, Y. W., et al. 2021, Phys. Rev. Lett., 126, 141101, doi: 10.1103/PhysRevLett.126.141101
Yuan & Narayan (2014) Yuan, F., & Narayan, R. 2014, Ann. Rev. Astron. Astrophys., 52, 529, doi: 10.1146/annurev-astro-082812-141003
Zdziarski et al. (2017) Zdziarski, A. A., Malyshev, D., Chernyakova, M., & Pooley, G. G. 2017, MNRAS, 471, 3657, doi: 10.1093/mnras/stx1846
Zdziarski et al. (2014) Zdziarski, A. A., Pjanka, P., Sikora, M., & Stawarz, Ł. 2014, MNRAS, 442, 3243, doi: 10.1093/mnras/stu1009

Appendix A Radiation Background

We model the radiation field of Cygnus X-1 from infrared to MeV as follows.

Following McConnell et al. (2002), we model the emission by the corona and disk region as three components, including i) blackbody radiation from a Shakura–Sunyaev thin disk. We compute the emission using agnpy with a disk luminosity $L_{\rm d}$ . The parameters of the inner and outer radii are selected to match the positions of the soft X-ray peak. ii) Inverse-Compton radiation by thermal electrons in the accretion flow, modeled as a power-law spectrum $dn_{\rm c,th}/d\epsilon\propto\epsilon^{-s_{\rm th}}\exp(-\epsilon/\epsilon_{\rm th,max})$ above a minimum energy $\epsilon_{\rm th,min}$ . iii) Inverse-Compton radiation by non-thermal electrons in the accretion flow, modeled as a second power-law spectrum, $dn_{\rm c,nth}/d\epsilon\propto\epsilon^{-s_{\rm nth}}\exp(-\epsilon/\epsilon_{\rm nth,max})$ starting from the maximum energy of the thermal component $\epsilon_{\rm th,max}$ . The third component is only included for the hard state. The geometry of the corona is modeled as a sphere, subject to soft photon flux from a source from the disk subtending a solid angle $\Omega\sim 2\pi$ and a source inside the corona subtending $\Omega\sim 4\pi$ .

Figure 3 compares our model to the X-ray and MeV gamma-ray spectral energy distribution of Cygnus X-1. We list the model parameters in Table 1. They are in agreement with the parameters obtained by fitting more sophisticated models to the X-ray data in the hard state and up to a factor of $\sim 2$ in the soft state (McConnell et al., 2002). The difference is mostly due to the much simpler model we use. For example, our model does not account for the Compton reflection. Moreover, we model the hard X-ray and MeV gamma-ray emission as single power-laws, while McConnell et al. (2002) computes the inverse Compton emission using ambient seed photons. Nonetheless, the details of the X-ray modeling barely matter here, as the high-energy neutrino and gamma-ray production is only sensitive to the energy density and spectral shape of the X-ray photons.

We model the stellar radiation as a blackbody spectrum with temperature $T_{*}=2.55\times 10^{4}$ K (Zdziarski et al., 2017) from a companion star with orbital separation $a_{*}=3.2\times 10^{12}\,\rm cm$ , with a radius $R_{*}=19\,\rm R_{\odot}$ (Pepe et al., 2015). However, due to the large distance between the companion star and the black hole, the energy density of the stellar radiation at the corona is negligible compared to the emission by the disk and the synchrotron radiation by secondary electrons.

Appendix B Optical Depth

Using the radiation fields described in the last section, we compute the optical depth to relativistic protons and high-energy $\gamma$ -rays due to various interaction processes. Figure 4 shows that in both states, the corona region is optically thick to $\gamma$ -rays above $\sim 10^{-1}$ GeV. The pionic $\gamma$ -rays are attenuated by the infrared through X-ray photons from the disk and cascade down to lower energies.

Appendix C Particle Transport Equations

We compute the neutrino and $\gamma$ -ray spectra by iteratively solving the following coupled transport equations.

$\displaystyle\frac{\partial N_{p}}{\partial t}+\frac{N_{p}}{\tau_{p,\rm esc}}+\frac{\partial}{\partial\gamma_{p}}\left[\left(\frac{d\gamma_{p}}{dt}\middle\|_{\rm syn}+\frac{d\gamma_{p}}{dt}\middle\|_{\rm BH}++\frac{d\gamma_{p}}{dt}\middle\|_{pp}+\frac{d\gamma_{p}}{dt}\middle\|_{p\gamma}\right)N_{p}\right]$	$\displaystyle=$	$\displaystyle Q_{p}^{\rm inj}$	(C1)
$\displaystyle\frac{\partial N_{e}}{\partial t}+\frac{N_{e}}{\tau_{e,\rm esc}}+\frac{\partial}{\partial\gamma_{e}}\left[\left(\frac{d\gamma_{e}}{dt}\middle\|_{\rm syn}+\frac{d\gamma_{e}}{dt}\middle\|_{\rm IC}\right)N_{e}\right]$	$\displaystyle=$	$\displaystyle Q_{e}^{p\gamma}+Q_{e}^{\rm BH}+Q_{e}^{pp}+Q_{e}^{\gamma\gamma}$	(C2)
$\displaystyle\frac{\partial N_{\gamma}}{\partial t}+\frac{N_{\gamma}}{\tau_{\gamma,\rm esc}}+N_{\gamma}c\left(a_{\gamma\gamma}+a_{\rm ssa}\right)$	$\displaystyle=$	$\displaystyle Q_{\gamma}^{p\gamma}+Q_{\gamma}^{pp}+Q_{\gamma}^{\rm IC}+Q_{\gamma}^{\rm e,syn}+Q_{\gamma}^{\rm p,syn}$	(C3)
$\displaystyle\frac{\partial N_{\nu}}{\partial t}+\frac{N_{\nu}}{\tau_{\nu,\rm esc}}$	$\displaystyle=$	$\displaystyle Q_{\nu}^{p\gamma}+Q_{\nu}^{pp}$	(C4)

Here $N_{i}\equiv(dN/dE)_{i}$ is the spectrum of a particle species $i=p,e,\gamma,\nu$ . For charged particles, the escape time $\tau_{\rm esc}$ is set to $t_{\rm conf}$ (equation 7) for a given proton or electron energy. For neutral particles, $\tau_{\rm esc}=t_{\rm cross}$ . The energy loss terms $d\gamma_{i}/dt|_{\rm process}$ and the injection terms $Q^{\rm process}_{j}$ denote the loss rate of particle species $i$ and the production of secondary particle $j$ due to one of the following process: synchrotron radiation (“syn”; Boettcher et al., 2012), Bethe-Heitler process (“BH”; Chodorowski et al., 1992), proton-proton interaction ( $pp$ ; Koldobskiy et al., 2021), photonpion production ( $p\gamma$ ; Kelner & Aharonian, 2008), inverse Compton radiation (“IC”; BLUMENTHAL & GOULD, 1970), and pair production (also called as pair annihilation, $\gamma\gamma$ ; Bottcher & Schlickeiser, 1997). For photons, the energy loss rates due to pair production and synchrotron self-absorption are computed using the efficiencies $a_{\gamma\gamma}$ (Coppi & Blandford, 1990) and $a_{\rm ssa}$ (Dermer & Menon, 2009). The equations are solved using a fully implicit finite difference scheme as in Stathopoulos et al. (2024) and evolved over a time period long enough to reach a steady-state solution.

The proton injection term $Q_{p}^{\rm inj}$ is assumed to be a single power-law from either the proton rest mass, in the turbulence scenario, or a break energy $\sim\sigma_{p}m_{p}c^{2}$ , in the reconnection scenario, to the maximum energy. The actual nonthermal proton spectrum could be more complicated than a power-law for several reasons. We have neglected the momentum diffusion terms in Equation C1 as the acceleration time is much shorter than the other timescales, though including such terms is needed to capture the broadening of the particle distribution due to momentum space diffusion (Becker et al., 2006; Stawarz & Petrosian, 2008). Moreover, the acceleration process could be more complicated than what has been assumed. For example, Mbarek et al. (2023) proposes an alternative acceleration picture to the turbulence and reconnection scenarios. They suggest that protons are injected into the corona with a hard spectrum and get re-accelerated to 1-10 TeV energies. Such a two-step acceleration is needed to explain the neutrino flux of NGC 1068, as neither relativistic magnetic reconnection nor turbulence seems to accelerate sufficient relativistic protons if their initial energy follows a thermal distribution.

Appendix D Gamma-ray spectrum

Figure 5 zooms into the spectral energy distribution of Cygnus X-1 in the gamma-ray energy range and compares our model to observations. The black curve indicates a state-averaged spectrum considering that on average the source spends $\sim 80\%$ and $\sim 20\%$ in hard and soft states, respectively (Grinberg et al., 2013). Our model explains the LHAASO observation above 20 TeV. We note that the emission between $\sim$ 10 GeV and $\sim$ 10 TeV would be further absorbed by the stellar photons (Albert et al., 2021), resulting in periodic modulation.

Table 1: Parameters of the X-ray emission model, assuming a black hole mass of

21.2\,M_{\odot}

and source distance

2.22\,\rm kpc

{ruledtabular}

State	$L_{\rm c,th}$	$s_{\rm th}$	$\epsilon_{\rm th,min}$	$\epsilon_{\rm th,max}$	$L_{\rm c,nth}$	$s_{\rm nth}$	$\epsilon_{\rm nth,max}$	$L_{\rm d}$
	$[\rm erg\,s^{-1}]$		[keV]	[keV]	$[\rm erg\,s^{-1}]$		[keV]	$[\rm erg\,s^{-1}]$
Hard	$2\times 10^{37}$	1.3	$0.3$	$150$	$5\times 10^{35}$	$2.2$	$2\times 10^{4}$	$4\times 10^{36}$
Soft	$1.2\times 10^{37}$	2.4	$2$	$2\times 10^{4}$	-	-	-	$1.2\times 10^{38}$