Supermassive Black Holes as Possible Sources of Ultrahigh-energy Cosmic Rays

The recent unprecedented discovery of extragalactic high-energy neutrinos has enabled pinpointing their source to blazars (IceCube Collaboration et al. 2018a, 2018b), which are supermassive black holes (SMBHs) at a distance of ∼1.75 Gpc with relativistic jets directed almost exactly toward us. It is generally believed that such neutrinos are tracers of ultrahigh-energy cosmic rays (UHECRs). UHECRs are the most energetic among particles detected on Earth, with an energy E > 10¹⁸ eV that is unreachable by the current most powerful particle accelerators, such as the Large Hadron Collider with a maximum energy <10¹³ eV per beam. Constituents of UHECRs were thought to be dominated by protons, as indicated by cosmic-ray fluorescence measurements (Abbasi et al. 2010, 2015), although recent observations are suggesting heavier constituents (Aab et al. 2017). For Galactic cosmic rays, one should observe anisotropy in the arrival direction dominantly on the Galactic plane. As observed by both the Pierre Auger Observatory (Pierre Auger Collaboration et al. 2017; Aab et al. 2018) in the southern hemisphere and the Telescope Array in the northern hemisphere (Abbasi et al. 2017), UHECRs with energy >10¹⁸ eV are extragalactic with very high confidence level. The spectrum of cosmic rays exhibits the presence of the so-called knee and ankle. Cosmic rays with energy up to ∼10^15.5 eV (knee) are generally believed to be produced in Galactic supernova explosions, while the significant lowering of flux between the knee and 10^18.5 eV (ankle) suggests a change in the source of such particles.

The flux of cosmic rays with energies >5 × 10¹⁹ eV is extremely low, which causes the main difficulty in unveiling their source and physics. In order to explain the highest energy cosmic rays, several exotic scenarios have been proposed, including extra dimensions, violation of Lorentz invariance (Bhattacharjee & Sigl 2000; Rubtsov et al. 2017), and the existence of new, exotic particles (Domokos & Kovesi-Domokos 1999). By analyzing the radio images of the blazar jets, Britzen et al. (2019) suggested that high-energy neutrino can possibly be explained by the collision of two jets. Among the astrophysical acceleration mechanisms for UHECRs, relativistic shocks in a plasma of relativistic jets have been previously considered among the most plausible (Blandford 2000). However, the recent results and estimates may indicate (Bell et al. 2018) that shock acceleration is not able to account for UHECR energies above 10²⁰ eV. Therefore, the production and acceleration mechanisms of UHECRs remain unclear.

Remarkably, an SMBH is the largest energy reservoir in the universe. Based on the irreducible mass for a black hole (Christodoulou & Ruffini 1971), it turns out that a rotating black hole has a maximum of 29% of its mass in rotational energy that is available for extraction and can be transformed into energy of accelerated particles. For an SMBH with typical mass of M = 10⁹ M_⊙ and dimensionless spin parameter a = 0.5, the available energy for extraction is of the order of E_BH ≈ 10⁷⁴ eV. It is therefore most pertinent to tap this enormous source as effectively and ultraefficiently as possible.

In this paper we invoke a novel and ultraefficient regime of the magnetic Penrose process (MPP) that electromagnetically extracts a black hole's rotational energy for accelerating cosmic-ray particles to ultrahigh energy beyond 10²⁰ eV.

1.1. Magnetosphere of an Astrophysical Black Hole

In realistic astrophysical scenarios, an arbitrary electromagnetic field that is present around a black hole is weak in the sense that its stress-energy tensor does not contribute to the spacetime geometry. This condition for a magnetic field of strength B and a black hole of mass M reads as follows (Gal'tsov & Petukhov 1978):

$$\begin{eqnarray}&&B\ll {B}_{{\rm{G}}}=\displaystyle \frac{{c}^{4}}{{G}^{3/2}{M}_{\odot }}\left(\displaystyle \frac{{M}_{\odot }}{M}\right)\approx {10}^{19}\displaystyle \frac{{M}_{\odot }}{M}\,\mathrm{Gauss}.\end{eqnarray} \tag{ 1 }$$

Here $G$ is the gravitational constant, $c$ is the speed of light and ${M}_{\odot }$ is mass of the Sun. The same value of 10¹⁹ holds for an electric field strength measured in statV/cm. Obviously, these conditions are satisfied in all known astrophysical scenarios. Therefore, the spacetime curvature around a black hole can be fully described by the standard Kerr metric. Hence, an arbitrary electromagnetic field surrounding a black hole can be considered as a test field in an axially symmetric Kerr spacetime background. This weakness of the electromagnetic field is compensated for by a large charge-to-mass ratio e/m for electrons and protons, whose motion will be essentially affected by magnetic fields already of the order of a few Gauss (Tursunov et al. 2016; Kološ et al. 2017).

The first vacuum solution of Maxwell equations in curved spacetime, namely the uniform magnetic field in the Kerr metric, was found by Wald (1974). Later, the effect of a plasma in the force-free approximation was considered in the well-known work of Blandford & Znajek (1977). Due to a lack of direct measurements of the exact shapes of field configurations around realistic black hole candidates, several numerical techniques have been employed that have shown strong connections between the shape of the magnetosphere and the characteristics of the accretion process (see, e.g., Punsly 2001; Meier 2012). In the presence of a plasma, the magnetosphere has a more complicated structure, as has been shown by full-fledged general relativistic magnetohydrodynamic simulations (see, e.g., Tchekhovskoy 2015; Janiuk et al. 2018; Nakamura et al. 2018; Porth et al. 2019). Despite that complexity, several analytical solutions to the black hole magnetosphere within the force-free approximation have been proposed (see, e.g., Gralla & Jacobson 2014; Grignani et al. 2019), where plasma is assumed to be in equilibrium with the electromagnetic field being magnetically dominated.

Although electromagnetic fields generated by a plasma, in general, have nonvanishing electric field components, it is usually expected that the electric field is effectively screened by the plasma, which makes the system magnetically dominated, that is, B² − E² > 0. This expectation is satisfied in the vacuum Wald (1974) solution everywhere outside the horizon, as well as in the force-free approach (Blandford & Znajek 1977; Gralla & Jacobson 2014). However, in the plasma-filled version of the Wald solution, B² − E² can turn negative within the ergosphere of the black hole (Komissarov 2004). Special interest is also paid in the literature to the boundary case B² = E² that occurs in a force-free plasma around a Kerr black hole (Menon & Dermer 2007; Brennan et al. 2013).

Further, we shall rely on the case in which the magnetic field contribution is dominant everywhere outside the horizon, and we use the Wald solution for our estimates, in which this condition is satisfied. Near the horizon, the strength of the induced electric field can be comparable to that of the magnetic field, while decreasing as the inverse square of the distance from the black hole. The electric field in this case is induced by a frame-dragging effect of twisting of magnetic field lines, which can be seen in the following example. Let us assume that the magnetic field is determined by at least one nonzero component of the four-vector potential so that ${A}^{\phi }\ne 0$, which corresponds to the axially symmetric configuration. We should also assume that A^t = 0 since any sufficient excess of charge in a plasma will be effectively screened. One can see that two covariant components of the electromagnetic potential are nonvanishing, namely ${A}_{t}={g}_{t\phi }{A}^{\phi }$ and ${A}_{\phi }={g}_{\phi \phi }{A}^{\phi }$. We give in detail in Section 4 the explicit form of these components in the Wald solution. Here we note that the Wald solution refers to a test field in otherwise empty space around a rotating black hole. The introduction of the plasma in this scenario does not alter spacetime geometry maintaining the test field condition. Plasma like the electromagnetic field would also share the same symmetries of axial symmetry and stationarity. Hence the Wald solution could be put forward without much hesitation even in the presence of plasma. Here we conclude that the Wald solution can be considered as a simple approximation to the realistic black hole magnetosphere and is well suited for both estimation purposes and understanding the leading-order contributions of the corresponding equations. This configuration has been effectively applied in the past for the explanation of various high-energy astrophysical phenomena (see, e.g., Kološ et al. 2017; Levin et al. 2018; Rueda et al. 2019; Ruffini et al. 2019; Tursunov et al. 2019; Stuchlík et al. 2020).

An extension of the Wald solution to moving black holes in binaries has been obtained by Morozova et al. (2014). The force-free approach to the similar problem of binary black holes has been studied by Alic et al. (2012) and Moesta et al. (2012).

The MPP was established in the mid-1980s by Wagh et al. (1985), Parthasarathy et al. (1986), and Bhat et al. (1985) as the process allowing the extraction of energy with efficiency exceeding η = 1 (for a review of early results, see Wagh & Dadhich 1989). In this section, we will show that the efficiency of the MPP, under certain conditions, can exceed η ∼ 10¹². We define efficiency of energy extraction in a standard manner as the ratio between gain and input energies. In particular, it appears that, depending on the initial setup, MPP can work in three basic regimes of efficiency: low, moderate, and ultrahigh. The latter case is able to provide ultrahigh energy for charged particles escaping from the vicinity of black holes in a straightforward manner, for characteristic values of magnetic field and even relatively moderate black hole spin. Below, we first discuss the original Penrose process (PP) and its relation to other competing mechanisms, giving brief historical remarks, and we derive its novel, ultraefficient regime at the end of this section.

The lower limit of MPP refers to the process originally discovered by Penrose (1969) in the absence of an external magnetic field. It is envisaged that a freely falling particle splits into two fragments inside the ergosphere, one of which can attain a negative energy relative to the observer at infinity, while the other, respecting energy conservation, comes out with an energy greater than that of the incident particle. Accretion of the negative-energy particle onto a black hole amounts to a negative energy flux, which is equivalent to extraction of energy from the black hole, and the only energy available for extraction is rotational. The maximum efficiency in this regime is only η_PP = 0.21 for an extremely rotating black hole. Moreover, as shown by Bardeen et al. (1972), PP requires the relative velocity between two fragments to be greater than 1/2c, and there is no conceivable astrophysical mechanism that can instantaneously accelerate particles to such high velocity. Thus, PP in the absence of electromagnetic interactions was a novel and purely geometric process, but it was not astrophysically viable as a power engine for high-energy sources.

PP was transformed into MPP by taking into account the presence of the magnetic field produced by the surrounding plasma dynamics. Now the energy required for a particle to ride on a negative-energy orbit could come from an electromagnetic interaction removing all constraints on relative velocity, and thus the process gets revived astrophysically (Wagh et al. 1985). Further, it was also shown by Parthasarathy et al. (1986) that its efficiency could exceed η > 1 for discrete particle accretion, a prediction that has been verified by fully relativistic MHD flow simulations in Narayan et al. (2014).

In a plasma setup, another process that could extract energy electromagnetically is the well-known Blandford & Znajek (1977) mechanism (BZ). It works on the principle that twisting of magnetic field lines by frame dragging produces a quadrupole electric potential difference between pole and equatorial plane, the discharge of which drives energy and angular momentum out from the hole. It is generally believed to be the leading mechanism for powering relativistic jets observed in a variety of black hole candidates. It was shown, however, that MPP is a more general process than BZ by Dadhich et al. (2018), because the latter requires the presence of a threshold magnetic field that is of the order of ∼10⁴ G, while MPP works for the entire range of the magnetic field. The latter could be thought of as the high magnetic field limit of the former. General relativistic MHD simulations by Lasota et al. (2014), Narayan et al. (2014), and Nemmen & Tchekhovskoy (2015) have shown that the energy extraction efficiency of this process is moderately high (η ≤ 10), but not ultrahigh, exceeding only a few hundred percent for a polarized plasma in a magnetic field. This is the moderate regime of MPP.

A third and the most efficient regime of MPP can accelerate charged particles to velocities much higher than one can hope to achieve in the above-described moderate regime, including the Blandford–Znajek mechanism (which by its setup uses charged matter only). Here a neutral particle is supposed to split into charged fragments in the ergosphere of a rotating black hole in the presence of an external magnetic field. As in the case of BZ, the twisting of magnetic field lines by frame dragging produces an electric field that can be associated with the electric charge of the black hole. In fact, in both the vacuum and plasma cases, the black hole acquires a net electric charge proportional to the black hole spin (see, e.g., Beskin 1997; Punsly 2001; Komissarov 2004; Levin et al. 2018; Zajaček et al. 2018; Zajacek & Tursunov 2019, and references therein).

A neutral particle can reach arbitrarily close to the horizon without being influenced by the electromagnetic field, and hence the split into charged fragments could occur very close to the horizon, and thereby an infalling charged fragment in addition to gravitational/geometric negative energy would have a very strong coulombic contribution, tremendously enhancing the quantum of energy being extracted (see a schematic representation of the model in Figure 1). This makes the process "ultra" efficient. For an idealized plasma or any other environment (containing charged matter only), as obtained for BZ, the point of splitting cannot occur very close to the horizon and hence cannot have the advantage of a tremendous gain of coulombic contribution by one of the charged fragments. This is why the efficiency of MPP in the moderate regime or BZ remains in the moderate range of an order of a few, as shown, for example, by Nemmen & Tchekhovskoy (2015), and hence could not reach the ultrahigh range. Although the proposed model is quite general, as a particular example we also consider beta decay of neutrons (which can appear, e.g., due to the nucleosynthesis process in the hot and dense plasma of an accretion disk; Janiuk 2014) in the dynamical environment of an SMBH, from which it follows that protons after decay can naturally reach an energy >10²⁰ eV for the characteristic values of magnetic field of order 10⁴ G and black hole mass 10⁹ M_⊙.

Zoom In Zoom Out Reset image size

Below we provide the main equations supporting the discussions given above. It is generally assumed in all studies of test fields around a rotating black hole that it shares symmetries of axial symmetry and stationarity. It should be noted that space around the black hole shares the black hole rotation in what is known as the "frame-dragging" phenomenon. Therefore the electromagnetic field and the plasma share these symmetry properties. This implies that the four-vector potential of the electromagnetic field A_μ has two nonvanishing covariant components A_t and A_ϕ (see discussion in Section 1.1 and Wald 1974). The presence of the magnetic field modifies the canonical four-momentum of charged test particles according to ${P}_{\mu }={{mu}}_{\mu }+{{qA}}_{\mu }$, where m, q, and u^μ are the mass, charge, and four-velocity of a test particle. The conserved components of the four-momentum, namely energy E and angular momentum L, can be written in the form

$$\begin{eqnarray}&&-E={P}_{t}={{mu}}_{t}+{{qA}}_{t},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&L={P}_{\phi }={{mu}}_{\phi }+{{qA}}_{\phi }.\end{eqnarray} \tag{ 3 }$$

The dynamics of charged particles around a Kerr black hole in the presence of a magnetic field have been widely studied in the literature (see, e.g., Stuchlík & Kološ 2016; Tursunov et al. 2016, 2018, and references therein). In addition to conservation of energy and angular momentum, the normalization of the four-velocity ${u}^{\alpha }{u}_{\alpha }=-\delta$ holds for both neutral and charged particles, where δ = 1 for massive particles and δ = 0 for massless particles. The energy of a particle is minimal at the equatorial plane, for which the four-velocity can be rewritten as ${u}^{\alpha }={u}^{t}(1,v,0,{\rm{\Omega }})$, where we denote by v = dr/dt and ${\rm{\Omega }}=d\phi /{dt}\equiv {u}^{\phi }/{u}^{t}$ the radial velocity and the angular velocity of the test particle as measured by an asymptotic observer, respectively. Substituting this into the normalization condition, we get the expression for Ω in the form

$$\begin{eqnarray}&&{\rm{\Omega }}=\displaystyle \frac{1}{D}\left(-\,{{Hg}}_{t\phi }\pm \sqrt{{u}_{t}^{2}\left({{Hg}}^{2}-{{Dg}}_{{rr}}{v}^{2}\right)}\right),\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&D=\delta \,{g}_{t\phi }^{2}+{u}_{t}^{2}\,{g}_{\phi \phi },\quad H=\delta \,{g}_{{tt}}+{u}_{t}^{2},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{g}^{2}={g}_{t\phi }^{2}-{g}_{\phi \phi }{g}_{{tt}},\quad {u}_{t}=-\left(E+{{qA}}_{t}\right)/m.\end{eqnarray} \tag{ 6 }$$

The sign in Equation (4) depends on whether the particle is corotating or counterrotating with respect to a locally nonrotating reference frame. Possible values of Ω are restricted by the limit of u^α tending to a null vector (Parthasarathy et al. 1986), that is, ${{\rm{\Omega }}}_{-}\leqslant {\rm{\Omega }}\leqslant {{\rm{\Omega }}}_{+}$, where Ω_± takes the form

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\pm }=\displaystyle \frac{1}{{g}_{\phi \phi }}\left(-{g}_{t\phi }\pm \sqrt{{g}_{t\phi }^{2}-{g}_{{tt}}{g}_{\phi \phi }}\right).\end{eqnarray} \tag{ 7 }$$

Let us now consider the split of a particle 1, not necessarily neutral, into two charged fragments 2 and 3 in the black hole ergosphere at the equatorial plane. The conservation laws before and after the split can be written in the form

$$\begin{eqnarray}&&{E}_{1}={E}_{2}+{E}_{3},{L}_{1}={L}_{2}+{L}_{3},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{q}_{1}={q}_{2}+{q}_{3},{m}_{1}\geqslant {m}_{2}+{m}_{3},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{m}_{1}{\dot{r}}_{1}={m}_{2}{\dot{r}}_{2}+{m}_{3}{\dot{r}}_{3},0={m}_{2}{\dot{\theta }}_{2}+{m}_{3}{\dot{\theta }}_{3},\end{eqnarray} \tag{ 10 }$$

where dots denote the derivative with respect to the proper time. If one of the particles after the split, for example particle 2, attains a negative energy, particle 3 comes out with an energy exceeding the energy of incident particle 1 at the expense of the rotational energy of the black hole. Conservation of the four-momentum at the splitting point leads to the following relation (Bhat et al. 1985):

$$\begin{eqnarray}&&{m}_{1}{u}_{1}^{\phi }={m}_{2}{u}_{2}^{\phi }+{m}_{3}{u}_{3}^{\phi }.\end{eqnarray} \tag{ 11 }$$

Substituting ${u}^{\phi }={\rm{\Omega }}\,{u}^{t}=-{\rm{\Omega }}\,Y/X$, where $Y=(E+q\,{A}_{t})/m$ and $X={g}_{{tt}}+{\rm{\Omega }}\,{g}_{t\phi }$, one can rewrite Equation (11) in the form

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{1}{m}_{1}{Y}_{1}\displaystyle \frac{{X}_{2}{X}_{3}}{{X}_{1}}={{\rm{\Omega }}}_{2}{m}_{2}{Y}_{2}{X}_{3}+{{\rm{\Omega }}}_{3}{m}_{3}{Y}_{3}{X}_{2}.\end{eqnarray} \tag{ 12 }$$

This equation leads to the final expression for the energy of escaping particle 3, which after several algebraic steps takes the following form:

$$\begin{eqnarray}&&{E}_{3}=\chi ({E}_{1}+{q}_{1}{A}_{t})-{q}_{3}{A}_{t},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&\chi =\displaystyle \frac{{{\rm{\Omega }}}_{1}-{{\rm{\Omega }}}_{2}}{{{\rm{\Omega }}}_{3}-{{\rm{\Omega }}}_{2}}\displaystyle \frac{{X}_{3}}{{X}_{1}},\quad {X}_{j}={g}_{{tt}}+{{\rm{\Omega }}}_{j}{g}_{t\phi },\end{eqnarray} \tag{ 14 }$$

where ${{\rm{\Omega }}}_{j}={(d\phi /{dt})}_{j}$ is an angular velocity of the jth particle, given by Equation (4) with the values limited by Equation (7).

Let us now define the efficiency of energy extraction as the ratio between gained and infalling energies:

$$\begin{eqnarray}&&\eta =({E}_{3}-{E}_{1})/{E}_{1}=-{E}_{2}/{E}_{1}.\end{eqnarray} \tag{ 15 }$$

If all particles are massive and charged, we obtain the following expression for efficiency:

$$\begin{eqnarray}&&\eta ={\eta }_{\mathrm{PP}}+\displaystyle \frac{{q}_{3}{A}_{t}-{q}_{1}{A}_{t}({\eta }_{\mathrm{PP}}+1)}{{m}_{1}{u}_{t1}+{q}_{1}{A}_{t}},\end{eqnarray} \tag{ 16 }$$

where all quantities are calculated at the point of split, and η_PP is the efficiency of the original Penrose (1969) process given by purely geometric factors. For the split close to the event horizon, η_PP takes the form

$$\begin{eqnarray}&&{\eta }_{\mathrm{PP}}=\displaystyle \frac{1}{2a}\left(\sqrt{2}\sqrt{1-\sqrt{1-{a}^{2}}}-a\right),\end{eqnarray} \tag{ 17 }$$

where a is the dimensionless spin of a black hole. This expression coincides with the results of Penrose (1969) and Bardeen et al. (1972).

Depending on whether the particles before and after the split are charged or neutral, one can distinguish three different regimes of efficiency. In the absence of electromagnetic fields or if all particles are neutral, the expression (16) turns into (17), given by purely geometric factors with the maximum value of 0.21 for an extremely rotating black hole (Penrose 1969).

If all particles are charged, one can see that in the astrophysically relevant conditions for elementary particles, such as electrons and protons, the following relation holds: $| \tfrac{q}{m}{A}_{t}| \gg | {u}_{t}|$. This inequality implies that in realistic conditions, the motion of a charged particle is dominated by the electromagnetic field. When a particle is neutral, its specific energy is given by E/m = −u_t, which is of the order of unity for a freely falling particle. When the particle is charged, the expression for specific energy $E/m=-{u}_{t}-\tfrac{q}{m}{A}_{t}$ is sufficiently influenced by the factor $\tfrac{q}{m}{A}_{t}$ that dominates the dynamics because of the very large value of the charge-to-mass ratio q/m for elementary particles. A more precise estimate can be obtained numerically (see, e.g., Kološ et al. 2017). Thus, the efficiency of MPP (16) in the moderate regime can be reduced to the following simple form:

$$\begin{eqnarray}&&{\eta }_{\mathrm{mod}.}\approx \displaystyle \frac{{q}_{3}}{{q}_{1}}-1.\end{eqnarray} \tag{ 18 }$$

The extraction of energy from the black hole occurs when ${q}_{3}\gt {q}_{1}$. Obviously, if particles before and after the split are charged, the efficiency remains moderate, reaching an order of a few in realistic conditions, because a plasma surrounding a black hole is usually considered to be neutral.

The situation changes dramatically if the incident particle is neutral (q₁ = 0) and splits into two charged fragments. In this case, the leading contribution to the efficiency is the third term on the right-hand side of Equation (16). The expression for efficiency in this case takes the form

$$\begin{eqnarray}&&{\eta }_{\mathrm{ultra}}={\eta }_{\mathrm{PP}}+\displaystyle \frac{{q}_{3}}{{m}_{1}}{A}_{t}\approx \displaystyle \frac{{q}_{3}}{{m}_{1}}{A}_{t}.\end{eqnarray} \tag{ 19 }$$

MPP in this regime is ultraefficient. The energy of the escaping particle according to Equation (15) is then given by

$$\begin{eqnarray}&&{E}_{3}=\left({\eta }_{\mathrm{ultra}}+1\right){E}_{1},\end{eqnarray} \tag{ 20 }$$

which can grow ultrahigh as we show below.

In general, the magnetic field has a complicated structure in the vicinity of the horizon, but in a small fraction of the space where the split occurs, one can consider the field to be approximately uniform. In this case, known as the Wald solution (Wald 1974), the expression for ultraefficiency (Equation (19)) takes the following form:

$$\begin{eqnarray}&&{\eta }_{\mathrm{ultra}}=\displaystyle \frac{1}{2}\left(\sqrt{\displaystyle \frac{{r}_{g}}{{r}_{\mathrm{ion}}}}-1\right)+\displaystyle \frac{{q}_{3}B\,a\,{r}_{g}}{2{m}_{1}{c}^{2}}\left(1-\displaystyle \frac{{r}_{g}}{2{r}_{\mathrm{ion}}}\right),\end{eqnarray} \tag{ 21 }$$

where r_ion is the ionization point (splitting point) of the neutral particle, and ${r}_{g}=2{GM}/{c}^{2}$ is the gravitational radius of a black hole. For quantitative estimates, let us consider a neutron beta decay,

$$\begin{eqnarray}&&{n}^{0}\to {p}^{+}+{W}^{-}\to {p}^{+}+{e}^{-}+{\bar{\nu }}_{{\rm{e}}},\end{eqnarray} \tag{ 22 }$$

in the vicinity of an SMBH having mass M, spin a, and magnetic field of strength B. Due to the large value of the charge-to-mass ratio for a proton, one can see from (21) that the leading contribution to the efficiency is given by the second term on the right-hand side, so the efficiency can be approximately written as ${\eta }_{\mathrm{ultra}}\approx {eB}\,a\,{r}_{g}/(2{m}_{{{\rm{n}}}^{0}}{c}^{2})$, where e is the charge of the proton and ${m}_{{{\rm{n}}}^{0}}$ is the neutron mass. From Equation (20) it follows that the energy of the escaping proton after neutron beta decay is determined by the relation

$$\begin{eqnarray}&&{E}_{{{\rm{p}}}^{+}}=\left({\eta }_{\mathrm{ultra}}+1\right){E}_{{{\rm{n}}}^{0}},\end{eqnarray} \tag{ 23 }$$

were ${E}_{{{\rm{n}}}^{0}}\approx {10}^{9}$ eV is the energy of a free neutron (mass-energy). Substituting here the relation for η_ultra given by Equation (21), one can estimate the energy of an escaping proton after beta decay of a free neutron as

$$\begin{eqnarray}&&{E}_{{{\rm{p}}}^{+}}=1.7\times {10}^{20}\mathrm{eV}\left(\displaystyle \frac{B}{{10}^{4}\,{\rm{G}}}\right)\left(\displaystyle \frac{M}{{10}^{9}{M}_{\odot }}\right)\left(\displaystyle \frac{a}{0.8}\right),\end{eqnarray} \tag{ 24 }$$

predicting the energy of proton E_p exceeding 10²⁰ eV for $M\sim {10}^{9}{M}_{\odot }$ and $B\sim {10}^{4}\,{\rm{G}}$. Here, we take the decay point at ${r}_{\mathrm{ion}}={r}_{g}$, far enough from the event horizon r_h = 0.8r_g that the high-energy particle is allowed to escape to infinity, avoiding infinite redshift. In the absence of a magnetic field, the energy-extracting action has to take place very close to the horizon to produce significant efficiency in the process. This is realistically very difficult to sustain and justify. The presence of the electromagnetic interaction has completely freed this constraint. In this estimate, we also neglect the effect of the antineutrino on the final energy of the escaping proton. In Figure 1 we depict results of numerical modeling of the ionization of a neutral particle skirting the inner edge of a Keplerian accretion disk for schematic purposes. The trajectory of an escaping high-energy particle after ionization of a freely falling neutral particle from the accretion disk is indicated by the blue color. It is important to note that the escaping particle after neutron beta decay is more likely a proton in the astrophysically favorable cases. This is because the Wald charge (or any black hole charge produced by twisting of magnetic field lines) is more likely to be positive in realistic cases (see, e.g., discussions in Wald 1974; Zajaček et al. 2018; Zajacek & Tursunov 2019). A more detailed numerical analysis of the process of acceleration is given in Section 4.

The dependence of energy of the escaping proton on the magnetic field for different black hole masses is given in Figure 2 (left). The proton energy versus black hole spin is shown in Figure 2 (right) for an SMBH of mass ${10}^{9}{M}_{\odot }$ and various values of magnetic field. One can see that the process does not require extreme or rapid rotation of the black hole.

Zoom In Zoom Out Reset image size

Constraints on the magnetic field and the SMBH's mass to produce UHECR particles are given in Figure 3, where several representative SMBH candidates are pointed out. The fitting of source candidates given in Figure 3 requires measurements of magnetic fields on the event horizon scales of nearby SMBHs (within approximately 100 Mpc from the Milky Way, due to the GZK cutoff effect) in addition to mass estimates. Currently, there are only a few sources for which magnetic fields are measured on such scales with reliable methods and precisions (we use the estimates obtained in Doeleman et al. 2012; Eckart et al. 2012; Kino et al. 2015; Baczko et al. 2016, for corresponding sources). One can also see from Figure 3 that an arbitrary SMBH candidate with mass in the range ${10}^{8}\mbox{--}{10}^{9}{M}_{\odot }$ that has been observed with relativistic jets can serve as a source of protons with energies over 10²⁰ eV, if we assume that the jets are produced in a Blandford–Znajek mechanism requiring a ${10}^{3}\mbox{--}{10}^{4}\,{\rm{G}}$ field (see bar labeled BZ).

Zoom In Zoom Out Reset image size

The proposed model also gives a relatively precise estimate of the maximum energy of protons produced by the SMBH in the center of our Galaxy, namely SgrA*, which has a highly ordered magnetic field (see, e.g., Eckart et al. 2012; Eatough et al. 2013; Morris 2015) that reaches 10–100 G on the event-horizon scales. The mass of SgrA* is estimated to be ≈4 × 10⁶ M_⊙ (Eckart et al. 2012; Parsa et al. 2017). Thus, the maximum energy of a proton produced after ionization near SgrA* reaches

$$\begin{eqnarray}&&{E}_{{p}^{+}}^{\,{\mathrm{SgrA}}^{* }}=5\times {10}^{15}\,\mathrm{eV}\left(\displaystyle \frac{B}{{10}^{2}\,{\rm{G}}}\right)\left(\displaystyle \frac{M}{4\times {10}^{6}{M}_{\odot }}\right)\left(\displaystyle \frac{a}{0.5}\right).\end{eqnarray} \tag{ 25 }$$

This energy remarkably coincides with the knee of the cosmic-ray energy spectra, above which the flux of particles is suppressed.

4.1. Escape along the Rotation Axis

Boundaries of motion for charged particles in the axially symmetric field configurations can be open to infinity along the rotation axis of the black hole, coinciding with the direction of the field lines. Such a collimated corridor around the axis of rotation of a black hole does not exist for neutral particles. In general, a charged particle is allowed to escape to infinity (along the magnetic field lines) if the resulting energy of the particle is greater than its rest energy at infinity. However, in some cases, the Lorentz γ factor of a charged particle along the rotating axis, that is, the direction of escape, remains around unity if the energy of the particle is ultrahigh. A high-energy charged particle produced at the black hole's equatorial plane can have most of its energy concentrated at the oscillatory (Larmor) energy mode in the equatorial plane, so the velocity along the axis of rotation (and escape direction) remains zero or moderate. In this section, we show that large values of the Lorentz γ factor of an escaping ultrahigh-energy particle from the inner regions of the black hole accretion disk may occur only in the presence of the induced charge of the black hole. In the asymptotically uniform magnetic field, this induced black hole charge is known as the Wald charge, arising from twisting of magnetic field lines by the black hole's rotation.

Let us denote the coordinate velocity v^μ, proper velocity u^μ, and Lorentz γ factor of the test particle as follows:

$$\begin{eqnarray}&&{v}^{\alpha }=\displaystyle \frac{{{dx}}^{\alpha }}{{dt}},\quad {u}^{\alpha }=\displaystyle \frac{{{dx}}^{\alpha }}{d\tau }=\gamma {v}^{\alpha },\quad \gamma =\displaystyle \frac{{dt}}{d\tau }.\end{eqnarray} \tag{ 26 }$$

In the asymptotic limit, that is, flat spacetime filled by the homogeneous magnetic field, one can find that the energy of a charged particle measured at infinity is nothing else but

$$\begin{eqnarray}&&{E}_{\infty }=E+{{qA}}_{t}^{\infty },\end{eqnarray} \tag{ 27 }$$

where E is the integral of motion defined in Equation (2), and ${A}_{t}^{\infty }$ is the asymptotic value of the time component of the four-vector potential A_μ of the electromagnetic field that causes the difference between E_∞ and E. In the case of a Kerr black hole immersed in a uniform magnetic field B, the asymptotic limit of A_μ reads as follows:

$$\begin{eqnarray}&&{A}_{\mu }^{\infty }=\left(-{Ba},0,0,\displaystyle \frac{1}{2}{{Bg}}_{\phi \phi }\right).\end{eqnarray} \tag{ 28 }$$

Thus, the Lorentz γ factor can be derived in the form

$$\begin{eqnarray}&&m\gamma ={{mu}}^{t}=\displaystyle \frac{{dt}}{d\tau }=E+{{qA}}_{t}^{\infty }={E}_{\infty }.\end{eqnarray} \tag{ 29 }$$

On the other hand, one can see that the energy E_∞ can be decomposed into the kinetic energy in the escape (vertical) direction, E_z, and the oscillatory (Larmor) energy, E_L, in the form (for details see Stuchlík & Kološ 2016)

$$\begin{eqnarray}&&{E}_{\infty }^{2}={E}_{{\rm{z}}}^{2}+{E}_{L}^{2}.\end{eqnarray} \tag{ 30 }$$

Ejection velocities u_z = u^z and v_z = v^z and the corresponding γ_z factor in the escape direction (coinciding with the rotation axis) can be found in the form

$$\begin{eqnarray}&&{{mu}}_{{\rm{z}}}={E}_{{\rm{z}}},\quad {v}_{{\rm{z}}}=\displaystyle \frac{{E}_{{\rm{z}}}}{{E}_{\infty }},\quad {\gamma }_{{\rm{z}}}=\displaystyle \frac{1}{\sqrt{1-{v}_{{\rm{z}}}^{2}}}=\displaystyle \frac{{E}_{\infty }}{{E}_{L}}.\end{eqnarray} \tag{ 31 }$$

In flat spacetime, both energies E_z and E_L are conserved, while in the black hole vicinity a transmutation between two energy modes can be observed that can maximize the γ_z factor along the black hole rotation axis (Stuchlík & Kološ 2016).

Let us now find the condition in which γ_z is maximal. This condition depends on the exact shape of the four-vector potential A_μ, which can have different forms depending on the stage of the black hole accretion. The solution of Maxwell equations for a uniform magnetic field B in the background Kerr black hole spacetime reads as (Wald 1974)

$$\begin{eqnarray}&&{A}_{t}=\displaystyle \frac{B}{2}\left({g}_{t\phi }+2{{ag}}_{{tt}}\right),\quad {A}_{\phi }=\displaystyle \frac{B}{2}\left({g}_{\phi \phi }+2{{ag}}_{t\phi }\right).\end{eqnarray} \tag{ 32 }$$

Here, the rotation of the black hole induces the electric field due to the frame-dragging effect that gives rise to the potential difference between the event horizon and infinity. We can see that the contravariant time component of A^μ is nonzero, being A^t = aB. Therefore, the solution (32) for A_μ in this form causes a selective accretion of charged particles of the same sign into the black hole. Selective accretion into the black hole occurs until A^t = 0, when the remaining nonvanishing component appears to be A^ϕ = B/2. Covariant components of A_μ at the final stage of the selective accretion have the following components:

$$\begin{eqnarray}&&{A}_{t}=\displaystyle \frac{B}{2}{g}_{t\phi },\quad {A}_{\phi }=\displaystyle \frac{B}{2}{g}_{\phi \phi }.\end{eqnarray} \tag{ 33 }$$

At this stage, the black hole accretes a charge equal to Q_W = 2aM B, which is known as the induced Wald charge (Wald 1974). The timescale of the selective accretion process is very short for astrophysical black holes. Moreover, an induced charge (in different form) should also arise in any other axially symmetric magnetic field configuration different from uniformity. Therefore, one can conclude that any astrophysical black hole candidate possesses nonzero electric charge that is gravitationally weak, but its effect on the charged particles cannot be neglected. Below, we will show that the induced charge of a black hole plays a crucial role in the effect of local acceleration of high-energy charged particles produced in the ionization of neutral matter from the inner regions of the black hole accretion flow. In the absence of the induced charge, high-energy particles can also be created within the MPP, but their γ factors along the escape direction (coinciding with the rotation and magnetic field axes) remain close to unity.

4.2. Numerical Analysis

In order to support the above discussions quantitatively, we numerically solve the problem of ionization of initially neutral particles for a particular set of initial conditions. The most general form of the equation of motion of charged particles can be written as

$$\begin{eqnarray}&&\displaystyle \frac{{{du}}^{\mu }}{d\tau }+{{\rm{\Gamma }}}_{\alpha \beta }^{\mu }{u}^{\alpha }{u}^{\beta }=\displaystyle \frac{q}{m}{F}_{\nu }^{\mu }{u}^{\nu }+{{ \mathcal F }}_{\mathrm{rad}.}^{\mu },\end{eqnarray} \tag{ 34 }$$

where ${F}_{\mu \nu }={\partial }_{\mu }{A}_{\nu }-{\partial }_{\nu }{A}_{\mu }$ is the tensor of the external electromagnetic field, and ${{ \mathcal F }}_{\mathrm{rad}.}^{\mu }$ is the radiation reaction force that has a sophisticated character, in general (see, e.g., Poisson 2004, and references therein). In the astrophysically relevant cases (when considering the motion of charged test particles, such as protons and ions around a magnetized Kerr black hole), the radiation reaction force can be simplified to the following form (see details in Tursunov et al. 2018):

$$\begin{eqnarray}&&{{ \mathcal F }}_{\mathrm{rad}.}^{\alpha }=\displaystyle \frac{2{q}^{3}}{3{m}^{2}}\left(\displaystyle \frac{{{DF}}_{\beta }^{\alpha }}{{{dx}}^{\mu }}{u}^{\beta }{u}^{\mu }\right.\left.+\,\displaystyle \frac{q}{m}\left({F}_{\beta }^{\alpha }{F}_{\mu }^{\beta }+{F}_{\mu \nu }{F}_{\sigma }^{\nu }{u}^{\sigma }{u}^{\alpha }\right){u}^{\mu }\right).\end{eqnarray} \tag{ 35 }$$

Radiation losses of cosmic rays due to synchrotron radiation in magnetic fields along the propagation distance are discussed in the following subsection. In order to parameterize the equations of motion (34), we introduce the following dimensionless parameter reflecting the relative influence of the Lorentz and gravitational forces:

$$\begin{eqnarray}&&{ \mathcal B }=\displaystyle \frac{{qGMB}}{2{{mc}}^{4}}.\end{eqnarray} \tag{ 36 }$$

Then we solve Equation (34) numerically, applying a Kerr black hole metric immersed in an external uniform magnetic field in the two limiting cases:

• 1.  
  uncharged black hole with A_μ given by Equation (32),
• 2.  
  black hole with induced Wald charge, $Q={Q}_{{\rm{W}}}\equiv 2{aMB}$, that is, with A_μ given by Equation (33).

The results of numerical modeling of MPP, that is, the ionization of a neutral particle near a black hole and the resulting fate of the charged particle, are shown in Figure 4 for uncharged and charged black hole cases. The description of the figure is given in the caption. Here we note that, in both cases, ionization of neutral particles leads to the production of high-energy charged particles, but escape of the particle to infinity (with large escape velocity) can be observed only in the case of a black hole with induced Wald charge (bottom row of Figure 4). Increasing the magnetic parameter ${ \mathcal B }$ leads to increasing the energy γ_z and, importantly, more narrow collimation of escaping charged particles. In the absence of the induced charge of the black hole, charged particles perform an oscillatory motion around magnetic field lines (Larmor type of motion), and this oscillatory energy cannot be effectively transformed into the translational kinetic energy in the perpendicular direction. Since the presence of a magnetic field always bounds the motion of charged particles in the equatorial plane, the final fate of the particle in the absence of an induced black hole charge leads to the collapse of the particle into the black hole. Therefore, the induced charge of the black hole arising from twisting of magnetic field lines plays the role of a local accelerator of high-energy charged particles produced in the energy extraction mechanisms, such as MPP.

Zoom In Zoom Out Reset image size

It is important to note that the ionization point can be located above the ergosphere of a rotating black hole. One can see this in Figure 5, where we plot trajectories of high-energy particles produced at various points of the equatorial plane. One can see that the energy of an ionized particle is larger when the ionization point is closer to the black hole. However, the maximal γ_z factor (corresponding to velocity in the vertical direction) is achieved when the ionization occurs close to the inner edge of the accretion disk (ISCO), although the differences between values of specific energy ${{ \mathcal E }}_{1}$ and γ_z factor of an escaping charged particle for ionization points around ISCO and below are not critical. In astrophysical conditions, the charged particles should escape to infinity throughout the disk along the funnels with lower matter density. The highest energy particles should originate from around ISCO.

Zoom In Zoom Out Reset image size

In numerical results presented in Figures 4 and 5, we have used positive values of the magnetic parameter ${ \mathcal B }\gt 0$, which correspond to the motion of positively charged particles, q > 0, in a magnetic field with the field lines oriented in the same direction as the axis of rotation of the black hole. This is the most astrophysically relevant scenario because, in general, magnetic field lines (generated by plasma dynamics corotating with a black hole) are supposed to share the axial symmetry of the black hole, at least in its vicinity. This leads to the positive sign of the induced charge of the black hole (see, e.g., Wald 1974; Zajaček et al. 2018).

In addition to a positively charged particle, ionization of a neutral particle leads to the creation of a negatively charged particle with ${ \mathcal B }\lt 0$, which always has negative energy with respect to the observer at infinity within the ergosphere. Outside the ergosphere, the energy of a negatively charged particle can be positive, but it is always lower than the initial energy of an incident neutral particle. For large values of magnetic parameter $| { \mathcal B }| \gg 1$, the energy of a negatively charged particle is always negative, due to the energy conservation law given by (8). Therefore, a negatively charged particle produced in the ionization of an initially neutral particle can never escape to infinity, remaining bounded or, in most cases, falling into the black hole. This neutralizes the "rotationally" induced electric field of the black hole, which is equivalent to the extraction of rotational energy from the black hole. In the opposite case of an antiparallel orientation of the magnetic field and rotational axes, the induced charge of the black hole is negative. This implies that an escaping high-energy particle (produced in the ionization of a neutral particle) has to be negatively charged.

In realistic scenarios, the magnetic parameter ${ \mathcal B }$ can be larger for several orders of magnitude than used in our numerical plots. For protons around a typical SMBH of mass M = 10⁹ M_⊙ and characteristic value of magnetic field B ∼ 10⁴ G, the parameter ${ \mathcal B }$ has the following value:

$$\begin{eqnarray}&&{{ \mathcal B }}_{\mathrm{SMBH}}\approx 2.3\times {10}^{11}\,\left(\displaystyle \frac{B}{{10}^{4}\,{\rm{G}}}\right)\left(\displaystyle \frac{M}{{10}^{9}{M}_{\odot }}\right).\end{eqnarray} \tag{ 37 }$$

For the best-known SMBH candidate located at the center of the Milky Way, we get the following estimate:

$$\begin{eqnarray}&&{{ \mathcal B }}_{{\mathrm{SgrA}}^{* }}\approx 9.4\times {10}^{6}\,\left(\displaystyle \frac{B}{100\,{\rm{G}}}\right)\left(\displaystyle \frac{M}{4\times {10}^{6}{M}_{\odot }}\right),\end{eqnarray} \tag{ 38 }$$

This implies that the energies, γ_z factors, and collimation of escaping charged particles in realistic conditions have to be also larger than demonstrated numerically in Figures 4 and 5. This gives rise to interpreting SMBHs as possible sources of UHECRs, as discussed in Section 2, with possible candidates given in Figure 3 in particular.

In addition to the induced electric field of the black hole, generated by the frame-dragging effect, an electric field can also appear that is due to charge separation in a plasma surrounding a black hole. The circular motion of the plasma of an accretion disk around the black hole in the presence of a magnetic field (with nonzero component of the magnetic field orthogonal to the orbital plane) necessarily leads to the separation of charges in the plasma and a resulting nonvanishing component of electric field of the plasma. In this case, high-energy charged particles escaping from the inner region of the accretion disk may have an additional component of the accelerating force. Such a configuration has already been applied in an investigation of the motion of Galactic center flare components by Tursunov et al. (2019). Possible acceleration of cosmic rays by the electric field of a plasma surrounding a black hole requires further investigation.

4.3. GZK Cutoff and Synchrotron Losses

Depending on a particle's type and energy, a primary UHECR can lose a large part of its energy in interactions with photons of the cosmic microwave background while propagating over distances comparable to the size of local cosmological structures. These interactions mainly appear as photopion production and force protons with energies above 5 × 10¹⁹ eV to lose a major part of their energy. Consequently, the spectrum of protons shows suppression of flux at these energies, which is known as the GZK cutoff (Greisen 1966; Zatsepin & Kuzmin 1966). Detection of UHECRs with energies beyond the GZK cutoff implies the location of sources within a distance of ∼100 Mpc if a primary particle is a proton.

On the other hand, inevitable interaction of a UHECR with a magnetic field along the trajectory can lead to synchrotron radiation loss. Although the suppression of energy of UHECRs in Galactic and intergalactic magnetic fields is relatively small, UHECRs can lose a sufficient amount of their energies in the source regions where magnetic fields can be considerably large. For an ultrarelativistic particle with charge q and mass m, the timescale of synchrotron loss is given by (Tursunov et al. 2018)

$$\begin{eqnarray}&&\tau \approx \displaystyle \frac{3\,{m}^{3}{c}^{5}}{{q}^{4}{B}^{2}f(r)},\qquad f(r)=1-\displaystyle \frac{2{GM}}{{{rc}}^{2}}.\end{eqnarray} \tag{ 39 }$$

Suppression of proton energy on the propagation distance in a magnetic field 10⁻⁵ G is demonstrated in Figure 6 for different values of initial energy. The cubic dependence of the decay timescale (39) on the particle mass implies that electrons lose their energy ∼10¹⁰ times faster than protons do. The characteristic timescale of synchrotron energy loss for high-energy electrons propagating in a magnetic field of 10⁴ G is of the order of ∼1 s; a similar timescale for protons is ∼10¹⁰ s. Therefore, for a typical SMBH with magnetic field of order 10⁴ G, the primary UHECRs are more plausibly protons or ions, while the decay timescales of electrons are too short for escape from the SMBH vicinity.

Zoom In Zoom Out Reset image size

We propose a mechanism that suggests SMBHs are sources of UHECRs. Employing the novel, ultraefficient regime of MPP and ionization of neutral particles, such as neutron beta decay near the horizon of a spinning black hole, we have shown that a proton's energy naturally exceeds 10²⁰ eV for an SMBH of 10⁹ M_⊙ and magnetic field 10⁴ G. We list the main advantages of the model as follows:

• 1.  
  It clearly predicts SMBHs as the source of the highest energy cosmic rays.
• 2.  
  It provides verifiable constraints on the mass and magnetic field of the SMBH candidate to produce UHECRs.
• 3.  
  It operates in viable astrophysical conditions for an SMBH with moderate spin and typical magnetic field strength in its vicinity.
• 4.  
  It does not require an extended acceleration zone for a particle to reach ultrahigh energy, nor the fine-tuning of accreting-matter parameters.
• 5.  
  The energy-extracting action can take place relatively far from the event horizon without risking the infinite redshift and ultraefficiency of energy extraction.
• 6.  
  The maximum energy of a proton in the process occurring at the Galactic center SMBH (10^15.5 eV) coincides with the knee of the cosmic-ray energy spectra.

The driving engine of the process is the presence of a gravitationally induced black hole charge that arises from the magnetic field twist that is due to black hole rotation in both vacuum and plasma surroundings. Numerically comparing the trajectories of charged particles in the absence and presence of the induced charge, we have shown that production of ultrahigh-energy particles after the ionization of neutral particles can be achieved in both cases. However, large velocities in the escaping "vertical" direction can be obtained only in the presence of the induced black hole charge.

We have shown that the ionization point should not necessarily be within the ergosphere of a rotating black hole, although the energy of an ionized particle decreases with increasing distance of the ionization point from the black hole. In fact, the maximum escape velocity of a charged particle is obtained near the innermost stable circular orbit that is still outside the ergosphere in many cases.

Our numerical results were obtained in the case of a vacuum magnetic field, given by the Wald solution. We have also shown that the process should work similarly in any axially symmetric magnetic field configuration (that shares the symmetries of a background Kerr spacetime metric at least near the black hole). Production of UHECRs in the plasma MHD case should be tested, and we give a clear prediction of similar results.

The described mechanism can, in principle, be applied to neutron stars, in which the lower masses are compensated for by large values of magnetic fields. Such studies, however, we shall leave for future investigations.

Since the synchrotron radiation loss of relativistic electrons is ∼10¹⁰ times faster than for protons, heavier constituents of UHECRs seem more plausible in this scenario.

The fit of candidate SMBHs with the proposed model also requires the measurement of magnetic fields on the event-horizon scale. Nowadays, the number of such precise measurements is very few and should increase with future global VLBI observations. Constraints on the source candidates with known SMBH masses and magnetic fields are given in Figure 3.

We believe that the proposed model of an SMBH as the power engine of UHECRs opens up new vistas for understanding this remarkable high-energy phenomenon, as well as for its applications in other similar high-energy settings.

A.T. thanks the International Mobility Project CZ.02.2.69/0.0/0.0/16_027/0008521 and is grateful to the 1st Physics Institute of the University of Cologne and the Inter-University Centre for Astronomy and Astrophysics in Pune for their kind hospitality. Z.S. acknowledges the Albert Einstein Centre for Gravitation and Astrophysics supported by the Czech Science Foundation grant No. 14-37086G. N.D. thanks Albert Einstein Institute, Golm for the summer visit. B.A. is supported in part by Project No. VA-FA-F-2-008 and No. MRB-AN-2019-29 of the Uzbekistan Ministry for Innovative Development, by the Abdus Salam International Centre for Theoretical Physics and by an Erasmus+ exchange grant between SU and NUUz.