Pair production and annihilation in advective accretion disks around black holes

We investigate the role of pair production and annihilation in the presence of radiation processes like bremsstrahlung and inverse-Comptonization, assuming no magnetic fields present in the system. The disk is assumed to be viscous, advective and in hydrostatic equilibrium. Photon-photon interactions leading to pair production is the most dominant process. Threshold for this process is, hv1hv2>∼(mec2)2 . We find that, for accretion rates lower than 0.8 M edd, we can neglect the presence of any pairs. Maximum number of pair production occurs in a region < 50r g. When the number density of pairs is optimal, a bump is observed near m e c 2 in the spectrum, which is a signature of the annihilated photons.


Introduction
Accretion is the best possible mechanism to explain the most luminous objects in the Universe, the Active Galactic Nuclei (AGN). Works by Hoyle & Littleton (1939) [1] and Bondi (1952) [2] laid the foundation for the theory of accretion. These theories shot into prominence when the enormous luminosities emitted from AGNs and X-ray binaries could only be explained with the help of accretion theory. Since then, accretion models have been improvised to match the observational evidences obtained from the new generation of telescopes. One such model was given by Shapiro, Lightman & Eardley [3]. This was the first accretion model which was composed of hot flows and could match the hard power law part of Cygnus X-1 (the first object to be realised as a black hole candidate). They predicted that the electron temperature in accretion disks to be approximately equal to or greater than the electron rest mass energy, kT e m e c 2 = 511 keV (where, k is the Boltzmann constant, T e is the electron temperature, m e is the mass of the electron and c is the speed of light). It was soon realised that the production of e + − e − pairs as well as its annihilation, is an important process and may play a significant role in accretion disc dynamics as well as in shaping the output spectrum.
In 1966, Jelley [4] pointed out that high energy gamma ray photons generated by luminous astrophysical sources on interaction with low-energy photons present in intergalactic and interstellar medium, could lead to the creation of electron-positron pairs. The reaction reads as : γ + γ → e + + e − . Herterich (1974) [5] showed that photon-photon pair production process has a threshold condition : E 1 E 2 (1 − cos θ pp ) ≥ 2(m e c 2 ) 2 , where E 1 = hν 1 and E 2 = hν 2 are the energies of the interacting photons and θ pp is the angle between them. He concluded that this process could significantly change the radiative output of the system as well as composition of the flow. In 1979, Liang [6] proposed one of the first models of an optically thin accretion disk solution including pair production. In this work he assumed the generation of e + − e − pairs through interaction of photons produced by unsaturated Comptonization and found an order of magnitude reduction in electron temperature because of enhanced cooling by these pairs. Unfortunately, detailed analysis was not presented. Also, the treatment was made for subrelativistic and trans-relativistic electron temperatures, assuming large optical depths, which allowed the photons to scatter efficiently. This scenario is not true in all astrophysical cases.   [7] studied the importance of pairs in AGNs and concluded that luminosities are reduced when there is sufficient pair production. Kusunose & Takahara (1988, 1990 [8,9,10] in a series of papers extensively discussed the effect of pairs in twotemperature accretion disks assuming pair equilibrium. These papers covered almost all the photon generation processes leading to the production of pairs : bremsstrahlung, synchrotron, inverse-Comptonization of bremsstrahlung and synchrotron photons (also see White & Lightman (1989) [11]) and external soft photons. But, the ion and electron energy advection term was neglected. Main conclusion of these papers were, for any given accretion rate, there existed two branches of solutions, one with a high pair to ion density ratio (defined as z) and one with a low z. The solution with a high z would generate unphysically high proton temperatures, halfheights and even long relaxation times. These were also unstable against thermal, secular and pair density perturbations. On the contrary, the solution with a low z was stable with negligible effect on the electron temperature or spectra. The authors also showed that a forbidden region will be formed when accretion rates exceed a certain critical value. In this region no steadystate solutions are possible. The reason for formation of this region is the enforcement of pair equilibrium, even in situations where pair production rate exceed the pair annihilation rate. This value of critical accretion rate is characteristic of the type of photon distribution.
Most of the initial works dealing with pair production either assumed a static plasma flow with the assumption of pair equilibrium [12,13,14] or a fixed efficiency of the system was considered [15]. Although these studies showed the importance of pair production in different accreting systems, they were not self-consistent. Shull (1979) [16] and Wandel, Yahil, and Milgrom (1984) [17] obtained self-consistent solutions but detailed work on the dynamics of the flow was not done and the radiative transfer part was not properly dealt with. They did not consider the effect of produced pairs back onto the system. Lightman (1982) [18] correctly quoted, that the problem of pair-production is highly non-linear; "pairs produce photons and photons produce more pairs". Yahel & Brinkman (1981) [19] and Yahel (1982) [20] investigated pair production in the most self-consistent way. But their methodology to obtain solutions were too sensitive, which inhibited them from getting a larger picture of the parameter space. Yahel & Brinkman (1981) concluded that there will be a characteristic spectral signature, a bump like feature near m e c 2 , due to e + − e − pairs, unless the number density of pairs is very low. Following them, more realistic self-consistent two-temperature work including pairs was done by Park & Ostriker (1989) [21] and Park (1990) [22] where they assumed spherical flows and implemented radiation hydrodynamic equations. Radiative processes considered were : line cooling, photoionizational heating, bremsstrahlung, Comptonization, pair creation and annihilation. They also considered the effects of generated radiation on the dynamics of the gas flow. The authors concluded that including pair processes in two-temperature solutions does not affect the accretion systems significantly. Advection-dominated disks soon came into existence and Kusunose & Mineshige (1996) [23] investigated the effect of pairs in these disks [24]. They concluded that pairs have negligible effect in these systems.
The discrepancies cited above were due to the fact that different works assumed different radiative processes and hence different radiation field producing pairs. In order to obtain a larger picture, we investigated this matter in greater details in the general-relativistic regime for one-temperature flows. There are three processes which could lead to pair production : photon-MHD-PP-2020 Journal of Physics: Conference Series 1640 (2020) 012022 IOP Publishing doi:10.1088/1742-6596/1640/1/012022 3 particle (electron, positron or proton), particle-particle and photon-photon interaction. Out of them photon-photon interaction is the most dominant process responsible for the production of e + − e − pairs [25]. Particle-particle interactions and photon-particle interactions are of less importance because their reaction cross-section decreases in the order of (1/137) (value of fine structure constant) and (1/137) 2 respectively. Hence, we can safely ignore them. In this work we assume that there are no magnetic fields present in the system. Thus, as a preliminary step, we include only bremsstrahlung and its inverse-Comptonization. The photons generated due to bremsstrahlung are in general, hard (of higher energy). If a photon's energy is lower than the energy of the ambient electrons, they can get upscattered to higher energies through the process of inverse-Comptonization. The probability for a photon to get scattered depends on the optical depth of the system. Thus, in the present work, the radiation field is made of a flat unscattered bremsstrahlung spectrum and a Wien tail. If these photons satisfy the criterion for pair production, e + − e − pairs would be produced [24]. Depending on the number density of the pairs produced, simultaneous annihilation will occur. After the annihilation process the radiation field will also contain these photons apart from the usual bremsstrahlung and inverse-Compton photon distribution.
Flows around BHs are trans-relativistic in nature. They are non-relativistic (kT < m e c 2 , adiabatic index=Γ ∼ 5/3) very far away from the horizon. As matter approaches the BH, it is compressed to smaller and smaller regions, with radiative processes playing significant role, the flow becomes mildly relativistic or relativistic (kT m e c 2 , Γ ∼ 4/3). The exact equation of state (EoS) for such a flow was given by Chandrashekar in 1939 [26]. This equation is computationally difficult to implement. Thus, we have used the Chattopadhyay & Ryu (2009) [27] EoS (CR EoS hereafter) for multispecies flow with variable adiabatic index. This EoS is analytical and exact which helped us to deal away with specifying any value of adiabatic index for each species at each point of the flow.
The paper is divided into the following sections. In section (2), we will cover the basic equations and assumptions used to model the flow. In section (3), we will discuss the methodology to find a self-consistent one-temperature transonic solution, including pair processes. We will show some important results in section (4) and finally conclude in section (5).

Assumptions and equations used to model the flow
We assume an advective viscous transonic accretion disk in hydrostatic equilibrium around a Schwarzschild BH. The treatment is in the pure general-relativistic regime where the line element in terms of spherical coordinate system is given by [28,29]: where, g µν s are the metric tensor components. The non-zero components are g tt = − 1 − 2 r ; g rr = 1 − 2 r −1 ; g θθ = r 2 ; g φφ = r 2 sin 2 θ. We have employed a system of units where, G = M = c = 1, such that the unit of length is in r g = GM BH /c 2 (where G = Gravitational constant, M BH = mass of BH) and the unit of velocity is in terms of light speed c. This unit system has been used throughout the paper unless mentioned otherwise. We assume the flow to be in steady state and around the equatorial plane, hence ∂/∂t = 0 and θ = 90 • , respectively. In addition, we assume fully ionized plasma and charge neutrality in the system. Accretion flows are governed by the following equations : (1) Continuity equation: (nu ν ) ;ν = S i − S o , where, n = number density of the particles and u µ s are the four-velocity components, and S i and S o are source and sink terms of pair production. In presence of pair production, the particle four-flux are not conserved. So, we define separately, the continuity equation for each particle (proton: n p , electron: n e − and positron: n e + ). (a) Proton: In the process of pair production, proton number density is not affected and hence the continuity equation can be integrated to obtain a form of conservation of proton mass flux or the proton-accretion rate of the system.
where, ρ p is the proton mass density and H is the half-height at each radius of the disk, the disk being in hydrostatic equilibrium. The form of half-height is obtained using the prescription given by Lasota (1994) [30,31].
(b) Positron: The continuity equation for positron can be written as [32] : where,ṅ C andṅ A are the creation and annihilation rates per unit volume, expressions of which will be discussed later.
(c) Electron : Charge neutrality demands that the number of positive charges equal the number of negative charges. This gives us the electron number where : (2) Conservation of energy-momentum : T µν ;ν = 0, where T µν is the energy-momentum tensor of a fully ionized viscous fluid and is given by : where, e, p and t µν are the local internal energy density, local isotropic gas pressure and viscous stress tensor, respectively. Projecting the four-divergence of T µν along the space direction, gives the relativistic Navier-Stokes equation, radial component of which is: and the integrated form of the azimuthal component gives : Here, L = hu φ and L 0 is the local bulk angular momentum defined at each radius and that defined at the horizon, respectively [33]. The specific angular momentum of the fluid is defined as Assuming shear to be the only factor giving rise to viscosity, we have t µν = −2ησ µν , where, η is the viscosity coefficient, defined as η = ρhν. Here,ν = α v ar(1 − v 2 ) is defined as the kinematic viscosity, with α v being the usual Shakura & Sunyaev viscosity parameter [34] and v = u r u r /(1 + u r u r ) is the radial three-velocity defined in the local co-rotating frame. Enthalpy and speed of sound in the medium is defined using h and a, respectively.
(3) First law of thermodynamics : u µ T µν ;ν = 0, which can be simplified using Eq.(5) to : where, ∆Q = Q + − Q − is the difference in the heating and cooling rate of the system.

Equation of state and the simplified equations of motion
As discussed before, we use the CR EoS which is given by, where,K = 1 + 1/χ, ξ = n p /n e − , χ = m e /m p and Θ = kT mec 2 is the non-dimensional form of temperature defined w.r.t the rest mass energy of the electron. Polytropic and adiabatic indices are defined as, N = 1 2 df dΘ and Γ = 1 + 1 N respectively. Using EoS, Eq.(9), we can simplify the first law of thermodynamics Eq.(8), to obtain : where, The term P is used to calculate the leptons internal energy consumed in producing pairs. If we simplify the azimuthal component of Navier-Stokes equation Eq.(7), we have : here, l = u φ . Now, we can obtain the differential equation for velocity from Eq.(6), using Eqs. (2)(3)(4) and (9)(10)(11), to obtain: where, The sound speed is defined as, a = Γp e+p . If we integrate Eq.(6) using Eq.(8) we obtain the generalised Bernoulli constant : which is a constant of motion even in the presence of dissipation. When the system is adiabatic and there is no dissipation E → E = −hu t , which is the canonical form of relativistic Bernoulli constant [31].

Pair production and annihilation processes
We assume pair production only through photon-photon interactions as mentioned before : γ + γ → e + + e − . The threshold condition is E 1 E 2 (1 − cos θ pp ) ≥ 2(m e c 2 ) 2 , suggesting that photons of energy > m e c 2 can only contribute to this process. Photon production mechanisms considered in this work are bremsstrahlung, inverse-Comptonization and the pair-annihilation photons.
Bremsstrahlung photons (ṅ br ) can be produced by e ± − p, e ± − e ± and e + − e − interactions. Their expressions are given in Svensson (1984) [14] and White & Lightman (1989) [11]. The photons generated through inverse-Comptonization (ṅ ic ) can be calculated using the forṁ n ic = f bṅbr , where f b is the fraction of emitted bremsstrahlung photons that are scattered upto the Wien peak before escaping the disk. Thus, it depends on the optical depth of the system as well as the temperature of the ambient electrons. It is calculated using the prescription given by Nakamura et. al. (1996) [35]. For computing the photons produced due to annihilation (ṅ ann )  6 we use the expression given by White & Lightman (1989) [11]. Thus, to calculate the total number of photons contributing to pair-production (n γ ), we solve the photon balance equation, which reads: where, t esc = H/c(1 + τ es ), is the photon escape time and τ es is the optical depth of the system. After we have information of the radiation fields present in the system, we identify all the possible combinations, which can lead to pair production. They are as follows : collisions between (1) photons of the same Wien peak (ṅ WW ) , (2) between Wien and flat unscattered bremsstrahlung photons (ṅ WF ), (3) between the photons belonging to the same unscattered flat part (ṅ FF ), (4) Wien photon and the annihilation photons (ṅ WA ), (5) flat part of bremsstrahlung photons with annihilation photons (ṅ FA ) and (6) between the same distribution of annihilation photons (ṅ AA ). Thus the creation rate can be written as : All these 6 scenarios have been included in this work, expressions of which are taken from White & Lightman (1989) [11] and Esin (1999) [24]. The annihilation rate can be written as [13] : where, r e is the electron radius.

Radiative processes
The radiative processes considered are as follows : • Electron is cooled via three processes : (1) bremsstrahlung (Q br ) radiation (through e ± − p, e ± − e ± and e + − e − collisions [11]) (2) inverse-Comptonization of bremsstrahlung photons (Q ic ) [11] and (3) annihilation (Q ann ) [12] • Viscous dissipation of energy leads to heating in the system. Shakura & Sunyaev α v parameter [34] controls the amount of this dissipation. • The pairs produced, have same thermal energies as the other particles present in the system, thus contributing to an overall heating term : Q pair−heat =< hν 1 + hν 2 >ṅ C . • If the temperature of the electrons are less than the energy of the ambient photons, then electrons can heat up via Compton heating (Q comp ) mechanism [29].
Thus to summarize : Q + = Q B + Q pair−heat + Q comp and Q − = Q br − Q ic − Q ann . In addition, we have implemented all the special and general relativistic transformations while computing the actual energy loss as well as while plotting the spectrum.

Sonic point conditions and shock conditions
Mathematically a sonic point is defined as the point where dv/dr has a 0/0 form. Hence, from Eq.(12), the sonic point conditions are: N = 0 and D = 0. To find the derivative of the velocity at the sonic point (dv/dr| c ) we need to use the L'Hospital rule. Presence of angular momentum and relativistic effects causes the multiplicity of sonic points. Depending on the flow parameters, an accretion solution may harbour one or more sonic points [29]. The sonic points are named according to its distance from the BH : inner (r ci ), middle (r cm ) and outer (r co ). Inner and outer sonic points are X-type sonic points and are physical in nature, i .e. flows pass through them, while O-type sonic points are unphysical in nature. Due to the multiplicity of sonic points, shocks can also be formed. Supersonic matter (after it has passed through r co ), due to the twin effect of centrifugal and thermal pressure could act as a barrier to the matter following it, inducing the formation of a centrifugal-force mediated shock transition. After the shock jump, matter becomes subsonic, but crosses the horizon supersonically after passing through r ci . The relativistic shock conditions [38]  v v 2 +pH] = 0, where, the square brackets denotes the difference of the quantities across the shock.

Methodology to obtain a solution 3.1. Finding a general transonic solution without pairs
Here, we discuss the method to find a general one-temperature advective transonic solution without any pair processes. The system we are dealing with, is still dissipative (viscous dissipation as well as radiative processes are included), hence sonic points are not known a priori. For a given a set of flow parameters (E, λ in , α v , β d ,Ṁ and M BH ), a solution may harbour one or more sonic points (see Section 2.4). Finding of sonic points is not trivial. The first step is to specify the boundary conditions. For this, we select a point asymptotically close to the horizon, say, r in = 2.001 (in terms of r g ), horizon being a co-ordinate singularity. The reason for this choice is that, near the horizon, the value of the variables are known, r → 2r g and v → c. Also, E → E = −hu t , because very near the horizon, gravity overpowers any process or interaction, preventing the matter to either interact with itself or its environment, making the region adiabatic in nature. We now give a step-by-step procedure for finding solutions: (1) Obtaining boundary conditions (v in , l in and L 0 ): We first assume any arbitrary value of Θ in . We evaluate the value of l in from the definition of specific angular momentum, λ in = −u φ /u t = −l in /(E/h in ), h in being a function of only Θ in . Simplifying E = −h in u t near the horizon, gives us an analytical expression for velocity near the horizon (v in ). Hence, supplying the value of E, h in (Θ in ) and l in (as obtained before), would give us the value of v in . The adiabaticity condition near the horizon reads E E, which implies γ φ exp(X f ) = 1, from Eq. (13). Calculating its derivative w.r.t r and then simplifying it, would give us a quadratic equation in (L in − L 0 ), (where, L in = l in h in ) from which we can find the value of L 0 .
(2) Finding sonic point (obtaining Θ in ): After finding v in , l in and L 0 for a given set of flow parameters and a given arbitrary Θ in , we integrate the equations of motion Eq.(3) and Eq.(10-12) from r in outwards, towards infinity using 4 th order Runga-Kutta integration technique. We check whether the solution obtained passes through the sonic point or not. If it does not, then we change Θ in , until the sonic point conditions are satisfied. We note that, although it is a pair-free iteration, we calculate the maximum positron density that could be obtained with the current radiative processes.
(3) Obtaining the transonic solution : As soon as we get the sonic point, we find the derivative dv/dr|c using L' Hospital rule and integrate further outwards until we reach infinity (around r → 10 4 r g ) or the velocity of the solution becomes too low. (4) Checking for other sonic points : Since angular momentum and relativistic effects induces the formation of multiple sonic points, we need to check whether any other sonic point exists for the same L 0 . For this we change Θ in by a larger value and repeat the steps described from 1-2. If there is another sonic point present then we follow step 3 to obtain a transonic solution. (5) Checking of shock transitions : When a certain set of flow parameters harbour multiple sonic points, we need to check whether the solutions passing through them satisfy the shock conditions (see, Section 2.4 for shock conditions). If, satisfied, the solution first passes through the outer sonic point (r co ) encounters a shock (r s ) and then passes through the inner sonic point (r ci ), entering the horizon supersonically. At the shock location, the velocity of the flow becomes subsonic while temperature and density jumps to higher values. The methodology described above is for a pair-free model. We name it as iteration 1. In iteration 2, we repeat steps 1-5, mentioned in the above section, using the positron number density calculated in iteration 1. But now, we include the pair creation and annihilation processes along with the processes included before in iteration 1. We continue this process of iterating, until the solutions converge.
To explain it further, we plot in Fig.(1a), Mach number, defined as, M = v/a and in Fig.(1b) the composition of the flow, defined using ξ = n p /n e − . The flow parameters used are, E = 1.001, λ in = 2.5,Ṁ = 0.6Ṁ edd , α v = 0.05 and M BH = 10M . Iteration 1 is denoted using red, dotted line. We can see that, the solution passes through a single sonic point (marked with a red star, in panel a, clearly visible in inset) and ξ = 1 for the first iteration (panel b), i.e. pair-free. As we perform iterations, there is a play between pair production and annihilation processes. The annihilation photons introduces an additional radiation field into the system, hence, there will be an increase in pair production and the interplay between these processes will go on. Coloured stars in panel (a) depicts the sonic point, which shifts to different locations with each iteration. We select the solution after 10 iterations or until the solutions converge, whichever is earlier. This is clearly visible from panels (a)  , in the presence of pair production and annihilation, for the given set of flow parameters, is shown in black, solid line. In panel (c), we plot the change in spectrum with iterations. Grey curves depict the emission due to annihilation, corresponding to each iteration, while the rest of the spectral emission comes from bremsstrahlung and its inverse-Comptonization. We see that, in iteration 1, there is no emission from annihilation, because of the absence of positrons. While, as we increase the number of iterations, the spectral shape changes depending on the pair processes and annihilation, ultimately converging to the spectrum shown in black, solid line. We investigate a general transonic one-temperature solution, including pair processes. The solution is obtained using the methodology elaborately described before. A typical solution is represented using Mach number (M ) and the covariant azimuthal component of four-velocity l or u φ plotted in Fig.(2a), while the composition of the flow ξ = n p /n e − = n p /(n p + n e + ) is plotted in Fig.(2b). Both of these have been plotted against log r. The different curves represent the different accretion rates of the system:Ṁ = 0.1Ṁ edd (red, dotted),Ṁ = 0.4Ṁ edd (green, dashed) andṀ = 0.8Ṁ edd (blue, solid). Other flow parameters are, E = 1.001, λ in = 2.60, α v = 0.01 and M BH = 10M . From the plot it is apparent that for these set of flow parameters only an outer sonic point exists. As we increase the accretion rate of the system the sonic point shifts inwards, towards the BH (clearly visible forṀ = 0.8Ṁ edd solution). This is mainly because of the increase in cooling of the system with the increase in accretion rate. We note that there is a minor decrease in the location of the sonic point for the solution represented by accretion rate 0.8Ṁ edd as compared to 0.4Ṁ edd . The solutions represented in panel (a) for these accretion rates is hence seen to be overlapped. Panel (b) interprets the presence of positrons in the flow. A decrease in ξ, suggests an increase in the number of positrons in the flow. We see that, at infinity and just near the horizon, ξ = 1, meaning that there are no positrons present in these regions. At infinity, radiative mechanisms are less effective due to low temperatures and low number densities. Thus, the radiation field giving rise to pair production is weak, leading to negligible production of positrons. While, near the horizon, photon trapping effects of the BH dominate, thereby reducing the number density of photons and hence, lowering of pair production rates occur. But, there could be ample amount of advected positrons and electrons. The annihilation rate is found to be higher in these regions reducing the number of positrons to ∼ 0. Now, as we increase the accretion rate of the system, matter supply becomes more, hence, there is more cooling in the system. Thus, the solution will take a higher temperature branch, to sustain the excessive cooling. This will, in turn, increase the pair production rate, contributing to the rise in positron number density in the flow and hence decreasing ξ (panel b). In panel (c), spectrums are plotted, where, with the increase in accretion rate, the system becomes more luminous. There is a bump in all the spectrums near ∼ m e c 2 , whose magnitude increases with the increase in accretion rate of the system. This bump is a signature of the amount of annihilation photons. More the accretion rate more would be the pairs produced and 10 hence more will be the photons generated through pair annihilation.

Conclusions
In this paper, we present one-temperature advective, viscous, transonic solutions, in the presence of cooling processes like bremsstrahlung, inverse-Comptonization as well as incorporating pair production and annihilation mechanisms. There was no fixed adiabatic index for the flow. CR EoS helped us in dealing with it. Pair production was assumed to be from the interactions between the photons generated through the above mentioned processes. We find that the number of pairs produced by bremsstrahlung, inverse-Compton and annihilation photons are not significant, since the minimum ξ is > 0.97 for an accretion rate ofṀ = 0.8Ṁ edd . For accretion rates lower than this, we can neglect the presence of any pairs. Also, the maximum number of pairs are produced in a region < 50r g . When the number density of pairs is not very low, like for accretion rates > 0.8Ṁ edd , there is a bump near m e c 2 in the spectrum. This is a spectral signature of the annihilated photons. We plan to carry out detailed analysis of this work in future.