Ultrahigh energy cosmic ray nuclei from remnants of dead quasars

TL;DR

This paper explores how dormant quasars could accelerate heavy nuclei to ultrahigh energies via the Blandford-Znajek mechanism, explaining observed cosmic ray spectra and correlating remnants with Telescope Array hotspots.

Contribution

It introduces a classification of quasar remnants based on photonuclear interactions affecting cosmic ray composition and energy spectrum.

Findings

Heavy nuclei can be accelerated in dormant quasars to ultrahigh energies.

Photonuclear interactions dominate energy losses, influencing cosmic ray composition.

Five nearby quasar remnants correlate with Telescope Array hotspots.

Abstract

We re-examine the possibility of ultrahigh energy cosmic rays being accelerated in nearby dormant quasars. We particularize our study to heavy nuclei to accommodate the spectrum and nuclear composition recently reported by the Pierre Auger Collaboration. Particle acceleration is driven by the Blandford-Znajek mechanism, which wires the dormant spinning black holes as Faraday unipolar dynamos. We demonstrate that energy losses are dominated by photonuclear interactions on the ambient photon fields. We argue that the local dark fossils of the past quasar activity can be classified on the basis of how source parameters (mass of the central engine and photon background surrounding the accelerator) impact the photonuclear interaction. In this classification it is possible to distinguish two unequivocal type of sources: those in which nuclei are completely photodisintegrated before escaping the acceleration region and those in which photopion production is the major energy damping mechanism. We further argue that the secondary nucleons from the photodisintegrated nuclei (which have a steep spectral index at injection) can populate the energy region below "the ankle" feature in the cosmic ray spectrum, whereas heavy and medium mass nuclei (with a harder spectral index) populate the energy region beyond "the ankle", all the way to the high energy end of the spectrum. In addition, we show that five potential quasar remnants from our cosmic backyard correlate with the hot-spot observed by the Telescope Array.

Tables

Table 1. Table 1: Relevant source parameters and normalization factor n 0 β superscript subscript 𝑛 0 𝛽 n_{0}^{\beta} , for β = − 4 , − 5 𝛽 4 5 \beta=-4,\,-5 .

Source	$L_{\mathrm{FIR}}/L_{\odot }$	$M_9$	$n_0^{-4}/\left(\mathrm{MeV}{\mathrm{cm}}^3\right)$	$n_0^{-5}/\left(\mathrm{MeV}{\mathrm{cm}}^3\right)$	$z$	References
NGC 2768		0.16	$7.8\times {10}^{23}$	$1.0\times {10}^{24}$	$0.0051$	Torres:2002bb ; Crocker:2008dk ; Temi:2009ue
NGC 2832		11.4	$3.0\times {10}^{19}$	$3.9\times {10}^{19}$	$0.0194$	Isola:2003jk ; Temi:2009ue ; Brown:2011xh
NGC 3610		0.05	$6.9\times {10}^{23}$	$8.9\times {10}^{23}$	$0.0066$	Torres:2002bb ; Temi:2009ue ; Temi:2003sj
NGC 3613		0.16	$4.8\times {10}^{23}$	$6.2\times {10}^{23}$	$0.0066$	Torres:2002bb ; Brown:2011xh ; Temi:2003sj
NGC 4125		0.28	$1.5\times {10}^{23}$	$1.9\times {10}^{23}$	0.0063	Torres:2002bb ; Temi:2009ue ; Temi:2003sj

Equations

$$L_{\mathrm{Edd}}=1.3\times 10^{47}\frac{M}{10^9{M}_{\odot }}erg/s\sim 10^2{L}_{\mathrm{galaxy}},$$

$${\dot{M}}_{\mathrm{Edd}}=\frac{L_{\mathrm{Edd}}}{\eta c^2}\simeq 22M_9{M}_{\odot }{\mathrm{yr}}^{-1}.$$

$$B_0\sim 10^4{M}_9^{-1/2}{\left(\frac{\dot{M}}{{\dot{M}}_{\mathrm{Edd}}}\right)}^{1/2}\mathrm{G}.$$

$$|\mathbf{J}|=a{J}_{\max }=a\frac{G{M}^2}{c},$$

$$r_S=\frac{2GM}{c^2}\simeq 3\times 10^{14}\left(\frac{M}{10^9{M}_{\odot }}\right)\mathrm{cm}.$$

$$r_h=r_g\left(1+\sqrt{1-a^2}\right),$$

$$r_h<r<r_{\max }=r_g\left(1+\sqrt{1-a^2{\cos }^2\theta }\right),$$

$$|\mathbf{\Omega }|=a\left(\frac{c}{2r_h}\right).$$

$$\mathbf{E}=-\mathbf{v}\wedge \mathbf{B}.$$

$$V$$

$$V\sim 10^{20}{M}_9{B}_4\mathrm{V},$$

$$\underset{r\to \infty }{\lim }{B}_{\hat{r}}=B_0\cos \theta ,\underset{r\to \infty }{\lim }{B}_{\hat{\theta }}=-B_0\sin \theta ,$$

$$\underset{r\to \infty }{\lim }{E}_{\hat{r}}=-\frac{B_{0}aM(3{\cos }^2\theta -1)}{r^2},\underset{r\to \infty }{\lim }{E}_{\hat{\theta }}=\mathcal{O}(r^{-4}).$$

$$E_{\max }=ZeV\sim 10^{11}Z{M}_9{B}_4\mathrm{GeV}\sim 10^{11}Z{M}_9^{1/2}\mathrm{GeV},$$

$$V_{\mathrm{eff}}\sim V{\left(\frac{\zeta }{r_h}\right)}^2.$$

$${\frac{dE}{dt}|}_{\mathrm{acc}}=Ze{V}_{\mathrm{eff}}\frac{c}{\zeta },$$

$${\tau }_{\mathrm{acc}}=E{\left[{\frac{dE}{dt}|}_{\mathrm{acc}}\right]}^{-1}=\frac{E}{Ze{V}_{\mathrm{eff}}}\frac{\zeta }{c}.$$

$${\gamma }_{\max }^{\mathrm{acc}}\sim 10^{11}\frac{Z}{A}{M}_9{B}_4{\left(\frac{\zeta }{r_h}\right)}^2,$$

$${\frac{dE}{dt}|}_{\mathrm{loss}}^{\mathrm{rad}}=\frac{2}{3}\frac{Z^2{e}^{2}c}{r_c^2}{\gamma }^4,$$

$${\tau }_{\mathrm{int}}^{\mathrm{rad}}=\gamma A{m}_p{\left[{\frac{dE}{dt}|}_{\mathrm{loss}}^{\mathrm{rad}}\right]}^{-1}.$$

$${\gamma }_{\max }^{\mathrm{rad}}\sim 10^{11}{\left(\frac{B_4}{Z}\right)}^{1/4}{M}_9^{1/2}{\left(\frac{r_c}{r_h}\right)}^{1/2}{\left(\frac{\zeta }{r_h}\right)}^{1/4}.$$

$${\gamma }_{\max }=\min \left\{{\gamma }_{\max }^{\mathrm{acc}},{\gamma }_{\max }^{\mathrm{rad}}\right\}.$$

$$n(\epsilon )=n_0\{\begin{array}{ll}(\epsilon /{\epsilon }_0{)}^{\alpha }& \epsilon <{\epsilon }_0\\ (\epsilon /{\epsilon }_0{)}^{\beta }& \text{otherwise},\end{array}$$

$$L_{\mathrm{FIR}}^{\mathrm{BH}}$$

$$n_0=\frac{L_{\mathrm{FIR}}^{\mathrm{BH}}}{4\pi r_g^{2}c}{\left\{\frac{{\epsilon }_0^{-\alpha }}{(\alpha +2)}\left[{\epsilon }_0^{\alpha +2}-{\epsilon }_{\min }^{\alpha +2}\right]-\frac{{\epsilon }_0^2}{(\beta +2)}\right\}}^{-1}.$$

$$\frac{1}{{\tau }_{\mathrm{int}}^{A\gamma }}=\frac{c}{2}{\int }_0^{\infty }d\epsilon \frac{n(\epsilon )}{{\gamma }^2{\epsilon }^2}{\int }_0^{2\gamma \epsilon }d{\epsilon }^{\prime }{\epsilon }^{\prime }\sigma ({\epsilon }^{\prime }),$$

$$\sigma ({\epsilon }^{\prime })=\pi {\sigma }_{\mathrm{res}}\frac{{\Gamma }_{\mathrm{res}}}{2}\delta ({\epsilon }^{\prime }-{\epsilon }_{\mathrm{res}}^{\prime }),$$

$$\frac{1}{{\tau }_{\mathrm{int}}^{A\gamma }(E)}$$

$$\frac{1}{\tau_{\rm int}^{A\gamma}(E)}=\frac{1}{\tau_{b}}\left\{\begin{array}[]{ll}\,(E_{b}/E)^{\beta+1}&~{}E\leq E_{b}\\ (1-\beta)/(1-\alpha)\left[\left(E_{b}/E\right)^{\alpha+1}-\left(E_{b}/E\right)^{2}\right]+\left(E_{b}/E\right)^{2}&~{}E>E_{b}\end{array}\right.\,,$$

$${\tau }_b=\frac{2E_b(1-\beta )}{c\pi {\sigma }_{\mathrm{res}}A{m}_p{\Gamma }_{\mathrm{res}}{n}_0}\mathrm{and}{E}_b=\frac{{\epsilon }_{\mathrm{res}}^{\prime }A{m}_p}{2{\epsilon }_0}.$$

Full text

I Introduction

The origin(s) of ultrahigh energy cosmic rays (UHECRs) remains a central question in high energy astrophysics Torres:2004hk ; Kotera:2011cp . Three principle observables drive the search for and characterization of cosmic ray sources: the energy spectrum, the nuclear composition, and the distribution of arrival directions. The cosmic ray energy spectrum encompasses a plummeting flux that drops from $10^{4}~{}{\rm m}^{-2}\,{\rm s}^{-1}$ at $E\sim 1~{}{\rm GeV}$ to $10^{-2}~{}{\rm km}^{-2}\,{\rm yr}^{-1}$ at $E\sim 10^{11}~{}{\rm GeV}$. Its shape is remarkably featureless and can be described by a broken power law: the spectrum steepens gradually from $E^{-2.7}$ to $E^{-3.3}$, then flattens to $E^{-2.6}$ near $10^{9.6}~{}{\rm GeV}$ forming a feature known as “the ankle”, and drops sharply at around $10^{10.5}~{}{\rm GeV}$ Abbasi:2007sv ; Abraham:2008ru ; Abraham:2010mj . The small variations of the spectral index can be interpreted either as a transition between cosmic ray populations, or as an imprint of cosmic ray propagation effects.

The simplest interpretation of the ankle is that above $10^{9.6}~{}{\rm GeV}$ a new population emerges which dominates the more steeply falling Galactic population of heavy nuclei. The extragalactic component can be dominated either by protons Bird:1993yi or heavies Allard:2005ha ; Allard:2005cx . Proton-dominance beyond the ankle is ultimately limited by the onset of photopion production on the cosmic microwave background, whereas dominance of a heavy composition is restricted by nucleus photodisintegration through the giant dipole resonance – the so called Greisen-Zatsepin-Kuz’min (GZK) suppression at around $10^{10.5}~{}{\rm GeV}$ Greisen:1966jv ; Zatsepin:1966jv . It has also been advocated that the ankle feature could be well reproduced by a proton-dominated power-law spectrum, where the ankle is formed as a dip in the spectrum from the energy loss of protons via Bethe-Heitler pair production Hillas:1967 ; Berezinsky:2002nc . In this case extra-galactic protons would already have started to dominate the spectrum somewhat beyond $10^{8.7}~{}{\rm GeV}$.

Optical observations of air showers with fluorescence telescopes or non-imaging Cherenkov detectors consistently find a predominantly light composition at around $10^{9}~{}{\rm GeV}$ Kampert:2012mx and the contribution of protons to the overall cosmic ray flux is $\gtrsim$ 50% in this energy range Aab:2014kda ; Aab:2014aea ; Aab:2016htd ; Aab:2016zth . Due the absence of a large anisotropy in the arrival direction of cosmic rays below the ankle Auger:2012an ; ThePierreAuger:2014nja , we can conclude that these protons must be of extra-galactic origin. At energies above $10^{10}~{}{\rm GeV}$, the high-statistics data from the Pierre Auger Observatory suggests a gradual increase of the fraction of heavy nuclei in the cosmic ray flux Aab:2014aea ; Aab:2014kda ; Aab:2016htd ; Aab:2016zth . Within uncertainties, the data from the Telescope Array (TA) is consistent with these findings Abbasi:2014sfa ; Abbasi:2015xga . Moreover, the Telescope Array has observed a statistically significant excess in cosmic rays with energies above $10^{10.7}~{}{\rm GeV}$ in a region of approximately 1150 square degrees centered on equatorial coordinates (R.A. = $146.7^{\circ}$, Dec. = $43.2^{\circ}$) – the TA hot spot Abbasi:2014lda . The absence of a concentration of nearby sources in this region of the sky favors a heavy nucleus hypothesis, whereby a few local sources within the GZK sphere can produce the hot spot through deflection in the Galactic ${\bm{B}}$-field. All in all, the apparent dominance of heavy nuclei in the vicinity of $10^{10}~{}{\rm GeV}$ strongly argues against the interpretations of the ankle given above.

One can (of course) accommodate the data with the addition of an ad hoc light extragalactic component below the ankle, with a steep injection spectrum Aloisio:2013hya . However, a more natural explanation of the entire spectrum and composition emerges while accounting for the “post-processing” of UHECRs through photodisintegration in the environment surrounding the source Unger:2015laa ; Farrar:2015ikt ; Globus:2015xga ; Anchordoqui:2014pca . In such models relativistic nuclei accelerated by a central engine to extremely high energies remain trapped in the turbulent magnetic field of the source environment. Their escape time decreases faster than the interaction time with increasing energy, so that only the highest energy nuclei can escape the source unscathed. In effect, the source environment acts as a high-pass filter on the spectrum of cosmic rays. All nuclei below the energy filter interact, scattering off the far-infrared photons in the source environment. These photonuclear interactions produce a steep spectrum of secondary nucleons, which is overtaken by the harder spectrum of the surviving nucleus fragments above about $10^{9.6}~{}{\rm GeV}$. These overlapping spectra could then carve an ankle-like feature into the source emission spectrum. The spectrum above the ankle exhibits a progressive transition to heavy nuclei, as the escape of non-interacting nuclei becomes efficient. In such models, the source evolution with redshift plays a paramount role on the source parameters. For models with positive redshift evolution (a greater number of sources per comoving volume at high redshift), both the imposition of a dominant light extragalactic component below the ankle and the source-ankle model favor a hard injection spectrum $\propto E^{-1}$. However, it has been recently noted that fitting the high end of the cosmic ray (CR) spectrum and composition with negative source evolution allows for softer injection spectra Taylor:2015rla .

On a separate track, it was recently pointed out that Fermi-LAT observations Ackermann:2014usa ; TheFermi-LAT:2015ykq severely constrain the distribution of UHECR sources Liu:2016brs . This is because on the way to Earth UHECRs interact with the radiation fields permeating the universe and give rise to energetic electron/positron pairs and photons, which in turn feed electromagnetic cascades resulting in a diffuse gamma-ray background radiation. The recent study in Liu:2016brs suggests that a local “fog” of UHECRs accelerated in nearby sources must exist so as to not overproduce the cascade gamma-ray flux beyond Fermi-LAT observations.

In consideration of these new results, it seems interesting to explore whether remnants of dead quasars, which are distributed preferentially in the low-redshift universe, could accelerate heavy nuclei up to the highest observed energies. Actually, nearby quasar remnants have long been suspected to be sources of UHECRs Boldt:1999ge ; Boldt:2000dx ; Torres:2002bb ; Neronov:2007mh . Particle acceleration proceeds via the Blandford-Znajek mechanism, which wires dormant spinning black holes as Faraday unipolar dynamos Blandford:1977ds ; Znajek . Detailed modeling with reasonable assumptions on source parameters suggests that protons can be efficiently launched at $E\gtrsim 10^{12}~{}{\rm GeV}$. Photopion production in collisions with ambient photons becomes a relatively important effect during the final phase of the acceleration process, and could damp the maximum attainable energy to values well below that imposed by the voltage drop. The large charge $Ze$ of heavy nuclei facilitates their acceleration to highest energy (scaling linearly with $Z$). However, heavy nuclei could also interact with the ambient photons and may suffer significant spallation. In this paper we study in detail the dominant processes for nucleus energy losses, while searching for plausible source parameters which may allow acceleration of heavy nuclei up to the highest observed energies.

The layout of the paper is as follows. In Sec. II we provide an assessment of what is known observationally about the evolution of quasar activity and then discuss the general framework use to relate the quasar population with massive dark objects at the center of local inactive galaxies. We also show that five nearby relics of such past activity correlate with the TA hot-spot. In Sec. III we review the generalities of particle acceleration near the event horizon of a spinning supermassive black hole and determine plausible source conditions for acceleration of UHECR nuclei. In Sec. IV we examine the various energy-loss mechanisms of one-shot acceleration, including radiative processes and photonuclear interactions with the ambient photons. We show that, in general, photonuclear energy losses are dominant. After that we propose a novel phenomenological model to explain the evolution of the UHECR spectrum and composition with energy, without fine-tuning. Finally, in Sec. V we present our conclusions.

II Dormant black holes after a shining past

The first extensive radio surveys of the 1950s revealed a new family of radio sources that, by the early 60s, had been associated with point-like optical counterparts Hazard . These quasi-stellar radio sources, or quasars, revealed very large systemic redshifts indicating cosmological distances Schmidt . This implied luminosities up to hundreds of times that of the Milky Way, apparently emitted from very small size-scales Smith , suggesting an extremely energetic emission mechanism. It has now been established that quasars are ultimately powered by the accretion of matter onto the $10^{7}-10^{10}M_{\odot}$ supermassive black holes (SMBHs) at cores of galaxies Zeldovich ; Salpeter:1964kb ; LyndenBell:1969yx ; Antonucci:1993sg ; Urry:1995mg .

Quasars and their relativistically-beamed counterparts, blazars, are the most luminous sub-class of active galactic nuclei (AGNs), which is the general name for an actively-accreting SMBH. In fact, these objects are the most luminous continuously-emitting objects in the Universe. Quasar luminosity is, to first order, constrained by the Eddington limit Eddington . This is the maximum luminosity of a spherically accreting object, determined by the point at which outward radiation pressure halts gravitational infall. The Eddington limit is proportional to the mass of the accreting object. Scaled for a quasar powered by a SMBH of $10^{9}~{}M_{\odot}$, this luminosity is:

[TABLE]

where $M_{\odot}$ is the solar mass and $L_{\rm galaxy}$ is the bolometric luminosity of large galaxies like the Milky Way. Note that accretion by the standard $\alpha$-disk, which is typically assumed for quasars, is also Eddington limited Shakura:1972te .

It is generally accepted that the number density of quasars peaked at redshift $z\sim 2$, when the universe was one-fifth of its present age Rees . Locally, essentially all galaxies with spheroidal components contain a SMBH Richstone:1998ky ; Ho:1998hz ; Kormendy:2004tc . The close relationships between the masses of these SMBHs and their host galaxy’s properties, including stellar velocity dispersion, spheroidal luminosity, and dynamical mass, indicates that the growth of SMBHs and their hosts are closely tied Kormendy:1995er ; Magorrian:1997hw ; Gebhardt:2000fk ; Ferrarese:2000se ; Haring:2004hr . These relationships also allow us to make a relatively direct determination of the local SMBH mass function Vika:2009ef ; McLure:2003em ; Graham:2007qb ; Marconi:2003tg ; Tamura:2005kc .

Another approach to deriving the local SMBH mass function is to integrate SMBH accretion over cosmic time based on the evolving high-$z$ quasar luminosity function Marconi:2003tg ; Tamura:2005kc ; Soltan:1982vf ; Chokshi ; Small92 ; Cao:2008pd ; Li:2011ik ; Shankar:2007zg . Such models are able to reproduce the “observered” low-$z$ mass function for quiescent SMBHs very successfully by assuming the typical radiative efficiencies for an $\alpha$-disk ($0.1\lesssim\eta\lesssim 0.2$), and reasonable Eddington fractions — quasar luminosities as a fraction of the Eddington limit (Eddington fraction) — of $0.2\lesssim L/L_{Edd}\lesssim 1$ Graham:2007qb ; Marconi:2003tg ; Cao:2008pd ; Li:2011ik ; Yu:2002sq . This general agreement supports for the idea that accretion in AGN mode is an important, and perhaps the principle mechanism by which local SMBHs grow Rees ; Small92 ; Marconi:2003tg ; Cao:2008pd . Thus, it is reasonable to assume that many of the more massive quiescent SMBHs in the local universe are ”remnant” quasars Rees ; Chokshi , and that these greatly outnumber active quasars in the modern epoch.

In their active phases, accreting SMBHs are expected to support strong magnetic fields. Assuming a rapidly-rotating black hole, which may be common in luminous quasars Cao:2008pd ; Volonteri:2007tu , this magnetic field is expected to survive the expiration of the quasar phase. We calculate a typical relic magnetic field strength for a rapidly rotating SMBH by assuming that it corresponds to maximum field strength acquired during its quasar phase. This in turn corresponds to its maximum accretion rate. To determine this accretion rate we assume the lower end of $\alpha$-disk radiative efficiency, $\eta=0.1$ (corresponding to a higher accretion rate per unit luminosity), and that the quasar radiates at the Eddington limit. The mass accretion rate needed to sustain the Eddington luminosity is

[TABLE]

where $M_{9}=M/(10^{9}M_{\odot})$. The accreting plasma is assumed to support an axisymmetric magnetic field configuration due to the generation of currents. Following Boldt:1999ge , we get an estimate of the characteristic magnetic field strength $B_{0}$ by assuming pressure equilibrium between the magnetic field and the in-falling matter Anchordoqui:2016dcp

[TABLE]

This magnetic field may be taken as a fiducial strength for quasars remnants at the upper end of the local SMBH mass and spin functions.

This represents the magnetic field strength developed by a maximally rotating black hole at the peak of its accretion phase and at the upper end of the mass function. In the following calculations we work with this fiducial strength, yielding. It will be seen that a relatively high field strength, $B_{0}\sim 10^{4}~{}{\rm G}$, is necessary to achieve a significant flux of $10^{11}~{}{\rm GeV}$ single-proton CRs around a $10^{9}\,M_{\odot}$ SMBH. However, when we consider black holes up to $10^{10}\,M_{\odot}$, and heavy nuclei CRs, the required field strength drops by up to two orders of magnitude. We also assume that a significant fraction of the maximum field strength is sustained as accretion rate drops into the remnant quasar phase. Further data on the spatial association of UHECRs with SMBHs will test this assumption.

Significant magnetic fields do appear to survive in the case of least one class of remnant quasar. BL Lacertae objects are AGNs of the blazar class – broadly, AGNs with jets that are closely aligned with our line of sight Urry:1995mg ; Blandford:1978 . BL Lacs show little or no signs of broad emission lines or thermal emission Urry:1995mg ; Marcha ; Giommi:2011sn . This implies that at least a subset of BL Lacs are starved of gas, and perhaps in a post-quasar phase. Under this interpretation, flat-spectrum radio quasars (FSRQs) become the progenitor population of low-accretion BL Lacs Cavaliere:2001gb ; Boettcher:2001jp ; Ajello:2011zi . This evolutionary argument is supported by BL Lacs’ negative cosmological evolution; their number densities rise at around the same epoch as the decline of FSRQs Ajello:2013lka . BL Lacs exhibit BZ-powered jets Maraschi:2002pp , indicating that significant magnetic fields may be sustained long after the cessation of significant accretion.

Hanny’s Voorwerp contains perhaps the first direct probe of quasar evolution Lintott:2009wq . Observations in the optical, ultraviolet, and X-ray indicate that this unusual object near the spiral galaxy IC 2497 contains highly ionized gas Jozsa:2009qm . The emission-line properties, and lack of X-ray emission from IC 2497, seem to indicate this spiral galaxy is a highly obscured AGN with a novel geometry arranged to allow photoionization of Hanny’s Voorwerp but not the galaxy’s own circumnuclear gas Rampadarath:2010yv , or alternatively the luminosity of the central source has decreased dramatically within the last $10^{5}$ yr Schawinski:2010up . Should this be the case, Hanny’s Voorwerp may exemplify the first detection of a quasar light echo. Interestingly, IC 2497 is inside the TA hot spot; see Fig. 1. However, the measured redshift $z=0.05$ of IC 2497, sets this dormant quasar outside our local GZK–horizon ($\sim 100$ Mpc). To sample sources inside the GZK–sphere we adopt a distance-limited ($z\lesssim 0.02$) compilation of quasar remnant candidates, capable of accelerating UHECRs Torres:2002bb ; Isola:2003jk . It is compelling that five of the objects contained in this sample are also inside the TA hot spot; see Fig. 1.

Centaurus A (Cen A) is a complex radio-loud source identified at optical frequencies with the galaxy NGC 5128. Cen A is the closest radiogalaxy to Earth (distance $\sim 3.4~{}{\rm Mpc}$) and has long been suspected to be a potential UHECR accelerator Cavallo:1978 ; Romero:1995tn . The Pierre Auger Collaboration has searched for anisotropies in the direction of Cen A scanning the energy threshold between $10^{10.6}~{}{\rm GeV}$ and $10^{10.9}~{}{\rm GeV}$ and counting events in angular radii ranging from $1^{\circ}$ to $30^{\circ}$ PierreAuger:2014yba . The strongest departure from isotropy (post-trial probability $\sim 1.4\%$) has been observed for $E>58~{}{\rm EeV}$ in a window of $15^{\circ}$: 14 events (out of a total of 155) have been observed in such an angular window while 4.5 are expected on average from isotropic distributions.

A statistical analysis to ascertain the significance of a possible correlation between UHECRs and SMBHs would require using data from Auger and TA. However, if the high energy end of the cosmic ray spectrum is the same in the north and in the south, there appears to be a difference in the energy scale between the two detectors Dawson:2013wsa , making such study unreliable at present. Future experiments that will observe both hemispheres using the same apparatus, such as the Probe Of Extreme Multi-Messenger Astrophysics (POEMMA), would make a robust statistical analysis more straightforward.111POEMMA has been recently selected by NASA for an in-depth probe mission concept study in preparation for the next deceadal survey Adams:2017fjh .

For convenience, we subdivide the SMBH sample into two subcomponents: one with masses $\gtrsim 10^{10}M_{\odot}$ and one with masses $\lesssim 10^{10}M_{\odot}$. The rationale for using $10^{10}M_{\odot}$ as the demarcation mass will become evident in Sec. IV. Roughly speaking, 16% of SMBHs have mass above $10^{10}M_{\odot}$ Isola:2003jk . Inside the TA hot spot one of six (five nearby) candidate sources is consistent with a mass$\gtrsim 10^{10}M_{\odot}$. Thus, we can assume that the ratios inside the hot spot are representative of the ratios of the universe at large.

III UHECR acceleration in quasar remnants

In realistic astrophysical situations involving quasars, the black hole itself is uncharged, and the gravity of the accretion disk is practically negligible. This means that the geometry of spacetime is described by the Kerr metric and therefore determined by two free parameters: the black hole mass $M$ and its angular momentum

[TABLE]

where $0\leq a\leq 1$ is the dimensionless spin parameter Kerr:1963ud . Before proceeding we note that for a Schwarzschild black hole, $a=0$ Schwarzschild:1916uq , whereas for a maximally spinning (a.k.a. extreme Kerr) black hole, $a=1$. As a matter of fact, the increase of the spin parameter would stop at $a\approx 0.998$ because photons emitted from the disk on retrograde paths are more likely to be captured by the hole than prograde ones, hence de-spinning the hole Thorne:1974ve .

For static black holes, the scale characterizing the event horizon $r_{h}$ is the Schwarzschild radius

[TABLE]

For spinning black holes, the (spherical) event horizon surface is defined by

[TABLE]

where $r_{g}=r_{S}/2$ is the gravitational radius. The ergosphere is the region described by the relation

[TABLE]

where $\theta$ is the angle to the polar axis; see e.g. Visser:2007fj ; Rieger:2011ch . Inside the ergosphere (which is ellipsoidal in shape), spacetime is pulled heavily towards the direction of the black hole rotation (frame dragging), and consequently no static observer can exist because the particles must co-rotate with the hole. The black hole’s angular velocity $\bm{\Omega}$ coincides with the angular velocity of the dragging of inertial frames at the horizon Rieger:2011ch

[TABLE]

It has long been known that a rotating black hole can radiate away its available reducible energy Penrose:1969pc ; Penrose:1971uk . In this context Blandford and Znajek (BZ) proposed a model of electromagnetic extraction of black hole’s rotational energy based on the analogy with the classical Faraday (unipolar induction) dynamo phenomenon Blandford:1977ds ; Znajek . In the BZ mechanism, the magnetic field $\bm{B}$ near the horizon taps the rotational energy of the black hole and generates powerful outflows of electromagnetic (Poynting) energy. BZ argue that spacetime frame dragging induces an electric field $\bm{E}$ that is strong enough to break the vacuum and establish an electron-positron force-free magnetosphere.

The particulars of the black hole magnetosphere could be very complicated, distinctively if the poloidal magnetic field near the horizon is misaligned with the black hole rotation axis Neronov:2007mh . For simplicity, herein we assume that the $\bm{B}$-field is aligned with the axis of rotation. This manageable system allows a transparent analytic calculation which captures the essence of most aspects of the BZ mechanism and leads to a correct order of magnitude.

Conduction electrons within the magnetosphere undergo collisions with the orbiting atoms, and this causes them to take up the rotational motion. At each point the electrons have a net drift velocity, ${\bf v}=\bm{\Omega}\wedge\bm{r}$, where $\bm{r}$ is the position vector of the point in question relative to an origin which is chosen to lie on the magnetic axis. The conduction electrons experience a magnetic force $-e{\bf v}\wedge\bm{B}$, which is mainly directed towards the central axis; as a result a negative charge appears in the core of the magnet, and a positive charge on its curved outer surface. In equilibrium an electrostatic field is set up, such that the total Lorentz force on the conduction electrons is zero. Namely, $\bm{F}=-e\bm{E}-e{\bf v}\wedge\bm{B}=\bm{0}$, and so

[TABLE]

We now calculate the potential difference $V$ corresponding to this electric field. It is straightforward to see by inspection of (9) that along the axis of rotation $\bm{E}=0$, and therefore $V=0$ at all points along this axis. For a maximally spinning hole, the potential drop between the central ($r=0$) and the marginal ($r=r_{\rm max}$) magnetic surfaces passing through the horizon is found to be

[TABLE]

where $B_{0}$ is the magnetic field strength. Now, recalling that $10^{4}~{}{\rm G}=1~{}{\rm T}=1~{}{\rm V\,s\,m^{-2}}$, we can rewrite (10) as

[TABLE]

where $B_{4}=B_{0}/10^{4}~{}{\rm G}$.

Far away from the horizon, the electromagnetic field can be expressed as Wald:1974np

[TABLE]

and

[TABLE]

The electric field has a quadrupole topology, which is distinctly seen in Fig. 2. The voltage difference between the horizon and $r=\infty$ is of order $V$. With such a huge voltage drop along field lines, one expects that electrons roaming in the vicinity of the horizon would be accelerated to huge energies and, upon colliding with stray photons and/or positrons, produce a cascade of $e^{+}e^{-}$ pairs. Very quickly therefore, the vacuum surrounding the black hole would be filled with a highly conducting plasma. The plasma-loaded field lines can serve as wires to complete a circuit between the pole and equator of the horizon. The electromotive force around this circuit is numerically comparable to (11). Thus, if a cosmic ray baryon can fully tap this potential acceleration up to extreme energies,

[TABLE]

would become possible. However, the charge density in the vicinity of accreting black holes could be so high that a significant fraction of this potential would be screened and so no longer available for particle acceleration. Therefore, it seems more appropriate to define an effective potential where the available length scale, the gap height $\zeta$, is explicitly taken into account. Following Rieger:2011ch , we take

[TABLE]

Accordingly, the characteristic rate of energy gain is found to be

[TABLE]

and so the acceleration timescale is given by

[TABLE]

In a realistic situation the energy losses of the accelerated particles limit the maximal energies to the values below the estimate of (14). It is this that we now turn to study. Before proceeding though, it is interesting to note that the axisymmetric stationary flow in the vicinity of the rotation-powered compact object could provide a source of ultrarelativistic nuclei with an extremely flat injection spectrum (typically $\propto E^{-1}$ Blasi:2000xm ; Kotera:2015pya ), in good agreement with the requirements to accommodate Auger observations Aab:2016zth ; Unger:2015laa . Moreover, a steepened injection spectrum would arise quite naturally if the medium surrounding the compact object were leaky at these late ages.

IV Mechanisms of energy dissipation

In the absence of energy losses the maximum Lorentz factor is given by

[TABLE]

where $E=\gamma Am_{p}$, with $m_{p}$ the proton mass and $A$ the nuclear baryon number. Within the potential drop, the nuclei follow the curved magnetic field lines and so emit curvature-radiation photons. The energy loss rate or total power radiated away by a single cosmic ray is Ochelkov

[TABLE]

where $r_{c}$ is the curvature radius of the magnetic field lines, and $\gamma$ is the Lorentz factor of the radiating particles. The characteristic cooling timescale is then given by

[TABLE]

Acceleration gains are balanced by radiative losses. In the absence of other damping mechanisms, the radiation reaction limit is Levinson:2000nx ; Levinson:2002ea

[TABLE]

All in all, direct electric field acceleration on the black hole magnetosphere would allow Lorentz factors up to

[TABLE]

Note that for our fiducial values $B_{4}\sim 1$, $M_{9}\sim 1$, and $r_{c}\sim r_{h}$, we have $\gamma_{\rm max}^{\rm acc}\sim\gamma_{\rm max}^{\rm rad}$ if $\zeta\sim r_{h}$. However, by taking $\zeta<r_{h}$ one can always suppress the emission of curvature photons. Moreover, as we show in what follows photonuclear interactions on the ambient photon background become the dominant mechanism for energy losses.

The estimation of the photon density in the vicinity of the SMBH is not straightforward. Studies for Cen A are quite extensive Whysong:2002ks ; Meisenheimer:2007zv , and the photon spectrum around the SMBH is consistent with a thermal origin: mostly dust emission heated by the nucleus Whysong:2002ks . For simplicity, in our calculations we approximate the photon spectrum as a broken power-law,

[TABLE]

which is also consistent with the data. Here $\varepsilon$ is the photon energy and the maximum of the density is at an energy of $\varepsilon_{0}$. The broken power-law spectrum allows a complete analytic treatment of photonuclear interactions, and the global result does not depend on the exact shape of the photon spectrum. Namely, the interaction times are comparable if the photon density is assumed to follow a (modified) black body spectrum Unger:2015laa . For more distance sources, it is not trivial to disentangle which photons are near the SMBH. As a first approximation, we use the galaxy-wide average far infra-red (FIR) flux reported in the NASA extragalactic database to normalize the spectrum and estimate the spectral indices NED . If a correlation becomes evident in the future a more detailed analysis would be required. By averaging over the thermal spectra of the sources studied in Wiklind we take $\alpha=2$ and $-5\leq\beta\leq-4$.

We normalize the spectrum to the FIR luminosity in the vicinity of the SMBH

[TABLE]

where $\varepsilon_{\rm min}=10^{-4}~{}{\rm eV}$ and $\varepsilon_{0}=10^{-2}~{}{\rm eV}$. There is a general consensus that the extended FIR emission along the dust lane of Cen A, $L_{\rm FIR}\sim 2\times 10^{10}L_{\odot}$, cannot only be attributed to dust heated by the active nucleus Marshall:1988 ; Eckart ; Wild:2000hq ; Stawarz:2006ku . Optical detections of luminous HII regions in relatively unoscured regions of the dark lane Graham:1981 suggest that recently formed stars embedded in the dust lane must be responsible for most of the broadly extended FIR emission. Only about 1% of the FIR emission comes from the active nucleus. With this in mind, we take $L_{\rm FIR}^{\rm BH}/L_{\rm FIR}\sim 10^{-2}$. The normalization factor is then

[TABLE]

In Table 1 we summarize the relevant properties of a subgroup of 5 candidate sources shown in Fig. 1. In the case of no-detection in the FIR we normalize the photon flux to the reported upper-limit. To determine the distance to these sources we adopt the usual concordance cosmology of a flat universe dominated by a cosmological constant, with $\Omega_{\Lambda}=0.692\pm 0.012$ and a cold dark matter plus baryon component $\Omega_{m}=0.308\pm 0.012$; the Hubble parameter as a function of redshift is given by $H^{2}(z)=H_{0}^{2}[\Omega_{m}(1+z)^{3}+\Omega_{\Lambda}]$, normalized to its value today, $H_{0}=100~{}h\,{\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1}$, with $h=0.678$ Olive:2016xmw .

To characterize the population of nearby quasar remnants we also study the case of NGC 5128. This case is special, as the photon spectrum peaks in the mid-infrared. We take $\epsilon_{0}=0.13~{}{\rm eV}$ and $L_{\rm mid-IR}/L_{\odot}=1.3\times 10^{8}$. The mass of the SMBH is somewhat uncertain Cappellari:2008db ; Neumayer:2007ic ; Marconi:2005qx . In our calculations we adopt $M_{9}=0.1$, yielding $n_{0}^{-4}=3.0\times 10^{23}~{}{\rm MeV^{-1}\,cm^{-3}}$ and $n_{0}^{-5}=3.9\times 10^{23}~{}{\rm MeV^{-1}\,cm^{-3}}$.

The interaction time for a highly relativistic nucleus propagating through an isotropic photon background with energy $\varepsilon$ and spectrum $n(\varepsilon)$, normalized so that the total number of photons in a box is $\int n(\varepsilon)d\varepsilon$, is given by Stecker:1969fw

[TABLE]

where $\sigma(\varepsilon^{\prime})$ is the photonuclear interaction cross section of a nucleus by a photon energy $\varepsilon^{\prime}$ in the rest frame of the nucleus.

We have found that for the considerations in the present work, the cross section can be safely approximated by the single pole of the narrow-width approximation,

[TABLE]

where $\sigma_{\rm res}$ is the resonance peak, $\Gamma_{\rm res}$ its width, and $\varepsilon^{\prime}_{\rm res}$ the pole in the rest frame of the nucleus. The factor of $1/2$ is introduced to match the integral (i.e. total cross section) of the Breit-Wigner and the delta function Anchordoqui:2006pe .

The mean interaction time can now be readily obtained substituting Eq. (27) into Eq. (26),

[TABLE]

Substituting (23) into (28) we finally find Unger:2015laa

[TABLE]

where

[TABLE]

The parameters characterizing the photodisintegration cross section are: $\sigma_{\rm res}\approx 1.45\times 10^{-27}~{}{\rm cm}^{2}\,A$, $\Gamma_{\rm res}=8~{}{\rm MeV}$, and $\epsilon^{\prime}_{\rm res}=42.65A^{-0.21}\,(0.925A^{2.433})~{}{\rm MeV},$ for $A>4$ ($A\leq 4$) Karakula:1993he . The parameters for the photopion production cross section are: $\sigma_{\rm res}\simeq 5.0\times 10^{-28}~{}{\rm cm}^{2}\,A$, $\Gamma_{\rm res}=150~{}{\rm MeV}$, and $\varepsilon^{\prime}_{\rm res}=(m_{\Delta}^{2}-m_{p}^{2})/(2m_{p})\simeq 340~{}{\rm MeV}$ Olive:2016xmw .

In Figs. 3, 4, 5, 6, 7, and 8 we compare the characteristic time scales of acceleration and energy losses for the various representative sources given in Table 1 and Cen A. We consider 56Fe and 28Si as fiducial nuclear species. In general, photonuclear interactions dominate the energy losses. By comparing the results for $\beta=-5$ and $\beta=-4$, we can see that there is almost no dependence on the photon spectral index. For $E\gtrsim 10^{11}~{}{\rm GeV}$, the accelerated nuclei reach the photopion production threshold, and the significance of this process becomes more important with increasing energy. The photopion production threshold corresponds to an energy-per-nucleon of $E/A\sim 10^{10}~{}{\rm GeV}$. Note that the onset of photopion production is above the required maximum energy $=10^{9.5}~{}{\rm GeV}$ for the fiducial model in Unger:2015laa . Therefore, we conclude that some quasar remnants are capable of launching heavy (and medium mass) nuclei up to the observed maximum energies.

Though Figs. 3, 4, 5, 6, 7, and 8 provide only a naïve illustration of particle acceleration with energy dissipation in the black hole dynamo it is evident that it possible to classify the one-shot acceleration mechanism in quasar remnants according to how source parameters (mass of the central engine and the ambient photon backgrounds) impact the photonuclear interaction. In the first type nuclei are completely photodisintegrated before escaping the acceleration region producing a flux of secondary nucleons, e.g. NGC 2768, 3610, 3613, 4125, and 5128. In the second type photopion production is the major energy damping mechanism and nuclei are able to escape without suffering significant spallation, e.g. NGC 2832 (which is near the center of the TA hot-spot). Combining the ideas put forward in Aloisio:2013hya ; Unger:2015laa we can now formulate a new hybrid model to explain the spectral shape of extragalactic cosmic rays, including the critical region of the ankle. Namely, the secondary nucleons produced in the photodisintegration process would have a soft spectral index at the source, and consequently can explain the energy region below the ankle. Note that after the photodisintegration process the secondary protons can still be accelerated to reach energies near the ankle. On the other hand, heavy and medium mass nuclei are emitted with a harder spectral index, and therefore (upon propagation to Earth) can populate the spectrum above the ankle, all the way to the GZK cutoff. We estimate that about 16% of the SMBH in the universe could accelerate UHECR nuclei via the BZ. At this stage, it is worthwhile to point out that for SMBHs associated with jets terminating in lobes (which can be detected in radio) like Cen A other cosmic ray acceleration mechanisms may be at play Biermann:1987ep ; Rachen:1992pg ; Rieger:2011rh .

It is also possible that the emitted nuclei could suffer additional photodisintegration while diffusing in the source environment, as described in Unger:2015laa . Photodisintegration of high-energy nuclei, $A+\gamma\to A^{\prime*}+X$, is followed by immediate photoemission from the excited daughter nuclei, $A^{\prime*}\to A^{\prime}+\gamma$ Anchordoqui:2006pd . The photodisintegration process in the source environment would then produce a flux of gamma-rays. A suggestion of a correlation between UHECRs and nearby gamma-ray emitting AGNs has been put forward in Nemmen:2010bp . Such a correlation would be consistent with our model. However, we should exercise caution since the results in Nemmen:2010bp are given with an a posteriori significance, and since the new AGN correlation analyses has not shown an increased in significance over time PierreAuger:2014yba .

In closing, we study the dependence of our conclusions on the magnetic field strength. As we have shown, the SMBHs which are capable of emitting UHECR nuclei without suffering significant spallation are those with masses $\gtrsim 10^{10}M_{\odot}$. We can then rewrite the acceleration time scale (17) keeping $B_{4}$ as a variable, i.e. we take

[TABLE]

In Fig. 9 we compare the acceleration and interaction time scales for different magnetic field strengths. As we anticipated in Sec. II, by simple inspection of Fig. 9 one can see that CR nuclei can be accelerated to the highest observed energies with $B_{4}\gtrsim 0.01$.

V Conclusions

We have investigated the possibility of UHECRs being accelerated in nearby supermassive black holes, which are dormant remnants of previously active quasars. We have shown that these fast spinning objects are latent dynamos that can accelerate heavy nuclei up to the maximum observed energies. We have also shown that energy losses are dominated by photonuclear interactions on the ambient photon fields. Armed with our findings we postulated a classification scheme for giant elliptical galaxies harboring supermassive black hole dynamos. In this classification we distinguish among two unequivocal type of cosmic ray sources according on how source parameters (mass of the central engine and the ambient photon backgrounds) impact the photonuclear interaction. In the first type nuclei are completely photodisintegrated before escaping the acceleration region producing a flux of secondary nucleons. In the second type photopion production is the major energy damping mechanism and nuclei are able to escape without suffering significant spallation. In the spirit of Aloisio:2013hya ; Unger:2015laa we conjectured that the spectral shape of extragalactic cosmic rays, including the pivotal region of the ankle, can be described by a smooth superposition of these two type of sources. Namely, the secondary nucleons produced in the photodisintegration process would have a soft spectral index at the source, and consequently can explain the energy region below the ankle. On the other hand, heavy and medium mass nuclei are emitted with a harder spectral index, and therefore (upon propagation to Earth) can populate the spectrum above the ankle, all the way to the GZK cutoff.

Acknowledgements.

We would like to acknowledge many useful discussions with our colleagues of the Pierre Auger Collaboration. We also thank Andrii Neronov for permission to reproduce Fig. 2. R.J.M., R.A.C., and J.J.G. are supported by the U.S. National Science Foundation (NSF) Grant No. AST-1153335. L.A.A. is supported by NSF Grant No. PHY-1620661 and by the National Aeronautics and Space Administration (NASA) Grant No. NNX13AH52G; he thanks the Department of Physics at UNLP for its hospitality. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the NSF or NASA.