Issue 
A&A
Volume 642, October 2020



Article Number  A209  
Number of page(s)  21  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/202037520  
Published online  23 October 2020 
Twotemperature solutions and emergent spectra from relativistic accretion discs around black holes
^{1}
Aryabhatta Research Institute of Observational Sciences (ARIES), Manora Peak, Nainital 263002, India
email: shilpa@aries.res.in, indra@aries.res.in
^{2}
Pt. Ravishankar Shukla University, Great Eastern Rd, Amanaka, Raipur, Chhattisgarh 492010, India
^{3}
IRFU/Service d’Astrophysique, Bat. 709 Orme des Merisiers, CEA Saclay, 91191 GifsurYvette Cedex, France
email: philippe.laurent@cea.fr
^{4}
Laboratoire Astroparticule et Cosmologie, Bâtiment Condorcet, 10 Rue Alice Domont et Léonie Duquet, 75205 Paris Cedex 13, France
Received:
17
January
2020
Accepted:
27
June
2020
Aims. We investigate a twotemperature advective transonic accretion disc around a black hole and analyse its spectrum in the presence of radiative processes such as bremsstrahlung, synchrotron, and inverseComptonisation. The aim is to link the emergent spectra with constants of motion of the accretion disc fluid, however, the number of unknowns in twotemperature theory exceeds the number of equations for a given set of constants of motion. We intend to remove the degeneracy using a general methodology and obtain a unique solution, along with its spectrum.
Methods. We used hydrodynamic equations (continuity, momentum, and energy conservation equation) to obtain sonic points and solutions. To solve these equations of motion we used the 4th order RungeKutta method. For the spectral analysis, general and special relativistic effects were taken into consideration. The system is, nonetheless, degenerate and we remove the degeneracy by choosing the solution with maximum entropy, as dictated by the second law of thermodynamics.
Results. We obtained a unique transonic solution for a given set of constants of motion. The entropy expression is a tool used to make a selection between the degenerate solutions. We found that Coulomb coupling is a weak energy exchange term, which allows protons and electrons to settle down into two different temperatures, justifying, hence, our study of twotemperature flows. The information of the electron flow allows us to model the spectra. We show that the spectra of accretion solutions depend on the associated constants of motion. At low accretion rates, bremsstrahlung is important. A fraction of the bremsstrahlung photons may be of higher energy than the neighbouring electrons, energising them through the process of Compton scattering. Synchrotron emission, on the other hand, provides soft photons, which can be inverseComptonised to produce a hard power law part in the spectrum. Luminosity increases with the increase in the accretion rate of the system, as well as with the increase in BH mass. However, the radiative efficiency of the flow has almost no dependence on the BH mass, but it sharply rises with the increase in the accretion rate. The spectral index, however, hardens with the increase in the accretion rate, while it does not change much with the variation in BH mass. In addition to the constants of motion, the value of the plasma beta parameter and magnitude of magnetic dissipation in the system also helps in shaping the spectrum. A shocked solution exists in twotemperature accretion flows in a limited region of the parameter space. We find that a shocked solution is always brighter than a solution without a shock.
Conclusions. An accreting system in twotemperature regime admits multiple solutions for the same set of constants of motion, producing widely different spectra. Comparing the observed spectrum with that derived from a randomly chosen accretion solution would give us a wrong estimation of the accretion parameters of the system. The form of entropy measurement we obtained helped us to remove the degeneracy of the solutions and allowed us to understand the physics of the system, shorn of arbitrary assumptions. In this work, we show how the spectra and luminosities of an accreting system depend on the constants of motion, producing solutions ranging from radiatively inefficient flows to luminous flows. An increase in BH mass quantitatively changes the system, making the system more luminous, while the spectral bandwidth also increases. A higher BH mass system spans the range from radio to gammarays. However, increasing the accretion rate around a BH of certain mass has little influence on the frequency range of the spectra.
Key words: hydrodynamics / accretion / accretion disks / black hole physics / shock waves / radiation mechanisms: general
© ESO 2020
1. Introduction
Accretion is the primary mechanism that may explain the observation of radiation coming from microquasars and active galactic nuclei (AGN). Stellar mass black holes (BH) and neutron stars are thought to reside in microquasars, whereas supermassive BHs reside in AGNs. The energy released due to the accretion of matter onto relativistic compact objects, such as neutron stars and BHs, is a fraction of the rest mass energy of the matter falling onto it. The study of the accretion process was initiated by Hoyle & Lyttleton (1939), where they studied accretion of matter onto a Newtonian star passing through the interstellar medium (ISM). In Bondi (1952) gave the first full analytical solution for spherical flows (also known as Bondi flows) around a static star. This theory is also relevant to the case of stellar winds. He showed that accretion and wind are characterised by a unique transonic solution having maximum entropy. But it took ten years, up until the discovery of quasars and Xray sources in the 1960s, that accretion phenomena gained in popularity. Salpeter (1964) and Zel’dovich (1964) extensively investigated accretion as a probable mechanism for driving and powering these luminous objects. It was concluded that the Bondi accretion model produced luminosities which are too low to explain the observations. As matter is radially falling and exhibiting short infall timescales compared to their cooling timescales, it leads to the low radiative efficiency of such flows. This finding led to the development of the famous Shakura & Sunyaev (1973) disc model (SSD) or the Keplerian disc model (Pringle & Rees 1972; Shakura & Sunyaev 1973; Novikov et al. 1973). In this model, it is assumed that matter is rotating in Keplerian orbits having negligible radial velocity an pressure gradient terms. The shear between differentially rotating matter gives rise to viscous stresses, which help in removing angular momentum outwards, allowing matter to spiral inwards and finally fall onto the central object. The time taken for inspiral allows for the matter to produce emission for a longer duration, unlike spherical flows. Since at every radius the angular momentum of the disc is Keplerian, the disc must therefore be geometrically thin. In other words, this implies that the heat generated due to viscous dissipation needs to be efficiently radiated away such that the angular momentum distribution remains close to the Keplerian one (hence the name “cool discs”, King 2012). In SSDs, matter and radiation are in thermal equilibrium, where each annulus of the disc emits a blackbody spectrum (or, depending on opacity, a modified version of it) that peaks at the temperature of the annulus of the disc. The composite spectrum is hence the sum of these blackbodies and is called the modified multicoloured blackbody spectrum. Although this model could successfully regenerate the thermal part of the spectrum from sources associated with BHs, it was unable to explain the nonthermal part. In addition, the assumption of Keplerian angular velocity at each annulus implied that the disc is arbitrarily terminated at the inner stable circular orbit (ISCO). Soon after, it was also concluded that these discs were thermally and secularly unstable (Pringle et al. 1973; Lightman & Eardley 1974; Artemova et al. 1996). In Thorne & Price (1975) argued that the instability present at the inner region of an SSD could expand into an optically thin gaspressure dominated region. Shapiro et al. (1976, hereafter SLE76) considered this geometrically thick and optically thin, puffedup region to be composed of protons and electrons described by two different temperature distributions. Using this model, they successfully reproduced the hard component part of Cygnus X1 spectrum from 8–500 keV. Unfortunately, the SLE model was also found to be thermally unstable (Piran 1978). If the disc is heated, it expands, reducing its number density and thereby its cooling rate. This makes the system even hotter, leading to a runaway thermal instability. Although the system proved to be unstable, this study has served as one of the cornerstones in twotemperature accretion theory.
The models discussed above have suffered from simplifying assumptions, for example, the cooling rate at each radius of the disc was equated with the heating rate and the advection term was not properly dealt with. In general, the heating and cooling rates do not need to be equal and some part of the heat could be advected inwards, along with the bulk motion of the system. Abramowicz et al. (1988) extensively investigated advection, in their “slim” optically thick accretion disc model, and found that the solutions obtained were thermally and viscously stable. The importance of advection was further demonstrated, using self similar solutions in the works of Narayan & Yi (1994), Abramowicz et al. (1995). Today, these discs are broadly classified as advectiondominated accretion flows (ADAFs) and Ichimaru (1977) was the first to propose it (for review also see BisnovatyiKogan & Lovelace 1997, 2001).
A general conclusion can be drawn from the above models that an accretion flow does not need to be Keplerian throughout but it could also be subKeplerian or a combination of both. In addition, the flow must be transonic. Matter that is very far away from the horizon is subsonic, whereas the BH boundary condition insists that matter should cross the horizon at the speed of light. Thus, accreting matter has to pass through at least one sonic point, or in other words, the BH accretion solution is necessarily transonic. Bondi (1952) in his seminal paper already highlighted the importance of transonicity for spherical accretion flows. Liang & Thompson (1980; hereafter LT80) argued that similarly to spherical flows, which are characterised by a single sonic point, rotating flows around BHs are characterised by multiple sonic points. Fukue (1987) extended their work and presented in detail the nature of sonic points in transonic rotating accretion flows. He concluded that even though a BH does not possess any hard surface, it can undergo a shock transition in the presence of multiple sonic points (also see Chakrabarti 1989).
All the works mentioned above (except SLE76) assume onetemperature accretion flows. Onetemperature flows are based on the fact that the timescale of the energy exchange process (like Coulomb coupling) between the ions and electrons is shorter or comparable to the dynamical time scale of the system, allowing the system to effectively settle down into a single temperature distribution (Le & Becker 2005; Chattopadhyay & Chakrabarti 2011; Kumar & Chattopadhyay 2014). Yet, in many astrophysical cases, infall timescales are generally much shorter than Coulomb collision timescales, that is, Coulomb coupling between the protons and electrons are not strong, allowing the two species to equilibrate to two different temperature distributions (Rees et al. 1982; Narayan & Yi 1995). Therefore, in addition to the advective transonic nature of an accretion flow around a BH, the gas is likely to be in the twotemperature regime. Also, the electrons are more prone to radiative cooling, compared to ions. This makes the electron temperature deviate largely from the protons especially in the inner regions of the accretion disc.
After the seminal paper of SLE76, a twotemperature assumption was largely used to model accretion flows as it could successfully reproduce the observed spectrum, electrons being the main radiators. Colpi et al. (1984) included advection and solved the energy equation (or, first law of thermodynamics) for spherical accretion flows (i.e. with no angular momentum), and assumed a freefall velocity field with radial dependence taking the form: v ∝ r^{−1/2}. The transonic nature of the flow was not taken into account, but emission processes relative to protons and electrons were discussed briefly. Narayan & Yi (1995; hereafter NY95) incorporated angular momentum and extensively discussed the nature of twotemperature, optically thin accretion discs. The equations of motion were solved under the selfsimilar assumption. It is to be noted that selfsimilarity in accretion solution around BH is plausible only at a large distance from the horizon and not near it. Additionally, a selfsimilar solution is not transonic. NY95 neglected the electron advection term and assumed the heating and cooling rates of electrons to be equal. There have been other notable works done in a twotemperature regime, where the issue of transonicity was bypassed, but the radiative transfer part was properly treated nonetheless. One such work is by Chakrabarti & Titarchuk (1995), where the spectrum was computed assuming twocomponents in the accretion flow: a Keplerian and another subKeplerian. An exact Comptonisation model was incorporated, following a discussion of the variation in the spectrum with the change in mass of BH and accretion rate. Mandal & Chakrabarti (2005) went a step further and included a nonthermal distribution of electrons along with the general thermal distribution inside the accretion flow. Similar to Colpi et al. (1984), in this paper the velocity field was assumed to be in freefall, but the shocked solution and its signature in the observed spectrum were qualitatively studied.
Nakamura et al. (1996) obtained the first global transonic, twotemperature accretion solutions in which the advection terms were considered without any simplifying assumptions. However, at the outer boundary, the authors assumed the ion temperature to be a fraction of the virial temperature and Coulomb coupling was equated with bremsstrahlung. Manmoto et al. (1997) extended the work of Nakamura et al. (1996) to calculate the spectrum. In addition, they slightly modified the outer boundary conditions, where the total gas temperature (and not the ion temperature) was assumed to be a fraction of the virial temperature and equated Coulomb coupling with the total cooling of electrons (bremsstrahlung, synchrotron, and inverseComptonisation). Recently, a more general transonic advective twotemperature accretion solution was obtained by Rajesh & Mukhopadhyay (2010), where the viscous stress was assumed to be proportional to the sum of ram pressure and gas pressure (see Chakrabarti 1996, for details on this particular viscosity prescription). In this work, a limited class of solutions were studied. This work was extended and a global class of transonic twotemperature solutions, including accretionshock, were investigated by Dihingia et al. (2017). The most striking feature of the transonic twotemperature studies mentioned above is that the solutions depend on the choice of inner or outer boundary conditions, however, based on singletemperature hydrodynamics, we know that transonic solutions are unique for a given set of constants of motion. We highlight this issue in greater detail below.
Hydrodynamic equations, even in the single temperature regime, admit an infinite number of solutions, but a transonic solution is physically favoured because it has the highest entropy among all possible global solutions (Bondi 1952). In addition, the location of the sonic point corresponds to a unique boundary condition. The number of hydrodynamic equations for flows in onetemperature and twotemperature regimes are exactly the same, but there is one more variable in case of twotemperature flows (i.e. the existence of different ion and electron temperature distributions). Obtaining a selfconsistent twotemperature flow would require solving basic hydrodynamic equations, but since there is one more flow variable in the twotemperature regime, the system is degenerate. We obtain a large number of transonic solutions for the same set of constants of motion, of the flow. A hint of the problem of degeneracy was reported briefly in LT80, but it was skirted by parameterising the temperatures of protons and electrons to some constant value, arguing that the coupling between these species is unknown. As mentioned earlier in this paper, several authors also assumed an arbitrary proton or electron temperature at the boundary where they started integrating to find possible solutions. A global treatment of twotemperature problem would require all the equations to be solved selfconsistently without taking recourse to any set of arbitrary assumptions on temperature values at any boundary. This problem of degeneracy was identified and a prescription for obtaining a unique transonic twotemperature solution was reported in Sarkar & Chattopadhyay (2019, hereafter SC19). However it was applied to flows having zero angular momentum. Spherical flows have single sonic points, which simplifies our problem, allowing us to focus only on the issue of degeneracy in a twotemperature regime. In SC19, we reported that for a given set of constants of motion, infinite transonic solutions exist, each having a sonic point property different from the rest. The question that arises is which solution should be chosen since nature generally disfavours degeneracy.
In a onetemperature regime, the transonic solution is unique (Le & Becker 2005; Becker et al. 2008) but for reasons cited above, a twotemperature solution is degenerate for the same set of constants of motion. Following Bondi (1952), we can look for the highest entropy solution to remove the degeneracy. However, the twotemperature energy equation within the adiabatic limit is not integrable (unlike in onetemperature regime, Kumar et al. 2013; Kumar & Chattopadhyay 2013, 2014; Chattopadhyay & Kumar 2016). This prevented us from obtaining an analytical expression for measuring entropy in the twotemperature regime. The integration of the energy equation is spoiled by the presence of the Coulomb coupling term. However keeping in mind the fact that near the horizon, gravity overpowers any other interaction or processes, we can neglect the Coulomb coupling term. Thus, as we reported in SC19 for the first time, we can use a form of entropy measurement which would only be applied near the horizon. Using the formula, we measured the entropies at the horizon for the transonic spherical twotemperature solutions obtained for a given set of constants of motion. We saw that the entropy was maximum for a particular solution. Based on the second law of thermodynamics, we selected the solution with maximum entropy as the solution that nature would prefer. This solved the problem of degeneracy in a twotemperature model. It may, however be noted that spherical flows were simple to handle thanks to the presence of single sonic points.
It may be noted that much like the flow speed, the thermal state of the flow around a BH is also transrelativistic in nature. This means that very far away from the horizon, matter is thermally nonrelativistic and as it approaches the BH, it could be subrelativistic or relativistic. Matter is referred as thermally relativistic when its thermal energy is comparable to or greater than its rest mass energy, kT/mc^{2} ≳ 1 (where k= Boltzmann constant, T= temperature, m= mass of the species), and its adiabatic index, (Γ) ∼4/3; it is nonrelativistic if its thermal energy is less than the rest mass energy (kT/mc^{2} < 1) and Γ ∼ 5/3. So, not only the temperature but also the mass of the constituent particles determine the relativistic nature of the flow. So, in twotemperature flows, where we consider two species with masses differing by ≳1000 times, an equation of state (EoS), with a fixed adiabatic index, is untenable. Chattopadhyay (2008) and Chattopadhyay & Ryu (2009), proposed an approximate EoS for multispecies flow (CR EoS, hereafter) that is analytical and computationally easy to handle. Although it is an approximation, it matches perfectly well (Vyas et al. 2015) with the relativistically perfect EoS obtained by Chandrasekhar (1939). In a later work, the authors extended the use of CR EoS for dissipative, relativistic accretion disc as well (Chattopadhyay & Kumar 2016). To incorporate the transrelativistic nature of protons and electrons, it is necessary to utilise a variable adiabatic index EoS.
In this paper, we would like to extend our previous Sarkar & Chattopadhyay (2019) model to rotating flows and, thereby, to remove the degeneracy of accretion solutions around BHs. Needless to say, choosing any one of the degenerate solutions would give us an incorrect picture of the accretion parameters of the BH. To handle the transrelativistic nature of these flows, we use the CR EoS, which freed us from having to specify the adiabatic indices of each species. Hydrostatic equilibrium along the disc thickness is maintained at each radius of the disc. The radiative processes with proper special and general relativistic corrections have been incorporated to compute the spectrum.
In this paper, we aim to analyse how the correct unique accretion solution depends on the constants of motion of the flow, along with the mass of the BH, and to discuss how these properties play an active role in shaping the spectrum. Further, we intend to study the relative contribution of various radiative processes on the emitted spectrum and also investigate which part of the disc is likely to contribute the most in the emitted spectrum. In addition, we want to check whether an accretion shock imparts any special spectral signature. Apart from this, we would like to study the dependence of radiative efficiency and the spectral index with the variation in the accretion rate of the BH and its mass, in the presence of different radiative processes.
This paper is divided into the following sections. In Sect. 2, we give a brief overview of the basic equations and assumptions used to model the flow. In Sect. 3, we discuss the solution procedure for finding a unique transonic twotemperature solution. We present the results in Sect. 4 and our conclusions in Sect. 5.
2. Twotemperature advective disc model: assumptions and governing equations
In this paper, our main intention is to obtain all possible accretion solutions onto a BH and compute the typical spectrum corresponding to each mode of accretion. The solution involving a rotating accretion flow is much more complicated than spherical accretion due to the presence of multiple sonic points. Therefore, the spectra are different for different kinds of solutions. Furthermore, twotemperature flows are degenerate, which makes the method for obtaining a unique solution very important. In the next subsections, we discuss the equations we used to model the twotemperature accretion flow. We also give an overview of the radiative processes we incorporated and the methodology we used to implement these processes in curved spacetime.
2.1. Equations of motion (EoM)
The spacetime metric around a nonrotating BH is described using a Schwarzschild metric:
where the metric tensors are expressed as, and g_{θθ} = g_{ϕϕ} = r^{2}, since an accretion flow is described around the equatorial plane. Here, r, θ, and ϕ are the usual spherical coordinates and t is the time coordinate, G= Gravitational constant, and M_{BH}= mass of the BH. We note that throughout the paper, we employ a system of units where the unit of length, velocity, and time are defined as r_{g} = GM_{BH}/c^{2}, c and r_{g}/c = GM_{BH}/c^{3}, respectively. All the variables used in the rest of the paper are written in this unit system unless mentioned otherwise. The BH system that we model is in a steady state and axissymmetric, therefore ∂/∂t = ∂/∂ϕ = 0. Moreover, at any radius, we assume that only the radial gradient of any quantity is dominant, therefore, ∂/∂θ = 0.
The radial component of the momentum balance equation is:
where e and p are the internal energy density and isotropic gas pressure, respectively, measured in the local fluid frame, and u^{μ}s are the components of fourvelocity. The mass accretion rate is obtained by integrating the conservation of four massflux:
where ρ = n(m_{p} + m_{e}) = local mass density of the flow, n is the particle number density, m_{p} and m_{e} are the mass of proton and electron, respectively, and H is the local halfheight of the disc. Here, Ṁ – the accretion rate, is a constant of motion throughout the flow. Halfheight is calculated assuming hydrostatic equilibrium along the vertical direction of the disc (Lasota 1994; Chattopadhyay & Chakrabarti 2011) which can be written as
Also, based on u_{μ}u^{μ} = −1, we obtain:
where, γ_{v} and γ_{ϕ} are the Lorentz factors in the radial and azimuthal directions, respectively and defined as and , where and v is the velocity in the local corotating frame. It can be shown that , where . The total Lorentz factor, therefore is, γ = γ_{v}γ_{ϕ}.
The first law of thermodynamics, or the energy balance equation, is and can be written as:
Here, ΔQ = Q^{+} − Q^{−}, where, Q^{+} and Q^{−} represent the rate of heating and cooling present in the flow for all the species, respectively. These rates are in units of ergs cm^{−3} s^{−1}, which are converted into geometric units before being applied in the above equation.
Coulomb coupling serves as an energy exchange process, which transfers energy between protons and electrons. In the singletemperature case, with Coulomb coupling being infinitely strong, protons and electrons equilibrate locally into a singletemperature distribution (Lanzafame et al. 1998; Lee et al. 2011, 2016; Chattopadhyay & Kumar 2016). However for twotemperature flows, it is not strong enough, hence, protons and electrons are allowed to thermalise at two different temperatures. In other words, the timescale required for protons and electrons to attain a thermal equilibrium and settle down into a single temperature is more than the timescale in which each of the two populations thermalise separately. To describe such flows, we need to use two separate energy equations, one for protons and another for electrons. These two energy equations are not independent and, as discussed, are connected by the Coulomb coupling term which acts as a cooling term for protons and a heating term for electrons if the proton temperature is higher than electron temperature.
If we integrate the equation of motion (Eq. (2)) with the help of the energy equation (Eq. (5)), we obtain the generalised Bernoulli constant, given by
where, h = (e + p)/ρ= specific enthalpy and . Here, represents the difference in the heating and cooling rates of the ith species. The generalised Bernoulli constant is conserved all throughout the flow, even in the presence of heating and cooling. The X_{f} term mainly comes up due to the presence of dissipation. In case of adiabatic flows, with no dissipation, X_{f} = 0 and E → ℰ = −hu_{t}, which is the canonical form of relativistic Bernoulli constant (Lightman et al. 1975; Chattopadhyay & Chakrabarti 2011).
2.2. EoS and the final form of the EoM
We ought to note that in this subsection, barred variables have been used to denote dimensional quantities and unbarred variables are, as before, nondimensional. In order to solve the above equations of motion, we need to supply an equation of state (EoS). In this paper, we use the ChattopadhyayRyu (CR) EoS given by Chattopadhyay (2008), Chattopadhyay & Ryu (2009). The form of CR EoS for multispecies flow is:
where, the summation is over ith species. Since, in this paper, we consider the accretion flow to be composed of protons and electrons (e^{−} − p^{+}) only, thus i represents these two species. Number density (n), mass density (ρ), and isotropic gas pressure (p), present in Eq. (7), can be represented in dimensional form in the following way:
where, η = m_{e}/m_{p} and . T_{i} is the temperature in units of kelvin, while is the non – dimensional temperature
defined with respect to the restmass energy of the respective ith species. The EoS, Eq. (7), can be simplified using Eqs. (8)–(10) to
where, and
Polytropic index and adiabatic index can be written as
The equation for halfheight (Eq. (4)), can be simplified to:
Here, λ = −u_{ϕ}/u_{t} is the specific angular momentum of the flow. Angular momentum plays a very important role in accretion disc physics because it can significantly modify the infall time scale. Using Eqs. (8)–(13), we can simplify the energy equation (Eq. (5)) and obtain two differential equations for temperature, one for protons and another for electrons. They are as follows:
respectively, where,
, and
If we simplify the radial component of the relativistic momentum balance equation (Eq. (2)) using Eqs. (8)–(15), we get the expression of gradient of threevelocity, which has the form:
where, , and
The effective sound speed (a) has been defined as
where,
2.3. Heating and cooling processes included in the flow
In the following subsections, we briefly discuss the processes that lead to the heating and cooling of plasma in the accretion flow.
2.3.1. Heating due to magnetic dissipation
The magnetic field in the medium surrounding the BH would be frozen into the highly conductive infalling plasma. As the matter falls inwards, its magnetic field strength would increase by 1/r^{2} and magnetic energy density (B^{2}/8π) by 1/r^{4}. Schwartzman (1971) argued that before the magnetic energy density exceeds the thermal energy density, turbulence and hydromagnetic instabilities would lead to reconnection of magnetic field lines. In other words, it means that the magnetic energy density is limited by equipartition with thermal energy density. This magnetic energy dissipated would heat up the matter, either protons or electrons or both, ensuring relativistic temperatures even far away from the BH. The expression for the dissipative heating rate, as given by Ipser & Price (1982) is
The above equation is a measure of the heating due to magnetic dissipation. However, there are uncertainties in the estimates of heating processes that are controlled by β_{d}. We used β_{d} = 0.001, unless stated otherwise, as a representative case. In Sect. 4.6, we varied the value of β_{d} and studied how the solutions depend on it. In this work, we assume that both protons and electrons can absorb the magnetic energy dissipated. Hence, we can write
where δ is the uncertainty parameter, which dictates the amount of heat absorbed by protons, the rest being absorbed by electrons. There is insufficient knowledge available in the literature relating to this issue. Thus, throughout our work, for simplicity, we consider δ = 0.5, which means that 50% of this heat would go to protons and the remaining 50% into electrons, unless otherwise specified.
In our present study, we investigate inviscid flows since the proper handling of the general relativistic (GR) form of viscosity in transonic flows is not trivial. The shear tensor in GR contains a derivative of u_{ϕ}, v and other terms (Peitz & Appl 1997, hereafter PA97). It is impossible to obtain a solution if all the terms of the shear tensor is considered. In PA97, the authors proposed an approximate form of shear tensor by neglecting all derivatives of v and then presented a limited class of solutions. Chattopadhyay & Kumar (2016) used the same form of viscosity but obtained the full range of solutions. They also computed the massloss from such advective accretion solutions. We envisage that the method to obtain viscous solution is not easy since the Bernoulli parameter of viscous flow has no analytical form. In addition, the sonic point is not known a priori and it needs to be obtained as a part of the eigenvalue of the solution. Moreover, the angular momentum on the horizon needs to be computed. And yet the solutions obtained are limited because the viscosity is still phenomenological and various terms of the relativistic version of the shear tensor has to be neglected in order to obtain a solution. The twotemperature regime further complicates the problem as pointed out above. Most of the works from the literature assume a Newtonian form of viscosity or the Shakura & Sunyaev (1973)αviscosity prescription (hereafter SS), which would not be appropriate for our GR model. So we avoided the use of any form of viscosity since the prime focus of this paper is to present a novel methodology for obtaining unique transonic twotemperature solutions for accretion discs around BHs. In addition, it is extensively shown in Figs. 2h–i; 3h–i of Chattopadhyay & Kumar (2016) that the specific angular momentum (λ = −u_{ϕ}/u_{t}) and bulk angular momentum (L = hu_{ϕ}) in the last few 100r_{g} is almost constant and subKeplerian. This is to be expected as gravity supersedes all other interactions near the BH horizon. In order to exhibit the qualitative effect of viscosity as a representative case, the onetemperature, viscous accretion disc solutions are presented in Appendix A, following the methodology of Chattopadhyay & Kumar (2016). We use two forms of viscosity, where the viscous stress tensor is given by (i) t_{rϕ} = −2η_{vis}σ_{rϕ}, (abbreviated as PA) and (ii) t_{rϕ} = −α_{vis}p (SS form of viscosity). The form of σ_{rϕ} of PA is adopted from Peitz & Appl (1997), Chattopadhyay & Kumar (2016), but presently, we assume the dynamical viscosity coefficient η_{vis} = ρhν_{vis}, instead of η_{vis} = ρν_{vis}, where ν_{vis} is the kinematic viscosity. In Fig. 14a, we show that for both the form of viscosities (PA and SS), λ ≈ constant for r ≲ 1000r_{g}. In Fig. 14b, we plot the heat dissipated by various processes. The PA form of viscosity is stronger than the SS type of viscosity, but it is generally much weaker than Q_{B}. It is comparable to Q_{B} only in a very narrow region. Close to the horizon, the angular momentum variation of the flow is quite small and the viscous heat dissipated is less than the magnetic heating; as a result, for simplicity, in this paper we study accretion in the weak viscosity limit. We concentrate on obtaining selfconsistent twotemperature, transonic, rotating accretion solutions by considering lowangularmomentum flows at the outer boundary and compute the spectra for such flows. We study how the accretion rate, angular momentum, and mass of the central black hole might affect the solution, as well as the emergent spectra that come from such solutions.
2.3.2. Coulomb coupling
As discussed above, Coulomb coupling (Q_{ep}) serves as an energy exchange process between protons and electrons. Therefore,
The expression for Coulomb coupling in cgs units (ergs cm^{−3}s^{−1}) is given by Stepney & Guilbert (1983),
where σ_{T} is the Thomson scattering crosssection, K_{i}(x)’s are the modified Bessel functions of second kind and ith order and ln Λ_{c} is the Coulomb logarithm which is set to be equal to 20.
Apprehensions have been expressed about the possibility that more efficient energy exchange processes might exist between ions and electrons, in addition to Coulomb coupling. In that case, the accretion flow may settle down into a single temperature distribution (Phinney 1981). Begelman & Chiueh (1988) used plasma waves and Sharma et al. (2007) used magnetorotational instability to increase the energy exchange between the two species inside the flow. However, some authors have raised doubts about the effectiveness of these processes (Blaes 2014; Abramowicz & Fragile 2013). In this paper, we ignore any type of collective effects and consider only Coulomb coupling as the main energy exchange process between the protons and electrons.
2.3.3. Inverse bremsstrahlung
Inverse bremsstrahlung () is a radiative loss term for the protons, the expression of which is given as (Boldt & Serlemitsos 1969),
2.3.4. Radiative mechanisms leading to cooling of electrons
The cooling of electrons could be caused by three basic cooling mechanisms (1) bremsstrahlung (Q_{br}), (2) synchrotron (Q_{syn}), and (3) inverseComptonisation (Q_{ic}). Emissivity due to bremsstrahlung (in ergs cm^{−3} s^{−1}) is given by Novikov et al. (1973),
We used thermal synchrotron radiation in our model, following the prescription of Wardziński & Zdziarski (2000). The emissivity is given by:
where, ν_{t} is the turnover frequency, above which the plasma is optically thin to synchrotron radiation and below which it is highly selfabsorbed by the electrons itself. For the calculation of ν_{t}, we need the information of magnetic field in the flow. For this purpose, we consider a stochastic magnetic field which is in partial or total equipartition with the gas pressure, as mentioned in Sect. 2.3.1. Thus, . We set β = 0.01 throughout this paper unless otherwise specified. In Sect. 4.5, we vary the value of β and discuss how the spectrum depends on it.
The soft photons generated through thermal synchrotron process could be inverseComptonised by the electrons present in the plasma. This is given by (Wardziński & Zdziarski 2000),
where, ζ is the enhancement factor. It is expressed as
Here, Γ_{inc} is the incomplete gamma function, is the spectral index which can be defined as:
It is the slope of the powerlaw photons that are generated due to inverseComptonisation at each radius. Therefore, the net spectral index (α) of the final inverseCompton spectrum is obtained from the contributions of all the values of α_{0} from each radius of the disc. In Eq. (26), , is the average amplification factor in energy of photon per scattering and P_{sc} = 1 − exp(−τ_{es}) is the probability that a photon is scattered. Here, τ_{es} is defined as the optical depth of the medium where electron scattering is important, the expression of which is given by (Turolla et al. 1986),
2.3.5. Compton heating
As discussed above, less energetic photons would cool the flow through the process of inverseComptonisation. However if the temperature of the electrons is less than the temperature of the photons present in the flow, then the electrons will gain energy via Compton scattering. This would lead to the Compton heating of the electrons. We assume that this has the same expression as that of inverseComptonisation, but the sign changes (Esin 1997), causing heating rather than cooling. Therefore,
2.3.6. Final expressions for ΔQ_{p} and ΔQ_{e}
So, to conclude, we have taken for heating and cooling of protons:
respectively. And for heating and cooling of electrons:
respectively.
Therefore, ΔQ_{p} = δQ_{B} − Q_{ep} − Q_{ib}, and ΔQ_{e} = (1 − δ)Q_{B} + Q_{ep} + Q_{comp} − Q_{br} − Q_{syn} − Q_{ic}.
In this paper, we ignore pion production and its contribution to the observed spectra, in addition to ignoring pair production arising from the interactions of high energy photons present inside the disc. Further on, in Sect. 4, based on posteriori calculations, we show that the contribution of both processes is not significant.
2.4. Entropy accretion rate expression
If we switch off the explicit heating and cooling of protons and electrons, the gradient of proton and electron temperatures becomes (using Eq. (5)):
Due to the presence of the Coulomb interaction term, we cannot integrate the above equation and obtain an analytical form^{1}. Hence, we cannot have a measure of entropy at every point of the flow.
However, an analytical expression is admissible only in regions where Q_{ep} is negligible. Such a region is near the horizon (r_{in}), where gravity overpowers any other interaction. The integrated form of Eqs. (31) and (32) are:
where κ_{1} and κ_{2} are the integration constants which are measures of entropy. The neutrality of the plasma implies n_{e in} = n_{p in} = n_{in}. The subscript “in” indicates quantities measured just outside the horizon. Therefore,
Thus, we can write
where,
The expression of the entropy accretion rate, using Eq. (3), can be written as
2.5. Sonic point conditions
The mathematical form of Eq. (16) suggests that at some point of the flow, where a = v, the denominator (𝒟) goes to 0. Then, for the flow to be continuous, the numerator (𝒩) also has to go to 0. This is called the sonic point of the flow. Sonic points exist whenever dv/dr = 𝒩/𝒟 = 0/0. Thus, the sonic point conditions are:
and
Here, the subscript “c” corresponds to the value of flow variables at the sonic point. The derivative at the sonic point dv/dr_{c}, is computed using the L’Hospital rule.
2.6. Shock conditions
The relativistic shock conditions or the RankineHugoniot conditions (Taub 1948) are:
Conservation of mass flux across the shock: [Ṁ] = 0
Conservation of energy flux: [Ė] = 0
Conservation of momentum flux: [,
where, Σ = 2ρH and W = 2pH are the vertically averaged density and pressure, respectively. The square brackets denote the difference of the quantities across the shock.
2.7. Observed spectrum
In this model, we include radiative processes such as bremsstrahlung, synchrotron, and inverseComptonisation, which give rise to emissions spanning the whole electromagnetic spectrum. This emission (measured in units of ergs s^{−1} Hz^{−1}) when plotted as a function of frequency (in units of Hz) gives us the spectrum. The spectrum is an observational tool that helps us in determining the intrinsic properties of any object (distant or nearby). Thus, the calculation of the correct spectrum is important. While obtaining a solution for a given set of flow parameters, that is, the mass of the BH, accretion rate etc, we have information on the emission coming from each radius; or, in other words, the spectrum at each radius is known and is computed as a function of the local v, ρ, and T. When the contribution from each radius is added, we get the total observed spectrum. The model presented in this paper is in the pure GR regime, so we take into account all the general and special relativistic effects in the observed spectrum. Below, we explain the methodology used to compute the spectrum.
Let us assume that the isotropic emissivity per frequency interval per unit solid angle in the fluid rest frame is j_{ν}. If we transform this emissivity to a local flat frame, then by using specialrelativistic transformations this becomes:
Here, θ′ is the angle which the velocity of the fluid element directed inwards makes with the line of sight.
We note that all the photons emerging from the disc do not need to reach the observer. Some would be captured by the BH due to its extreme gravity. The amount of emission captured by the BH, depends on its distance from the BH. The expression to calculate this was given by Zel’dovich & Novikov (1971):
where θ^{*} is the angle within which photons will be captured by the BH and hence lost.
Now if we integrate the emissivity expression over the whole volume of the disc and on all solid angles, taking into account θ^{*}, we get the luminosity of the system as a function of frequency and, hence, the spectrum. We also account for the gravitational redshift, which introduces a factor of into the observed frequency. Furthermore, if we want to calculate the bolometric luminosity of the system, we need to integrate the frequencydependent luminosity over all the frequencies. For more details on calculating the spectrum, see Shapiro (1973). In the total spectrum, there are signatures of all the emission processes; this is discussed extensively in the results section. Bremsstrahlung emission always comes in the highfrequency part of the spectrum. When the accretion rate of the system is low, the contribution from bremsstrahlung emission is visible in the spectrum. Synchrotron emission is characterised by the turnover/absorption frequency (ν_{t}). InverseComptonisation, on the other hand, is identified as a powerlaw part in the spectrum, following the relation F_{ν} ∝ ν^{α}, α as the spectral index.
3. Solution procedure
Accretion discs around BHs are transonic in nature and may possess multiple sonic points (LT80; Fukue 1987; Chakrabarti 1989). The nature of the sonic point is also dictated by the slope of the solution at the sonic point. If the slope (i.e. dM/dr_{c}, here M = v/a is the Mach number) admits two real roots at the sonic point, then a solution can actually pass through it. These types of sonic points are termed Xtype or saddletype. The accretion solution corresponds to the negative slope, while the excretion solution corresponds to the positive slope. However, if both the roots of the slope are imaginary or complex, then matter cannot pass through them. These sonic points are called Otype (imaginary slope) and spiraltype (complex), respectively. The combined effect of the flow parameters like E, λ, & Ṁ, determines the number of sonic points formed inside an accretion flow, as well as the topology of the solution. An accretion flow with low values of λ admits only one outer sonic point (r_{co}; located at larger distance from the BH), while those with higher values of λ admit only inner sonic points (r_{ci}; closer to the BH). In the intermediate λ range, the accretion disc may admit a maximum of three sonic points : inner (r_{ci}), middle (r_{cm}), and outer (r_{co}), where r_{ci} and r_{co} are Xtype, while the middle sonic may be spiral or Otype depending on whether the system is dissipative or nondissipative, respectively (for further information on sonic points see, Holzer 1977; LT80; Ferrari et al. 1985; Fukue 1987; Chakrabarti 1989). Flows with high value of E, generally admit only one r_{ci}. On the other hand, Ṁ controls radiative cooling and modifies the thermal state of the flow. This, in turn, modifies the range of E and λ, which allows for the formation of multiple sonic points.
In Sect. 3.1, we describe the method for finding sonic points. In Sect. 3.2, we show that the transonic solutions are degenerate and we present an extensive discussion of how to remove the degeneracy.
3.1. Method for obtaining sonic points in twotemperature flows
In a single temperature regime, the location of the sonic point and its property is unique for a given set of constants of motion. In dissipative systems, sonic points are not known a priori and they are obtained selfconsistently by integrating the equations of motion. So, presently, in the two temperature regime, we follow exactly the same procedure to solve the equations, as done in the single temperature realm. We need to select some fixed boundary from where we can start integrating dv/dr, dΘ_{p}/dr and dΘ_{e}/dr to find the sonic point. As r → 2r_{g}: v → c and E is defined at the horizon, the latter being a constant of motion. However, we cannot start integration from r = 2r_{g} because of coordinate singularity on the horizon. Therefore, we select a point asymptotically close to the horizon (r_{in} = 2.001r_{g}). It is here, where gravity overpowers any other processes or interactions, therefore, infall timescales are much smaller than any other timescales. In other words, at r = r_{in}, X_{f} → 0 and E → ℰ = −hu_{t} (from Eq. (6)). Simplifying this, we obtain an expression of v_{in} in terms of E, λ, r_{in}, Θ_{pin}, and Θ_{ein}. We list the procedure for obtaining a transonic solution below:

For a given set of values of E, λ & Ṁ, we start integration from a point asymptotically close to the horizon at r_{in} = 2.001 (in units of r_{g}).

As r_{in} → 2, E → ℰ = −h_{in}u_{t} = h_{in}(1 − 2/r_{in})^{1/2}γ_{v}γ_{ϕ}. Here h_{in} = h_{in}(Θ_{pin}, Θ_{ein}) is the specific enthalpy at r_{in}.

We supply Θ_{pin}.

We also supply Θ_{ein}. Then v_{in} is obtained as a function of the flow parameters and can be expressed as

Using the values of Θ_{pin}, Θ_{ein}, and v_{in} we integrate Eqs. (14)–(16), from r = r_{in} to outwards. As we integrate, we simultaneously check the sonic point conditions (Eqs. (38) and (39)).

If a sonic point is not found, we supply another value of Θ_{ein} and repeat steps 4 and 5 until the sonic point conditions are satisfied.

Once a sonic point is found, we integrate the equations of motion from sonic point to larger distances and obtain the full, global, transonic twotemperature accretion solution.

After we locate one sonic point, we change Θ_{ein} again and repeat steps 4–7 to check if any other sonic point exists. If it is found, we obtain its corresponding transonic solution.
We note that if our supplied guess value of Θ_{pin} is unphysical, then even by iterating with all possible values of Θ_{ein}, the sonic point conditions can never be satisfied. This is basically the modified version of the methodology adopted by Le & Becker (2005), who obtained accretion solutions for a dissipative flow in the single temperature regime.
We illustrate the procedure to find transonic solutions, enlisted above, in Figs. 1a–c. All three panels in this figure plots Mach number (M = v/a) vs log r. The accretion disc parameters are λ = 2.5, E = 1.000045, and Ṁ = 0.001 Ṁ_{Edd} around a BH of 10 M_{⊙}. We would like to point out that all these accretion disc parameters and the BH mass chosen are for representative purposes only. These have been varied and their effect on the solution were studied later.
Fig. 1. Method for finding sonic points. Solutions are presented in terms of Mach number M (=v/a) vs log r plot. Θ_{pin} = 7.162 × 10^{−2} for all iterations. Panel a: iterations to obtain inner sonic point r_{ci} (black circle) and panel b: iterations to obtain outer sonic point r_{co} (black star). Various branches plotted are multivalued branches (MB; green dasheddot), transonic (TS; red dashed), and supersonic (SB; blue dasheddot). Respective Θ_{ein}s are mentioned inside the panels. Panel c: full set of transonic solutions. Global accretion solution (red solid) through r_{co} and accretion solution through r_{ci} (red, dashed) is plotted. Equatorial global wind (through r_{co}) and nonglobal wind (through r_{ci}) are represented using red dotted curve. The accretion disc flow parameters used are λ = 2.5, E = 1.000045, Ṁ = 0.001 Ṁ_{Edd}, and M_{BH} = 10 M_{⊙}. 
In Fig. 1a, we present the method to obtain inner sonic point or r_{ci} and the transonic solution through it. Following steps 1, 2, and 3, we supply Θ_{pin} = 7.162 × 10^{−2} for the aforementioned values of E, λ, and Ṁ. Following step 4, we start by supplying a high value of Θ_{ein} and obtain v_{in}. Then we integrate the equations of motion (step 5). For higher values of Θ_{ein}, we obtain a multivalued branch (MB) of solutions. We plot one such MB solution (green dasheddot) corresponding to Θ_{ein} = 11.948. Clearly, the MB solutions are not correct. We reduce the value of Θ_{ein} (step 6) and repeat the whole procedure up to step 5. We observe that as we reduce Θ_{ein}, the MB solutions approach the transonic solution, that is, they would shift rightward, but in all probability, we would overshoot the transonic solution and end up with a purely supersonic branch (SB) solution (i.e. when v > a or M > 1 at all r). We plot a representative case of a purely SB solution (blue dasheddot) corresponding to Θ_{ein} = 7.731. When the solutions corresponding to various Θ_{ein} suddenly shift from an MB solution to SB, then we know that the Θ_{ein} corresponding to a transonic solution (TS) lies in between these two values of Θ_{ein}. We iterate using the electron temperature at r_{in} within the range 7.731 < Θ_{ein} < 11.948 and obtain the transonic solution (TS) (red dashed) and the sonic point is at r_{ci} = 5.186 (black circle). Then, by following step 7, we obtain the complete transonic solution from r_{in} to a large distance through r_{ci}. In Fig. 1b, we present the procedure for finding the existence of the outer sonic point for the same set of flow parameters. Following step 8, we reduce Θ_{ein} by a comparatively large value, such that we start obtaining MB type solutions similar to the ones we obtained while trying to locate r_{ci}. In this panel, we present an example of an MB solution (green dasheddot) for Θ_{ein} = 6.571. We repeat steps 4–6 and check when the solution jumps from MB to SB. In this instance, it corresponds to Θ_{ein} = 7.106 (blue dasheddot). We iterate on Θ_{ein} between the two limits, 6.571 < Θ_{ein} < 7.106, until we obtain a transonic solution through r_{co} = 3644.9 (black star). No other sonic point exists for these flow parameters (E, λ, Ṁ). In Fig. 1c, we plot the complete set of transonic solutions, accretion (red dashed and red solid) as well as equatorial wind solutions (red dotted) for the disc parameters λ = 2.5, E = 1.000045 and Ṁ = 0.001 Ṁ_{Edd} around a BH of 10 M_{⊙}. The global accretion solution connecting infinity to the horizon is represented by a red solid line, while a red dashed line represents an accretion solution which is not global. We should remember that not all set of disc parameters produce multiple sonic points and this point is discussed in detail in the following sections.
3.2. Presence of degeneracy in twotemperature transonic solutions : Method for removing it to obtain unique transonic solutions, invoking the second law of thermodynamics
In the last section, we laid down the procedure to obtain transonic solution by supplying a guess value of Θ_{pin} and iterating with various values of Θ_{ein}, until we get the sonic point. This has been elaborately discussed in steps 1–8 in Sect. 3.1. Now, if we choose a different value of Θ_{pin} and again follow the steps 1–8 of Sect. 3.1 for the same set of disc parameters (E, λ & Ṁ), we obtain a different transonic solution with distinctly different sonic point properties. This means that for a given set of constants of motion, twotemperature EoMs admit multiple transonic solutions. Since the number of EoMs (accretion rate, momentum, and energy equation) are less than the number of flow variables (density, velocity components, and two temperatures), even the transonic solutions become degenerate as a result. From the second law of thermodynamics, it is clear that out of all the possible solutions, only the solution with the highest entropy () would be expected to be favoured by nature (Sarkar & Chattopadhyay 2019). In dissipative systems, entropy is not conserved, so we measure entropy in the region near the event horizon (see, Eq. (37) of Sect. 2.4). For each transonic solution corresponding to a given Θ_{pin}, we compute . If the computed is plotted with respect to Θ_{pin}, then there is a clear maxima in . Following the second law of thermodynamics, the solution (corresponding to a Θ_{pin}) with maximum entropy () is the physically plausible solution. Hence, we are able to constrain the degeneracy and obtain a unique transonic twotemperature solution for a given set of constants of motion.
In Figs. 2a–g, we illustrate the methodology for obtaining a unique twotemperature transonic solution. We choose the accretion disc flow parameters: E = 1.0015, λ = 2.6, Ṁ = 0.02 Ṁ_{Edd}, and M_{BH} = 10 M_{⊙}. As described above, we supply a Θ_{pin} or equivalently T_{pin} and then we iterate on Θ_{ein} (or T_{ein}) to obtain a transonic solution. The entropy accretion rate corresponding to the transonic solution for the particular T_{pin}, is plotted in Fig. 2g. The black solid line is for those T_{pin} whose solution passes through outer sonic point (r_{co}) and black dotted line is for those passing through inner sonic point (r_{ci}). From this figure, we select few T_{pin}s (points marked “a”–“f”) which are (a) 3.1 × 10^{11} K, (b) 5.605 × 10^{11} K, (c) 6.04 × 10^{11} K, (d) 6.460 × 10^{11} K, (e) 6.554 × 10^{11} K and (f) 7.0 × 10^{11} K and plot their solutions in terms of M vs log r, presented in Figs. 2a–f, respectively. The global accretion solutions are represented by solid curves, dotted are the wind types, while the dashed curves are accretion solutions which are not global. There is a range of T_{pin} where both inner and outer sonic points are present and is the multiple sonic point regime. For example, “b”, “c” and “d” has both inner (marked circle) and outer sonic points (marked triangle). For point “e”, the solutions passing through inner and outer sonic points (magenta triangle and magenta circle), have almost the same value of (see inset). A corresponding solution is plotted in Fig. 2e, which shows that the global solution (connecting horizon to large distances) passes through r_{co}. Now, corresponding to point “a” (orange triangle, in panel g), the solution passes only through r_{co} (Fig. 2a), while for point “f” (blue circle), the solution passes through r_{ci} (Fig. 2f). However, the entropy is exactly same for both these points. This means that the solution character can be completely different even if has the same value. Similarly another pair, “b” and “d”, also has the same entropy but different T_{pin}. Also, their corresponding solutions are significantly different, although both solutions lie in the multiple sonic point regime. This shows that not only all solutions presented in the figure have the same E, Ṁ, λ, but they may even have the same and yet, the solutions are completely different.
Fig. 2. Left: M vs log r plot for various values of T_{pin}. (a) T_{pin} = 3.1 × 10^{11} K, (b) T_{pin} = 5.605 × 10^{11}K, (c) T_{pin} = 6.04 × 10^{11} K, (d) T_{pin} = 6.460 × 10^{11} K, (e) T_{pin} = 6.554 × 10^{11} K and (f) T_{pin} = 7.0 × 10^{11} K. Global solutions are represented by solid lines. Panel g: vs T_{pin}. Solid black curve is for the solutions passing through outer sonic point, while dotted black curve is for solutions passing through inner sonic point. Panels a–f: solutions corresponding to the points marked in right panel g. The disc flow parameters are E = 1.0015, λ = 2.6 and Ṁ = 0.02 Ṁ_{Edd}. The spacetime is described by a BH of mass 10 M_{⊙}. 
In order to drive the point all the way home, we present , & L (luminosity) for all degenerate solutions in Table 1. Some solutions can be about four times more luminous than other solutions. It is evident from the figure and table that degeneracy in twotemperature model is a serious problem. An observational parameter such as L is quite different for different degenerate solutions. Any random choice from the pool of degenerate solutions would provide us with completely wrong information about the system as well as an erroneous spectrum. Hence, the removal of degeneracy is important. Although all the solutions presented above have the same energy, angular momentum, and accretion rate, only one of them possesses the highest entropy (highest ). In this particular case, the highest entropy solution is the one corresponding to point “c” (red triangle) in –T_{pin} curve (Fig. 2g), and the correct, unique twotemperature accretion solution is represented in Fig. 2c (red solid).
3.3. The stability of the highest entropy transonic solutions
In this section, we investigate the stability of the unique transonic twotemperature solution selected from the available set of degenerate solutions. We note that the proposed unique solution is that of the highest entropy and, based on the second law of thermodynamics, nature should prefer it. Therefore, the solution should be stable. However, because there is a degeneracy of solutions, we need to study the stability of the proposed unique solution as well. We provide a qualitative analysis for the stability of the unique twotemperature solution.
Let us assign and T_{pinmax} = T_{pin} for which is maximum. Let us further define δT_{pin} which is the difference between adjacent higher and lower T_{pin} and as the gradient of . Now the change in is given by
It is clear that Δ = 0 for any extrema of curve, but the solution is said to be stable if T_{pin} moves away from the value T_{pinmax} and the system adjusts automatically to regain its old value. We prefer the graphical method to investigate the stability and the technique is similar to the first derivative test for obtaining local extrema (Melo 2014). We plot Δ vs T_{pin} in Fig. 3. At T_{pin} < T_{pinmax}, the figure shows that Δ > 0. Then from Eq. (42), . So the system tends to go to higher , which we denote using a rightward arrow. Similarly, for T_{pin} > T_{pinmax}, Δ < 0 which implies . Since, based on the second law of thermodynamics, a physical system is not amenable to decreasing its entropy, therefore for T_{pin} > T_{pinmax}, the system would tend to come back to T_{pinmax}. We represent this in the figure by using a leftward arrow. In other words, entropy can only increase (), if T_{pin} → T_{pinmax} from either side of T_{pinmax}. Hence, a solution corresponding to T_{pin} = T_{pinmax} must be stable. Therefore, in addition to the fact that T_{pinmax} corresponds to a solution with maximum entropy, we can conclude that the solution is also stable.
Fig. 3. Stability analysis of the unique transonic twotemperature solution with maximum entropy. The flow parameters used are same as Fig. 2. is plotted against variation of T_{pin}. The arrows indicate that Δ converge at T_{pin} = T_{pinmax} (blue dot) and is the stable equilibrium solution. This T_{pin} is the solution with maximum entropy marked “c” in Fig. 2g. 
4. Results
Twotemperature accretion solutions are parameterised by E, λ, Ṁ. In addition, β_{d} and β controls heating and cooling. Since we quote accretion rates in terms of the Eddington rate, therefore, the information on M_{BH} also enters the solution. In this section, we study in detail the twotemperature accretion solutions and discuss their spectral properties. We only analyse solutions with maximum entropy, selected from the available degenerate group of solutions.
4.1. General twotemperature solutions
In Fig. 4, we study a typical twotemperature transonic advective accretion disc solution. The parameters used are, E = 1.000045, λ = 2.5, Ṁ = 0.001 Ṁ_{Edd}, and M_{BH} = 10 M_{⊙}. In Fig. 4a, we plot vs T_{pin}. The solution with T_{pin} = 6.0 × 10^{11} K has the maximum entropy (marked with green triangle) and the corresponding M vs log r (green solid line) is plotted in Fig. 4b. The global solution passes through an outer sonic point whose position is r_{co} = 3040.182 (black star). The radial threevelocity (v; cyan solid) in corotating frame and flow velocity in the azimuthal direction (v_{ϕ}; blue dotted) is plotted in Fig. 4c. Matter far away from the horizon has negligible velocity in the radial as well as in the azimuthal directions. As it approaches the BH (r → r_{g}), v → c, thus satisfying the BH boundary condition. On the other hand, v_{ϕ} increases with the decrease of r, but is maximised at r ∼ 3r_{g} and finally goes to zero on the horizon. This is mainly because near the horizon infall timescale is much shorter than any other timescales. The strong gravity does not allow the matter enough time to rotate in the azimuthal direction. In Fig. 4d, we plot the number density (in units of cm^{−3}) as a function of radius. The number density increases with the decrease in radius, as it should be for a convergent flow. Values for T_{p} (red dotted) and T_{e} (orange solid) are plotted in Fig. 4e. The cooling processes are dominated by electrons compared to protons. They are, however, coupled by a Coulomb coupling term which acts as an energy exchange term between the protons and electrons, as previously discussed. This term is weak, which allows protons and electrons to equilibrate into two different temperatures (T_{p} and T_{e}), unlike in the case of singletemperature flows where Coulomb coupling term is assumed to be very efficient, allowing protons and electrons to attain a single temperature. Figure 4f shows that adiabatic indices of both protons (red dotted) and electrons (orange solid) vary with the flow. This justifies our use of CR EoS. Γ_{p} ∼ 1.66 and Γ_{e} ∼ 1.60 at large distances away from the BH, hence, both the species are thermally nonrelativistic. When the flow approaches the BH, Γ_{e} becomes relativistic i.e. Γ_{e} ∼ 1.33 near the horizon. We can see that Γ_{p} does not vary much but it becomes mildlyrelativistic near the horizon owing to the higher mass of protons. In Fig. 4g, we prove that the generalised Bernoulli constant is a constant of motion throughout the flow, even in the presence of dissipation.
Fig. 4. (a) plotted against T_{pin}. The entropy for inner sonic point solutions (black, dotted) and outer sonic points (black, solid) are presented. T_{pin} marked with green triangle corresponds to maximum entropy solution. Flow variables plotted are (b) M (green, solid), (c) v (cyan, solid) and v_{ϕ} (blue, dotted), (d) log n (magenta, solid), (e) T_{p} (red, dotted), T_{e} (orange, solid), (f) Γ_{p} (red, dotted), Γ_{e} (orange, solid), (g) E (brown, solid) as functions of log r. The flow parameters are E = 1.000045, λ = 2.5, Ṁ = 0.001 Ṁ_{Edd} and M_{BH} = 10 M_{⊙}. Panel b: the sonic point is marked with a black star. 
4.1.1. Emissivities and spectral properties
In Fig. 5, we present the heating and cooling rates for the solution plotted in Fig. 4. In Fig. 5a, we plot the heating terms. The red solid line represents the heating due to magnetic dissipation. This amount of heat is assumed to be equally distributed among protons and electrons. Blue dotted line represents the Compton heating of electrons. It is mainly the hard bremsstrahlung photons present inside the flow which leads to the heating of the electrons, which is due to their higher energy as compared to electrons. In Fig. 5b, we plot Coulomb Coupling term (yellow solid line) and inverse bremsstrahlung (green dotted line). The strongest heating term is represented by Q_{B}.
Fig. 5. Emissivity vs log r plot for the flow presented in Fig. 4, shown here in the top three panels. Bottom panel: spectrum of the accretion flow. 
In Fig. 5c, we plot the emissivities of all the cooling processes for electrons: bremsstrahlung (red dotted line), synchrotron (green dashed line), and inverseComptonisation (blue dasheddotted line). For the present set of disc parameters, at the outer boundary of the disc, the temperatures are nonrelativistic and, therefore, bremsstrahlung emission dominates over all other processes. In the inner regions of the accretion disc, that is, close to the BH horizon, synchrotron and inverseComptonisation becomes important and exceeds bremsstrahlung. However, inverseComptonisation is less than synchrotron emission, mainly due to the low accretion rate of the flow. The total cooling of electrons is represented by a black solid line. In Fig. 5d, we plot the spectrum for the accretion flow (black solid line). It is plotted by summing up the contributions of all emission processes at each radius. General and special relativistic frame transformations from the fluid rest frame to the observer frame are taken into account while computing the spectra, including photon capture and photon bending effect due to the presence of strong gravity. This is widely discussed in Sect. 2.7. Spectrum of each emission process is also plotted. Bremsstrahlung is shown as a red dotted line, synchrotron by a green dashed line, and inverseComptonisation by a blue dasheddotted line. The overall luminosity of the system is low, namely, L = 2.536 × 10^{29} ergs s^{−1}, with an radiative efficiency of η_{r} = 1.957 × 10^{−5}. It may be noted that efficiency is defined as η_{r} = L/(Ṁc^{2}). The spectral index is α = 1.744.
4.2. Contributions of different regions of the accretion disc to the overall spectrum
In Fig. 5d, we show the contribution of all the emission processes in the total broad band continuum spectrum. However, in the following, we aim to investigate the contribution of various regions of an accretion disc in the overall spectrum. We chose a set of flow parameters: E = 1.0002, λ = 2.48, Ṁ = 0.05 Ṁ_{Edd}, and M_{BH} = 10 M_{⊙}. In Fig. 6a, we plot the Mach number M of the accretion flow and in Fig. 6b, we plot T_{p} (black solid) and T_{e} (black dotted) as a function of log r. In panels a and b, we indicate various regions with vertical lines which represents accretion disc section from 2 − 3r_{g} (magenta, dotted), 3 − 5r_{g} (blue, dashed), 5 − 8r_{g} (green, shortdasheddotted), 8 − 10r_{g} (brown, longdasheddotted), 10 − 100r_{g} (orange, dasheddoubledotted) and 100 − 1000r_{g} (red, dashedtripledotted). The spectra from all these regions are separately over plotted in Fig. 6c, the colourcoding of the spectra matches the region from which they are computed. The black curve represents the overall spectrum for the disc parameters stated above. The contribution from the region r = 10^{3} − 10^{4}r_{g} is too low in the overall spectrum and is not plotted, therefore, in order to avoid cluttering the figure. The spectrum computed from the region 2 − 3r_{g} is low in spite of the high values of n and T_{e} since a significant number of photons emitted from that region are captured by the BH. Most of the high energy emission is contributed by accreting matter from the region between 3 − 5r_{g} and 5 − 8r_{g} and the lowenergy end of the spectra from this region is always around and above 10^{12} Hz. We tabulated these details and other spectral properties in Table 2. We can conclude from the table that ∼90% of the emission comes from a region < 10r_{g} of the accretion disc. The lower energy part of the spectrum is mostly contributed by the outer part of the disc. Since we only considered the advective disc, the spectrum is hard and the radiative efficiency for this particular set of disc parameters is less than 1%.
Fig. 6. (a) M and (b) log T vs log r, and (c) total spectrum (black solid) plus contribution from various length scales of the accretion disc, 2 − 3r_{g} (magenta, dotted), 3 − 5r_{g} (blue, dashed), 5 − 8r_{g} (green, shortdasheddotted), 8 − 10r_{g} (brown, longdasheddotted), 10 − 100r_{g} (orange, dasheddoubledotted) and 100 − 1000r_{g} (red, dashedtripledotted). Flow parameters are E = 1.0002, λ = 2.48, Ṁ = 0.05 Ṁ_{Edd}, and M_{BH} = 10 M_{⊙}. 
4.3. Dependence of accretion solutions and corresponding spectra on energy and angular momentum
In Figs. 7 and 8, we investigate the dependence of accretion solutions and the corresponding spectra on E and λ for a 10 M_{⊙} BH with Ṁ = 0.01 Ṁ_{Edd}. In these figures, E increases from left to right and the values are E = 1.0005, 1.001, 1.003, and 1.01; whereas, as we go from top to bottom, λ increases as λ = 2.40, 2.55, 2.70, and 2.85. In short, E changes along the row while λ changes along the column. Lowangularmomentum flows (λ = 2.40) behave as Bondi flow, possessing a single sonic point, through which the global solution passes (see Figs. 7a1–a4), irrespective of the value of E. As angular momentum increases, the rotation head of the specific energy (E) of the flow plays a significant role inside the system and multiple sonic points form in an appreciable section of the parameter space. For λ = 2.55 (Figs. 7b1–b3), multiple sonic point exists in a large range of E. In Figs. 7b1–b2, the global solution (blue, solid) passes through the outer sonic point; whereas in Fig. 7b3, the solution (blue, solid) harbours a shock and passes through both inner and outer sonic points. In Fig. 7b4, only a single sonic point exists. This is mainly due to the fact that for flows with higher energy, the distribution of sound speed a is generally higher compared to flows with lower E. Therefore, the flow can only become transonic when v increases significantly, which can happen only very close to the BH. For low values of E, the sound speed distribution a(r) is comparatively low. Therefore the sonic points form further out because the flow becomes transonic whenever the infall velocity v(r) attains moderately high values.
Fig. 7. Variation of solutions, M as a function of log r with variation of E and λ. From left to right: the specific energy increases as E = 1.0005, 1.001, 1.003, and 1.01. From top to bottom: the angular momentum increases as λ = 2.40, 2.55, 2.70, and 2.85. The other parameters are Ṁ = 0.01 Ṁ_{Edd} and M_{BH} = 10 M_{⊙}. 
Fig. 8. Variation of the spectrum with E and λ. The set of values for E and λ and other parameters are same as in Fig. 7. 
Angular momentum has a different effect on the flow structure. If we increase λ, then the distribution of v_{ϕ}(r) increases. Higher values of v_{ϕ} restrict the increase of v to moderate values, except near the horizon. So for flows with higher λ, the sonic points shift towards the BH. For even higher λ ≥ 2.85, only inner sonic point exists irrespective of the value of E (see Figs. 7d1–d4). In Figs. 7c1–c2 which are for λ = 2.70, shocks form even at low energies. If one compares with single temperature accretion discs (Chattopadhyay & Chakrabarti 2011; Kumar & Chattopadhyay 2014, 2017; Chattopadhyay & Kumar 2016), it is clear that multiple sonic points form in a much smaller range of the energyangular momentum parameter space of a twotemperature accretion disc, as is shown in Figs. 7a1–d4.
Figure 8 shows the corresponding spectra, which span from 10^{12} − 10^{22} Hz. As a general trend, with the increase in λ of the system, luminosity increases since matter gets enough time to radiate. However the spectral shape and slope (arising because of inverseComptonisation) remains roughly the same, except for the solutions which harbours shock. The spectral slope is flatter in case of shocked solutions. With the increase in E, thermal energy of the system increases, emission is hence higher. The spectral shape and slope is visibly changed. Bremsstrahlung emission, the broad peak in the higher frequency range, is increased with the increase in E (left to right), while angular momentum seems to have little effect on this particular radiative process. Since, in this case, we are dealing with a flow with a low accretion rate, the spectrum is relatively soft as inverseComptonisation is not significant.
4.4. Shocked solution, spectra and the parameter space
In Figs. 7b3,c1,c2, the accretion solutions admit stable shocks. As previously discussed, for low λ, a flow admits only one sonic point (Figs. 7a1–a4). As λ increases, the flow possess multiple sonic points. A flow can pass through both the sonic points only when the shock conditions are satisfied (see, Sect. 2.6). With the increase in λ, the centrifugal term increases and the twin effect of the centrifugal and the thermal term can restrict the infalling matter, leading to a centrifugal pressure mediated shock transition. In the following section, we analyse the shocks present in twotemperature accretion flows.
In Fig. 9a, we present a typical shock solution for the parameters E = 1.002, λ = 2.58, Ṁ = 0.2 Ṁ_{Edd}, and M_{BH} = 10 M_{⊙}. For these parameters, the flow possess multiple sonic points. The blue dashed line is for the accretion solution passing through the outer sonic point r_{co}. When this solution becomes supersonic, it encounters a shock at r_{sh} = 20.952. Then it jumps to the subsonic branch and enters the BH supersonically after crossing through the inner sonic point r_{ci}. The global solution is represented with a red solid line. The compression ratio (, ± implies post and preshock quantities, respectively) is 1.459 in this case. In Fig. 9b, we plot the number density as a function of r. At the shock, there is an increase in number density of both protons and electrons equally. This leads to increased cooling in the system which is evident from Fig. 9c, where we plot emissivities of various cooling processes related to electrons. The corresponding spectrum is plotted in Fig. 9d. In Figs. 9c–d, bremsstrahlung is represented using dotted green line, synchrotron in dashed yellow, and inverseComptonisation in dasheddotted magenta, while the total cooling is represented by solid black line. The accretion rate of the system is high, so the spectrum is mainly dominated by inverseComptonisation, especially in the postshock region.
Fig. 9. Typical shocked solution (a), with its corresponding number density (b), emissivities, (c) and spectrum (d), is presented. The parameters taken are E = 1.002, λ = 2.58, Ṁ = 0.2 Ṁ_{Edd}, and M_{BH} = 10 M_{⊙}. 
In Fig. 9d, the total spectrum of the system is plotted in black, while superimposed upon it is the spectrum (blue solid) of the shockfree solution (blue dashed of panel a). The luminosity of the system is 2.831 × 10^{36} ergs s^{−1}, which corresponds to an efficiency (η_{r}) of 1.09%, while for a shockfree branch (blue dashed), the luminosity would be 1.293 × 10^{36} ergs s^{−1} and η_{r} = 0.50%. Because of the shock, the luminosity and, hence, the efficiency of the system doubles. However, it seems that there is no special spectral signature of shock in accretion flow, except that the luminosity of the powerlaw part of the spectrum increases. This is also evident in Figs. 8b3, c1, and c2.
The bounded region in Eλ space in Fig. 10, represents the shock parameter space, i.e. a flow with E, λ values from the bounded region for the given Ṁ that is to undergo a stable shock transition. Each bounded region or shockparameter space is characterised by different accretion rates: Ṁ = 0.01 (green, solid), Ṁ = 0.1 (blue, dashed), and Ṁ = 1.0 (red, dotted) around a 10 M_{⊙} BH. We can see that as Ṁ increases, the parameter space decreases and shifts to the lower angular momentum side. High values of Ṁ imply higher rates of cooling and, therefore, much hotter flow at the outer boundary, which can accrete and form the disc. Hence, even for lower λ, the centrifugal term in conjunction with the thermal term can resist the infall to produce an accretion shock. That is why for higher Ṁ, the shock parameter space shifts to the lower λ values. For low Ṁ, the shock parameter space is almost similar, it significantly changes only in the presence of high accretion rates. More interestingly, it is clear that an accretion flow with high accretion rate may also harbour accretion shocks.
Fig. 10. Shock parameter space for Ṁ = 0.01 Ṁ_{Edd} (green, solid), 0.10 Ṁ_{Edd} (blue, dashed), and 1.00 Ṁ_{Edd} (red, dotted) around a 10 M_{⊙} BH. 
4.5. Dependence of spectrum on β
The magnitude of stochastic magnetic field inside the flow is controlled by β. Any change to it would lead to a change in synchrotron emission from electrons and, eventually, a change in the radiation due to inverseComptonisation. Hence, the spectra that an observer would see significantly depends on the value of β. In Fig. 11, we plot the change in spectra with change in β for the flow parameters E = 1.003, λ = 2.54, Ṁ = 0.1 Ṁ_{Edd}, and M_{BH} = 10 M_{⊙}. We varied β: 0.002 (blue, solid), 0.01 (green, dashed), and 0.02 (red, dotted). The present flow has high accretion rate where cooling is more pronounced. Even for low β, the power law signature in the spectrum arising due to inverseComptonisation, is hard. However, the dominant emission comes from bremsstrahlung, as can be inferred from its bump at higher frequency regime. As we increase β, the bump feature vanishes. This is mainly due to the fact that with the increase in β, synchrotron emission and, hence, inverseComptonisation increases more as compared to bremsstrahlung, which is independent of the magnitude of the magnetic field in the flow. The synchrotron turnover frequency also shifts to higher frequencies with the increase in β.
Fig. 11. Change in spectra with increase in β = 0.002 (blue, solid), 0.01 (green, dashed), and 0.02 (red, dotted). Other parameters used are E = 1.003, λ = 2.54, and Ṁ = 0.1 Ṁ_{Edd} in an accretion disc around M_{BH} = 10 M_{⊙}. 
4.6. Dependence of solutions and spectra on β_{d}
Figure 5a shows that magnetic dissipation is a more efficient heating mechanism compared to Compton heating and to the Coulomb heating of electrons. In Figs. 12a1–b2, we vary β_{d}, which controls magnetic dissipation. We compare the solutions and resultant spectra for β_{d} = 0.013 (blue, solid), 0.015 (green, dashed), and 0.017 (red, dotted). For Figs. 12a1,a2, we chose a higher ratio between magnetic and gas pressure, β = 0.2, and accretion rate Ṁ = 1.0 Ṁ_{Edd}. For Figs. 12b1,b2, we selected Ṁ = 1.5 Ṁ_{Edd} and β = 0.15. For both the cases, we have E = 1.001, λ = 2.61 and M_{BH} = 10 M_{⊙}. The sonic point, luminosities, and spectral index of the accretion flows are given in Table 3. Evidently, luminosity and, hence, efficiency, decreases with increasing β_{d}, but increases with increasing Ṁ. If the dissipative heating is higher (i.e. higher β_{d}), then matter with lower temperature at large r may achieve the same E. Therefore, the accretion flow would have an overall lower temperature and would emit less. That is exactly, what is observed in Figs. 12a2,b2, where the luminosity goes down with the increase in β_{d}. The spectrum is decisively hard for superEddington accretion rates. It means that no single flow parameter can dictate whether the spectrum is hard or soft; instead all the flow parameters together contribute for the final outcome. However, it is clear that β_{d} does influence the emitted spectra and luminosity significantly.
Fig. 12. Accretion solutions (a1, b1) and their corresponding spectra (a2, b2) plotted for a flow with E = 1.001, λ = 2.61 around M_{BH} = 10 M_{⊙}. Various curves are for β_{d} = 0.013 (blue, solid), β_{d} = 0.015 (green, dashed) and β_{d} = 0.017 (red, dotted). The accretion rates and ratio of magnetic to gas pressure are chosen as Ṁ = 1.0 Ṁ_{Edd}, β = 0.2 (a1, a2), and Ṁ = 1.5 Ṁ_{Edd}, β = 0.15 (b1, b2). 
4.7. Possibility of pair production and pion production
Up to this point, we have investigated how various factors can affect the twotemperature solutions and resulting spectra. However, it may be noted that we have ignored pair production from particleparticle interactions in the accretion disc or from accretion disc radiations. We have also ignored the production of gammarays due to high energy interactions like pion decay. We assumed that these processes will not significantly affect the solutions. In Appendix B, we investigate the pair production processes a posteriori. We compare the number densities of protons n_{p} with positrons n_{e+} (Figs. B.1a1, b1) generated through photon interactions produced in accretion discs as well as compare the total emissivity Q_{total} with pair annihilation emissivity Q_{ann} (Figs. B.1a2, b2). We consider two sets of accretion disc parameters (1) Ṁ = 1.0, β = 0.2 (Figs. B.1a1, a2) and (2) Ṁ = 1.5, β = 0.15 (Figs. B.1b1, b2). The rest of the parameters are common for both the cases, namely, β_{d} = 0.013, E = 1.001, and λ = 2.61. These two accretion disc cases are described around a BH of 10 M_{⊙}. After the a posteriori calculations, elaborately discussed in Appendix B, we can conclude that, n_{e+} ≪ n_{p} and Q_{ann} ≪ Q_{total}.
In Appendix C, we compute the production of pions (π^{0}) a posteriori and the gamma rays emitted due to its decay. We plot log T_{p} vs log r in Figs. C.1a1, b1 and the corresponding spectra in Figs. C.1a2, b2. We study the generation of pions and gamma ray photons for two cases (1) by varying accretion rate (Ṁ = 0.01 : red, dotted, 0.1 : green, dashed, and 1.0 : blue, solid) around a BH of M_{BH} = 10 M_{⊙} (Figs. C.1a1, b1) and (2) by varying mass of the BH (M_{BH} = 10^{2} : blue, solid, 10^{4} : green, dashed, and 10^{6} : red, dotted), keeping the accretion rate Ṁ = 0.1 constant. The other disc parameters are E = 1.0007, λ = 2.61, β = 0.01, & β_{d} = 0.001. Luminosity for higher Ṁ is higher and so is the gammaray produced by decay of pions. The same trend is observed when we increase the BH mass. However, the gamma ray luminosity is always < 10^{−5} times that of the total luminosity. An elaborate discussion on these two cases is presented in Appendix C. We can conclude that the consideration of pair production or pion decay does not affect the accretion solutions or the spectra significantly.
4.8. Dependence on Ṁ and M_{BH}
In Fig. 13a, we plot the continuum spectra for Ṁ = 0.1 Ṁ_{Edd} (blue, solid), Ṁ = 0.6 Ṁ_{Edd} (green, dashed), and Ṁ = 1.2 Ṁ_{Edd} (red, dotted) from a disc around a BH of 10M_{⊙}. The disc becomes brighter as Ṁ increases and even the efficiency also increases. The spectra also becomes harder, mainly because the inverseCompton output increases with the increase in number density of hot electrons inside the flow. However, the range of frequency ν on which the spectrum is distributed do not increase appreciably with the increase in accretion rate. The corresponding spectral properties are presented in Table 4.
Fig. 13. Spectra from (a) M_{BH} = 10 M_{⊙} for different accretion rates Ṁ = 0.1 Ṁ_{Edd} (blue, solid), Ṁ = 0.6 Ṁ_{Edd} (green, dashed), and Ṁ = 1.2 Ṁ_{Edd} (red, dotted); (b) Ṁ = 0.1 Ṁ_{Edd} but around M_{BH} = 10 M_{⊙} (blue, solid), M_{BH} = 10^{3} M_{⊙} (magenta, dashed), and M_{BH} = 10^{6} M_{⊙} (brown, dotted). Other disc parameters are E = 1.001 and λ = 2.4. 
In Fig. 13b, we plot spectra from discs with the same accretion rate Ṁ = 0.1 Ṁ_{Edd}, but around different M_{BH} which are: 10 M_{⊙} (blue, solid), 10^{3} M_{⊙} (magenta, dashed), and 10^{6} M_{⊙} (brown, dotted). The more massive the black hole, the brighter the disc since absolute accretion rate increases. In addition, the spectrum spans a larger range of ν, with significant emission from radio to γ rays. It may be noted that higher M_{BH} results in a more broadband spectra, therefore, the disc becomes more luminous but the spectral index does not change much. The spectral properties are presented in Table 5.
4.9. Luminosity, efficiency, and spectral index of twotemperature flows
In Fig. 14a, we calculate the luminosities and in Fig. 14b, we plot the efficiency of the accretion of matter onto BHs of different masses (10 M_{⊙}, 10^{3} M_{⊙} and 10^{6} M_{⊙}) as a function of the accretion rate (Ṁ). The sizes of the circles are in order of increasing value of the BH mass. The parameters used are E = 1.001 and λ = 2.4. We note that luminosity rises steeply with the increase in accretion rate of the system for all BH masses. The more matter is supplied, the greater the conversion of it into energy. However, at higher accretions rates, luminosities approach asymptotic values. Radiation emitted by the accretion disc is the effect of conversion of gravitational energy released in the act of accretion into electromagnetic radiation. So as Ṁ increases, emission increases due to an increased supply of matter. However, it cannot emit more than the energy obtained from the accretion process and, therefore, it reaches a ceiling around η ∼ 10%. It is apparent from Fig. 14b, that the efficiency is slightly affected by the BH mass. The spectral index (α) is represented as colour bar over both Figs. 14a,b. It changes visibly with the increase in accretion rate of the system (also, see Fig. 13a and Table 4) but does not change much with the change in BH mass (also, see Fig. 13b and Table 5).
Fig. 14. (a) Variation of bolometric luminosity (in ergs s^{−1}) and (b) efficiency (in %) as a function of Ṁ (in units of Eddington rate, Ṁ_{edd}). Colour bar indicates the spectral index (α). The BHs of different masses: 10 M_{⊙} (small circle), 10^{3} M_{⊙} (medium circle), and 10^{6} M_{⊙} (largest circle) are represented with increasing sizes of the circles. The parameters used here are E = 1.001 and λ = 2.4. 
5. Discussion and conclusions
In this paper, we study solutions for twotemperature accretion discs around nonrotating BHs. We note that the spin of the BH may play an important role in jet generation via a process called BlandfordZnajek mechanism (Blandford & Znajek 1977), however, accretion is still the primary mechanism for explaining the observed luminosities. A proper twotemperature accretion solution is the best way to obtain the spectra from such systems. Twotemperature equations produce a degenerate set of solutions even when they have the same set of disc parameters such as the generalised Bernoulli parameter (E), accretion rate (Ṁ) and angular momentum (λ). For the given set of disc parameters, a choice of proton temperature may produce a transonic solution through the outer sonic point, while some other choice of the temperature would produce a solution through inner sonic point and yet another would produce solutions which undergo shock transition. The resulting radiation also varies accordingly. In fact some solutions may exhibit four times more luminosity compared to some other solutions (see Fig. 2 and Table 1). Therefore, this degeneracy issue is serious and calls for urgent attention.
In this work, we lay down the methodology to obtain a unique twotemperature solution using the principles of the second law of thermodynamics. We stated that the solution with the highest entropy near the horizon is the correct solution. Since the proposed correct transonic solution is the one with the highest entropy, therefore it is warranted that these solutions should be stable for the relevant boundary conditions. In fact, the collective wisdom of the community on accretion solutions has led to an expectation that close to the horizon, accretion should be transonic. However, for a given set of accretion disc parameters such as Ṁ, E, and λ, we do have a large number of transonic twotemperature solutions and the question of stability of the solutions arises. In Sect. 3.3, we show that the gradient of entropy of the flow with the proton temperature, that is, , is such that it tends to push the solution towards the temperature corresponding to the highest entropy solution. In other words, if our proposed solution is perturbed, then would automatically try to restore the solution to the one corresponding to the highest entropy.
Once we have established the method for obtaining the unique twotemperature solution in a rotating disc, we investigate the effects of various disc parameters on twotemperature accretion solutions. We obtained all possible solutions depending on E and λ for a given Ṁ and M_{BH} (Fig. 7) and, in addition, we also plot the emitted spectrum (Fig. 8). This also shows Ṁ or M_{BH} do not alone determine the emitted spectrum or even the luminosity. Depending on E and λ, the solution changes and, thus, so do the spectrum and luminosity. The constants of motion are uniquely linked to the obtained spectrum.
There are indeed shocked accretion solutions even in the twotemperature regime. The shocked solutions are more luminous because in the postshock region, inverseComptonisation becomes effective, the intensity of the powerlaw photons increases, as compared to a shockfree solution (Fig. 9). We also show that accretion flow can harbour steady shocks in a small but significant patch of the energy angular momentum parameter space (see Fig. 10). However, in general, the shock strength in the twotemperature flow is less than that in a onetemperature flow. Moreover, we did not find any particular spectral signature for the presence of shock, only that the shocked solution is more luminous than shockfree ones. However, it has been found in cases where the accretion rate of the flow is low, and where a weak bremsstrahlung feature is visible (in the high frequency end of the spectrum) in a shockfree solution, that this feature disappears in a shocked solution (see Figs. 8c1,c2).
The radiative properties of a BH system depend on β and β_{d}, along with E, λ, Ṁ, and M_{BH}. For low values of β and β_{d}, the radiative efficiency is around a few percent, but for higher values, the efficiency can easily surpass ten percent, even for the same accretion rate and mass of the BH. We also showed that the spectra becomes broadband if the mass of the central BH considered is higher, it also becomes more luminous but the spectral index remains roughly the same. Whereas, with the increase in accretion rate of the BH, the bandwidth of spectra remains the same, while the luminosity and the spectral index change significantly.
We did not consider pair production from the radiation of the accretion disc, nor did we consider particle production due to high energy interaction of the protons. We show, based on a posteriori calculations, that pair production is negligible and, therefore, the contribution of the pair annihilation emission in the total spectrum of the disc is also negligible. Similarly, we showed, with the help of a posteriori estimate, that the gammaray production from pion decay is also negligible.
In this work, we do not consider the viscosity of the flow, but invoking our results from our previous works on viscous accretion solutions in the single temperature regime (and also Appendix A), we argue that since the angular momentum is almost constant in the inner part of the disc and the viscous heating is much weaker than the magnetic dissipation, we neglect viscosity to ease our computation. Considering the singletemperature viscous flow as the representative case, we may conclude that there would not be any qualitative change resulting from a consideration of viscosity, although a possible quantitative effect cannot be ruled out. In fact, since we show magnetic dissipation to be a very efficient heating process, we could incorporate most heating processes by tuning the parameter β_{d}.
In conclusion, it is absolutely necessary to obtain a unique transonic twotemperature solution. To interpret the relevant observations, the hydrodynamics of the system must be properly handled. Selecting an arbitrary solution would mislead us. We allowed the second law of thermodynamics to dictate and select the solution, without taking recourse to any assumption, so that consistency is maintained. In addition to spherical flow (Sarkar & Chattopadhyay 2019), the accretion disc solution also corresponds to the highest entropy solution.
In a singletemperature regime, the absence of Coulomb coupling makes it easier to integrate the corresponding equation and obtain an analytical measure of entropy for all r (Kumar et al. 2013).
Acknowledgments
The authors acknowledge the anonymous referee for helpful suggestions which improved the quality of the paper. SS acknowledges Mr. Kuldeep Singh for help in python plotting.
References
 Abramowicz, M. A., & Fragile, P. C. 2013, LRR, 16, 1 [Google Scholar]
 Abramowicz, M. A., Czerny, B., Lasota, J. P., & Szuszkiewicz, E. 1988, ApJ, 332, 646 [NASA ADS] [CrossRef] [Google Scholar]
 Abramowicz, M. A., Chen, X., Kato, S., Lasota, J. P., & Regev, O. 1995, ApJ, 438, L37 [NASA ADS] [CrossRef] [Google Scholar]
 Artemova, I. V., BisnovatyiKogan, G. S., Bjoernsson, G., & Novikov, I. D. 1996, ApJ, 456, 119 [NASA ADS] [CrossRef] [Google Scholar]
 Becker, P. A., Das, S., & Le, T. 2008, ApJ, 677, L93 [CrossRef] [Google Scholar]
 Begelman, M. C., & Chiueh, T. 1988, ApJ, 332, 872 [NASA ADS] [CrossRef] [Google Scholar]
 BisnovatyiKogan, G. S., & Lovelace, R. V. E. 1997, ApJ, 486, L43 [NASA ADS] [CrossRef] [Google Scholar]
 BisnovatyiKogan, G. S., & Lovelace, R. V. E. 2001, New Astron. Rev., 45, 663 [NASA ADS] [CrossRef] [Google Scholar]
 Blaes, O. 2014, SSRv, 183, 21 [Google Scholar]
 Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433 [NASA ADS] [CrossRef] [Google Scholar]
 Bondi, H. 1952, MNRAS, 112, 195 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Boldt, E., & Serlemitsos, P. 1969, ApJ, 157, 557 [NASA ADS] [CrossRef] [Google Scholar]
 Chakrabarti, S. K. 1989, ApJ, 347, 365 [NASA ADS] [CrossRef] [Google Scholar]
 Chakrabarti, S. K. 1996, ApJ, 464, 664 [NASA ADS] [CrossRef] [Google Scholar]
 Chakrabarti, S. K., & Titarchuk, L. G. 1995, ApJ, 455, 623 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1939, An Introduction to the Study of Stellar Structure (Chicago: Univ. of Chicago Press) [Google Scholar]
 Chattopadhyay, I. 2008, AIP Conf. Ser., 1053, 353 [NASA ADS] [CrossRef] [Google Scholar]
 Chattopadhyay, I., & Chakrabarti, S. K. 2011, Int. Journ. Mod. Phys. D, 20, 1597 [NASA ADS] [CrossRef] [Google Scholar]
 Chattopadhyay, I., & Kumar, R. 2016, MNRAS, 459, 3792 (CK16) [NASA ADS] [CrossRef] [Google Scholar]
 Chattopadhyay, I., & Ryu, D. 2009, ApJ, 694, 492 [NASA ADS] [CrossRef] [Google Scholar]
 Colpi, M., Maraschi, L., & Treves, A. 1984, ApJ, 280, 319 [CrossRef] [Google Scholar]
 Colpi, M., Maraschi, L., & Treves, A. 1986, ApJ, 311, 150 [CrossRef] [Google Scholar]
 Dahlbacka, G. H., Chapline, G. F., & Weaver, T. A. 1974, Nature, 250, 37 [NASA ADS] [CrossRef] [Google Scholar]
 Dihingia, I. K., Das, S., & Mandal, S. 2017, MNRAS, 475, 2164 [CrossRef] [Google Scholar]
 Eilek, J. A. 1980, ApJ, 236, 664 [CrossRef] [Google Scholar]
 Esin, A. A. 1997, ApJ, 482, 400 [NASA ADS] [CrossRef] [Google Scholar]
 Esin, A. A. 1999, ApJ, 517, 381 [NASA ADS] [CrossRef] [Google Scholar]
 Ferrari, A., Trussoni, E., Rosner, R., & Tsinganos, K. 1985, ApJ, 294, 397 [NASA ADS] [CrossRef] [Google Scholar]
 Fukue, J. 1987, PASJ, 39, 309 [NASA ADS] [Google Scholar]
 Gould, R. J., & Schréder, G. P. 1967, PhRv, 155, 1404 [Google Scholar]
 Holzer, T. E. 1977, JGr, 82, 23 [NASA ADS] [CrossRef] [Google Scholar]
 Hoyle, F., & Lyttleton, R. A. 1939, Proc. Cam. Phil. Soc., 35, 405 [NASA ADS] [CrossRef] [Google Scholar]
 Ichimaru, S. 1977, ApJ, 214, 840 [NASA ADS] [CrossRef] [Google Scholar]
 Ipser, J. P., & Price, R. H. 1982, ApJ, 255, 654 [NASA ADS] [CrossRef] [Google Scholar]
 King, A. 2012, MmSAI, 83, 466 [Google Scholar]
 Kolykhalov, P. I., & Syunyaev, R. A. 1979, SvA, 23, 189 [Google Scholar]
 Kumar, R., & Chattopadhyay, I. 2013, MNRAS, 430, 386 [NASA ADS] [CrossRef] [Google Scholar]
 Kumar, R., & Chattopadhyay, I. 2014, MNRAS, 443, 3444 [NASA ADS] [CrossRef] [Google Scholar]
 Kumar, R., & Chattopadhyay, I. 2017, MNRAS, 469, 4221 [NASA ADS] [CrossRef] [Google Scholar]
 Kumar, R., Singh, C. B., Chattopadhyay, I., & Chakrabarti, S. K. 2013, MNRAS, 436, 2864 [NASA ADS] [CrossRef] [Google Scholar]
 Lanzafame, G., Molteni, D., & Chakrabarti, S. K. 1998, MNRAS, 299, 799 [NASA ADS] [CrossRef] [Google Scholar]
 Lasota, J. P. 1994, Theory of Accretion Discs, eds. W. J. Dushl, J. Frank, F. MeyerHofmeister, & W. M. Tscharnuter (Dordrecht: Kluwer), 2, 341 [CrossRef] [Google Scholar]
 Le, T., & Becker, P. A. 2005, ApJ, 632, 476 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, S.J., Ryu, D., & Chattopadhyay, I. 2011, ApJ, 728, 142 [CrossRef] [Google Scholar]
 Lee, S.J., Chattopadhyay, I., Kumar, R., Hyung, S., & Ryu, D. 2016, ApJ, 831, 33 [NASA ADS] [CrossRef] [Google Scholar]
 Liang, E. P. T., & Thompson, K. A. 1980, ApJ, 240, L271 [CrossRef] [Google Scholar]
 Lightman, A. P., & Eardley, D. M. 1974, ApJ, 187, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Lightman, A. P., Press, W. H., Price, R. H., & Teukolksy, S. 1975, Problem Book in Relativity and Gravitation (New Jersey: Princeton University Press) [Google Scholar]
 Mandal, S., & Chakrabarti, S. K. 2005, A&A, 434, 839 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Manmoto, T., Mineshige, S., & Kusunose, M. 1997, ApJ, 489, 791 [NASA ADS] [CrossRef] [Google Scholar]
 Melo, C. A. H. 2014, Appl. Math. Comput., 232, 1025 [Google Scholar]
 Nakamura, K. E., Kusunose, M., Matsumoto, R., & Kato, S. 1996, PASJ, 48, 761 [CrossRef] [Google Scholar]
 Narayan, R., & Yi, I. 1994, ApJ, 428, L13 [NASA ADS] [CrossRef] [Google Scholar]
 Narayan, R., & Yi, I. 1995, ApJ, 452, 710 [NASA ADS] [CrossRef] [Google Scholar]
 Novikov, I. D., & Thorne, K. S. 1973, Black Holes, eds. B. S. Dewitt, & C. Dewitt (New York: Gordon and Breach), 343 [Google Scholar]
 Pringle, J. E., & Rees, M. J. 1972, A&A, 21, 1 [NASA ADS] [Google Scholar]
 Pringle, J. E., Rees, M. J., & Pacholczyk, A. G. 1973, A&A, 29, 179 [NASA ADS] [Google Scholar]
 Peitz, J., & Appl, S. 1997, MNRAS, 286, 681 (PA97) [NASA ADS] [CrossRef] [Google Scholar]
 Phinney, E. S. 1981, ESASP, 161, 337 [Google Scholar]
 Piran, T. 1978, ApJ, 221, 652 [NASA ADS] [CrossRef] [Google Scholar]
 Rees, M. J., Begelman, M. C., Blandford, R. D., & Phinney, E. S. 1982, Nature, 295, 17 [NASA ADS] [CrossRef] [Google Scholar]
 Rajesh, S. R., & Mukhopadhyay, B. 2010, MNRAS, 402, 961 [CrossRef] [Google Scholar]
 Salpeter, E. E. 1964, ApJ, 140, 796 [CrossRef] [Google Scholar]
 Sarkar, S., & Chattopadhyay, I. 2019, Int. J. Mod. Phys. D, 28, 1950037 (SC19) [CrossRef] [Google Scholar]
 Schwartzman, V. F. 1971, Sov. Astr.AJ, 15, 377 [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
 Shapiro, S. L. 1973, ApJ, 180, 531 [NASA ADS] [CrossRef] [Google Scholar]
 Shapiro, S. L., Lightman, A. P., & Eardley, D. M. 1976, ApJ, 204, 187 (SLE76) [NASA ADS] [CrossRef] [Google Scholar]
 Sharma, P., Quataert, E., Hammett, G. W., & Stone, J. M. 2007, ApJ, 667, 714 [NASA ADS] [CrossRef] [Google Scholar]
 Stepney, S., & Guilbert, P. W. 1983, MNRAS, 204, 1269 [NASA ADS] [Google Scholar]
 Taub, A. H. 1948, Phys. Rev., 74, 328 [NASA ADS] [CrossRef] [Google Scholar]
 Svensson, R. 1982a, ApJ, 258, 321 [NASA ADS] [CrossRef] [Google Scholar]
 Svensson, R. 1982b, ApJ, 258, 335 [NASA ADS] [CrossRef] [Google Scholar]
 Svensson, R. 1984, MNRAS, 209, 175 [NASA ADS] [CrossRef] [Google Scholar]
 Thorne, K. S., & Price, R. H. 1975, ApJ, 195, L101 [NASA ADS] [CrossRef] [Google Scholar]
 Turolla, R., Nobili, L., & Calvani, M. 1986, ApJ, 303, 573 [CrossRef] [Google Scholar]
 Vyas, M. K., Kumar, R., Mandal, S., & Chattopadhyay, I. 2015, MNRAS, 453, 2992 (VKMC15) [NASA ADS] [CrossRef] [Google Scholar]
 Wardziński, G., & Zdziarski, A. 2000, MNRAS, 314, 183 [NASA ADS] [CrossRef] [Google Scholar]
 Weaver, T. A. 1976, Phys. Rev. A, 13, 1563 [NASA ADS] [CrossRef] [Google Scholar]
 Zel’dovich, Y. B. 1964, Sov. Phys. Dok., 9, 195 [Google Scholar]
 Zel’dovich, Y. B., & Novikov, I. D. 1971, Relativistic Astrophysics Stars and Relativity (Chicago: University of Chicago Press), 1 [Google Scholar]
 Zdziarski, A. A. 1985, ApJ, 289, 514 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Effects of viscosity in the system
Here, we discuss the effect of viscosity as (1) a mechanism to remove angular momentum outwards and (2) a source of heating in the system. The handling of viscosity is not trivial and applying it to twotemperature flows might further complicate the scenario and divert us from the question at hand, which is to find the unique transonic twotemperature accretion solutions for rotating flows. Presently, we recall the physics of viscosity in a relativistic but a singletemperature disc like Chattopadhyay & Kumar (2016) and show that near the horizon viscosity would have marginal effect. And in the outer region it would have a more significant effect, but since that region contributes less in the spectrum, so for our purpose, we can neglect it without compromising on the qualitative aspect of the present work. In viscous onetemperature flows we have azimuthal component of radialmomentum equation, the integrated form of which is:
where, L_{0} is the bulk angular momentum at the horizon. L is the local bulk angular momentum expressed as L = hu_{ϕ} = hl. The specific angular momentum is defined as λ = −u_{ϕ}/u_{t} = −l/u_{t}. The dynamical viscosity coefficient is η_{vis} = ρhν_{vis}. Here, ν_{vis} = α_{vis}ar(1 − v^{2})^{2}, is the kinematic viscosity, α_{vis} being the Shakura & Sunyaev viscosity parameter. is the r − ϕ component of the shear tensor which has been evaluated using the expression given in Peitz & Appl (1997), Chattopadhyay & Kumar (2016). We abbreviate this form of viscosity as PA. Using the same procedure followed in Chattopadhyay & Kumar (2016), we find solutions for E = 1.0005, α_{vis} = 0.01, λ_{0} = 2.60, Ṁ = 0.01 Ṁ_{Edd} and M_{BH} = 10 M_{⊙}. Also, we investigate another case of onetemperature flows using the same set of parameters but assuming a form of viscosity which is generally followed in nonrelativistic accretion disc equations and is given by t_{rϕ} = −2η_{vis}σ_{rϕ} = −αP. We abbreviate it as the SS form of viscosity.
We plot in Fig. A.1a the angular momentum distribution of the system as a function of radius where each curve represent PA form of viscosity (green, solid) and SS viscosity (magenta, dotted). It is evident from the figure that angular momentum has been transported outwards (for both the cases) due to the presence of viscosity. Within a distance r ∼ 1000r_{g}, angular momentum is almost constant, similar to the case of inviscid flows. The SS form of viscosity is weaker so it is less efficient in removing angular momentum and hence angular momentum remains almost constant ≲3 × 10^{4}r_{g}. So, neglecting viscosity within these regions does not affect the system qualitatively. Viscosity in addition to removing angular momentum, also heats up the system. We plot in Fig. A.1b the heat that is dissipated due to the presence of the PA form of viscosity (green, solid) and SS form of viscosity (magenta, dotted). We see that using SS viscosity is inefficient with regard to heating up the system and it is always 3 orders of magnitude less than PA viscosity. Very far away from the BH, it is 9 orders of magnitude less than the former. We also compare the heating due to magnetic dissipation (Q_{B}, blue, dashed), which is the source of heating in this work of twotemperature inviscid flows. In Fig. A.1b, Q_{B} has been calculated using the velocity, temperature, and pressure of the corresponding onetemperature flow and using β_{d} = 0.02. We see that Q_{B} is an efficient source of heating in the system. It always supersedes the heating due to the presence of GR viscosity, except in a very small region near the BH, where it becomes comparable. Thus, our assumption of taking Q_{B} as a source of heating in the absence of viscosity suffices for our problem. If viscosity was present, then the total heating would not have changed much, except in a very narrow region.
Fig. A.1. (a) Distribution of specific angular momentum (λ) and (b) heating in the system as a function of radius (log r), when viscosity assumed is relativistic (green, solid) and Newtonian (magenta, dotted). Panel b: we also compare the heating due to magnetic dissipation, Q_{B} (blue, dashed), assuming β_{d} = 0.02. The flow parameters are E = 1.0005, α_{vis} = 0.01, λ_{0} = 2.60, Ṁ = 0.01 Ṁ_{Edd}, and M_{BH} = 10 M_{⊙}. 
Appendix B: Estimation of electronpositron pair production in an advective twotemperature accretion disc
In this section, we estimate how many pairs can be produced using the twotemperature accretion solution as the background solution. There are three processes which could lead to pair (electron and positron) production in accretion discs, namely: photonparticle (electron, positron, or proton) interaction, particleparticle interaction, and photonphoton interaction. Svensson (1982a,b, 1984) gave a detailed analysis of the effect of electronpositron pairs present in relativistic and mildly relativistic plasmas. From these papers, it is apparent that photonphoton interaction is the dominant process responsible for the generation of pairs in accretion discs around BH. The reason behind particleparticle interactions and photonparticle interactions to be of less importance is their excessively small reaction crosssections, which are of the order of 1/137 (value of fine structure constant) and (1/137)^{2}, respectively. The photonphoton pair production process has a threshold condition, which is E_{1}E_{2}(1 − cos θ_{pp}) ≥ 2(m_{e}c^{2})^{2}, where E_{1} and E_{2} are the energies of the photons and θ_{pp} is the angle between these two photons. The advective twotemperature accretion disc solution take into account synchrotron emission which can produce ample amount of soft photons.
These photons are too soft to satisfy the criterion for pair production (Esin 1999). However, these photons, after interacting with highenergy electrons, can get upscattered to high energies, contributing to pair production. Also, the bremsstrahlung emission process produce ample amount of hard photons. Thus, at any particular radius, we can assume the radiation field to be made of a flat bremsstrahlung spectrum which is flat with a high energy cutoff and a Comptonisation spectrum, which is the sum of cutoff power law and the Wien tail (Gould & Schréder 1967; Zdziarski 1985; Esin 1999). We used the formula given by (see Eq. (B1) of the paper Svensson 1984) to compute the rate of photon photon pair production from two power law (PL) photon distribution with an exponentially cutoff, Wien photon interaction (W–W) and power law photons with Wien photons (PL–W). The pair density is estimated a posteriori, by using the temperature profile and velocity profile of the twotemperature solution of this paper. The positron number density is computed from
Here, n_{e+} is the positron number density, S^{±} are the source and the sink terms or pair production and annihilation rate, respectively. S^{±} rates are adopted from Svensson (1984, 1982a, respectively). Since the accretion disc studied in this paper is composed of e^{−} − p^{+} fluid, we first integrate Eq. (B.1) with only the S^{+} term to compute the maximum possible n_{e+} produced. With this distribution of n_{e+}, we compute annihilation rate. We iterate few times till the solution converges. To estimate the production of electronpositron pairs (i.e. e^{−} − e^{+}), we chose two sets of accretion disc parameters presented in the manuscript since the spectrum for this disc parameters is hard and it is possible to obtain significant hard photons which satisfy pairproduction threshold conditon mentioned above. We compute pair production for the case of β_{d} = 0.013 of Fig. 12, and present in Figs. B.1a1, a2 the same two cases Ṁ = 1.0, β = 0.2 and (b1, b2) Ṁ = 1.5, β = 0.15. In the upper panels (a1, b1), we compare the proton number density n_{p} (black, solid), estimated positron number density n_{e+} when annihilation rate is ignored (blue, dotted) and the ones for which both production and annihilation rates are considered (green, dashed). The blue curve is the maximum possible number of positrons that can be produced in the disc and n_{e+} ≪ n_{p}. In the bottom panels (a2, b2), we plot the corresponding emissivity Q_{ann} obtained due to the annihilation of pairs (red, dotted) and compared that with the total emissivity Q_{tot} (brown, solid). The estimated number density of positrons is negligible and the contribution to the total emissivity is negligibly small compared to the total emissivity obtained from radiative processes like, synchrotron, bremsstrahlung, and Comptonisation. This a posteriori estimation justifies our assumption of not considering pair production in the present work.
Fig. B.1. a1, b1: comparison of number density of protons n_{p} (black, solid), positron number densities n_{e+} (without annihilation, blue, dotted), and with both production and annihilation rates (green, dashed). a2, b2: comparison of emissivities of the total radiative cooling Q_{tot} and emissivities due to annihilation of pairs Q_{ann}. For two sets of accretion disc parameters, (a1, a2) Ṁ = 1.0, β = 0.2 and (b) Ṁ = 1.5, β = 0.15. The other parameters are β_{d} = 0.013, E = 1.001, λ = 2.61 and M_{BH} = 10 M_{⊙}. 
Appendix C: Estimation of the gammaray emission by pion interaction
In this section, we discuss whether pion production leads to a viable amount of cooling in the disc and whether it exhibits any observational signature. The reactions leading to pion (π^{±}, π^{0}) production by protonproton interactions are as given below (Eilek 1980):
The threshold temperature for these reactions is 290 MeV. For negative pions temperatures of > 2 GeV are required. Thus, it is can be assumed that negligible π^{−} will be present in the disc. The π^{0} further decay into gammaray photons (Kolykhalov & Syunyaev 1979):
and π^{+} decays into muon and muon neutrinos, which further decays into electron neutrinos, muon antineutrinos and positrons, respectively.
We restrict our study to neutral pions π^{0} since we are interested to study the gamma ray emissivities obtained from an accreting Schwarzschild BH. The rate of π^{0} production is given by (in units of cm^{−3} s^{−1}):
Here, (units of cm^{3} s^{−1}) is the velocityweighted crosssection, which was evaluated for π^{0} by Dahlbacka et al. (1974), assuming experimental cross sections for pion production and a relativistic MaxwellBoltzmann distribution for protons. This was further investigated by Weaver (1976) and Kolykhalov & Syunyaev (1979), who computed for π^{0} as well as for π^{+}. This function is strongly dependent on the proton temperature. In Colpi et al. (1986) obtained a bestfit to these curves, the expression of which is given in Eqs. (10) of their paper and the same form is used to compute the emissivity of γrays. It is clear from Eq. (C.2) that each π^{0} decays into two photons, therefore, the number of photons produced per unit time per unit volume is . We analysed a posteriori the total gammaray luminosity as measured by an observer at infinity using the methodology adopted in Colpi et al. (1984).
We checked the production of γ rays by varying the accretion rate of the system from Ṁ = 0.01 (red, dotted) to 0.10 (green, dashed) and 1.0 (blue, solid) for a 10 M_{⊙} BH (see Figs. C.1a1, a2). Since this emission is crucially dependent on the proton temperature, we plot log T_{p} as a function of log r in Fig. C.1a1 for the different accretion rates. The set of disc parameters used are λ = 2.61 and E = 1.0007. For Ṁ = 0.01, & 0.1, the accretion solution does not undergo shock transition (see, Fig. 9) and the T_{p} distribution of the global accretion solution is similar. However, for the same set of E, λ but Ṁ = 1.0, there is a stable accretion shock and the T_{p} jumps at the shock location. The corresponding spectra (Fig. C.1a2) for the three values of Ṁ show a marked difference where the shocked accretion solution is more luminous and the spectrum is harder (solid, blue), and becomes less luminous and softer for lower Ṁ. The gammaray emission (calculated a posteriori) is represented in grey colour. As a result of π^{0} decay, the contribution in the high energy regime increases. Quantitatively, the gammaray luminosity for different accretion rates are related by the following relation L_{γph} (Ṁ = 1.0) ≃ 200 L_{γph}(Ṁ = 0.1) and L_{γph} (Ṁ = 0.1) ≃ 100 L_{γph} (Ṁ = 0.01). Although L_{γph} (Ṁ = 1.0) is much higher than that compared to lower Ṁs but compared to the overall luminosity for each Ṁ, L_{γph} is pitiably low. We chose Ṁ = 0.1 and the same values of E and λ and then studied the gammaray production for accretion discs onto different masses of central BH. We plot the proton temperature distribution T_{p} vs r in loglog scale and the corresponding spectra in Figs. C.1b1, b2, for central BH masses M_{BH} = 10^{2} (blue, solid), 10^{4} (green, dashed) and 10^{6} (red, dotted). The M_{BH} affects the system quantitatively but not qualitatively. While the luminosity increases and spectra become more broad band with the increase in M_{BH}, but the efficiency of gammaray emission (L_{γph}/(Ṁc^{2})) remains almost same ∼10^{−8}. In both cases (change in accretion rate and mass of BH), the fractional change in luminosity is always < 10^{−5}. This analysis justifies our assumption of neglecting pion production leading to emission of γrays.
Fig. C.1. a1, a2: dependence on the accretion rate, Ṁ = 0.01 (red, dotted), 0.1 (green, dashed), and 1.0 (blue, solid). (a1) log T_{p} as a function of log r and (a2) log νL_{ν} with log ν. (b1, b2) Dependence on mass of BH, M_{BH} = 10^{2} (blue, solid), 10^{4} (green, dashed), 10^{6} (red, dotted) in the production of γ rays. (b1) log T_{p} as a function of log r and in (b2) the corresponding spectra are plotted. The gamma ray emission is presented in grey. 
All Tables
All Figures
Fig. 1. Method for finding sonic points. Solutions are presented in terms of Mach number M (=v/a) vs log r plot. Θ_{pin} = 7.162 × 10^{−2} for all iterations. Panel a: iterations to obtain inner sonic point r_{ci} (black circle) and panel b: iterations to obtain outer sonic point r_{co} (black star). Various branches plotted are multivalued branches (MB; green dasheddot), transonic (TS; red dashed), and supersonic (SB; blue dasheddot). Respective Θ_{ein}s are mentioned inside the panels. Panel c: full set of transonic solutions. Global accretion solution (red solid) through r_{co} and accretion solution through r_{ci} (red, dashed) is plotted. Equatorial global wind (through r_{co}) and nonglobal wind (through r_{ci}) are represented using red dotted curve. The accretion disc flow parameters used are λ = 2.5, E = 1.000045, Ṁ = 0.001 Ṁ_{Edd}, and M_{BH} = 10 M_{⊙}. 

In the text 
Fig. 2. Left: M vs log r plot for various values of T_{pin}. (a) T_{pin} = 3.1 × 10^{11} K, (b) T_{pin} = 5.605 × 10^{11}K, (c) T_{pin} = 6.04 × 10^{11} K, (d) T_{pin} = 6.460 × 10^{11} K, (e) T_{pin} = 6.554 × 10^{11} K and (f) T_{pin} = 7.0 × 10^{11} K. Global solutions are represented by solid lines. Panel g: vs T_{pin}. Solid black curve is for the solutions passing through outer sonic point, while dotted black curve is for solutions passing through inner sonic point. Panels a–f: solutions corresponding to the points marked in right panel g. The disc flow parameters are E = 1.0015, λ = 2.6 and Ṁ = 0.02 Ṁ_{Edd}. The spacetime is described by a BH of mass 10 M_{⊙}. 

In the text 
Fig. 3. Stability analysis of the unique transonic twotemperature solution with maximum entropy. The flow parameters used are same as Fig. 2. is plotted against variation of T_{pin}. The arrows indicate that Δ converge at T_{pin} = T_{pinmax} (blue dot) and is the stable equilibrium solution. This T_{pin} is the solution with maximum entropy marked “c” in Fig. 2g. 

In the text 
Fig. 4. (a) plotted against T_{pin}. The entropy for inner sonic point solutions (black, dotted) and outer sonic points (black, solid) are presented. T_{pin} marked with green triangle corresponds to maximum entropy solution. Flow variables plotted are (b) M (green, solid), (c) v (cyan, solid) and v_{ϕ} (blue, dotted), (d) log n (magenta, solid), (e) T_{p} (red, dotted), T_{e} (orange, solid), (f) Γ_{p} (red, dotted), Γ_{e} (orange, solid), (g) E (brown, solid) as functions of log r. The flow parameters are E = 1.000045, λ = 2.5, Ṁ = 0.001 Ṁ_{Edd} and M_{BH} = 10 M_{⊙}. Panel b: the sonic point is marked with a black star. 

In the text 
Fig. 5. Emissivity vs log r plot for the flow presented in Fig. 4, shown here in the top three panels. Bottom panel: spectrum of the accretion flow. 

In the text 
Fig. 6. (a) M and (b) log T vs log r, and (c) total spectrum (black solid) plus contribution from various length scales of the accretion disc, 2 − 3r_{g} (magenta, dotted), 3 − 5r_{g} (blue, dashed), 5 − 8r_{g} (green, shortdasheddotted), 8 − 10r_{g} (brown, longdasheddotted), 10 − 100r_{g} (orange, dasheddoubledotted) and 100 − 1000r_{g} (red, dashedtripledotted). Flow parameters are E = 1.0002, λ = 2.48, Ṁ = 0.05 Ṁ_{Edd}, and M_{BH} = 10 M_{⊙}. 

In the text 
Fig. 7. Variation of solutions, M as a function of log r with variation of E and λ. From left to right: the specific energy increases as E = 1.0005, 1.001, 1.003, and 1.01. From top to bottom: the angular momentum increases as λ = 2.40, 2.55, 2.70, and 2.85. The other parameters are Ṁ = 0.01 Ṁ_{Edd} and M_{BH} = 10 M_{⊙}. 

In the text 
Fig. 8. Variation of the spectrum with E and λ. The set of values for E and λ and other parameters are same as in Fig. 7. 

In the text 
Fig. 9. Typical shocked solution (a), with its corresponding number density (b), emissivities, (c) and spectrum (d), is presented. The parameters taken are E = 1.002, λ = 2.58, Ṁ = 0.2 Ṁ_{Edd}, and M_{BH} = 10 M_{⊙}. 

In the text 
Fig. 10. Shock parameter space for Ṁ = 0.01 Ṁ_{Edd} (green, solid), 0.10 Ṁ_{Edd} (blue, dashed), and 1.00 Ṁ_{Edd} (red, dotted) around a 10 M_{⊙} BH. 

In the text 
Fig. 11. Change in spectra with increase in β = 0.002 (blue, solid), 0.01 (green, dashed), and 0.02 (red, dotted). Other parameters used are E = 1.003, λ = 2.54, and Ṁ = 0.1 Ṁ_{Edd} in an accretion disc around M_{BH} = 10 M_{⊙}. 

In the text 
Fig. 12. Accretion solutions (a1, b1) and their corresponding spectra (a2, b2) plotted for a flow with E = 1.001, λ = 2.61 around M_{BH} = 10 M_{⊙}. Various curves are for β_{d} = 0.013 (blue, solid), β_{d} = 0.015 (green, dashed) and β_{d} = 0.017 (red, dotted). The accretion rates and ratio of magnetic to gas pressure are chosen as Ṁ = 1.0 Ṁ_{Edd}, β = 0.2 (a1, a2), and Ṁ = 1.5 Ṁ_{Edd}, β = 0.15 (b1, b2). 

In the text 
Fig. 13. Spectra from (a) M_{BH} = 10 M_{⊙} for different accretion rates Ṁ = 0.1 Ṁ_{Edd} (blue, solid), Ṁ = 0.6 Ṁ_{Edd} (green, dashed), and Ṁ = 1.2 Ṁ_{Edd} (red, dotted); (b) Ṁ = 0.1 Ṁ_{Edd} but around M_{BH} = 10 M_{⊙} (blue, solid), M_{BH} = 10^{3} M_{⊙} (magenta, dashed), and M_{BH} = 10^{6} M_{⊙} (brown, dotted). Other disc parameters are E = 1.001 and λ = 2.4. 

In the text 
Fig. 14. (a) Variation of bolometric luminosity (in ergs s^{−1}) and (b) efficiency (in %) as a function of Ṁ (in units of Eddington rate, Ṁ_{edd}). Colour bar indicates the spectral index (α). The BHs of different masses: 10 M_{⊙} (small circle), 10^{3} M_{⊙} (medium circle), and 10^{6} M_{⊙} (largest circle) are represented with increasing sizes of the circles. The parameters used here are E = 1.001 and λ = 2.4. 

In the text 
Fig. A.1. (a) Distribution of specific angular momentum (λ) and (b) heating in the system as a function of radius (log r), when viscosity assumed is relativistic (green, solid) and Newtonian (magenta, dotted). Panel b: we also compare the heating due to magnetic dissipation, Q_{B} (blue, dashed), assuming β_{d} = 0.02. The flow parameters are E = 1.0005, α_{vis} = 0.01, λ_{0} = 2.60, Ṁ = 0.01 Ṁ_{Edd}, and M_{BH} = 10 M_{⊙}. 

In the text 
Fig. B.1. a1, b1: comparison of number density of protons n_{p} (black, solid), positron number densities n_{e+} (without annihilation, blue, dotted), and with both production and annihilation rates (green, dashed). a2, b2: comparison of emissivities of the total radiative cooling Q_{tot} and emissivities due to annihilation of pairs Q_{ann}. For two sets of accretion disc parameters, (a1, a2) Ṁ = 1.0, β = 0.2 and (b) Ṁ = 1.5, β = 0.15. The other parameters are β_{d} = 0.013, E = 1.001, λ = 2.61 and M_{BH} = 10 M_{⊙}. 

In the text 
Fig. C.1. a1, a2: dependence on the accretion rate, Ṁ = 0.01 (red, dotted), 0.1 (green, dashed), and 1.0 (blue, solid). (a1) log T_{p} as a function of log r and (a2) log νL_{ν} with log ν. (b1, b2) Dependence on mass of BH, M_{BH} = 10^{2} (blue, solid), 10^{4} (green, dashed), 10^{6} (red, dotted) in the production of γ rays. (b1) log T_{p} as a function of log r and in (b2) the corresponding spectra are plotted. The gamma ray emission is presented in grey. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.