Prepared for submission to JCAP

Neutrino and pair creation in

reconnection-powered coronae of
accreting black holes
arXiv:2410.12638v2 [astro-ph.HE] 17 Oct 2024

D. Karavola,1 M. Petropoulou,1,2 D. F. G. Fiorillo,3 L. Comisso,4

L. Sironi4,5
1 National and Kapodistrian University of Athens, University Campus Zografos, GR 15784,
Athens, Greece
2 Institute of Accelerating Systems & Applications, University Campus Zografos, GR 15784,

Athens, Greece
3 Deutsches Elektronen-Synchrotron DESY, Platanenallee 6, 15738 Zeuthen, Germany
4 Department of Astronomy and Columbia Astrophysics Laboratory, Columbia University,

New York, NY 10027, USA

5 Center for Computational Astrophysics, Flatiron Institute, New York, NY, 10010, USA

E-mail: [email protected], [email protected], [email protected],

[email protected], [email protected]

Abstract. A ubiquitous feature of accreting black hole systems is their hard X-ray emission
which is thought to be produced through Comptonization of soft photons by electrons and
positrons in the vicinity of the black hole, in a region with optical depth of order unity. The
origin and composition of this Comptonizing region, known as the corona, is a matter open for
debate. In this paper we investigate the role of relativistic protons accelerated in black-hole
magnetospheric current sheets for the pair enrichment and neutrino emission of AGN coronae.
Our model has two free parameters, namely the proton plasma magnetization σp , which
controls the peak energy of the neutrino spectrum, and the Eddington ratio λEdd (defined
as the ratio between X-ray luminosity LX and Eddington luminosity LEdd ), which controls
the amount of energy transferred to secondary particles. For sources with λEdd ≳ λEdd,crit
(where λEdd,crit ∼ 10−1 for σp = 105 or ∼ 10−2 for σp = 107 ), proton-photon interactions
and γγ annihilation produce enough secondary pairs to achieve Thomson optical depths
τT ∼ 0.1 − 10. Additionally, we find that the neutrino luminosity scales as L2X /LEdd for
λEdd ≲ λEdd,crit , while it is proportional to LX for higher λEdd values. We apply our model to
four Seyfert galaxies, including NGC 1068, and discuss our results in light of recent IceCube
observations.
Contents

1 Introduction 1

2 Model 2

3 Analytical Estimates 5
3.1 Neutrino luminosity and peak energy 5
3.2 Calorimetric limit for pγ interactions 6
3.3 Secondary pair number density 6
3.3.1 The γγ annihilation channel 7

4 Numerical results 8

5 Neutrino predictions for Seyfert galaxies 15

6 Summary and Discussion 18

A Analytical expression for Lν+ν̄ 24

B Other channels for pair production 25

B.1 The photomeson (pγ) production channel 25
B.2 Proton-photon (Bethe-Heitler) pair production 25

C Proton distribution with a post-break slope of 2 26

D γ − γ optical depth 27

1 Introduction

Active Galactic Nuclei (AGN) are the most powerful steady emitters of non-thermal radiation.
They are likely powered by accretion onto a central supermassive black hole. The accretion
disk feeding the black hole is responsible for the emission of optical-ultraviolet (OUV) radi-
ation, while at higher energies from keV to hundreds of keV, in the X-ray band, usually a
power-law emission is observed; a typical AGN spectrum is shown, e.g., in figure 1 of Ref. [1].
The origin of the X-rays is traditionally attributed to a compact region close to the black
hole, the so-called corona, where a population of hot electrons can Comptonize soft OUV
photons, leading to the observed power law [2, 3]. From the point of view of energetics, this
moves the question to what is the mechanism that keeps heating the electrons in the corona,
as they would otherwise cool due to strong radiative losses.
One possible mechanism for maintaining the high temperature of electrons in the corona
is magnetic reconnection, which has also been proposed for stellar-mass accreting black holes
(e.g. [4, 5]). Recently, Ref. [6] pointed out that if current sheets form in accreting systems with
high density of seed photons, inverse Compton cooling would keep the electron population
at low temperatures (kTe ≪ 100 keV). These “cold” electrons are found in plasmoids that
move with mildly relativistic speeds and are formed in current sheets. Interestingly, radiative
particle-in-cell (PIC) simulations have shown that the bulk motions of cold plasma in the

–1–
current sheet can mimic a quasi-thermal electron distribution with an effective temperature
close to 100 keV [7–9]. Additionally, a large fraction of the dissipated magnetic energy is
converted into X-ray radiation. In this context, the magnetic field strength in the coronal
region can be inferred by the observed X-ray flux.
In the last few years, the question of coronal AGN emission has taken on a new rel-
evance, after the detection of high-energy neutrinos from NGC 1068 [10], a Seyfert galaxy
at ≃ 10.1 Mpc. The integrated all-flavor neutrino luminosity is Lν+ν̄ ≃ 1042 erg/s, while
a comparable TeV γ-ray emission is not observed [11]. This suggests as a site for neutrino
emission a compact region where γ-rays would be attenuated, so that the corona emerges as a
natural candidate [12]. This raises new questions on the nature of the non-thermal processes
powering the corona, since not only leptons but also hadrons need to be energized to power
neutrino emission.
A scenario which incorporates magnetic reconnection powering hadron acceleration was
recently proposed by some of us in Ref. [13]. In this scenario, a large-scale reconnection layer
forming in the black hole magnetospheric region (see sketch in figure 1) accelerates protons
to energies of tens of TeV. The subsequent pγ interactions with the coronal X-ray photons
trigger an electromagnetic cascade. The latter results in a production of pairs maintaining
optical thickness (Thomson optical depth of the corona ≳1), while also powering a neutrino
signal consistent with the IceCube observations. While this study was motivated by the
IceCube results of NGC 1068, and as such it was tailored to this source, it remains an open
question to understand what features are likely generic to many AGN coronae, and how do
the expectations for the neutrino signal vary among different sources.
In this work, we address this question, by conducting a comprehensive study of the
neutrino emission and pair density expected from the coronal region of AGN in the pres-
ence of magnetospheric current sheets. The predictions for the electromagnetic and neutrino
signal depend on only three parameters, namely the black hole mass, the Eddington ratio
of the accreting system (defined using the X-ray luminosity of the corona), and the proton
plasma magnetization. We perform numerical calculations of the electromagnetic cascade
and hadronic neutrino emission, comparing the results with the analytical expectations, to
provide the scalings of the expected signal with these three parameters. In this way, we are
able to make predictions on the expected neutrino signal from other AGN cores similar to
NGC 1068.
This paper is structured as follows. In Section 2 we present a detailed description of our
coronal modelf, and in Section 3 we provide analytical expressions for the predicted neutrino
luminosity and secondary pair density. In Section 4 we present the numerical results of our
parameter study using the code ATHEνA. In Section 5 we test our model on individual
Seyfert galaxies and compare the numerical results for the neutrino spectra to those observed
by IceCube. We conclude in Section 6 with a summary of our results and a discussion.

2 Model

We set up a model describing the coronal environment in which the X-ray emission of the
corona is powered by magnetic dissipation at a large-scale magnetic reconnection layer (with
typical size R of a few gravitational radii). Motivated by theory and simulations [6, 7], we
express the broadband X-ray luminosity of the corona LX as
c
LX = ηX SP A = ηX Erec BA (2.1)
4π

–2–
Figure 1: AGN sketch showing the central massive black hole surrounded by an accretion
disk. Perpendicularly to the disk a relativistic jet can be formed. A corona can be shown
close to the central object. The numbers in the sketch refer to: (1) magnetospheric current
sheet, (2) circumnuclear medium, (3) turbulent accretion disk (4) jet and (5) jet boundary
current sheet.

where ηX ∼ 0.5 is the fraction of the magnetic energy dissipated to the X-ray photon field,
SP is the Poynting energy flux, A is the surface of the reconnection layer (taken to be 4πR2 ),
B and Erec = βrec B are respectively the magnetic field and electric field strengths in the
upstream region of the layer, with βrec ∼ 0.1 being the reconnection speed in units of the
speed of light c [14–16]. Then, the magnetic field can be expressed as
 1
LX 2
B= . (2.2)
ηX cβrec R2
We assume that X-ray photons are produced via Comptonization of low-energy seed photons
(e.g. from the accretion disk) by electron-positron pairs in the corona that do not necessarily
have a hadronic origin. We do not attempt to model in detail the spectral formation of the
X-ray coronal emission, but we discuss potential improvements in this regard in Section 6.
Instead, we adopt an observationally motivated X-ray photon spectrum, in this regard a power
law with photon index sX ∼ 2 extending from EX,min = 0.1 keV to EX,max = 100 keV,

nX (E) = nX,0 E −sX , EX,min ≤ E ≤ EX,max , (2.3)

where nX (E) is the differential number density of photons in the corona. The X-ray spectra
of AGN suggest that the optical depth for Comptonization is moderate (i.e., τT ∼ 0.1 − 10,
[3, 17, 18]). We can therefore estimate the pair density that is needed to produce the assumed
X-ray spectrum, as n± = τT /(σT R), where σT ≃ 6.65·10−25 cm2 is the Thomson cross section.
The magnetization of the pair plasma is then computed as
B2 σT LX ℓX
σe = 2
= 3
= , (2.4)
4πn± me c 4πηX βrec τT Rme c ηX βrec τT
where ℓX = σT LX /4πRme c3 is the X-ray compactness of the corona. The latter can also be
expressed in terms of the Eddington luminosity of an accreting black hole with mass Mbh ,

–3–
LEdd ≃ 1.3 · 1045 Mbh /(107 M⊙ ) erg s−1 , as:
LX 1 mp 1 mp
ℓX = = λEdd , (2.5)
LEdd R̃ me R̃ me
where R̃ = R/Rg , Rg = GMbh /c2 is the black-hole gravitational radius, and λEdd is the
Eddington ratio (defined using the broadband X-ray luminosity and not the bolometric lumi-
nosity as commonly done in AGN studies). The expected pair magnetization of the coronal
plasma is estimated as
−1 −1 −1
σe ≈ 3700 λEdd,−2 τT,−1 ηX,−0.3 βrec,−1 R̃−1 (2.6)
where we introduced the notation qx = q/10x .
Driven by PIC simulation results of relativistic magnetic reconnection (σe ≫ 1) we
expect that any protons, if present in the upstream region of the layer, will be injected to the
acceleration process, and energize at a rate dEp /dt = eβrec Bc, where Ep = γp mp c2 and γp
is the proton Lorentz factor [14] (for a qualitative discussion see also [15, 16]). We consider
that a fraction ηp of the Poynting luminosity is transferred to non-thermal protons. Similarly
to Eq. 2.1, the bolometric energy injection rate into relativistic protons can be written as
ηp
Lp = ηp SP A = LX . (2.7)
ηX
We consider that protons are accelerated in a prior phase in the upstream region, before
injected in the coronal environment. As a result, we assume that the differential spectrum of
protons injected into the corona per unit time by the acceleration process is also described
by a broken power law [14–16],
(
−s +1
d2 Np γp,brp γp−1 , 1 ≤ γp ≤ γp,br
= Ṅp,0 −s (2.8)
dγp dt γp p , γp,br < γp ≤ γmax .

where Ṅp,0 is a normalization factor, γp,br is the break Lorentz factor, and sp is the post-
break power-law slope, which is taken to be 2 or 3 (the choice and effect of these values will be
discussed in Section 4). The post-break slope can be affected by the strength of the guide field,
see e.g. Ref. [19]. The break Lorentz factor γp,br is dictated by energy conservation arguments,
namely the injected luminosity in the accelerated protons, Lp is at most comparable to the
Poynting flux, SP , times the area of the current sheet. Assuming that all protons from the
upstream plasma are accelerated by the reconnection process, it follows that γp,br ≃ σp [20],
where the proton magnetization value is defined as
B2 n± me
σp = 2
= σe , (2.9)
4πnp mp c n p mp
with np being the proton number density. The maximum energy of the proton distribution is
set either by the size of the corona (Hillas confinement criterion, Ep,max = eBR [21]), or it is
radiation-limited (i.e., the energy loss rate balances the acceleration rate). Since, in the cases
we examine, the radiation-limited Lorentz factor is the lowest one, Ep,max is determined by
the balance between proton cooling and acceleration.
Using Eq. 2.8 we may express the differential power of protons injected into the corona
as
mp c2 ≤ Ep ≤ Ep,br
(
d2 Np 1,
Lp (Ep ) ≡ Ep = L0  Ep −sp +1 (2.10)
dtdEp Ep,br , Ep,br < Ep ≤ Ep,max .

–4–
where Ep,br = σp mp c2 ≫ Ep,min , and L0 is a normalization that relates to Lp as
(   −1
−1 E
Ep,br 1 + ln Ep,max , sp = 2
L0 ≈ Lp p,br (2.11)
(2Ep,br )−1 , sp = 3.

3 Analytical Estimates

In this section we derive analytical expressions for the neutrino energy spectrum, and the
total secondary pair density produced in the coronal region due to proton-photon and photon-
photon interactions, as a function of the black hole mass and the coronal X-ray luminosity. We
also examine for what parameter values the proton power transferred to secondary populations
is saturated (calorimetric limit). Our analytical expressions, which are derived under several
simplifying assumptions, will be compared against numerical results in Section 4.

3.1 Neutrino luminosity and peak energy

High-energy neutrinos are produced through interactions of protons with X-ray photons from
the corona. We can estimate the neutrino luminosity as

3
Eν Lν+ν̄ (Eν ) ≈ fpγ (Ep )Ep Lp (Ep )|Ep ≈20Eν , (3.1)
8
where fpγ is the photomeson production efficiency for a proton of energy Ep , and is defined
as
 
−1 R
fpγ (Ep ) = min 1, tpγ (Ep ) . (3.2)
c

The corona is said to be optically thick to photomeson production processes when the first
term in the curly bracket dominates, and optically thin otherwise. For the pγ energy-loss
timescale tpγ , we can use Eq. 9 from Ref. [13] and write the second term in the curly brackets
as

R LX 1 min(Ep , Ep,∗ ) λEdd,−2 min(Ep , Ep,∗ )

t−1
pγ (Ep ) = 4σ̂pγ 3 2
≃ 0.3 (3.3)
c me c 9 ln(10)R̃Rg ϵr mp c R̃ 25 TeV

where ϵr = 390 is the energy threshold of the interaction (in units of me c2 ), σ̂pγ ≈ 70 µb is
the effective cross section [22], and Ep,∗ is the characteristic energy of protons interacting at
threshold with the lowest energy coronal photons,

ϵr mp me c4  ϵ
X,min
−1
Ep,∗ = ≃ 930 TeV . (3.4)
2ϵX,min 0.1 keV

When Ep > Ep,∗ the photomeson production efficiency becomes approximately independent
of the proton energy. Substitution of expressions 3.2 and 3.3 into Eq. 3.1 leads to
 
3 λEdd,−2 min (Ep , Ep,∗ )
Eν Lν+ν̄ (Eν ) ≈ min 1, 0.3 Ep Lp (Ep )|Ep ≈20Eν (3.5)
8 R̃ 25 TeV

–5–
or, after using Eqs. 2.10 and 2.11 and sp = 3,

Eν

 
 , Eν ≤ Eν,br
3ηp λEdd,−2 min (Eν , Eν,∗ )  Eν,br
Eν Lν+ν̄ (Eν ) ≈ LX min 1, 1.2
Eν −1
 
16ηX R̃ 5 TeV
, Eν > Eν,br



Eν,br
(3.6)
where Ep ≈ 20Eν , Ep,∗ ≈ 20Eν,∗ and Eν,br ≈ 5σp,5 TeV. When evaluated at Eν = Eν,br ,
Eq. 3.6 provides an approximation for the bolometric neutrino luminosity, which is de-
rived in Appendix A). For coronae optically thin to photomeson production, we expect
Eν,br Lν (Eν,br )/LX ∝ λEdd , or to be constant otherwise.

3.2 Calorimetric limit for pγ interactions

When the efficiency of pγ interactions becomes equal to unity, then all the available energy
of the parent proton population is transferred to secondary particles, such as leptons and
neutrinos. Any further increase in the number density of target photons, for fixed proton
luminosity, does not lead to an increase of the secondary injection rate. This saturation
indicates the calorimetric limit of protons in the corona (equivalently the corona becomes
optically thick to pγ).
Starting from the condition fpγ (E)|E=min(Ep,br ,Ep,∗ ) = 1 and using Eq. 3.2 we find that
the calorimetric limit is achieved when the Eddington ratio of the source is

 −1
min(Ep,br , Ep,∗ )
λEdd,crit ≃ 0.03R̃ . (3.7)
25 TeV

This critical value of the Eddington ratio scales as ∝ σp−1 , as long as Ep,br < Ep,∗ , since
Ep,br ∝ σp . Alternatively, we can define a critical neutrino energy, Eν,crit , that indicates if
the proton population in the corona is calorimetric to pγ interactions for a given value of
λEdd . Using the relation Ep ≈ 20Eν and re-arranging Eq. 3.7 we find

R̃
Eν,crit ≃ 4 TeV . (3.8)
λEdd,−2

where we assumed that Eν,br ≤ Eν,∗ (this holds for σp ≤ 106 ).

If Eνcrit ≲ Eν,br , the neutrino production at its peak is already calorimetric, meaning
that the majority of the proton energy distribution is transferred to secondaries.

3.3 Secondary pair number density

In order to estimate the total density of secondary pairs created in the corona environment
indirectly through hadronic interactions, we take into account three pair creation channels,
as outlined below1 :
1
The symbol γ without subscript is used to indicate the target photons for each interaction.

–6–
 
γ−γ annihilation


 p+ + π 0 → γTeV + γTeV −−−−−−−−−−→ e+ + e−
or




n + π → µ+ + ν̄µ
0 +



(i) p+ + γ → or
p + π + π → µ+ + ν̄µ + µ− + νµ
+ − +








µ+ (µ− ) → e+ (e− ) + ν̄e (νe )



(ii) p+ + γ → p+ + e+ + e−
cool
(iii) γMeV + γ → e+ + e−

The first (pγ) channel injects pairs in the corona indirectly via the attenuation of TeV
photons produced in neutral pion decays, generated by relativistic protons, and directly via
the decay of charged pions. Additionally, Bethe-Heitler proton-photon interactions (second
channel) act as a direct pair production channel. Lastly, the third (γγ annihilation) channel
injects pairs through the attenuation of MeV photons resulting from the electromagnetic
cascade (namely, photons emitted by the relativistic pairs injected in channels (i) and (ii)).
Because of radiative losses the injected relativistic pairs cool down reaching Lorentz factors
of γe ∼ 1, forming the pool of non-relativistic pairs.
Simple physical arguments show that the γγ annihilation channel vastly dominates over
the others in terms of pair production. pγ injects an energy content in TeV photons which is
nearly immediately reprocessed into a comparable energy content in MeV photons. However,
for the same energy content, the number density of TeV photons injected is much smaller than
the number density of MeV photons, by about 6 orders of magnitude. Therefore, the amount
of pairs that can be injected by pγ pair production directly is roughly suppressed by 6 orders of
magnitude compared to the γγ injection channel. Bethe-Heitler injection is also subdominant
compared to the γγ annihilation channel, because the energy injected directly into pairs by
the former channel is less than the energy ending up to MeV photons. In Appendix B, we
explicitly evaluate the amount of pairs injected by pγ and Bethe-Heitler, confirming these
order-of-magnitude estimates. In the following, we present in detail the predictions for the
γγ annihilation channel.

3.3.1 The γγ annihilation channel

The annihilation of photons with energies Eγ ∼ 1 − 10 MeV (henceforth, soft γ rays) will
lead to the creation of pairs. Assuming that each particle takes half of the energy of the
γ-ray photon, E± ≈ Eγ /2, and that all of the soft γ-ray luminosity is transferred to pairs2 ,
LMeV ≈ Lγγ± , we can calculate the injection rate of relativistic pairs as

Lγγ
±
ṅ± = , (3.9)
E± V

2
This implies that the corona is optically thick to soft γ-rays. In Appendix D we present numerical results
for the γγ opacity of the corona supporting our assumption.

–7–
where V = 4πR3 /3 is the corona volume. Because the cooling timescale of the γγ injected
pairs is short compared to the particle escape timescale from the corona3 , which is taken to
be R/c, Eq. 3.9 can also be used for the injection rate of non-relativistic pairs. Then, the
number density of non-relativistic pairs can be estimated as
R 3LMeV
n± = ṅ± = . (3.10)
c 2πR2 cEγ
The next step is the estimation of the MeV photon luminosity. The process that dom-
inates the production of soft γ-ray photons depends on the system parameters, like σp (this
will become clearer by our numerical results presented in the next section). One possibil-
ity is that the synchrotron emission of Bethe-Heitler pairs dominates the MeV range (see
for example figure 2). Using approximate expressions for the relevant effective cross sec-
tions, it can be shown that Lν+ν̄ /LBH,syn ≈ 33 (nX,pγ /nX,BH ), where nX,pγ and nX,BH are
the target photon number densities for photomeson and Bethe-Heitler production processes
respectively (for more details see Ref. [23]). When nX (E) ∝ E −2 , as assumed here, then
Lν+ν̄ /LBH,syn ≈ 33 (EX,BH /EX,pγ ) ≈ 0.23 (i.e. of order unity), with EX,BH = 2me c2 /γp
and EX,pγ = mπ c2 /γp ≃ 290me c2 /γp being the target photon energies for on-threshold
Bethe-Heitler and photomeson interactions respectively. Similar estimates can be made if
the emission is dominated by the radiation of secondary pairs from photomeson interactions.
Therefore, regardless of the origin of soft γ-ray photons, we can use the approximate relation
3
LMeV ≈ Lν+ν̄ ≈ fpγ (Ep,br )Ep,br Lp (Ep,br ). (3.11)
8
Using Eqs. 3.10 and 3.11 the number density of the pairs produced by the interaction of
photons with energy Eγ with the coronal ones can be estimated as:
λ2 λEdd,−2 min(Ep,br , Ep∗ )
   
9ηp fpγ LX −3 Edd,−2 10 MeV
nγγ
± = = 5.4 · 1012
cm min 1, 0.3 .
32πηX Eγ R2 c R̃2 LX,43 Eγ R̃ 25 TeV
(3.12)
where we used 10 MeV as the nominal value for the γ-ray photon energy, where the corona
becomes optically thick to γγ annihilation (for more details, see Appendix D).
The total pair density can be estimated by the sum of the pγ (Eq. B.4), Bethe-Heitler
(Eq. B.8) and the γγ annihilation (Eq. 3.12) terms. However, inspection of the three afore-
mentioned expressions shows that the attenuation of soft γ-ray photons (1 − 10 MeV) is the
most important contributor to the density of secondary pairs, in agreement with the qualita-
tive arguments presented above.

4 Numerical results
In this section we present the numerical results of our parametric study, which was performed
using the leptohadronic numerical code ATHEνA [24, 25]. The code solves a system of cou-
pled partial differential equations for all stable particle species inside a spherical, magnetized
region, namely electrons and positrons (pairs), protons, neutrons4 , photons, and neutrinos.
All parameter values are listed in Table 1.
3
In most cases synchrotron cooling dominates the radiative losses of pairs, and ctsyn /R =
6πme c2 /(σT B 2 R)(γe /10)−1 ≪ 1. This is also validated by our numerical results in section 4.
4
Neutrons have a decay time of ∼15 minutes, in their own rest frame. However, they are considered stable
particles, in our calculation, due to the fact that in the lab rest frame they have a decay timescale much longer

–8–
 Table 1: Model parameters. The top three are the key parameters that we vary.

Parameter Symbol (unit) Value range

Black hole mass Mbh (M⊙ ) 105 − 109
Broadband (0.1-100 keV) X-ray luminosity LX (erg s−1 ) 1040 − 1044
Proton plasma magnetization σp 103 − 107
Proton post-break slope sp 2 or 3
Maximum proton Lorentz factor γp,max 108
Size of reconnection layer R (Rg ) 3
Reconnection rate βrec 0.1
X-ray photon index sX 2
X-ray to magnetic energy density ratio ηX 0.5
Non-thermal proton to Poynting luminosity ratio ηp 1/3

Figure 2: Decomposition of the spectral energy distribution (SED) of the hadronic cascade
emission (black solid line) for σp = 105 (left panel) and σp = 107 (right panel). The black
dashed line represents the all-flavor neutrino energy distribution. The dash-dotted grey line
indicates the corona X-ray spectrum. Other parameter used here are: sp = 2, LX = 1043 erg
s−1 and Mbh = 107 M⊙ (λEdd = 10−2 ).

In figure 2 we show indicative electromagnetic (solid lines) and all-flavor neutrino (dashed
lines) spectra for the AGN coronal region as described by our model for two values of σp .
Black solid line shows the photon spectrum escaping the corona with all leptonic and hadronic
processes contributing to it, while yellow solid line represents the same spectrum without
accounting for γγ annihilation. Moreover, each of the other colored lines represents one
individual hadronic spectral component: the orange one refers to the emission from secondary
pairs produced in photomeson interactions, the blue one indicates the Bethe-Heitler pair
emission component, and the magenta one corresponds to the proton synchrotron spectrum.
We observe that the cascade spectra escaping the corona (solid black lines) in both panels
are similar. The aforementioned characteristic is because the escaping photon spectrum is
shaped by photon-photon pair creation, which has no direct dependence on σp . In particular,
high-energy photons (E ≳ 10 MeV) interact with lower energy ones creating relativistic
electron-positron pairs that cool radiatively mainly due to synchrotron losses, resulting in an

than the escape time, R/c.

–9–
Figure 3: Characteristic timescales as a function of the proton Lorentz factor for the pa-
rameters of figure 2. Solid lines represent the timescales for the acceleration process (black),
proton synchrotron cooling (magenta), photomeson production (orange), and Bethe-Heitler
pair production (blue). The X-ray photons of the corona are considered as targets for the
two latter processes. The dotted black line represents the Hillas-limited Lorentz factor, while
magenta and orange dotted lines represent the proton Lorentz factor at which the acceleration
timescale is equal to the synchrotron and photomeson ones, respectively.

extended and almost flat spectrum up to 10 MeV (black line). However, the contribution
of the individual hadronic components to the total spectrum is vastly different in the two
cases. For σp = 105 , protons that carry most of the energy of the population, i.e., those with
energy Ep,br = σp mp c2 , interact with X-ray photons of energy 0.1 keV close to the threshold
energy for Bethe-Heitler pair production. For such interactions the proton loss timescale due
to Bethe-Heitler pair production is the shortest, and also comparable to the photomeson loss
timescale – see figure 3. As a result, the energy transferred from protons with energy Ep,br
to secondary leptons is similar for the two interactions, hence their radiative output will be
comparable – see figure 2. The Bethe-Heitler synchrotron spectrum has a well defined peak
at ∼ 1 MeV and is produced by pairs injected with γe ≈ σp (for more details about the
Bethe-Heitler spectra, see [23]). Meanwhile, pion decays produce on average more energetic
pairs with γe ≈ 0.05 (mp /me )σp (see Eq. B.3 and [26]) that emit photons of energy ∼ 10 GeV.
The proton synchrotron spectrum is subdominant, as the relevant loss timescale is longer
than the other processes. The relative importance of proton synchrotron radiation, Bethe-
Heitler pair production and photomeson production in the cooling of protons changes for
higher proton break energies. When σp = 107 , protons cool preferentially via the photomeson
production process, which has the shortest timescale. Note that the photomeson timescale
becomes almost independent of the proton energy for interactions away from the threshold
value of the interaction (see also Eq. 3.4), as opposed to the Bethe-Heitler timescale (its
effective cross section increases logarithmically with interaction energy [25]). Moreover, the
proton synchrotron timescale becomes comparable to the photomeson loss timescale. There-
fore, for σp = 107 we expect that the synchrotron spectra from protons and pairs from pion
decays will have a comparable contribution to the total emission, while the Bethe-Heitler syn-

– 10 –
Figure 4: Left panel: Energy distribution of protons at injection (transparent solid lines)
and at steady state (bold solid lines). Steady state neutron energy distributions are also
shown (transparent dashed lines). The color bar refers to the proton magnetization value,
which determines the peak of the injected proton energy distribution (see Eq. 2.8). Right
panel: Spectral energy distribution of cascade photons (solid lines) and neutrinos of all flavors
(dashed lines) emitted from the coronal region. Other parameters used: LX = 1043 erg s−1
and Mbh = 107 M⊙ (λEdd = 10−2 ).

chrotron component should be subdominant, as indeed shown in the right panel of figure 2.
We note that in order to show the Bethe-Heitler synchrotron component separately, we re-run
ATHEνA only with electron synchrotron and Bethe-Heitler pair production on. As a result,
protons do not lose part of their energy due to proton synchrotron or photopion production
(as they normally do – see black and yellow curves). As a result, there is more energy available
for the Bethe-Heitler produced pairs, resulting in the blue curve being slightly higher than
the yellow one (which includes all the interactions apart from γγ annihilation).
We study next the impact of σp on the photon cascade and neutrino emission from the
AGN coronal region. Figure 4 shows our numerical results for LX = 1043 erg s−1 , Mbh =
107 M⊙ , sp = 3, and varying σp in the range 103 −107 . We do not explore even lower σp values
as these would lead to peak neutrino production at energies ≪ 1 TeV, where observations are
dominated by the neutrino atmospheric background. In the left panel we present the injected
and steady-state proton distributions, as well as the steady-state neutron distribution in the
corona. In the right panel we display the electromagnetic and neutrino spectra emitted from
the coronal region. As σp increases the cascade spectrum becomes more luminous, because
the production rate of secondary pairs is higher; the energy lost by the proton population
increases with σp as shown in the left panel (see also figure 3). Moreover, the electromagnetic
spectra show a sharp cutoff at about 10 MeV due to γγ pair production on the highest energy
X-ray photons (100 keV) for all σp values. Similarly, the dip of the spectrum at about 10 GeV
is caused by the attenuation of γ-ray photons by the 0.1 keV photons of the corona that
have larger number densities (for more details about the attenuation of energetic photons, see
Appendix D). Such attenuation leads to creation of secondary leptonic population which later
on radiate through synchrotron, producing broad and almost flat (in ELE representation)
electromagnetic spectra (E ∈ [10−3 , 107 ] eV, see figure 2). Given that in our model γγ
pair creation of the aforementioned populations occurs for all σp values, AGN coronae with
different luminosity or black-hole masses will be characterized by similar cascade spectra as
those shown in figure 4.

– 11 –
Figure 5: Fraction of the proton (solid dark blue line), photon (dash-dotted orange line),
neutrino (dashed light blue line), and neutron (dotted red line) bolometric energy density
over the total energy density of all four species as a function of σp . The corona luminosity
in all panels is LX = 1043 erg s−1 but the Eddington ratio changes from 10−1 in panel (a) to
10−4 in panel (d) in decrements of 10.

The peak of the neutrino energy distribution, which is determined by protons with
energy Ep,br = σp mp c2 , moves to higher energies as σp increases, as shown on the right
panel of figure 4. Moreover, the peak neutrino luminosity increases for higher σp values (see
Eq. 3.6), for the same reason the cascade emission becomes overall more luminous: for higher
σp values, more energy is carried by protons of the distribution that fulfill the energy threshold
condition for pion production, resulting in a larger fraction of proton energy transferred to
secondary particles (neutrinos, photons, and pairs). As the coronal proton population tends
to become calorimetric to photomeson interactions, the cascade and neutrino luminosities
tend to saturate to a constant value. Furthermore, as σp increases, the overall shape of the
neutrino spectrum changes, with the extent of the power law above the peak neutrino energy
decreasing. The latter is due to the fact that as σp increases the break energy of the proton
energy distribution moves towards its upper energy limit, which is set by radiative losses.
To demonstrate the impact of the Eddington ratio, λEdd , we performed a series of nu-
merical runs for an accreting system with LX = 1043 erg s−1 and Mbh = 106 − 109 M⊙ ,
and computed the energy densities of protons, and secondary particles at steady state in the
corona for σp = 103 −107 . Our results are presented in figure 5. We find that as λEdd decreases
(from panel a to panel d), the amount of energy transferred from the proton population to
the secondary particles, namely photons, neutrons and neutrinos, also decreases for all values

– 12 –
Figure 6: Non-relativistic secondary pair number density (left column) and all-flavor (bolo-
metric) neutrino luminosity over corona X-ray luminosity (right column) as a function of
the Eddington ratio (defined with respect to the broadband X-ray luminosity). Results for
σp = 105 and 107 are shown in the upper and lower panels, respectively.

of σp . In the system with the lowest λEdd (panel d), protons do not lose a significant amount
of their energy and thus, the system is optically thin to photohadronic interactions for all
σp values. On the contrary, in the system with the highest λEdd (panel a) protons transfer
most of their energy to secondary populations, with most of it going to neutrons and photons.
Furthermore, by comparing panels (a) and (b), we observe that the σp value at which protons
start to significantly lose energy increases almost linearly with λ−1
Edd (see also Eq. 3.7). For
example, in panel (a) where λEdd = 10−1 , the protons start to lose a noticeable fraction of
their energy (∼ 0.5 of their total energy) at σp ≈ 105 , while in panel (b) where λEdd = 10−2 ,
the same thing happens for σp ≈ 106 , see Eq. (3.7).
We examine next whether proton-photon interactions and γγ annihilation in the coronal
region are able to create enough pairs to result in Thomson optical depths of 0.1 to 10. In
order to do so, we show, in the left-hand side panels of figure 6, the secondary non-relativistic
pair number density, n± , computed numerically as a function of the Eddington ratio, λEdd ,
for two values of σp . Colored markers represent the numerical results of our parametric study
for cases of (i) varying compact object mass, Mbh = 105 − 109 M⊙ (mass range indicated
by the colorbar) and (ii) varying coronal luminosity, LX = 1040 − 1044 erg s−1 (luminosity
values plotted with different markers – see legend of figure 6). Furthermore, the color-shaded
bands represent the pair density values that correspond to τT = 0.1 − 10 for each black
hole mass. Changes in either one of these parameters affect the number density of X-ray

– 13 –
Figure 7: All-flavor (bolometric) neutrino luminosity over corona X-ray luminosity as a
function of Thomson optical depth due to secondary pairs for cases with LX = 1043 erg s−1
and varying black hole mass (see colorbar), for two values of σp (see legend).

photons in the corona, and in turn the opacities for proton-photon interactions and γγ pair
production. When the secondary pair density is plotted against λEdd it can be described by
a broken power-law function, with the break occurring when protons become calorimetric in
the corona (see Eq. 3.7 and vertical dotted line). For comparison purposes, we overplot the
analytical expression for n± , as given by Eq. 3.12, when LX (dashed green lines) or Mbh
(dashed black lines) is varying separately. We observe that our analytical approximation
for the pair density in the coronal region, as given by Eq. 3.12, is in good agreement with
the numerical results, verifying our assumptions that the dominant source of pairs is γγ
annihilation of energetic photons produced indirectly via proton-photon interactions. The
above result is valid both for the case where we vary the mass of the central black hole
while we keep the broadband coronal luminosity fixed at LX = 1043 erg s−1 (color-changing
crosses) and the scenario in which the coronal broadband luminosity varies while the black
hole mass is fixed to Mbh = 107 M⊙ (dark green shape-changing symbols). We also highlight
that our estimation of the calorimetric limit describes well enough the change of slope of the
dashed lines, marking the transition from an optically thin to an optically thick corona to pγ
interactions. In addition, we observe that the numerical n± values fall within their respective
color-shaded bands for λEdd ≳ λEdd,crit , while for λEdd ≲ λEdd,crit the density of secondary
pairs in the corona is not high enough as to result in optical depths of order unity. Therefore,
the dominant source of pairs in the coronae of black holes accreting at low Eddington ratios
cannot be of hadronic origin.
Furthermore, the right-hand side panels of figure 6 show that the fraction of the coronal
luminosity transferred to neutrinos is proportional to λEdd for an optically thin source and
becomes constant for an optically thick one, as described in detail in section 3.1. Moreover,

– 14 –
the fraction of luminosity transferred to neutrinos is the same for sources with the same λEdd
(and σp , even though the black hole masses and broadband X-ray luminosities are different
(see, for example, overlapping green squares with red crosses in right panels of figure 6). By
comparing the upper and lower panels on the right-hand side of the figure we observe that the
fraction of LX transferred to neutrinos is higher for higher σp values (and fixed λEdd ). This
trend can be understood as follows. For lower σp values, such as 105 , protons that carry most
of the population energy, i.e. those with energies of Ep,br , are close to threshold for Bethe-
Heitler interactions with the coronal photons (see also proton timescale plot in figure 3). As a
result, a non-negligible amount of the energy available in the proton population is channelled
to secondary Bethe-Heitler pairs, leaving less energy available for neutrino production. The
aforementioned condition is not taken into account in our analytical calculations (see instead
Eq. 6 in Ref. [27] for a more accurate treatment), resulting in the over-prediction of Lν+ν̄ /LX
for σp = 105 (see top right panel of figure 6). Nonetheless, the analytical predictions of the
ratio (see Appendix A) for the optically thin regime, i.e. for λEdd ≲ λEdd,crit , are in very good
agreement with the numerical results, and with the analytical approximation as derived by
[13] (solid black line).
We can combine the information shown in figure 6 into a diagram showing the bolometric
neutrino luminosity over the broadband X-ray luminosity as a function of the Thomson optical
depth attributed to secondary pairs (see figure 7). We observe that as λEdd increases (see
colorbar) both the bolometric neutrino luminosity and the Thomson optical depth increase.
Interestingly, for those cases with 0.1 ≲ τT,sec ≲ 10, the bolometric neutrino luminosity is
≳ 10−2 LX , with the highest values expected when σp ∼ 107 .

5 Neutrino predictions for Seyfert galaxies

The coronal region of AGN has been proposed as a production site of high-energy neutrinos
since the late seventies [e.g. 28, 29]. The IceCube Collaboration has recently presented com-
pelling evidence of high-energy neutrino emission (significance of 4.2σ) from NGC 1068 [30],
a nearby Seyfert II galaxy located at a distance of ∼10.1 Mpc (for a review, see [31]). The
best-fit neutrino spectrum is soft, dΦν /dEν ∝ Eν−3.2 , extending from 1.5 to 15 TeV. These
experimental results can be reproduced by the model presented in section 2, as demonstrated
in detail in Ref. [13]. More recently, the IceCube Collaboration performed a stacking analy-
sis of 27 additional Seyfert galaxies, which are part of the BAT AGN Spectroscopic Survey
(BASS) [32], in search of a neutrino signal [33]. They report excesses of neutrinos (at 2.7σ
significance compared to the background expectation) from two Seyfert galaxies, NGC 4151
and CGCG 420-015. Furthermore, NGC 4945 is a close-by (located at 3.8 Mpc) obscured
Seyfert II galaxy and a very bright emitter in the infrared energy regime, with that radia-
tion assumed to be the reprocessed UV and optical emission originating from the center of
the source, near the massive black hole. The aforementioned source, also, emits in the hard
X-rays, being comparably luminous to NGC 4151 [34]. Such characteristics make NGC 4945
an interesting target for future instruments, such as IceCube-Gen2 [35] and KM3Net [36].
We apply the model as described in section 2 with the difference that the maximum
proton energy is not fixed at 108 GeV. Instead, the latter is set by the minimum value
between the Hillas criterion [21], and the equilibrium between the characteristic acceleration
time and the characteristic cooling time (see for example figure 3). For NGC 1068, NGC 4151,
CGCG 420015 and NGC 4954 the aforementioned limits are 9.0 · 109 GeV, 4.4 · 1010 GeV,
1.4 · 1010 GeV and 1.6 · 1010 GeV, respectively.

– 15 –
 Table 2: Properties of Seyfert galaxies that are likely associated with high-
energy neutrinos.

Name Distance L2−10 keV L†X Mbh λEdd Refs.

(Mpc) (erg s−1 ) (erg s−1 ) (M⊙ )
NGC 1068 10.1 8.5 · 1042 3.6 · 1043 6.7 · 106 5 · 10−2 [I]
NGC 4151 15.8 1.0 · 1044 4.3 · 1044 3 · 107 ‡ 10−1 [II]
CGCG 420-015 128.8 2.0 · 1042 8.8 · 1042 2 · 108 5 · 10−4 [III]
NGC 4945 3.8 1.2 · 1042 5 · 1042 4 · 106 4 · 10−2 [IV]
†
Based on the intrinsic 2-10 keV X-ray luminosity by [37] and recalculated to the 0.1-100
keV band assuming a photon index sX = 2.
‡
Average value from different methods: Mbh = (3.76±1.15)·107 M⊙ [from stellar dynamics,
38], Mbh = 3.0+0.8 7 +0.75 7
−2.2 · 10 M⊙ [from gas dynamics, 39], and Mbh = 2.45−0.55 · 10 M⊙ [from
reverberation mapping, 40].
Parameter referencing: [I]: distance by [31] and mass by [31] [II]: mass by [40] and distance
by [41] [III]: mass and distance by [42] [IV]: mass by [43] and distance by [44].

In figure 8 we present the neutrino spectra predicted by our model for the galaxies of
Table 2 and σp = 103 − 107 . The latter is the main free parameter of the model for a given
source. We also overplot, for each source, the best-fit (all flavor) neutrino luminosity and
the 68% statistical uncertainty (grey-shaded bands) based on a power-law model as obtained
by IceCube5 [33]. Our model can satisfactorily describe the neutrino luminosity observed
by IceCube for the cases of NGC 1068 and NGC 4151. However, the preferable proton
magnetization value which better describes the observation data, differs between sources. For
NGC 4151, the neutrino spectra obtained for σp = 103 and 104 enclose the IceCube spectrum
with its statistical uncertainty. A value closer to 103.5 would better describe the observational
data. An accurate determination of σp would require fitting the IceCube neutrino events under
the assumption that the signal is described by our model, a procedure that lies beyond the
scope of this paper. For NGC 1068 the neutrino spectrum that is closer to the observed one
is obtained for σp = 105 in agreement with [13]6 , and (ii) the slightly different values of ηp
adopted. In panel (c) we display the case of CGCG 420-015 for which the model can predict
neutrinos with energies above 1 TeV, but with a luminosity, which cannot be fine-tuned, many
orders of magnitude lower than the one inferred by IceCube. We attribute the low predicted
neutrino luminosity to the fact that the corona of this source is less compact (has an X-ray
luminosity comparable that of the other sources and the largest layer among the three due to
larger black hole mass, thus the lowest λEdd and ℓX ). As a result, the corona is optically thin
to photomeson interactions, and the efficiency of energy transfer from protons to neutrinos is
low.
In figure 9 we show the ratio of the peak neutrino luminosity over the broadband X-ray
luminosity, as produced by our model, and in comparison to observational data for NGC 1068,
NGC 4151 and CGCG 420-015. Filled black symbols represent the numerical values obtained
by evaluating the model neutrino luminosity spectrum at Eν,br . For NGC 1068, NGC 4151
and CGCG 420-015 we use σp values (105 , 103 , and 106 respectively). The respective model
neutrino spectra are close to the power-law IceCube fits (grey shaded area in figure 8) in terms
5
We convert the reported energy fluxes to luminosities by using the luminosity distances listed in Table 2.
6
Ref. [13] find a neutrino component for NGC 1068, with a morphology very similar to the one of figure 8.
Minor differences between the two model spectra arise due to (i) the different σp values used, with [13] fine
tuning the latter to better match the minimum energy of the IceCube bow-tie

– 16 –
Figure 8: Numerical results of all-flavor neutrino spectra for NGC 1068 (panel a), NGC 4151
(panel b), CGCG 420-015 (panel c) and NGC 4945 (panel d), computed for a range of σp
values (see color bar). The bold solid line represents the model neutrino spectrum that is the
closest one (from those displayed) to the IceCube data (grey-shaded bands by [33]), while the
yellow area represent the best fit and 1σ deviation from Ref. [30].

of the observed energy range and/or neutrino luminosity. Open symbols refer to observational
values that were computed using the maximum (peak) value of the best-fit power-law model
for the neutrino flux, Φ(Eν ) ∝ E −γ , acquired by [33], the luminosity distances, and black
hole masses listed in Table 2. The vertical bars on the open symbols indicate the range of
minimum and maximum fluxes obtained when using the 1σ contours of the power-law fits
presented in figure 8 of Ref. [33]. Differences between the open and black filled points originate
from deviations of the predicted neutrino spectrum from the power-law model assumed by
IceCube [33]. Moreover, the colored curves in figure 9 represent the analytical approximation
of Eq. 3.6 for a coronal environment that is optically thin to pγ interactions (fpγ < 1 – see
section 3 and Eq. 3.7), while the black dashed line describes the system when it is in the
optically thick regime.
Comparing our numerical results against the analytical expressions, we notice that
NGC 1068 is closer to the optically thick regime, while the other two sources have low
photomeson production efficiency (optically thin regime). In particular, CGCG 420-015 is
a source with low Eddington ratio. Its neutrino luminosity as inferred by IceCube is com-

– 17 –
Figure 9: Ratio of peak neutrino luminosity and broadband X-ray corona luminosity as a
function of the Eddington ratio. Colored curves are computed using the analytical approxi-
mation of Eq. 3.6 estimated on the approximate peak neutrino energy, Eν,br , for various σp
values (see color bar). The horizontal dashed black line represents the calorimetric limit (Eq.
3.6). The various symbols represent the numerical results (filled black) shown in figure 8 and
observational values (open) which were acquired by [33] and Table 2. The vertical bars on
the open symbols are computed from the 1σ contours of the power-law fits presented in [33],
as explained in the text.
.

parable to the X-ray luminosity of the source, exceeding the theoretical prediction even for
the calorimetric limit. In other words, despite CGCG 420-015 having a X-ray luminosity
similar to that of NGC 4151, the larger black hole mass implies a much wider emission region
which is much less prone to photohadronic emission. Moreover, it is interesting to note that
while NGC 1068 and NGC 4151 have similar Eddington ratios, our model can accommodate
the difference in the observed values of Eν,br Lν (Eν,br )/LX as a result of different σp values.
Therefore, if one can observationally determine LX , Mbh , and Eν,br Lν (Eν,br ), one can use this
diagram as a guideline to infer the proton magnetization in the coronal environment. Because
σp also relates to the peak energy of the neutrino spectrum, as shown in figure 8, the inferred
σp value can be then cross checked against neutrino spectra, whenever available.

6 Summary and Discussion

We introduced a model for reconnection-powered AGN coronae. We assumed that a fraction of

the dissipated magnetic energy via reconnection is used to power the X-ray coronal emission,
whose spectrum was described by a power law between 0.1 keV and 100 keV with photon
index 2. We also considered that protons existing in the coronal environment are accelerated
in magnetic reconnection regions near the central black hole. The number density distribution
of the injected hadrons into the corona is described by a broken power law, with slopes of −1
and −3 respectively before and after the break energy σp mp c2 , where σp is the proton plasma

– 18 –
magnetization in the corona. These relativistic protons interacting with coronal photons
produce secondary pairs and high-energy neutrinos. The pairs can contribute to enrich the
lepton population in the corona, while the neutrinos can provide a novel signal visible at
Earth. Based on the aforementioned model, we performed a parameter scan of the coronal X-
ray luminosity, black hole mass, and proton plasma magnetization. For each parameter set we
computed the high-energy neutrino spectrum and estimated the number density of secondary
pairs produced in the corona using the leptohadronic code ATHEνA. We complemented our
parametric study with analytical calculations that showcase the impact of the three main
physical parameters.
Describing the differential photon density spectrum of the corona as a power law with
slope −2 can be considered a simplification which, however, does not significantly alter our
results. In Ref. [45], the authors found, by fitting observational data of 228 AGN, that the
coronal photon index, sX is proportional to log10 (Lbol /LEdd ), with Lbol being the bolometric
AGN luminosity. The aforementioned dependence can be considered a weak one and, thus,
the parameter space that creates for the sX is very narrow, with sX ∈ [1, 3]. For sX > 1, most
of the coronal photons, in number, are found at the lower limit of the distribution (e.g. 0.1
keV) and are the dominant targets for photo-hadronic interactions. As a result, changes in
sX would not affect the overall picture of the corona apart from its luminosity; this can be
altered, at most, by a factor of 9 compared to the one obtained for the nominal value sX = 2.
Our results, both numerical and analytical, highlight the importance of the Eddington
ratio parameter, λEdd . Such a dependence results in the initially 3D parameter space of our
model (LX , Mbh , σp ) being reduced to a 2D parameter space (λEdd , σp ), with the X-ray
luminosity LX only fixing the overall energetics of the source. More specifically, we find that
for λEdd ≲ λEdd,crit (see Eq. 3.7) the system is optically thin to pγ interactions. As a result,
the latter are not capable of sustaining a high-enough pair density in the corona to account
for the inferred Thomson scattering depths of AGN coronae. Therefore, for λEdd ≲ λEdd,crit ,
the percentage of the coronal luminosity, LX , transferred to neutrinos is proportional to the
Eddington ratio, λEdd , while that aforementioned percentage becomes constant otherwise (see
e.g. Eq. A.4 and figure 6), corresponding to a calorimetric regime where all of the protons
dump their energy in secondary particles.
In this work we have not specified the origin of the pairs responsible for the Comp-
tonization of soft photons to the X-ray regime. Nonetheless, the pair number density in the
corona is expected to be the sum of hadronic-generated (secondary) pairs and of primary
pairs (i.e., those generated by other channels). For the parametric search we performed in
section 4 we found 10−7 ≲ τT,sec ≲ 10, with the highest values obtained in systems with the
lightest black holes and highest proton plasma magnetizations (see figure 7). Assuming that
the total Thomson optical depth (with both primary and secondary contributions) is τT ∼ 1,
we estimated the pair plasma magnetization for the cases we explored (see e.g. Eq. 2.6) and
found that 10−5 ≲ σe /σp ≲ 1.
Moreover, we comment on the dependence of our results on the post-break slope sp of
the accelerated proton population, which was taken to be 3 so far. 3D PIC simulations of
reconnection in the relativistic regime 1 ≪ σe < σp , which is of relevance to our work, have
not yet been performed. There is evidence, however, that the post-break slope is related
to the strength of the non-reconnecting magnetic field, known as the guide field. 3D PIC
simulations in electron-positron plasmas indicate steeper power laws for stronger guide fields
(see e.g. Figs. 3 and 4 in [19]). For sp = 2 the energy of the hadronic population is equally
distributed in logarithmic energy to protons with Ep ≥ Ep,br = σp mp c2 . As a result, all

– 19 –
protons above the break play an important part in pair creation and neutrino production.
We performed numerical runs for the same parameters as those used in figure 4 except for
sp = 2 (Appendix C). We find that, in this case, the change of σp does not significantly
affect the shape and luminosity of the electromagnetic spectrum or the luminosity and peak
position of the neutrino energy distribution (see figure 10).
In this study we assumed that the X-ray non-thermal emission of the corona is fixed
(i.e., it does not depend on the hadronic interactions and does not evolve with time). Under
this hypothesis, we computed the expected density of secondary pairs created by proton-
photon and photon-photon interactions. In cases where hadronic interactions are responsible
for producing a high density of secondary pairs, the assumption of a fixed coronal field has
to be modified, and a time-dependent framework should be adopted. We envision a situation
where accelerated protons start interacting with an ambient soft thermal radiation field, e.g.
the UV bump of the accretion disk, and produce secondary pairs which will contribute to the
pool of non-relativistic pairs in the corona. Such pairs can Comptonize soft photons to higher
energies. As a result, a fraction of the injected energy into secondary pairs can be transferred
to X-ray photons, that will obtain a power-law spectrum dictated by the pair density in the
corona [46]. Therefore protons will be interacting with an evolving radiation field. In this
new description of the coronal environment, it is not clear whether a steady state can be
reached, and if its emission will be characterized by a hard (photon index ≲ 2 ) non-thermal
X-ray photon field. We intend to improve our model by accounting for an evolving coronal
radiation field that will be dynamically coupled to the relativistic proton population of the
corona in a separate work.
We have, additionally, performed some numerical runs for XRB-like environments, with
the mass of the compact object being Mbh = 1 − 10 M⊙ and a broadband X-ray luminosity
of LX = 1037 − 1038 erg s−1 . The aforementioned parameters yield λEdd values equivalent
to the AGN case of LX = 1043 erg s−1 and Mbh = 106 M⊙ (discussed in Section 4). We
find that such XRB-like systems behave in a similar way as the AGN studied in this work.
More specifically, electromagnetic and neutrino spectra are shaped like the ones presented in
sections 4 and 5, but they are less luminous than those of AGN coronae; their luminosities
scale with the X-ray luminosity of the system. We attribute the similarities to the fact that
both XRB and AGN corona we studied have similar compactnesses (in terms of photons and
protons).
An interesting aspect of our model is that the neutrino spectrum from an AGN corona
can be determined if basically two parameters are known, the ratio λEdd and σp (if sp is
specified). The constraining nature of the model and its simplicity make it appealing for
the calculation of the diffuse neutrino flux from Seyfert galaxies. Contrary to other works
where NGC 1068 is used as a prototype, our model suggests that NGC 1068 may not be
representative of the population in terms of neutrino production (see e.g. figure 9). In fact,
we find that our model can correctly represent the neutrino signal inferred from NGC 1068
and NGC 4151 at IceCube only if the former has a larger σp . For the recently reported
neutrino excess at CGCG 420-015, we find that the coronal energetics is much smaller than
the one suggested by the excess, due to the larger black hole mass which widens the corona
and makes it less optically thick. Thus, should this excess be confirmed as a signal, it would
suggest a non-coronal origin for the neutrinos.
Our two-parameter model can be extended to predict the neutrino flux from the AGN
corona population. Being a superposition of sources with different properties (i.e., λEdd and
σp ), the diffuse neutrino spectrum of the population is expected to be broader than the one

– 20 –
computed using a single source as a prototype (see e.g. [47] and [48] in a different context).
A comparison of our model prediction against the most updated IceCube measurement of the
diffuse flux will be studied in a separate work.

Acknowledgments
We thank Chengchao Yuan for useful comments on this manuscript. M.P. acknowledges
support from the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the “2nd
call for H.F.R.I. Research Projects to support Faculty members and Researchers" through
the project UNTRAPHOB (Project ID 3013). DFGF is supported by the Alexander von
Humboldt Foundation (Germany) and, when this work was begun, was supported by the
Villum Fonden (Denmark) under Project No. 29388 and the European Union’s Horizon
2020 Research and Innovation Program under the Marie Sklodowska-Curie Grant Agreement
No. 847523 “INTERACTIONS.” L.C. acknowledges support from NSF grant PHY-2308944
and NASA ATP grant 80NSSC22K0667. L.S. acknowledges support from DoE Early Career
Award DE-SC0023015 and NASA ATP 80NSSC24K1238. This work was also supported
by a grant from the Simons Foundation (MP-SCMPS-00001470) to L.S., and facilitated by
Multimessenger Plasma Physics Center (MPPC, NSF PHY-2206609 to L.S.).

References
[1] J.S. Collinson, M.J. Ward, H. Landt, C. Done, M. Elvis and J.C. McDowell, Reaching the peak
of the quasar spectral energy distribution - II. Exploring the accretion disc, dusty torus and host
galaxy, MNRAS 465 (2017) 358 [1610.04221].
[2] R.A. Sunyaev and L.G. Titarchuk, Comptonization of X-Rays in Plasma Clouds - Typical
Radiation Spectra, A&A 86 (1980) 121.
[3] F. Haardt and L. Maraschi, A Two-Phase Model for the X-Ray Emission from Seyfert
Galaxies, ApJ 380 (1991) L51.
[4] A.M. Beloborodov, Plasma Ejection from Magnetic Flares and the X-Ray Spectrum of Cygnus
X-1, ApJ 510 (1999) L123 [astro-ph/9809383].
[5] B.F. Liu, S. Mineshige and K. Shibata, A Simple Model for a Magnetic Reconnection-heated
Corona, ApJ 572 (2002) L173 [astro-ph/0205257].
[6] A.M. Beloborodov, Radiative Magnetic Reconnection Near Accreting Black Holes, ApJ 850
(2017) 141 [1701.02847].
[7] L. Sironi and A.M. Beloborodov, Kinetic Simulations of Radiative Magnetic Reconnection in
the Coronae of Accreting Black Holes, ApJ 899 (2020) 52 [1908.08138].
[8] N. Sridhar, L. Sironi and A.M. Beloborodov, Comptonization by reconnection plasmoids in
black hole coronae I: Magnetically dominated pair plasma, MNRAS 507 (2021) 5625
[2107.00263].
[9] N. Sridhar, L. Sironi and A.M. Beloborodov, Comptonization by reconnection plasmoids in
black hole coronae II: Electron-ion plasma, Monthly Notices of the Royal Astronomical Society
518 (2022) 1301.
[10] IceCube Collaboration, R. Abbasi, M. Ackermann, J. Adams, J.A. Aguilar, M. Ahlers et al.,
Evidence for neutrino emission from the nearby active galaxy NGC 1068, Science 378 (2022)
538 [2211.09972].
[11] V.A. Acciari et al., Constraints on Gamma-Ray and Neutrino Emission from NGC 1068 with
the MAGIC Telescopes, ApJ 883 (2019) 135 [1906.10954].

– 21 –
[12] K. Murase, Hidden Hearts of Neutrino Active Galaxies, ApJ 941 (2022) L17 [2211.04460].
[13] D.F.G. Fiorillo, M. Petropoulou, L. Comisso, E. Peretti and L. Sironi, TeV Neutrinos and Hard
X-Rays from Relativistic Reconnection in the Corona of NGC 1068, ApJ 961 (2024) L14
[2310.18254].
[14] A. Chernoglazov, H. Hakobyan and A.A. Philippov, High-Energy Radiation and Ion
Acceleration in Three-dimensional Relativistic Magnetic Reconnection with Strong Synchrotron
Cooling, Oct., 2023.
[15] H. Zhang, L. Sironi, D. Giannios and M. Petropoulou, The origin of power-law spectra in
relativistic magnetic reconnection, Feb., 2023.
[16] H. Zhang, L. Sironi and D. Giannios, Fast particle acceleration in three-dimensional relativistic
reconnection, The Astrophysical Journal 922 (2021) 261.
[17] C. Ricci, L.C. Ho, A.C. Fabian, B. Trakhtenbrot, M.J. Koss, Y. Ueda et al., BAT AGN
Spectroscopic Survey - XII. The relation between coronal properties of active galactic nuclei and
the Eddington ratio, MNRAS 480 (2018) 1819 [1809.04076].
[18] P.O. Petrucci, D. Gronkiewicz, A. Rozanska, R. Belmont, S. Bianchi, B. Czerny et al.,
Radiation spectra of warm and optically thick coronae in AGNs, A&A 634 (2020) A85
[2001.02026].
[19] G.R. Werner and D.A. Uzdensky, Nonthermal Particle Acceleration in 3D Relativistic
Magnetic Reconnection in Pair Plasma, ApJ 843 (2017) L27 [1705.05507].
[20] L. Comisso, Concurrent Particle Acceleration and Pitch-angle Anisotropy Driven by Magnetic
Reconnection: Ion-electron Plasmas, ApJ 972 (2024) 9 [2405.18227].
[21] A.M. Hillas, The Origin of Ultra-High-Energy Cosmic Rays, ARA&A 22 (1984) 425.
[22] A.M. Atoyan and C.D. Dermer, Neutral Beams from Blazar Jets, ApJ 586 (2003) 79
[astro-ph/0209231].
[23] D. Karavola and M. Petropoulou, A closer look at the electromagnetic signatures of
Bethe-Heitler pair production process in blazars, May, 2024.
[24] S. Dimitrakoudis, A. Mastichiadis, R.J. Protheroe and A. Reimer, The time-dependent
one-zone hadronic model - first principles, A&A 546 (2012) A120.
[25] A. Mastichiadis, R.J. Protheroe and J.G. Kirk, Spectral and temporal signatures of
ultrarelativistic protons in compact sources, astro-ph/0501156.
[26] M. Petropoulou, S. Dimitrakoudis, P. Padovani, A. Mastichiadis and E. Resconi, Photohadronic
origin of γ -ray BL Lac emission: implications for IceCube neutrinos, MNRAS 448 (2015)
2412 [1501.07115].
[27] A. Mastichiadis and M. Petropoulou, Hadronic X-Ray Flares from Blazars, ApJ 906 (2021)
131 [2009.12158].
[28] D. Eichler, High-energy neutrino astronomy: a probe of galactic nuclei?, ApJ 232 (1979) 106.
[29] V.S. Berezinskii and V.L. Ginzburg, On high-energy neutrino radiation of quasars and active
galactic nuclei, MNRAS 194 (1981) 3.
[30] R. Abbasi et al., Evidence for neutrino emission from the nearby active galaxy NGC 1068,
Science 378 (2022) 538 [2211.09972].
[31] P. Padovani, E. Resconi, M. Ajello, C. Bellenghi, S. Bianchi, P. Blasi et al., Supermassive black
holes and very high-energy neutrinos: the case of NGC 1068, arXiv e-prints (2024)
arXiv:2405.20146 [2405.20146].

– 22 –
[32] C. Ricci, B. Trakhtenbrot, M.J. Koss, Y. Ueda, I. Del Vecchio, E. Treister et al., BAT AGN
Spectroscopic Survey. V. X-Ray Properties of the Swift/BAT 70-month AGN Catalog, ApJS
233 (2017) 17 [1709.03989].
[33] R. Abbasi, M. Ackermann, J. Adams, S.K. Agarwalla, J.A. Aguilar, M. Ahlers et al., IceCube
Search for Neutrino Emission from X-ray Bright Seyfert Galaxies, arXiv e-prints (2024)
arXiv:2406.07601 [2406.07601].
[34] E. Aguilar-Ruiz, N. Fraija, J.C. Joshi, A. Galvan-Gamez and J.A. de Diego, Cosmic rays,
neutrinos and GeV-TeV gamma rays from Starburst Galaxy NGC 4945, Physical Review D
104 (2021) 083013.
[35] M.G. Aartsen, R. Abbasi, M. Ackermann, J. Adams, J.A. Aguilar, M. Ahlers et al.,
Icecube-gen2: the window to the extreme universe, Journal of Physics G: Nuclear and Particle
Physics 48 (2021) 060501.
[36] S. Adrián-Martínez, M. Ageron, F. Aharonian, S. Aiello, A. Albert, F. Ameli et al., Letter of
intent for km3net 2.0, Journal of Physics G: Nuclear and Particle Physics 43 (2016) 084001.
[37] C. Ricci, B. Trakhtenbrot, M.J. Koss, Y. Ueda, I. Del Vecchio, E. Treister et al., BAT AGN
Spectroscopic Survey. V. X-Ray Properties of the Swift/BAT 70-month AGN Catalog, The
Astrophysical Journal Supplement Series 233 (2017) 17.
[38] C.A. Onken, M. Valluri, J.S. Brown, P.J. McGregor, B.M. Peterson, M.C. Bentz et al., The
Black Hole Mass of NGC 4151. II. Stellar Dynamical Measurement from Near-infrared Integral
Field Spectroscopy, ApJ 791 (2014) 37 [1406.6735].
[39] E.K.S. Hicks and M.A. Malkan, Circumnuclear Gas in Seyfert 1 Galaxies: Morphology,
Kinematics, and Direct Measurement of Black Hole Masses, ApJS 174 (2008) 31 [0707.0611].
[40] G. De Rosa, M.M. Fausnaugh, C.J. Grier, B.M. Peterson, K.D. Denney, K. Horne et al.,
Velocity-resolved Reverberation Mapping of Five Bright Seyfert 1 Galaxies, ApJ 866 (2018) 133
[1807.04784].
[41] W. Yuan, M.M. Fausnaugh, S.L. Hoffmann, L.M. Macri, B.M. Peterson, A.G. Riess et al., The
Cepheid Distance to the Seyfert 1 Galaxy NGC 4151, ApJ 902 (2020) 26 [2007.07888].
[42] M. Koss, B. Trakhtenbrot, C. Ricci, I. Lamperti, K. Oh, S. Berney et al., BAT AGN
Spectroscopic Survey. I. Spectral Measurements, Derived Quantities, and AGN Demographics,
ApJ 850 (2017) 74 [1707.08123].
[43] L.J. Greenhill, J.M. Moran and J.R. Herrnstein, THE DISTRIBUTION OF H2O MASER
EMISSION IN THE NUCLEUS OF NGC 4945, .
[44] I.D. Karachentsev, R.B. Tully, A. Dolphin, M. Sharina, L. Makarova, D. Makarov et al., The
Hubble flow around the CenA / M83 galaxy complex, The Astronomical Journal 133 (2007)
504.
[45] B. Trakhtenbrot, C. Ricci, M.J. Koss, K. Schawinski, R. Mushotzky, Y. Ueda et al., BAT AGN
Spectroscopic Survey (BASS) – VI. The γX–L/LEdd relation, Monthly Notices of the Royal
Astronomical Society 470 (2017) 800.
[46] A.P.L. George B. Rybicki, Radiative processes in astrophysics, Haruard-Smithsonian Center for
Astrophysics.
[47] P. Padovani, R. Gilli, E. Resconi, C. Bellenghi and F. Henningsen, The neutrino background
from non-jetted active galactic nuclei, A&A 684 (2024) L21 [2404.05690].
[48] M. Petropoulou, S. Coenders, G. Vasilopoulos, A. Kamble and L. Sironi, Point-source and
diffuse high-energy neutrino emission from Type IIn supernovae, MNRAS 470 (2017) 1881
[1705.06752].

– 23 –
[49] P.S. Coppi and R.D. Blandford, Reaction rates and energy distributions for elementary
processes in relativistic pair plasmas, MNRAS 245 (1990) 453.
[50] M.J. Chodorowski, A.A. Zdziarski and M. Sikora, Reaction Rate and Energy-Loss Rate for
Photopair Production by Relativistic Nuclei, ApJ 400 (1992) 181.
[51] S.R. Kelner and F.A. Aharonian, Energy spectra of gamma-rays, electrons and neutrinos
produced at interactions of relativistic protons with low energy radiation, Physical Review D 78
(2008) 034013.

A Analytical expression for Lν+ν̄

We compute the bolometric neutrino luminosity in the optically thin regime by integrating
Lν+ν̄ (Eν ) given in Eq. 3.1.
We first examine the case with Eν,br ≃ 5 σp,5 TeV ≤ Eν,∗ ≃ 46.5 TeV, that holds for
σp ≲ 106 . Then, there are three energy ranges of interest: (i) Eν ≤ Eν,br , (ii) Eν,br < Eν ≤
Eν,∗ , and (iii) Eν,∗ < Eν ≤ Eν,max ≃ 500 TeV γp,max,8 . We obtain,

" ! #
2
Eν,min
 
(thin) Eν,br 1 Eν,∗ Eν,∗
Lν+ν̄ =A 1− 2 + ln +1− . (A.1)
25 TeV 2 Eν,br Eν,br Eν,max

where Eν,min ≪ Eν,br is the energy of the less energetic neutrinos that are produced by protons
interacting near the threshold with coronal photons of energy EX,max , and

18 ηp λEdd,−2
A= LX . (A.2)
16 ηX R̃

For proton magnetizations σp > 106 , the energy regimes of interest for the integration become
(i) Eν ≤ Eν,∗ , (ii) Eν,∗ ≤ Eν < Eν,br , and (iii) Eν,br < Eν ≤ Eν,max , and the bolometric
luminosity reads
" ! #
2
Eν,min
(thin) Eν,∗ Eν,∗ Eν,∗ Eν,br
Lν+ν̄ =A 1− 2
+1− +1− . (A.3)
25 TeV 2Eν,br Eν,∗ Eν,br Eν,max

In the optically thick regime the result is simply written as

(thick) 3 ηp
Lν+ν̄ ≈ LX (A.4)
8 ηX

where we assumed Eν,min ≪ Eν,br ≪ Eν,max .

The expression derived above is valid when the timescale of the photopion production,
tpγ is significantly shorter than the timescales of all the other hadronic interactions, such as
proton synchrotron or Bethe-Heitler pair production. If the aforementioned condition is not
satisfied, the proton population loses a fraction of its energy due to other interactions apart
from photopion production. As a result the amount of energy transferred to neutrinos is
reduced in comparison to cases where tpγ ≪ tp,syn and tpγ ≪ tBH .

– 24 –
B Other channels for pair production
B.1 The photomeson (pγ) production channel
We consider the production of neutral pions (π 0 ) through interactions of protons with X-ray
photons. Neutral pions of energy Eπ0 = kpγ Ep ≈ 0.2Ep will instantaneously decay into two
γ-ray photons. The average γ-ray photon energy (in the lab frame) is taken to be half of that
of the parent pion,
 
1 Ep,br
Eγ = Eπ0 ≈ 0.1Ep,br ≈ 2.5 TeV . (B.1)
2 25 TeV
The TeV photons from pion decay may interact with low-energy photons, i.e., EEγ ≥
 −1
E
2(me c2 )2 ⇒ E ≥ 0.2 eV 25 p,brTeV , through photon-photon annihilation to produce electron-
positron pairs. The interactions of TeV photons with X-ray coronal photons (EX ≥ 0.1 keV)
will take place far from the energy threshold of γγ pair production, where the cross section
will be about three orders of magnitude lower than its maximum value (close to the thresh-
old) [49]. Nonetheless, photons produced by the secondary pairs themselves via synchrotron
radiation will provide a sufficient number of low-energy targets in order for the TeV photons to
interact close to the threshold. This will become clear by our numerical results that account
for all photon populations in the corona (see section 4).
We therefore assume that all the luminosity injected into pions is transferred to γ-ray
photons, which in turn will be attenuated, thus transferring their luminosity into secondary
pairs. We may express the pair energy injection rate as
1
Lpγ
± ≈ Lπ 0 ≈ fpγ Ep,br Lp (Ep,br ) (B.2)
2
where fpγ is defined in Eq. 3.2.
Regarding charged pion production, we also assume that the pion is created with 20
per cent the energy of the parent proton. The pion will decay, resulting in the creation of
3 neutrinos and one lepton. If the pion energy is equally distributed to all the secondary
particles, then the lepton energy upon production reads,
pγ→π ± →ee 1
E± ≈ Eπ± = 0.05Ep,br . (B.3)
4
Summing up the contributions of the neutral pion and charged pion decays, and using Eq. 3.10
we estimate the secondary pair number density as
Ep,br −1 λEdd,−2
2
η′
   
pγ 7 −3 p λEdd,−2 min (Ep,br , Ep,∗ )
n± ≈ 0.13 · 10 cm min 1, 0.3 .
ηX 25 TeV R̃2 LX,43 R̃ 25 TeV
(B.4)
In the optically thick limit, we find npγ 2 −1 −2
± ∝ λEdd LX R̃ ∝ LX R−2 whereas in the optically
thin limit the dependence on coronal parameters is stronger, npγ 3 −1 −3
± ∝ λEdd LX R̃ ∝ L2X R−3 .

B.2 Proton-photon (Bethe-Heitler) pair production

Another hadronic channel for pair production is proton-photon (Bethe-Heitler) pair produc-
tion. The bolometric luminosity transferred from the proton population to the secondary pair
population can be written as
LBH
± = fBH Ep,br Lp (Ep,br ) (B.5)

– 25 –
where fBH denotes the efficiency of this process, which is written as

fBH = min {1, σ̂BH nt R} (B.6)

where σ̂BH ≃ 8 × 10−31 cm2 is the maximum effective
cross section accounting for the
inelasticity of the interaction [50], nt = 3Lt / 4πR2 cϵt is the number density of target
photons for Bethe-Heitler, and ϵt is the target photon energy. The latter is expressed
as ϵt ≃ 16me c2 /γp ≈ 0.32 keV(25 TeV/Ep,br ) assuming interactions close to the thresh-
old [23, 25].
We also assume that the created pairs have a Lorentz factor equal to the one of the parent
proton (this Lorentz factor corresponds to the peak of the energy distribution of injected pairs
as shown in [23]),
BH me
E± = Ep,br . (B.7)
mp

Substitution of expressions B.5, B.6 and B.7 into Eq. 3.10 yields

Ep,br −1 λEdd,−2
2
η′
    
BH 8 −3 p λEdd,−2 Ep,br
n± = 1.2 · 10 cm min 1, 0.1 (B.8)
ηX 25 TeV R̃2 LX,43 R̃ 25 TeV

C Proton distribution with a post-break slope of 2

In this work we describe the proton differential number distribution as a broken power law
with a break at γp = σp , with σp being the proton magnetization parameter (see Eq. 2.8).
The pre-break slope is fixed at −1, while the post-break slope value we used was set to −3
(or sp = 3 equivalently). We compute the electromagnetic and neutrino spectra for the same
parameters as those used in figure 4, with the only change being the post-break slope of the
proton differential number distribution, which is, now, sp = 2.

Figure 10: Same as figure 4 but for a post-break slope of the proton distribution of sp = 2.

Our results are presented in figure 10. We observe that the cascade spectrum does not
highly depend on σp value when sp = 2. In this case the energy is equally distributed (in
logarithm) to protons with γp ≥ σp . Combining the latter with the almost constant value
of the photopion effective cross section [51], the protons above the break of the distribution
can find and interact efficiently with the luminous coronal photons. The almost universal
photon spectrum, especially at energies where γγ annihilation is most important as a source

– 26 –
of pairs (∼ 10 MeV, see also section 3.3), suggests that the secondary pair density will also
not depend strongly on σp .
Furthermore, we highlight that the neutrino spectrum is not significantly affected by
the change of σp . We notice that the peak neutrino luminosity changes only by a factor of 3
when σp varies by 4 orders of magnitude. Most importantly, the peak energy of the neutrino
spectrum, also, remains unchanged when we vary σp , as it is controlled by the highest energy
protons of the distribution. Therefore, if the post-break slope of the proton distribution is
sp = 2, both the neutrino peak energy and the neutrino post-break slope are insensitive to
σp , as opposed to the sp = 3 case (see figure 4). This is an important qualitative difference
that could be used to distinguish between possible post-break slopes of the accelerated proton
population.
Finally, we note that changing sp will not significantly alter the behaviour of the system
in the n± − λEdd or Lν+ν̄ /LX − λEdd parameter spaces, since both quantities depend on the
broadband corona luminosity (see Eqs. 3.6 and 3.12 and figure 6).

D γ − γ optical depth
In section 4 we presented the numerical results of our model. In particular, in the right panel
of figure 4, we show the cascade electromagnetic spectrum of the coronal region, noting that
its morphology is similar for all σp values, as a result of γγ annihilation for photons with
energies ≥ 10 MeV. In order to validate this we numerically calculate the optical depth of the
γγ annihilation process, as
(ϵγ ϵ′ )2 − 1
Z
dϵ′ n(ϵ′ ) ′


τγγ (ϵγ ) ≈ 0.625σT R ln ϵ γ ϵ (D.1)
2/ϵγ (ϵγ ϵ′ )3

where we used the approximate γγ cross section from Ref. [49] and with ϵγ = Eγ /(me c2 ),
R being the radius of the coronal region, σT being the Thomson cross section and n(ϵ)
representing the differential number distribution of target photons.
The results of the τγγ calculation for the cases of figure 4, are shown in figure 11. The
solid lines in the top left panel represent the optical depth as a function of the energy of
photons that are annihilated. In the same panel, the dashed lines show the cascade electro-
magnetic spectrum for each of the cases (in a ϵ2 n(ϵ) versus ϵ representation), while the dashed
lines in the bottom right panel refer to the number distribution of the photon population in
the corona. Different colors correspond to different σp values (see color bar). We note that
the shape of the photon number distribution is reflected to the optical depth with the peak
at 10−10 MeV being responsible for the very high opacity value at 1010 MeV and thus the
absorption of all the luminosity produced at such energies (mostly due to photomeson – see
appendix B). Therefore, the low-energy photons from the electromagnetic cascade can atten-
uate the TeV to PeV γ rays, even in the absence of other ambient dense low-energy photons
(e.g. from the disk). We, also, see a second peak in the photon distribution at energies of
∼ 10−3 MeV that corresponds to the minimum energy photons of the X-ray corona. This
secondary peak in the target photon number density is reflected to τγγ at energies ∼ 103 MeV.
Furthermore, the vertical dotted line indicates the photon energy where the system goes from
the optically thin to the optically thick regime. The aforementioned transition happens at
∼ 10 MeV and we notice that photons with higher energies are attenuated almost completely.
Moreover, in the right panel of figure 11, we show the time evolution of the optical depth
for one of the cases displayed on the left panel, the one with σp = 105 . Results are shown

– 27 –
Figure 11: Top panel: Solid lines represent the optical depth for γγ annihilation process as a
function on the high energy photon that is absorbed for the cases of figure 4. Colorbar shows
the proton magnetization value, while dotted lines represent the cascade electromagnetic
spectrum for each of the cases, normalized to its maximum. Vertical dotted line in both panels
represent τγγ = 1 and we see that photons of higher energies are significantly attenuated.
Bottom panel: Dashed lines show the number density of the photon distribution acting as
targets for γγ, for the same cases as the ones in the top panel. We note how the morphology of
the target photon number distribution reflects in the shape of τγγ . Right panel:Time evolution
of the optical depth for the case of σp = 105 , presented in the left panels. The time span
shown is one light-crossing time with a step between the snapshots of 0.01 R/c.

for a time interval of one light crossing time, R/c, since this is long enough for a steady state
to be reached. We see that initially the system is optically thin to photons with PeV energy
values, as there are no photons below 100 eV initially present in the corona. However, as time
progresses and photopion production contributes to the TeV-PeV γ-ray emission, secondary
pairs radiate low energy photons acting as γγ targets to the former, and the corona becomes
optically thick to very high-energy γ rays within one light crossing time.

– 28 –