Blandford-Znajek jets in MOdified Gravity

General relativity (GR) will be imminently challenged by upcoming experiments in the strong gravity regime, including those testing the energy extraction mechanisms for black holes. Motivated by this, we explore magnetospheric models and black hole jet emissions in Modified Gravity (MOG) scenarios. Specifically, we construct new power emitting magnetospheres in a Kerr-MOG background which are found to depend non-trivially on the MOG deformation parameter. This may allow for high-precision tests of GR. In addition, a complete set of analytic solutions for vacuum magnetic field configurations around static MOG black holes are explicitly derived, and found to comprise exclusively Heun's polynomials.


Introduction
In recent years new astrophysical observations have given striking confirmations to the prediction of General Relativity (GR).The detection of gravitational waves by the LIGO-Virgo-KAGRA collaboration provided the first direct observations of coalescing binary systems of black holes and neutron stars [1][2][3], whereas the images revealed by the Event Horizon Telescope (EHT) collaboration confirmed the presence of supermassive black holes harboured in the nuclei of galaxies [4,5], and gave support to the hypothesis that the astrophysical jets associated to Active Galactic Nuclei (AGNs) are powered by spinning black holes via the Blandford-Znajek (BZ) mechanism [6][7][8].
The new observational opportunities at our disposal pave the way to the detection of strong-gravity effects in a variety of astrophysical systems, where potential deviations from GR are expected to manifest.In this regard, gravitational waves detectors and telescope arrays such as the EHT not only constitute a source of new observational discoveries, but also an invaluable ground where to make comparison between precise theoretical predictions of GR and alternative theories of gravity.For instance, the first multimessenger observation of a neutron star merger/GRB event, GW170817/GRB170817A [9], not only confirmed binary neutron star mergers as progenitor of short gamma-ray bursts, but also allowed to impose stringent constraints that ruled out all the theories which are not consistent with the weak form of the equivalence principle [10].Similarly, the comparison with EHT observations of the M87* shadow permitted to exclude many models of "black hole mimickers" as well as black holes with additional massive scalar hairs as candidates for the central object [7].
Among the candidates that survived observations, the Scalar-Tensor-Vector Gravity (STVG) [11] -also referred to as MOdified Gravity (MOG) in the literature -provides a covariant extension of GR in which electromagnetic and gravitational perturbations propagate with the same speed, thus being consistent with the GW170817/GRB170817A event [12], and admitting spinning black hole solutions [13] whose shadows are qualitatively indistinguishable from the one of M87* within the present EHT precision [7].From a theoretical perspective, the interest for MOG theory stems from the fact that it constitutes a covariant generalisation of GR based on an action principle that offers phenomenological explanations to the galactic rotation curves [14,15], the bullet cluster phenomenon [16] and cluster dynamics [17], without the necessity of introducing dark matter.In light of this, it is urgent to develop new theoretical signatures that can help to distinguish GR from MOG in present and future astrophysical observations.One of the most characteristic features of spinning black holes is the possibility of extracting rotational energy via the Penrose process [18,19], whose electromagnetic manifestation, namely the BZ mechanism [6], is considered to be the primal engine behind the emission of relativistic jets in AGN.Despite its great importance in the field of relativistic astrophysics, the inherently complex dynamics behind the BZ mechanism constituted a major obstacles for theorists to clarify details about its physics, and only in the last two decades progress in this directions were made either by means of numerical simulations [20][21][22][23][24], and by means of analytic studies [25][26][27][28][29][30][31][32][33][34].
Recently there have been efforts to understand how deviations from GR can manifest in the physics of relativistic jets.Examples range from studies on emission mechanisms in regular black hole metrics [35], Kerr-Sen black holes [36], quadratic and cubic theories of gravity [37,38], and on signatures from extra dimensions [39].A first attempt to investigate how MOG deformations reflect on jet dynamics has been conducted in [40], where the Blandford & Payne mechanism [41] was assumed to cause the jet launching and particle trajectories around Kerr-MOG black holes were computed.
In this work we consider an alternative path with respect to [40].We explore the Kerr-MOG scenario by regarding the jet as being powered by the black hole, with the BZ mechanism responsible for the extraction of energy and angular momentum, and by focusing on the dynamics of force-free magnetospheres around MOG black holes.
Indeed, the comparison between EHT observations and General Relativistic Magneto-HydroDynamics (GRMHD) simulations revealed that the power emitted in the launching region of M87* is consistent with the BZ mechanism [7,8].More specifically, numerical results from GRMHD simulations [23,42] showed that the rate of energy extracted in the BZ mechanism from a black hole surrounded by a razor-thin accretion disc can be generally expressed as where the quantity (2πΨ H ) represents the total flux threading the event horizon, Ω H is the black hole angular velocity, and κ is a factor taking into account the geometry of the magnetic field configuration (κ = 2π 3 • 1 4π 2 ≈ 0.053 for a monopolar field).The function f (Ω H ) is approximately f (Ω H ) ≈ 1 for slowly-spinning black holes, reproducing the widely used quadratic scaling for ĖBZ computed originally in [6], whereas in the high-spin regime its explicit expression is typically given in terms of an expansion in powers of Ω H . Currently, the most accurate expression for f (Ω H ), obtained in the case of a monopolar magnetosphere around a Kerr black hole, was computed analytically up to orders O(Ω 6 H ) in [34], improving previous estimates derived by means of numerical fit with GRMHD simulations [42] which truncated the series at O(Ω 4 H ). It is crucial to stress that the authors in [42] gave indications that the function f (Ω H ) remains the same even for different magnetic field geometries.Being the BZ mechanism dependent on both the configuration of the magnetosphere surrounding a black hole and the theory of gravity considered, one expects that the power extracted from black holes in alternative scenarios can constitute a signature to reveal deviations from standard GR results [35,43,44].In what follows we argue that such deviations are mainly encoded in the function f (Ω H ). The main purpose of this work is indeed to derive an explicit expression for f (Ω H ) in the MOG case beyond the leading order approximation, and show that this function depends in a non-trivial manner on the MOG deformation parameter that accounts for deviations from GR.In doing this we apply the standard BZ perturbation theory to construct analytic models of black hole magnetospheres and to explicitly compute the energy extracted via the BZ mechanism, along the lines of [30] and [34].
More generally, one expects that the expression of f (Ω H ) is characteristic of the underlying theory of gravity on which the BZ mechanism is set to operate.Hence, further investigations on the BZ mechanism have the potential to produce new theoretical predictions and observational signatures that would enable current and future horizon-scale measurements done by the EHT collaboration to distinguish GR from alternative theories of gravity, test the Kerr paradigm and clarify details concerning the structure of magnetic fields around spinning black holes.
While future high-precision polarimetric observations with EHT are expected to put more stringent constraints on the magnetic flux and magnetic field geometry threading the black hole [7], an accurate knowledge of f (Ω H ) becomes crucial if one aims to use power estimates in order to discriminate between jet models in GR and modified theories of gravity.As mentioned earlier, in alternative scenarios one might expect that f (Ω H ) will depend on one or more additional deformation parameters, considerably increasing the parameter space of the theory.Therefore, it is important to pursue an analytical derivation of f (Ω H ).This is also motivated by the fact that the role of f (Ω H ) becomes relevant in the high-spin regime where GRMHD simulations are computationally expensive [45], and analytic results can give important insights for the development of future numerical codes.
The present work is structured as follows.In Sec. 2 we review some of the main features of STVG and of its Kerr-MOG solution for spinning black holes, whereas in Sec. 3 we collect standard results and equations that captures magnetospheric dynamics in stationary and axisymmetric backgrounds, which find application in the case under study.In Sec. 4 vacuum magnetic field configurations around static MOG black holes are classified.This will serve as the starting point for the perturbative construction of a spinning monopolar magnetosphere around Kerr-MOG, which we detail in Sec. 5.An explicit expression for the power and angular momentum extracted in the BZ mechanism, together with a first-order expression for the factor f (Ω H ) in terms of the MOG deformation parameter is given in Sec.6, where we also make comparisons with the usual BZ theory in GR.We conclude the paper with a discussion in Sec. 7.
Where not explicitly specified we adopt geometrised units by setting G N = c = 1 for the Newton's constant and the speed of light, and signature (−, +, +, +) for the metric.

MOdified Gravity (MOG)
We provide in this section a short review of the STVG theory, aimed at clarifying its main features and the motivation behind its construction.We first list the various fields characterising this theory and write the action explicitly.We then focus on presenting the Kerr-MOG black hole, a stationary and axisymmetric solution of the STVG equations of motion which, for the rest of the paper, will be regarded as the background to construct magnetospheric models and study the BZ mechanism.

Scalar-Tensor-Vector Gravity action
STVG has been developed with the aim of providing a fully relativistic generalisation of GR, in which the weak-field modified gravity effects become appreciable over galactic scales, thus explaining astrophysical phenomena, as for example the galaxy rotation curves, without postulating the existence of non-baryonic dark matter.
The main idea of STVG consists in enhancing the gravitational coupling constant G to a dynamical scalar field whose asymptotic value exceeds the Newton constant, G N = 6, 67 × 10 −11 Nm 2 /kg 2 .In order to compensate for the increased value of G at scales comparable to the solar system, where GR provides an accurate description of the physics we observe, a repulsive Yukawa interaction is included in the theory.This is achieved by adding a Proca massive vector field ϕ µ to the action, so that, at short scales, the value of G reduces to G N and STVG effectively reconciles with GR and its Newtonian limit without violating any known local observations.One of the main advantages of MOG theory is that it descends from an action principle.The action consists of four contributions [11,15,17] where with R being the Ricci scalar computed from the dynamical metric g µν .The additional contributions S ϕ and S s are respectively associated to the vector and scalar fields of the theory, explicitly given by1 (2.3) In the expressions above the quantity B µν ≡ 2∂ [µ ϕ ν] represents the field strength for ϕ µ (x), coupled to the dynamical scalar field µ(x), whose vacuum state plays the role of an effective mass for the vector field, and set the scale at which the deviations from GR become appreciable.As mentioned earlier, in STVG the effective gravitational coupling is promoted to a dynamical scalar field G(x), whose value is in general allowed to vary in space and time.
The quantities V (ϕ), V (G) and V (µ) denote possible self-interactions of the fields.Finally the term S M represents the action contribution associated to source matter fields, which are assumed to couple with both the metric g µν and the vector field ϕ µ .

Kerr-MOG black holes
A vacuum static and spherically symmetric black hole solution for STVG has been found in [11,46] by solving the field equations associated to the action (2.1).The solution was later generalised to the stationary and axisymmetric case in [13], and denoted as the Kerr-MOG black hole.In obtaining such solutions one neglects self-interaction potentials by setting V (ϕ) = V (G) = V (µ) = 0, and fix constant values for the scalar fields G and µ.In particular, one assumes that the asymptotic value of the enhanced gravitational constant, The asymptotic solution for the vector field ϕ µ reveals that its time component is a Yukawatype potential [11] whose charge K is postulated [13] to acquire a gravitational character by means of the relation in order to recover the Newtonian value for G in the proximity of a source mass M .Assuming a minimal coupling between test particles and the vector field, it was shown that the repulsive Yukawa term (2.5) adds to the Newtonian gravitational potential, leading to a modified Newtonian dynamics in the weak field regime of the MOG background [11,15].The mass value µ can be fixed in such a way that effects of deviations from GR -namely the cut-off scale for the Yukawa potential -manifest at distances of kiloparsecs from the gravitational source.The value µ ≈ 2.6 × 10 −28 eV, corresponding to a lengthscale ℓ = µ −1 ≈ 23.8 kpc, proved to correctly reproduce the galaxy rotation curves [14,15], the bullet cluster [16] and cluster dynamics [17], but can be effectively neglected in the derivation of black hole solutions [13].Under this approximation the vector field leads to an additional repulsive Coulomb contribution to the gravitational potential.
The Kerr-MOG spacetime [13] is a stationary and axisymmetric three-parameters family of solution for STVG, which describes a black hole characterised by an angular momentum and an additional gravitational charge reflecting the presence of the additional vector field.The spacetime metric g µν and the vector field ϕ µ can be expressed in Boyer-Lindquist coordinates x µ = (t, r, θ, φ) in terms of the line element ds 2 = g µν dx µ dx ν , namely where, adopting units for which G N = 1, one has with a being the angular momentum per unit mass.Accordingly, the Arnowitt-Deser-Misner (ADM) [47] mass, angular momentum and gravitational charge for the Kerr-MOG black hole [48] are defined as It is important to remark that, while the parameter M only serves as an integration constant of the STVG field equations, it is M α that corresponds to the gravitational mass of the Kerr-MOG black hole, in accordance with the ADM Hamiltonian formulation [48].We notice that the presence of a non-vanishing MOG deformation parameter α, defined in Eq. (2.4), leads to a difference between M and the ADM mass M α .
In the following it will be convenient to introduce a dimensionless spin parameter, given by ϵ and for the rest of the paper we assume a, or analogously ϵ, to be non-negative.In the limit ϵ → 0 the metric (2.7) describes a static MOG black hole [11,13], whereas for α → 0 it smoothly reduces to the Kerr solution of standard GR.The Kerr-MOG metric (2.7) admits an inner and an outer horizon, denoted here as r ± respectively, solutions of ∆ = 0, and respectively located at It is possible to associate an angular velocity to the event horizon r + , given by 12) The boundary of the ergoregion, or ergosurface, can simply be defined as the locus of spacetime points where the asymptotic timelike Killing field ∂ t becomes null, i.e.where g tt = 0. Consequently, the ergosurface r e (θ) is found at By requiring the gravitational charge K in Eq. (2.9) and the event horizon position r + as defined in Eq. (2.11) to be real quantities, in accordance with the cosmic censorship conjecture, one obtains physical bounds on the MOG deformation parameter α [48,49] 0 In particular extremality, is defined when the equality holds in the second relation above, meaning that the two horizons collapse in a single spherical surface of radius r ± = M α = |a| √ 1 + α.Hence, a non-vanishing value for the deformation parameter α implies that, for the same ADM mass and angular momentum, the size of the horizon for a Kerr black hole is always bigger than the one of its MOG counterpart.Notice also that, unlike in the Kerr metric, having that the specific angular momentum equates the ADM mass is not possible in Kerr-MOG.In order to make a comparison with the situation in GR, we label respectively with M ⋆ and J ⋆ the ADM mass and angular momentum of a GR Kerr black hole.Equating the ADM mass with its Kerr-MOG counterpart means M ⋆ = M α and it is immediate to notice that, while for extreme Kerr J ⋆ = M2 ⋆ , one has J = M 2 α / √ 1 + α < J ⋆ in Kerr-MOG.Thus, by spinning up a Kerr and a Kerr-MOG black holes of same ADM masses one would reach extremality in the MOG scenario before than in the GR case.

Black hole magnetospheres dynamics
In this section we review conventions and equations that will be pivotal for extending the BZ model to the case of Kerr-MOG black holes.In particular, we assume that a vacuumbreaking process similar to that outlined by Blandford and Znajek [6] exists so that at equilibrium a plasma-filled magnetosphere, represented by a probe electromagnetic field F µν sustained by a current density of the plasma j µ , surrounding the Kerr-MOG black hole establishes.As a first exploration of the magnetospheric problem in a MOG background, we also assume that an eventual coupling between electromagnetic fields and the MOG vector field ϕ µ is negligible 2 .This is to avoid that the Kerr-MOG background can act as a source for additional electromagnetic fields, and guarantees that Maxwell's equations are still given by ∇ [ρ F µν] = 0 and ∇ ν F µν = j µ , with ∇ µ being the covariant derivative defined with respect to the Kerr-MOG metric in Eq. (2.7).

Force-Free Electrodynamics (FFE)
The minimum non-trivial level of description for magnetospheres around compact objects like black holes and neutron stars is correctly captured by the theory of Force-Free Electrodynamics (FFE) [6,51,52], a regime of magnetohydrodynamics in which most of the energy resides in the electromagnetic sector of the system and the inertia, as well as the thermal degrees of freedom, of the plasma can be effectively neglected [53].Under these assumptions the magnetic fields dominates the energy and momentum balance of the system, and the stress-energy tensor respects the force-free approximation T µν ≃ T EM µν .The FFE system of equations is therefore given by the Maxwell's equations ∇ [ρ F µν] = 0, and ∇ ν F µν = j µ supplemented with the force-free condition By means of the non-homogeneous Maxwell's equations, and assuming the field to be magneticallydominated, F µν F µν > 0 [21,51,52], it is possible to formally eliminate the current j µ from the dynamical description and summarise the FFE equations as Accordingly, FFE is a non-linear regime of Maxwell's electrodynamics which is only specified by the electromagnetic field configurations, whereas the current j µ can be regarded as a secondary quantity that descends from the fields, which only serves to sustain the magnetic fields and screen the longitudinal component of the electric fields.
We consider magnetospheres around a Kerr-MOG black hole that share the same symmetries of the background, i.e. that are stationary and axisymmetric [6,54].As usual in magnetohydrodynamics, stationary and axisymmetric flows are characterised by a set of integrals of motion taking constant values along the magnetic field lines [55].More specifically, the magnetic flux across a circular loop of radius r sin θ piercing the rotational axis of the black hole, Ψ = Ψ(r, θ), is constant along the projection of the field lines on the poloidal plane (r, θ).In the force-free limit the electromagnetic field is totally determined by two additional integrals of motion [52], respectively representing the poloidal current I = I(r, θ) flowing through a loop around the rotational axis, and Ω = Ω(r, θ), namely the angular velocity of the magnetic field lines.These two are subject to integrability conditions which descend from Eq. (3.2), and imply that I ≡ I(Ψ) and Ω ≡ Ω(Ψ) [6,34,53].
As anticipated, a generic stationary and axisymmetric force-free field in the Kerr-MOG background is written solely in terms of the integrals of motion, and reads [52] where the one-form is defined to be always orthogonal to the four-velocity of a particle which is co-rotating together with the field lines, with an angular velocity Ω(Ψ).
The magnetospheric problem amounts to find explicit functional expressions for the field variables Ψ(r, θ), I(Ψ) and Ω(Ψ) on a given background.

Grad-Shafranov equation and critical surfaces
Due to the structural analogy between Kerr and Kerr-MOG it is possible to obtain similar relations that relate the force-free field variables in the stationary axisymmetric case.In this section, therefore, we limit to record the standard set of relations used in the construction of magnetospheric models.For more detailed derivations and discussions we refer the reader to previous work in the literature [6,34,52,55].
As typical for stationary and axisymmetric flows in magnetohydrodynamics, the field variables of the magnetospheric problem are put in relation altogether by means of the Grad-Shafranov equation [53,55].In the Kerr-MOG spacetime this can be compactly written in a semi-covariant manner [34] as The equation above, also called stream equation, emerges by combining the r and θ components of the force-free condition, Eq. (3.2), in the Kerr-MOG background, together with the field structure (3.4) [34].Besides specifying boundary conditions, which allow to select a particular solution, the Grad-Shafranov equation must also be supplemented with regularity conditions at its critical surfaces [24,56,57].Similar to the case of the standard Kerr spacetime in Boyer-Lindquist coordinates, for Kerr-MOG black holes the magnetospheric structure is featured with the presence of four critical surfaces located at the event horizon, at asymptotic infinity and at two light surfaces.
In order to demand regularity of the electromagnetic field at the horizon and at infinity, the following two Znajek conditions [21,24,30,34] must be imposed where the labels +, ∞ have been respectively used to evaluate quantities as the radial coordinate approaches the horizon and infinity ( e.g.I + ≡ I| r=r + and I ∞ ≡ I| r=∞ ), and with Ω H being the black hole angular velocity, as given in Eq. (2.12).The light surfaces are defined as the locus of points where the velocity of an observer co-rotating with the field lines becomes null, thus satisfying the condition g µν η µ η ν = 0 [52].Such an equation generally admits two real and distinct solutions which respectively define the Inner/Outer Light Surface (ILS/OLS).At both light surfaces the following Robin-type condition [34,52,56,57] dubbed reduced stream equation, must hold.This ensures the smoothness of the magnetic flux Ψ when the magnetic field lines cross these critical surfaces.Detailed information about the properties of the ILS and OLS in the case of the Kerr background can be found in [21], and they directly extend also to the Kerr-MOG case due to the structural analogy between the two metrics.While the ILS is a closed surface that is always comprised between the event horizon and the ergosphere, and therefore is a distinguishing feature of black hole magnetospheres, the OLS is characterised by an open topology, that constitutes the black hole analogue of pulsar's light cylinder [24].
We conclude this section by stressing that, as in the GR case, Eq.s (3.6), (3.7) and (3.8), supplemented with appropriate boundary conditions which specify a particular topology for the magnetic field configuration, constitute the complete set of equations needed to address the magnetospheric problem.The authors of [30] first showed the crucial role of the light surfaces in the construction of a consistent solution for a monopolar magnetosphere around a Kerr black hole by exploiting perturbative techniques enhanced with a matched asymptotic expansion scheme.In [34] such a methodology has been improved and used to derive semi-analytical results leading to the computation of new perturbative contributions in the power extracted by means of the BZ mechanism, that solved previous discrepancies between numerical simulations and analytic approaches.The results of [34] are pivotal for a possible understanding of the non-perturbative structure of the BZ theory.For a comprehensive discussion on the Grad-Shafranov equation and the regularity conditions at the critical surfaces in the context of Kerr black holes we refer the reader to [34].

Vacuum solutions in static MOG backgrounds
The BZ perturbative approach [6,30,34] consists in building perturbations around vacuum solutions in static black hole backgrounds.It is therefore essential to derive and classify the vacuum electromagnetic field configurations surrounding a static MOG black hole, so as to use them as a starting point to extend the BZ perturbation theory to the case of a Kerr-MOG background.In this section we therefore focus on static MOG black holes which we can describe by setting the spin parameter ϵ = 0 in Eq. (2.7).The metric, thus, reads with ∆ = (r − r+ )(r − r− ) and the two horizons in the static case are located at r± = M α (1 ± (1 + α) −1 ).It worths noticing that, as a consequence of the fact that the "charge" appearing in the Kerr-MOG metric is tied to the gravitational mass, no extreme limit for a static black hole exists in MOG, as opposed to the Reissner-Nordström solution in GR.
Vacuum solutions, characterised by j µ = 0, can be constructed by assuming I(Ψ) = Ω(Ψ) = 0 and by demanding the magnetic flux Ψ to obey the Grad-Shafranov equation (3.6) in the static MOG background.The latter reduces to the homogeneous equation LΨ(r, θ) = 0, where the Laplace operator in the static MOG background reads By assuming a solution of the form Ψ(r, θ) ∼ R (ℓ) (r)Θ ℓ (θ) it is possible to separate the variables, with the radial harmonics R (ℓ) (r) and the angular harmonics Θ ℓ (θ) being eigenfunctions of two independent Sturm-Liouville problems [58], respectively defined by In the rest of the paper we consider split-field configurations [52], and only focus on the domain θ = [0, π/2].Solutions in the entire space can directly be obtained by reflection across the equator θ = π/2.In what follows we obtain analytic expressions for the solutions of the two Sturm-Liouville problems (4.3); in particular, the magnetic flux solution of the vacuum equation is constructed by combining the radial and angular harmonics, according to Ψ(r, θ) where an integration constant c 0 is included and the summation on the index ℓ is implicit.
The regular solutions we consider are subject to the following boundary conditions [24,30,34] with the first two conditions demanding regularity of the field at θ = 0, and the last being a normalisation condition on the total magnetic flux passing through a surface encompassing the event horizon.Later in this section we show that an additional boundary condition in the asymptotic region is needed in order to specify the topology of the vacuum magnetic field [29].

Angular eigenfunctions
Eigenfunctions for the angular part of the Laplace operator (4.3) have been extensively studied in the literature [58], and here we limit to summarise their main properties.Angular harmonics obeying L (ℓ) θ [Θ ℓ ] = 0 where ℓ is a positive integer, reduce to a degenerate case of the hypergeometric functions [59].One of the solutions is characterized by its irregularity at θ = 0. On the other hand, the second solution takes the form of a truncated polynomial that terminates after a finite number of terms.These two types of solutions represent distinct mathematical representations of the behavior of the system.The choice between these solutions will be based on the regularity condition at θ = 0, Eq. (4.5).
The regular solution Θ ℓ can be directly expressed in terms of Gegenbauer polynomials [58,60].More specifically where the C α ℓ (x) are orthogonal polynomials in the interval x ∈ [−1, 1] with respect to the weight function Below we list the explicit expression for ℓ = 0, 1, 2, 3, which are the only relevant for the present discussion The first solution, in particular, is related to the monopole field configuration that we will mostly concentrate on in this paper.

Radial eigenfunctions
In order to study the radial part of the Laplace operator, it is convenient to redefine r in terms of a dimensionless radial coordinate This shows that the radial equation r [R (ℓ) ] = 0 is actually a second-order linear Ordinary Differential Equation (ODE) of the kind where ∆(w) = r2 + (w − 1)(w − w a ) with w a < 1, and This makes manifest the fact that Eq. (4.10) is in general characterised by four regular singular points, located at w = 0, w a , 1, ∞, with 1 > w a ≥ 0 due to the coordinate choice adopted.
It is immediate to notice that the Laplacian operator in the standard GR case (α → 0) is recovered by simply taking the limit r− → 0 and r+ → 2M , in such a way that w a → 0.
Every linear second-order ODE with at most 4 regular singularities in the complex plane can always be reduced to the Heun's Equation [61] where the complex number q is called the accessory parameter, and with each of the 4 regular singularities w = 0, 1, w a , ∞ related to a pair of characteristic exponents, respectively given by (0, 1 − γ), (0, 1 − δ), (0, 1 − ε) and (λ, χ).In order to guarantee that infinity constitutes a regular singular point, the parameter ε is subject to the Fuchsian constraint [61] so that the Heun's equation is specified by the set of 6 parameters w a , q, δ, γ, λ, χ.We notice that the Heun's equation and its confluent versions have often found applications in black hole physics (see Ref.s [62][63][64][65][66][67] for a partial literature).
In the case considered here, Eq. (4.10) directly reduces to the Heun's equation (4.12) upon identifying As mentioned above, in the limit w a → 0 the standard Schwarzschild case in GR can be recovered and, accordingly, the Heun's equation reduces to the Papperitz-Riemann equation, whose solutions are given in terms of hypergeometric functions [58].
Similarly to what one does in the standard Schwarzschild case, we look for two families of polynomials, which arise as Frobenius solutions of the Heun's equation (4.12).There are 2 × 4 × 6 Frobenius solutions at each singular point, for a total of 192 Frobenius' solutions [68].The symbol H l (w a , q, λ, χ, δ, γ; w) is typically used to label the solution at w = 0 with characteristic exponent 0. All the remaining solutions can be determined by acting with Möbius and indeces transformations on the Heun's equation (4.12).Heun's polynomials in particular arise as particular solutions which are regular at three singular points [69].One also demands that the correct solutions are capable of reproducing the eigenfunctions for the Laplacian operator in the standard Schwarzschild background upon taking the limit w a → 0. In particular, the two classes of solutions respecting all the aforementioned characteristics can be written as follows for ℓ ≥ 1 We stress that this are new and original results obtained in this work, representing radial harmonics for vacuum electromagnetic fields around black holes with a Reissner-Nordström metric structure.One can check explicitly that for w a → 0 the standard Schwarzschild vacuum fields in terms of hypergeometric functions are recovered [58].More generally, the function U ℓ (w a , w) admits a truncated power series representation, akin to whose coefficients a n (w a ) obey the two-terms recurrence relation where From the series representation above one can explicitly compute the polynomials.For later convenience we list below the expressions of the polynomials for both U ℓ and V ℓ corresponding to ℓ = 1, 2: We notice that at the event horizon, w → 1, the functions U ℓ are properly defined, whereas V ℓ manifests a logarithmic divergence Viceversa, for w → ∞ the asymptotic values are Hence the vacuum radial functions U ℓ are convergent at the horizon and divergent at infinity, whereas the opposite holds for the funtions V ℓ .

Asymptotically monopolar static field
As an example of vacuum magnetic field solutions around static MOG background (4.1), we can consider configurations with an isotropic distribution of magnetic field lines at infinity.A monopolar magnetosphere is specified by the following asymptotic boundary condition [29] lim Among the exact solutions we derived, there exists a monopole magnetic field solution given by Ψ(r, θ) where the integration constant has been chosen in order to respect the normalisation condition of the flux in Eq. (4.5).
A monopolar magnetosphere around a static MOG black hole is thus indistinguishable from one surrounding a Schwarzschild black hole in GR.However, as we detail in the next section, this is no longer true for a monopole magnetosphere around spinning black holes, with the static solution constructed here that will serve as the starting point to extend the BZ perturbative procedure to the Kerr-MOG case.For comparison in the same plot the case of a vertically uniform field (dashed black) around a Schwarzschild black hole in GR -whose horizon is depicted in gray -is reported.The plot assumes α = 10 to magnify the distorsion of the magnetic field lines in the vicinity of the non-spinning MOG black hole.

Asymptotically vertical static field
Before moving to consider the case of a stationary MOG background we present here another interesting topology for a vacuum field.Following the definition given in [29], a vertical topology [70] is defined by the condition that Ψ remains finite for r → ∞ and the product r • θ is kept fixed.This can be more easily visualised in cartesian coordinates, r = √ x 2 + z 2 and cos θ = z/ √ x 2 + z 2 , with the asymptotic boundary condition that now reads By considering vacuum solutions of the type Ψ(r, θ) = U ℓ (r)Θ ℓ (θ) which are regular at the event horizon, and by converting to cartesian coordinates, it is possible to observe that for large values of z one has Ψ(x, z) ∼ x 2 z ℓ−1 .Thus, the only vacuum solution which is regular both at the event horizon and at infinity and that is consistent with the asymptotically vertical boundary condition (4.24), corresponds to ℓ = 1 and explicitly reads where an integration constant was fixed by means of the normalisation condition in Eq. (4.5).We notice that for α → 0 one recovers a vertical magnetic field configuration in the Schwarzschild background with profile given by Ψ = r 2 /(2M ) 2 sin 2 θ [70].
An interesting result of this paper is that, in contrast with the GR case, vacuum magnetic fields surrounding static MOG black holes, and characterised by a vertical asymptotic profile, are not uniformly vertical.In Fig. 1 we present an illustrative comparison between the asymptotically vertical magnetic field for a Schwarzschild and a non-spinning MOG black hole with the same ADM mass.The derivation obtained here shows that, interestingly, deviations from GR not only affect the dynamics at large scales, but can also be observed in strong-gravity magnetostatic configurations.This might suggest that also the geometrical factor κ in the BZ formula (1.1) for the energy extracted can receive modification by the presence of the MOG parameter α, as opposed to the case of a monopole magnetosphere where only the factor f (Ω H ) will change.

Force-free Kerr-MOG magnetospheres
In this section we generalise the BZ approach [6,30,34] for the construction of a splitmonopole magnetosphere in the Kerr-MOG background in the regime of slow rotation, namely for small values of the black hole spin ϵ.At each order in perturbation theory we provide explicit expressions for the magnetospheric field variables, and we comment about the regularity of the solution across all the critical surfaces that characterise the magnetospheric problem.

Leading order solution
Since in the limit ϵ → 0 the Kerr-MOG metric (2.7) reduces to the static MOG metric (4.1), a perturbative solution of the Grad-Shafranov equation (3.6) for small values of ϵ can be constructed order by order in an ϵ expansion, by starting from a vacuum field in the static MOG background, F µν ∼ O(ϵ 0 ), with an associated current which is assumed to scale as j µ ∼ O(ϵ).At the leading order in the expansion for small ϵ then the field and the vector current automatically satisfy the force-free constraint, F µν j ν ∼ O(ϵ).These assumptions are satisfied as long as the force-free field variables scale as In this section we only consider the split-monopole configuration which has been derived in the previous section after specifying the boundary conditions (4.5) and (4.22).
With the assumptions (5.1), the condition g µν η µ η ν = 0 that determine the light surfaces reveal the following scaling in the spin parameter ϵ In other words, at the leading order, the ILS coincide with the event horizon, whereas the OLS is located at infinity.The non-perturbative scaling of the OLS is at the core of the Matched Asymptotic Expansion scheme that is needed to consistently construct higher-order solutions in the BZ perturbation theory [30,34].In the case of a split-monopole field in Kerr the necessity of the MAE scheme becomes evident at the fourth order in perturbation theory [30,34].We expect the same to be true in Kerr-MOG.However, in order to disentangle possible modified gravity effects from specific configurations of the magnetosphere [37] it turns out to be sufficient to truncate the series at the third order in the perturbative expansion, so that in the present work we do not need to deal with the MAE scheme.
All the quantities will be normalised to the black hole gravitational mass M α , which is independent from the expansion parameter ϵ and will facilitate the comparison with results for Kerr black holes with same ADM mass as the Kerr-MOG black hole under consideration.

First order in the small spin regime
The ansatz for the field variables expansion at the first subleading order is , where a factor proportional to the ADM mass of the black hole is included so as to have dimensionless perturbative coefficients, and the integrability conditions (3.3) have been used to constrain the dependence of i 1 and ω 1 .

Magnetic flux
The function ψ 1 obeys the source-less stream equation Lψ 1 = 0, and the unique solution which is regular at both the horizon and at infinity is the trivial solution ψ 1 = 0.

Poloidal current and angular velocity of the field lines
By means of the Znajek condition (3.7) at the event horizon one obtains (5.5) The asymptotic Znajek condition (3.7) instead gives By comparing the two conditions above one obtains explicit expressions for the current and the angular velocity at this order in perturbation theory (5.7) Notice that in the limit α → 0, these two quantities reproduce the results known in the literature for the Kerr spacetime, namely i 1 = 1/4 Θ 1 (θ) and ω 1 = 1/4.
The expressions in Eq.s (5.5) and (5.6) can also be derived by demanding regularity of the stream equation respectively at the ILS and at OLS, by means of (3.8) [34,71].As already explained, for ϵ = 0, the ILS and the OLS are located at the event horizon and at infinity respectively.Turning on the spin parameter ϵ, it is possible to compute perturbatively the new location of the ILS and OLS.For the first correction we find By evaluating the reduced stream equation (3.8) at the new location for the ILS and OLS, one gets an equation for i 1 and ω 1 which is precisely solved by the expressions given in Eq. (5.7).Thus, all the quantities defined at this order in perturbation theory are regular at all the critical surfaces.

Second order in the small spin regime
Moving to the next perturbative order, the expansion now reads where the dependence on i 2 and ω 2 can be directly inferred from the integrability conditions (3.3).

Magnetic flux
Expanding the stream equation (3.6) up to O(ϵ 2 ) and using the known expressions for ψ 0 , ω 1 and i 1 , one obtains the following equation for the function ψ 2 (r, θ) where we recall that r± = M α 1 ± (1 + α) −1 label the outer/inner horizon positions in the case of a static MOG black hole (see Sec. 4).
By assuming the solution to be separable, ψ 2 (r, θ) = R 2 (r)Θ 2 (θ), one can project this equation on Θ 2 (θ) by using the orthogonality conditions which characterise the angular harmonics Θ ℓ (θ).This produces a non-homogeneous differential equation for the radial part which can be conveniently written in terms of the dimensionless radial coordinate w = r/r + and w a = r− /r + , introduced in Eq. (4.9), as ( 5.11) This equation is the same type of Heun's equation that we studied in the vacuum static case (See Eq. (4.12)), with ℓ = 2 and with an additional non-homogeneous term.It is possible to determine analytically a particular solution to this equation which is regular at both the event horizon (w = 1) and at infinity after properly adding the homogeneous solutions U 2 (w a ; w) We illustrate in Fig. 2 the behaviour of the function R 2 (w) for some representative values of the deformation parameter α.Notice that by taking the limit w a → 0, and after using the inversion formula for the dilogarithms, one has which is precisely the result known in the Kerr spacetime.The function R 2 (w) is smooth in the asymptotic region which trivially reproduces the result known in the Kerr spacetime R ∞ 2 = M 4r in the limit ω a → 0 [6,34].This means that the solution at O(ϵ 2 ) extends all the way up to the asymptotic region, smoothly crossing the OLS.
By making use of the inversion relation one can also infer the behaviour at the static horizon w = 1, i.e. r = r+ = 1 + (1 + α) −1 .This simply reads (5.15) Remarkably, the limit w a → 0 is finite and reproduces the value obtained expanding up to second order in the spin parameter in the case of a monopole solution for the Kerr magnetosphere, R H 2 = 6π 2 −49 72 [34].We also recall that explicitly one has w a = 1+α− √ 1+α 1+α+ √ 1+α for the Kerr-MOG spacetime.

Poloidal current and angular velocity of the field lines
The Znajek condition at the horizon, expanded at the order ϵ 2 , gives (5.16) Similarly, the Znajek condition at infinity at the order ϵ 2 gives By comparing the two relations above, one immediately has for consistency that The corrections to the positions of the light surfaces can be obtained by demanding that r erg (θ) ≥ r ILS (θ) ≥ r + which, together with Eq. (5.18), leads to the following explicit expressions for the light surface positions at order ϵ 2 (5.19)

Third order in the small spin regime
By including terms of order ϵ 3 , the expansion of the force-free field variables reads Again, the integrability conditions I = I(Ψ) and Ω = Ω(Ψ), Eq. (3.3), dictate the structure of the expansion at the third order in perturbation theory.

Magnetic flux
Similarly to what we saw at the first perturbative order, the function ψ 3 obeys a sourceless stream equation Lψ 3 = 0.This means that the trivial solution ψ 3 = 0 is the only one consistent with the boundary conditions of the split-monopole configuration.

Poloidal current and angular velocity of the field lines
Expanding the Znajek condition at the horizon (3.7) up to order ϵ 3 , one obtains a relation between i 3 and ω 3 of the type (5.21) The asymptotic Znajek condition at the third perturbative order, instead, simply produces By comparing the two equations above one obtains the following expression for i 3 (θ) while for ω 3 (θ) we have In the limit α → 0 Eq.s (5.23) and (5.24) correctly reproduce the results for the Kerr magnetosphere, namely [34] i where in the GR case R H 2 = 6π 2 −49 72 .
The position of ILS now read whereas the OLS location is (5.27) Notably, the limit α → 0 yields a finite result that correctly reproduces the light surfaces locations for the Kerr magnetosphere [34] at the third perturbative order in the spin parameter.It is possible to verify by direct substitutions that i 3 and ω 3 as defined in Eq.s (5.23) and (5.24), satisfy the stream equation when evaluated at the locations (5.26) and (5.27).
The solution obtained, hence, is regular at the horizon and at infinity, implying regularity at the ILS and OLS as well.

Spinning split-monopole in Kerr-MOG and bunching of field lines
We conclude this section by presenting plots of the analytic solution we perturbatively derived for a split-monopole magnetosphere in a slowly rotating Kerr-MOG background up to third order in the spin parameter ϵ.
It appears evident from the left panel of Fig. 3 that the magnetic field lines are smooth at both the ILS and OLS, whose respective locations are given analytically in Eq.s (5.26) and (5.27).It is possible to observe that, because of the enhanced gravitational attraction of a Kerr-MOG black hole, the light surfaces lie closer to the event horizon compared to the case of the magnetosphere for a GR Kerr black hole of the same ADM mass.This is displayed in the right panel of Fig. 3, where the fractional deviations of the positions of the light surfaces with respect to the GR case is plotted as a function of the polar angle.
In Fig. 4 the ratio between the angular velocity of the magnetic field lines Ω(Ψ) and the angular velocity of the black hole Ω H is plotted with respect to the polar angle θ.The comparison with the GR case (α = 0) makes it clear that a monopolar magnetosphere around a Kerr-MOG black hole spins faster when compared to the Kerr case.
Numerical [42] as well as analytical [72] studies of magnetospheres in Kerr spacetime lead to the observation that the magnetic field lines tend to bunch up towards the rotational axis θ = 0 when the black hole is in the high-spinning regime.With the magnetospheric solution derived here, we are now able to investigate whether this bunching of field lines also occur around a Kerr-MOG black hole.
To this aim one can compute the contravariant component of the radial magnetic field which, according to Eq. (3.4), is related to the magnetic flux through B r = ∂ θ Ψ/ √ −g.By converting the expansion in the spin parameter ϵ into an expansion for a dimensionless angular velocity at the horizon ω H = M α Ω H , the second-order accurate expression for B r reads which directly reduces to the expression known in the Kerr metric when α → 0 [42]., plotted together with the magnetospheric surfaces of interest: red for the ILS, purple for the OLS, gray for the ergosphere and black for the event horizon.The magnetic field lines correspond to curves of constant Ψ, depicted in blue for ϵ = 0.9 and as dotted lines in the static case ϵ = 0.The plot has been obtained by fixing α = 0.23.Right Panel: plot for the fractional deviation ∆ X = (r X (α) − r X (0))/r X (0) of the ILS and OLS positions in the case ϵ = 0.9 and α = 0.23, with respect to the GR case α = 0.The negative values of the fractional deviations for both the ILS and the OLS indicate that the critical surfaces are closer to the black hole in the Kerr-MOG case.   .The colours represent different spin parameters, whereas dot-dashed and solid lines distinguish between the GR case and the maximal value of α consistent with (2.14).The increasing of the value for B r at θ = 0 is considered as a signature of the bunching of field lines.
In Fig. 5 we present a plot for B r evaluated at the horizon r = r + , obtained by varying the MOG deformation parameter α and the angular velocity of the Kerr-MOG black hole ω H .In particular, we present a comparison between the GR case, α = 0, and the case in which α takes its maximum value allowed for a given black hole spin, i.e. α * = ϵ −2 − 1, (see Eq. (2.14)).The former case is depicted in the plot with dot-dashed curves, for which ω H = ϵ/(2 + 2 √ 1 − ϵ 2 ), whereas for α = α * solid curves are adopted and one has ω * H = ϵ/(1 + ϵ 2 ).All possible intermediate values for α, thus, lie within dot-dashed and solid curves of the same colour in Fig. 5.
The figure makes immediate to observe that increasing the angular velocity of the Kerr-MOG black hole and the MOG deformation parameter leads to two competitive effects.For fixed black hole spin, indeed, B r H at θ = 0 increases as the deformation parameter α approaches the maximum value.In other words, MOG deviations from GR contribute in a positive manner to the bunching of the field lines towards the rotational axis of the black hole.As one spins the black hole up, the maximum value for B r H at θ = 0 is reached for ω ≈ 0.42, namely ϵ ≈ 0.54.After this, the value of B r H decreases until it reaches ω H → 1/2 (corresponding to α = α * and ϵ → 1), where the window of allowed values for α narrows down to α → 0. In this limit, clearly, the value of B r H attains the corresponding value for an extreme Kerr black hole.

Blandford-Znajek mechanism in Kerr-MOG
We assume that the jet is black-hole powered, and that most of its energy is extracted by means of the BZ mechanism [6].The perturbative solution we derived for the magnetosphere can be therefore exploited to compute the rate of energy and angular momentum extracted from the Kerr-MOG black hole.
In particular, the power and the angular momentum per unit of time extracted at the horizon can be computed by means of the following integrals [34,52] Ė(r + ) = 2π π 0 Ω(r + , θ)I(r + , θ)∂ θ Ψ(r + , θ)dθ , By making use of the expansions of the field variables in the spin parameter ϵ one can write the expressions above as a series expansion in ϵ as follows where inside the parenthesis we have explicitly written the dependence on the coefficients of the expansion for Ψ, I and Ω that contribute to the integrals in Eq. ( 6.1) at the specific order in the ϵ expansion.By means of Eq. (2.12) it is possible to trade the expansion in the dimensionless spin parameter ϵ for an expansion in the angular velocity of the black hole Ω H , and arrive at the more familiar expressions for the BZ rate of energy and angular momentum extraction at the horizon of a Kerr-MOG black hole As mentioned in the introduction, the prefactors are related to the characteristic geometrical quantity κ = 2π 3 • 1 4π 2 ≈ 0.053 of a monopolar magnetosphere (that, as already inferred in sec.4.3, remains unaltered with respect to the GR case), whereas the deviation functions and the explicit expressions for C E,L (α) are given by We recall that R H 2 (α) is given in Eq. (5.15) as a function of α, which in the GR case reduces to R H 2 = 6π 2 −49 72 .In other words, the expression for the deviation function f E α (Ω H ) above is a highly non-linear function of the deformation parameter α which, in the limit α → 0, correctly reproduces H ), as previously obtained for Kerr black holes [34,42].An analogous argument holds for the function f L α (Ω H ). It is interesting to notice that, at this order in perturbation theory, the deviation functions for the energy and angular momentum extraction rates share the same dependence on the deformation parameter α, and differ only by means of numerical factors.
It is immediate to recognise that the leading term in Eq. (6.3) is precisely the same that one would have obtained for the BZ mechanism derived in [6] for the magnetosphere -25 - 3 2+α , displayed in the legend together with the corresponding value of α.The dotted line stands for the curve of all the maximum values for the power extracted from a Kerr-MOG black hole at extremality.The region between the dashed and the dotted curve is the region accessible for the MOG case, according to the physical constraints (2.14).More specifically, the green and red areas respectively correspond to estimated range of values for α in the case of supermassive [74,75] and stellar mass [76] MOG black holes.
surrounding a Kerr black hole.Such a degeneracy was first noted in [37] 3 to affect the leading order term of the power and angular momentum extracted for the BZ mechanism in alternative theories of gravity.Effects of modified gravity and of specific configurations for a black hole magnetospheres can therefore be isolated and disentangled only by considering subleading contributions in the BZ mechanism, that enters at order Ω 2 H in the factor f E α (Ω H ).
The plot in Fig. 6 displays the function f E α (Ω H ) according to Eq. (6.4) as a function of the black hole angular velocity Ω H and for specific values of the MOG parameter α.In particular, the values in the parameter space have been chosen to be consistent with previous estimates for the deformation parameter that can be found in the literature [77].For stellar mass Kerr-MOG black holes [76] derived an upper limit α < 0.1 (light red zone in Fig. 6).For supermassive Kerr-MOG black holes the deformation parameter lies in the range α ∈ [0.03, 2.47] (light green zone in Fig. 6), with the upper limit obtained in [74] to reproduce the rotational curves of white dwarf galaxies, and the lower limit derived in [75] to study globular cluster velocity dispersion.As it is clear from the picture if we consider a Kerr and a Kerr-MOG black hole of the same ADM mass and angular velocity Ω H , the power extracted at the horizon is reduced in the Kerr-MOG case compared to the result for a Kerr black hole.
H ). The deviations from the GR factor f E 0 are of order ≤ −5% in the case of stellar mass black hole, whereas they can attain ≤ −14% in the case of supermasssive black holes.
In Fig. 7 we plot the relative deviation of f E α with respect to its GR limit f E 0 .Interestingly, the relative deviations become more relevant in the region of the parameter space which corresponds to MOG black holes in the supermassive regime, which precisely constitute the primary candidates for EHT observations.Our results in Eq.s (6.3), (6.4) and (6.5) (illustrated in Fig.s 6 and 7) show that by combining high-precision estimates of the jet power with independent measurements of the black hole spin or angular frequency, it is possible to probe the metric of astrophysical black holes and possibly put constraints on deformation parameters [35].In the present work, specifically, we focused on the Kerr-MOG scenario, and obtained a non-degenerate expression for the BZ power emitted at order O(Ω 4 H ) without making assumptions on the magnitude of the deformation parameter α.This constitutes an advancement with respect to the current literature about the BZ mechanism in alternative theories of gravity, which either truncated the expression for the power emitted at the leading order [36,73] or exploited a double expansion in both small spin and small deformation parameters to derive next-to-leading order results [37,38].

Concluding remarks
The main purpose of this work is to study the BZ mechanism around Kerr-MOG black holes [13].In order to accomplish this goal several intermediate results have been achieved.
More specifically, in Sec. 4 we analytically classified all vacuum static magnetic field configurations around non-spinning black holes in MOG in terms of angular harmonics and radial Heun's polynomials [61].It is important to stress that these results are solely based on the singularity structure of the Laplacian operator in metrics akin to the Reissner-Nordström metric.We therefore envision that the solutions here derived can also be useful in studying magnetic field configurations around electrically charged black holes in the test field limit [78].We explicitly showed that, while the solution for a static monopolar vacuum field in MOG is indistinguishable from its GR counterpart, the case with vertical asymptotic topology is qualitatively different when compared to the case of a Schwarzschild black hole in the stronggravity region.We expect that this difference can reflect in the geometrical factor κ present in the expression for the energy extracted in the BZ mechanism, though further investigations are needed in this direction.
In Sec. 5 we considered the BZ perturbative approach [6,30,34] in order to construct the first analytical model for a spinning monopolar magnetosphere in a Kerr-MOG background, up to the third order in the perturbative expansion.At each order in perturbation theory we proved the smoothness of the solution across all the critical surfaces characterising the magnetospheric problem, and we studied how the presence of a MOG deformation parameter contributes in a positive manner to the bunching of the field lines towards the rotational axis of the black hole.
Having an analytical description of black hole magnetospheres is important and interesting in its own right.First of all since our understanding of the energy extraction mechanism is incomplete, and only through an analytical model one can attain a deeper understanding.Moreover, the analytical model are complementary to the numerical simulations.In the analytical model one can directly obtain the dependence on the key parameters such as the black hole angular velocity, whereas the simulations can only cover one set of parameters at a time.In addition, the analytic solution derived here can be beneficial to adapt GRMHD codes which exploits force-free approximation and stationarity, and to perform numerical simulations of black hole magnetospheres in the Kerr-MOG background.
Finally, in Sec.6 the explicit expression for the power extracted at the horizon of a Kerr-MOG black hole in the BZ mechanism was computed.We showed that its expression, Ė = 2π 3 Ω 2 H f α (Ω H ), is formally similar to the one obtained in the Kerr black hole background.In fact, in the case of monopolar magnetospheres, only the function f α (Ω H ), that accounts for deviations from a quadratic dependence on the angular velocity Ω H , allows to distinguish the MOG case from the standard GR case.As an important result, we showed with an explicit example that the expression for f (Ω H ) depends on the specific theory of gravity on which the BZ mechanism is set to operate.For the Kerr-MOG case, we derived f α (Ω H ) up to orders O(Ω 4 H ), as explicitly given in Eq. (6.4), and the subregion of the parameter space within which the function can vary is depicted in Fig. 6.In Fig. 7 we also showed that the fractional deviation of f α (Ω H ) from the expression it takes for the standard Kerr case is more relevant in the range of the MOG parameter α which characterises supermassive black holes.
As already emphasized, the analytical approach is relevant in order to obtain a clear understanding of the physics of black hole magnetospheres.Moreover, our model and the expression we obtained for f α (Ω H ) can provide analytic support for the construction of novel GRMHD simulations that take into account the MOG deformation parameter, and which can be used to constrain future high-precision horizon-scale observations from the EHT collaboration.In the context of GR it is known that in order to reproduce the numerical data for the power emitted in the BZ mechanism for black holes in the high spin regime further subleading corrections in Ω H are needed in the expression of f (Ω H ) [34,42].Given that GRMHD simulations in the high-spinning regime for black holes are computationally expensive [45], and that a complete knowledge of the BZ mechanism in modified theories of gravity would require an entire scanning of the parameter space, enhanced by the presence of one or more deformation parameters, analytic models as the one proposed here and higher-orders extensions are expected to provide precious information to overcome these issues.We leave the construction of additional subleading corrections for future works.Finally, while this research focused on the specific case of the MOG scenario, it would be extremely interesting to extend the analysis of the BZ mechanism to theory-agnostic backgrounds, such as the Konoplya-Rezzolla-Zhidenko metric [79].A first step in this direction was taken in [73], even though the analysis was limited to the leading order contribution for the power emitted, which cannot be used to distinguish GR from alternative theories of gravity due to a degeneracy among the spin and deformation parameters, see [37] and our discussion in Sec. 6.We plan to investigate theory-agnostic backgrounds in future projects.The model constructed here should be considered synergetic to future theory-agnostic studies which can use our results to make comparison with the specific MOG case.

Figure 1 .
Figure 1.Asymptotically vertical magnetostatic field lines (blue) around a static MOG black hole.For comparison in the same plot the case of a vertically uniform field (dashed black) around a Schwarzschild black hole in GR -whose horizon is depicted in gray -is reported.The plot assumes α = 10 to magnify the distorsion of the magnetic field lines in the vicinity of the non-spinning MOG black hole.

Figure 2 .
Figure 2. Plot for the function R 2 (w), obtained by varying the MOG parameter α.The dashed line represents the function R 2 (w) in the GR case, Eq. (5.13).The dotted vertical line marks the position of the static event horizon at w = 1, namely r = r+ .

Figure 3 .
Figure 3. Left Panel: the magnetic flux Ψ for a monopolar configuration around a Kerr-MOG black hole, plotted together with the magnetospheric surfaces of interest: red for the ILS, purple for the OLS, gray for the ergosphere and black for the event horizon.The magnetic field lines correspond to curves of constant Ψ, depicted in blue for ϵ = 0.9 and as dotted lines in the static case ϵ = 0.The plot has been obtained by fixing α = 0.23.Right Panel: plot for the fractional deviation ∆ X = (r X (α) − r X (0))/r X (0) of the ILS and OLS positions in the case ϵ = 0.9 and α = 0.23, with respect to the GR case α = 0.The negative values of the fractional deviations for both the ILS and the OLS indicate that the critical surfaces are closer to the black hole in the Kerr-MOG case.

Figure 4 .
Figure 4. Angular distribution for the velocity of the magnetic field lines Ω(Ψ) in the monopolar case, for four different values of the MOG parameter, α = 0, 0.1, 0.15, 0.23, and for ϵ = 0.9.Notice that the value α ≈ 0.23 approximately corresponds to the maximal value of the MOG parameter when the spin of the Kerr-MOG black hole is fixed to ϵ = 0.9.Viceversa α = 0 corresponds to a Kerr black hole.The plot has been obtained by converting the expansion in the spin parameter ϵ into an expansion in the black hole angular velocity Ω H , Ω = ω 1 Ω H + ω 3 Ω 3 H + O(Ω 4 H ).

Figure 5
Figure5.The colours represent different spin parameters, whereas dot-dashed and solid lines distinguish between the GR case and the maximal value of α consistent with(2.14).The increasing of the value for B r at θ = 0 is considered as a signature of the bunching of field lines.

Figure 6 .
Figure 6.Plot for the function f α , given in Eq. (6.4) and characterising the BZ mechanism in the MOG background, as a function of the dimensionless horizon angular velocity ω H = M α Ω H .The curves are obtained by varying the MOG parameter α for specific reference values.In particular, the dashed black line corresponds to the GR case, α = 0.Because of the physical bound (2.14), the curves in the MOG case truncate at the maximum spin value ϵ * = 1/ √ 1 + α, ω * =