Inflow-outflow around black holes

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

Rajiv Kumar Department of Astronomy, Xiamen University, Xiamen, Fujian 361005, China Wei-Min Gu Department of Astronomy, Xiamen University, Xiamen, Fujian 361005, China Jiujiang Research Institute of Xiamen University, Jiujiang 332000, China

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
journal: ApJ

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 analytical study of the 2D disk with outflows has been started with a relaxation of vertical hydrostatic equilibrium in the disk by Narayan & Yi (1995a). They used self-similar ADAF solutions (Narayan & Yi, 1994) in the radial direction with symmetry conditions on the equatorial plane and solved the flow variables along the polar direction. However, they could not get actual outflow solutions because they assumed mass accretion rate is independent of radial distance and thus polar velocity, . Since they have found in their solutions that the Bernoulli parameter is positive, therefore the outflows may form close to a rotation axis. Subsequently, Xu & Chen (1997) have included the non-zero in their study and found accretion and ejection solutions. After that, the theoretical studies of the 2D disk with self-similar solutions along the radial directions have been done by many authors with the outflows (Xue & Wang, 2005; Gu et al. , 2009; Jiao & Wu, 2011; Gu, 2012, 2015) and without outflows (Habibi et al. , 2017). The similar studies have also done in magneto-hydrodynamics regime (Samadi & Abbassi, 2016; Mosallanezhad et al. , 2016; Zeraatgari et al. , 2018). An one more analytical study for the outflows with self-similar solutions in one-dimensional flow has done by Blandford & Begelman (1999) with assuming the mass accretion rate varies with some power of the radial distance. These outflow solutions are known as adiabatic inflow-outflow solutions (ADIOS). Further, they have presented their work with a family of two-dimensional self-similar solutions for the outflows (Blandford & Begelman, 2004). Simultaneously, many numerical simulations for the investigations of accretion-ejection have been also done with optically thin hot flows by Stone et al. (1999); Igumenshchev & Abramowicz (1999, 2000); Yuan et al. (2012a, b); Bu et al. (2016a, b), a review by Yuan & Narayan (2014) and with optically thick, hyper accreting flows by Ohsuga et al. (2005); Okuda et al. (2007); Yang et al. (2014); Jiang et al. (2014); Jiao et al. (2015); Jiang et al. (2017). There are a few more models for explanation of jet generation by the extraction of rotational energy of Kerr BHs (Blandford & Znajek, 1977), by anchoring of matter with magnetic lines (Blandford & Payne, 1982) and by shocked generated extra-thermal gradient force in the post-shock region (Molteni et al. , 1996a, b; Kumar & Chattopadhyay, 2013; Chattopadhyay & Kumar, 2016; Lee et al. , 2016; Kumar & Chattopadhyay, 2017).

The self-similar solution gets popularity for the analytical studies of the 2D disk with/without outflows in both hydrodynamics (HD) and magneto-hydrodynamics (MHD) regime because it simplifies fluid equations in the radial direction, which makes ODEs with only polar derivatives or independent of radius and radial derivatives. Gu (2012, 2015) have mentioned that the outflows form naturally from advective accretion disk for both optically thin and thick gas medium. Moreover, the hot flows with bremsstrahlung, synchrotron emissivity and Comptonization of soft photons can give rise to the hard spectrum of the BHCs (Chakrabarti & Titarchuk, 1995; Narayan & Yi, 1995b; Yuan & Narayan, 2014). Since the bipolar jets have been seen during hard state with radiatively inefficient flows around the BHCs (Remillard & McClintock, 2006). Therefore, we assumed inner part of the disk is hot advective and radiatively inefficient accretion flows for the study of jet generations around the BHs. We believed on the two zone configuration of the accretion disk (Esin et al. , 1997; Das & Sharma, 2013), one inner part hot sub-Keplerian advective radiative inefficient accretion flows (RIAFs) and other outer part, which is geometrical thin and cool Keplerian optically thick (Shakura-Sunyaev disk). For time being, we did not consider radiative emissivities in the flow since this work mainly focus on the study of jet generation in the 2D disk and full consideration of radiative advective flows leave for future work.

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,


the components of Navier-Stokes equation,


and the energy equation,


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,


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:


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,


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


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


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,


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,


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,


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




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


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


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,


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.

4.1 ADAFs solutions

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

Figure 1: Variations of the radial flow variables with radial distance, . Panels are showing variation of (a, d, g), (b, e, h) and (c, f, i). The panels (a-c) are plotted for parameters, with different (dotted red), (dashed black) and (long-dashed blue). The panels (d-f) are plotted for energy parameter with different (dotted, red), (dashed, black) and (long-dashed, blue). The panels (g-i) are plotted for viscosity parameter with different (dotted, red), (dashed, black) and (long-dashed, blue). The second and third columns are plotted for same and . The solid curve (cyan color) represented the Keplerian angular momentum distribution in panels (a, d, g).

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 .

Figure 2: The 2D disk structure represented with inflow density contour and velocity vector field, which is representing velocity direction of (or ) with a magnitude of (or ). A first, second and third rows are plotted for different , viscosity parameters () and , respectively. All the panels are plotted with corresponding same disk parameters of the Figure (1) with and . The first row has and other two rows have . Here solid black, dotted green, ‘+’ symbol blue and open square symbol black curves are representing disk surface, outflow sonic surface, inflow sonic surfaces and inflow supersonic region, respectively. Color bar at top of figures represents inflow density variation. Here the axes and density are having unit of and , respectively.

Figure (2) is representing the 2D disk structures with velocity vector fields and density contours. In a first row of the Figure (2) is plotted with different (panel 2a), (panel 2b) and (panel 2c), which are corresponding radial input solutions of the first column of the Figure (1). Disk thickness is high for the panel (2a) in the row. The disk thickness decreases with decreasing towards panel (2c), which is radiation-dominated for . Since the flow is accelerated more with decreasing , therefore, the gas formed the disk surface at lower latitude and the radiation dominated flows may have strong outflows. Moreover, the inflow matter becomes supersonic close to the horizon after ‘+’ symbol blue color line (where and ). Here, the disk surface ( and ) and outflow sonic surface ( and ) are represented with solid black line and dotted green line, respectively. The disk surface (solid line) separated the inflow and outflow regions in the disk. The velocity vectors are mostly directed towards the BH in the inflow region and in the outflow region, they are going out. Here, we found that the gas or radiation pressure dominated flows are having the outflows. Which are consistent with simulations for both the gas and radiation pressure supported flows. The case with advection-dominated and can resemble for high accretion rate flows with high luminous BH sources. The case with advection-dominated and can resemble for low accretion rate with low luminous BH sources. In a second row of the Figure (2) is plotted with different viscosity parameters (panel 2d), (panel 2e) and (panel 2f) for the input parameters corresponding to solutions of second column of the Figure (1). In panel (2d), which is plotted with we get two kind of the sonic surfaces above the disk surface (solid black), one for the outflow, where (dotted green) and other for the fail outflows (‘+’ blue), when the solutions fail to make sonic transition in the outflow and matter velocity again becomes . Corresponding to these two surface regions, the detail variations of flow variables along direction are presented below in Figure (3). In the second row, we increased the value of viscosity parameters and panel (2e) is plotted with then we got smooth inflow and outflow surfaces. The outflow strength and region are also increased. If we further increased the viscosity, (panel 2f) then we got almost no outflows with 2D disk inflow structure. Typical behavior of the flow variables of the second row panels are presented in Figure (4) with some fixed radius. In these cases, the outflows are strongly depends on the viscosity and with changing viscosity we can get no outflows, weak and strong outflows. They may explain various states of the BHCs, since the viscous time and cooling time scales are changed with changing viscosity or mass accretion rate parameter (Das & Sharma, 2013). Both the parameters are also changed the distribution of angular momentum in the flow (Kumar & Chattopadhyay, 2014), so the flowing matter can becomes Keplerian or sub-Keplerian. Now in last row, we are represented 2D structure corresponding to solutions of last column of the Figure (1). Here we found that the outflow region and strength depend on the transition radius and both are increased with increasing . The depends on the . For small , the outflow region is small and strength is also weak due to low local energy of the flow, when we compared with panels (2g), (2h) and (2i) in the row, which are having transition radius (), () and (), respectively. If we compare with the panel (2e), which has () and the panel (2i) then they give almost same disk structures because both have small flow energy difference. So, for the small advective disk (), the outflow region and strength are much affected by changing the ADAF disk size and accordingly the Keplerian disk size is also changed as explained by Esin et al. (1997). However in Esin et al. (1997) paper the ADAF size is decreased by increasing mass accretion rate, which also changed the distribution of local energy of the flow as similarly did it in the present paper.

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.

Figure 3: Variations of flow variables with polar angle, . Panels are showing variations of radial velocity, (a), evaporation velocity, (b), azimuthal velocity, (c), mass density, (d), gas pressure, (e) and temperature, (f). These solutions are drawn with different radii, (solid, red), (dotted, blue) and (dashed, black). All curves are drawn from the Figure (2d). Here all the velocities, , and are having units of , , and , respectively.

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

Figure 4: Variations of flow variables with polar angle, . Panels are showing variations of (a), (b), (c), (d), (e) and (f). These solutions are drawn for fixed radius, with different viscosity parameters, (solid, red), (dotted, blue) and (dashed, black) and other parameters are corresponding from the second row of the Figure (2).

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.

Figure 5: Variations of flow variables with polar angle, . Panels are showing variations of (a), (b), (c), (d), (e) and (f). These solutions are drawn with different radii, (solid, red), (dotted, blue) and (dashed, black) and the disk parameters are same from the panel (c) of Figure (2).

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

4.2.1 effect of and

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.

Figure 6: The density contours and velocity fields are plotted with same disk parameters corresponding to the first panel, panel (d) and panel (f) of the Figure (2). Here figures have plotted with different values of as mentioned in each panel.

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.

Figure 7: These solutions are plotted with changing viscosity parameter, (solid red), (dotted blue) and (dashed black) at same and other parameters are same as in the panel (f) of Figure (2).

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

Figure 8: The density contours and velocity fields are presented with different and other disk parameters are same as in the panel (e) of Figure (6). Here value of is mentioned in both panels with viscosity parameters.

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


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.

Figure 9: Variation of with (panel a) and (panel b) are presented. In the panel (a), curves are presented with fixed and other disk parameters are (solid red), (dotted blue) and (dashed black). In the panel (b), curves are presented with fixed and other disk parameters are same as the panel (a). Panels (c) and (d) are presented variations of with at fixed for parameters . The panel (c) is plotted with different, (solid red), (dotted blue) and (dashed black) for fixed . The panel (d) is plotted with different, (solid red), (dotted blue) and (dashed black) for fixed .

In Figure (9), we are presented the variations of outer boundary of the outflows () with (panel 9a), (panel 9b) and Bernoulli parameter with in panels 9(c and d) and other details are written in the caption. The two curves, solid line (red) and dotted line (blue) in the panel (9a) are represented with different flow constant of motion and , respectively. So, they are have different for and for . The curve with lower is has small outflow region and which more clear towards lower values of in the panel (9a). Here is higher means more matter going out from the disk or higher mass outflow rate. Again in the same panel (9a), another curve with dashed (black) line is plotted with the same as the curve dotted (blue) but both have different and . Other flow parameters are same for both the curves. We found that the outflows are high with higher for same . Since high rises more temperature and kinetic energy, therefore local specific energy of the flow is increased as seen in the Figure (1f). Here outflow region is increased with decreasing but disk thickness is decreased as seen in the Figure (6). For high , we did not find the outflows but has the inflow 2D disk structure as seen in the Figure (2f). In the panel (9b), the outflow region is increased with increasing and . The two curves, solid and dashed line are become maximum around and decreased with further increasing . Since for , the gas pressure variation along the radial direction becomes almost flat on the equatorial plane. In panels (9c), curves with variations of the are plotted for same solutions of the Figure (7). The solid red curve () is decreased with decreasing and gave no outflow solution as in the (2f). Other two curves ( and ), the are increased with decreasing at high latitude and gave outflows. In panel (9d), the