Pair Plasma Cascade in Rotating Black Hole Magnetospheres with a Split Monopole Flux Model

An electron-positron cascade in the magnetospheres of Kerr black holes (BHs) is a fundamental ingredient to fueling the relativistic γ-ray jets seen at the polar regions of galactic supermassive BHs (SMBHs). This leptonic cascade occurs in the spark gap region of a BH magnetosphere where the unscreened electric field parallel to the magnetic field is present; hence, it is affected by the magnetic field structure. A previous study explored the case of a thin accretion disk, representative of active galactic nuclei. Here we explore the case of a quasi-spherical gas distribution, as is expected to be present around the SMBH Sgr A* in the center of our Milky Way galaxy, for example. The properties and efficiency of the leptonic cascade are studied. The findings of our study and the implications for SMBH systems in various spectral and accretion states are discussed. The relationships and scalings derived from varying the mass of the BH and background photon spectra are further used to analyze the leptonic cascade process to power jets seen in astronomical observations. In particular, one finds the efficiency of the cascade in a quasi-spherical gas distribution peaks at the jet axis. Observationally, this should lead to a more prominent jet core, in contrast to the thin disk accretion case, where it peaks around the jet–disk interface. One also finds the spectrum of the background photons plays a key role. The cascade efficiency is maximum for a spectral index of 2, while harder and softer spectra lead to a less efficient cascade.


INTRODUCTION
The environment surrounding rotating supermassive black holes (SMBH) are not only fully general relativistic but can contain an extreme astrophysical plasma environment.These systems are surrounded by a plasmarich magnetosphere and emit highly energetic γ-ray jets emanating from near the horizon.Blandford and Znajek (BZ) (Blandford & Znajek 1977) were the first to theorize how the magnetospheric plasma could power the mcsitarz@ku.edumedvedev@ku.edualexlford@ku.edu* Adjunct Research Support γ-ray jets by converting the rotational energy to electromagnetic energy via the Poynting vector ( ⃗ S).This process can be expressed using the spin down luminosity (Koide et al. 2002) where Ω F denotes the angular frequency of the magnetic field lines rotating about the BH, ω H = ac/r s r H defines the angular frequency, r s = 2GM c −2 is the Schwarzschild radius, r H is the horizon radius r H = GM/c 2 + [(r s /2) 2 − a 2 ] 1/2 , a = J/M c is the spin parameter of a BH, and B ⊥ is the perpendicular component of the magnetic field with respect to the BH surface.(Blandford 1993), (Blandford & Levinson 1995), (Levinson & Blandford 1996), and others (Beskin et al. 1992;Hirotani & Okamoto 1998;Levinson & Globus 2016;Globus & Levinson 2013) further proposed a mechanism to fill the magnetosphere of a rotating, uncharged black hole with plasma that could then undergo the BZ mechanism.A similar mechanism works in a pulsar magnetosphere, see (Michel 1975;Cheng et al. 1986a,b;Beskin et al. 1986;Sturner et al. 1995;Lovelace et al. 2006;Sulkanen & Lovelace 1990;Higgins & Henriksen 1997, 1998;Sturrock 1971) and the references therein.These mechanisms for powering jets often assume that the magnetospheric plasma is force free: where ρ e , ⃗ E, ⃗ j, c, and ⃗ B are the charge density, electric field, current density, speed of light, and magnetic field (respectively).For a force free plasma, the magnetic pressure (P B = B 2 /2µ 0 ) must also greatly exceed the plasma pressure (P P = nk B T ).This allows nonmagnetic forces to be neglected.The region where e ± plasma can be produced is called the "spark gap" or just the "gap." If the conditions are right in the gap, an e ± plasma cascade occurs, filling the force free magnetosphere with plasma which in turn powers the γ-ray jets through the BZ mechanism.This cascade exists in a cyclical lifespan described in the following paragraph.
Assume the Kerr BH (KBH) is effectively in a plasma with a magnetic field that threads the horizon.Also assume that there is a bath of soft photons originating from the accretion disk surrounding the KBH.The KBH now can act as a conductor, the spin of the BH creates a motional electric field with frame dragging effects.If the component of the electric field parallel to the magnetic field is non-zero then particle acceleration can occur.Electrons and positrons accelerated by this parallel electric field have a chance to Compton up-scatter the soft photons existing around a BH (produced by an accretion flow, for example).The up-scattered γ-ray photons have a chance to interact with each other and produce and electron-positron pair: provided their total energy in the center-of-mass frame exceeds 2m 2 c 2 .With two new charged particles, this process of e ± acceleration, Compton up-scattering, and γγ pair production repeats.
This process has two threshold parameters that can dictate the processes efficiency and productivity.The size of the region of plasma (referred to now as the gap) must be large enough.The larger the gap, the more region is available for a particle to accelerate and scatter, the higher the probability of scattering.Second, the photons must have enough energy to exceed the rest mass of the secondary particles.If the mass-energy threshold is not met, the pair production will not occur.This work will extend the study of (Ford et al. 2018) by sweeping the parameter space of BH mass and spectral index using a different magnetic flux function than that of previous works: which represents the force-free "split-monopole" configuration.This new model more accurately reenacts what occurs in nature which can further our understanding of various systems seen in observations without the use of assumptions or limiting the types of environments covered by the previous flux model.Section 2 will introduce the mathematical framework of the study is conducted in along with the cascade mechanism expressions and system of equations that are solved along with the boundary conditions imposed on the system (Hirotani & Okamoto 1998;Hirotani & Pu 2016).§3 details the techniques involved in solving the rigid system of equations that govern the BL processes.The data of the study is presented in §4, followed by a results discussion and conclusion in §5.

Metric and Coordinate System
Given a KBH with mass M , angular momentum J, and charge Q, we assume the force free magnetosphere surrounding the surface is stationary and asymmetric.To study this system, we employ the Kerr metric (Kerr 1963) expressed through the Boyer-Lindquist coordinate system (Boyer & Lindquist 1967).This system of coordinates describes space-time in the vicinity of a BH.The coordinate transform for an oblate spheroid is as follows: with the time dimension still expressed by t.
In the Kerr Black Hole scenario, charge Q = 0.This constraint simplifies the Boyer-Lindquist coordinates while also allowing for a dimensionless spin parameter, a, to be defined.Black hole characteristics angular momentum J and charge Q are constrained by mass M , which may take any positive value.The constraint equation simplifies to This can be rearranged to express limits on the value In natural units (c = G = 1), J/M is labeled the Kerr Parameter.This work will employ a unit system in which the spin parameter a ≡ J/M c (found by multiplying the original expression by GM c −2 ).With a defined, the horizon radius can be re-expressed using the gravitation radius r g = GM/c 2 and spin parameter, . With these parameters, the Kerr metric under the Boyer-Lindquist coordinates system can be properly expressed using coordinates (t, r, θ, ψ).
The Kerr metric can be succinctly expressed via tensor notation ds 2 ≡ g tt dt 2 +2g tψ dtdψ+g ψψ dψ 2 +g rr dr 2 +g θθ dθ 2 .(9) Here we can note that the term dtdψ implies a timespace coupling in the plane of rotation as long as J > 0 holds.The field tensors can be defined in terms of length scales ∆ (called the discriminant), Σ, A, r H , and a: Using the above expressions, the line element ds 2 can be fully expressed in terms of the Boyer-Lindquist coordinates and BH base characteristics (shown in full in Appendix A).From the full expression, it can be condensed using the spin parameter, radius r s , length scales ρ 2 , ∆, Σ 2 , and scalar functions introduced by (Price & S. 1986) α, ω, and ω.These new condensing expression are defined as: The line element can now be expressed in its final form

Black Hole Magnetosphere Field Equations
To describe the electromagnetic fields and phenomena that occur in the magnetosphere of a KBH, we apply the 3 + 1 split rule of electrodynamics laws (Price & S. 1986); where 4D space-time is split into global time t and 3D curved space.
Before moving on, the fields inside the gap must be defined according to the force free condition on the gap.We recall the force free condition for the magnetosphere (ρ e ⃗ E + j/c × ⃗ B = 0) where we now define as the charge density in the plasma contributed by the positrons and electrons.This then recovers the starting Poisson equation For a force free magnetosphere that follows the degenerate condition ⃗ E • ⃗ B = 0, see (Macdonald 1984), the toroidal magnetic field can be written as where I(Ψ) is the current flowing back towards the BH surface with constant flux Ψ(r, θ) (the actual form of the flux function is discussed below).The poloidal field is defined as.
From this point, we use the notation for fields ⃗ F T = ⃗ F ⊥ and ⃗ F P = ⃗ F ∥ , interchangeably in relation to the surface of the KBH.To find the electric field components, we recall the gap has a frozen-in condition.The electric field can either be calculated through or more simply Here we define the rotational velocity of the field lines measured by a ZAMO (zero angular momentum observer) by where the angular velocity of the field lines Ω F = dϕ/dt and angular velocity of BH ω = (dτ /dt) ZAM O .From these prescriptions, the electric field components recovered are and This gives a Poisson equation of the degenerate, forcefree magnetosphere as This is defined as as the Goldreich-Julian charge density (Goldreich & Julian 1969).Rewriting the electric field as the Poisson equation inside the gap is finally For a tensor notation based discussion of ρ GJ please see (Hirotani & Pu 2016).
The last expression to be defined is the flux Ψ(r, θ).Previous works used a double split monopole defined by This flux function was used to mime the existence of a thin accretion disk in the vicinity of a horizon.In the current work we instead use the split-monopole configuration.This configuration naturally occurs in numerical simulations when the accreting gas is hot.This corresponds to a geometrical thick accretion flow, such as advection or convection dominated accretion flow (ADAF and CDAF).Such accretion flows are generally x-ray bright (in contrast to thin disk accretion), so they are likely present in many AGNs.Splitting the monopole into two hemispherical planes (Michel 1973) can be expressed as As one can see, a discontinuity arises at the equatorial plane (θ = π/2) where a current sheet is present.For this work, θ will be swept from 0 to π/2, avoiding the discontinuity and assuming symmetry in the magnetosphere.
The field configuration derived above is chosen for the following reasons.First, this depicts a simplified model of an accreting BH where the jets are located in the polar regions and the disk is associated with the equatorial plane.Second, inside the gap, the field will always be radial.Therefore, the gap and the plasma production within it is nicely captured by the above configuration.
While the gap in the magnetosphere has been mentioned and the fields inside defined, the null surface has yet to be fully explained.There is a surface that exists in the magnetosphere where ρ GJ = 0.This null area has very strong ⃗ E || parallel to the magnetic field lines.There exists a surface of charge deficit around this electric field.As described above, an e ± cascade in this region is required to maintain the force free condition on the gap that will power the BH jets.

e ± Cascade Process
With the global field equations established, we can shift focus to the gap itself and the equations that govern the cascade system.Recall that ⃗ E can no longer be screened out near the gap due to insufficient plasma, giving rise to to the important ⃗ E || field.This phenomena generates the flowing, charge separated plasma seen in the magnetosphere.A schematic of the system is seen in (Fig. 1).While the poloidal magnetic field, ⃗ B p , (Eq.16) may be oblique in general, it can be constrained to the full perpendicular case for this discussion.In this circumstance, the acceleration of charges by the parallel electric field (Eq.21) is at its most efficient and therefore the cascade system itself is at its most efficient.Further limiting the system into one dimension, define the x coordinate to be perpendicular to the "null surface" and zero at the center of the gap.with r 0 being the null surface of the gap.For a spherical gap, with x increasing outwardly, r 0 is the radius at which ρ GJ is zero.With the dimensional reduction, the Poisson equation (Eq.24) now reads With the gap (r 0 ) considerably smaller than H (less than 1% of r H ), the Poisson equation can be adjusted by expanding ρ GJ (x, θ) about x = 0 This recovers the expansion coefficient A θ = ∂ r ρ GJ (x = 0, θ), which is on the order of (Ω F B/2πce).
Inside the gap, E || is accelerating the e ± charges and Compton scattering them off of UV-photons from the accretion disk, making the charges rapidly loose momentum.Though, longitudinal motion is still steady due to E || .From this set up, the equation of motion for a charge can be written using the Lorentz factor Γ, the Thompson scattering cross section σ T , and the energy density of the background radiation field U b ; The Thompson scattering cross section is defined by The second term in the equation of motion (Eq.30) represents Compton drag, which may be overestimated in a general sense due to simplification.This simplification still holds for pair production.At the boundaries of the gap, E || is expected to be zero.Except for at this condition (everywhere else E || > 0), the RHS of (Eq.30) cancels out in the first order where Γ reaches a terminal value The newly created e ± have perpendicular momentum P ⊥ .This momentum is gradually lost as the charges scatter off background radiation photons in a drag length defined by This drag length is smaller than the Compton mean free length1 where the charges scatter off background soft photons to produce the desired γ-ray photons required for pair production.Using these lengths, and the length of the half gap, H, the following relation is defined From this relation, two major approximations are deduced: charges migrate in one direction and the charges have a monoenergetic motion with respect to Γ(x).This relation holds due to X-ray and UV photons contributing to the drag term.Despite the various non-simplified drag terms, the charges can still produce the necessary γ-ray photons with energies up to Γm e c 2 .These γ-ray photons are emitted from inverse Compton scattering (Beskin et al. 1992), not curvature radiation or synchrotron radiation.These new γ-ray photons are free to pair produce charges by collisions with background photons.
For a γ-ray photon of energy ϵ γ m e c 2 to undergo pair production it must collide with a background photon with energy ϵ s m e c 2 .The photons must satisfy where µ is the cosine of the collision angle θ.The minimum energy required for production occurs when the photons collide head on (µ = −1) with the most energetic soft photon collision occurring with ϵ s = ϵ max .Only γ-rays at energies greater than the above minimum energy (ϵ −1 max ) can contribute to pair production.

Basic Physical System Equations
With the physical process of the cascade described, the equations governing the charges and fields can now be examined in the cascade context.To begin, it is necessary and helpful to define the motion of the electrons and positrons in the gap.The motion of the charges is set by the direction of the current as it moves towards the polar regions.For this system, the electrons will migrate inward toward the BH while the positrons migrate outward away from the BH.This set-up is accurate in lower latitudes (small distance away from BH) and appears opposite in high latitudes.But, the treatment of physics holds in each scenario as the signs of the various terms in (Eq.28) resolve themselves and create the same effect.This research will consider low latitudes only.The speed of these mobile charges can be defined as with defined for brevity.The speed is shared by all charges due to the monoenergetic assumption.From this, the continuity equation (Berestetskii et al. 1989) then can be defined as F ± are the Boltzmann equations for the γ-ray photons, defined later in the section (Eq.45).η p (x) is defined as the angle averaged pair production redistribution function expressed via Here, σ p representing the cross section of pair production in collisions between photons of energies m e c 2 ϵ s and m e c 2 ϵ γ moving angle µ with respect to each other; only when the energy relation is met (Eq.35) does the above expression have a non-vanishing value.With the ν being defined as, n ± (x) and F ± (x, ϵ γ ) represent the number density of outwardly/inwardly moving particles and γ-ray photons in the ± x-direction in a non-dimensional energy interval ϵ γ ≈ ϵ γ + dϵ γ , respectively.Finally, dN s /dϵ s is the number density of soft background photons in a non-dimensional energy interval ϵ s ≈ ϵ s + dϵ s .The γ-ray photons created by inverse Compton scattering are highly beamed in the same direction as the charges, e ± , which is to say that the distribution functions can be fully described by F ± .
Current is conserved as it is carried along the field lines.From one combination of the continuity equation (Eq.38) which recovers E || must go to zero at the borders of the gap, which occurs when j 0 = j crit .j 0 is now defined as the current density along a field line, out-flowing from the gap at a constant rate.When j 0 takes the value, within order of magnitude (Price & S. 1986), of j 0 ≈ then energy and angular momentum can be effectively extracted from the KBH.Another combination of the continuity equation recovers In place of the original continuity equation (Eq.38), the two combination of the continuity equation presented above (Eq.43 and Eq.44) will be used.
The single dimensional motion of γ-ray photons obeys a Boltzmann equation of the form Above, η c (x) is the Compton redistribution function with σ KN as the Klein-Nishina (Klein & Nishina 1928) cross section defined by (Lightman & Rybicki 1979) There is the assumption that energy transfer from e ± with Lorentz factor Γ to a photon with incident energy m e c 2 ϵ s is approximately m e c 2 Γ 2 ϵ s in the Compton redistribution function.This holds in its simplified form, but can be more complex if need be.
The migrating e ± 's and γ-rays are described by the differential equations (Eqs. 29,30,44,45) while the number density functions n ± are related by (Eq.43 and Eq.45) (independent of each other).This sums to a system of five differential equations to be solved.
For the system of five equations to be solved, there must first be assumptions about the background radiation field made first.The spectral number density of background radiation per unit interval of ϵ s can be modeled by a power law This C(α) term is a decreasing function with respect to α expressed via with epsilon terms defining the cutoff of the spectrum.U b is the background radiation field's energy density.
2.5.Non-Gray Analysis of the γ-Ray Distribution Examining the expression for η p (Eq. 39), no gray approximation that can be attributed to it to solve the Boltzmann equation (Eq.45) due to its ϵ γ dependence.To rectify this, ϵ γ is split into bins so that η p is approximated in each bin.To do this, let ξ i and ξ i−1 be the upper and lower limits of the i th bin (limits are sufficiently close).Using this, the right hand side of (Eq.44) using a summation instead of solving the following integral types where and In place of integral (44) the following is used (53) with χ number of energy bins.Then (Eq.45) can be integrated over an energy bin between each limit to form where The system of equations now composes of (Eqs.29, 30, 53, 54).

Gap Boundary Conditions
To fully study the physical processes and environment within the gap, it is helpful to consider the case(s) when the functions of E || , Γ, n + and f ± i have symmetric conditions.These boundary conditions are defined in a way so symmetries allow the conditions to be set at the center of the gap (x = H) and at the edge of the gap (x = 0), allowing only integration over half of the gap.This gap width should only be within a few percents of the black hole radius.The gap boundary is defined when the parallel electric field vanishes.Using (Eq.28) at x = H and E || = 0, a smooth curve resulting in dE || dx = 0 at the half gap width can be found with (Eq.43) and To begin, E || should not change sign within the gap and the function should vanish at the boundaries of the gap itself.With H << r 0 (Eq.27), we may assume that E || is an even function with respect to x.This condition also applies to Γ; Next, the assumption of functional symmetry is imposed on n ± (particles) and F ± (γ-ray photons) The consequences of the symmetric functions include essential requirements and features of the pair production cascade.Therefore, (Eqs.29, 30, 53, 54) are solved within 0 ≤ x ≤ H.
With these functional conditions imposed the derivation of the full boundary conditions at x = 0 and x = H can begin.First, we derive the conditions for the inner boundary.From (Eqs.29, 30, 57) which is equivalent to is found for x = 0. From (Eqs.43 and 58), x = 0 also recovers 2n Then, following the second symmetric function, (Eqs.52 and 59) for each value of ϵ γ .
The boundary at x = H is formulated to be an free boundary that ensures E || smoothly vanishes as it approaches H.To accomplish this, any inwardly propagating particle -in this scenario the e − 's -should not enter the gap from x > H. Combining this condition with (Eq.43) at x = H results with Additionally, the charge density distribution must remain continuous at the outer boundary through (Eqs.29 and 64) to recover 1 4π Similar to the inner boundary, no γ-ray photons may come into the gap from the outside (all up-scattered photons must be created inside the gap), represented mathematically via From the above derivations, there is now (2m + 5) total boundary conditions for (2m + 3) total unknown functions for E || (x), Γ(x) (which is contained in the β(x) function), n + , and f ± i (x) for (i = 1, 2, ..., m) and two uncalculated constants H and j 0 .These are formed from the conditions E(H) = 0 and E ′ x (H) = 0, so j 0 now plays the role as an eigenvalue of the boundary value problem.
It is sufficient for the investigation into the pair production to consider γ-ray photons satisfying which is a dimensionless energy parameter.Below the minimum energy threshold of ϵ min = ϵ −1 max , photons fail to contribute to pair production.The lowest energy bin may be defined as β 0 m e c 2 with β 0 = ϵ −1 max .These may be chosen based on the study at hand and will be numerically defined in the following section.The underlying goal is then to seek solutions to the boundary conditions by way of the shooting method of solving boundary value problems.This method involves taking the boundary value problem and reducing it to an initial value problem at different conditions until the solution that satisfies the original problem is met.For this study, the shooting algorithm begins at x = 0 and finishes at the outer boundary x = H.
For a realistic model of a BH, the current flowing along each field line is determined more by the global system physics rather than the micro-scale physics inside the gap.(Price & S. 1986) shows that the load connecting wind/jet to the unipolar induction is the reasoning behind this.We define the "null surface" to be where the gap center is located and where ρ GJ vanishes.This definition imposes a symmetrical center point in the gap, and allows j 0 to be treated as an eigenvalue in a standard boundary value problem (facilitating the shooting method) and focusing on finding j 0 with the outer bound of x = +H.Following this, the shooting method scheme can be listed as follows: • Very small values of j 0 (with |n + − n − | is also very small), (Eq.29) can be approximated with • The above quadratic form function cannot solve 1 4π • As j 0 increases, the RHS will also increase (monotonically) with x in (Eq.29) • As j 0 increases more, continues to decrease until it vanishes at some j cr (U b , α).
• Above ȷ cr , the second condition can not be satisfied, no matter what initial conditions are given.
To summarize thus far: the system seeks a solution that satisfies the second condition by adjusting j 0 to j cr (U b , α).This resulting system is then molded by 29 different equations for 31 different unknown variables which are then integrated under 31 boundary conditions.A table summarizing these boundary conditions can be seen in Appendix C.
To study the e ± cascade seen in the magnetospheres of KBHs, a C ++ code (Hirotani & Okamoto 1998;Ford et al. 2018) is used to solve the system of equations for a given set of initial conditions.To conduct this analysis, a set of initial conditions was set and solved for.This control initial condition set can be found in Appendix B8.From here, a single parameter is incrementally changed ("swept") and the new system is then solved.This generates a new set of initial conditions the code then solves for and the process repeats; set, change, solve, set.This incremental change must be small enough that the code can start its process using the previous set of initial conditions and the new parameter.
This code relies on the shooting method of solving a boundary value problem (BVP) using a set of initial parameters defined at the top of the code (gap width, Lorentz factors, energy bins, etc.).The shooting method takes the multi-dimensional BVP and simplifies it to a single initial value problem (IVP).From the IVP, the code takes the solution to the equations at one boundary (x = 0) and attempts to find the solution that solves the equations at the other bound (x = H) by "shooting" different solution guesses to the boundary.Once the solution is found, namely j 0 , the remaining parameters are solved for and the parameter under analysis is incrementally changed, and the code solves the system again.If the change in the swept parameter is too big, the guess of solution will be too far off to solve, and the system of equations will collapse.
While this code accurately solves the mathematics within the gap, the scope of the code stops there.Despite the importance of the solutions found in this gap, it is also important to visualize how this plasma fuels the jets at the poles.Various studies, including ones done by (Crinquand et al. 2020;Ptitsyna & Neronov 2016;Kisaka et al. 2020;Chen et al. 2018) use PIC codes to model the system globally, from soft photon bath to jets.These studies are just as important as they allow us to study the state of plasma and gap and view their effects on the jets themselves; which we view through observational astronomy.

ANALYSIS OF CHANGING FLUX MODEL
One of the main goals of this study is to alter the flux model used from the simplified model seen in previous studies (Ford et al. 2018) to a more realistic model motivated by simulations (Crinquand et al. 2020).The old flux expression, sin(θ) 2 , represents a toy model of accretion.Caveat is "pathology" at large distances, represented by a diagonal line at large r.The updated model, (1 − cos(θ)), has no caveat present.Physically, this model represents sites with no accretion (MW type systems), accretion from a thick, hot disk (AGN systems), and far away, thin disks.Pictorially, we can see the null surface reaching outward from the surface of the KBH for the toy model.Mathematically, the difference between the flux models θ dependence of A θ ≡ ∂ r ρ GJ .This difference stems from the expansions coefficient's (A) distribution of field around the null surface, x = 0 (within first order).For the updated model, the null surfaces stretch at equal latitudes around the BH.This new flux model more accurately represents the magnetic field structure -split monopole.Simulations show the split monopole configuration forms near KBHs, so it is a natural force free configuration.A table of the initial conditions (ICs) used in this paper are found in Appendix B.
The differences in the flux model can be examined through global visualizations of the solutions to many of the equations presented in the section above.Looking at distinct theoretical solutions for ρ GJ , stagnation point (flow of in flowing or out flowing material is zero), light cylinders (regions of warped spacetime so severs that there is no sub-luminal path out), and ergosphere (region in which a observer must rotate with the KBH and can no longer be a stationary observer), initial differences in the models can be seen, but are not complete.Building on these solutions, we can view more differences in the flux models by the plasma density ρ GJ (Fig. 2) and the parallel electric field (Fig. 3).These plots more greatly illustrate the major differences in field and plasma density distribution.We see two very distinct model structures in both quantities; in structure, magnitude, and distribution.The effects these differences have on the structure of the gap will be examined in the following subsections.

Flux Transition
To study this model, the code must first be set up for the new flux conditions.As stated above, the code must sweep a parameter from the control initial conditions to the desired point, as opposed to plugging in the new model directly.To achieve this, the flux models are written in a summation form with each other as a function of percentage, with η representing the percent of the old model and ι representing the percent of the new model: where we define ι = η-1 to restrict the free variables to only one.The code begins with η = 1 and incrementally decreases until η = 0, where the new flux model is fully implemented.This change of flux is done at the axis of rotation of the KBH θ = 0.
Plotting the three main gap characteristics viewed in this study (half gap width, Lorentz factor at gap center, and outgoing particle flux), we can examine the transition of fluxes at the polar cap of the KBH.As a note: the plots shown in this paper that are labeled "normalized" are normalized to the initial value in that specific parameter sweep.This allows a progression trend to be seen with respect to the starting point.
We can immediately see in (Fig. 4) the size of the gap gradually increases as the flux shifts to the new, realistic model.With the increase in gap size, there is also a decrease of Lorentz factor within the center of the gap.This can be from the expanded space in which particles decrease motion and loose energy.The most drastic change in trend is shown in the normalized outgoing particle flux, where we see a trend that leads to almost 100% decrease in flux at the polar cap.The loss of energy in particles from the Lorentz trend hints that they can no longer scatter of photons, decreasing the flux.While these trends can show general behavior or trend, they are not indicative of the while system, as they are only at θ = 0.
Viewing the comparison of the old model to the new model (Fig. 5) in the same KBH system (that is only the flux has changed from the control ICs) gives a more complete view of what should be expected with the split monopole model.Immediately, it can be seen that all three plots show an oppositely directed trend, which is what is expected in this switch.It can also be noted the lack of severity in the trends when compared to the sin(θ) 2 model.This depicts a more evenly structured gap and a less drastic change in emission based on viewing angle to the axis of rotation.
The observed opposite trends for the quasi-spherical gas distribution studied here and the thin accretion disk studied in a previous paper (Ford et al. 2018) indicates the importance of the magnetic field geometry.The former (quasi-spherical) case show the highest cascade efficiency at high latitudes, toward the jet axis.The latter, in contrast, indicates highest efficiency around the jet-disk interface.This can be a direct prediction of our study and a possible observational diagnostic of the field geometry at astrophysical sources.

Varying Mass
To fully investigate the differences in the flux models, direct comparisons from parameters studied in the old flux model must be made.For this study, KBH mass was varied and then examined as a function of angle θ.For the old flux model, data from (Ford et 2018) was used to study KBH's of 10 6 M ⊙ , 10 7 M ⊙ , and 10 8 M ⊙ .For the updated mass model, 10 7 M ⊙ , 10 8 M ⊙ , and 10 9 M ⊙ was used.To begin, the half gap width of the gap is examined in (Fig. 6).In the unnormalized plots, the higher the mass, the larger the gap is in a general sense.Continuing in the unnormalized trends, not only is there no drastic slope change as the gap approaches the equator, but the direction of the trend itself opposite, with the new model constantly decreasing.Viewing the trend as a ratio from the pole to the equator, we see many differences in the models.Starting with the obvious, instead of a increasing exponential we see a decreasing quasi-linear relation, with plateaus at the x-axis limits.We also see more definite spacing between the masses as opposed to the tight grouping of the sin(θ) model.
The outgoing flux is discussed (Fig. 7) next.We see for both models that the higher the mass of the KBH, the less outgoing flux overall.The updated model sees as small general increase within the mass trends versus the exponential decrease of the toy model, in the unnormalized figures.Under normalization, the old model again displays a bunched exponential tend with the decrease of flux as the angle increases.The new model displays a bunch to spread relationships between the masses, with the lowest mass being between the highest and lowest mass.The trend shows a sin(x 2 ) function, see Appendix D.
Next, the Lorentz factor of the gap is examined in (Fig. 8).Beginning with the unnormalized plots, we see an inverse relationship between increasing mass and decreasing Lorentz factor.For the new model, as the gap reaches the equatorial plane where the gap was smaller, the Lorentz factor slightly increases, as opposed to the exponential decrease seen in the old model.This inverse relationship between Lorentz factors (and by extension γ-ray energies) and gap width can be attributed to the E || field within the gap.The bigger the gap present, the smaller electric field is needed to accelerate particles to the cascade condition.When the cascade begins, it begins to screen the field out.Viewing the normalized trends, the updated model continues to show a more complex trend as well as bunching near the pole of the KBH.These trends show a counter-intuitive cos(cos(x)) dependence that encapsulates the distinct double inflec-tion line see Appendix D.
Following the Lorentz Factor, the value A θ = ∂ r ρ GJ is studied (Fig. 9).This factor shows an inverse dependence on the magnitude of the trend, but the behavior of the trend is not affected by the varying mass.In the normalized plots, the new flux continues to show a polynomial type trend that increases as θ increases around the BH, compared to the decreasing exponential given by the sin 2 (θ) flux model.The unnormalized trends again show a steadily increasing behavior opposite the linear to decreasing trends of the old flux.
The last three figures (Figs. 10,11,12) discussed show trends at discrete polar angles (with the axis of rotation set to θ = 0).These angles are (from closest to axis of rotation to equator): π/58, π/4, π/3, and π/2.The most striking observation is seen in the length of the trends as the mass increases.This is due to effect of Gauss's law at the outer boundary position.
x = H (outer boundary) is determined by the condition that the Goldreich-Julian charge density ρ GJ (x = H) matches the actual charge density ρ, which is given by the current density J divided by c.Here, J is conserved along each flow line because of stationality.Because boundary conditions dictate no particle injections at boundaries, the leptons are composed of either electrons or positrons at the boundary, which means that the real charge density ρ is given by J c at x = H.Therefore, if the gap half width H is measured in gravitational radius, H is solely determined as a function of the dimensionless conserved current density, j = J ρ GJ (x=H)c .Subsequently, if we compare at the same dimensionless current j, the actual gap width (e.g., in cm) increases with increasing BH mass, because the gravitational radius increases accordingly.The second overarching observation is the relative angular symmetry within the three trends.While the previously discussed trends had a more direct relation to the θ dependence of the flux model, the following trends including E || (x), ρ(x), and voltage V (x) appear to show a minimal dependence on θ.Despite a small sample size, the figures show minimal difference in magnitude and overall trend behavior across the entire calculated area for a given plotted value.This is a primary observation, as the entire calculated range (0 − π 2 ) is only represented by 4 polar angles.Finally, the difference in trends between the flux models.For the first three θ values in (Figs. 10,11,12), there are two trend lines per mass, where the square data point represents the old flux model and the triangle data point represents the masses in the new flux model.For each θ value where both fluxes are plot- ted, the difference between the two fluxes for 10 8 M ⊙ and 10 9 M ⊙ is noticeable, but minimal relative to the difference generated by the difference in mass.This includes the previously discussed observations (trend length along the x-axis and angular symmetry) The mass containing the largest differential between fluxes is the lowest mass.However, the trend lines for 10 7 M ⊙ is consistently approximately one magnitude apart.
The first trend (Fig. 10) depicts the plasma density (ρ) as a function of x in cm from the center of the gap at x = 0 to the gap edge at x = H, where H is the half gap width.Along with the observations previously discussed, the trends show a direct relation between mass of the BH and the steepness of the rise in plasma density from center to edge of the gap.This difference in slope can be seen at two points on the trend.As the lines begin, there is a steep increase in plasma density right after x = 0, after which is an area of inflection where the trend levels out and continues the remaining way to the edge of the gap, x = H.It can be seen that the higher the mass of the BH, the lower these two slopes are in relation to lower BH masses.Recall, that the boundary of the gap x = H is defined as the region where plasma density ρ is equal to ρ GJ .The null surface x = 0 is defined as the region ρ GJ = 0.As the charge density (ρ GJ ) increases (rapidly from 0 at first, then gradually as x increases), plasma density rises as well until the two meet.The steepness in trends can be explained by the relative distance the plasma density has before it meets the charge density.With a thinner gap, the density gradient must be larger to cover the same two magnitudes the higher mass systems who posses a smaller gradient over a much larger distance.
Next, the electric field E || (in statvolt/cm) (Fig. 11) as a function of gap coordinate center (in cm) x = 0 to edge x = H is plotted.The electric field peaks, for each mass, just after the null surface and exponentially drops as it approaches the edge.This is a result of the electric field being screened as it approaches the edge of the gap.Also recall, that by definition, the electric field must be zero at the edge.This explains each mass tending to zero as x approaches H.The difference in gradient of the masses is explained by this concept.Because the lower masses have higher starting field magnitude, they must reach zero in a shorter distance.
This can be shown if we recall ) is plotted on the left with the normalized trend above the unnormalized trends.We again see the same behaviors and inter-model relationships as before.This is coupled with the intuitive understanding that as the gap widens, the Lorentz factor will decrease.
and (Hirotani & Pu 2016), where ρ GJ denotes the Goldreich-Julian charge density, Ω denotes the angular frequency of the magnetic field lines, ω does the frame-dragging angular frequency at position x in the gap, B the magnetic field strength at point x, and r g = GM/c 2 the gravitational radius.After some algebra, (Eqs.72 and 73) can be related through Note that dimensionless gap width H/r g cancels out in the RHS.Thus, E ∥ ∝ M −5/4 , which is shown by the curves in (Fig. 11).As E || (x) (cgs) increases, so does Γ(x) for e − 's.Which in turn creates in an increased pair production rate, resulting in more efficient screening of E || (x) near the outer edge at x = H.This results in a smaller gap H for lower BH masses M , this confirming the previous statement.
Finally, the voltage drop in statvolt between the center of the gap (in cm) and a point in the gap is plotted and is plotted on the left with the normalized trend above the unnormalized trends.The unnormalized trends are plotted on the logarithmic scale.We see that for both flux models, the mass is inversely related to the trend's magnitude, but does not affect how the trend behaves.In contrast, the flux function does affect the trend behavior.
studied 12.The voltage in the gap is calculated via where dx is the x-coordinate (the center of the gap is x = 0) interval.Again, some symmetries and relations can be immediately deduced.The relative difference between the old and new flux modes (with the one exception) and the difference between the different angles plotted both continue to be minimal.We see that as E || decreases as x increases in (Fig. 11), the voltage's gradient correspondingly remains positive throughout the gap.The voltage increases with a steepness corresponding to mass (and therefore gap width) with steeper trends being seen at lower mass (thinner gap width).
Looking at all of the trends from the (1 − cos(θ)) model, we examine the relationships between the three main parameters within the model.The trends for A θ and voltage shows some changes based on flux, but mainly the difference can be explained by the previous values or the varying mass.Despite the wildly different functional expressions for the trends, we can derive behaviors in the gap seen across the three masses.We can see a direct relationship between the Lorentz factor and the outgoing flux from the gap.As the ZAMO moves down from the pole to the equator, more radiation is produced.2 is the equator.The square boxes in the legend represent color, while the shapes on the left of the legend denote which flux is being plotted.By definition, the electric field reach zero at the edge.The difference in gradient of the trends is explained by lower masses having a higher starting field magnitude, thus, they must reach zero in a shorter distance.The divergence seen at low mass can be attributed to a numerical issue due to the stiffness of the system of equations.The first challenges to the extended spectral study of e ± cascade in the magnetosphere of black holes comes with changing the spectral number density of ambient soft photons of the system of equations.As expressed in an earlier section, the system of equations that governs this system is rigid, and pushes back any changes to the equations when the system is stable.Viewing the spectral number density as a function of photon energy for different α in (Fig. 13) illuminates the drastic difference of function between integer spectral index numbers.The remaining panels show the spectral energy density as a function of α for the minimum, average, and maximum input spectral energy.The steep curves and sharp inflection is the source of the challenge when it comes to sweeping the index variable.The abrupt change in not ideal for the shooting method and this system of equations, causing the calculations require immense computational resources.
While changing the spectral index has stretched computational resources and the analysis of different spectral states incident on the KBH system is ongoing, the main trends as a function of α at the axis of rotation have been recovered (Fig. 14).While half gap width and Lorentz factor show a manageable change in parameter, the outgoing flux has posed certain numerical challenges.
The change in trend of the half gap width, Lorentz factor, and particle flux as a function of spectral index illustrates the drastic inflections encountered when the code is solving the BVP, shown in (Fig. 14).These inflections are the main challenge when solving for different spectral states, as the code must begin at the known solution (α = 2) to solve for the next spectral state.The different colors of the line signify the various data sets used to create a semi-continuous trend.The code runs for such a extended period of time that a single sweep is not feasible at this time.We emphasize that the α = 2 spectrum yields the most favorable conditions for the pair cascade.This is seen from the minimum Lorentz factor and minimum gap width.Softer (α > 2) and harder spectra (α < 2) reduce the cascade efficiency.

RESULTS AND DISCUSSION
In this paper we explored the role of the magnetic field configuration around a spinning BH on the electronpositron cascade and its properties in particular.Two models, representing a thin accretion disk and a quasispherical gas distribution were compared.We have shown that the new magnetic flux model for a KBH e ± cascade has not only shown the expected oppositely  directed trends, but delivered a more complex mathematical relationship as a function of polar angle.In comparison with the original flux model of sin 2 (θ), (1 − cos(θ)) allows a more realistic treatment of accretion disks far way from the KBH, close disks that are filled with thick, hot material, and systems that have no disk at all such as the Milky Way.This paper details the beginning of an extended study on this new flux model as well as the first computational challenges seen when beginning a study on the spectral states of KBH systems.This study suggests that the parameter trends as a function of θ are not simple monotonic functions, but more complex expressions that have a greater variance withing the differing parameter sets, such as the varying mass shown in this work.The intuitive nature of behaviors in the gap remain; the γ-ray photon flux correlates with the max Lorentz factor and anti-correlates with the gap width.
Particularly interesting results of this study concern the efficiency of the pair cascade.First, the efficiency of the pair cascade in a quasi-spherical gas distribution case studied here increases with latitude and peaks at the jet axis.In contrast, in the thin disk accretion case, highest efficiency shown around the jet-disk interface.This finding is prediction of our study and a possible observational diagnostic of the field geometry at astrophysical sources.Second, the the spectrum of background photons is also a key factor in the leptonic cascade phenomenon.The efficiency is maximum for a spectra index of two.Both harder and softer spectra result is a less efficient leptonic cascade.
Follow up studies are underway examining different parameters and will be readily compared to the results found in (Ford et al. 2018).The physical and spectral parameter space available can help shed light on how these systems utilize the BZ mechanism to fuel powerful general relativistic jets seen in observations.This work was supported by NSF grant PHY-2010109.

Figure 1 .
Figure 1.Schematic depiction of the spark gap and SMBH system.Generally H << r0, unless the cascade is not efficient of present at all.

Figure 2 .
Figure 2. Visual comparison of the charge density environment of the KBH.sin(θ) 2 is shown on the left and (1 − cos(θ)) on the right.In arbitrary units, ρGJ is plotted around the KBH with blue signifying positive density and red negative.The new model shows a less layered distribution of charge with lesser magnitude areas being seen outside the null surface, away from the KBH.The differences in the null surface distribution are seen in the red solid line to follow the less complex distribution of charge as the new model shows a single shell around the KBH while the old model (left) shows a farther reaching trend.

Figure 3 .
Figure 3. Visual comparison of the parallel electric field of the KBH.sin(θ) 2 is shown on the left and (1 − cos(θ)) on the right.In arbitrary units, ⃗ E || is plotted around the KBH with green signifying positive field and red negative.Again, the new model shows a less complex distribution of field.The same contours are plotted for each model, and with this we see the relative bunching of the field in the old model as well as the less circular shapes of the new model's contours as it folds around the event horizon.

Figure 4 .
Figure 4.The transition between the two flux models can be viewed as a function of percentage of the old model.Though these figures only show the axis of rotation, it does give a insight to how the gap will change over all angles.We see the gap widen as the new model takes over, this widening leads to more space for particles to decelerate.This shows up in the Lorentz factor of the particles steadily decreasing.With this decrease in velocity, the outgoing flux follows, as seen in the exponential decrease of outgoing flux.

Figure 5 .
Figure 5.After the new flux model is fully integrated, the two expressions can be directly compared at the axis of rotation (θ = 0).The normalized trends are shown on the right while the unnormalized data is plotted on the left.As expected, the direction of the new model's trend is opposite of that of the old model.The new model also shows less severity as the model sweeps from the pole to the equator.

Figure 6 .
Figure 6.Comparison of gap widths in old and new models (left and right panels, respectively).The old model (sin 2 (θ)) is plotted on the left with the normalized trend above the unnormalized trends.We can immediately see the opposite and more complex nature of the (1 − cos(θ)) model.also note that the old models trends are more bunched in the normalized plots while the new model has variance.

Figure 7 .
Figure 7.Comparison of γ-ray fluxes in old and new models (left and right panels, respectively).The old model (sin 2 (θ))is plotted on the left with the normalized trend above the unnormalized trends.We see the same behaviors and inter-model relationships as before.Interestingly, the total γ-ray flux increases toward the equator, in contrast to the old model.

Figure 8 .
Figure8.Comparison of peak Lorentz factor in old and new models (left and right panels, respectively).The old model (sin 2 (θ)) is plotted on the left with the normalized trend above the unnormalized trends.We again see the same behaviors and inter-model relationships as before.This is coupled with the intuitive understanding that as the gap widens, the Lorentz factor will decrease.

Figure 9 .
Figure 9.Comparison of A θ = ∂rρGJ in old and new models (left and right panels, respectively).The old model (sin 2 (θ))is plotted on the left with the normalized trend above the unnormalized trends.The unnormalized trends are plotted on the logarithmic scale.We see that for both flux models, the mass is inversely related to the trend's magnitude, but does not affect how the trend behaves.In contrast, the flux function does affect the trend behavior.

Figure 10 .
Figure 10.Comparison of plasma density (ρ) as a function of gap coordinate (x in cm) from the center of the gap at x = 0 to the gap edge at x = H, where H is the half gap width, for three different BH masses (line colors) for both flux models (shape of point) when applicable.The panels represent different, discrete angle θ values around the BH where θ = 0 is the axis of rotation and θ = π 2 is the equator.The square boxes in the legend represent color, while the shapes on the left of the legend denote which flux is being plotted.By definition, the electric field reach zero at the edge.The difference in gradient of the trends is explained by lower masses having a higher starting field magnitude, thus, they must reach zero in a shorter distance.The divergence seen at low mass can be attributed to a numerical issue due to the stiffness of the system of equations.

Figure 11 .
Figure 11.Comparison of electric field (E || ) (in statvolt/cm) as a function of gap coordinate (x in cm) for three different BH masses (line colors) for both flux models (shape of point) when applicable.The square boxes in the legend represent color, while the shapes on the left of the legend denote which flux is being plotted.The panels represent different, discrete angle θ values around the BH where θ = 0 is the axis of rotation and θ = π 2 is the equator.

Figure 12 .
Figure 12.Comparison of voltage (V (x) = x 0 E ∥ (x)dx) (in statvolt) as a function of gap coordinate (x in cm) for three different BH masses (line colors) for both flux models (shape of point) when applicable.The panels represent different, discrete angle θ values around the BH where θ = 0 is the axis of rotation and θ = π 2 is the equator.The square boxes in the legend represent color, while the shapes on the left of the legend denote which flux is being plotted.The H/2 voltage peak is another result of the electric field vanishing at the boundaries.With the inverse relation ship between electric field strength and distance, as the Electric field tends to zero, the coordinates continue to rise, giving way to a parabolic fit to the voltage patterns in the gap

Figure 13 .
Figure 13.In the top right, the spectral number density as a function of photon energy is plotted for different levels of spectral index (α) in log-log scale.The functional trend as a function of energy across the differing indices illuminates the challenges of sweeping the system of equations.The remaining three plots show the spectral number density as a function of alpha for the minimum, average (semi-log), and maximum (semi-log) photon energy used in this study.The function's most drastic changes all center around the initial condition of α = 2, creating a challenge to sweep to other spectral indices.

Figure 14 .
Figure 14.Depictions of the three main parameter trends (half width, Lorentz factor, and outgoing flux) as a function of α.The different colors represent different data sets from the simulation as spectral index is swept from α = 2.The sharp inflection seen around α = 2 is believed to be the main challenge.

Table 2 .
(Ford et al. 2018)tions for the SMBH system cascade equations in the gap, recreated from(Ford et al. 2018)