# Constraining the baryon loading factor of AGN jets: implication from the
gamma-ray emission of the Coma cluster

Xin-Yue Shi1,2, Yi Zhang1,2, Ruo-Yu Liu1,2 and Xiang-Yu Wang1,2 1School of
Astronomy and Space Science, Xianlin Road 163, Nanjing University, Nanjing
210023, China
2Key laboratory of Modern Astronomy and Astrophysics (Nanjing University),
Ministry of Education, Nanjing 210023, People’s Republic of China

###### Abstract

High-energy cosmic rays (CRs) can be accelerated in the relativistic jets of
Active Galactic Nuclei (AGNs) powered by supermassive black holes. The baryon
loading efficiency onto relativistic CR baryons from the accreting black holes
is poorly constrained by observations so far. In this paper, we suggest that
the $\gamma$-ray emission of galaxy clusters can be used to study the baryon
loading factor of AGN jets, since CRs injected by AGN jets are completely
confined in the galaxy clusters and sufficiently interact with intra-cluster
medium via hadronic process, producing diffuse $\gamma$-rays. We study the
propagation of CRs in the galaxy clusters and calculate the radial
distribution of the gamma-rays in the galaxy cluster with different injection
rates from AGNs. By comparison with the $\gamma$-ray flux and upper limits of
the Coma cluster measured by $Fermi$-LAT and VERITAS, we find the upper limit
of the average baryon loading factor (defined as the efficiency with which the
gravitational energy is converted into relativistic particles) to be
$\eta_{p,\mathrm{grav}}<0.1$. The upper limit is much lower than that required
to account for diffuse neutrino flux in the conventional blazar models.

Galaxies: clusters: individual: Coma – Cosmic rays – $\gamma$-rays: galaxies:
clusters

##  1 introduction

Powerful relativistic jets generated by the accretions of the supermassive
black holes (SMBHs), i.e., the active galactic nuclei (AGNs), are believed to
be one of the promising acceleration sites of ultra-high energy cosmic rays
(UHECRs) (Biermann, 1988; Takahara, 1990; Rachen & Biermann, 1993; Berezinsky
et al., 2006; Dermer et al., 2009) and high-energy neutrinos (Mannheim et al.,
1992; Atoyan & Dermer, 2001; Murase et al., 2014). The baryon loading factor
$\eta_{p,\mathrm{rad}}$ is usually defined as the the ratio between total CR
luminosity ($L_{p}$) and the bolometric radiation luminosity of the jet
($L_{\gamma}$), i.e., $\eta_{p,\mathrm{rad}}=L_{p}/L_{\gamma}$. The baryon
loading factor is very relevant to the UHECR and neutrino production
processes. There have been some discussions about the baryon loading factor in
AGN jets based on the lepto-hadronic model of blazars (e.g., Böttcher et al.,
2013; Xue et al., 2019a) as well as on the potential blazar neutrinos events
measured by IceCube (IceCube Collaboration et al., 2018a, b; Eichmann et al.,
2018; Taboada & Stein, 2019; Giommi et al., 2020; Paliya et al., 2020;
Rodrigues et al., 2021). The accelerated protons or nuclei in blazar jets can
interact with the radiation fields, which produces high-energy pions that
eventually decay into photons and neutrinos. However, due to the constraints
on the emission of the electromagnetic cascades by the multi-wavelength
observations, the required power of relativistic protons is significantly
higher than the observed radiative luminosity and even the Eddington
luminosity of the SMBH (e.g. Keivani et al., 2018; Ansoldi et al., 2018; Gao
et al., 2019; Cerruti et al., 2019; Oikonomou et al., 2019; Rodrigues et al.,
2019; Zhang et al., 2020), unless multiple radiation zones along the jet are
considered (Liu et al., 2019; Xue et al., 2019b, 2021)

In addition, if blazars can account for the diffuse neutrinos detected by
IceCube, $\eta_{p,\mathrm{rad}}\sim 3-300$ is required with the lower limit
for the case of a flat CR spectrum (Murase et al., 2014). However, there
exists a large uncertainty about this baryon loading factor so far (Kadler et
al., 2016; Righi et al., 2017; Palladino et al., 2019). In the $p\gamma$ model
of blazars, the energy of relativistic protons is mainly lost through
photohadronic interactions. As the cooling timescale is much longer than the
escape timescale in the AGN jets(Böttcher et al., 2013; Murase et al., 2014;
Gao et al., 2017), the photohadronic interactions are very inefficient (Sikora
et al., 2009; Sikora, 2011), and thus CRs readily escape from the the jets and
host galaxy. It is usually difficult to estimate the baryon loading factor of
AGNs using the radiation of AGNs (Protheroe & Szabo, 1992; Berezhko, 2008).

Galaxy clusters are the largest gravitationally bound objects in the universe.
In the galaxy clusters, CR protons injected by AGN jets are fully confined by
the turbulent magnetic fields in cosmological time (Völk et al. 1996;
Berezinsky et al. 1997, see Brunetti & Jones 2015 for a recent review). Hence,
the confined CR protons can sufficiently interact with intra-cluster medium
(ICM) via $pp$ collision to produce $\gamma$-rays and neutrinos. Meanwhile, as
galaxy clusters are large CR reservoirs which contain CRs injected in the
past, the diffuse $\gamma$-rays in the present can reflect the total injected
energy of CRs in history. Thus, galaxy clusters are favorable environment for
constraining the efficiency of the baryon loading factor. The gamma-ray
emission from individual galaxy cluster was searched for a long time with the
Fermi Large Area Telescope (LAT; (Han et al., 2012; Ackermann et al., 2016),
but so far only a $\gamma$-ray source towards the Coma cluster (Abell 1656)
has been detected by Fermi-LAT (Xi et al., 2018; Adam et al., 2021; Baghmanyan
et al., 2021).

The Coma cluster, located at a distance of $\sim$100 Mpc (z = 0.023), is one
of the nearest and most massive galaxy clusters(Smith et al., 1998; Kubo et
al., 2007). While the mass-to-energy conversion efficiency, denoted by
$\epsilon$, under the standard, radiatively efficient accretion disk model
(Shakura & Sunyaev, 1976) is $\approx 0.1$, a higher efficiency
$\epsilon\approx 0.3$ is also possible for thin-disk accretion on to a extreme
Kerr black hole(Thorne, 1974). Therefore, we adopt a intermediate value of
$W_{\mathrm{g}}\approx 0.2M_{\mathrm{BH,tot}}c^{2}\approx 2.1\times 10^{64}\
\mathrm{erg}$ to estimate the total gravitational energy extracted from black
holes in the Coma cluster, where $M_{\mathrm{BH}}=5.8\times 10^{10}M_{\odot}$
is the total black hole mass(Ensslin et al., 1998). As jets are powered by
accretions of supermassive black holes, the total energy of AGN jets should be
less than or of the same order of magnitude as the releasing gravitational
potential energy. Using the amount of the total gravitation energy, we define
a new baryon loading factor ($\eta_{p,\mathrm{grav}}$) as the efficiency with
which the gravitational energy is converted into relativistic particles in the
present paper.

Multi-wavelength observations and analysis ranging from low frequency radio
wavelengths to $\gamma$-rays centered in the direction of the Coma cluster
have been reported(Kent & Gunn, 1982; Ajello et al., 2009; Brunetti et al.,
2012; Xi et al., 2018; Bonafede et al., 2021). The Coma cluster hosts a giant
radio halo(Giovannini et al., 1993), and extended soft thermal X-ray emission
is observed by the ROSAT all-sky survey(Briel et al., 1992). Using 9 years of
$Fermi$ LAT data, Xi et al. (2018) reported the discovery of $\gamma$-ray
emission from the Coma cluster with an extended spatial structure. The
integral energy flux of $\gamma$-ray emission in the energy range of $0.2-300$
GeV is $\sim 2\times 10^{-12}\ \mathrm{erg\ cm^{-2}\ s^{-1}}$ with a
relatively soft spectral index of $\approx$ -2.7(Xi et al., 2018). The
detection is later confirmed by the $Fermi$-LAT Collaboration(Abdollahi et
al., 2020; Ballet et al., 2020), who found a source in the direction of the
Coma cluster (named as 4FGL J1256.9+2736 in the 4FGL-DR2 catalog), as well as
by other groups (Adam et al., 2021; Baghmanyan et al., 2021). Observations
with VERITAS $\gamma$-ray detector, High Energy Stereoscopic System telescopes
(H.E.S.S.) and other imaging atmospheric Cherenkov telescopes have also
provided upper limits on the very-high-energy (VHE) $\gamma$-ray flux of Coma
cluster (Perkins et al., 2006; Aharonian et al., 2009; Arlen et al., 2012;
Aleksić et al., 2012). We will use these observations to constrain the baryon
loading factor of AGN jets in the Coma cluster.

This paper is organized as follows. In Section 2, we study the propagation of
CRs in the Coma cluster and calculate the radial distribution, taking into
account effects of two different injection histories. The observed
$\gamma$-ray profile and the total flux of the Coma cluster are calculated in
Section 3. The baryon loading factors are given in section 4. We discuss the
results and draw some conclusions in Section 5.

##  2 The radial density profile of CRs in the Coma cluster

Galaxy clusters are able to confine CRs for cosmological timescale (Völk et
al., 1996; Berezinsky et al., 1997). The virial radius of the Coma cluster is
$r_{\mathrm{vir}}\sim 3$ Mpc (Łokas & Mamon, 2003; Kubo et al., 2007). As the
typical scale of magnetic field fluctuations $l_{\mathrm{c}}$ is about $1-10$
% of the virial radius (Brunetti & Jones, 2015), we have $l_{\mathrm{c}}\sim$
0.3 Mpc. The magnetic field strength of the Coma cluster can be derived from
the Faraday rotation measures, $B\sim 5\mu\mathrm{G}$ (Feretti et al., 1995;
Bonafede et al., 2010). If the Larmor radius $r_{\mathrm{L}}$ of a proton or
nucleus with change $Z$ is smaller than the coherence length
($r_{\mathrm{L}}<l_{\mathrm{c}}$), which is satisfied with energy
$E_{\mathrm{c}}\lesssim 10^{21}\ Z(B/5\ \mathrm{\mu G})(l_{\mathrm{c}}/0.3\
\mathrm{Mpc})\ \mathrm{eV},$ the propagation of CRs in the cluster turbulent
magnetic field is in the diffusive regime (Kotera & Lemoine, 2008).

In the diffusive regime, the diffusion coefficient is defined as

$\displaystyle D_{\mathrm{cl}}$ $\displaystyle=$
$\displaystyle\frac{1}{3}\left(\frac{B}{\delta
B}\right)^{2}cr_{\mathrm{L}}^{2-w}l_{\mathrm{c}}^{w-1}$ (1)
$\displaystyle\approx$ $\displaystyle 8.3\times
10^{31}\left(\frac{l_{\mathrm{c}}}{0.3\
\mathrm{Mpc}}\right)^{2/3}\left(\frac{E_{p}}{1\ \mathrm{PeV}}\right)^{1/3}$
$\displaystyle\times\left(\frac{BZ}{5\ \mathrm{\mu G}}\right)^{-1/3}\
\mathrm{cm^{2}s^{-1}},$

where we have assumed $B\sim\delta B$, and $w=5/3$ is the spectral index for
Kolmogorov diffusion. The diffusion timescale of the CR in the Coma cluster is

$\displaystyle t_{\mathrm{diff}}$ $\displaystyle\approx$
$\displaystyle\frac{r_{\mathrm{vir}}^{2}}{2D_{\mathrm{cl}}}\approx
16.4\left(\frac{r_{\mathrm{vir}}}{3\
\mathrm{Mpc}}\right)^{2}\left(\frac{l_{\mathrm{c}}}{0.3\
\mathrm{Mpc}}\right)^{-2/3}$ (2) $\displaystyle\times\left(\frac{E_{p}}{1\
\mathrm{PeV}}\right)^{-1/3}\left(\frac{BZ}{5\ \mathrm{\mu G}}\right)^{1/3}\
\mathrm{Gyr}.$

For CRs with $E_{p}\lesssim 0.4\ \mathrm{PeV}$ the diffusion timescale is
longer than the Hubble time $t_{\mathrm{H}}\sim 14\ \mathrm{Gyr}$, so these CR
protons are confined in the cluster by magnetic fields.

As CRs can not escape the galaxy cluster, the radial distribution of
accumulated CRs in the Coma cluster depends on the injection history of CRs,
which determines the $\gamma$-ray production rate at present. We consider two
reference cases of CR injection history as follows:

  * I

. AGN redshift evolution
The AGN activity changes significantly as a function of redshift, which
increases with redshift up to $z\approx 2$ and then flattens(Hasinger et al.,
2005; Ajello et al., 2012). As the CR injection is expected from the whole
population of AGNs in the Coma cluster, the redshift evolution of AGN
emissivities needs to be taken into consideration. We generally assume that
the injection rate in the Coma cluster follows the AGN luminosity density
injected in the universe as function of redshift (or look-back time). The
luminosity functions of both X-ray-selected (Hasinger et al., 2005) and Fermi-
selected (Ajello et al., 2012) FSRQ samples have the similar profile which
peaks around redshift 1-2, the latter shows a steeper decline after the
redshift peak. In this reference case, we use the luminosity density of Fermi-
selected sample to account for the injection rate, $Q_{p}(t)\propto
L(t_{\mathrm{lb}})$, where $t_{\mathrm{lb}}$ is the look-back time.

  * II.

Constant injection rate
In this case, we consider a simple constant injection rate for a comparison,
$Q_{p}(t)=const.$. The injection time is assumed to last for $\sim 10$ Gyr.

The injected CRs at the sources are assumed to have a power-law spectrum,
$Q_{p}(E_{p},t)\propto E_{p}^{-\alpha}$. For the moment, we assume the
spectrum following the one expected in the strong shocks, i.e., $\alpha\simeq
2$.

Neglecting the energy loss of particles and assuming that the particle
injection occurs at $r=0$, the probability that the CRs reaches radius $r$
after time $t$ can be given by (Aharonian, 2004)

$p(E_{p},r,t)=\frac{e^{-r^{2}/\left(4D_{\mathrm{cl}}(E_{p})t\right)}}{8\pi^{3/2}\left(D_{\mathrm{cl}}(E_{p})t\right)^{3/2}}.$
(3)

Therefore, the radial density distribution of the CRs can be written as

$\frac{dN_{p}(E_{p},r)}{dE_{p}dV}=\int_{0}^{t_{\mathrm{max}}}Q_{p}(E_{p},t)p(E_{p},r,t)dt.$
(4)

Further assuming the total injection energy of protons equals to the total
gravitational energy extracted from the central black holes, we can calculate
the radial density distribution in the Coma cluster. Figure 1 shows the
obtained radial density distribution of the CRs. As lower energy CRs diffuse
more slowly than higher energy ones, the density distribution of 10 TeV CRs is
steeper than that of 1 PeV. The time evolution of the injection rate also
influences the density profile of the CRs. In Figure 2, the contributions of
different time intervals are illustrated separately using different color
scales.

As redshift evolution reaches a peak at the epoch $z\sim 1-2$ (Hasinger et
al., 2005; Ajello et al., 2012), if the injection rate in the cluster follows
the redshift evolution, the contributions of earlier injected CRs, which are
currently diffuse to farther radius, are larger than that in the case of the
constant injection rate. Therefore, for the case that the injection rate
follows the redshift evolution case, the radial density distribution is
flatter. Different redshift evolutions of luminosity density have slight
effect on the distribution compared to the constant case.

Figure 1: The density profile of CRs in the Coma cluster assuming all the
gravitational energy converted into CR energy (i.e.,
$\eta_{p,\mathrm{grav}}=1$) and a flat CR spectrum (i.e., $\alpha=2$). The
constant and redshift evolutions injections of 10 TeV (solid) and 1 PeV
(dashed) CRs are shown in red and blue lines respectively.  Figure 2: The
solid line shows the density profile of the CRs ($E_{p}=10$ TeV) in the Coma
cluster, which is same as Figure 1. The contributions of different starting
and ending times are illustrated separately by dashed lines (with different
color scales). Left: A constant injection rate, time interval $\Delta t=2$
Gyr. Right: The injection rate followed the redshift evolution of the AGN
X-ray luminosity, time interval $\Delta t=t_{\mathrm{lb,max}}/5$.

##  3 $\gamma$-ray production

The confined CRs interact with the ICM in the cluster via the $pp$
interactions. The products of $pp$ interaction are charged and neutral pions
($\pi^{\pm}$ and $\pi^{0}$). Charged pions quickly decay into muons, then into
electrons/positrons and neutrinos. Neutral pions decay into $\gamma$-ray
photons, $\pi^{0}\rightarrow\gamma+\gamma$. In this section, we focus on the
production of $\gamma$-ray photons and calculate the radial profile of the
$\gamma$-ray flux of the Coma cluster.

Following the calculation in Kelner et al. (2006), the $\gamma$-ray emissivity
can be estimated as

$\displaystyle J_{\gamma}(E_{\gamma},{\bf r})$ $\displaystyle=$
$\displaystyle\frac{dN_{\gamma}}{dE_{\gamma}dVdt}$ (5) $\displaystyle=$
$\displaystyle cn_{\mathrm{ICM}}({\bf
r})\int_{E_{\gamma}}^{\infty}\sigma_{pp}(E_{p})\frac{dN_{p}(E_{p},{\bf
r})}{dE_{p}dV}$ $\displaystyle\times
F_{\gamma}(\frac{E_{\gamma}}{E_{p}},E_{p})\frac{dE_{p}}{E_{p}}$

where $\sigma_{pp}(E_{p})$ is the total inelastic cross section of $pp$
interactions, $F_{\gamma}(E_{\gamma}/E_{p},E_{p})$ is the spectrum of the
secondary $\gamma$-ray in a single collision and $n_{\mathrm{ICM}}({\bf r})$
is the density profile of thermal gas in the ICM.

We use the classical “$\beta$ model” of density profile to describe the
thermal gas in the ICM(Cavaliere & Fusco-Femiano, 1976), which is

$n_{\mathrm{ICM}}(r)\approx
n_{\mathrm{ICM}}(0)\left[1+\left(r/r_{\mathrm{c}}\right)^{2}\right]^{-3\beta/2}.$
(6)

The parameters of ICM density profile can be measured from its bremsstrahlung
emission in the X-ray waveband. For the Coma cluster, the parameters in Eq.6
are

$\displaystyle n_{\mathrm{ICM}}(0)=3.42\times 10^{-3}\ \mathrm{cm^{-3}},$
$\displaystyle r_{\mathrm{c}}=290\ \mathrm{kpc}\sim
0.1r_{\mathrm{vir}},\beta=0.75,$ (7)

which are obtained from the X-ray all-sky survey of $Roentgen$ $Satellite$
($ROSAT$) position sensitive proportional counter (PSPC)(Briel et al., 1992).

Figure 3 shows the geometry between the Coma cluster and Earth with a
spherical coordinate system and the original point locates at the center of
the Coma cluster. The luminosity distance of Coma cluster is $R\approx 103$
Mpc, so the distance $D$ between the element volume $dV$ with coordinate
$(r,\theta,\varphi)$ and the Earth is
$D=\sqrt{R^{2}+r^{2}-2Rr\mathrm{sin}\theta\mathrm{cos}\varphi}$, and the
viewing angle is $\mathrm{sin}\psi\approx
r\mathrm{sin}\theta\mathrm{cos}\varphi/D$. Assuming that the CR diffusion and
the ICM density is the spherical symmetry distribution, the $\gamma$-ray
intensity profile average over the solid angle can be calculated as

$\Phi(E_{\gamma},\psi)=\frac{dN_{\gamma}}{dE_{\gamma}d\Omega
dt}=\frac{1}{2\pi(1-\mathrm{cos}\psi_{\mathrm{vir}})}\int\frac{J_{\gamma}(E_{\gamma},{\bf
r})}{4\pi{D}^{2}}dV,$ (8)

where $dV=r^{2}\mathrm{sin}\theta drd\theta d\varphi$,
$2\pi(1-\mathrm{cos}\psi_{\mathrm{vir}})$ is the solid angle of the
$\gamma$-ray flux, $\psi_{\mathrm{vir}}$ is the angular distance of the virial
radius.

Figure 3: The diagrammatic sketch of the geometry between the Coma Cluster and
Earth. The luminosity distance between the Coma Cluster and Earth is $R\approx
103$ Mpc. The distance between the element volume $dV$ with the coordinate
$(r,\theta,\varphi)$ and Earth is $D$, and the viewing angle is $\psi$.

We illustrate the gamma-ray morphologies as a function of the viewing angle
$\psi$ at different photon energies in Figure 4. As shown in Figure 4, the
predicted $\gamma$-ray intensity profile follows both the radial density
distribution of CRs and the gas. In both cases of injection history considered
here, the $\gamma$-ray intensity profiles show flat cores and steep decline at
large radius.

Figure 4: The normalized radial profile of $\gamma$-ray flux in terms of
viewing angle $\psi$ at 10 GeV (red), 100 GeV (blue) and 1 TeV (green). The
flux is normalized by fixing the flux to be unity at 10 GeV at the center of
Coma cluster. The results of constant (dashed line) and redshift evolution
(solid line) injection rates are illustrated separately.

##  4 baryon loading factor in the Coma cluster

###  4.1 time-averaged baryon loading factor

In this section, using the total releasing gravitational potential energy
$W_{\mathrm{g}}$ as denominator, we define the time-averaged baryon loading
factor of the Coma cluster as

$\eta_{p,\mathrm{grav}}\approx\frac{W_{p,\mathrm{tot}}}{W_{\mathrm{g}}},$ (9)

where $W_{p,\mathrm{tot}}$ is the total historical injected energy of the CRs
and $W_{\mathrm{g}}$ is the total releasing gravitational potential energy
from the central black hole, $W_{\mathrm{g}}\approx 0.2M_{\mathrm{BH}}c^{2}$.

The total mass of central black hole of AGNs in the Coma cluster which is
proportional to the total gravitational energy $W_{g}$ in Eq.9 can also
linearly effect the results. Using $M-\sigma$ relationship, the central black
hole mass was suggested to be $\sim 10^{11}M_{\odot}$(Kormendy & Ho, 2013)
given $\sigma\sim 1000\ \mathrm{km\ s^{-1}}$ in the Coma cluster. However, no
such a massive black hole has been found at the cluster center. The central BH
mass of one of the most massive supergiant elliptical galaxy NGC 4889 is
measured to be $M_{\mathrm{BH}}=2.1\times 10^{10}M_{\odot}$(McConnell et al.,
2011). The ratio of the black hole mass to that of the spheroidal component of
the stellar population of the host galaxy is estimated as
$\eta_{\mathrm{BH}}=0.002-0.006$(Kormendy & Richstone, 1995; Wang & Biermann,
1998). In this paper, we use a conservative value $\eta_{\mathrm{BH}}=0.002$,
corresponding to $M_{\mathrm{BH}}\approx 5.8\times 10^{10}M_{\odot}$ (Ensslin
et al., 1998) to give the constraint.

The $\gamma$-ray photons produced via the $pp$ interaction will contribute to
the total gamma-ray flux from the cluster, we thus can use the observed
$\gamma$-ray flux to restrict the total injection energy of protons. The
absorption of gamma-ray photons by extragalactic background light (EBL)
photons is not important for gamma rays below 10 TeV from a source at 100 Mpc.
On the other hand, due to the streaming instability, CRs of energy less than
tens of GeV may be dissipated by self-excited Alfvén waves (Kulsrud & Pearce
1969; Skilling 1971, see Appendix A for detailed calculations in the Coma
cluster), so we use the highest energy bin of $\approx$ 55 GeV in spectrum
energy distributions of the extended emission component of the Coma cluster
for the single radio model from Xi et al. (2018) to restrict the total
injection energy of CRs at 95% confidence level.

Using this upper limit, we estimate the maximum allowed total injected energy
of CR protons above 1 GeV and the corresponding baryon loading
$\eta_{p,\mathrm{grav}}$, for two different injection histories, as well as
the power-law slope of the injection spectrum $\alpha$. Since the expected
pionic gamma-ray flux at a certain energy depends on both these two
parameters, the obtained constraints on the two parameters are coupled. The
relationships between $\eta_{p,\mathrm{grav}}$ and $\alpha$ of protons for two
reference cases are shown separately in Figure 5. The best constraint is
$\eta_{p,\mathrm{grav}}\lesssim 0.02$ for the constant injection case and
$\lesssim 0.04$ for the redshift evolution case when $\alpha=2.1$. For softer
spectral indices, the upper limit increases to $\eta_{p,\mathrm{grav}}\lesssim
0.1$ for the redshift evolution case. In Figure 6 we compare the observed
$\gamma$-ray flux with the theoretical prediction with the parameters
$\eta_{p,\mathrm{grav}}\lesssim 0.04$ (for the case of redshift evolution) and
$\alpha=2.1$. For the sake of illustration, the energy range below the the
highest energy bin of $\approx$ 55 GeV is shown in dotted line and is not used
for comparison, because the CR distribution at the lower energy range could be
modified by the streaming instability.

Note that the streaming instability is less important at higher energy where
the particle number density is lower. A more conservative constraint can be
obtained by using the 99% C.L. integral flux upper limit for $E_{\gamma}>220$
GeV in the Coma cluster core region ($\psi=0.4^{\circ}$) with the 18.6 hr
VERITAS observation (Arlen et al., 2012)111HESS also had a measurement on the
Coma cluster, obtaining a higher flux upper limit (Aharonian et al., 2009),
and hence is less constraining.. Arlen et al. (2012) provide three values of
the upper limit under different assumptions of the gamma-ray index, i.e.,
-2.1, -2.3, and -2.5. We here take the flux upper limit with -2.3 (which is
$\sim 4.7\times 10^{-12}\ \mathrm{ph\ cm^{-2}\ s^{-1}}$) for comparison with
our theoretical prediction, noting that the difference from the given upper
limits for -2.1 and -2.5 is only at the level of about 10%. The obtained
constraint on $\eta_{p,\mathrm{grav}}$ with the VERITAS observation are shown
in Figure 5. For softer spectral indices, the upper limit increases up to
$\eta_{\mathrm{p}}\lesssim 1$ in the redshift evolution case. We also show the
maximum $\gamma$-ray spectra allowed by the VERITAS observation for $\alpha=2$
and $\alpha=2.5$ in Figure 6.

In Fig. 5, we see that the obtained $\eta_{p,\mathrm{grav}}$ with the constant
CR injection rate is systematically smaller by approximately a factor of 2
than that obtained with considering the AGN redshift evolution. This is
because the amount of recently injected cosmic rays is increased in the
constant injection case with respect to the redshift evolution case. The
recently injected CRs have not diffused to the large radius, and hence the
radial distribution of CR density is steeper in the case of constant
injection. Since the gas density at smaller radius is higher, this would
enhance the overall gamma-ray production efficiency and thus leads to a
stricter constraint on the baryon loading factor. Note that a negative
evolution of blazars is observed by Fermi-LAT and X-ray selected sample, which
is mainly in those low-luminosity, high-synchrotron–peaked BL Lac objects
(Beckmann et al., 2003; Ajello et al., 2014). If the negative evolution is
employed in the calculation, the obtained constraint on
$\eta_{p,\mathrm{grav}}$ would become even stricter than that in the constant
injection case, as it would further enhance the amount of recently injected
CRs. However, this class of blazars constitutes only a small fraction of total
blazars and is thus not representative.

In addition to blazar jets, the obtained constraints on the baryon loading
factor may also apply to AGN-driven outflows, which have been suggested as
high-energy particle accelerators as well (e.g. Jiang et al., 2010; Faucher-
Giguère & Quataert, 2012; Wang & Loeb, 2015; Lamastra et al., 2016; Liu et
al., 2018), since the energy of those powerful outflows essentially originate
from the gravitational energy of matter fall into SMBHs.

Figure 5:  The relationship between upper limit of baryon loading factor
$\eta_{p,\mathrm{grav}}$ and power-law index $\alpha$ of injected CRs for the
constant (solid line) and redshift evolution (dashed line) case. The index
range is $2\leq\alpha\leq 2.5$. Results of $Fermi$-LAT (orange) and VERITAS
(blue) are shown separately in different colors.  Figure 6:  Spectrum energy
distribution and upper limits for $Fermi$-LAT (red) and VERITAS (blue)
observations of the Coma cluster. The $Fermi$-LAT data shows the extended
emission component of the Coma cluster for the single radio model with
integration viewing angle $\psi\approx 1.25^{\circ}$, corresponding to the
virial radius of the Coma cluster, and the upper limits are at 95% confidence
level. The 99% confidence level upper limit above an energy threshold of 220
GeV with $\alpha=2.3$ and $\psi\approx 0.4^{\circ}$ is shown in the VERITAS
data. The $\gamma$-ray spectra in our model are illustrated by dashed line in
the figure for comparisons. The best constraint using the $Fermi$-LAT data for
the redshift evolution case is $\eta_{p,\mathrm{grav}}=0.04$ and $\alpha=2.1$
(red, energy range below the highest energy bin is illustrated by dotted
line). Assuming that the injection spectrum index of CRs $\alpha=2$ (dashed)
and $\alpha=2.5$ (dotted), $\gamma$-ray spectra using constraints from VERITAS
observation for the redshift evolution case are also illustrated.

###  4.2 radiation related baryon loading factor

The baryon loading factor defined in Murase et al. (2014) is
$\eta_{p,rad}=L_{p}/L_{\gamma}$, where $L_{p}$ is the total CR luminosity,
$L_{\gamma}$ is corresponding to the bolometric radiation luminosity of the
jet. In this section, we use the total radiation energy $W_{\gamma}$ to obtain
the baryon loading $\eta_{p,rad}$ for comparison.

Assuming that the total radiation energy is proportional to the black hole
mass, we first use the total released gravitational energy to normalize the
integral radiation energy, $\eta_{\gamma}=W_{\gamma}/W_{\mathrm{g}}$.
Therefore, the baryon loading factor becomes to
$\eta_{p,\mathrm{rad}}=W_{\mathrm{p,tot}}/W_{\gamma}=\eta_{p,\mathrm{grav}}/\left\langle\eta_{\gamma}\right\rangle$,
where $\left\langle\eta_{\gamma}\right\rangle$ is the average fraction of the
electromagnetic energy to the gravitational energy from observed blazars.

The sample is taken from Ghisellini et al. (2014), which is composed by 191
FSRQs and 26 BL Lac objects. Because of the lack of measurements of the black
hole mass for BL Lac objects, we only use the FSRQs from the sample to
calculate $\eta_{\gamma}$ (see Appendix B for detailed calculations).

The baryon loading factors with different injection indexes are shown in
Figure 7, $\eta_{p,\mathrm{rad}}\lesssim 1$ ($Fermi$-LAT), $\lesssim 10$
(VERITAS). Note that in the estimation of modified baryon loading factor
$\eta_{p,\mathrm{rad}}$, the same normalization with released gravitational
energy of black hole mass is used in $\eta_{p,\mathrm{grav}}$ and
$\eta_{\gamma}$. In the calculation of
$\left\langle\eta_{\gamma}\right\rangle$, the central black hole mass of each
blazar in our sample is adopted. However, in the calculation of
$\eta_{p,\mathrm{grav}}$ in Section 4.1, as $W_{p,\mathrm{tot}}$ accounts for
the contributions of total CR injection history in the Coma cluster, the
released gravitational energy of black hole mass is estimated by the total
mass of central black holes in the Coma cluster.

Figure 7:  Same as Figure 5, but only
$\eta_{p,\mathrm{rad}}=\eta_{p,\mathrm{grav}}/\left\langle\eta_{\gamma}\right\rangle$
for the redshift evolution case is shown. The dashed line shows results using
the lower 95% confidence range of $\left\langle\eta_{\gamma}\right\rangle$.

##  5 Discussions and Conclusions

Galaxy clusters can effectively confine the CRs in cosmological times, so CRs
can sufficiently interact with ICM to produce $\gamma$-ray and neutrino
radiation. In this paper, taking into account effects of the injection history
of AGN jets, we have studied the propagation and distribution of CRs in the
Coma cluster and obtained constraints on the average baryon loading factor
using the $\gamma$-ray observations. The upper limits of the average baryon
loading factor are $\eta_{p,\mathrm{grav}}\sim 0.01$ and
$\eta_{p,\mathrm{grav}}\sim 0.1$, respectively, from the $Fermi$-LAT and
VERITAS observations for various cosmic-ray power-law indexes. We also use the
integral radiation energy to obtain the upper limits on the conventional
baryon loading factor $\eta_{p,\mathrm{rad}}$, which are
$\eta_{p,\mathrm{rad}}\sim 1$ ($Fermi$-LAT) and $\eta_{p,\mathrm{rad}}\sim 10$
(VERITAS), respectively. If such a constraint can be generalized to all the
AGN in the universe, one may conclude that blazars cannot be the major sources
of the diffuse neutrino background measured by IceCube, when comparing this
upper limit to the theoretically required one (e.g. Murase et al., 2014;
Palladino et al., 2019).

Gamma-ray emission can be also produced by shocks associated with galaxy
merger processes. The Coma cluster may be undergoing the merging
process(Tribble, 1993; Gurzadyan & Mazure, 2001). Shocks associated with
galaxy merger processes would also accelerate particles to relativistic
energies and produce $\gamma$-ray emission(Colafrancesco & Blasi, 1998; Ryu et
al., 2003). In this scenario, the $\gamma$-ray profile is significantly
different from the morphologies predicted by the central injection model in
Figure 4 (Planck Collaboration et al., 2013; Ackermann et al., 2016). Hence
the gamma-ray profile can be used to distinguish the two scenarios for the
gamma-ray emission. Note that, as our main purpose is to give the constraint
on the upper limit of baryonic loading efficiency, the $\gamma$-ray emission
contributed by galaxy merger does not affect our result and can only reduce
the upper limits.

In this paper, we considered the total gravitational energy of the accretion
matter as the upper limit of the energy budget for accelerated protons in AGN
jets or outflows. We note that an analysis of luminous blazars’s spectral
energy distribution based on the conventional one-zone model shows that the
jet/outflow’s power can be larger than the gravitational power of the
accretion matter(Ghisellini et al., 2014). The general relativistic
magnetohydrodynamic (GRMHD) simulations of the radiative-inefficient accretion
flow disks also shows that the jet efficiency is allowed to be greater than
unity for some parameters (McNamara et al., 2000). This may not be
unreasonable: the accretion process can amplify the magnetic field, which
means that the gravitational potential energy of the falling matter can be
stored in the magnetic energy and be released in the certain specific state.
As a result, the instant jet/outflow power may exceed the accretion power. An
observational constraint for the typical AGN phase lifetime is $\sim 10^{5}$
yr(Schawinski et al., 2015), and the typical timescale of the blazar flare
duration is from days to months. Both timescales are much shorter compared to
the total growth time $10^{7}-10^{9}$ yr of AGNs (Fabian & Iwasawa, 1999; Yu &
Tremaine, 2002), or the Salpeter time. It implies that black holes grow via
many such short bursts, so the time-averaged baryon loading factor
$\eta_{p,\mathrm{grav}}$ and $\eta_{p,\mathrm{rad}}$ obtained in this paper
might not apply to one short, specific status of AGNs, such as an intense
flare.

Note that as the population of field galaxies is much larger than that inside
clusters, the employed redshift evolution, which is obtained from the entire
FSRQ population, may not be representative for AGNs in the clusters. Among
many previous studies that are based on Chandra X-ray survey and the Sloan
Digital Sky Survey, some suggested significant suppression of the AGN fraction
in clusters than in the field (Haggard et al., 2010; Hwang et al., 2012;
Haines et al., 2012; Hwang et al., 2012; Ehlert et al., 2013, 2014), but some
did not(Martini et al., 2013; Melnyk et al., 2013; Koulouridis et al., 2014).
Such inconsistency among different studies is probably due to the
observational bias(Xue, 2017). If the AGN inside the cluster has a stronger
redshift evolution than that in the field, the constraint on the baryon
loading factor would be relaxed, while if the AGN redshift evolution is weaker
in the cluster than in the field, the constraint would be stricter.

As the acceleration timescale in AGN jets is much shorter than the cooling
timescales in the energy range in our study, and CRs can escape faster than
they efficiently interact in jets, we use a power-law spectrum to describe the
CRs leaking into the cluster. However, if CRs are still confined in the AGN
jets when the acceleration shuts off, they will loss energy via adiabatic
cooling due to the expansion of the jet. This would reduce the energy of
cosmic rays leaking into the cluster and further increase the upper limit on
the baryon loading factor.

The recent IceCube data shows that the high-energy astrophysical muon-neutrino
spectrum is consistent with a single power law function (Abbasi et al., 2021).
They also test the possible existence of an additional astrophysical component
following certain source-class specific flux prediction. These tests give
independent constraints on the normalization of different theoretical models,
which may be converted to the constraints on the baryon loading factors. For
example, the baryon loading factor $\eta_{p,\mathrm{rad}}$ of AGN inner jets
model from Murase et al. (2014) is $3-300$ and is reduced to about half of the
original value $\sim 1-150$ based on IceCube’s constraint, which is generally
in agreement with our result.

## Acknowledgements

We thank the anonymous referee for constructive comments, and Huirong Yan and
Chong Ge for helpful discussions. This work is supported by NSFC grants
U2031105, 11625312 and 11851304, the National Key R & D program of China under
the grant 2018YFA0404203, and China Manned Spaced Project (CMS-CSST-2021-B11).

##  Appendix A CR streaming in the Coma cluster

CR protons in clusters can stream along magnetic field lines down the CR
gradient. As CR stream in the plasma, a slight anisotropy in the CRs can
generate unstable growth in the Alfvén waves, which is the so-called streaming
instability(Kulsrud & Pearce, 1969; Skilling, 1971). The resulting wave growth
rate can be written as (Kulsrud & Pearce, 1969)

$\Gamma_{\mathrm{cr}}\sim\Omega_{0}\frac{n_{\mathrm{cr}}}{n_{\mathrm{i}}}\frac{v_{\mathrm{D}}}{v_{\mathrm{A}}}$
(A1)

where $\Omega_{0}=eB/(m_{p}c)$ is the cyclotron frequency, $n_{\mathrm{cr}}$
is the number density of CR protons, $n_{i}$ is the ion density in the plasma.
The streaming speed $v_{\mathrm{D}}$ can be written as
$v_{\mathrm{D}}=D(E_{p})\nabla n_{cr}/n_{cr}$, and the Alfvén speed is
$v_{\mathrm{A}}=B/\sqrt{4\pi n_{\mathrm{i}}m_{p}}$.

The wave would slow down CR diffusion and limit streaming speeds to the speed
of the waves. However, there are a variety of damping mechanisms during the
scattering, including nonlinear Landau damping in the collisionless limit and
wave damping by background magneto-hydrodynamic turbulence. The non-linear
Landau damping rate is (Kulsrud, 2005)

$\Gamma_{\mathrm{NL}}\approx
0.3\frac{\Omega}{\mu}\frac{v_{\mathrm{i}}}{c}\left(\frac{\delta
B}{B}\right)^{2},$ (A2)

where $\Omega=\Omega_{0}/\gamma_{p}$, $v_{\mathrm{i}}$ is the velocity of
thermal ions. The turbulence damping rate is (Farmer & Goldreich, 2004)

$\Gamma_{\mathrm{tur}}\sim\frac{v_{\mathrm{A}}}{\sqrt{r_{\mathrm{g}}l_{\mathrm{A}}}},$
(A3)

where $r_{\mathrm{g}}=\gamma_{p}m_{p}c/(eB)$,
$l_{\mathrm{A}}=L_{\mathrm{MHD}}/M_{\mathrm{A}}^{3}$, the injection scale in
the ICM of Coma cluster is $L_{\mathrm{MHD}}=300$ kpc, and the Alfvénic Mach
number $M_{\mathrm{A}}=10$.

Figure 8 shows radius distributions of the growth and damping rates in the ICM
of the Coma cluster for different proton energies. The loading factor
$\eta_{p,\mathrm{grav}}$ = 0.01 and $\alpha=2$ for constant injection case is
adopted in the figure. When the energy of injected protons $E_{p}=50$ GeV, the
growth of resonance wave exceeds the damping rate in the inner region of the
cluster which is also the main region of $\gamma$-ray production. When the
energy of injected protons reaches 500 GeV (the corresponding energy of
product gamma-ray photons is $\sim$50 GeV), the streaming instability of CR
transportation in the ICM is completely suppressed by damping mechanisms, so
we can use the flux upper limit of the highest-energy bin of the Fermi-LAT’s
data to give the constraints.

Figure 8: The growth and damping rate of CR protons in the ICM of Coma
cluster. The density profile of injection protons is constant injection with
power-law index $\alpha=$2 and the total injection energy is limited by Fermi
observations.

##  Appendix B calculation of the radiation energy fraction $\eta_{\gamma}$

In this appendix, we use a sample taken from Ghisellini et al. (2014) of 217
blazars (composed by 191 FSRQs and 26 BL Lac objects) which have been detected
in the $\gamma$–ray band by the $Fermi$-LAT and have been spectroscopically
observed in the optical band. For blazars, using one–zone leptonic model, the
bolometric luminosity $L_{\mathrm{bol}}$ can be established by
multi–wavelength data(Ghisellini & Tavecchio, 2010; Ghisellini et al., 2014).
The bulk Lorentz factor $\Gamma$ is determined by the viewing angle $\psi$,
$\sin\psi\sim 1/\Gamma$. Having the bolometric jet luminosity
$L_{\mathrm{bol}}$ and the bulk Lorentz factor $\Gamma$, the absolute
radiative power $P_{\mathrm{rad}}$ can be calculated as
$P_{\mathrm{rad}}=2fL_{\mathrm{bol}}/\Gamma^{2},$ where the factor 2 accounts
for two jets, $f$ is the numerical factor for external Compton process (4/3)
or synchrotron and the self Compton emission (16/5).

The Salpeter time or e-folding time of an accreting black hole is
$t_{\mathrm{Sal}}\approx 5\times 10^{7}\ \mathrm{yr}$, which gives the
characteristic timescale for the black hole mass to increase by one e-fold, if
quasar’s luminosity equals to the Eddington limit(Salpeter, 1964). If the
characteristic growth time of black holes is comparable to the Salpeter time,
and the growth of black hole mass occurs mainly during the AGN phases, the
integral radiation energy can be estimated as
$W_{\gamma}=P_{\mathrm{rad}}t_{\mathrm{Sal}}$.

The fraction of integral radiation energy to the total released gravitational
energy from black holes is defined as
$\eta_{\gamma}=W_{\gamma}/(0.2M_{\mathrm{BH}}c^{2})$, where $M_{\mathrm{BH}}$
is the central black hole mass. For FSRQs, the central black hole mass can be
estimated through the virial H$\beta$, MgII and CIV broad emission lines by
reverberation mapping(Ghisellini et al., 2014). The histogram of
$\eta_{\gamma}$ from 191 FSRQs in the sample is shown in Figure 9. The
distribution is fitted with a log-normal distribution. The average value
$\left\langle\eta_{\mathrm{\gamma,FSRQ}}\right\rangle\approx 0.03$, with a 95%
confidence interval $\eta_{\mathrm{\gamma,FSRQ}}\in[0.007,1.5]$. As the BL Lac
objects have no central black hole mass estimation because of the lack of high
S/N broad-line measurements, we use the average black hole mass of the FSRQs
to give a simple estimation of 26 BL Lac objects in the sample,
$\mathrm{log}_{10}(M_{\mathrm{BH,avg}}/M_{\odot})\approx 8.5$, then
$\eta_{\mathrm{\gamma,BLL}}\approx 0.03$, which is in good agreement with that
obtained with the FSRQs.

Figure 9:  The histogram of $\eta_{\gamma}$ for the 191 blazars in our sample.
The distribution (orange solid line) is fitted with a log–normal Gaussian
distribution with average value
$\left\langle\mathrm{log}_{10}(\eta_{\gamma})\right\rangle\approx-1.48$
($\left\langle\eta_{\gamma}\approx 0.03\right\rangle$) and width
$\sigma=0.67$.

## References

  * Abbasi et al. (2021) Abbasi, R., Ackermann, M., Adams, J., et al. 2021, arXiv:2111.10299
  * 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. (2016) Ackermann, M., Ajello, M., Albert, A., et al. 2016, ApJ, 819, 149. doi:10.3847/0004-637X/819/2/149
  * Adam et al. (2021) Adam, R., Goksu, H., Brown, S., et al. 2021, A&A, 648, A60. doi:10.1051/0004-6361/202039660
  * Aharonian (2004) Aharonian, F. A. 2004, Very High Energy Cosmic Gamma Radiation: A Crucial Window on the Extreme Universe. Edited by AHARONIAN FELIX A. Published by World Scientific Publishing Co. Pte. Ltd., . ISBN #9789812561732. doi:10.1142/4657
  * Aharonian et al. (2009) Aharonian, F., Akhperjanian, A. G., Anton, G., et al. 2009, A&A, 502, 437
  * Ajello et al. (2009) Ajello, M., Rebusco, P., Cappelluti, N., et al. 2009, ApJ, 690, 367. doi:10.1088/0004-637X/690/1/367
  * Ajello et al. (2012) Ajello, M., Shaw, M. S., Romani, R. W., et al. 2012, ApJ, 751, 108. doi:10.1088/0004-637X/751/2/108
  * Ajello et al. (2014) Ajello, M., Romani, R. W., Gasparrini, D., et al. 2014, ApJ, 780, 73. doi:10.1088/0004-637X/780/1/73
  * Aleksić et al. (2012) Aleksić, J., Alvarez, E. A., Antonelli, L. A., et al. 2012, A&A, 541, A99. doi:10.1051/0004-6361/201118502
  * Ansoldi et al. (2018) Ansoldi, S., Antonelli, L. A., Arcaro, C., et al. 2018, ApJ, 863, L10. doi:10.3847/2041-8213/aad083
  * Arlen et al. (2012) Arlen, T., Aune, T., Beilicke, M., et al. 2012, ApJ, 757, 123. doi:10.1088/0004-637X/757/2/123
  * Atoyan & Dermer (2001) Atoyan, A. & Dermer, C. D. 2001, Phys. Rev. Lett., 87, 221102. doi:10.1103/PhysRevLett.87.221102
  * Baghmanyan et al. (2021) Baghmanyan, V., Zargaryan, D., Aharonian, F., et al. 2021, arXiv:2110.00309
  * Ballet et al. (2020) Ballet, J., Burnett, T. H., Digel, S. W., et al. 2020, arXiv:2005.11208
  * Beckmann et al. (2003) Beckmann, V., Engels, D., Bade, N., et al. 2003, A&A, 401, 927. doi:10.1051/0004-6361:20030184
  * Berezhko (2008) Berezhko, E. G. 2008, ApJ, 684, L69. doi:10.1086/592233
  * Berezinsky et al. (1997) Berezinsky, V. S., Blasi, P., & Ptuskin, V. S. 1997, ApJ, 487, 529. doi:10.1086/304622
  * Berezinsky et al. (2006) Berezinsky, V., Gazizov, A., & Grigorieva, S. 2006, Phys. Rev. D, 74, 043005. doi:10.1103/PhysRevD.74.043005
  * Biermann (1988) Biermann, P. L. 1988, High Energy Astrophysics, 134
  * Bonafede et al. (2010) Bonafede, A., Feretti, L., Murgia, M., et al. 2010, A&A, 513, A30. doi:10.1051/0004-6361/200913696
  * Bonafede et al. (2021) Bonafede, A., Brunetti, G., Vazza, F., et al. 2021, ApJ, 907, 32. doi:10.3847/1538-4357/abcb8f
  * Böttcher et al. (2013) Böttcher, M., Reimer, A., Sweeney, K., et al. 2013, ApJ, 768, 54. doi:10.1088/0004-637X/768/1/54
  * Briel et al. (1992) Briel, U. G., Henry, J. P., & Boehringer, H. 1992, A&A, 259, L31
  * Brunetti et al. (2012) Brunetti, G., Blasi, P., Reimer, O., et al. 2012, MNRAS, 426, 956. doi:10.1111/j.1365-2966.2012.21785.x
  * Brunetti & Jones (2015) Brunetti, G. & Jones, T. W. 2015, Magnetic Fields in Diffuse Media, 557
  * Cavaliere & Fusco-Femiano (1976) Cavaliere, A., & Fusco-Femiano, R. 1976, A&A, 500, 95
  * Cerruti et al. (2019) Cerruti, M., Zech, A., Boisson, C., et al. 2019, MNRAS, 483, L12. doi:10.1093/mnrasl/sly210
  * Colafrancesco & Blasi (1998) Colafrancesco, S. & Blasi, P. 1998, Astroparticle Physics, 9, 227. doi:10.1016/S0927-6505(98)00018-8
  * Dermer et al. (2009) Dermer, C. D., Razzaque, S., Finke, J. D., et al. 2009, New Journal of Physics, 11, 065016. doi:10.1088/1367-2630/11/6/065016
  * Eichmann et al. (2018) Eichmann, B., Rachen, J. P., Merten, L., et al. 2018, J. Cosmology Astropart. Phys, 2018, 036. doi:10.1088/1475-7516/2018/02/036
  * Ehlert et al. (2013) Ehlert, S., Allen, S. W., Brandt, W. N., et al. 2013, MNRAS, 428, 3509. doi:10.1093/mnras/sts288
  * Ehlert et al. (2014) Ehlert, S., von der Linden, A., Allen, S. W., et al. 2014, MNRAS, 437, 1942. doi:10.1093/mnras/stt2025
  * Ensslin et al. (1998) Ensslin, T. A., Wang, Y., Nath, B. B., et al. 1998, A&A, 333, L47
  * Fabian & Iwasawa (1999) Fabian, A. C. & Iwasawa, K. 1999, MNRAS, 303, L34. doi:10.1046/j.1365-8711.1999.02404.x
  * Farmer & Goldreich (2004) Farmer, A. J. & Goldreich, P. 2004, ApJ, 604, 671. doi:10.1086/382040
  * Faucher-Giguère & Quataert (2012) Faucher-Giguère, C.-A. & Quataert, E. 2012, MNRAS, 425, 605. doi:10.1111/j.1365-2966.2012.21512.x
  * Feretti et al. (1995) Feretti, L., Dallacasa, D., Giovannini, G., et al. 1995, A&A, 302, 680
  * Gao et al. (2017) Gao, S., Pohl, M., & Winter, W. 2017, ApJ, 843, 109. doi:10.3847/1538-4357/aa7754
  * Gao et al. (2019) Gao, S., Fedynitch, A., Winter, W., et al. 2019, Nature Astronomy, 3, 88. doi:10.1038/s41550-018-0610-1
  * Giommi et al. (2020) Giommi, P., Padovani, P., Oikonomou, F., et al. 2020, A&A, 640, L4. doi:10.1051/0004-6361/202038423
  * Ghisellini & Tavecchio (2010) Ghisellini, G. & Tavecchio, F. 2010, MNRAS, 409, L79. doi:10.1111/j.1745-3933.2010.00952.x
  * Ghisellini et al. (2014) Ghisellini, G., Tavecchio, F., Maraschi, L., et al. 2014, Nature, 515, 376. doi:10.1038/nature13856
  * Giovannini et al. (1993) Giovannini, G., Feretti, L., Venturi, T., et al. 1993, ApJ, 406, 399. doi:10.1086/172451
  * Gurzadyan & Mazure (2001) Gurzadyan, V. G. & Mazure, A. 2001, New A, 6, 43. doi:10.1016/S1384-1076(01)00043-4
  * Haines et al. (2012) Haines, C. P., Pereira, M. J., Sanderson, A. J. R., et al. 2012, ApJ, 754, 97. doi:10.1088/0004-637X/754/2/97
  * Han et al. (2012) Han, J., Frenk, C. S., Eke, V. R., et al. 2012, MNRAS, 427, 1651. doi:10.1111/j.1365-2966.2012.22080.x
  * Haggard et al. (2010) Haggard, D., Green, P. J., Anderson, S. F., et al. 2010, ApJ, 723, 1447. doi:10.1088/0004-637X/723/2/1447
  * Hasinger et al. (2005) Hasinger, G., Miyaji, T., & Schmidt, M. 2005, A&A, 441, 417
  * Hwang et al. (2012) Hwang, H. S., Park, C., Elbaz, D., et al. 2012, A&A, 538, A15. doi:10.1051/0004-6361/201117351
  * IceCube Collaboration et al. (2018a) IceCube Collaboration, Aartsen, M. G., Ackermann, M., et al. 2018, Science, 361, eaat1378. doi:10.1126/science.aat1378
  * IceCube Collaboration et al. (2018b) IceCube Collaboration, Aartsen, M. G., Ackermann, M., et al. 2018, Science, 361, 147. doi:10.1126/science.aat2890
  * Jiang et al. (2010) Jiang, Y.-F., Ciotti, L., Ostriker, J. P., et al. 2010, ApJ, 711, 125. doi:10.1088/0004-637X/711/1/125
  * Kadler et al. (2016) Kadler, M., Krauß, F., Mannheim, K., et al. 2016, Nature Physics, 12, 807. doi:10.1038/nphys3715
  * Kelner et al. (2006) Kelner, S. R., Aharonian, F. A., & Bugayov, V. V. 2006, Phys. Rev. D, 74, 034018. doi:10.1103/PhysRevD.74.034018
  * Kent & Gunn (1982) Kent, S. M. & Gunn, J. E. 1982, AJ, 87, 945. doi:10.1086/113178
  * Kimura et al. (2015) Kimura, S. S., Murase, K., & Toma, K. 2015, ApJ, 806, 159. doi:10.1088/0004-637X/806/2/159
  * Kubo et al. (2007) Kubo, J. M., Stebbins, A., Annis, J., et al. 2007, ApJ, 671, 1466. doi:10.1086/523101
  * McConnell et al. (2011) McConnell, N. J., Ma, C.-P., Gebhardt, K., et al. 2011, Nature, 480, 215. doi:10.1038/nature10636
  * Keivani et al. (2018) Keivani, A., Murase, K., Petropoulou, M., et al. 2018, ApJ, 864, 84. doi:10.3847/1538-4357/aad59a
  * Kormendy & Richstone (1995) Kormendy, J. & Richstone, D. 1995, ARA&A, 33, 581. doi:10.1146/annurev.aa.33.090195.003053
  * Kormendy & Ho (2013) Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511
  * Kotera & Lemoine (2008) Kotera, K. & Lemoine, M. 2008, Phys. Rev. D, 77, 023005. doi:10.1103/PhysRevD.77.023005
  * Koulouridis et al. (2014) Koulouridis, E., Plionis, M., Melnyk, O., et al. 2014, A&A, 567, A83. doi:10.1051/0004-6361/201423601
  * Kulsrud & Pearce (1969) Kulsrud, R. & Pearce, W. P. 1969, ApJ, 156, 445. doi:10.1086/149981
  * Kulsrud (2005) Kulsrud, R. M. 2005, Plasma physics for astrophysics, by R.M. Kulsrud Princeton series in astrophysics. Princeton, NJ: Princeton University Press, 2005
  * Lamastra et al. (2016) Lamastra, A., Fiore, F., Guetta, D., et al. 2016, A&A, 596, A68. doi:10.1051/0004-6361/201628667
  * Liu et al. (2018) Liu, R.-Y., Murase, K., Inoue, S., et al. 2018, ApJ, 858, 9. doi:10.3847/1538-4357/aaba74
  * Liu et al. (2019) Liu, R.-Y., Wang, K., Xue, R., et al. 2019, Phys. Rev. D, 99, 063008. doi:10.1103/PhysRevD.99.063008
  * Łokas & Mamon (2003) Łokas, E. L. & Mamon, G. A. 2003, MNRAS, 343, 401. doi:10.1046/j.1365-8711.2003.06684.x
  * Mannheim et al. (1992) Mannheim, K., Stanev, T., & Biermann, P. L. 1992, A&A, 260, L1
  * Martini et al. (2013) Martini, P., Miller, E. D., Brodwin, M., et al. 2013, ApJ, 768, 1. doi:10.1088/0004-637X/768/1/1
  * McNamara et al. (2000) McNamara, B. R., Wise, M., Nulsen, P. E. J., et al. 2000, ApJ, 534, L135. doi:10.1086/312662
  * Melnyk et al. (2013) Melnyk, O., Plionis, M., Elyiv, A., et al. 2013, A&A, 557, A81. doi:10.1051/0004-6361/201220624
  * Murase et al. (2014) Murase, K., Inoue, Y., & Dermer, C. D. 2014, Phys. Rev. D, 90, 023007. doi:10.1103/PhysRevD.90.023007
  * Oikonomou et al. (2019) Oikonomou, F., Murase, K., Padovani, P., et al. 2019, MNRAS, 489, 4347. doi:10.1093/mnras/stz2246
  * Paliya et al. (2020) Paliya, V. S., Böttcher, M., Olmo-García, A., et al. 2020, ApJ, 902, 29. doi:10.3847/1538-4357/abb46e
  * Palladino et al. (2019) Palladino, A., Rodrigues, X., Gao, S., et al. 2019, ApJ, 871, 41. doi:10.3847/1538-4357/aaf507
  * Perkins et al. (2006) Perkins, J. S., Badran, H. M., Blaylock, G., et al. 2006, ApJ, 644, 148. doi:10.1086/503321
  * Planck Collaboration et al. (2013) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2013, A&A, 554, A140. doi:10.1051/0004-6361/201220247
  * Protheroe & Szabo (1992) Protheroe, R. J. & Szabo, A. P. 1992, Phys. Rev. Lett., 69, 2885. doi:10.1103/PhysRevLett.69.2885
  * Rachen & Biermann (1993) Rachen, J. P. & Biermann, P. L. 1993, A&A, 272, 161
  * Righi et al. (2017) Righi, C., Tavecchio, F., & Guetta, D. 2017, A&A, 598, A36. doi:10.1051/0004-6361/201629412 ’
  * Rodrigues et al. (2019) Rodrigues, X., Gao, S., Fedynitch, A., et al. 2019, ApJ, 874, L29. doi:10.3847/2041-8213/ab1267
  * Rodrigues et al. (2021) Rodrigues, X., Heinze, J., Palladino, A., et al. 2021, Phys. Rev. Lett., 126, 191101. doi:10.1103/PhysRevLett.126.191101
  * Ryu et al. (2003) Ryu, D., Kang, H., Hallman, E., et al. 2003, ApJ, 593, 599. doi:10.1086/376723
  * Salpeter (1964) Salpeter, E. E. 1964, ApJ, 140, 796. doi:10.1086/147973
  * Schawinski et al. (2015) Schawinski, K., Koss, M., Berney, S., et al. 2015, MNRAS, 451, 2517. doi:10.1093/mnras/stv1136
  * Shakura & Sunyaev (1976) Shakura, N. I. & Sunyaev, R. A. 1976, MNRAS, 175, 613. doi:10.1093/mnras/175.3.613
  * Sikora et al. (2009) Sikora, M., Stawarz, Ł., Moderski, R., et al. 2009, ApJ, 704, 38. doi:10.1088/0004-637X/704/1/38
  * Sikora (2011) Sikora, M. 2011, Jets at All Scales, 275, 59. doi:10.1017/S1743921310015644
  * Smith et al. (1998) Smith, D. A., Herter, T., & Haynes, M. P. 1998, ApJ, 494, 150. doi:10.1086/305175
  * Skilling (1971) Skilling, J. 1971, ApJ, 170, 265. doi:10.1086/151210
  * Taboada & Stein (2019) Taboada, I. & Stein, R. 2019, The Astronomer’s Telegram, 12967
  * Takahara (1990) Takahara, F. 1990, Progress of Theoretical Physics, 83, 1071. doi:10.1143/PTP.83.1071
  * Thorne (1974) Thorne, K. S. 1974, ApJ, 191, 507. doi:10.1086/152991
  * Tribble (1993) Tribble, P. C. 1993, MNRAS, 263, 31. doi:10.1093/mnras/263.1.31
  * Völk et al. (1996) Völk, H. J., Aharonian, F. A., & Breitschwerdt, D. 1996, TeV Gamma-ray Astrophysics. Theory and Observations, 279
  * Wang & Biermann (1998) Wang, Y. & Biermann, P. L. 1998, A&A, 334, 87
  * Wang & Loeb (2015) Wang, X. & Loeb, A. 2015, MNRAS, 453, 837. doi:10.1093/mnras/stv1649
  * Xi et al. (2018) Xi, S.-Q., Wang, X.-Y., Liang, Y.-F., et al. 2018, Phys. Rev. D, 98, 063006. doi:10.1103/PhysRevD.98.063006
  * Xue et al. (2019a) Xue, R., Liu, R.-Y., Wang, X.-Y., et al. 2019, ApJ, 871, 81. doi:10.3847/1538-4357/aaf720
  * Xue et al. (2019b) Xue, R., Liu, R.-Y., Petropoulou, M., et al. 2019, ApJ, 886, 23. doi:10.3847/1538-4357/ab4b44
  * Xue et al. (2021) Xue, R., Liu, R.-Y., Wang, Z.-R., et al. 2021, ApJ, 906, 51. doi:10.3847/1538-4357/abc886
  * Xue (2017) Xue, Y. Q. 2017, New A Rev., 79, 59. doi:10.1016/j.newar.2017.09.002
  * Yu & Tremaine (2002) Yu, Q. & Tremaine, S. 2002, MNRAS, 335, 965. doi:10.1046/j.1365-8711.2002.05532.x
  * Zhang & Li (2017) Zhang, B. T. & Li, Z. 2017, J. Cosmology Astropart. Phys, 2017, 024. doi:10.1088/1475-7516/2017/03/024
  * Zhang et al. (2020) Zhang, B. T., Petropoulou, M., Murase, K., et al. 2020, ApJ, 889, 118. doi:10.3847/1538-4357/ab659a


#  Gradient vector fields of codimension one on the 2-sphere with at most ten
singular points

Svitlana Bilun, Bohdana Hladysh, Alexandr Prishlyak and Vladislav Sinitsyn

###### Abstract

We describe all possible (469) topological structures of сodimension one
gradient vector fields on the 2-sphere with at most ten singular points. To
describe structures, we use a graph whose edges are one-dimensional stable
manifolds. The saddle-node singularity is specified by selecting a pair of
vertices-edge or edge-face, and the saddle connection is specified by a
T-vertex.

##  1 Introduction

In this paper, we consider gradient vector fields on a sphere. Since the
function increases along each trajectory, the field has no cycles and
polycycles. In general position, a typical gradient field (a vector field of
the codimension 0) is a Morse field (Morse-Smale field without closed
trajectories). In typical one-parameter families of gradient vector fields,
two types of bifurcations are possible: a saddle-node and a saddle connection.
The corresponding vector fields at the time of the bifurcation are fields of
codimension 1. In our case, they completely determine the topological type of
the bifurcation. To classify Morse fields, a cell complex structure (diagram)
is often used, in which cells of dimension n are stable manifolds of singular
points with Morse index equal to n. We apply this approach to the
classification of vector fields of codimension 1 on the 2-sphere. The cell
structure on the 2-sphere is determinated by the imbedded graph, which is the
1-skeleton.

Without loss of generality, we assume that the number of singular points does
not increase under bifurcation (as the parameter increases). The saddle-node
bifurcation is defined by a pair of cells corresponding to the singular points
participating in the bifurcation. We mark this pair on the diagram with a
green arrow or a green triangle. A saddle-node bifurcation in the diagram
corresponds to a point of degree 3 (T-vertex), where two edges (half-edges)
are opposite and the third is perpendicular to them.

Then, the separatrix that connects the saddle with the node (source or sink)
contracts to a point under the saddle-node bifurcation.

Selection of pairs of cells of neighboring dimensions is used in the discrete
Morse theory [8].

Topological invariants of functions were constructed in the papers of Kronrod
[17] and Reeb [41] for oriented maniofolds, in [20] for non-orientable two-
dimensional manifolds and in [7, 15, 16] for manifolds with boundary, in [35]
for non-compact manifolds.

In general, Morse vector fields (Morse-Smale vector fields without closed
orbits) are gradient field of Morse functions. If we fix the value of
functions in singular points the field determinate the topological structure
of the function [20, 43]. Therefore, Morse–Smale vector fields classification
is closely related to the classification of the functions.

Topological classification of smooth function on closed 2-manifolds was also
investigated in [6, 16, 15, 35, 34, 39, 20, 24, 50, 47, 2, 42, 44], on
2-manifolds with the boundary in [13, 16, 14] and on closed 3-manifolds in
[33].

In [18, 21, 22, 23, 40, 1, 29, 32, 45, 49, 30, 31, 38, 28], the
classifications of flows on closed 2- manifolds and [5, 19, 32, 29, 25, 37,
38] on manifolds with the boundary were obtained. Topological properties of
Morse-Smale vector fields on 3-manifolds was investigated in [46, 48, 36, 24,
49, 26, 27, 11, 4, 3].

The topological invariants of graphs and their embeddings in the 2-manifolds
can be founded in [23, 10, 12, 9].

In view of the bijection between vector fields and flows on closed manifolds,
we identify these two concepts.

The purpose of this paper is to describe all possible topological structures
of the gradient flow structure of dimension 1 on a sphere with no more than
four saddles (a saddle-node point is also considered a saddle).

##  2 Typical one-parameter bifurcations of gradient flows on a sphere

Typical vector fields on compact 2-manifolds are Morse-Smale fields. Among the
gradient fields are Morse fields or Morse-Smale gradient-like fields that do
not contain closed trajectories. They satisfy three properties:

1) singular points are nondegenerate;

2) there are no saddle connections;

3) $\alpha$-limit ($\omega$-limit) set of each trajectory is a singular point.

In typical one-parameter field families, one of these conditions is violated.
Violation of the first condition leads to a saddle-node bifurcation, and the
second to the appearance of a saddle connection. The third condition cannot be
violated, because according to the Poincaré-Bendixon theorem, the
$\alpha$-limit ($\omega$-limit) set of every trajectory on the sphere is
either a singular point, or a cycle, or a polycycle. Since gradient fields
have no cycles and polycycles, this set is a singular point.

###  2.1 Saddle-node

The saddle-node bifurcationis shown in Fig. 1 if the node is the source.

Figure 1: saddle-nod bifurcation

It can be described by the equation $V(x,y,a)=\\{x,y^{2}+a\\}$ if the node is
the source and the equation $V(x,y,a)=\\{-x,-y^{2}a\\}$ if the node is a sink.
Here, $a$ is a parameter. If $a<0$ we get the flow before the bifurcation, if
$a>0$ we get the flow after the bifurcation, and if $a=0$ we get the flow at
the moment of the bifurcation (the flow of codimensionality 1). The saddle-
sink bifurcation can be obtained from the saddle-source rotation of the
direction on all trajectories, that is, by taking the reversed flow. Note that
the saddle-node bifurcation corresponds to a typical one-parameter-method
catastrophe (deformation) of functions: $f(x,y,a)=x^{2}+y^{3}+ay$.

Consider a stable manifold (red separatrix in the figures), for a singular
saddle-source point. It cannot start at this particular point because the
gradient field contains no loops. This means that a stable saddle manifold
participating in a saddle-source bifurcation does not form a loop. Therefore,
in order to determine the saddle-source bifurcation on the cellular partition
given by the stable multiples of the Morse field before the bifurcation, it is
necessary to select a 1-cell that is not a loop and one of its ends. In the
following figures, we will do this with the help of a blue arrow or blue
triangles at the ends of the ribs.

Similarly to the above, a saddle-sink bifurcation is possible if the unstable
manifold of the saddle before the bifurcation does not form a loop. This means
that the corresponding 1-cell is contained in the boundary of two different
2-cells. To set such a bifurcation, you need to select a pair: 1 cell and one
of the two adjacent 2-cells. In the pictures, we will do this with the help of
a blue arrow (or a blue triangle) directed from the 1-cell towards the
corresponding 2-cell.

Note that the bifurcation of the saddle flow in Morse’s theory corresponds to
the operation of reducing critical points or reducing complementary handles.

###  2.2 Saddle connection

The bifurcation of the saddle connection is shown in Fig. 2. It can be
described by the equation $V(x,y,a)=\\{x^{2}-y^{2}-1,-2xy+a\\}$.

Figure 2: saddle connection

We describe one of the possible situations when it appears in typical one-
parameter families of gradient fields. Let $p,q$ be saddle critical points of
the function $f$, $f(p)<f(q)$,
$[\frac{f(p)+f(q)}{2}-\varepsilon,\frac{f(p)+f(q)}{2}+\varepsilon]$ does not
contain critical values of $f$, $L$ of the line component of level
$f^{(}-1)(\frac{f(p)+f(q)}{2})$, $u$ is the separatrix of $p$, $v$ is the
separatrix of $q$, which intersect $L$. Then according to Morse theory,
$f^{-1}([\frac{f(p)+f(q)}{2}-\varepsilon,\frac{f(p)+f(q)}{2}+\varepsilon])$ is
homeomorphic to the cylinder $S^{1}\times[0,1]$ and in the coordinates $(s,t)$
($s$–polar angle on the circle) on it, the trajectories are given by segments
$s=\text{const}$. Consider the twisting of a cylinder given by the formula
$(s,t)\to(s+a,t)$ ($a$– deformation parameter). At the same time, the
trajectories turn into helical lines. With a continuous change of the
parameter $a$, the point of intersection of the trajectory $u$ with the circle
$frac{f(p)+f(q)}{2}+\varepsilon$ continuously rotate along it. The
intersection point of the trajectory $v$ with this circle remain unchanged.
According to the intermediate value theorem, there is a moment in time $a$ at
which these points coincide.Then the bifurcation of the saddle connection
occur at this point of time $a$. Note that the construction described above
can be implemented by changing the Riemannian metric (changing angles
considered perpendicular), without changing the function itself.

In the figuress, the bifurcation of the saddle connection will be indicated by
a rib that is perpendicularly glued to the middle of the other rib (in the
form of the letter T). A glued edge (1-cell) containing a separatrix
connection be called perpendicular. Therefore, a graph containing a T-vertex
and a perpendicular edge corresponds to a flow of codimensionality 1.

In Morse’s theory, bifurcation of the saddle connection corresponds to the
operation of adding handles (sliding one handle over another handle of the
same index).

###  2.3 Distinguishing graph

Let’s construct a graph embedded in a sphere, the vertices of which are
sources, saddle-node points, and T-vertices, and the edges are stable one-
dimensional manifolds. On the graph, we select a T-vertex, if it is present
and has a perpendicular edge in it. And in the case of a saddle-source, we
select the corresponding vertex-edge pair, and in the case of a saddle-sink,
we select an edge-face pair. A spherical graph with a selected element
(T-vertex and edge, vertex-edge or edge-face pair) is called a distinguishing
graph of a vector field or simply a field graph. Two such graphs are
equivalent if there is an isomorphism of them as spherical graphs that
preserves the selected element.

###### Theorem 1

Two vector fields of codimensionality 1 on the sphere are topologically
equivalent if and only if their distinguishing graphs are equivalent. Each
spherical graph with a selected pair is a distinguishing graph of a vector
field of codimensional 1.

Proof. The necessity follows directly from the construction of the
distinguishing graph. To prove sufficiency, let’s construct a separatrix
diagram for the graph. To do this, consider a dual graph. If there is no
T-vertex, then the selected pair defines a half-edge on the initial or dual
graph. Let’s squeeze it into a point. The remaining half-edges are flow
separatrixes. Their direction is determined from the vertices of the initial
graph to the vertices of the dual graph. So, we got a separatrix diagram of a
vector field, which specifies its topological type. In the case of a T-vertex,
for the perpendicular edge, we still construct a double edge, and for the
other two edges incident to the T-vertex, we construct an edge that is a
continuation of the perpendicular edge. Therefore, the T-vertex turn into a
vertex of valence 4 (saddle). Further considerations are the same as in the
case without a T-vertex.

Since the separatrix diagram of a graph is defined uniquely, such a graph
defines a vector field.

##  3 Fields with one and two saddles.

Let us describe the structure of all possible vector fields of
codimensionality 1 with one and two saddles. For saddle-node bifurcations, we
define the flows to the bifurcations by their graphs. The number of ribs is
equal to the number of saddles. There are two graphs with one edge: 1)
$G^{1}_{1}$ – with two vertices (segment – Fig. 3 1)), 2) $G^{1}_{2}$ – with
one vertex (loop – Fig.3. 7) ).

Figure 3: Flows of codimensionality 1 with one and two saddles. 1)-6) saddle-
source, 7)-12) saddle-sink, 13)-16) saddle connection. Graphs 1) $G^{1}_{1}$,
2) $G^{2}_{1}$, 3) $G^{2}_{2}$, 5) $G^{2}_{3}$, 7) $G^{1}_{2}$ 9) $G^{2}_{4}$.

Only four graphs are possible with two edges: 1) $G^{2}_{1}$ – a pair of
multiple edges, 2) $G^{2}_{2}$ – a chain of length 2, 3) $G^{2}_{3}$ – a
segment with a loop, 4) $G^{2}_{4}$ – a pair of loops.

On $G^{1}_{1}$ there is a single (with accuracy to homeomorphism) saddle-
source bifurcation (see Fig.3. 2). Saddle-sink is not possible because
$G^{1}_{1}$ does not split the plane. On $G^{1}_{2}$, on the contrary, there
is a saddle-sink bifurcation, and a saddle-source bifurcation does not exist
because the only edge is a loop. Note that $G^{1}_{2}$ is dual to $G^{1}_{1}$,
the corresponding vector fields are obtained from each other by multiplication
by $-1$. When the direction of the vector field changes, the source turns into
a sink, so the saddle-source bifurcation turns into a saddle-sink. Bifurcation
of the saddle connection for flows with one saddle is not possible, since two
saddles take part in it.

On $G^{2}_{1}$, all possible saddle-source pairs are symmetric, so there is
only one such bifurcation 2). Given the symmetry of $G^{2}_{2}$, there are two
saddle-node bifurcations on it: 3) – with an extreme vertex, 4) – with a
central vertex. On the graph $G^{2}_{3}$, the loop cannot participate in the
saddle-source bifurcation. For another segment edge, there are two options for
choosing its end. So, we have two saddle-source bifurcations 8) and 9).
Saddle-source bifurcation is not possible on a pair of loops. If we consider
the reverse flows, then the resulting 5 saddle-source bifurcations will turn
into saddle-sink bifurcations on the dual graphs shown in Fig.3. 8)-12).

If both saddles participate in the bifurcation of the saddle connection, then
one of them will be normal, and the other - perpendicular. If the normal edge
is a segment, then there are two options for the perpendicular edge: 13) and
14). If the regular edge is a loop, then we also have two options 15) and 16).

We went through all the possible options, and therefore it is fair

###### Theorem 2

The sphere has the following structures of codimensionality 1 flows:

  * •

two flow structures with three singular points (1 and 7);

  * •

there are no streams with four points;

  * •

ten structures with five points (2–6 and 8–12);

  * •

four structures with six points (13–16).

##  4 Fields with three saddles

In the following, we will describe only saddle-source bifurcations, since they
specify all possible saddle-sink bifurcations on dual graphs. Let us describe
all possible graphs with three edges, see Fig.4.

Figure 4: Fields of codimensionality 1 with three saddles: the saddle-source
case  Figure 5: Flows of codimensionality 1 with three saddles: a peculiarity
of the saddle connection

Let a graph contain four vertices, then it is a tree. There can be only two
such trees: the chain $G^{3}_{1}$ and the Y-graph $G^{3}_{1}$. If the graph
has three vertices, then they are connected by a chain. Let’s number them so
that the second vertex is the middle one on this chain. Let us first consider
the case when the third edge is not incident to the second vertex. If it
connects the first and third vertices, then we get a triangle $G^{3}_{3}$, and
if it has one vertex at its ends, then the graph $G^{3}_{4}$ is two chains
with a loop at the end. If the third edge starts from vertex 2, then we have
two options: 1) it ends at the extreme vertex – graph $G^{3}_{5}$ or at vertex
2 – graphs $G^{3}_{6}$ and $G^{3}_{7}$. Now let’s write all binary graphs in
reverse order. We will obtain graphs $G^{3}_{8}$ – $G^{3}_{1}2$ with two
vertices and graphs $G^{3}_{13}$ and $G^{3}_{14}$ with one vertex.

In Fig.5, each of the graphs $G^{3}_{k}$ is depicted above its number $k$).

All possible saddle-source bifurcations of these graphs are depicted with blue
arrows.

In Fig.5, all possible bifurcations of the saddle connection are shown with
three saddles.

So, fair

###### Theorem 3

The sphere has the following structures of codimensionality 1 flows:

  * •

56 structures with seven singular points;

  * •

20 structures with eight singular points.

##  5 Flows with four saddles

Graphs with four edges (see Fig.6) are of the following types:

a) with five vertices (1–3),

b) with four vertices (4–12),

c) with three vertices (13–26),

d) with two vertices (dual to 4–12),

e) with one vertex (dual to 1–3).

Figure 6: Vector fields with four saddles: the saddle-node case

Each graph of type e) is dual to a graph of type a). Graphs of type d) are
dual to graphs of type b). Therefore, only graphs with five, four or three
vertices and all their possible saddle-node bifurcations using blue triangles
are shown in Fig.6.

Saddle-node bifurcations give inverse saddle-node bifurcations on a dual
graph. Therefore, in order to find the total number of flows of codimensional
1 of this type, we multiply the number of fields (64) given by graphs 1–12 by
2 and add the number of fields (89) given by graphs 13–26. Finally we have:

###### Theorem 4

On the sphere there are 217 different codimensional 1 vector field structures
with nine singular points.

Next, we consider bifurcations of the saddle connection in flows with four
saddles (10 singular points). Two cases are possible for an edge perpendicular
to a T-vertex: 1) it does not break the graph, 2) it breaks the graph.

In the first case, if it is removed, we obtain a connected graph with three
edges (if two edges that are incident T-vertices are connected into one edge).
The removed edge is defined by a perpendicular vector in the T-vertex and one
of the corners of the face in which it is directed. We add one to these
numbers, which corresponds to the case when the perpendicular edge has a
vertex of degree 1 at one of its ends. All possible such cases (130) are shown
in Fig.7.

Figure 7: Vector fields with four saddles: a case of the saddle connection

Consider the case when a perpendicular edge divides the graph into two
subgraphs. Let’s glue two edges at the vertex T into one and denote this
subgraph by $K$, and the other subgraph by $L$. Then $K$ can contain one, two
or three edges. In the latter case, the graph $L$ has no edges and only one
vertex, which we considered earlier. If the graph $K$ has one edge
($G^{1}_{1}$ or $G^{1}_{2}$), then in each of them the perpendicular vector is
determined uniquely (up to homeomorphism), that is, we have 2 options. On
graphs with two edges, there are 8 different possibilities to choose the
corner at the vertices ($G^{2}_{1}$ – 1, $G^{2}_{2}$ – 2, $G^{2}_{3}$ – 3,
$G^{2}_{4}$ – 2) . Therefore, $2\times 8=16$ fields of this type are possible.

If $K$ has two edges, then there are 7 ways for the perpendicular vector to
the edges: $G^{2}_{1}$ – 1, $G^{2}_{2}$ – 1, $G^{2}_{3}$ – 3, $G^{2}_{4}$ – 2.
The graph $L$ contains one edge and in each of two such graphs the angle is
specified uniquely. Therefore, we have $7\times 2=14$ fields of this type.

So, fair

###### Theorem 5

There are 160 different flow structures of codimensionality 1 with ten
singular points on the sphere.

## Conclusion

All possible structures of flows of codimensionality 1 (typical bifurcations)
on a two-dimensional sphere with no more than four saddles (10 singular
points) were found: with three singular points – 2 structures, with four –
none, with five points – 10, with six – 4, with seven – 56, with eight – 20,
with nine – 217 and with ten points – 160 structures. We hope that the
research carried out in this paper can be extended to other surfaces: spheres
with holes, surfaces of a larger genus, in particular non-orientable ones.

## References

  * [1] O. Akchurin, S. Bilun, and A. Prishlyak.  Three-color graph as the 1-skeleton of the 2-sphere triangulation.  arXiv preprint arXiv:2209.05737, 2022. doi:10.48550/ARXIV.2209.05737.
  * [2] S. Bilun and A. Prishlyak.  The closed Morse 1-forms on closed surfaces.  Visn., Mat. Mekh., Kyv. Univ. Im. Tarasa Shevchenka, 2002(8):77–81, 2002.
  * [3] S. Bilun and A. Prishlyak.  Visualization of Morse flow with two saddles on 3-sphere diagrams.  arXiv preprint arXiv:2209.12174, 2022. doi:10.48550/ARXIV.2209.12174.
  * [4] S. Bilun, A. Prishlyak, and A. Prus.  Morse flows with fixed points on the boundary of 3-manifold.  arXiv preprint arXiv:2209.04019, 2022. doi:10.48550/arXiv.2209.04019.
  * [5] S. Bilun, M. Hrechko, O. Myshnova, A. Prishlyak.  Structures of optimal discrete gradient vector fields on surface with one or two critical cells  arXiv preprint arXiv:2303.07258, 2023. doi:10.48550/arXiv.2303.07258.
  * [6] S. Bilun, A. Prishlyak, S. Stas, A. Vlasenko.  Topological structure of Morse functions on the projective plane  arXiv preprint arXiv:2303.03850, 2023. doi:10.48550/arXiv.2303.03850.
  * [7] A.V. Bolsinov and A.T. Fomenko.  Integrable Hamiltonian systems. Geometry, Topology, Classification.  A CRC Press Company, Boca Raton London New York Washington, D.C., 2004\.  724 p.
  * [8] R. Forman. Morse theory for cell complexes. Advances in Mathematics. 134 (1): 90–145. 1998. doi:10.1006/aima.1997.1650.
  * [9] J. Gross and T. Tucker.  Topological graph theory.  Wiley, 1987.
  * [10] F. Harary.  Graph Theory.  Westview Press, 1969.
  * [11] Ch. Hatamian and A. Prishlyak.  Heegaard diagrams and optimal morse flows on non-orientable 3-manifolds of genus 1 and genus 2.  Proceedings of the International Geometry Center, 13(3):33–48, 2020\. doi:10.15673/tmgc.v13i3.1779.
  * [12] P. Hilton and S. Wylie.  Homology theory: An introduction to algebraic topology.  Cambridge University Press, 1968.
  * [13] B. I. Hladysh and A. O. Pryshlyak.  Functions with nondegenerate critical points on the boundary of the surface.  Ukrainian Mathematical Journal, 68(1):29–41, 2016. doi:10.1007/s11253-016-1206-5.
  * [14] B. I. Hladysh and A. O. Pryshlyak.  Deformations in the General Position of the Optimal Functions on Oriented Surfaces with Boundary.  Ukrainian Mathematical Journal, 71(8):1173–1185, 2020. doi:10.1007/s11253-019-01706-8.
  * [15] B.I. Hladysh and A.O. Prishlyak.  Topology of functions with isolated critical points on the boundary of a 2-dimensional manifold.  SIGMA. Symmetry, Integrability and Geometry: Methods and Applications, 13:050, 2017. doi:0.3842/SIGMA.2017.050.
  * [16] B.I. Hladysh and A.O. Prishlyak.  Simple morse functions on an oriented surface with boundary.  Журнал математической физики, анализа, геометрии, 15(3):354–368, 2019. doi:10.15407/mag15.03.354.
  * [17] A.S. Kronrod.  Functions of two variables.  Russian Mathematical Surveys, 5:24–134, 1950.
  * [18] Z. Kybalko, A. Prishlyak, and R. Shchurko.  Trajectory equivalence of optimal Morse flows on closed surfaces.  Proc. Int. Geom. Cent., 11(1):12–26, 2018. doi:10.15673/tmgc.v11i1.916.
  * [19] M. Losieva and A. Prishlyak.  Topology of morse–smale flows with singularities on the boundary of a two-dimensional disk.  Pr. Mizhnar. Heometr. Tsentr, 9(2):32–41, 2016. doi:10.15673/tmgc.v9i2.279.
  * [20] D.P. Lychak and A.O. Prishlyak.  Morse functions and flows on nonorientable surfaces.  Methods of Functional Analysis and Topology, 15(03):251–258, 2009\.
  * [21] A.A. Oshemkov and V.V. Sharko.  Classication of morse-smale flows on two-dimensional manifolds.  Matem. Sbornik, 189(8):93–140, 1998.
  * [22] M.M. Peixoto.  On the classication of flows of 2-manifolds.  Dynamical Systems (Proc. Symp. Univ. of Bahia, Salvador, Brasil, 1971), pages 389–419, 1973.
  * [23] A.O. Prishlyak.  On graphs embedded in a surface.  Russian Mathematical Surveys, 52(4):844, 1997. doi:10.1070/RM1997v052n04ABEH002074.
  * [24] A.O. Prishlyak.  Morse–smale vector fields without closed trajectories on-manifolds.  Mathematical Notes, 71(1-2):230–235, 2002. doi:10.1023/A:1013963315626.
  * [25] A.O. Prishlyak.  On sum of indices of flow with isolated fixed points on a stratified sets.  Zhurnal Matematicheskoi Fiziki, Analiza, Geometrii [Journal of Mathematical Physics, Analysis, Geometry], 10(1):106–115, 2003.
  * [26] A.O. Prishlyak.  Complete topological invariants of morse–smale flows and handle decompositions of 3-manifolds.  Fundamentalnaya i Prikladnaya Matematika, 11(4):185–196, 2005.
  * [27] A.O. Prishlyak.  Complete topological invariants of morse-smale flows and handle decompositions of 3-manifolds.  Journal of Mathematical Sciences, 144:4492–4499, 2007.
  * [28] A. Prishlyak and L. Di Beo.  Flows with minimal number of singularities on the Boy’s surface  Proceedings of the International Geometry Center, 15(1):32–49, 2020\.
  * [29] A. Prishlyak and M. Loseva.  Topological structure of optimal flows on the girl’s surface.  Proceedings of the International Geometry Center, 15(3-4):184–202, 2022.
  * [30] A. Prishlyak, A. Prus, and S. Huraka.  Flows with collective dynamics on a sphere.  Proc. Int. Geom. Cent, 14(1):61–80, 2021. doi:10.15673/tmgc.v14i1.1902.
  * [31] A. Prishlyak and M. Loseva.  Topology of optimal flows with collective dynamics on closed orientable surfaces.  Proceedings of the International Geometry Center, 13(2):50–67, 2020\. doi:10.15673/tmgc.v13i2.1731.
  * [32] A. Prishlyak and A. Prus.  Morse-smale flows on torus with hole.  Proc. Int. Geom. Cent., 10(1):47–58, 2017. doi:10.15673/tmgc.v1i10.549.
  * [33] A.O. Prishlyak.  Equivalence of morse function on 3-manifolds.  Methods of Func. Ann. and Topology, 5(3):49–53, 1999.
  * [34] A.O. Prishlyak.  Conjugacy of morse functions on surfaces with values on a straight line and circle.  Ukrainian Mathematical Journal, 52(10):1623–1627, 2000. doi:10.1023/A:1010461319703.
  * [35] A.O. Prishlyak.  Morse functions with finite number of singularities on a plane.  Meth. Funct. Anal. Topol, 8:75–78, 2002.
  * [36] A.O. Prishlyak.  Topological equivalence of morse–smale vector fields with beh2 on three-dimensional manifolds.  Ukrainian Mathematical Journal, 54(4):603–612, 2002.
  * [37] A.O. Prishlyak.  Topological classification of m-fields on two-and three-dimensional manifolds with boundary.  Ukrainian Mathematical Journal, 55(6):966–973, 2003.
  * [38] A.O. Prishlyak and M.V. Loseva.  Optimal morse–smale flows with singularities on the boundary of a surface.  Journal of Mathematical Sciences, 243:279–286, 2019.
  * [39] A.O. Prishlyak and K.I. Mischenko.  Classification of noncompact surfaces with boundary.  Methods of Functional Analysis and Topology, 13(01):62–66, 2007\.
  * [40] A.O. Prishlyak and A.A. Prus.  Three-color graph of the morse flow on a compact surface with boundary.  Journal of Mathematical Sciences, 249(4):661–672, 2020. doi:10.1007/s10958-020-04964-1.
  * [41] G. Reeb.  Sur les points singuliers d’une forme de pfaff complétement intégrable ou d’une fonction numérique.  C.R.A.S. Paris, 222:847—849, 1946.
  * [42] V.V. Sharko.  Functions on manifolds. Algebraic and topological aspects., volume 131 of Translations of Mathematical Monographs.  American Mathematical Society, Providence, RI, 1993.
  * [43] S. Smale.  On gradient dynamical systems.  Ann. of Math., 74:199–206, 1961.
  * [44] С.В. Білун, А.О. Гагай, О.О. Ворончук, M.В. Лосєва.  Деформації замкнених 1-форм Морса на замкнених поверхнях.  Праці Міжнародного геометричного центру, 6(4):52–56, 2013\.
  * [45] В.М. Кузаконь, В.Ф. Кириченко, and О.О. Пришляк.  Гладкi многовиди. Геометричнi та топологiчнi аспекти.  Працi Iнcтитуту математики НАН України.—2013.—97.—500 с, 2013.
  * [46] А.О. Пришляк.  Векторные поля Морса–Смейла с конечным числом особых траекторий на трехмерных многообразиях.  Доповiдi НАН України, (6):43–47, 1998.
  * [47] А.О. Пришляк.  Сопряженность функций Морса.  Некоторые вопросы совр. математики. Институт математики АН Украины, Киев, 1998.
  * [48] А.О. Пришляк.  Топологическая эквивалентность функций и векторных полей Морса—Смейла на трёхмерных многообразиях.  Топология и геометрия. Труды Украинского мат. конгресса, pages 29–38, 2001.
  * [49] А.О. Пришляк.  Векторные поля Морса–Смейла без замкнутых траекторий на трехмерных многообразиях.  Математические заметки, 71(2):254–260, 2002\.
  * [50] О.О. Пришляк.  Топологія многовидів.  Київський університет, 2015.

Taras Shevchenko National University of Kyiv

Svitlana Bilun Email address<EMAIL_ADDRESS>Orcid ID: 0000-0003-2925-5392

Bohdana Hladysh Email<EMAIL_ADDRESS>Orcid ID: 0000-0001-8935-1453

Alexandr Prishlyak Email address<EMAIL_ADDRESS>Orcid ID:
0000-0002-7164-807X

Vladislav Sinitsyn Email<EMAIL_ADDRESS>Orcid ID:
0009-0000-0928-0676