Numerical models of neutrino and gamma-ray emission from magnetic reconnection in the core of radio-galaxies

TL;DR

This paper models neutrino and gamma-ray emissions from magnetic reconnection in radio-galaxies' cores, suggesting a partial explanation for VHE data and predicting neutrino fluxes below IceCube detection levels.

Contribution

It introduces a combined numerical approach to simulate particle acceleration and emission in the accretion flow of radio-galaxies, specifically applying it to Centaurus A.

Findings

CR injection can partially reproduce VHE data

Neutrino emission peaks at ~10^16 eV

Neutrino flux is two orders below IceCube limits

Abstract

Non-blazar radio-galaxies emitting in the very-high-energy (VHE; >100 GeV) regime offer a unique perspective for probing particle acceleration and emission processes in black hole (BH) accretion-jet systems. The misaligned nature of these sources indicates the presence of an emission component that could be of hadronic origin and located in the core region. Here we consider turbulent magnetic reconnection in the BH accretion flow of radio-galaxies as a potential mechanism for cosmic-ray (CR) acceleration and VHE emission. To investigate if this scenario is able to account for the observed VHE data, we combine three numerical techniques to self-consistently model the accretion flow environment and the propagation of CRs plus electromagnetic cascades within the accretion flow zone. Here we apply our approach to the radio-galaxy Centaurus A and find that injection of CRs consistent with magnetic reconnection power partially reproduce the VHE data, provided that the accretion flow makes no substantial contribution to the radio-GeV components. The associated neutrino emission peaks at $\sim10^{16}$ eV and is two orders of magnitude below the minimum IceCube flux.

Tables

Table 1. Table 1: Model parameters of the calculated VHE SED shown in Fig. 4. All the models correspond to a mass accretion rate of ⟨ M ˙ a ​ c ​ c ⟩ = 1.3 × 10 − 3 delimited-⟨⟩ subscript ˙ 𝑀 𝑎 𝑐 𝑐 1.3 superscript 10 3 \langle\dot{M}_{acc}\rangle=1.3\times 10^{-3} M ⊙ yr -1 and a CR injection spectrum d ​ N C ​ R / d ​ ε ∝ ϵ − κ ​ exp ⁡ { − ε / ε c ​ u ​ t } proportional-to 𝑑 subscript 𝑁 𝐶 𝑅 𝑑 𝜀 superscript italic-ϵ 𝜅 𝜀 subscript 𝜀 𝑐 𝑢 𝑡 dN_{CR}/d\varepsilon\propto\epsilon^{-\kappa}\exp\{-\varepsilon/\varepsilon_{cut}\} , with ε m ​ i ​ n = 10 13 subscript 𝜀 𝑚 𝑖 𝑛 superscript 10 13 \varepsilon_{min}=10^{13} eV and ε c ​ u ​ t = 5 × 10 16 subscript 𝜀 𝑐 𝑢 𝑡 5 superscript 10 16 \varepsilon_{cut}=5\times 10^{16} eV.

Model	$T_p/T_e$	$W_{CR}$ [erg s-1]	$\kappa$	$W_{CR}/W_{rec}$	$W_{CR}/⟨\dot{M_{acc}}⟩c^2$	$W_{esc}$ [erg s-1]
m1	125	$2\times {10}^{43}$	1.3	13.2	2.7$\times {10}^{-1}$	1.9$\times {10}^{42}$
m2	190	$4\times {10}^{42}$	1.0	2.4	5.4$\times {10}^{-2}$	1.2$\times {10}^{41}$
m3	240	$9\times {10}^{41}$	1.0	0.3	1.2$\times {10}^{-2}$	2.6$\times {10}^{40}$

Equations

$$W_{rec}=1.52\times 10^{42}f\left(\frac{{\dot{M}}_{acc}}{\text{mbox}{M}_{\odot }\text{mbox}{yr}^{-1}}\right)\left(\frac{T_p}{T_e}\right)\text{mbox}{ergs}^{-1},$$

Full text

1 Introduction

Misaligned radio-galaxies emitting at very-high-energies (VHE¿100 GeV; Centaurus A [1], M87 [2], IC 310 [3], NGC 1275 [4], 3C 264 [5], PKS 0625-354 [6]111The nature of PKS 0625-354 and IC 310 present features of both radio-galaxies and blazars; see also [7]. ) offer a unique perspective for probing particle acceleration and emission processes in black hole (BH) accretion-jet systems. The misaligned nature of these sources challenge the interpretation of the VHE data with a single zone synchrotron self-Compton (SSC) jet-emission scenario, as Doppler boosting factors of only a few are allowed to model the observed spectrum. Multi-zone SSC scenarios can explain the SED including the VHE part, but more free parameters have to be introduced. (see e.g., [8], [1]). Larger Doppler factors can be applied assuming a bent jet scenario, i.e., that VHE emission is mainly produced within the zone where the jet is more aligned to the observer line of sight (see e.g., [9]). Alternatively, the VHE emission component of radio galaxies could be of is of hadronic origin222Hadronic processes in astrophysics have recently been favoured by simultaneous neutrino and $\gamma$-ray detections in the direction of the blazar TXS 0506+056 ([10]) and two more multi-messenger associations of sources not yet identified [11]. and in the core region ([12], [13], [14], [15], [16]).

Here we consider the scenario where the VHE emission of radio-galaxies is produced by the interactions of CRs with the BH accretion flow, where CRs are accelerated by turbulent magnetic reconnection. We then investigate the conditions of the accretion flow as well as the distribution of injected CRs required to reproduce the observed data. Magnetic reconnection has previously been discussed as an efficient CR acceleration mechanism in turbulent and magnetised plasmas ([17], [18], [19], [20], [21], [22], [23]). A description of the phenomenology and analytical and numerical studies of turbulent magnetic reconnection in BH accretion flows and jets can be found in the contributions by [24] and [25] in these Proceedings.

In this work, we combine three numerical techniques to model the accretion flow environment and propagation of CRs to self-consistently investigate the emission and absorption of VHE $\gamma$-rays as well as the production of neutrinos within the accretion flow zone. We then compare the CR power required to reproduce the VHE SED with the magnetic reconnection power of the accretion flow considering the analytical model of [21] for magnetic reconnection in magnetically dominated accretion flows (MDAF). Here we focus our analysis on Centaurus A (Cen A), the closest radio-galaxy emitting in the VHE regime. However, our approach is aimed to be applied in different sources displaying VHE emission and radiative inefficient accretion flows (RIAFs; [26], [27]).

In the next section we describe the numerical accretion flow simulation that we employ to obtain the gas density, magnetic, and photon fields where CRs propagate and interact. In Section 3 we describe the Monte Carlo simulation of CRs plus electromagnetic cascading that we employ to account for gamma-ray and and neutrino emission. Finally, we discuss our results in Section 4.

2 The numerical model for background accretion flow

We adopt a numerical general relativistic (GR) magneto-hydrodynamic (MHD) RIAF approach together with GR synchrotron + inverse Compton radiative transfer to model the gas density, magnetic and photon fields in the BH accretion flow of Cen A. To do this, we employ the axi-symmetric harm code [28], together with the radiative transfer grmonty code [29]. The accretion flow is simulated within a spatial domain of 40 $R_{g}$, where the accretion is triggered by magneto-rotational-instability (MRI) on a torus in initial equilibrium having a poloidal magnetic field with a maximum plasma beta $\beta=25$. We use the dimensionless BH spin parameter $a=0.94$, the gas specific heat ratio $\gamma=4/3$ and a 256$\times$256 resolution for the $R$ and $\theta$ spherical coordinates.

Different initial magnetic field configurations are expected to produce different dynamical effects in the accretion flow. For instance, a multi-loop poloidal initial magnetic field can produce, after enough evolution time, suppression of MRI and very efficient jets, forming the so-called magnetically arrested disc [30], [31]). By contranst, an accretion flow solution with turbulent magnetic field can be obtained with an initial magnetic field porfile with a single poloidal loop which produce a less efficient jet [32],[33]. For the case of an initial toroidal magnetic field, a spontaneous dipolar magnetic flux emerge producing transient relativistic jets [31]. These different flow behaviours are also expected to produce different outcomes in the soft photon field emission as well as in the CRs emission. We will explore these possibilities in a future work.

For the simulation employed here, we adopt the mass of $M_{BH}=5\times 10^{7}$ M⊙ for the central BH in Cen A to define the spatial and temporal scales333 The spatial and temporal scales are defined as $R_{g}=GM_{BH}/c^{2}$ and $R_{g}/c$, respectively, where $G$ is the gravitational constant, $c$ the speed of light and $M_{BH}$ is the mass of the super massive BH of Cen A.. Thus, we define the gas density scale $\rho_{0}$ fixing the mass accretion rate with the value of $\dot{M}_{acc}=1.3\times 10^{-3}$ M⊙ yr${}^{-1}\simeq 1.17\times 10^{-3}$ MEdd (consistent with the accretion rate obtained by [34] for Cen A). Given the gas density scale, the three-magnetic field scale is defined as $B_{0}=c\sqrt{4\pi\rho_{0}}$.

For a given snapshot, the photon field is obtained calculating the radiation flux with the grmonty code at different radius and polar angles within the accretion flow zone. Here we assume a constant proton-to-electron temperature ratio $T_{p}/T_{e}$ along all the spatial domain444 A more appropriate model for the electron temperature should possibly consider its dependence on the plasma magnetisation, see e.g., [35], [36], [33]., where electrons follow a relativistic thermal energy distribution. In Fig. 1 we show the gas number density and magnetic field maps of the simulation snapshot at the integration time of $t=3600R_{g}/c$. We note that in a time interval $\Delta t\simeq 400R_{g}/c$ centered in the snapshot showed in Fig. 1, the accretion flow maintains a turbulent magnetic field, with an overall behaviour with no drastic changes. Thus, we choose a sequence of snapshots within this time interval, to mimic a quasi-continuous injection and emission of CRs (see the next section).

In Fig. 2 we plot the associated photon field map obtained with $T_{p}/T_{e}=190$. We use the $T_{p}/T_{e}$ ratio as a free parameter to obtain different accretion flow models. These models have the same gas density and magnetic field, but different photon field profiles that correspond to the values of $T_{p}/T_{e}=125$, 190, and 240 (see also Table 1)555 The temperature of thermal protons is obtained assuming a $\gamma$-law equation of state $p=(\gamma-1)u$, where p is the gas pressure and $u$ the gas internal energy and thus, $T_{p}$ $\propto u/\rho$. These variables are given by the GRMHD numerical simulation output and therefore are function of the spherical radius $R$ and the polar angle $\theta$. For the proton-to-electrom temperature values assumed in this work, the electron temperature $T_{e}$ take values in the range of $10^{8-10}$ K, within 20 $R_{g}$. The spectral energy distribution of these photon field profiles are represented in Fig. 4 with the histograms in the energy region $<10^{7.5}$ eV.

3 Monte Carlo simulations of CR emission and $\gamma$-$\gamma$+IC cascading

If magnetic reconnection efficiently accelerates CRs in the accretion flow of Cen A, they could potentially produce gamma-rays and neutrinos due to hadronic interactions with the accretion flow environment. Once they are produced, neutrinos escape practically with no absorption. On the other hand, $\gamma$-rays are susceptible to be absorbed and develop electromagnetic cascades due to the soft radiation field of the accretion flow environment. Here we assess if the CR power required to reproduce the VHE data of Cen A, can be produced by the magnetic reconnection power in the accretion flow.

We consider the magnetic field, gas density and photon field given by the GRMHD model described in the previous section, as the target environment for hadronic interactions. We then employ the Monte Carlo CRPropa 3 code [37] to simulate: (i) proton-proton interactions, (ii) photomeson interactions, (iii) $\gamma$-$\gamma$ pair production, (iv) inverse Compton scattering of secondary leptons produced in the previous interactions, and (v) synchrotron cooling to account for energy losses of charged particles. However, to save computational resources we do not follow the photons produced by this last process666 We note that this synchrotron photons have energies ¡ 100 GeV, and in this work we are mainly interested in modelling the VHE region of the SED. Tracking the synchrotron photons to low energies is computationally expensive for CRPropa simulations, where the trajectory of every particle generated within the accretion flow zone is followed..

The mean-free-paths (MFPs) of particles and photons used in the CRPropa 3 simulations are calculated with the gas density and the photon field obtained from the GRMHD simulation. This accretion flow simulation is axi-symmetric and thus, the calculated MFPs of particles and photons vary with the radial and polar positions along their propagation. We show in Fig. 3 the interaction MFPs corresponding to the accretion flow snapshot of Fig. 1 and 2 at different radii from the central BH. In this background environment, CRs with energies $\sim 10^{15-17}$ eV would produce gamma-rays and neutrinos by photon-pion processes and proton-proton interactions. On the other hand, CRs with energies $\lesssim 10^{14}$ eV would produce $\gamma$-rays and neutrinos mainly by proton-proton interactions. The strongest absorption by pair creation occurs for gamma-rays with energies of $\sim 10^{14}$ eV. IC scattering of secondary leptons contribute with $\gamma$-ray emission principally for energies ¡$10^{13}$ eV.

If the CR trajectories were ballistic within the accretion flow, proton-proton as well as photo-pion interactions were highly inefficient as the MFPs for this interactions exceed the spatial domain of the considered model (40 $R_{g}$) for most of the energy range (see Fig. 3). We note however that CRs are readily trapped by the magnetic field of background environment (see Fig. 1), which favours the efficiency of CR interactions777 Consider for instance the Larmor radius of a PeV proton in the presence of a magnetic field of $10^{2}$ G: $R_{L}\simeq 3.3\times 10^{6}(E_{p}$/GeV)(G/$B)$ cm $\simeq 10^{10}$ cm $\simeq 10^{-2}R_{g}.$.

To reproduce the VHE data of Cen A, we consider three accretion flow models, corresponding to the same mass accretion rate but with different photon field profiles, which are obtained using $T_{p}/T_{e}=125$, 190, and 240. For each model, we emulate a continuous injection of CRs by simulating four burst-like injection of CRs, each one within the GRMHD snapshots at $t_{i}=$3440, 3520, 3600, and 3680 $ct/R_{g}$ (similarly to the CR simulations in [27]). We then calculate the observed fluxes of gamma-rays as $\nu F_{\nu}=(4\pi R_{s}^{2})^{-1}\epsilon^{2}\dot{N}_{\epsilon}/\Delta\epsilon\,\,(\mbox{erg}\,\mbox{s}^{-1}\,\mbox{cm}^{-2})$, where $R_{s}=3.8$ Mpc is the distance from CenA to us, $\Delta\epsilon$ is the size of the energy bin, and $\dot{N}_{\epsilon}$ is the number of photons per unit time within the energy bin $\epsilon$ that leave the spherical detection boundary at $40R_{g}$. The rate of escaping photons is calculated as $\dot{N}_{\epsilon}=\frac{1}{\Delta t}\sum\limits_{i=1}^{4}N_{\epsilon,i},$ where $N_{\epsilon,i}$ are the number of photons within the energy bin $\epsilon$, arriving at the detection boundary (at $40R_{g}$), and produced by the burst-like simulation of CRs injected within the snapshot at the time $t_{i}$. $\Delta t=320R_{g}$ is the time interval during which the code register the arrivals of particles/photons at the spherical boundary, produced by all the burst-like CR injection events.

In all background models, we only consider CR protons for simplicity, injecting them with a power-law energy distribution with exponential cut-off $dN/d\epsilon\propto\epsilon^{-\kappa}\exp\{-\epsilon/\epsilon_{cut}\}$. According to the analytical calculations of [15], the acceleration rate by magnetic reconnection is balanced by CR energy losses at $\sim 3\times 10^{17}$ eV for parameters corresponding to Cen A (see their Fig.3). On the other hand, [38] show that a similar acceleration mechanism driven by reconnection can accelerate CRs up to 10 PeV in typical active galactic nuclei hosting hot accretion flows. We then choose the intermediate value of $\epsilon_{cut}=5\times 10^{16}$ eV for the cut-off energy in our models. For the minimum energy of CR injection we choose the value of $\epsilon_{min}=10^{13}$ eV, which gives the best results to reproduce, or partially reproduce the VHE data.

The calculated SED of these models is shown in Fig. 4 and the associated neutrino emission in Fig. 5. The power-law index $\kappa$, the power of injected CRs $W_{CR}$, and the power of CRs that escape the spherical detection boundary $W_{esc}$ (at 40$R_{g}$), are listed Table 1 for each emission model. In this table we also compare the CR power $W_{CR}$ with the magnetic reconnection power $W_{rec}$ and with the accretion power $\langle\dot{M}_{acc}\rangle c^{2}$. The magnetic reconnection power is calculated as a function of the accretion rate and the proton-to-electron temperature of the accretion plasma, following the analytical model of [21]:

[TABLE]

where $f$ is a combination of dimensionless factors888We define $f\equiv A\Gamma^{-1}(11.5\alpha^{10/3}+14.89)^{1/2}$, being $A$ the ratio of the height to the radius of the magnetic reconnection zone, $\Gamma^{-1}=[1+(v_{A}/c)^{2}]^{1/2}$ the relativistic correction factor of the Alfven velocity, and $\alpha$ the viscosity parameter. See [21] for details. and here we adopt the value $f=5$.

4 Summary and Discussion

We investigate the scenario where the VHE of radio-galaxies is produced by turbulent magnetic reconnection in the BH accretion flow. To do this, we simulate the interaction of CRs with the BH accretion flow combining (i) numerical GRMHD, (ii) leptonic radiative transfer, and (iii) Monte Carlo propagation of CR plus electromagnetic cascades within the accretion flow zone which, based on analytical and numerical results, subject to fast reconnection driven by turbulence (see [25] and [24] for details). We then compare the power of CRs required to reproduce the VHE data, with the magnetic reconnection power given by the analytical model of [21].

Here we apply our approach to model the VHE emission of Cen A and consider three accretion flow models corresponding to a fixed value of the accretion rate and different values of the $T_{p}/T_{e}$ ratio giving different photon field profiles (see Table 1 and Fig. 4). Less power of CR injection is required to match data points in the VHE regime for background models with lower photon field densities (see Table 1.). This is naturally expected as weaker radiation fields imply less absorption of VHE $\gamma$-rays. From the models presented in this work, only the model m3 is consistent with the magnetic reconnection power (i.e., satisfying the condition $W_{CR}/W_{rec}<1$). This emission model matches only two points of the VHE data and partially contributes to the rest of the VHE tail. The leptonic soft radiation of this model m3 do not contribute substantially to the radio-GeV emission, which is in agreement with the picture where the radio-GeV data of Cen A is produced in the sub-parsec jet by SSC of a single electron population ([44], [45], [46]). Furthermore, the neutrino flux associated to the model m3 is well below the upper limits predicted by the IceCube extragalactic diffuse emission, which is more realistic for a single source in comparison with the neutrino fluxes of the models m2 and m3. This is simply because the contribution of photo-pion interactions to neutrino production is diminished as the radiation field is lower. We then conclude that CRs accelerated by magnetic reconnection are able to partially contribute to the current VHE data of Cen A with $\gamma$-rays produced within the accretion flow, provided that this accreting plasma makes no substantial contribution to the radio-GeV observed emission.

As described in the previous section, the neutrino and VHE radiation fluxes are calculated by simulating hadronic interactions and electromagnetic cascading inside a spherical boundary at 40 $R_{g}$ from the BH. Beyond this detection boundary no further emission and absorption processes are considered. However, further interactions of the CRs that escape this detection boundary might contribute to the observed gamma-ray and neutrinos fluxes. Particularly, interactions with the circumnuclear disk (CND) that surrounds the core of Cen A ([47], [48], [49]), might appear as core emission if these interactions take place along our line of sight. The CND is observed with a projected size of $\sim$200$\times$400 pc [34], and has regions with densities of $\sim 10^{3-5}$ cm*-3* ([50]). The mean-free-path of CRs travelling inside the CND could then be of $\sim 10^{3}$ pc and thus, a non-negligible fraction the CRs escaping the accretion zone might produce $\gamma$-rays within this gaseous structure. This propagation effect will be investigated in a forthcoming paper to complete the study on the VHE emission of Cen A initiated in this proceeding.

5 Acknowledgments

We acknowledge support from the Brazilian agencies FAPESP (grant 2013/10559-5) and CNPq (grant 308643/2017-8). The simulations presented in this lecture have made use of the computing facilities of the GAPAE group (IAG-USP) and the Laboratory of Astroinformatics IAG/USP, NAT/Unicsul (FAPESP grant 2009/54006-4). RAB is supported by the FAPESP grant 2017/12828-4 and JCRR by the FAPESP grant 2017/12188-5.