© 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 super-massive 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 X-ray 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 non-thermal 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 gas-pressure dominated region. Shapiro et al. (1976, hereafter SLE76) considered this geometrically thick and optically thin, puffed-up 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 X-1 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 two-temperature 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 advection-dominated accretion flows (ADAFs) and Ichimaru (1977) was the first to propose it (for review also see Bisnovatyi-Kogan & 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 sub-Keplerian 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 one-temperature accretion flows. One-temperature 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 two-temperature 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 two-temperature 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 free-fall 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 two-temperature, optically thin accretion discs. The equations of motion were solved under the self-similar assumption. It is to be noted that self-similarity in accretion solution around BH is plausible only at a large distance from the horizon and not near it. Additionally, a self-similar 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 two-temperature 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 two-components in the accretion flow: a Keplerian and another sub-Keplerian. 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 non-thermal 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 free-fall, but the shocked solution and its signature in the observed spectrum were qualitatively studied.

Nakamura et al. (1996) obtained the first global transonic, two-temperature 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 inverse-Comptonisation). Recently, a more general transonic advective two-temperature 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 two-temperature solutions, including accretion-shock, were investigated by Dihingia et al. (2017). The most striking feature of the transonic two-temperature studies mentioned above is that the solutions depend on the choice of inner or outer boundary conditions, however, based on single-temperature 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 one-temperature and two-temperature regimes are exactly the same, but there is one more variable in case of two-temperature flows (i.e. the existence of different ion and electron temperature distributions). Obtaining a self-consistent two-temperature flow would require solving basic hydrodynamic equations, but since there is one more flow variable in the two-temperature 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 two-temperature problem would require all the equations to be solved self-consistently 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 two-temperature 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 two-temperature 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 one-temperature regime, the transonic solution is unique (Le & Becker 2005; Becker et al. 2008) but for reasons cited above, a two-temperature 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 two-temperature energy equation within the adiabatic limit is not integrable (unlike in one-temperature 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 two-temperature 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 two-temperature 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 two-temperature 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 trans-relativistic in nature. This means that very far away from the horizon, matter is thermally non-relativistic and as it approaches the BH, it could be sub-relativistic or relativistic. Matter is referred as thermally relativistic when its thermal energy is comparable to or greater than its rest mass energy, kT/mc² ≳ 1 (where k= Boltzmann constant, T= temperature, m= mass of the species), and its adiabatic index, (Γ) ∼4/3; it is non-relativistic if its thermal energy is less than the rest mass energy (kT/mc² <  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 two-temperature 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 trans-relativistic 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 trans-relativistic 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 two-temperature solution. We present the results in Sect. 4 and our conclusions in Sect. 5.

2. Two-temperature 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, two-temperature flows are degenerate, which makes the method for obtaining a unique solution very important. In the next sub-sections, we discuss the equations we used to model the two-temperature accretion flow. We also give an overview of the radiative processes we incorporated and the methodology we used to implement these processes in curved space-time.

2.1. Equations of motion (EoM)

The space-time metric around a non-rotating BH is described using a Schwarzschild metric:

$$\begin{array}{c}\hfill \mathrm{d}{s}^2=g_{\mathit{tt}}{c}^2\mathrm{d}{t}^2+g_{\mathit{rr}}\mathrm{d}{r}^2+g_{\theta \theta }\mathrm{d}{\theta }^2+g_{\varphi \varphi }\mathrm{d}{\varphi }^2,\end{array}$$(1)

where the metric tensors are expressed as, $-g_{\mathit{tt}}=g_{\mathit{rr}}^{-1}=(1-2G{M}_{\mathrm{BH}}/c^{2}r)$ and g_θθ = g_ϕϕ = r², 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²,  c and r_g/c = GM_BH/c³, 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 axis-symmetric, 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:

$$\begin{array}{c}\hfill u^r\frac{\mathrm{d}{u}^r}{\mathrm{d}r}+\frac{1}{r^2}-(r-3)u^{\varphi }{u}^{\varphi }+(g^{\mathit{rr}}+u^r{u}^r)\frac{1}{e+p}\frac{\mathrm{d}p}{\mathrm{d}r}=0,\end{array}$$(2)

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 four-velocity. The mass accretion rate is obtained by integrating the conservation of four mass-flux:

$$\begin{array}{c}\hfill \dot{M}=4\pi \rho H{u}^{r}r,\end{array}$$(3)

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 half-height of the disc. Here, Ṁ – the accretion rate, is a constant of motion throughout the flow. Half-height is calculated assuming hydrostatic equilibrium along the vertical direction of the disc (Lasota 1994; Chattopadhyay & Chakrabarti 2011) which can be written as

$$\begin{array}{c}\hfill H=\sqrt{\frac{p{r}^3}{\rho {\gamma }_{\varphi }^2}}·\end{array}$$(4)

Also, based on u_μu^μ = −1, we obtain:

$$\begin{array}{c}\hfill -u_t=\sqrt{(1-\frac{2}{r})}{\gamma }_v{\gamma }_{\varphi },\end{array}$$

where, γ_v and γ_ϕ are the Lorentz factors in the radial and azimuthal directions, respectively and defined as ${\gamma }_v=\sqrt{1/(1-v^2)}$ and ${\gamma }_{\varphi }=\sqrt{1/(1-v_{\varphi }^2)}$, where $v_{\varphi }=\sqrt{-u_{\varphi }{u}^{\varphi }/u_t{u}^t}$ and v is the velocity in the local co-rotating frame. It can be shown that $v^2={\gamma }_{\varphi }^2{v}_{\widehat{r}}^2$, where $v_{\widehat{r}}=\sqrt{-u_r{u}^r/u_t{u}^t}$. The total Lorentz factor, therefore is, γ = γ_vγ_ϕ.

The first law of thermodynamics, or the energy balance equation, is $u_{\mu }{T}_{;\nu }^{\mu \nu }=\mathrm{\Delta }Q$ and can be written as:

$$\begin{array}{c}\hfill u^r[\left(\frac{e+p}{\rho }\right){\rho }_{,r}-e_{,r}]=\mathrm{\Delta }Q.\end{array}$$(5)

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⁻³ s⁻¹, 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 single-temperature case, with Coulomb coupling being infinitely strong, protons and electrons equilibrate locally into a single-temperature distribution (Lanzafame et al. 1998; Lee et al. 2011, 2016; Chattopadhyay & Kumar 2016). However for two-temperature 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

$$\begin{array}{c}\hfill E=-h{u}_t\exp (X_f),\end{array}$$(6)

where, h = (e + p)/ρ= specific enthalpy and $X_f=\int \frac{\mathrm{\Delta }{Q}_{\text{p}}+\mathrm{\Delta }{Q}_{\text{e}}}{\rho h{u}^r}\mathrm{d}r$. Here, $\mathrm{\Delta }{Q}_i=Q_i^{+}-Q_i^{-}$ 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 sub-section, barred variables have been used to denote dimensional quantities and unbarred variables are, as before, non-dimensional. In order to solve the above equations of motion, we need to supply an equation of state (EoS). In this paper, we use the Chattopadhyay-Ryu (CR) EoS given by Chattopadhyay (2008), Chattopadhyay & Ryu (2009). The form of CR EoS for multispecies flow is:

$$\begin{array}{c}\hfill \overline{e}={\displaystyle \sum _i{\overline{e}}_i=\sum _i[{\overline{n}}_i{m}_i{c}^2+{\overline{p}}_i\left(\frac{9{\overline{p}}_i+3{\overline{n}}_i{m}_i{c}^2}{3{\overline{p}}_i+2{\overline{n}}_i{m}_i{c}^2}\right)],}\end{array}$$(7)

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:

$$\begin{array}{c}\hfill \overline{n}={\displaystyle \sum _i{\overline{n}}_i={\overline{n}}_p+{\overline{n}}_e=2{\overline{n}}_e,}\end{array}$$(8)

$$\begin{array}{c}\hfill \overline{\rho }={\displaystyle \sum _i{\overline{n}}_i{m}_i={\overline{n}}_{\mathrm{e}}{m}_{\mathrm{e}}+{\overline{n}}_{\mathrm{p}}{m}_{\mathrm{p}}={\overline{n}}_{\mathrm{e}}{m}_{\mathrm{e}}(1+\frac{1}{\eta })={\overline{n}}_{\mathrm{e}}{m}_{\mathrm{e}}\stackrel{\sim }{K},}\end{array}$$(9)

$$\begin{array}{c}\hfill \overline{p}={\displaystyle \sum _i{\overline{p}}_i=\sum _i{\overline{n}}_{i}k{T}_i={\overline{n}}_{\mathrm{e}}k(T_{\mathrm{e}}+T_{\mathrm{p}})={\overline{n}}_{\mathrm{e}}{m}_{\mathrm{e}}{c}^2({\mathrm{\Theta }}_{\mathrm{e}}+\frac{{\mathrm{\Theta }}_{\mathrm{p}}}{\eta }),}\end{array}$$(10)

where, η = m_e/m_p and $\stackrel{\sim }{K}=1+1/\eta$. T_i is the temperature in units of kelvin, while ${\mathrm{\Theta }}_i=\frac{k{T}_i}{m_i{c}^2}$ is the non – dimensional temperature

defined with respect to the rest-mass energy of the respective ith species. The EoS, Eq. (7), can be simplified using Eqs. (8)–(10) to

$$\begin{array}{c}\hfill \overline{e}={\overline{n}}_{\mathrm{e}}{m}_{\mathrm{e}}{c}^2(f_{\mathrm{e}}+\frac{f_{\mathrm{p}}}{\eta })=\frac{\overline{\rho }{c}^{2}f}{\stackrel{\sim }{K}},\end{array}$$(11)

where, $f_i=1+{\mathrm{\Theta }}_i\left(\frac{9{\mathrm{\Theta }}_i+3}{3{\mathrm{\Theta }}_i+2}\right)$ and $f=f_{\text{e}}+\frac{f_{\text{p}}}{\eta }\cdot$

Polytropic index and adiabatic index can be written as

$$\begin{array}{c}\hfill N_i=\frac{\mathrm{d}{f}_i}{\mathrm{d}{\mathrm{\Theta }}_i}\mathrm{and}{\mathrm{\Gamma }}_i=1+\frac{1}{N_i}·\end{array}$$(12)

The equation for half-height (Eq. (4)), can be simplified to:

$$\begin{array}{c}\hfill H=\sqrt{\frac{[r^3-{\lambda }^2(r-2)]}{\stackrel{\sim }{K}}({\mathrm{\Theta }}_{\mathrm{e}}+\frac{{\mathrm{\Theta }}_{\mathrm{p}}}{\eta })}.\end{array}$$(13)

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:

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\mathrm{\Theta }}_{\mathrm{p}}}{\mathrm{d}r}=-\frac{2{\mathrm{\Theta }}_{\mathrm{p}}}{2N_{\mathrm{p}}+1}(\mathcal{A}+\frac{1}{v(1-v^2)}\frac{\mathrm{d}v}{\mathrm{d}r})-\mathbb{P}\eta ,\end{array}$$(14)

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\mathrm{\Theta }}_{\mathrm{e}}}{\mathrm{d}r}=-\frac{2{\mathrm{\Theta }}_{\mathrm{e}}}{2N_{\mathrm{e}}+1}(\mathcal{A}+\frac{1}{v(1-v^2)}\frac{\mathrm{d}v}{\mathrm{d}r})-\mathbb{E},\end{array}$$(15)

respectively, where,

$\mathcal{A}=-\frac{r}{r(r-2)}-\frac{3r^2-{\lambda }^2}{2[r^3-{\lambda }^2(r-2)]},\mathbb{P}=\frac{2\mathrm{\Delta }{Q}_{\text{p}}\tilde{K}}{\rho u^r(2N_{\text{p}}+1)},$, and $\mathbb{E}=\frac{2\mathrm{\Delta }{Q}_{\text{e}}\stackrel{\sim }{K}}{\rho u^r(2N_{\text{e}}+1)}.$

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 three-velocity, which has the form:

$$\begin{array}{c}\hfill \frac{\mathrm{d}v}{\mathrm{d}r}=\frac{\mathcal{N}}{\mathcal{D}},\end{array}$$(16)

where, $\mathcal{N}=-\frac{1}{r(r-2)}+\frac{{\lambda }^2{\gamma }_{\varphi }^2(r-3)}{r^4}+a^2\mathcal{A}+\frac{{\mathrm{\Gamma }}_{\text{p}}{N}_{\text{p}}\mathbb{P}+{\mathrm{\Gamma }}_{\text{e}}{N}_{\text{e}}\mathbb{E}}{h\stackrel{\sim }{K}}-\frac{\Delta Q}{\rho h{u}^r}$, and $\mathcal{D}=\frac{v}{1-v^2}\left(1-\frac{a^2}{v^2}\right).$

The effective sound speed (a) has been defined as

$$\begin{array}{c}\hfill a^2=\frac{\mathcal{G}}{h\stackrel{\sim }{K}},\end{array}$$(17)

where, $\mathcal{G} = \frac{2{\mathrm{\Gamma }}_{\text{p}}{N}_{\text{p}}{\mathrm{\Theta }}_{\text{p}}}{\eta (2N_{\text{p}}+1)}+\frac{2{\mathrm{\Gamma }}_{\text{e}}{N}_{\text{e}}{\mathrm{\Theta }}_{\text{e}}}{(2N_{\text{e}}+1)}\cdot$

2.3. Heating and cooling processes included in the flow

In the following sub-sections, 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² and magnetic energy density (B²/8π) by 1/r⁴. 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

$$\begin{array}{c}\hfill {\overline{Q}}_{\mathrm{B}}=\frac{3u^{r}c}{2r{r}_{\mathrm{g}}}\frac{B^2}{8\pi }=\frac{3u^{r}c}{2r{r}_{\mathrm{g}}}{\beta }_{\mathrm{d}}\overline{p}=\frac{3u^{r}c}{2r{r}_{\mathrm{g}}}{\beta }_{\mathrm{d}}{\overline{n}}_{\mathrm{e}}k(T_{\mathrm{p}}+T_{\mathrm{e}})\mathrm{ergs}{\mathrm{cm}}^{-3}{\mathrm{s}}^{-1}.\end{array}$$(18)

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

$$\begin{array}{c}\hfill Q_{\mathrm{p}}^{+}=\delta Q_{\mathrm{B}}\mathrm{and}{Q}_{\mathrm{e}}^{+}=(1-\delta )Q_{\mathrm{B}},\end{array}$$(19)

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 mass-loss 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 two-temperature 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 two-temperature 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 sub-Keplerian. 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 one-temperature, 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ϕ = −α_visp (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 self-consistent two-temperature, transonic, rotating accretion solutions by considering low-angular-momentum 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,

$$\begin{array}{c}\hfill Q_{\mathrm{p}}^{-}=Q_{\mathrm{e}}^{+}=Q_{\mathrm{ep}}.\end{array}$$(20)

The expression for Coulomb coupling in cgs units (ergs cm⁻³s⁻¹) is given by Stepney & Guilbert (1983),

$$\begin{array}{cc}\hfill {\overline{Q}}_{\mathrm{ep}}& =\frac{3}{2}\frac{m_{\mathrm{e}}}{m_{\mathrm{p}}}{\overline{n}}_{\mathrm{e}}{\overline{n}}_{\mathrm{p}}{\sigma }_{T}ck\frac{T_{\mathrm{p}}-T_{\mathrm{e}}}{K_2(1/{\mathrm{\Theta }}_{\mathrm{e}})K_2\left(1/{\mathrm{\Theta }}_{\mathrm{p}}\right)}\text{ln}{\mathrm{\Lambda }}_c\hfill \\ \hfill & \times [\frac{2{({\mathrm{\Theta }}_{\mathrm{e}}+{\mathrm{\Theta }}_{\mathrm{p}})}^2+1}{{\mathrm{\Theta }}_{\mathrm{e}}+{\mathrm{\Theta }}_{\mathrm{p}}}{K}_1\left(\frac{{\mathrm{\Theta }}_{\mathrm{e}}+{\mathrm{\Theta }}_{\mathrm{p}}}{{\mathrm{\Theta }}_{\mathrm{e}}{\mathrm{\Theta }}_{\mathrm{p}}}\right)+2K_0\left(\frac{{\mathrm{\Theta }}_{\mathrm{e}}+{\mathrm{\Theta }}_{\mathrm{p}}}{{\mathrm{\Theta }}_{\mathrm{e}}{\mathrm{\Theta }}_{\mathrm{p}}}\right)],\hfill \end{array}$$(21)

where σ_T is the Thomson scattering cross-section, 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 magneto-rotational 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 (${\overline{Q}}_{\mathrm{ib}}$) is a radiative loss term for the protons, the expression of which is given as (Boldt & Serlemitsos 1969),

$$\begin{array}{c}\hfill {\overline{Q}}_{\mathrm{ib}}=1.4\times {10}^{-27}{\overline{n}}_{\mathrm{e}}^2\sqrt{\frac{m_{\mathrm{e}}}{m_{\mathrm{p}}}{T}_{\mathrm{p}}}.\mathrm{ergs}{\mathrm{cm}}^{-3}{\mathrm{s}}^{-1}.\end{array}$$(22)

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) inverse-Comptonisation (Q_ic). Emissivity due to bremsstrahlung (in ergs cm⁻³ s⁻¹) is given by Novikov et al. (1973),

$$\begin{array}{c}\hfill {\overline{Q}}_{\mathrm{br}}=1.4\times {10}^{-27}{\overline{n}}_{\mathrm{e}}^2\sqrt{T_{\mathrm{e}}}(1+4.4\times {10}^{-10}{T}_{\mathrm{e}}).\end{array}$$(23)

We used thermal synchrotron radiation in our model, following the prescription of Wardziński & Zdziarski (2000). The emissivity is given by:

$$\begin{array}{c}\hfill {\overline{Q}}_{\mathrm{syn}}=\frac{2\pi }{3}\frac{{\nu }_t^3}{r{r}_g}{m}_{\mathrm{e}}{\mathrm{\Theta }}_{\mathrm{e}},\end{array}$$(24)

where, ν_t is the turnover frequency, above which the plasma is optically thin to synchrotron radiation and below which it is highly self-absorbed 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, $B=\sqrt{8\pi \beta \overline{p}}$. 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 inverse-Comptonised by the electrons present in the plasma. This is given by (Wardziński & Zdziarski 2000),

$$\begin{array}{c}\hfill {\overline{Q}}_{\mathrm{ic}}=\zeta {\overline{Q}}_{\mathrm{syn}},\end{array}$$(25)

where, ζ is the enhancement factor. It is expressed as

$$\begin{array}{c}\hfill \zeta =3\phi {\left(\frac{x_t}{{\mathrm{\Theta }}_{\mathrm{e}}}\right)}^{{\alpha }_0-1}[{\mathrm{\Gamma }}_{\mathrm{inc}}(1-{\alpha }_0,\frac{x_t}{{\mathrm{\Theta }}_{\mathrm{e}}})+\frac{6{\mathrm{\Gamma }}_{\mathrm{inc}}({\alpha }_0)P_{\mathrm{sc}}}{{\mathrm{\Gamma }}_{\mathrm{inc}}(2{\alpha }_0+3)}].\end{array}$$

Here, Γ_inc is the incomplete gamma function, $x_t=\frac{h{\nu }_t}{m_e{c}^2},$ $\phi =\frac{[1+{(2{\mathrm{\Theta }}_{\text{e}})}^2]}{[1+10{(2{\mathrm{\Theta }}_{\text{e}})}^2]}, {\alpha }_0$ is the spectral index which can be defined as:

$$\begin{array}{c}\hfill {\alpha }_0=-\frac{\ln P_{\mathrm{sc}}}{\ln A}·\end{array}$$(26)

It is the slope of the power-law photons that are generated due to inverse-Comptonisation at each radius. Therefore, the net spectral index (α) of the final inverse-Compton spectrum is obtained from the contributions of all the values of α₀ from each radius of the disc. In Eq. (26), $A=1+4{\mathrm{\Theta }}_{\text{e}}+16{\mathrm{\Theta }}_{\text{e}}{}^2$, 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),

$$\begin{array}{c}\hfill {\tau }_{\mathrm{es}}=0.4{[1+{(2.22T_{\mathrm{e}}\times {10}^{-9})}^{0.86}]}^{-1}\overline{\rho }H{r}_{\mathrm{g}}.\end{array}$$(27)

2.3.5. Compton heating

As discussed above, less energetic photons would cool the flow through the process of inverse-Comptonisation. 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 inverse-Comptonisation, but the sign changes (Esin 1997), causing heating rather than cooling. Therefore,

$$\begin{array}{c}\hfill {\overline{Q}}_{\mathrm{e}}^{+}={\overline{Q}}_{\mathrm{comp}}.\end{array}$$(28)

2.3.6. Final expressions for ΔQ_p and ΔQ_e

So, to conclude, we have taken for heating and cooling of protons:

$$\begin{array}{c}\hfill Q_{\mathrm{p}}^{+}=\delta Q_{\mathrm{B}},\mathrm{and}{Q}_{\mathrm{p}}^{-}=Q_{\mathrm{ep}}+Q_{\mathrm{ib}},\end{array}$$(29)

respectively. And for heating and cooling of electrons:

$$\begin{array}{c}\hfill Q_{\mathrm{e}}^{+}=(1-\delta )Q_{\mathrm{B}}+Q_{\mathrm{ep}}+Q_{\mathrm{comp}},\mathrm{and}{Q}_{\mathrm{e}}^{-}=Q_{\mathrm{br}}+Q_{\mathrm{syn}}+Q_{\mathrm{ic}}\end{array}$$(30)

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)):

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\mathrm{\Theta }}_{\mathrm{p}}}{\mathrm{d}r}=\frac{{\mathrm{\Theta }}_{\mathrm{p}}}{N_{\mathrm{p}}}\frac{1}{n_{\mathrm{p}}}\frac{\mathrm{d}{n}_{\mathrm{p}}}{\mathrm{d}r}+\frac{Q_{\mathrm{ep}}\eta \stackrel{\sim }{K}}{\rho u^r{N}_{\mathrm{p}}}\mathrm{and}\end{array}$$(31)

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\mathrm{\Theta }}_{\mathrm{e}}}{\mathrm{d}r}=\frac{{\mathrm{\Theta }}_{\mathrm{e}}}{N_{\mathrm{e}}}\frac{1}{n_{\mathrm{e}}}\frac{\mathrm{d}{n}_{\mathrm{e}}}{\mathrm{d}r}-\frac{Q_{\mathrm{ep}}\stackrel{\sim }{K}}{\rho u^r{N}_{\mathrm{e}}}·\end{array}$$(32)

Due to the presence of the Coulomb interaction term, we cannot integrate the above equation and obtain an analytical form¹. 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:

$$\begin{array}{cc}& n_{\mathrm{e}\mathrm{in}}={\kappa }_1\exp \left(\frac{f_{\mathrm{e}\mathrm{in}}}{{\mathrm{\Theta }}_{\mathrm{e}\mathrm{in}}}\right){\mathrm{\Theta }}_{\mathrm{e}\mathrm{in}}^{\frac{3}{2}}{(3{\mathrm{\Theta }}_{\mathrm{e}\mathrm{in}}+2)}^{\frac{3}{2}}\hfill \end{array}$$(33)

$$\begin{array}{cc}\hfill & n_{\mathrm{p}\mathrm{in}}={\kappa }_2\exp \left(\frac{f_{\mathrm{p}\mathrm{in}}}{{\mathrm{\Theta }}_{\mathrm{p}\mathrm{in}}}\right){\mathrm{\Theta }}_{\mathrm{p}\mathrm{in}}^{\frac{3}{2}}{(3{\mathrm{\Theta }}_{\mathrm{p}\mathrm{in}}+2)}^{\frac{3}{2}},\hfill \end{array}$$(34)

where κ₁ and κ₂ 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,

$$\begin{array}{c}\hfill n_{\mathrm{in}}^2=n_{\mathrm{e}\mathrm{in}}{n}_{\mathrm{p}\mathrm{in}}\Rightarrow n_{\mathrm{in}}=\sqrt{n_{\mathrm{e}\mathrm{in}}{n}_{\mathrm{p}\mathrm{in}}}.\end{array}$$(35)

Thus, we can write

$$\begin{array}{c}\hfill n_{\mathrm{in}}=\kappa \sqrt{\exp \left(\frac{f_{\mathrm{e}\mathrm{in}}}{{\mathrm{\Theta }}_{\mathrm{e}\mathrm{in}}}\right)\exp \left(\frac{f_{\mathrm{p}\mathrm{in}}}{{\mathrm{\Theta }}_{\mathrm{p}\mathrm{in}}}\right){\mathrm{\Theta }}_{\mathrm{e}\mathrm{in}}^{\frac{3}{2}}{\mathrm{\Theta }}_{\mathrm{p}\mathrm{in}}^{\frac{3}{2}}{(3{\mathrm{\Theta }}_{\mathrm{e}\mathrm{in}}+2)}^{\frac{3}{2}}{(3{\mathrm{\Theta }}_{\mathrm{p}\mathrm{in}}+2)}^{\frac{3}{2}}},\end{array}$$(36)

where, $\kappa =\sqrt{{\kappa }_1{\kappa }_2.}$

The expression of the entropy accretion rate, using Eq. (3), can be written as

$$\begin{array}{cc}\hfill {\dot{\mathcal{M}}}_{\mathrm{in}}& =\frac{\dot{M}}{4\pi \kappa (m_{\mathrm{e}}+m_{\mathrm{p}})}\hfill \\ \hfill & =\frac{4\pi n_{\mathrm{in}}(m_{\mathrm{e}}+m_{\mathrm{p}})H_{\mathrm{in}}{u}_{\mathrm{in}}^r{r}_{\mathrm{in}}}{4\pi \kappa (m_{\mathrm{e}}+m_{\mathrm{p}})}\hfill \\ \hfill & =\left[\sqrt{\exp \left(\frac{f_{\mathrm{e}\mathrm{in}}}{{\mathrm{\Theta }}_{\mathrm{e}\mathrm{in}}}\right)\exp \left(\frac{f_{\mathrm{p}\mathrm{in}}}{{\mathrm{\Theta }}_{\mathrm{p}\mathrm{in}}}\right){\mathrm{\Theta }}_{\mathrm{e}\mathrm{in}}^{\frac{3}{2}}{\mathrm{\Theta }}_{\mathrm{p}\mathrm{in}}^{\frac{3}{2}}{(3{\mathrm{\Theta }}_{\mathrm{e}\mathrm{in}}+2)}^{\frac{3}{2}}{(3{\mathrm{\Theta }}_{\mathrm{p}\mathrm{in}}+2)}^{\frac{3}{2}}}\right]\hfill \\ \hfill & \times H_{\mathrm{in}}{u}_{\mathrm{in}}^r{r}_{\mathrm{in}}.\hfill \end{array}$$(37)

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:

$$\begin{array}{cc}& -\frac{1}{r_{\mathrm{c}}(r_{\mathrm{c}}-2)}+\frac{{\lambda }_{\mathrm{c}}^2{\gamma }_{{\varphi }_{\mathrm{c}}}^2(r_{\mathrm{c}}-3)}{r_{\mathrm{c}}^4}+a_{\mathrm{c}}^2{\mathcal{A}}_{\mathrm{c}}+\frac{{\mathrm{\Gamma }}_{\mathrm{pc}}{N}_{\mathrm{pc}}{\mathbb{P}}_{\mathrm{c}}+{\mathrm{\Gamma }}_{\mathrm{ec}}{N}_{\mathrm{ec}}{\mathbb{E}}_{\mathrm{c}}}{h_{\mathrm{c}}\stackrel{\sim }{K}}\hfill \\ \hfill & -\frac{\mathrm{\Delta }{Q}_{\mathrm{c}}}{{\rho }_{\mathrm{c}}{h}_{\mathrm{c}}{u}_{\mathrm{c}}^r}=0\hfill \end{array}$$(38)

and

$$\begin{array}{c}\hfill \frac{v_{\mathrm{c}}}{1-v_{\mathrm{c}}^2}(1-\frac{a_{\mathrm{c}}^2}{v_{\mathrm{c}}^2})=0\Rightarrow v_{\mathrm{c}}=a_{\mathrm{c}}.\end{array}$$(39)

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 Rankine-Hugoniot conditions (Taub 1948) are:

Conservation of mass flux across the shock: [Ṁ] = 0

Conservation of energy flux: [Ė] = 0

Conservation of momentum flux: [$[\Sigma h{\gamma }_v^2{v}^2+W]=0$,

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 inverse-Comptonisation, which give rise to emissions spanning the whole electromagnetic spectrum. This emission (measured in units of ergs s⁻¹ Hz⁻¹) 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 special-relativistic transformations this becomes:

$$\begin{array}{c}\hfill j{\prime }_{\nu \prime }=j_{\nu }\frac{1-v^2}{{(1-v\cos \theta \prime )}^2}\mathrm{and}\nu \prime =\nu \frac{\sqrt{1-v^2}}{(1-vcos\theta \prime )}·\end{array}$$(40)

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):

$$\begin{array}{c}\hfill |\cos {\theta }^{\ast }|=\sqrt{\frac{27}{4}{\left(\frac{2}{r}\right)}^2\left(\frac{2}{r}\right)+1},\end{array}$$(41)

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 $\sqrt{1-2/r}$ into the observed frequency. Furthermore, if we want to calculate the bolometric luminosity of the system, we need to integrate the frequency-dependent 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 high-frequency 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). Inverse-Comptonisation, on the other hand, is identified as a power-law 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 X-type or saddle-type. 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 O-type (imaginary slope) and spiral-type (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 X-type, while the middle sonic may be spiral or O-type depending on whether the system is dissipative or non-dissipative, 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 two-temperature 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 self-consistently 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:

1. 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).

2. As r_in → 2, E → ℰ = −h_inu_t = h_in(1 − 2/r_in)^1/2γ_vγ_ϕ. Here h_in = h_in(Θ_pin, Θ_ein) is the specific enthalpy at r_in.

3. We supply Θ_pin.

4. We also supply Θ_ein. Then v_in is obtained as a function of the flow parameters and can be expressed as

   $$v_{\mathrm{in}}={[1-\frac{(1-2/r_{\mathrm{in}})}{{\mathcal{E}}^2}\frac{r_{\mathrm{in}}^3}{\{r_{\mathrm{in}}^3-{\lambda }^2(r_{\mathrm{in}}-2)\}}{h}_{\mathrm{in}}^2]}^{1/2}.$$

5. 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)).

6. 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.

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

8. 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 rci (black circle) and panel b: iterations to obtain outer sonic point rco (black star). Various branches plotted are multivalued branches (MB; green dashed-dot), transonic (TS; red dashed), and supersonic (SB; blue dashed-dot). Respective Θeins are mentioned inside the panels. Panel c: full set of transonic solutions. Global accretion solution (red solid) through rco and accretion solution through rci (red, dashed) is plotted. Equatorial global wind (through rco) and non-global wind (through rci) are represented using red dotted curve. The accretion disc flow parameters used are λ = 2.5, E = 1.000045, Ṁ = 0.001 ṀEdd, and MBH = 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⁻² 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 multi-valued branch (MB) of solutions. We plot one such MB solution (green dashed-dot) 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 dashed-dot) 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 dashed-dot) 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 dashed-dot). 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 two-temperature 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, two-temperature 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 ($\dot{\mathcal{M}}$) 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 ${\dot{\mathcal{M}}}_{\text{in}}$. If the computed ${\dot{\mathcal{M}}}_{\text{in}}$ is plotted with respect to Θ_pin, then there is a clear maxima in ${\dot{\mathcal{M}}}_{\text{in}}$. Following the second law of thermodynamics, the solution (corresponding to a Θ_pin) with maximum entropy (${\dot{\mathcal{M}}}_{\text{in}}$) is the physically plausible solution. Hence, we are able to constrain the degeneracy and obtain a unique transonic two-temperature solution for a given set of constants of motion.

In Figs. 2a–g, we illustrate the methodology for obtaining a unique two-temperature 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 ${\dot{\mathcal{M}}}_{\text{in}},$ 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_pins (points marked “a”–“f”) which are (a) 3.1 × 10¹¹ K, (b) 5.605 × 10¹¹ K, (c) 6.04 × 10¹¹ K, (d) 6.460 × 10¹¹ K, (e) 6.554 × 10¹¹ K and (f) 7.0 × 10¹¹ 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 ${\dot{\mathcal{M}}}_{\text{in}}$ (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 ${\dot{\mathcal{M}}}_{\text{in}}$ 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 ${\dot{\mathcal{M}}}_{\text{in}}$ and yet, the solutions are completely different.

Fig. 2. Left: M vs log r plot for various values of Tpin. (a) Tpin = 3.1 × 1011 K, (b) Tpin = 5.605 × 1011K, (c) Tpin = 6.04 × 1011 K, (d) Tpin = 6.460 × 1011 K, (e) Tpin = 6.554 × 1011 K and (f) Tpin = 7.0 × 1011 K. Global solutions are represented by solid lines. Panel g: ${\dot{\mathcal{M}}}_{\text{in}}$ vs Tpin. 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 space-time is described by a BH of mass 10 M⊙.

In order to drive the point all the way home, we present $T_{\text{pin}},{\dot{\mathcal{M}}}_{\text{in}},r_{\text{c}},v_{\mathrm{c}}$, & 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 two-temperature 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 ${\dot{\mathcal{M}}}_{\text{in}}$). In this particular case, the highest entropy solution is the one corresponding to point “c” (red triangle) in ${\dot{\mathcal{M}}}_{\text{in}}$–T_pin curve (Fig. 2g), and the correct, unique two-temperature 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 two-temperature 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 two-temperature solution.

Let us assign ${\dot{\mathcal{M}}}_{\text{in|max}}=\max ({\dot{\mathcal{M}}}_{\text{in}})$ and T_pin|max = T_pin for which ${\dot{\mathcal{M}}}_{\text{in}}$ is maximum. Let us further define δT_pin which is the difference between adjacent higher and lower T_pin and $\mathrm{\Delta }=(\mathrm{d}{\dot{\mathcal{M}}}_{\text{in}}/\mathrm{d}{T}_{\text{pin}})$ as the gradient of ${\dot{\mathcal{M}}}_{\text{in}}$. Now the change in ${\dot{\mathcal{M}}}_{\text{in}}$ is given by

$$\begin{array}{c}\hfill \delta {\dot{\mathcal{M}}}_{\mathrm{in}}=\left(\frac{\mathrm{d}{\dot{\mathcal{M}}}_{\mathrm{in}}}{\mathrm{d}{T}_{\mathrm{pin}}}\right)\delta T_{\mathrm{pin}}=\mathrm{\Delta }\delta T_{\mathrm{pin}}.\end{array}$$(42)

It is clear that Δ = 0 for any extrema of ${\dot{\mathcal{M}}}_{\text{in}}-T_{\text{pin}}$ curve, but the solution is said to be stable if T_pin moves away from the value T_pin|max 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_pin|max, the figure shows that Δ >  0. Then from Eq. (42), $\delta {\dot{\mathcal{M}}}_{\text{in}}>0$. So the system tends to go to higher ${\dot{\mathcal{M}}}_{\text{in}}$, which we denote using a rightward arrow. Similarly, for T_pin >  T_pin|max, Δ <  0 which implies $\delta {\dot{\mathcal{M}}}_{\text{in}}<0$. Since, based on the second law of thermodynamics, a physical system is not amenable to decreasing its entropy, therefore for T_pin >  T_pin|max, the system would tend to come back to T_pin|max. We represent this in the figure by using a leftward arrow. In other words, entropy can only increase ($\delta {\dot{\mathcal{M}}}_{\text{in}}>0$), if T_pin → T_pin|max from either side of T_pin|max. Hence, a solution corresponding to T_pin = T_pin|max must be stable. Therefore, in addition to the fact that T_pin|max corresponds to a solution with maximum entropy, we can conclude that the solution is also stable.

Fig. 3. Stability analysis of the unique transonic two-temperature solution with maximum entropy. The flow parameters used are same as Fig. 2. $\mathrm{\Delta }=(\mathrm{d}{\dot{\mathcal{M}}}_{\text{in}}/\mathrm{d}{T}_{\text{pin}})$ is plotted against variation of Tpin. The arrows indicate that Δ converge at Tpin = Tpin|max (blue dot) and is the stable equilibrium solution. This Tpin is the solution with maximum entropy marked “c” in Fig. 2g.

4. Results

Two-temperature 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 two-temperature 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 two-temperature solutions

In Fig. 4, we study a typical two-temperature 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 ${\dot{\mathcal{M}}}_{\text{in}}$ vs T_pin. The solution with T_pin = 6.0 × 10¹¹ 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 three-velocity (v; cyan solid) in co-rotating 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⁻³) 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 single-temperature 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 non-relativistic. 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 mildly-relativistic 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) ${\dot{\mathcal{M}}}_{\text{in}}$ plotted against Tpin. The entropy for inner sonic point solutions (black, dotted) and outer sonic points (black, solid) are presented. Tpin 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) Tp (red, dotted), Te (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 MBH = 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 inverse-Comptonisation (blue dashed-dotted line). For the present set of disc parameters, at the outer boundary of the disc, the temperatures are non-relativistic 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 inverse-Comptonisation becomes important and exceeds bremsstrahlung. However, inverse-Comptonisation 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 inverse-Comptonisation by a blue dashed-dotted line. The overall luminosity of the system is low, namely, L = 2.536 × 10²⁹ ergs s⁻¹, with an radiative efficiency of η_r = 1.957 × 10⁻⁵. It may be noted that efficiency is defined as η_r = L/(Ṁc²). 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, short-dashed-dotted), 8 − 10r_g (brown, long-dashed-dotted), 10 − 100r_g (orange, dashed-double-dotted) and 100 − 1000r_g (red, dashed-triple-dotted). The spectra from all these regions are separately over plotted in Fig. 6c, the colour-coding 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³ − 10⁴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 low-energy end of the spectra from this region is always around and above 10¹² 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 − 3rg (magenta, dotted), 3 − 5rg (blue, dashed), 5 − 8rg (green, short-dashed-dotted), 8 − 10rg (brown, long-dashed-dotted), 10 − 100rg (orange, dashed-double-dotted) and 100 − 1000rg (red, dashed-triple-dotted). Flow parameters are E = 1.0002, λ = 2.48, Ṁ = 0.05 ṀEdd, and MBH = 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. Low-angular-momentum 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 MBH = 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 energy-angular momentum parameter space of a two-temperature accretion disc, as is shown in Figs. 7a1–d4.

Figure 8 shows the corresponding spectra, which span from 10¹² − 10²² 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 inverse-Comptonisation) 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 inverse-Comptonisation 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 two-temperature 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 ($R=u_{-}^r/u_{+}^r$, ± implies post and pre-shock 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 inverse-Comptonisation in dashed-dotted 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 inverse-Comptonisation, especially in the post-shock 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 MBH = 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 shock-free solution (blue dashed of panel a). The luminosity of the system is 2.831 × 10³⁶ ergs s⁻¹, which corresponds to an efficiency (η_r) of 1.09%, while for a shock-free branch (blue dashed), the luminosity would be 1.293 × 10³⁶ ergs s⁻¹ 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 power-law 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 shock-parameter 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 inverse-Comptonisation. 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 inverse-Comptonisation, 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, inverse-Comptonisation 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 MBH = 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 super-Eddington 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 MBH = 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 two-temperature solutions and resulting spectra. However, it may be noted that we have ignored pair production from particle-particle interactions in the accretion disc or from accretion disc radiations. We have also ignored the production of gamma-rays 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 (π⁰) 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² : blue, solid, 10⁴ : green, dashed, and 10⁶ : 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 gamma-ray produced by decay of pions. The same trend is observed when we increase the BH mass. However, the gamma ray luminosity is always < 10⁻⁵ 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 inverse-Compton 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) MBH = 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 MBH = 10 M⊙ (blue, solid), MBH = 103 M⊙ (magenta, dashed), and MBH = 106 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³ M_⊙ (magenta, dashed), and 10⁶ 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 two-temperature 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³ M_⊙ and 10⁶ 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 electro-magnetic 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), 103 M⊙ (medium circle), and 106 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 two-temperature accretion discs around non-rotating BHs. We note that the spin of the BH may play an important role in jet generation via a process called Blandford-Znajek mechanism (Blandford & Znajek 1977), however, accretion is still the primary mechanism for explaining the observed luminosities. A proper two-temperature accretion solution is the best way to obtain the spectra from such systems. Two-temperature 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 two-temperature 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 two-temperature 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, $\mathrm{d}{\dot{\mathcal{M}}}_{\text{in}}/\mathrm{d}{T}_{\text{pin}}$, 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 $\mathrm{d}{\dot{\mathcal{M}}}_{\text{in}}/\mathrm{d}{T}_{\text{pin}}$ 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 two-temperature solution in a rotating disc, we investigate the effects of various disc parameters on two-temperature 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 two-temperature regime. The shocked solutions are more luminous because in the post-shock region, inverse-Comptonisation becomes effective, the intensity of the power-law photons increases, as compared to a shock-free 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 two-temperature flow is less than that in a one-temperature flow. Moreover, we did not find any particular spectral signature for the presence of shock, only that the shocked solution is more luminous than shock-free 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 shock-free 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 broad-band 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 gamma-ray 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 single-temperature 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 two-temperature 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.

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

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 two-temperature flows might further complicate the scenario and divert us from the question at hand, which is to find the unique transonic two-temperature accretion solutions for rotating flows. Presently, we recall the physics of viscosity in a relativistic but a single-temperature 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 one-temperature flows we have azimuthal component of radial-momentum equation, the integrated form of which is:

$$\begin{array}{c}\hfill -\rho u^r(L-L_0)=2\eta {\sigma }_{\varphi }^r,\end{array}$$(A.1)

where, L₀ 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 = α_visar(1 − v²)², is the kinematic viscosity, α_vis being the Shakura & Sunyaev viscosity parameter. ${\sigma }_{\varphi }^r$ 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, λ₀ = 2.60, Ṁ = 0.01 Ṁ_Edd and M_BH = 10 M_⊙. Also, we investigate another case of one-temperature flows using the same set of parameters but assuming a form of viscosity which is generally followed in non-relativistic 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⁴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 two-temperature inviscid flows. In Fig. A.1b, Q_B has been calculated using the velocity, temperature, and pressure of the corresponding one-temperature 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, QB (blue, dashed), assuming βd = 0.02. The flow parameters are E = 1.0005, αvis = 0.01, λ0 = 2.60, Ṁ = 0.01 ṀEdd, and MBH = 10 M⊙.

Appendix B: Estimation of electron-positron pair production in an advective two-temperature accretion disc

In this section, we estimate how many pairs can be produced using the two-temperature accretion solution as the background solution. There are three processes which could lead to pair (electron and positron) production in accretion discs, namely: photon-particle (electron, positron, or proton) interaction, particle-particle interaction, and photon-photon interaction. Svensson (1982a,b, 1984) gave a detailed analysis of the effect of electron-positron pairs present in relativistic and mildly relativistic plasmas. From these papers, it is apparent that photon-photon interaction is the dominant process responsible for the generation of pairs in accretion discs around BH. The reason behind particle-particle interactions and photon-particle interactions to be of less importance is their excessively small reaction cross-sections, which are of the order of 1/137 (value of fine structure constant) and (1/137)², respectively. The photon-photon pair production process has a threshold condition, which is E₁E₂(1 − cos θ_pp) ≥ 2(m_ec²)², where E₁ and E₂ are the energies of the photons and θ_pp is the angle between these two photons. The advective two-temperature 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 high-energy 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 cut-off and a Comptonisation spectrum, which is the sum of cut-off 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 cut-off, 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 two-temperature solution of this paper. The positron number density is computed from

$$\begin{array}{c}\hfill n_{e^{+}}=\frac{1}{\mathit{rH}}\int (S^{+}-S^{-})rH\mathrm{d}r.\end{array}$$(B.1)

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 electron-positron 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 pair-production 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 np (black, solid), positron number densities ne+ (without annihilation, blue, dotted), and with both production and annihilation rates (green, dashed). a2, b2: comparison of emissivities of the total radiative cooling Qtot and emissivities due to annihilation of pairs Qann. 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 MBH = 10 M⊙.

Appendix C: Estimation of the gamma-ray 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 (π^±,  π⁰) production by proton-proton interactions are as given below (Eilek 1980):

$$\begin{array}{cc}& p+p\to p+n+{\pi }^{+}\hfill \\ \hfill & p+p\to p+p+{\pi }^0\hfill \\ \hfill & p+p\to d+{\pi }^{+}.\hfill \end{array}$$(C.1)

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 π⁰ further decay into gamma-ray photons (Kolykhalov & Syunyaev 1979):

$$\begin{array}{c}\hfill {\pi }^0\to 2{\gamma }_{\mathrm{ph}}\end{array}$$(C.2)

and π⁺ decays into muon and muon neutrinos, which further decays into electron neutrinos, muon anti-neutrinos and positrons, respectively.

$$\begin{array}{cc}& {\pi }^{+}\to {\mu }^{+}+{\nu }_{\mu }\hfill \\ \hfill & {\mu }^{+}\to {\nu }_e+e^{+}+{\overline{\nu }}_{\mu }\hfill \end{array}$$(C.3)

We restrict our study to neutral pions π⁰ since we are interested to study the gamma ray emissivities obtained from an accreting Schwarzschild BH. The rate of π⁰ production is given by (in units of cm⁻³ s⁻¹):

$$\begin{array}{c}\hfill {\mathcal{R}}_{{\pi }^0}=\frac{n^2}{2}{⟨\overline{\sigma }\overline{v}⟩}_{{\pi }^0}.\end{array}$$(C.4)

Here, ${⟨\overline{\sigma }\overline{v}⟩}_{{\pi }^0}$ (units of cm³ s⁻¹) is the velocity-weighted cross-section, which was evaluated for π⁰ by Dahlbacka et al. (1974), assuming experimental cross sections for pion production and a relativistic Maxwell-Boltzmann distribution for protons. This was further investigated by Weaver (1976) and Kolykhalov & Syunyaev (1979), who computed $⟨\overline{\sigma }\overline{v}⟩$ for π⁰ as well as for π⁺. This function is strongly dependent on the proton temperature. In Colpi et al. (1986) obtained a best-fit 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 π⁰ decays into two photons, therefore, the number of photons produced per unit time per unit volume is $n^2{⟨\overline{\sigma }\overline{v}⟩}_{{\pi }^0}$. We analysed a posteriori the total gamma-ray 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 gamma-ray emission (calculated a posteriori) is represented in grey colour. As a result of π⁰ decay, the contribution in the high energy regime increases. Quantitatively, the gamma-ray 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 gamma-ray production for accretion discs onto different masses of central BH. We plot the proton temperature distribution T_p vs r in log-log scale and the corresponding spectra in Figs. C.1b1, b2, for central BH masses M_BH = 10² (blue, solid), 10⁴ (green, dashed) and 10⁶ (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 gamma-ray emission (L_γph/(Ṁc²)) remains almost same ∼10⁻⁸. In both cases (change in accretion rate and mass of BH), the fractional change in luminosity is always < 10⁻⁵. 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 Tp as a function of log r and (a2) log νLν with log ν. (b1, b2) Dependence on mass of BH, MBH = 102 (blue, solid), 104 (green, dashed), 106 (red, dotted) in the production of γ rays. (b1) log Tp 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 rci (black circle) and panel b: iterations to obtain outer sonic point rco (black star). Various branches plotted are multivalued branches (MB; green dashed-dot), transonic (TS; red dashed), and supersonic (SB; blue dashed-dot). Respective Θeins are mentioned inside the panels. Panel c: full set of transonic solutions. Global accretion solution (red solid) through rco and accretion solution through rci (red, dashed) is plotted. Equatorial global wind (through rco) and non-global wind (through rci) are represented using red dotted curve. The accretion disc flow parameters used are λ = 2.5, E = 1.000045, Ṁ = 0.001 ṀEdd, and MBH = 10 M⊙.

In the text

Fig. 2. Left: M vs log r plot for various values of Tpin. (a) Tpin = 3.1 × 1011 K, (b) Tpin = 5.605 × 1011K, (c) Tpin = 6.04 × 1011 K, (d) Tpin = 6.460 × 1011 K, (e) Tpin = 6.554 × 1011 K and (f) Tpin = 7.0 × 1011 K. Global solutions are represented by solid lines. Panel g: ${\dot{\mathcal{M}}}_{\text{in}}$ vs Tpin. 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 space-time is described by a BH of mass 10 M⊙.

In the text

Fig. 3. Stability analysis of the unique transonic two-temperature solution with maximum entropy. The flow parameters used are same as Fig. 2. $\mathrm{\Delta }=(\mathrm{d}{\dot{\mathcal{M}}}_{\text{in}}/\mathrm{d}{T}_{\text{pin}})$ is plotted against variation of Tpin. The arrows indicate that Δ converge at Tpin = Tpin|max (blue dot) and is the stable equilibrium solution. This Tpin is the solution with maximum entropy marked “c” in Fig. 2g.

In the text

Fig. 4. (a) ${\dot{\mathcal{M}}}_{\text{in}}$ plotted against Tpin. The entropy for inner sonic point solutions (black, dotted) and outer sonic points (black, solid) are presented. Tpin 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) Tp (red, dotted), Te (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 MBH = 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 − 3rg (magenta, dotted), 3 − 5rg (blue, dashed), 5 − 8rg (green, short-dashed-dotted), 8 − 10rg (brown, long-dashed-dotted), 10 − 100rg (orange, dashed-double-dotted) and 100 − 1000rg (red, dashed-triple-dotted). Flow parameters are E = 1.0002, λ = 2.48, Ṁ = 0.05 ṀEdd, and MBH = 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 MBH = 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 MBH = 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 MBH = 10 M⊙.
In the text