Recent advances in numerical simulations of magnetically arrested accretion onto supermassive black holes have shed light on the formation and dynamics of magnetospheric current sheets near the black hole horizon. By considering the pair magnetization
An accreting supermassive black hole (SMBH) is the powerhouse of an active galactic nucleus (AGN). M87 is a massive elliptical radio galaxy at a distance of 16.8 Mpc and has a black hole mass about 6.5 • 10 9 M ⊙ [1]. The non-thermal electromagnetic emission emerging from the vicinity of its SMBH, M87*, is variable, occasionally exhibiting flares in X-rays and very high-energy (VHE, > 0.35 TeV) γ-rays [2][3][4], on timescales as short as a few times the light crossing time (1 -20t cr ) of the black hole's gravitational radius. The production of non-thermal radiation from the vicinity of the SMBH conveys the potential for particle acceleration to high energies, thereby suggesting the plausibility of proton acceleration to ultra-high energies. Notably, M87* has been postulated as a potential accelerator of Ultra High Energy Cosmic Rays (UHECRs) [5,6]. The properties of the accelerated particle population, such as their distribution in energy and maximal energy reached, are determined by the dominant particle acceleration mechanism at work and the physical conditions in the accelerator region. Recently, the Event Horizon Telescope (EHT) has provided insights into the innermost regions of the accretion flow surrounding M87*. By analyzing these images alongside predictions from general relativistic magnetohydrodynamic (GRMHD) models, the
Collaboration suggests a scenario in which a magnetically arrested disk (MAD) exists around a black hole rotating with a moderate to high spin parameter [7].
High-resolution three-dimensional (3D) GRMHD simulations of accretion onto rotating SMBHs have revealed transient and non-axisymmetric magnetospheres during periods of magnetic flux decay at the horizon, leading to significant drops in mass accretion rates and the emergence of thin equatorial current sheets separating the two sides of the polar jet [8]. The intricate dynamics of these current sheets, including plasmoid-mediated reconnection and the injection of reconnection-heated plasma into the accretion disk and jet boundary, can provide a dynamic environment where particles can accelerate and eventually power outbursts of electromagnetic radiation, commonly known as flares. However, challenges persist in understanding the triggering mechanisms behind large flux eruption events and how they are impacted by the system's physical parameters, such as the spin of the SMBH, or when the system is not in the MAD regime, and others.
Magnetic reconnection in plasmas with high magnetization (σ ≫ 1), defined as twice the ratio of the magnetic energy density to the plasma enthalpy density, is an efficient process of magnetic energy dissipation. A significant fraction of the dissipated energy is used to accelerate particles to relativistic energies ∼ σmc 2 [9][10][11]. The physics of reconnection can only be captured from first principles utilizing fully kinetic particle-in-cell (PIC) simulations. A large body of work on 2D simulations of relativistic reconnection in pair plasmas has demonstrated that particles are accelerated into power-law distributions, dN/dγ ∝ γ -p , that extend up to ∼ σ and with a power-law slope p that depends on the plasma magnetization [11,12]. Particles reaching the highest energies were injected into the acceleration process by interactions with the non-ideal electric fields at the so-called X-points of the current sheet, regions where the conditions of ideal MHD break down [13]. A secondary acceleration that could push particles to γ ≫ σ was also found to operate for particles trapped in compressing plasmoids, albeit on much longer timescales [14,15]. Recent large-box 3D PIC simulations of reconnection in pair plasmas have revealed that particles with γ > σ can escape from plasmoids (flux ropes in 3D) and can accelerate to even higher energies while meandering between the two sides of the reconnection layer in the upstream (inflow) region [16,17]. These particles will terminate their acceleration phase by reentering into the flux ropes where the acceleration is not efficient or escape from the system. The primary goal of this paper is to examine the role of magnetospheric current sheets in the context of magnetically arrested accretion onto SMBHs, with a particular focus on M87*. We explore the evolution of relativistic particle populations accelerated in current sheets during magnetic reconnection, and investigate their implications for the production of non-thermal radiation and the potential acceleration of ultra-high energy protons. To achieve this goal, we employ a two-step approach that combines theoretical modeling and numerical calculations. Motivated by the findings of recent 3D PIC simulations, we develop a model for non-thermal emission from the reconnection region, accounting for the presence of two distinct particle populations ("free" and "trapped", located at the reconnection upstream and flux ropes, respectively) [17], as well as the interplay between these two populations and secondary pairs injected into the system through γγ pair creation. The main parameters of the model are the pair magnetization in the upstream region of the current sheet, the JCAP12(2024)009 length of the current sheet, and the mass accretion rate. Subsequently, we utilize numerical methods to calculate the non-thermal radiation spectra arising from these pair populations within a one-zone leptonic model, considering interactions with both non-thermal and disk photon fields. By applying our model to the specific case of M87*, we aim to provide insights into the observed phenomena and discuss the implications of our findings on magnetospheric current sheets and their potential role in particle acceleration near SMBHs.
This paper is structured as follows. In section 2 we outline our theoretical framework and present the key parameters of our model. In section 2.1 we describe the pairs that are in the free acceleration phase while in section 2.2 we discuss the pair dynamics in the downstream region analytically. In section 2.3 we showcase the numerical approach to the problem. In section 3 we apply our model to M87*. In section 3.1 we present the numerical results of the pair and photon spectrum. Using these numerical results in section 3.2 we determine which combinations of the key parameters can reproduce some of the observations in the X-rays. In section 3.3 we check whether the production of ultra-high energy protons in M87* is feasible. In section 3.4 we discuss what is the effect of the pair creation in the system. Finally, we present the conclusions of this work in section 4.
Recent 3D high-resolution general relativistic MHD simulations of accretion onto a black hole through a magnetically arrested disk [MAD, 18,19] have revealed the formation of short-lived current sheets in the black hole magnetospheric region [8]. The reconnection layers have a typical length of a few gravitational radii of the supermassive black hole (SMBH), i.e. l ≲ 10r g , where r g = GM/c 2 ≃ 1.6 • 10 14 M/(10 9 M ⊙ ) cm, and their lifespan is of the order of 10 l/c for the largest current sheets [8].
One can estimate the strength of the magnetic field threading the black hole horizon by noting that in the MAD regime the dimensionless magnetic flux threading the black hole, ϕ BH = Φ BH / Ṁ cr 2 g , cannot exceed an approximate value of 50 [19]. This translates to the relation
where Ṁ is the accretion rate onto the SMBH, and Φ BH = 4πr 2 H B 0 is the magnetic flux threading the black hole horizon of radius r H = r g f (a s ), where f (a s ) = (1 + 1 -a 2 s ) and a s is the dimensionless spin of the SMBH. Introducing the dimensionless accretion rate, ṁ = Ṁ / ṀEdd , where ṀEdd = L Edd /(η c c 2 ) is the Eddington accretion rate and η c is a matter-to-luminosity conversion factor [20], we can estimate the magnetic field strength as
where we have introduced the notation Q x = Q/10 x and M is expressed in units of the solar mass. In what follows, we assume that B 0 is the typical magnetic field strength in the upstream region of the magnetospheric current sheets. Recent 3D kinetic simulations of reconnection in pair plasmas [17,21] have shown that the total pair distribution in the reconnection region is composed of the so-called "free"
particles that undergo active acceleration in the inflow region, by meandering between the two sides of the reconnection layer, and the "trapped" particles that end up inside plasmoids (flux tubes) and do not undergo further acceleration. These can cool due to radiative losses and eventually get advected out from the reconnection region. We develop an analytical model to describe the steady-state populations of free and trapped particles and then complement our analysis with numerical calculations of non-thermal radiation accounting for physical processes, like inverse Compton scattering and photon-photon pair production.
We assume that cold pairs dominate the upstream plasma. The pair magnetization σ e , which is one of the main free parameters of the model, is then defined as
where n e ± is the pair density in the upstream. The injection rate of pairs to the downstream region can be written as,
where η rec = v rec /c ∼ 0.06 is the quasi-steady state reconnection rate achieved in 3D PIC simulations of reconnection [17,21], A = πR 2 is the area of a cylindrical current sheet of radius R ≡ Rr g , and the factor of 2 accounts for the inflow of particles from both sides of the sheet. We assume that a fraction ζ of the pairs entering the layer will become "free" after being first accelerated to a Lorentz factor γ inj ∼ σ e ≫ 1 within the reconnection layer.
The free pairs may leave the upstream region of active acceleration and become eventually trapped in the plasmoids of the current sheet, where they continue cooling due to radiative losses. The remaining (1 -ζ) fraction of the pairs entering the layer will only experience the injection stage within the current sheet (e.g. at X-points [13,22]), without ever becoming free in the upstream, before getting trapped in plasmoids. We adopt ζ = 0.06 as a typical value that is motivated by simulations [17] (for more details, see appendix A).
In what follows we compute the distribution functions of the free and trapped pairs, and the resulting synchrotron spectra using a semi-analytical approach.
We first consider the emission produced by pairs in the free acceleration phase. In the free phase pairs can be accelerated by the ideal electric field and lose energy due to synchrotron radiation 1 . The kinetic equation governing their differential number distribution, N free (γ) ≡ dN free /dγ, is written as [23]:
1 In the analytical model we ignore inverse Compton losses (ICS) due to disk radiation, but these are included in the numerical calculations. Nevertheless, the magnetic energy density is typically much larger than the disk photon radiation density, making ICS losses negligible. For instance, using parameters relevant to M87*, we find
where L soft ≈ ν EHT Lν EHT ≃ 10 41 erg s -1 , ν EHT = 230 GHz, R soft = 5rg and, B0 = 100G.
-4 -
where γ inj ∼ σ e . Eq. (2.6) is a simplified version of a Boltzmann-type equation for the particle distribution function defined in a 6-dimensional phase space F (⃗ x, ⃗ p, t). In this simplified form, we do not account for the orientation of momenta or the spatial location of the particles. Effectively, it can be seen as the equation obtained after averaging over these variables. Using eqs. (2.2) and (2.3) the free particle injection rate can be expressed in terms of the main model parameters as
The corresponding energy injection rate can be estimated as
where L BZ ≃ πr 2 g B 2 0 c/24 is the power extracted from a maximally rotating black hole through the Blandford -Znajek process [19,24].
The acceleration and synchrotron loss rates in eq. (2.6) are defined respectively as γacc = η rec eB 0 κ z m e c ≡ β a , (2.9)
where κ z ∼ 1 is the particle velocity along the z direction of the electric current in units of c (in what follows we set κ z = 1 and do not show it explicitly). We note that a pitch angle of π/2 was assumed in eq. (2.10) because free particles move preferentially along the z direction and therefore perpendicular to the magnetic field direction. Finally, t fr esc (γ) is the escape timescale from the free acceleration phase,
(2.11)
as found in PIC simulations [17]. In figure 1 we present the relevant timescales for an electron or a positron while also including the advection timescale defined in section 2.2.1 where we assume R = 1. The solution of eq. (2.6) is [17,23]:
where s free = t acc /t fr esc ≈ 1, and the characteristic timescale τ is given by
In the above expressions γ rad = β a /β s is the synchrotron radiation-limited Lorentz factor of pairs, which is found by solving γacc = -γsyn , and it can also be expressed as
JCAP12(2024)009 The quantity τ (γ) (see eq. (2.13)) represents the characteristic timescale over which the system reaches a steady state for a given Lorentz factor γ. τ (γ) is an increasing function of γ, indicating that pairs with higher Lorentz factors take longer to reach a steady state due to significant energy losses and larger acceleration timescale. The condition for which eq. (2.12) is valid can be rewritten in the following form γ ≤ γ rad tanh 2βαt γ rad + ln γ rad +γ inj γ rad -γ inj . This expression defines the Lorentz factor γ at time t where there is a balance between acceleration and energy losses.
From this point on we restrict our analysis to the regime of σ e < γ rad ∼ 10 6 because the free particle channel would not exist otherwise [16,17]. We will return to this point later in the Discussion.
The pair distribution of eq. (2.12) can be approximated by a power law for γ ≪ γ rad ,
(2.15)
The power-law approximation fails to capture the exact solution given by eq. (2.12) only at γ ≃ γ rad , underestimating it by a factor of 2 -3. Using the approximation of eq. (2.15) we can compute the total energy of the free population, which reads
where we used eqs. (2.7), (2.9), (2.14), (2.15), and we assumed that γ inj ≪ γ rad while deriving the expression in the second line of the equation.
Meanwhile, the total available energy in the system (i.e., the dissipated magnetic energy in the current sheet's lifetime) can be estimated as
where where the factor of 2 accounts for the energy inflow from both sides of the current sheet, and T = 15(R/c) ≃ 0.9 RM 9 days is the assumed lifetime of the current sheet motivated by the results from 3D GRMHD simulations [8], and E rec = η rec B 0 is the ideal electric field in the upstream region of the current sheet. The energy ratio is then given by
where we used eq. (2.2). Only a small fraction of the available system energy is carried by the free particles, as their maximal energy is strongly limited by the synchrotron losses and they cannot tap all the available energy. This can also be verified by comparing γ rad in eq. (2.14) with the maximum Lorentz factor attained within the acceleration region of length l = 2R and electric field E rec ,
Using eqs. (2.14) and (2.19), the energy ratio of eq. ( 2.18) can be simply expressed as
(2.20)
Because of the hard power-law distribution of free pairs (s free = 1) their synchrotron spectrum will peak roughly at the photon energy corresponding to the burnoff limit
where B cr = m 2 e c 3 /(eℏ) ≃ 4.4 • 10 13 G is the Schwinger magnetic field strength. The corresponding peak synchrotron luminosity at the target photon energy can be estimated as
where we used eqs. (2.7) and (2.8). The synchrotron luminosity of the free pairs and the peak synchrotron photon energy are independent of the uncertain pair magnetization in the magnetospheric region. The reason for this outcome is that the differential number density of free pairs dN free /dγ is a power-law with an index of -1 which correspond to an equal number of pairs in each Lorentz factor. Since the magnetic field is independent of σ e the luminosity emitted by this channel will also be independent of the pair magnetization. Another way of thinking about this is that the luminosity that goes into the free pairs L free e,inj is a fraction of the Poynting luminosity which is also independent of σ e .
The peak synchrotron photons may undergo pair production (close to the energy threshold for photon annihilation) when they interact with less energetic synchrotron photons of a energy
The corresponding synchrotron luminosity is
where we have assumed that the synchrotron spectrum is produced by a monoenergetic pair distribution at γ = γ rad . The optical depth for the attenuation of the peak synchrotron photons can be approximated as
where n l is the number density of target photons of energy ϵ l inside the volume of the assumed cylindrical current sheet of height 2η rec R, and σ T /4 is the approximate peak value for the cross section [25]. The energy injected per unit time into secondary pairs from γγ production can be then written as
which is orders of magnitude lower than the free pair injection luminosity, is independent of σ e , and scales quadratically with ṁ. Secondary pairs will be injected into the system with
and their steady state distribution can be obtained from eq. (2.6) after using the appropriate source term (Q sec inj ≈ 2τ γγ L pk syn /ϵ syn ), and neglecting the acceleration term. Therefore, γγ pair production is not expected to alter the total pair density in the system or impact the overall photon spectrum, unless more low-energy photons are available, in addition to the ones provided by the free particles (as we assumed so far). In the next section, we will investigate whether the synchrotron radiation of trapped pairs may provide these additional target photons.
Trapped pairs can be classified into two distinct groups. The first group consists of free pairs that escape the (upstream) region of active acceleration and become trapped in the plasmoids of the reconnection region. These pairs are injected into the trapped phase with Lorentz factor from γ inj ∼ σ e up to γ rad , but with a softer power law than the free pairs, as we describe below. The second group of trapped pairs consists of particles that were accelerated in the current sheet up to σ e (at e.g. X-points [13]) and then became trapped in plasmoids where no further fast acceleration takes place. Consequently, we can express the total injection rate to the trapped population (integrated over γ) as follows,
-8 -
where Q trap X corresponds to the injection rate of trapped particles that went through the injection phase of acceleration but did not experience a free phase of acceleration, and Q trap fr indicates the injection rate of pairs from the free to the trapped phase.
The differential injection rate of pairs from the free to the trapped phase can be expressed as
where we used eq. (2.15) to derive the right-hand side of the expression. These pairs that are trapped in the reconnection region may escape on the advection timescale of plasmoids, i.e.
is the Alfvén speed (see also appendix B of ref. [21]). Meanwhile, they experience synchrotron losses in a magnetic field with comparable strength to the upstream region [26].
The steady-state trapped pair distribution originating from the free phase can be obtained by solving eq. (2.6) without the acceleration term, after using eq. (2.29) as the injection term,foot_0 and after replacing t fr esc (γ) with t adv (see also ref. [17]),
(2.30) where γ syn cool is the Lorentz factor of pairs losing half of their initial energy due to synchrotron cooling within the advection timescale, and it is given by,
The first branch of eq. (2.30) describes the spectrum in the energy range where the injection from the free channel takes place, i.e. from γ inj to γ rad . The second branch is practically zero, unless γ syn cool < γ rad and is more relevant when γ syn cool ≪ γ inj ; this case is also known as the fast-cooling regime. The condition for fast cooling, γ syn cool ≪ γ inj , translates to 2 ≪ σ e η c,-1 ( ṁ-5 R 0 ) -1 . If β s t adv γ= t adv /t syn (γ) ≃ 2γ ṁ-5 R 0 /(f 4 (α s )η c,-1 ) ≫ 1 (where we used eq. (2.10)), the steady-state distribution can be approximated by a broken power law,
The total number of trapped particles originating from the free phase is given by
where the expression in the second line has been derived in the limit of 1 ≪ γ inj ≪ γ rad . The total number of trapped particles originating from the free phase scales as γ -1 inj ∝ σ -1 e , since it is determined by the number of free pairs whose injection rate scales also as σ -1 e (see eq. (2.7)). In this limit, the total energy of the trapped pair population reads,
and is larger than the energy of the free particle population (see eq. (2.16)) for lower magnetization values (higher number densities).
The rate at which pairs are injected directly into the e.g., X-point phase of acceleration is parameterized as (1-ζ)Q tot e,inj . These pairs upon acceleration obtain a power-law spectrum that extends from γ ∼ 1 to γ ∼ σ e with a slope of ∼ -1 in the limit of high magnetizations [13,27,28]. For simplicity, we call these pairs N trap X although their spectrum is not only determined by acceleration at the X-points. The differential injection rate then reads
where the normalization factor 1/ ln(σ e ) is found by requiring that σe
(
The steady-state distribution of X-trapped particles can be computed as [29]
(2.37)
The total energy of this pair channel will be comparable to E trap fr,tot since the injection rates for both pair channels are a fraction of Q tot e,inj and they share the same escape timescale, t adv , from the system.
For the numerical calculation of the non-thermal radiation from the free and trapped pairs in the reconnection region we utilize the code LeHaMoCfoot_1 [30]: a versatile time-dependent lepto-hadronic modeling code that includes synchrotron emission and absorption, inverse Compton scattering (on thermal and non-thermal photons), γγ pair production, photopair (Bethe-Heitler pair production) and photopion production processes, as well as proton-proton inelastic collisions. For this work, we will utilize the leptonic module of the code.
The code has been designed to describe the time evolution of the particle distribution within a spatially homogeneous spherical source. Particles that are injected into the source are assumed to instantaneously fill the whole volume, while all particles escape from the source on the same timescale (i.e., their radial position from the center of the source is not considered JCAP12(2024)009 when computing their geometric escape). Moreover, the emissivities of all processes are angle-averaged assuming isotropic distributions for charged particles and photons.
Therefore we do not account for the different spatial distributions of the free and trapped populations and their emitted photons in our numerical model (all are assumed to occupy the same region since we can not distinguish the upstream and downstream regions). Nonetheless, their spatial differences are encoded in the fact that free particles undergo active acceleration, while trapped particles do not. Moreover, with the adopted numerical framework we cannot directly study potential anisotropic radiative signatures [31]. To ensure some equivalency between the physical current sheet and the numerical emitting region (in regards to densitydependent interaction rates), we consider a sphere of effective radius R eff ≃ 0.5Rη 1/3 rec,-1 that has the same volume as the cylindrical current sheet of radius R and height 2η rec R (motivated by the width of the largest plasmoids in the layer, see figure 2). Finally, to capture the realistic escape of trapped pairs, which happens along the current sheet on the plasma advection timescale, we set their escape timescale equal to R/c. Secondary pairs are mainly produced from the interaction of synchrotron photons from the free population with the other photons in the system. We assume that they are injected within the system's volume, do not undergo the free acceleration phase, and are considered part of the trapped population. Since they are trapped in the midplane, we assume they cool via synchrotron, inverse Compton losses, and escape the system in R/c with R being the half-length of the current sheet.
In this section, we apply our numerical model for radiation from magnetospheric current sheets to the SMBH M87*. We first compute the non-thermal radiation from the free and trapped populations, accounting also for inverse Compton scattering and γγ pair production, for different values of the upstream pair magnetization. We then set limits on the latter parameter using multi-wavelength observations of M87*, and discuss the implications of our findings for the pair enrichment of the upstream plasma and UHE proton acceleration (section 3.1).
Parameter Symbol Value [Units] Distance d M87 16.4 [Mpc] * Black hole mass M 6.5 • 10 9 [M ⊙ ] Dimensionless mass accretion rate ṁ 10 -6 -10 -5** Current sheet half-length R 1 [r g ] Effective radius R eff 0.5 [r g ] Reconnection rate η rec 0.06 Upstream magnetic field B 0 70 -221 [G] † Pair magnetization σ e 10 -10 6 ‡ Fraction of free particle injection ζ 0.06 Multiple of r g R 1 Matter-to-luminosity conversion factor η c 0.1 Spin of the SMBH α s 1 * Adopted from [32]. ** Motivated by modeling of polarimetric observations [7]. † Estimated using eq. (2.2) for the range of ṁ values listed here. ‡ Motivated by the regime of interest σ e < γ rad ∼ 10 6 (see eq. (2.14))
. Such low values imply that the magnetospheric region, where current sheets form, has pair densities much greater than the Goldreich-Julian value. This could materialize if some plasma is channeled from the disk to the layer in the dynamically evolving environment of MAD, or/and if pairs are generated by γγ absorption as described in [33].
Table 1. Model parameters and values used for the application to M87*.
We compute the steady-statefoot_2 distributions of free and trapped pairs, as well as the emitted photon spectra using the code LeHaMoC. All results are obtained at T = 15R eff /c, assuming that this is the relevant lifetime of the current sheet. The parameters are listed in table 1.
In figure 3 we illustrate the decomposition of the pair spectrum under two distinct accretion rate scenarios ( ṁ = 10 -6 on the left and ṁ = 10 -foot_3 on the right), for two different values of the magnetization, σ e = 10 3 (top) and σ e = 10 6 (bottom). The red dashed line in all plots represents the population of free pairs N free (γ) accelerated in the upstream region from γ inj ≈ σ e to γ rad with a power law index ∼ -1 as described in section 2.1. The acceleration process leads to a notable spike in the distribution at γ rad , indicating the accumulation of pairs in this energy range. The pile-up at the radiation-limited Lorentz factor is expected to grow with the lifetime of a current sheet. The small variation in the position of γ rad JCAP12(2024)009 luminosity and the 10 MeV photon luminosity depend linearly on ṁ, see eqs. (2.22) and (2.26).
For lower σ e values we expect the production of more low energy photons compared to cases with higher σ e (see also figure 4). These photons are upscattered to higher energies, i.e., above 10 MeV up to TeV energies (see figure 4). The attenuation of these photons is responsible for the production of secondary pairs above γ > 10 up to γ ∼ 10 6 . We find that the energy that goes into the secondary pairs from the attenuation of TeV photons is more significant for lower σ e values because of the higher energy of target photons.
Figure 4 presents a decomposition of the photon spectrum emitted by different pair populations for the same parameters used for figure 3. We also show with grey markers multi-wavelength observations of M87* [34,35]. The grey solid line in figure 4 represents the photon field from the inner accretion flow of M87* as assumed in our numerical model (which is included as a source of soft photons for inverse Compton scattering). For energies between ∼ 4 • 10 -3 eV and 4 eV where no data are available, we assume that the photon spectrum is a power law whose index is benchmarked by UV observations. 5 The magenta-shaded region is JCAP12(2024)009 constructed using X-ray observations taken over a 12-year-long period with Chandra from the core region of M87 [36]. In general, the X-ray spectrum observed by Chandra is harder when brighter. The red dashed line represents the photon spectrum produced by free pairs N free (γ) via synchrotron and inverse Compton scattering. The synchrotron emission peaks at 10 MeV, originating from pairs with γ ≃ γ rad , and has a peak luminosity that is in agreement with the analytical estimate of eq. ( 2 c,-1 f -2 (a s ), we observe the contribution to the total emission from N trap fr , where the spectrum is flat in νL ν since these photons are produced by cooled pairs with power law index ∼ -3, while at lower energies the spectrum has both N trap fr and N trap X contributions. Synchrotron photons from free pairs are mainly absorbed by the photons originating from N trap fr and serve as the source of secondary pairs N sec γγ , as discussed in the previous paragraph 6 . Because of the dependence of the secondary energy injection rate on ṁ (see eq. (2.26)), their synchrotron emission will be more luminous for higher accretion rates (compare magenta dotted lines in left and right panels of figure 4). The synchrotron spectra of secondary pairs are also limited by synchrotron self-absorption at photon energies below ∼ 10 -3 eV, and their contribution to the total luminosity is negligible compared to the emission from free-trapped pairs above ϵ syn (σ e ). Synchrotron photons can be upscattered by the same pairs to higher energies 7 . Because of the low energy density of the target photons compared to the energy density of the magnetic field the produced luminosity of the inverse Compton scattered photons is typically ∼ 1 -2 orders of magnitude lower than the synchrotron luminosity depending on σ e . The photon spectrum produced in the GeV band by Compton-scattered photons and at GHz frequencies by synchrotron radiation cannot simultaneously explain the experimental data by Fermi-LAT, H.E.S.S. (upper limits), MAGIC, and VERITAS as shown by the grey points in the right panel of figure 5. The origin of this emission may originate from different regions or be attributed to protons.
The left panel of figure 5 depicts the overall distribution of free pairs (solid lines) and trapped pairs (dashed lines) at steady state obtained for different initial values of the pair magnetization (σ e ), while maintaining a constant value of ṁ = 10 -6 . Therefore, the upstream magnetic field is fixed in all cases, and by changing σ e we effectively change the upstream number density of pairs. At this stage, we do not account for the contribution of secondary pairs into the definition of σ e . Dotted lines indicate secondary pairs generated by γγ absorption. The bump of the secondary pair spectrum at γ ∼ 10 remains unaffected because these pairs originate by the attenuation of the 10 MeV photons whose luminosity and energy are independent of σ e as shown in eq. (2.22). Nonetheless, the number of more from the SMBH. 6 In appendix B we discuss the dependence of the γγ pair production optical depth on ṁ and σe. 7 Each pair population "sees" all non-thermal photons (produced by all pair populations) as targets for inverse Compton scattering. We also take into account the disk photons, whose energy density is computed within a sphere of radius R soft > R eff as U soft = 3 L soft /(4πR 2 soft c) ≃ 0.04 erg cm -3 . Here, L soft ≈ ν EHT Lν EHT ≃ 10 41 erg s -1 , ν EHT = 230 GHz, and R soft = 5rg.
-15 -JCAP12(2024)009 energetic secondary pairs (γ ≫ 10) increases with decreasing σ e due to the increase in the inverse Compton scattered photon luminosity (see right panel). On the right panel, we demonstrate the resulting photon spectra, showing the free and trapped pair contributions separately. In summary, higher magnetization values lead to harder X-ray spectra (right panel) and higher number densities of secondary pairs relative to the initial pair density of the system (left panel).
Neglecting at the moment the contribution of secondary pairs in the total number density of plasma (i.e., σ e is defined based on the primary pair density), we can use X-ray observations to constrain σ e for a given ṁ. We remind that the non-thermal emission of the current sheet is transient, lasting about a day, assuming a lifetime of T = 15(R eff /c) ≃ 2.6 days. A comparison to observations obtained on similar timescales is, therefore, the most meaningful.
By considering 4.7 ks X-ray observations from July 31st, 2007 to March 28th, 2019 with Chandra/ACIS [36] (see the magenta shaded area in figure 5) we can determine which combinations of ṁ and σ e result in spectra within the observed range. We adopt the range of accretion rate values provided by ref. [7]. In figure 6 we present the photon index.foot_4 map of the synchrotron pair spectrum in the 0.2-10 keV energy range (see color bar) for various combinations of ṁ and σ e . For σ e ≲ 10 4 we find soft X-ray spectra, as the emission is dominated by the cooled part of the trapped pair population. Instead, for higher magnetizations the spectra become harder, reaching the asymptotic value of 3/2. The non-hatched region indicates solutions that fall inside the observed range of X-ray spectra (see the magenta-shaded region in figure 5). Similarly, we show the 0.2-10 keV X-ray luminosity map as predicted by our model for a combination of σ e and ṁ values. Inspection of both JCAP12(2024)009 panels shows that only a small portion of the considered parameter space is consistent with the spread of X-ray observations from the core region of M87. For σ e ≲ 10 5 , only the lowest accretion rates are viable producing a flat spectrum (Γ ≃ 2), while accretion rates as high as 10 -5 are plausible for σ e ∼ 10 6 .
Ultra-high energy acceleration in the vicinity of SMBHs in AGN has been widely discussed (for a recent review, see [6]). Here, we examine if acceleration by the reconnection-driven electric field in magnetospheric current sheets can push protons to energies above 10 18 eV.
To determine whether acceleration at magnetospheric current sheets in M87* can push protons to ultra-high energies, one must compare the acceleration timescale with the various loss timescales relevant for relativistic protons. Knowledge of the target photon density and spectrum is needed for photomeson and photopair (Bethe-Heitler) production processes.
We consider first the disk photons. Apart from the EHT measurement flux at 230 GHz, whose origin is spatially resolved on a radial scale ≲ 7r g , all other fluxes are likely produced at larger distances in the accretion flow. We consider the next two limiting cases: (i) a high photon-density scenario, where all the observed flux emanates from the same region as the 230 GHz flux, and (ii) a low photon-density scenario, where just the 230 GHz photons from the inner accretion flow act as targets for proton cooling. In addition to the photons from the accretion flow, we take into account the non-thermal photons emitted by the free and trapped pairs in the reconnection region (as described in section 3.1).
The characteristic proton timescales (for their definitions, see appendix C) are plotted in figure 7 against the proton Lorentz factor for typical parameter values (see table 1). The yellow-and red-colored bands indicate respectively the range of synchrotron loss and acceleration timescales obtained by using the lowest (solid lines) and highest (dashed lines) values of the mass accretion rate, as inferred from the modeling of polarimetric observations [7].
factors shown in the figure . For σ e = 10 4 photopion losses on the non-thermal photons from pairs are dominating the proton energy losses, limiting their acceleration to energies E p ≲ 1 EeV (crossing of the dotted blue line and red band). The radiation-limited proton energy, in this case, is lower than the maximum energy that can be achieved in a current sheet of length l = 2R and electric field E rec = 0.1η rec,-1 B 0 ,
Higher magnetizations, e.g. σ e > 10 4 , would lead to lower number densities of target photons (see figure 5), shifting the blue dotted line upwards, and making synchrotron losses the relevant limiting mechanism of proton acceleration.
In figure 8 we show the maximum energy achieved by protons, namely min[E
p,max ] where E (rad) p,max denotes the radiation limited energy of protons due of cooling, across various magnetization values σ e and for two distinct mass accretion rates ṁ. It is evident that E (rad) p,max is contingent upon both parameters, while E (acc) p,max is independent of σ e (see grey solid and dashed-dotted line in figure 8). For magnetization σ e ≲ 10 3 we observe that the maximum energy achieved by protons, which is limited by photopion energy losses, is almost constant because the target photon field does not differ much for these values of magnetization (see e.g. figure 5). The dependence of E (rad) p,max on ṁ for σ e ≲ 10 3 is driven by the fact that the trapped pairs and their synchrotron target photon density are higher for larger accretion rates as higher ṁ values correspond to a higher magnetic field (refer to eq. (2.2)). Therefore the acceleration is limited to lower energy values for larger ṁ. For 10 3 ≲ σ e ≲ 10 5 , proton acceleration is also limited by photopion energy losses. The rising trend of E p,max in the range 10 3 ≲ σ e ≲ 10 5 is attributed to the change of the minimum injection energy of trapped pairs ∼ σ e m e c 2 , which is reflected upon the distribution of target photons produced by these pairs (see the spectral break as σ e increases in figure 5). As σ e increases further, for both values of ṁ, we find a constant E p,max . For the low accretion rate, protons gain the maximum energy available into the system as given by eq. (3.1) (see black dashed-dotted line in figure 8). Conversely, for higher accretion rates and thus stronger magnetic fields, synchrotron losses become the limiting factor in proton acceleration (see blue dashed line in figure 8).
In figure 5 we showed that, for fixed ṁ, the number of secondary pairs injected at γ ≲ 10 becomes progressively larger than the number of trapped pairs for higher initial values of σ e . Therefore, secondary pairs may play a crucial role in the self-regulation of plasma magnetization (see also [37,38]). By employing eqs. (2.22), (2.23), and (2.27) and assuming that the target photon field for the attenuation of 10 MeV photons 9 is constant for all magnetizations we can show that the ratio between the secondary pair density (n sec γγ ) and the number density n e ± corresponding to the initial value of σ e (see eq. (2.3)) is given by
We also numerically explored the non-thermal radiation arising from the free and trapped pair populations in a magnetospheric current sheet for parameters relevant to M87* and different initial magnetization values, accounting for synchrotron emission and absorption, inverse Compton scattering of the disk and synchrotron photons, and photon-photon pair production. For σ e ≲ 10 4 the emission from the trapped pair population overshoots the optical/UV and X-ray flux measurements, unless ṁ ∼ 10 -6 (or even lower). The emission from the free population of particles is, however, independent of σ e and peaks at ∼ 10 MeV energies, as shown in figure 5. A transient magnetospheric current sheet may therefore power a MeV photon flare with a hard spectrum and day-long duration (comparable to the lifetime of a current sheet with a length of a few gravitational radii). The luminosity of the flare depends on the reconnection rate, which is well constrained from PIC simulations, the fraction of particles participating in the free phase of acceleration, which can be benchmarked with PIC simulations, and on source parameters, like the length of the current sheet, the accretion rate and black hole mass and spin (see eq. (2.22)). Note that the MeV flare luminosity does not depend on the unknown plasma magnetization in the magnetospheric region. Therefore, a MeV γ-ray monitoring instrument, like AMEGO-X [39], would be ideal for testing the particle acceleration scenario in magnetospheric current sheets.
Our analysis also revealed that non-thermal radiation from pairs may limit the acceleration of protons in current sheets due to photopion losses. The proton radiation-limited energy depends on σ e and ṁ through the low-energy (< 0.1 eV for E p ≃ 10 18 eV) photon number density when photopion losses are dominant, but also through the magnetic field whenever synchrotron losses dominate. The radiation-limited proton energy generally increases with increasing σ e till it reaches the maximum energy possible by the reconnection electric field along the whole length of the current sheet. The latter can be as high as 20 EeV for the largest magnetospheric current sheets (length of 10 r g ) in M87*. The length of the current sheet does not impact the target photon energy density. Therefore, the photomeson loss timescale will be the same for different sizes of the current sheet. Even though the available potential energy for proton acceleration is larger in larger current sheets, their acceleration is limited by synchrotron losses. By utilizing eq. (C.2) and eq. (3.1) we can show that E
p,max is the synchrotron radiation-limited proton energy. Therefore for larger current sheets and σ e > 10 5 the acceleration will be limited by the synchrotron losses to 3 EeV M
rec,-1 f (a s ) energies. Sgr A* has been discussed as a PeV particle accelerator [5]. We therefore discuss the implications of our model for particle acceleration in magnetospheric current sheets for Sgr A*. In contrast to M87*, Sgr A* possesses a lower mass M Sgr A * ≃ 4 • 10 6 M ⊙ [40], and accretes at an even lower rate, i.e., ṁ ∈ (10 -8 , 10 -7 ) (assuming η c = 0.1) [40]. Assuming that accretion onto the black hole happens in the MAD regime, the magnetic field strength in the magnetospheric region of Sgr A*, as determined by eq. (2.2), is comparable with the one in M87*. Given that the length of the current sheets is a multiple of r g , the electric potential drop along a magnetospheric current sheet in Sgr A* will be lower than in M87*. The upper limit on the energy achieved by protons in Sgr A* is found, using eq. (3.1), to be
This suggests that radiation losses will not restrict proton acceleration, unlike in M87*, enabling protons to tap the system's maximum energy potential. Therefore, magnetospheric current sheets in Sgr A* can push protons to maximum energies of a few PeV.
In the analyses carried out in sections 2 and 3, we incorporated the influence of the supermassive black hole's spin, denoted as a s , across all evaluated physical properties. Assuming that current sheets may also form for slowly rotating black holes and that the MAD relation in eq. (2.1) equally applies, we can assess the effects of spin on the energy spectrum and particle acceleration in magnetospheric current sheets. For example, the magnetic field intensity in the upstream region B 0 , the synchrotron limited Lorentz factor γ rad , and consequently, the synchrotron luminosity of free pairs are all influenced by f (a s ), which is defined as 1 + 1 -a 2 s (see eqs. (2.2), (2.14), and (2.22)). Given that f ranges from 1 to 2 for maximally spinning and non-rotating black holes, respectively, the MeV peak synchrotron luminosity could be lower by a factor of 10 for a black hole with a s ∼ 0.1. Additionally, the highest energy that protons can attain in the acceleration zone scales as f -2 (a s ), with smaller spins corresponding to lower proton energies.
For parameters relevant to M87* we showed that the trapped pair population cools due to synchrotron losses to γ ≪ γ inj ≃ σ e (see eq. (2.31)). Nonetheless, our model is based on the results of 3D PIC simulations without radiative losses. Synchrotron losses of plasma trapped in flux tubes can lead to reduced internal pressure, thus decreasing the transverse size of plasmoid structures in the reconnection region, as demonstrated in 3D radiative PIC simulations [31]. In this case, the acceleration timescale of free particles, or the escape timescales of free particles from the region of active acceleration could differ from those adopted here. Moreover, the ratio ξ between the free and the trapped population might be affected. To robustly address this issue, a detailed investigation of the free phase of acceleration in the regime of strong synchrotron cooling (i.e., γ syn cool ≪ σ e ≪ γ rad ) is needed. The numerical calculations of the non-thermal emission were performed in a single-zone framework (see section 2.3). According to this, all particle populations (free and trapped) occupy the same volume wherein the magnetic field and photon fields are homogeneous. In reality, the current sheet is a much more complex system, with trapped particles residing mainly in flux tubes and free particles moving above and below the current sheet, sampling the fields in the upstream region close to the sheet. Pairs entering the free phase are doing quasiperiodic deflections between the two sides of the reconnection layer with their distance from the midplane being < 0.1l [16]. Therefore in zeroth approximation, we can assume that the system is evolving in the same one-zone blob and the photon fields produced by the distinct pair populations coexist within the same region.
We estimated the maximum energy that protons can achieve in magnetospheric current sheets, assuming the same rate of acceleration as for the free electrons. To make further predictions for the electromagnetic (e.g., proton synchrotron radiation) and neutrino signatures of the relativistic hadronic population, one would have to know the number density of protons in the magnetospheric region and the shape of their distribution. While it is unlikely that the plasma in the magnetospheric region of an accreting black hole consists of electrons and protons, it is possible for a small number of protons (compared to the pairs) to be present. Specifically, in our model, we have operated under the assumption that the upstream JCAP12(2024)009 magnetization is governed by the pair number density; in other words, the rest mass energy density of the plasma is controlled by the pairs. Let us consider a proton population and adopt a power-law description for their number density distribution, represented as n p (γ) = k p γ -p , where k p and p denote the normalization and power-law index respectively. Assuming that the power-law index would be p ∼ 1, and the distribution would extend up to γ (syn) p,max (as indicated in eq. (C.2)). We can then set an upper limit on the number density of protons in the region by considering that the energy density in relativistic protons cannot exceed the magnetic energy density available for dissipation. This leads to the following constraint on the proton density. We can demonstrate that by equating the magnetic field energy density with the proton energy density, the ratio of pairs to protons density in this system becomes n p /n e ± < 0.25 • 10 -7 σ e,4 (18.4 + ln(γ
p,max,8 , where we normalized the maximum proton Lorentz factor to 10 8 (motivated by figure 8). Consequently, if protons were accelerated into a hard power-law up to the radiation-limited Lorentz factor, their number density should be negligible compared to the pair density.
Throughout this work, we have assumed that the magnetization in the upstream region remains below the synchrotron radiation-limited Lorentz factor of pairs γ rad (σ e < γ rad ), in which case the scenario of particle acceleration during the free phase can operate. Recently [41] have discussed the production of TeV flares from pairs created via γγ pair production, in the current sheets of M87*. Despite the similarities to our work, there is a major difference regarding the regime of interest; ref. [41] focuses on cases where the initial magnetization of the upstream plasma (i.e. before pair enrichment) is much higher than γ rad because the magnetospheric region is assumed to have a pair density equal to the Goldreich-Julian value (i.e. very low-density plasma region threaded by strong magnetic fields). In this regime, particles are mainly accelerated in the downstream reconnection region, with acceleration occurring predominantly in regions where the magnetic field strength is less than the electric field strength, as detailed in [12]. There, particle Lorentz factors may exceed γ rad due to anisotropic effects (i.e. small pitch angles) and potentially attain energies up to a few σ e m e c 2 , contingent upon the strength of synchrotron cooling. In this scenario, γγ pair production takes place on synchrotron photons emitted by a single particle population. The calculations of multi-wavelength spectra are then performed assuming a steady state where σ e is lower than its initial value as a result of pair enrichment. However, the pair density is calculated using optical depth arguments and does not incorporate dynamic effects. This might underestimate the total pair enrichment as shown in appendix D.
In this work, we have operated under the assumption of a one-zone approach, wherein both pair populations N free and N trap occupy the same spatial volume and emit nonthermal photons isotropically as discussed in section 2.3. However, to better describe real-world conditions, future investigations should aim to model pair production in scenarios where photons emitted by different pair populations have different geometrical distributions. Moreover, a more accurate calculation of pair production should incorporate interaction angles of photons and account for the production of pairs at different distances from the current sheet. Such an approach would offer a more nuanced understanding of secondary pair production rates and pair enrichment of the current sheet. Additionally, incorporating anisotropic effects into pair distributions is crucial, especially when pairs are accelerated JCAP12(2024)009 above γ rad (i.e. γ rad < σ e ), as it can influence cooling processes [42] and the resulting spectral characteristics. Finally in 3.3 we discussed the maximum energy that a single proton can achieve in the current sheet of M87*. It is essential to integrate realistic proton distributions and to delve into the intricate dynamics of proton acceleration in the context of 3D reconnection of pair-proton plasmas. An understanding of whether the free phase of acceleration applies to protons and the factors governing the broken power law distributions observed in recent simulations [31] is essential. Finally, radiation losses as shown may play a role in the formation of the distribution. All these are essential to quantify the expected photon and neutrino emission arising from a proton population within the current sheet.
We presented a model of particle acceleration at current sheets formed in the magnetospheric region of a rotating SMBH. The particle distribution in the reconnection region consists of two distinct populations that are coined "free" and "trapped". Free electrons and positrons are accelerated while meandering between both sides of a transient magnetospheric current sheet and can produce synchrotron-powered MeV flares. The luminosity of such flares, which is independent of the upstream plasma magnetization, is about 20 percent of the Blandford-Znajek power, and the duration of the flares is determined by the lifetime of the current sheets. Trapped electrons and positrons produce lower energy synchrotron radiation, which can extend down to a few MeV energies, and can produce X-ray counterparts to the MeV flares. Secondary pairs created via the attenuation of ≳ MeV photons by lower energy synchrotron photons may enrich the plasma, thus reducing the initial upstream magnetization by a factor of 100 within the lifetime of the current sheet. Low-energy synchrotron photons from trapped pairs are important targets for pion production and can limit proton acceleration to EeV energies for parameters relevant to M87*. N trap fr can be approximated by a broken power law with slopes ∼ -2 and ∼-3 below and above σ e respectively (see eq. (2.32)). The corresponding synchrotron spectrum will then be L syn (ϵ) ∝ ϵ -1/2 and ∝ ϵ -1 below and above ϵ syn (σ e ) respectively. In the context of M87* the luminosity of target photons can be approximated as L l ≈ L syn (ϵ l /ϵ syn (σ e )) 1/2 if ϵ syn (σ e ) > ϵ l ⇔ σ e > 10 5.2 , or L l ≈ L syn otherwise.
The optical depth for γγ pair creation in these two scenarios can then be written as assuming that L syn is given by (2.22) and that photons occupy a spherical volume of radius R eff . The latter assumption is made in order to compare it with the numerical results of LeHaMoC.
In figure 10 we plot the optical depth computed numerically using LeHaMoC for two values of ṁ (all other parameters are the same as in table 1), which agrees well with the analytical estimate. The inclusion of photons produced by trapped pairs as targets for γγ absorption increases the optical depth of the interaction and thus the pair creation rate within the system's volume.
The proton synchrotron energy loss timescale is defined as,
The energy loss time scale of a high energy proton with Lorentz factor γ p due to photopion production is given by [43],
where εth = 145 MeV, n γ (ϵ) ≡ dN γ /dV dϵ denotes the differential density of photons in the source frame of reference, while all barred quantities are measured in the proton's rest frame.
The energy-dependent cross-section and the inelasticity (fraction of energy transferred per collision to secondaries) of the process are represented by σ pπ and k pπ , respectively. The energy loss timescale for Bethe-Heitler pair production of a high energy proton with Lorentz factor γ p is given by [44] t pγ,ee ≈ γ p dγ p dt
-1 pγ,ee ≃ 3σ T cα f m e 8πγ p m p ∞ 2 dk n γ k 2γ p ϕ(k) k 2 -1 , (C.4)
where α f is the fine structure constant, σ T is the Thomson cross-section, and k = 2γ p ϵ/(m e c 2 ) is the photon energy in the proton's rest frame in units of the electron rest mass energy.
We observe that for certain combinations of ṁ and σ e the density of secondary pairs generated in the system through γγ absorption may exceed the primary density in a steady state (see figure 9). In order to assess whether pair creation influences the system dynamics, we conducted numerical calculations using the radiative code LeHaMoC. In these numerical experiments, we allowed the magnetization σ e of the system to vary over time. The recalculation of σ e was done every 10R eff /c in order to let the system approach its new steady state after each modification of the magnetization. First, we can calculate analytically the new magnetization after each pair enrichment event using eq. (3.2). If we initiate the numerical calculation with magnetization σ e,0 , which corresponds to a number density n e ± ,0 , then after 10R eff /c when the system will be close to a steady state, the total number density of pairs will be, n e ± ,1 = n e ± ,0 (1 + η r (σ e,0 )) (D.1)
where η r ≡ n sec γγ /n e ± which is a function of σ e (see eq. (3.2)). Modifying the magnetization accordingly we find that, σ e,1 = σ e,0 n e ± ,0 n e ± ,1 (D.2)
By repeating the same procedure we find that the magnetization in the step i + 1 of the modification will be given by,
In the analytical model, we neglect the injection of secondary pairs below γinj from γγ pair production. Their contribution will be assessed in the next section through numerical calculations.
https://github.com/mariapetro/LeHaMoC/
We note that the trapped pair population reaches a steady state only asymptotically. In practice, after ∼ 10R eff /c the densities reach 99.9954 percent of their asymptotic values.
All flux measurements, besides the EHT measurement at 230 GHz, which is spatially resolved, should be considered as upper limits on the flux in the reconnection region as they are likely produced at larger distances
The photon index Γ is defined as the slope of the differential photon energy spectrum and is calculated by performing a linear fit to the quantity log 10 (νFν ) ∝ (-Γ + 2) log 10 (ν) as derived by the model in the desired energy range
For the number density of the target photon field we assume that is given from n l = 3L pk syn /(4πR 2 cϵ l ) where in the calculation we assumed that R = Rrg.
The value at which the initial magnetization σe,0 does not get regulated by pair enrichment depends on the system parameters.