Inflow-outflow around black holes

The 2D disk structure with advective transonic inflow-outflow solutions around black holes

Abstract

We solved analytically viscous two-dimensional (2D) fluid equations for accretion and outflows in spherical polar coordinates () and obtained explicitly flow variables in and directions around black holes (BHs). We investigated global transonic advection-dominated accretion flow (ADAF) solutions in direction on an equatorial plane with using Paczyński-Wiita potential. We used radial flow variables of ADAFs with symmetric conditions on the equatorial plane, as initial values for integration in direction. In the study of 2D disk structure, we used two-azimuthal components of viscous stress tensors namely, and . Interestingly, we found that the whole advective disk is not participating in outflow generation and the outflows form close to the BHs. Normally, outflow strength increased with increasing viscosity parameter (), mass-loss parameter () and decreasing gas pressure ratio (). Outflow region increased with increasing , for and decreasing for . The is effective in angular momentum transportation at high latitude and outflows collimation along an axis of symmetry since it changes polar velocity () of the flow. The outflow emission is also affected by the ADAF size and decreased with decreasing it. Transonic surfaces formed for both inflows (, very close to BH) and outflows (). We also explored no outflows, outflows and failed outflows regions, which mainly depend on the viscosity parameters.

accretion, accretion disks – black hole physics – hydrodynamics
1

1 Introduction

An accretion disk is associated with many astrophysical objects, e.g., compact objects (black holes, neutron stars, and white dwarfs) and young stellar objects. Accreting gas onto these objects can generate radiation and bipolar outflows/jets due to an extraction of its gravitational energy. These objects are often associated with non-relativistic to relativistic bipolar jets. Especially, the relativistic jets have been observed around accreting black hole candidates (BHCs), for instance, active galactic nuclei (AGNs) and black hole ray binaries (BHXBs). The jets from the AGN M87 are emerged from an extremely small central region of a source within (Junor et al. , 1999) but recent observation shows that even smaller region less than (Doeleman et al. , 2012), where is a Schwarzschild radius. AGNs and BHXBs are believed to harbor supermassive BHs and stellar mass BHs at the center, respectively. denotes mass of the Sun. Moreover, the BHXBs are also showing typically two type of spectral states in their observations, one high soft state, which is radiatively efficient and dominated by thermal radiation with black body spectrum in soft X-ray regime and second low hard state, which is radiatively inefficient and dominated by non-thermal radiation with some power law spectrum in the hard X-ray regime (Remillard & McClintock, 2006). These two states are also connected with many intermediate states and interestingly, bipolar jets and quasi-periodic oscillations are associated with the hard state in the BHXBs (Gallo et al. , 2003; Fender et al. , 2004). However, such changes in the spectral states are yet to be observed for the AGNs. Since the time scales of AGNs and BHXBs can be scaled by the mass of the BHs, but inner boundary conditions are same, so the basic physics of both kind of objects can be similar (McHardy et al. , 2006). In this context, there are a several theoretical and numerical studies on accretion processes with Keplerian/sub-Keplerian flows and that can play an important role in generations of soft spectrum (Shakura & Sunyaev, 1973; Novikov & Thorne, 1973; Abramowicz et al. , 1988) and hard spectral state (Sunyaev & Titarchuk, 1980; Chakrabarti & Titarchuk, 1995; Narayan & Yi, 1995b; Esin et al. , 1997), hard/soft state transitions (Wu et al. , 2016), and the outflows from the accretion disks around BHs (Narayan & Yi, 1995a; Molteni et al. , 1996a, b; Igumenshchev & Abramowicz, 1999, 2000; Ohsuga et al. , 2005; Okuda et al. , 2007; Yang et al. , 2014; Das et al. , 2014; Bu et al. , 2016a, b; Lee et al. , 2016; Kumar & Chattopadhyay, 2017; Jiang et al. , 2017). The mechanism for the jet generation and evolution of it is still not much clear and a topic of active research in the fields of theory and observations.

The present paper is based on the study of structure of accretion disk in the 2D flow with outflows and extension of previous studies (Narayan & Yi, 1995a; Xu & Chen, 1997; Blandford & Begelman, 1999; Xue & Wang, 2005; Xie & Yuan, 2008; Jiao & Wu, 2011). In this paper, there are following things, which are new from the previous 2D analytical studies as one by one, since we want to study inflow-outflow structure close to the BH, therefore we used pseudo-Newtonian potential (Paczyński & Wiita, 1980). Which incorporates general relativistic effects very close to the BH. Second, we considered two-azimuthal components of viscous stress tensor (Stone et al. , 1999; Yuan et al. , 2012a) from out off nine-components (Xue & Wang, 2005). Since it is mostly believed that the azimuthal component of magnetic stress is more important in the angular momentum transfer by the study of magneto-rotational instability (MRI) simulations (Balbus & Hawley, 1998). The present study is axisymmetric with the 2D HD rotating flow so we used anomalous shear stress, which can approximates the magnetic stress and following Stone et al. (1999); Yuan et al. (2012a) for assuming two-azimuthal components of the viscosity are non-zero. Their effects on the disk structures have also discussed in the simulation by Yang et al. (2014). Third, we have investigated global transonic ADAF solutions on the equatorial plane and their flow variables used as the boundary conditions for the integration of differential equations in direction. When doing this, the fluid ODEs are still dependent on the radius of the disk and other radial derivatives, unlike with using the self-similar assumptions (Xue & Wang, 2005; Jiao & Wu, 2011). In present study, our main interest is to investigate the inflow-outflow structure close to the BHs with changing various flow parameters, namely, disk viscosity parameters (), grand specific energy (), gas pressure ratio () and mass-loss parameter () in the fluid flows. The structure of this paper is in next section 2, model fluid equations and assumptions, section 3, solution procedure, section 4, numerical results and in the last section, conclusions of our study.

2 Model fluid Equations and Assumptions

We considered viscous hydrodynamic fluid equations for advective accretion-outflow solutions with steady-state and axisymmetric in the spherical polar coordinates (). We assumed pseudo-Newtonian geometry (Paczyński & Wiita, 1980) around the Schwarzschild BHs. For time being, we are ignoring magnetic field in the accretion and outflows. We represented the viscous fluid equations and the flow variables in geometrical unit system and chosen , where, and are mass of the BH, universal gravitational constant and speed of the light, respectively. Therefore, units of a length, flow velocity ( or sound speed), energy, specific angular momentum, mass, density, pressure and time are , , , , , , and , respectively. We also assumed that the two-components of viscous stress tensor are effective in plane, which are and as following Stone et al. (1999). Thus, the conserved form of the fluid equations in the 2D become as,
the continuity equation,

 1r2∂∂r(r2ρvr)+1rsinθ∂∂θ(ρvθsinθ)=0, (1)

the components of Navier-Stokes equation,

 vr∂vr∂r+vθr∂vr∂θ−v2θ+v2ϕr+1ρ∂P∂r−Fr=0, (2)
 vr∂vθ∂r+vθr∂vθ∂θ+vrvθr−v2ϕcotθr+1rρ∂P∂θ=0, (3)
 vr∂vϕ∂r+vθr∂vϕ∂θ+vϕr(vr+vθcotθ)=1ρr[1r2∂∂r(r3trϕ)+∂tθϕ∂θ+2tθϕcotθ], (4)

and the energy equation,

 ρ[vr∂ϵ∂r+vθr∂ϵ∂θ−Pρ{vrρ∂ρ∂r+vθrρ∂ρ∂θ}]=fQ+, (5)

where is viscous heating rate and is advection factor (Narayan & Yi, 1995a). For simplicity we assumed is fixed. Since the values of should not be arbitrary, therefore for brevity we used only for highly advective flow, inspite of radiation-dominated or gas-dominated flow. However, should be determined with relevant radiation mechanisms. () is total pressure, is gas pressure and is radiation pressure, which could be due to blackbody emissivity for optically thick medium (Abramowicz et al. , 1988) or bremsstrahlung and synchrotron emissivity for optically thin medium (Narayan & Yi, 1995b). is dimensionless temperature of the fluid, , where , , and are the Boltzmann constant, mean molecular weight of the gas, mass of the proton and mass of the electron, respectively. We assumed for fully ionized flow. is central attractive force around the BH. The specific internal energy (Kato et al. , 2008; Jiao & Wu, 2011) is

where is known as adiabatic index and defined as ratio between heat capacities. is effective , is polytropic index and is the gas pressure ratio. The two-azimuthal components of the viscous stress tensors are written as,

 τrϕ=η1(∂vϕ∂r−vϕr)  and  τθϕ=η2r(∂vϕ∂θ−vϕcotθ), (7)

where and are coefficients of viscosity and is Keplerian angular velocity on the equatorial plane. The and are the Shakura-Sunyaev viscosity parameters. The flow variables in the plane are defined (Xue & Wang, 2005) as:

 Mass density  ρ(r,θ)=ρ=ρ1(θ)ρ2(θ=π/2,r),
 Polar velocity \deleted{or evaporation velocity}  vθ(r,θ)=vθ=vθ1(θ)vθ2(θ=π/2,r),
 Azimuthal velocity  vϕ(r,θ)=vϕ=vϕ1(θ)vϕ2(θ=π/2,r),
 Fluid temperature  Θ(r,θ)=Θ=Θ1(θ)Θ2(θ=π/2,r), (8)

where the flow variables with ’’ in brackets are represented variation along the direction for a given and they are called as polar flow variables and corresponding derivatives will be polar flow derivatives. The flow variables with ’’ in brackets are represented variation along the radial direction and they are called as radial flow variables and corresponding derivatives will be radial flow derivatives. Here , we are following same as in the previous studies (Xue & Wang, 2005; Jiao & Wu, 2011) and corresponding radial derivative. Using above definitions in equations (1-5) then we get ordinary differential equations (ODEs) of the 2D flows,

 ρvrr[2r+1ρ2dρ2dr+1vr2dvr2dr]+ρvθ[1ρ1dρ1dθ+1vθ1dvθ1dθ+cotθ]=0 (9)
 vr1vrdvr2dr+vθvr2rdvr1dθ−v2θ+v2ϕr+Θ~tβρ2dρ2dr+Θ1~tβdΘ2dr−Fr=0, (10)
 vrrvθ1dvθ2dr+vθvθ2dvθ1dθ+vrvθ−v2ϕcotθ+Θβ~t1ρ1dρ1dθ+Θ2β~tdΘ1dθ=0 (11)
 rvϕ1vrdvϕ2dr+vϕ(vr+vθcotθ)+vθvϕ2dvϕ1dθ= Θ1vϕ1ρ2r2d(r2τrϕe)dr+α2vϕ2Θβ~tΩKr[τθΘ1dΘ1dθ+τθρ1dρ1dθ+d2vϕ1dθ2+vϕ1+τθcotθ] (12)
 vrΘ[NeffΘ2dΘ2dr−1ρ2dρ2dr]+vθΘr[NeffΘ1dΘ1dθ−1ρ1dρ1dθ]=β~tfQ+, (13)

where and is effective polytropic index. We have solved above equations (9-13) explicit way and following similar methodology as used in papers (Xue & Wang, 2005; Jiao & Wu, 2011). Since we are avoiding self-similar solution definitions along the radial direction, therefore firstly, we have to find out the radial flow variables with corresponding derivatives of the ADAF on the equatorial plane (detail equations are presented in appendix A), then get polar flow variables using symmetric boundary conditions on the equatorial plane and finally integrate above equations along the polar direction. Before doing so, we are assuming some symmetric properties with boundary conditions in the next subsection.

2.1 Boundary conditions for inflow-outflow

In order to solve ODEs (9-13) in direction, so we used symmetric boundary conditions at from the rotation axis. Which are obtained from the reflection symmetry and following the previous studies (Xue & Wang, 2005; Jiao & Wu, 2011),

 vθ1(π/2)=0=dρ1(π/2)dθ=dΘ1(π/2)dθ=dvr1(π/2)dθ=dvϕ1(π/2)dθ;  ρ1(π/2)=1. (14)

Since is an evaporation velocity for the generation of outflows so we assumed before the outflow at , it is zero but becomes non-zero immediately, when matter goes upward from the equatorial plane. Therefore we used is non-zero on the equatorial plane and represented below in equation (18). Moreover, the outflows are started from the equatorial plane so here we assumed total flow density at is equal to the inflow density means , which implies at all the radius. Here the is changing with the radius and also depends on and mass accretion rate () as represented in below equation (15). But the disk structure is independent of and . By using above definitions, we obtained explicitly the fluid equations in the pure radial direction at (appendix A). Since gas can evaporate from the accretion disk to infinity (Narayan & Yi, 1995a; Esin et al. , 1997; Gu, 2015), therefore we assumed mass loss in the continuity equation (A1), which defined as (Blandford & Begelman, 1999),

 ˙Min=−4πr2ρevrecosθe=˙Mb(rrb)s, (15)

where is exponent and called as the mass-loss parameter, is radial distance from the BHs when the disk started outflows from the equatorial plane and other quantities, and have denoted in appendix (A). According to the simulation paper by Ohsuga et al. (2005), is not a constant in the disk but average value has estimated around . Since there are limitations in the analytical approach, therefore, we assumed ’’ as a parameter and a constant for a particular solution. is the mass accretion rate at radius . Here, and is dimensionless mass accretion rate. is the Eddington mass accretion rate in the geometrical unit. Here, , corresponds to constant accretion rate means no mass loss from the disk. Since we want to study the outflows, therefore ’’ should be greater than zero. Now the equation (15) after differentiation can be written as,

 2r+1ρedρedr+1vredvredr=sr. (16)

Since we assumed that the radial components of flow variables and its derivatives are same for all values of polar angle at or above the equatorial plane for a particular radius. Therefore, the equation (9) with the help of the equation (16) becomes,

 vrs+vθ[1ρ1dρ1dθ+1vθ1dvθ1dθ+cotθ]=0. (17)

On solving\deletedfull model fluid equations (10-13) with equation (17), there are still many intricacies, so we made one more simplification on (the equatorial plane), i.e., we assumed double derivative of azimuthal velocity is vanishing as one of possibility since number of unknown problem the fluid equations (10-13, 17) with (14) at , we still need one more boundary condition in order to get flow variables. So we assumed from following as the equation (14). Thus, the polar flow variables on are estimated from equations (10-13, 17) with using equation (14) and after some simplifications, we get,

 aev2r1(90)+bevr1(90)−Fr=0,  Θ1(90)=x0x2vr1(90),  vϕ1(90)=√x3x4vr1(90)  and  dvθ1(90)dθ=−svrvθ2, (18)

where , ,  ,  ,  ,  ,  . Here, , , , and corresponding radial derivatives, , , , calculated from transonic ADAF solutions (Narayan et al. , 1997; Lu et al. , 1999) on the equatorial plane from equations (A10-A12). Here, subscript ’’ denotes values of the flow variables on the equatorial plane. In next section we will discuss solution procedure to find critical points and ADAF solutions.

3 Solution Procedure

Since the BH accretion is necessarily transonic because of the nature of gravity around central objects. Therefore, we first define and find out the critical point of the accretion flow in following subsections.

3.1 Critical point conditions

The critical point is a point of discontinuity of differential equation and mathematical is defined as form. So, the critical point conditions are obtained from the equation (A11),

 N=0 ⟹ (v2ϕe)crc+Frc+2(a2se)crc+fNeq(Λ+e)c=0 (19)

and

 D=0 ⟹ (v2re)c−(a2se)c=0. (20)

Here subscript ’’ denotes the flow quantities at the critical point and the radial velocity gradient at critical points obtained by lHospital rule. We found critical points by satisfying equations (19-20) together, with the help of integration of equations (A10-A12), for given set of parameters ( and ), detail explanations in appendix (B). We integrated the differential equations (A10-A12) from horizon to outward with the help of equation (A8). For this, we used a very nice technique for calculation of asymptotic flow variables very close to the horizon, which is describing in next subsection.

3.2 Method to find asymptotic flow variables

For hunting of the critical point location, we used a methodology as described in many papers (Becker et al. , 2008; Kumar & Chattopadhyay, 2013, 2014; Kumar et al. , 2014; Chattopadhyay & Kumar, 2016). Using Frobenius expansion for calculation of asymptotic value of for the differential equation (A12), the expression is

 λe=λ0+ζ(r−rS)Λ,   r→rS, (21)

where and are constants and to be determined by equation (A12) with using the equation (21), we get,

 limr→rSdλedr=limr→rS[2λer−γeffvreΩKζ(r−rS)Λα1a2s], (22)

Here, we assumed for limit and . Where, is free-fall velocity and . The value of will be obtained by iterations with satisfying the conditions (19, 20). With using expressions of and in the equation (22) then above equation can be written as,

 limr→rSγeffδζ(r−rS)Λα1a2se√2r(r−rS)3/2=2λ0rS (23)

When eliminating all terms from the above equation then we require . So, we get, . For a choice of value, we obtained a value of then we calculated flow variables very close to the horizon, say . Now we can obtain values of and at with the help of equations (21) and (A8) then we can integrate outward fluid equations (A10-A12) from . The detail method for finding critical points (CP) and disk structure are described in the appendix (B). \deletedWhen we combined equations (21) and (A8) with the value of and expression of at . Thus we got a polynomial in or . Now, supplying the parameters (or ), , , and then we solved the polynomial for at for first choice . Once obtained, other quantities and easily get with the help of the and equation (21). We now can integrate differential equations (A10-A12) outward from by using , and and simultaneously, checked the sonic point equations (19-20). If sonic conditions are not satisfied then we reduced the value of and repeat the whole procedure again from the solving polynomial for to checking sonic point conditions. This solution procedure repeated till satisfy sonic conditions. When ensuring it then we obtained critical point location () for given flow parameters. Once obtained, we integrated equations (A10-A12) outward along the radial direction for a given with other disk parameters. Then we investigated outer boundaries of ADAF solution (Narayan et al. , 1997; Lu et al. , 1999) by changing with repeating whole above procedure. Once obtained for the ADAF solution, we simultaneously integrated ODEs ( the radial equations A10-A12 and the polar equations 10-13 with 17) along the radial direction (inward and outward) and along the polar direction with using obtained radial flow variables and symmetric boundary conditions on the equatorial plane, we got complete 2D disk structure of the fluid flows.

4 Numerical results

We solved analytically 2D fluid equations with assuming explicitly radial fluid equations on the equatorial plane and first integrated along the radial direction, say at then immediately at same , we solved along the polar direction from the equatorial plane (details in the appendix B). Since there are many analytical studies on the 2D disk structure with using ADAF self-similar assumptions on the equatorial plane (Narayan & Yi, 1995a; Xu & Chen, 1997; Blandford & Begelman, 1999; Xue & Wang, 2005; Jiao & Wu, 2011). Therefore, we investigated the ADAF solutions on the equatorial plane for the calculations of the polar flow variables. So we would first represent the ADAF solutions in coming subsection and later next subsection with complete inflow-outflow solutions. In present work, we used both extreme values of or , one , where depends on , which may change the disk flow variables and structure with changing . And other , where for any value of . Here the mass inflow density and pressure of the gas have been calculated with the mass accretion rate and for all the solutions of this paper.

We used flow parameters to find the transonic accretion solutions on the equatorial plane, which are and .

In Figure (1), we represented typical ADAF solutions as previously shown by Narayan et al. (1997); Lu et al. (1999). Here, which are plotted with different values of in a first column, values of viscosity parameter () in a second column, and values of the grand specific energy () of the flow in a third column, which changes . Here is outer boundary of the ADAF or assumed transition radius from the Keplerian to the sub-Keplerian flows of the two zone configuration of the disk. The distribution of specific angular momentum (), bulk velocity () with sound speed () and the Bernoulli parameter are plotted in panels 1(a, d, & g), panels 1(b, e, & h) and panels 1(c, f, & i), respectively. In panel (1a), the value of is lowest for and increases with decreasing , when keeping other parameters are fixed. Since is lower for , therefore and are higher in panel (1b). The is lower for higher in panel (1c), since the is low, which is not compensated by high values of and . In the second column of the Figure (1), we changed and kept other parameters same. In panel (1d), values of is higher for . Since angular momentum transported less for lower viscosity when or is same. Therefore and are lower in panel (1e) and is also low in panel (1f), which is not compensated by higher . The third column of the Figure (1), we plotted curves with different means changing and keeping other parameters fixed. Here, the distribution of are almost same close to the BH and higher when approaching to for lower in panel (1g). As the expected variation of and are low for corresponding high around in panel (1h). The is lower for lower in panel (1i), so this may indicates for shorter , the possibility of the outflows may weak. Here, for different and are and , respectively. We get same power law scaling with the radius for and as in the Narayan et al. (1997) and almost independent of but the scaling for is changing significantly with (or ) as in the Figure (1h), e.g., for , for and for (as seen in Narayan et al. (1997)). Interestingly, these scaling rules are also same for the total and () on the equatorial plane, when calculating the 2D structures. Moreover, all the sub-Keplerian solutions have in the vicinity of the BHs but around have and mass density will be higher at since , so before transition radius, we believe that flow was Keplerian. Since the sub-Keplerian flows are showing positive in the intermediate values of , therefore may give rise outflows (Narayan & Yi, 1995a). Therefore we used these sub-Keplerian hot flows for the generation of the outflows and investigated the 2D disk structures as presented in next subsection.

4.2 Inflow-outflow solutions

Here outflow solutions above the equatorial plane in directions are calculated only up to a sonic surface when outflow Mach number () becomes equal to one, where is the total velocity of fluid and is a sound speed. Since the sonic surface rises a kind of discontinuity in the analytical integration of differential equations, therefore, integrations are invalid after the sonic surface without taking any proper methodology to solve the discontinuity. Therefore the fate of these outflows after the sonic surface is unknown in this study but we can predict that these disks may have strong outflows on the basis of transonic nature of the outflows and will reach very large after crossing the sonic surface. Here we used parameters and are non-zero when ODEs integrated along the polar directions. All the 2D disk figures with velocity vectors and density contours are plotted up to the radius size (), where outflows are started to generate from the disk. Here we are redefined a few flow variables in their physical units, e.g., the flow density , the gas pressure and flow temperature . Here, , , where and are accretion rate in unit of the Eddington accretion rate and mass of the BH in unit of the solar mass, respectively and .

In the Figure (2), the inflow matter formed the sonic surface very close to the BH () and hence matter enters into the BH supersonically. Moreover, the inflow region away from the BH () also becomes supersonic before the outflows are started to generate from the disk, which is marked by (open square) black symbol in the second and third columns of the Figure (2), except panel (f). But the inflow matter is always subsonic very close to the equatorial plane before inner critical point or sonic surface (). These supersonic regions are having locally higher density, lower , high and high rising , so total local velocity is high therefore arrow length is large in the same regions. The more detail of variations of the flow variables, we will present in Figure (5). Interestingly, the inflow matter is supersonic before making outflows. As the matter is moving inward, the flow variables are changed very fast and the flow becomes subsonic. This supersonic to subsonic transition of the inflow matter along the radial direction may give hint for the possibility of occurrence of shock transition in the flow, although this transition is not much sharp as accretion shocks studied in the literature (Fukue, 1987; Chakrabarti, 1989; Becker et al. , 2008; Kumar et al. , 2013, 2014; Chattopadhyay & Kumar, 2016; Lee et al. , 2016; Kumar & Chattopadhyay, 2017). The outflows occurred close to the BH because the thermal pressure and rotation velocity are increasing very fast and the local energy (as the profile of in the Figure 1c, f and i) becomes sufficient to generate bipolar outflows. We get disk surface with changing slope at every radius, which is unlike to the previous studies of 2D disk structure with self-similar assumptions (Jiao & Wu, 2011), they got the inflow disk surface with a constant slope. We also investigated no outflows, outflows and failed outflows regions of the 2D flow, which are depending on the disk parameters.

The above descriptions of the Figure (2) have based on the observations of the velocity fields and the outflow size. Now for more detail study of these figures, we plotted typical flow variables along direction with some fixed radius.

The Figure (3) is presented for three radii along the direction with same disk parameters, which are corresponding to the panel (2d). These are the typical outflow solutions with low viscosity values and we found variations of the flow variables are different from the previous analytical studies (Xu & Chen, 1997; Xue & Wang, 2005; Jiao & Wu, 2011) but have some basic qualitative similarities with the outflow solutions. The basic properties of the outflows are the radial velocity should becomes at some ’’ above the equatorial plane and other flow variables behavior may vary with boundary conditions on the equatorial plane. Here we get failed outflow solutions (solid red, ), which is multi-valued solution means has same value at various , outflow solution (dotted blue, ), which is also multi-valued and radially out outflow solution (dashed black, ) means at high latitude. Since these outflows are mainly driven by combinations of centrifugal force and gradient of pressures (gas or radiation) force, and behavior of temperature and angular velocity vary with the radial distance, therefore outflows and disk structure changed with the radius. If we see panel (3a), initially the for inflow and at some ’’ becomes zero, which gives the disk surface and gives outflow region. In panel (3b), the polar velocities () are increasing with decreasing ’’ but dashed black curve again turn back and approaches towards zero because is started to decreasing at high latitude (panel 3c). The is higher for lower radius solutions, this behavior maybe due to corresponding higher rotation velocity as in the panel (3c). The gas density (panel 3d) and pressure (panel 3e) are monotonically decreasing towards axis due to the expansion of gas above the equatorial plane. \deletedHere and are normalized with minimum values of and on the equatorial plane, respectively. The behavior of temperature is also not monotonic, it is decreasing and increasing toward axis in panel (3f) due to multi-valued nature of . And nature of is mostly depends on the and therefore on the viscosity. Although the dependence of multi-valued nature of velocities are very complicated but mostly depend on the viscosity, which we will see in the next figure. The solution corresponding to (solid red) is failed outflow (or fail to make transonic outflow solution and becomes again less than zero) due to very fast decreasing (panel 3c). Although is increasing but did not produce sufficient pressure gradient force to maintain positive and resulting matter falls back towards the BH. The solution corresponding to also feels more gravity than other solutions with higher values of . If we compare this kind of solution with no outflows analytical MHD solutions but for all (Zeraatgari et al. , 2018), in this case also decreases vary fast at high latitude. So plays a key role in generating the outflows. The black dashed curve behaves almost similar ways as the bipolar accretion outflow solutions are represented by Xu & Chen (1997).

The Figure (4) is drawn for fixed radius with different viscosity parameters and used same parameters corresponding to the second row of the Figure (2). In panel (4a), is high with high viscosity in the inflow region since is low (Figure 1d) therefore is also low in panel (4c). The is again high in the outflow region corresponding to same value of due to high acceleration for high viscosity, if we compare curves with solid red () and dotted blue (). Since is also high for high viscosity (panel 4b), therefore the total outflow velocity () is high in both curves. The (panel 4c) and (panel 4f) are monotonically increasing and decreasing, respectively for higher viscosity as compared with multi-valued curve for low viscosity (solid red). Since higher viscosity makes flow more hotter as is higher (panel 4f) and also transports more angular momentum therefore, somehow which makes smooth variation of (panel 4c). The solution corresponding to (dashed black) is not producing outflow due to low , as all three solutions have same gravity pull at . The (panel 4d) and (panel 4e) are decreasing smoothly with ’’. Here is higher for low viscous solution and therefore is also higher, since inflow is low. Interestingly, the gas density and pressure are not becoming zero at high latitude, specifically in high viscosity solutions because integrations are terminated at the outflow sonic surfaces.

Three curves of Figure (5) are plotted from three different regions of the Figure (2c), which are an inner (no outflow), middle (outflow) and outer (supersonic inflow ) regions. The solution for no outflow (solid red, ) region is close to the BH and experience more gravity and the combined fluid centrifugal force and pressure gradient force are not sufficient to defend gravity. So, the matter is not able to cross the disk surface and always, (panel 5a). The next middle with outflow region, solution (dotted blue, ) has appropriate forces, which make positive and give outflows. From outer part of the disk, the solutions from this region are showing very different behavior. If we look at dashed black curve (for ), as matter start expanding subsonically upward with increasing and along direction. Same time (panel 5a) and (panel 5f) are decreasing or increasing, simultaneously and resulting flow becomes supersonic and subsonic in the inflow region. In the same region is also increasing or decreasing (panel 5d). In this region, velocities, density and gas pressure gradients are high from the solutions of other two regions and also cooler. This kind of inflow supersonic regions above the equatorial plane as seen in the second and third columns of the Figure (2) except panel (f), which are surrounded by the black square symbol line. This kind of regions are not found with the low (panel 2d) or high (panel 2f) viscosity and low (panel 2g) solutions. Here the disk viscosity parameters roughly categorized as is low and is high.

The outflow solutions have a good qualitative agreement with radial self-similar adopted 2D HD (Xue & Wang, 2005; Jiao & Wu, 2011) and MHD (Samadi & Abbassi, 2016; Mosallanezhad et al. , 2016) flows. Since the disk vertical thickness of our solutions depends on the radius and flow parameters. So there is a possibility that the disk thickness can matched with previous studies for some suitable flow parameters. If we compare results for 2D disk structure with no outflow (Narayan & Yi, 1995a; Zeraatgari et al. , 2018) then the our disk vertical thickness is low since the integration is terminated at the sonic surface. For no outflow disk structure (Figure 2f), the thickness increased with the increasing viscosity. Moreover, the nature of flow variables profiles along the direction are mostly consistent with simulation by Yang et al. (2014).

effect of τθϕ and s

Now we are studying in this subsection effects on the 2D disk structure with variation of for and . Here we are taking three cases, one from the Figure (2a) which is plotted with , second, which has low viscosity with weak outflows from the Figure (2d) and third, high viscosity with no outflows from the Figure (2f). All the cases are presented in Figure (6) with variation of and keeping other parameters fixed for each cases.

The first case of the Figure (2a) is represented with two viscosity parameters, (panel 6a) and (panel 6b). The panel (6a) has disk structure with no outflows because more transfer of angular momentum due to high and as a result decreased much at high latitude and flow variables variation became similar to the dashed black curve of the Figure (4). So, all the matter will fall supersonically onto the BH after crossing the sonic surface. In the panel (6b) is plotted with , which gave outflow solutions and outflow region also increased from the Figure (2a). The second case from the Figure (2d) is represented with (panel c) and (panel d) of the Figure (6). When we compared with the Figure (2d), the outflow region and strength are increased in panel (6c) but decreased in panel (6d). In the panel (6c), is very small from the value of , so velocity vectors are almost parallel to the equator and the matter seems going back at high latitude. The last case of the Figure (2f) is again drawn with two different viscosities, (panel 6e) and (panel 6f). In both panels 6(e and f), outflow region is increased with decreasing from the Figure (2f). From this study, we can say that disk structure depends on and also depends on value of viscosity parameters. Since for low viscosity (), say , the outflows are increasing with increasing and for high viscosity (), say , the outflows are increasing with decreasing , when keeping fixed.

In Figure (7), we are represented variation of three velocities and flow temperature with at fix for different values, which are taken from the Figures (2f, 6e and 6f). For (solid red) is not showing the outflow. Since variation and values of is less (panel 7c) and is also decreasing above the equatorial plane (panel 7d), therefore combine effect of the outflows driving forces are not enough to make positive. When we decreased (dotted blue) and (dashed black) then the is high and rising faster at high latitude (panel 7c). So becomes positive and gives outflow. The similar behavior of has also found in the simulation with the inclusion of (Yang et al. , 2014), which decreases at high latitude. Moreover, (panel 7a) and (panel 7b) are increasing faster with decreasing , so the outflow strength is also increased.

Now, we are changed value of with keeping other parameters fixed and studied effects on the disk structure. Here we are taken case of the Figure (6e) with changing and presented in Figure (8).

In both panels of the Figure (8), the outflow region and strength are increased with increasing . Since the outflows are much affected by the viscosity parameters and mass-loss parameter. Therefore, we want to see the variation of local energy of the inflow-outflow and definition of the local energy of the flow is

 B(r,θ)=B=v2r2+v2ϕ2+v2θ2+h+Φ. (24)

This is modified Bernoulli energy parameter for the 2D flow. Which is similar to the local energy defined in the appendix (A) as the Bernoulli parameter on the equatorial plane, when is zero. Here, is specific enthalpy.