A 2.5-dimensional viscous, resistive, advective magnetized accretion-outflow coupling in black hole systems: a higher order polynomial approximation

Outflows and jets are ubiquitous in nature. They are observed both in the local universe, mostly in black hole (BH) X-ray binaries (BHXRBs) which are believed to harbor stellar mass BHs called microquasars (Mirabel & Rodríguez 1994, 1998; Eikenberry et al. 1998), and in powerful extragalactic radio sources (Blandford & Rees 1974; Begelman et al. 1984; Mirabel 2003) where well-collimated outflows or jets emerge continuously from the nuclear region of their host active galaxies (active galactic nuclei, AGNs) or quasars harboring supermassive BHs. The accreting hot plasma around BHs powered by extreme gravity of the central object results in the formation of an outflow/jet which extracts mass, angular momentum and energy from the inner regions of the accretion flow. The outflows in microquasars are observed only in the low-hard state of BHXRBs (Fender et al. 2004; Rushton et al. 2010) which are radiatively inefficient. Radiatively inefficient accretion flows (RIAFs) are hot gas pressure dominated systems which are geometrically thick (h(r)/r≳0.1) and optically thin, where h(r) is the scale-height of the accretion region, with the dynamics of the flow strongly sub-Keplerian and advection dominated (Narayan & Yi 1994, 1995; Yuan et al. 2003). RIAFs occur when the mass accretion rate $\dot{M}$ is very low (presumably with $\dot{M}\,\lesssim {10}^{-3}\,{\dot{M}}_{\rm{Edd}}$), where ${\dot{M}}_{\rm{Edd}}$ is the Eddington accretion rate or the accretion rate corresponding to the Eddington luminosity. Outflows/jets are not likely to be observed in the high-soft state of BHXRBs (Rushton et al. 2010 and references therein), which are believed to be powered by a geometrically thin and optically thick radiation pressure dominated Keplerian accretion disk (Shakura & Sunyaev 1973; Frank etal. 2002).

Theoretically speaking, it has been argued (Narayan & Yi 1994; Ghosh & Mukhopadhyay 2009; Ghosh et al. 2010; Bhattacharya et al. 2010) that a geometrically thick advective accretion flow has a strong tendency to drive bipolar outflows due to high thermal energy content of the hot gas. They may be additionally propitious to propel outflows/jets because their vertical thick structure enhances the large-scale poloidal component of the magnetic field, which plays a critical role in launching strong and collimated outflows (Meier 1999; Ballantyne & Fabian 2005).

Apart from the low-hard state of BHXRBs which power jets, at the other end of the spectrum, strong outflows and jets are observed in low excitation radio galaxies (LERGs) harboring supermassive BHs. LERGs, a more general case of low luminous AGNs (LLAGNs), seem to be accreting gaseous plasma directly from the hot X-ray emitting phase of the interstellar medium (ISM) or from the hot X-ray halos surrounding the galaxy or from the hot phase of the intergalactic medium (IGM) quasi-spherically in a radiatively inefficient mode with a near Bondi rate (Allen et al. 2006; Hardcastle et al. 2006, 2007; Balmaverde et al. 2008; Ghosh & Banik 2015). LERGs thus resemble the low-hard state of BHXRBs, having geometrically thick and optically thin gas pressure dominated advective accretion flow, accreting hot gas at a high sub-Eddington accretion rate. This strongly advective radiatively inefficient accretion paradigm (RIAF) or hot mode accretion having considerable geometrical thickness is more prone to emanate outflows/jets and is very conducive to propel matter vertically outwards out of the accreting region.

However, with the increase in $\dot{M}$ as ${10}^{-3}{\dot{M}}_{\rm{Edd}}\ll \dot{M}\,\lesssim {10}^{-2}\,{\dot{M}}_{\rm{Edd}}$, the flow tends to be more centrifugally dominated and becomes moderately advective, with the central BH accreting relatively cold gas as compared to the hot mode accretion. Incidentally, the moderately ad-vective accretion flow does not occur in a radiatively inefficient mode, but where a considerable amount of both gas and radiation pressure seem to be present in the system, rendering the flow to have a moderate optical depth. Geometrically, the inner advective region would then be relatively thinner than that corresponding to RIAFs. The moderately advective accretion paradigm may also be susceptible to ejecting outflows. The difference between this paradigm and that with RIAF, however, may lie in the acceleration mechanisms to eject bipolar outflows and jets, which we would eventually investigate in this study. Nonetheless, with the increase in $\dot{M}$ as the flow tends to become more rotationally/centrifugally dominated with lesser geometrical thickness, the efficacy of the disk to eject outflows diminishes. Beyond $\dot{M}\gt {10}^{-2}{\dot{M}}_{\rm{Edd}}$, the flow would eventually tend towards a Keplerian nature (Yuan & Narayan 2014; paper II (in preparation)). A geometrically thin Keplerian accretion disk plausibly fails to account for the launching and acceleration of outflows and jets (Ghosh & Mukhopadhyay 2009; Ghosh et al. 2010; Bhattacharya et al. 2010; also see Czerny & You 2016 for further discussion).

Extensive work has been pursued on the origin of outflows/jets, since the seminal work of Blandford & Payne (1982) in studying accretion powered hydro-magnetic outflows, which we focus upon in this study. Physical understanding of accretion powered hydromagnetic outflows has either been performed in the stationary self-similar approximation (Blandford & Payne 1982; Pudritz & Norman 1986; Wardle & Koenigl 1993; Contopoulos 1996; Ferreira & Pelletier 1995; Ferreira 1997; Casse & Ferreira 2000; Narayan et al. 2007) in the quasi-analytical regime to demonstrate the importance of the poloidal component of the magnetic field to launch outflowing matter from the Keplerian accretion disk, or through magnetohydrodynamic (MHD) simulations in both nonrelativistic as well as in relativistic regimes (for details, see the introduction in Ghosh & Mukhopadhyay 2009; Bhattacharya et al. 2010 and references therein; also see Pudritz et al. 2007 and references therein).

In most of these studies the authors remain focused mainly on the launching of outflows/jets from the geometrically thin Keplerian disk. The formation of the accretion powered outflow and jet is directly related to the efficacy of extraction of angular momentum and energy from the magnetized accretion flow. However, the exact mechanism by which the radial accretion is diverted into strong outflows and plausible jets still remains theoretically elusive. Notwithstanding, jet launching is completely an MHD process. Accreting material diffuses across magnetic field lines threading the accretion region, and is then lifted upwards by MHD forces which then couples to the field and becomes accelerated magnetocentrifugally. However, if the accreting system is strongly gas pressure dominated, it may happen that the gas pressure gradient would play a more contributory role in lifting the plasma vertically outwards along with the help of magnetic forces. In addition, turbulent reconnection of magnetic field lines may lead to flux annihilation (Bisnovatyi-Kogan & Lovelace 2000). Magnetic energy dissipates through turbulent magnetic reconnection, which may also power the outflow/jet emission (de Gouveia dal Pino & Lazarian 2005; de Gouveia Dal Pino et al. 2010; Giannios et al. 2010 and references therein).

Most studies on accretion disk and outflow/jet have evolved separately, assuming these two to be apparently dissimilar objects. In light of both deeper theoretical understanding and observational inferences (Merloni et al. 2003; Fender & Belloni 2004; Fender et al. 2004), it is evident that the dynamics of outflow and the underlying accretion are strongly correlated (for details see Ghosh & Mukhopadhyay 2009; Ghosh et al. 2010; Bhattacharya et al. 2010 and references therein). Outflows and jets observed in AGNs and X-ray binaries (XRBs) can only originate in an accretion powered system, where the accreting plasma around gravitating objects like BHs acts as a source, whereas the outflow and then the jet acts as one of the possible sinks (Bhattacharya et al. 2010). The implicit coupling between accretion and outflow is then essentially governed by conservation laws; conservation of matter, momentum and energy. The outflowing matter carries away mass, angular momentum and energy extracted from the accreting plasma (Pudritz & Norman 1986). We do not intend to investigate the physics of jet formation and its launching mechanism which is altogether a different field of research, but would like to focus entirely on the inter-correlating dynamics of the accretion and outflow within the coupled accretion-outflow region, and the conditions/criteria for jet launching. Any proper understanding of the dynamics and the conditions of jet launching should necessarily require a robust understanding of the dynamics of the magnetized advective accretion region coupled to outflow, governed explicitly by conservation laws. The relevant dynamical solutions at the accretion-outflow coupled surface (the surface from where outflow decouples from radially inward accretion flow) would then act as boundary conditions at the base of the jet.

Accretion disk-outflow/jet coupling has been studied on a number of occasions from an observational angle, in accretion powered systems (Fender & Belloni 2004; Yuan et al. 2005; Körding et al. 2007; Neilsen & Lee 2009; Soleri et al. 2010; Miller-Jones et al. 2011, 2012; Miller et al. 2013; King et al. 2013). From a theoretical perspective, few self-similar studies have been attempted in the context of accretion-outflow coupling, both in a non-magnetized and in a magnetized regime (see Pudritz et al. 2007 and references therein; Ghosh & Mukhopadhyay 2009; Ghosh et al. 2010; Bhattacharya et al. 2010; Blandford & Begelman 1999; Soria et al. 1997; Samadi & Abbassi 2016; Xue & Wang 2005; Jiao & Wu 2011; Jiao et al. 2015). Notwithstanding, theoretically, it is still difficult to construct a reasonably satisfactory and definitive model of the magnetized accretion-outflow/jet coupled region, owing to the complicated geometry and inconclusive understanding of the inflow-outflow coupled region. On the other hand, few simulations on disk-outflow/jet coupling have also been peformed (Nishikawa et al. 2005; McKinney & Narayan 2007b,a). In these simulations, how the matter gets deflected from the equatorial plane has been studied largely in the Keplerian regime. Casse & Keppens (2004) performed an advective, resistive MHD simulation of the accretion-ejection structure with the inclusion of the energy equation, however neglecting the viscosity and radiative loss from the system. Nonetheless, it is still difficult to simultaneously simulate the accretion and outflow regions because the time scales of accretion and outflow are in general very different.

In one of our previous works on accretion-induced outflow (Ghosh & Mukhopadhyay 2009), we made an attempt to formulate a unifying scheme of inflow and outflow in a 2.5-dimensional advective paradigm. In the present work, we endeavor to develop a robust viscous, resistive and advective MHD accretion-induced outflow model in the 2.5-dimensional paradigm, restricting ourselves to the non-relativistic regime, in the context of accretion powered hydromagnetic outflows/jets, focusing entirely on the inter-correlating dynamics of the accretion and outflow within the coupled accretion-outflow region, without aspiring to explore the mechanism of launching and ejection of outflows and jets. A complete 2.5-dimensional viscous, resistive, advective global MHD numerical solution of such a system is left for further work, which is beyond the scope at present; we confine our treatment to quasi-analytical/quasi-numerical power law self-similarity (e.g., Narayan & Yi 1994) in a quasi-stationary configuration, by upholding the conservation equations. All physical quantities are scaled as powers in r and z according to their dimensions, in the limit of higher order polynomial expansion. We will perform a critical analysis of the accretion-outflow region and provide a complete quasi-analytical family of solutions. Although the quasi-analytical self-similar solutions are approximate, they, however, can provide a strong physically intuitive picture of accretion dynamics coupled with the outflow, as well as physical criteria/conditions to eject outflows/jets.

In the next section, we present the formulation of our model. Section 3 describes the quasi-analytical/quasi-numerical procedure to solve the model equations of the accretion-induced outflow. In Section 4, we evaluate the coefficients of our self-similar solutions and analyze them. In Section 5, we investigate the nature and behavior of the family of solutions for accretion-induced outflow within the bounded accretion coupled outflow region. Finally, we end up in Section 6 with a summary and discussion.

We formulate the accretion-outflow coupled model by considering a 2.5-dimensional viscous, resistive, advective accretion flow geometry in the MHD regime, in a quasi-stationary state. The vertical flow is explicitly included in the system. We adopt the cylindrical coordinate system (r, φ, z) to describe a quasi-stationary, mean axisymmetric accretion flow. As we have incorporated outflow in our system, within the accretion-outflow coupled region in an advective paradigm, all the dynamical flow variables; namely, radial velocity (v_r), azimuthal velocity (v_φ), specific angular momentum (λ), vertical velocity or outflow velocity (v_z), isothermal sound speed (c_s), mass density (ρ), thermal pressure (P) and magnetic field components (B_r, B_z, B_φ), vary in both r and z. The dynamical equations are vertically integrated over an arbitrary scale-height h(r) from −h to +h. Here h(r) is not a hydrostatic disk-scale height but a photospheric surface within which accretion and outflow are coupled. Above h(r), the outflow decouples from the accretion flow, gets further accelerated in the hot nonthermal magnetized corona and finally forms a relativistic jet. We focus on this accretion-outflow coupled region within which the flow is mostly bounded. We only include the rφ component of turbulent stress, which is responsible for radial transport of angular momentum outwards (angular momentum gets transported due to the diffusion of turbulent eddies). Vertical transport of angular momentum occurs mainly through large-scale magnetic stresses, where the outflowing matter magnetically extracts or removes angular momentum. The accreting mass is assumed to be much less compared to that of the central object, and hence the flow is not self-gravitating. As the accretion flow is associated with turbulence, we express all the dynamical variables in terms of the mean and fluctuating parts, generically represented as

$$\begin{eqnarray}&&F=\overline{F}+{F}^{^{\prime} },\end{eqnarray} \tag{ 1 }$$

where $\bar{F}$ is then the appropriate mean value and F' is the fluctuation corresponding to that variable. We adhere to the conventional statistical averaging approach of Favre (Favre decomposition) used for compressible turbulent flows, where the mean value of the flow velocity is actually the density weighted mean velocity; and corresponding to the rest of the variables like density, pressure and magnetic field, the mean values are the conventional mean. This will then define a mean flow in the disk, plausibly associated with some kind of time average; for a quasi-stationary flow, one can conceive of a time average, with time scale much larger than the time scales of turbulent fluctuations. As the turbulence dominates the transport process in the accretion flow, we neglect microscopic viscosity and resistivity, in a usual way. Statistical averaging of MHD equations generates a large number of turbulent correlation terms. In the present study, we restrict ourselves to first order turbulent correlation and neglect second order and higher order correlation terms.

The turbulence is defined in terms of mean Reynolds and Maxwell stress described through correlations given by

$$\begin{eqnarray}\begin{array}{ll}{\bar{t}}_{ij} & =\overline{{t}_{ij}^{ {\mathcal R} }}+\overline{{t}_{ij}^{ {\mathcal M} }}\\ & \Rightarrow -[\overline{\rho {v}_{i}^{^{\prime} }{v}_{j}^{^{\prime} }}- \ (\overline{\frac{{B}_{i}^{^{\prime} }{B}_{j}^{^{\prime} }}{4\pi }- \ {\delta }_{ij}\frac{{B}^{^{\prime} }{}^{2}}{8\pi }})],\end{array}\end{eqnarray} \tag{ 2 }$$

where ${\bar t}_{ij}$ is the net turbulent stress. The Reynolds and Maxwell stresses are commonly parameterized in terms of kinetic and magnetic turbulent viscosities, ${\nu }_{ij}^{ {\mathcal R} }$ and ${\nu }_{ij}^{ {\mathcal M} }$ respectively, as

$$\begin{eqnarray}&&\overline{\rho {v}_{i}^{^{\prime} }{v}_{j}^{^{\prime} }}=-{\nu }_{ij}^{ {\mathcal R} }\overline{\rho }\overline{{s}_{ij}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&-(\overline{\frac{{B}_{i}^{^{\prime} }{B}_{j}^{^{\prime} }}{4\pi }-{\delta }_{ij}\frac{{B}^{^{\prime} }{}^{2}}{8\pi }})=-{\nu }_{ij}^{ {\mathcal M} }\overline{\rho }\overline{{s}_{ij}},\end{eqnarray} \tag{ 4 }$$

where ${\bar{s}}_{ij}=\frac{\partial {\bar{v}}_{i}}{\partial {x}_{j}}+\frac{\partial {\bar{v}}_{j}}{\partial {x}_{i}}-\frac{2}{3}\nabla \cdot \bar{v}{\delta }_{ij}$ is the strain tensor. The turbulent viscosities are parameterized through an α prescription as

$$\begin{eqnarray}&&{\nu }_{ij}^{ {\mathcal R} }\sim {\alpha }_{ij}^{ {\mathcal R} }{\bar{c}}_{\rm{s}}h,\,\,\,\,\,\,\,\,{\nu }_{ij}^{ {\mathcal M} }\sim {\alpha }_{ij}^{ {\mathcal M} }{\bar{c}}_{\rm{s}}h,\end{eqnarray} \tag{ 5 }$$

where ${\nu }_{ij}={\nu }_{ij}^{ {\mathcal R} }+{\nu }_{ij}^{ {\mathcal M} }$ is the net turbulent viscosity and ${\alpha }_{ij}={\alpha }_{ij}^{ {\mathcal R} }+{\alpha }_{ij}^{ {\mathcal M} }$ is the net turbulent viscosity parameter. With these parameterizations, the coupled MHD accretion-outflow dynamical equations for mean flow in quasi-stationary state are as follows:

• (a)  
  Mass transfer:

  $$\begin{eqnarray}&&\frac{\partial (\bar{\rho }{\bar{v}}_{j})}{\partial {x}_{j}}=0.\end{eqnarray} \tag{ 6 }$$

  We define the net mass flow rate, which is a constant, through an integro-differential equation as

  $$\begin{eqnarray}&&\displaystyle {\int }_{-h}^{+h}\displaystyle {\int }_{r}\displaystyle {\int }_{0}^{2\pi }\left[\frac{1}{r}\frac{\partial }{\partial r}(r\bar{\rho }{\bar{v}}_{r})+\frac{\partial }{\partial z}(\bar{\rho }{\bar{v}}_{z})\right]rd\varphi drdz=-\dot{M}\end{eqnarray} \tag{ 7 }$$

  where the first term is the signature of the radial accretion flow and the second term is attributed to outflow. If we discard ${\bar{v}}_{z}$ (neglecting outflow), Equation (7) reduces to a height integrated continuity equation of the accretion flow, and where $\dot{M}$ would then be the usual mass accretion rate.
• (b)  
  Momentum transfer:The momentum balance equation in the mean field MHD is given by

  $$\begin{eqnarray}\begin{array}{ll}\frac{\partial \left(\bar{\rho }{\bar{v}}_{i}{\bar{v}}_{j}\right)}{\partial {x}_{j}}= & -\bar{\rho }\frac{\partial {\varphi }_{G}}{\partial {x}_{i}}-\frac{\partial \bar{p}}{\partial {x}_{i}}\\ & +\frac{\partial }{\partial {x}_{j}}\left(\frac{{\bar{B}}_{i}{\bar{B}}_{j}}{4\pi }-{\delta }_{ij}\frac{{\bar{B}}^{2}}{8\pi }\right)\\ & +\frac{\partial }{\partial {x}_{j}}\left(\overline{{t}_{ij}^{ {\mathcal R} }}+\overline{{t}_{ij}^{ {\mathcal M} }}\right).\end{array}\end{eqnarray} \tag{ 8 }$$

  Using Equation (6) and integrating Equation (8) vertically, the radial momentum balance equation is given by

  $$\begin{eqnarray}&&\displaystyle {\int }_{-h}^{+h}\left[\bar{\rho }{\bar{v}}_{r}\frac{\partial {\bar{v}}_{r}}{\partial r}-\bar{\rho }\frac{{\bar{\lambda }}^{2}}{{r}^{3}}+\bar{\rho }{\bar{v}}_{z}\frac{\partial {\bar{v}}_{r}}{\partial z}\right.\\ &&+\bar{\rho }{F}_{\rm{Gr}}+\frac{\partial \bar{P}}{\partial r}+\frac{1}{4\pi }\left(\frac{\bar{{B}_{\varphi }}}{r}\frac{\partial }{\partial r}(r\bar{{B}_{\varphi }}\right)\\&&+\left.\bar{{B}_{z}}\frac{\partial \bar{{B}_{z}}}{\partial r}-\bar{{B}_{z}}\frac{\partial \bar{{B}_{r}}}{\partial z})\right]dz=0.\end{eqnarray} \tag{ 9 }$$

  In deriving this we have used divergence criteria of the magnetic field. F_Gr is the radial component of the gravitational force. In a similar fashion, we write the azimuthal momentum balance equation as

  $$\begin{eqnarray}&&\begin{array}{l}\displaystyle {\int }_{-h}^{+h}\left(\bar{\rho }\frac{{\bar{v}}_{r}}{r}\frac{\partial \bar{\lambda }}{\partial r}+\bar{\rho }\frac{{\bar{v}}_{z}}{r}\frac{\partial \bar{\lambda }}{\partial z}\right)dz=\\ \frac{1}{{r}^{2}}\frac{\partial }{\partial r}\displaystyle {\int }_{-h}^{+h}\left({r}^{2}{\bar{t}}_{r\varphi }\right)dz\\ +\displaystyle {\int }_{-h}^{+h}\frac{1}{{r}^{2}}\frac{\partial }{\partial r}\left(\frac{{r}^{2}\bar{{B}_{r}}\bar{{B}_{\varphi }}}{4\pi }\right)dz+\frac{\bar{{B}_{\varphi }}\bar{{B}_{z}}}{4\pi }{|}_{-h}^{+h}.\end{array}\end{eqnarray} \tag{ 10 }$$

  The last term in the right hand side of Equation (10) is the magnetic torque that acts at the accretion-outflow surface, and which helps in transporting the angular momentum vertically outwards. This term is responsible for mass loss in the wind. Next we derive the vertical momentum balance equation which is obtained from the z component of Equation (8).

  $$\begin{eqnarray}&&2\displaystyle {\int }_{0}^{h}\left[\bar{\rho }{\bar{v}}_{r}\frac{\partial {\bar{v}}_{z}}{\partial r}+\bar{\rho }{\bar{v}}_{z}\frac{\partial {\bar{v}}_{z}}{\partial z}+\bar{\rho }{F}_{Gz}+\frac{\partial \bar{P}}{\partial z}\right.\\ &&+\left.\frac{\partial }{\partial z}\left(\frac{{\bar{B}}_{\varphi }^{2}+{\bar{B}}_{r}^{2}}{8\pi }\right)-\frac{{\bar{B}}_{r}}{4\pi }\frac{\partial {\bar{B}}_{z}}{\partial r}\right]dz=0,\end{eqnarray} \tag{ 11 }$$

  where F_Gz is the vertical component of gravitational force. Equation (11) is integrated from 0 to h due to the reflection symmetry of all the dynamical variables across the accretion-outflow coupled surface. Equation (11) contains information on the outflow dynamics within the accretion-outflow region. The matter starts to accelerate vertically outwards from just above the equatorial plane of the accretion region. If there is no outflow, then v_z = 0, and if we neglect the magnetic pressure and magnetic stresses, Equation (11) reduces to the well known hydrostatic equilibrium condition in the disk, and the usual hydrostatic disk scale-height can be obtained.
• (c)  
  Divergence condition:

  $$\begin{eqnarray}&&\displaystyle {\int }_{-h}^{+h}\left[\frac{1}{r}\frac{\partial }{\partial r}(r{\bar{B}}_{r})+\frac{\partial {\bar{B}}_{z}}{\partial z}\right]dz=0.\end{eqnarray} \tag{ 12 }$$

  The divergence condition determines the symmetry property of magnetic field components. Whether the radial and vertical components of the magnetic field will follow odd and even symmetry or vice-versa in the z direction can be ascertained from the above equation.
• (d)  
  Magnetic induction:The turbulent magnetic induction equation is derived from the mean field MHD theory (e.g., Krause & Raedler 1980). Following the usual procedure and neglecting the dynamo effect, the steady state induction equation in tensorial form is given by

  $$\begin{eqnarray}&&{\epsilon }_{ijk}{\epsilon }_{kmn}\frac{\partial }{\partial {x}_{j}}({\bar{v}}_{m}{\bar{B}}_{n})-{\epsilon }_{ijk}\frac{\partial }{\partial {x}_{j}}\left({\epsilon }_{kmn}{\nu }_{ml}^{ {\mathcal R} }\frac{\partial {\bar{B}}_{n}}{\partial {x}_{l}}\right)\,=0,\end{eqnarray} \tag{ 13 }$$

  where ${\nu }_{ml}^{ {\mathcal R} }$ is the kinetic part of the turbulent viscosity. The above equation has been written in the most general form considering an anisotropic turbulence. Splitting Equation (13) into radial and azimuthal directions, we obtain after vertical integration

  $$\begin{eqnarray}&&\begin{array}{l}\displaystyle {\int }_{-h}^{+h}\left[\frac{\partial }{\partial z}({\bar{v}}_{r}{\bar{B}}_{z}-{\bar{v}}_{z}{\bar{B}}_{r})\left.\\ +\right.\frac{\partial }{\partial z}\left({\nu }_{zz}^{ {\mathcal R} }\frac{\partial {\bar{B}}_{r}}{\partial z}-{\nu }_{rr}^{ {\mathcal R} }\frac{\partial {\bar{B}}_{z}}{\partial r}\right)\right]dz=0,\end{array}\end{eqnarray} \tag{ 14 }$$

  and

  $$\begin{eqnarray}&&\begin{array}{r}\displaystyle {\int }_{-h}^{+h}\left[\frac{\partial }{\partial z}({\bar{v}}_{\varphi }{\bar{B}}_{z}-{\bar{v}}_{z}{\bar{B}}_{\varphi })-\frac{\partial }{\partial r}({\bar{v}}_{r}{\bar{B}}_{\varphi }\right.\\ \left.-{\bar{v}}_{\varphi }{\bar{B}}_{r})+\frac{\partial }{\partial z}\left({\nu }_{zz}^{ {\mathcal R} }\frac{\partial {\bar{B}}_{\varphi }}{\partial z}\right)\right.\\ +\left.\frac{\partial }{\partial r}\left({\nu }_{rr}^{ {\mathcal R} }\frac{\partial {\bar{B}}_{\varphi }}{\partial r}+{\nu }_{\varphi \varphi }^{ {\mathcal R} }\frac{{\bar{B}}_{\varphi }}{r}\right)\right]dz=0.\end{array}\end{eqnarray} \tag{ 15 }$$

  We do not show the vertical component of the induction equation here as it is similar to that of the radial equation, and contains the same information regarding magnetic dynamics of the flow. Note that turbulent diffusion in the induction equation arises only from the kinetic part of the turbulent stress tensor through Reynolds stress. This is attributed to the mean field approximation, where we split the quantities into the mean and turbulent parts.
• (e)  
  Energy conservation:The Poynting flux S_j is given by

  $$\begin{eqnarray}&&{S}_{j}=\frac{c}{4\pi }{\epsilon }_{jlm}{E}_{l}{B}_{m}.\end{eqnarray} \tag{ 16 }$$

  Using mean field MHD and the induction equation, and by discarding microscopic resistivity, we write

  $$\begin{eqnarray}\begin{array}{ll}\frac{\partial {\bar{S}}_{j}}{\partial {x}_{j}}= & \frac{{\bar{v}}_{j}}{c}{\epsilon }_{jlm}{\bar{B}}_{l}{\bar{J}}_{m}+\frac{{\bar{J}}_{j}}{c}{\epsilon }_{jlm}\bar{{v}_{l}^{^{\prime} }{B}_{m}^{^{\prime} }}\\ & +\frac{{\bar{v}}_{j}}{c}{\epsilon }_{jlm}\bar{{B}_{l}^{^{\prime} }{J}_{m}^{^{\prime} }}+\frac{{\bar{B}}_{j}}{c}{\epsilon }_{jlm}\bar{{J}_{l}^{^{\prime} }{v}_{m}^{^{\prime} }},\end{array}\end{eqnarray} \tag{ 17 }$$

  where J is the current density. In contrast to the other correlation terms, the last term in the right hand side of Equation (17) is not a first order correlation term, but rather a second order or higher order correlation term (Pessah et al. 2006a); we omit this term in our present study. Presuming that the major contribution to heat generation in the accretion disk is due to turbulent stress, which can be considered to be reasonably true in a global sense (as assumed in the standard picture of accretion disk theory), we write the energy conservation equation, using Equation (17) and neglecting the kinetic dynamo effect, as

  $$\begin{eqnarray}&&\begin{array}{l}\frac{\partial }{\partial {x}_{j}}\left[\bar{\rho }{\bar{v}}_{j}\left(\frac{{\bar{v}}^{2}}{2}+\frac{8-3\beta }{2\beta }\frac{{\bar{P}}_{g}}{\bar{\rho }}+{\varphi }_{G}\right)-{\bar{v}}_{i}\bar{{t}_{ij}^{ {\mathcal R} }}\right]\\ \,\,\,\,\,\,\,\,\,\,\,\,+\frac{{\bar{v}}_{j}}{c}{\epsilon }_{jlm}{\bar{B}}_{l}{\bar{J}}_{m}+\frac{{\bar{v}}_{j}}{c}{\epsilon }_{jlm}\bar{{B}_{l}^{^{\prime} }{J}_{m}^{^{\prime} }}\\ \,\,\,\,\,\,\,\,\,\,\,\,-\frac{1}{4\pi }{\epsilon }_{jlm}{\nu }_{lk}^{ {\mathcal R} }\frac{\partial {\bar{B}}_{m}}{\partial {x}_{k}}{\epsilon }_{jlm}\frac{\partial {\bar{B}}_{m}}{\partial {x}_{l}}+\frac{\partial {F}_{j}^{r}}{\partial {x}_{j}}=0,\end{array}\end{eqnarray} \tag{ 18 }$$

  where $\beta ={\bar{P}}_{g}/({\bar{P}}_{g}+{\bar{P}}_{r})$; the ratio of gas pressure to the total pressure in the accretion-outflow coupled region is assumed to be constant. ${\bar{P}}_{r}$ is radiation pressure in the gas-radiation mixture. The radiation field has been assumed to be locally isotropic. The ratio of specific heat of the fully ionized gas and radiation are taken as γ_g = 5/3 and γ_r = 4/3 respectively. The 'effective ratio of specific heat' Γ is related to β through Γ = (8 – 3β)/(6 – 3β). In deriving Equation (18), we have neglected the turbulent thermal conductivity. The last term in the above equation represents transport of radiative flux. Multiplying Equation (8) by ${\bar{v}}_{j}$ and using the continuity equation, the vertically integrated energy budget for accretion-induced outflow is thus obtained below

  $$\begin{eqnarray}&&\begin{array}{l}\displaystyle {\int }_{-h}^{+h}[\frac{3}{2}(2-\beta )\bar{\rho }{\bar{v}}_{r}\frac{\partial {\bar{c}}_{\rm{s}}^{2}}{\partial r}-{\bar{v}}_{r}{\bar{c}}_{\rm{s}}^{2}\frac{\partial \bar{\rho }}{\partial r}]dz\\ +\displaystyle {\int }_{-h}^{+h}[\frac{3}{2}(2-\beta )\bar{\rho }{\bar{v}}_{z}\frac{\partial {\bar{c}}_{\rm{s}}^{2}}{\partial z}-{\bar{v}}_{z}{\bar{c}}_{\rm{s}}^{2}\frac{\partial \bar{\rho }}{\partial z}]dz\\ =\displaystyle {\int }_{-h}^{+h}[\bar{\rho }{\nu }_{r\varphi }^{ {\mathcal R} }{\bar{s}}_{r\varphi }^{2}+\frac{1}{4\pi }{\epsilon }_{jlm}{\nu }_{lk}^{ {\mathcal R} }\frac{\partial {\bar{B}}_{m}}{\partial {x}_{k}}{\epsilon }_{jlm}\frac{\partial {\bar{B}}_{m}}{\partial {x}_{l}}]\\ dz-2{F}^{r+}. \% \end{array}\end{eqnarray} \tag{ 19 }$$

  The left hand side of the above equation is the signature of advection of energy flux in radial and vertical directions due to accretion and outflow respectively. The first and second terms on the right hand side express turbulent viscous heating due to the rφ component of the stress tensor and the turbulent Ohmic dissipation or Joule's heat loss respectively in the accretion flow. Turbulent Ohmic dissipation symbolizes a resistive flow due to which constant annihilation of the magnetic flux occurs. The last term of the equation is the flux of radiation that is escaping from the accretion-outflow surface. Similar to that of the induction in Equation (13), turbulent diffusion in the energy equation also arises only through the Reynolds stress tensor.The net heat flux generated in the accretion flow is defined by

  $$\begin{eqnarray}\begin{array}{ll}{q}^{+}= & \displaystyle {\int }_{-h}^{+h}\left(\bar{\rho }{\nu }_{r\varphi }^{ {\mathcal R} }{\bar{s}}_{r\varphi }^{2}\right.\\ & +\left.\frac{1}{4\pi }{\epsilon }_{jlm}{\nu }_{lk}^{ {\mathcal R} }\frac{\partial {\bar{B}}_{m}}{\partial {x}_{k}}{\epsilon }_{jlm}\frac{\partial {\bar{B}}_{m}}{\partial {x}_{l}}\right)dz.\end{array}\end{eqnarray} \tag{ 20 }$$

  Defining f as a constant cooling factor scaled through a relation

  $$\begin{eqnarray}&&{q}^{+}-2{F}^{r+}=(1-f){q}^{+},\end{eqnarray} \tag{ 21 }$$

  the final form of the heat flux equation is obtained by coupling Equations (19), (20) and (21). A constant f implies that the radiative heat loss from the accretion surface has a linear proportionality with the net heat flux generated in the system. Although f should vary radially, however, within a small inner accretion region where the accretion and outflow are coupled, this presumption is an acceptable approximation. The parameter f reflects the extent to which the accretion flow is advective. f → 1 represents an accretion flow with efficient cooling or which is radiatively efficient. At the other extremity, f → 0 epitomizes a strongly advective or advection dominated system which is radiatively inefficient.Equations (7), (9), (10), (11), (12), (14), (15) and (19) exhibit the dynamical behavior of accretion-induced outflow in a resistive MHD paradigm for the mean flow. In the next section, we explore a technique to solve these eight coupled partial differential equations using Equations (5), (20) and (21), rigorously.

The equations describing a turbulent magneto-fluid in the last section are extremely complicated and it is beyond our scope to solve them numerically. The only procedure left in this context is to explore some kind of approximate analytical method or quasi-analytical method to solve them. As a first approximation, we follow a power law self-similar approach (see Narayan & Yi 1994) with a polynomial expansion to solve the equations for the mean flow, to obtain the class of solutions. For the present purpose we invoke a generalized n^th degree polynomial expansion for all the mean quantities, where the mean flow variables are functions of both the radial and vertical coordinate. We restrict ourselves to the Newtonian paradigm as power law self-similar solutions are valid only in the limit of Newtonian approximation. The generalized Newtonian potential at any (r, z) in the accretion-outflow region is written in the form of a power series

$$\begin{eqnarray}\begin{array}{ll}{\varphi }_{G}(r,z)= & -GM\left({r}^{-1}-\frac{1}{2}{r}^{-3}{z}^{2}\right.\\ & \left. +\frac{3}{8}{r}^{-5}{z}^{4}-\cdots \right),\end{array}\end{eqnarray} \tag{ 22 }$$

where M is the mass of the central object. The gravitational force in the radial and vertical directions is then written in the form of a polynomial expansion as shown below.

$$\begin{eqnarray}&&{F}_{{\rm{Gr}}}(r,z)=GM\displaystyle \sum _{n=0}^{\infty }\left(\begin{array}{c}-3/2\\ n\end{array}\right){r}^{-2-2n}{z}^{2n}.\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&{F}_{Gz}(r,z)=GM\displaystyle \sum _{n=0}^{\infty }\left(\begin{array}{c}-3/2\\ n\end{array}\right){r}^{-3-2n}{z}^{2n+1},\end{eqnarray} \tag{ 24 }$$

where the parentheses on the right hand side of Equations (23) and (24) denote the usual binomial coefficient. F_Gr and F_Gz have even and odd symmetry in the z direction, respectively. We seek a polynomial expansion with a similar form for all mean flow variables, where the flow of matter in the accretion region is being considered to have reflection symmetry about the equatorial plane. All the hydrodynamical variables will have even symmetry except those directly related to F_Gz. The mean outflow velocity ${\bar{v}}_{z}$ will then have an odd symmetry in the z direction. Consequently, the radial component of the magnetic field ${\bar{B}}_{r}$ and the vertical component of the magnetic field B_z will have even and odd symmetry configurations, respectively, which is required from the divergence condition of the magnetic field. The odd symmetry configuration of the magnetic field has been used previously on other occasions in the context of outflows/jets from the accretion flows (e.g., Lovelace et al. 1987; Samadi & Abbassi 2016). The flow velocities, angular momentum and density are then written in the following polynomial form

$$\begin{eqnarray}&&\begin{array}{l}{\bar{v}}_{r}(r,z)=\displaystyle \sum _{n=0}^{\infty }{v}_{r2n}{r}^{a-2n}{z}^{2n},\\ {\bar{v}}_{\varphi }(r,z)=\displaystyle \sum _{n=0}^{\infty }{v}_{\varphi 2n}{r}^{b-2n}{z}^{2n},\\ \bar{\lambda }(r,z)=\displaystyle \sum _{n=0}^{\infty }{v}_{\varphi 2n}{r}^{b-2n+1}{z}^{2n},\\ {\bar{v}}_{z}(r,z)=\displaystyle \sum _{n=0}^{\infty }{v}_{z(2n+1)}{r}^{c-2n}{z}^{2n+1},\\ {\bar{c}}_{s}(r,z)=\displaystyle \sum _{n=0}^{\infty }{c}_{s2n}{r}^{d-2n}{z}^{2n},\\ \bar{\rho }(r,z)=\displaystyle \sum _{n=0}^{\infty }{\rho }_{2n}{r}^{e-2n}{z}^{2n}.\end{array}\end{eqnarray} \tag{ 25 }$$

Similarly, the components of the magnetic field can be expanded as

$$\begin{eqnarray}&&\begin{array}{l}{\bar{B}}_{r}(r,z)=\displaystyle \sum _{n=0}^{\infty }{B}_{r2n}{r}^{i-2n}{z}^{2n},\\ {\bar{B}}_{\varphi }(r,z)=\displaystyle \sum _{n=0}^{\infty }{B}_{\varphi 2n}{r}^{j-2n}{z}^{2n},\\ {\bar{B}}_{z}(r,z)=\displaystyle \sum _{n=0}^{\infty }{B}_{z(2n+1)}{r}^{k-2n}{z}^{2n+1},\end{array}\end{eqnarray} \tag{ 26 }$$

where v_r2n, v_φ2n, v_{z(2n + 1)}, C_s2n, ρ_2n, B_r2n, B_φ2n and B_{z(2n + 1)} are dimensionless coefficients which will be evaluated from MHD conservation equations.

We determine the exponents a, b, c, d, e, i, j and k by self comparison of various terms in the model equations. Substituting the solutions from Equations (25) and (26) in the MHD conservation equations and comparing the exponents of r and z, we obtain a = −1/2, b = −1/2, c = −3/2, d = −1/2, e = −3/2, i = −5/4, j = −5/4 and k = −9/4.

Using the above value of exponents and using the polynomials in Equations (25) and (26), Equations (9), (10), (11), (12), (14), (15) and (19) can be expanded in the power of aspect ratio (h/r). Using a simple technique given by the form

$$\begin{eqnarray}\begin{array}{ll} & \displaystyle \sum _{n=0}^{\infty }{u}_{n}{x}^{n}\displaystyle \sum _{n=0}^{\infty }{v}_{n}{x}^{n}\displaystyle \sum _{n=0}^{\infty }{w}_{n}{x}^{n}\\ = & \displaystyle \sum _{n=0}^{\infty }{x}^{n}\displaystyle \sum _{m=0}^{\infty }\displaystyle \sum _{l=0}^{\infty }{u}_{n-m}{v}_{m-l}{w}_{l},\end{array}\end{eqnarray} \tag{ 27 }$$

the above polynomial equations then can be written in a generic form as

$$\begin{eqnarray}&&{A}_{0}{\left(\frac{h}{r}\right)}^{0}+{A}_{1}{\left(\frac{h}{r}\right)}^{2}+{A}_{2}{\left(\frac{h}{r}\right)}^{4}+\cdots =0,\end{eqnarray} \tag{ 28 }$$

where, A₀, A₁, A₂ and . . . are zeroth order, first order, second order and higher order coefficients which are nonlinear functions of v_r2n, v_φ2n, v_{z(2n + 1)}, C_s2n, φ_2n, B_r2n, B_φ2n and B_{z(2n + 1)} corresponding to n = 0, n = 1, n = 2, . . ., respectively. Equation (28) is a linear combination in powers of h/r which are linearly independent.

If we neglect all the terms of order ≥ (h/r)², and only keep the zeroth order term, the generic Equation (28) after neglecting the magnetic field contribution will reduce to algebraic equations of Narayan & Yi (1994). To exemplify, in the appendix, we have shown the polynomial expansion of the integro-differential mass transfer expression, Equation (7), and the radial momentum balance expression, Equation (9), explicitly. Even if the accretion flow has considerable thickness, h in general would always be less than r. Further, magnetic stresses will compress or squeeze the accretion region. Considering the expression in Equation (28) up to the term (h/r)² would then be a reasonable approximation. Restricting the expansion up to (h/r)² corresponding to n = 1 and neglecting the terms with orders ≥ (h/r)⁴ in the generic expression (28), we equate A₀ and A₁ to zero, respectively. Extending this to MHD conservation equations and assuming isotropic turbulence, after rigorous algebra, we will then have fifteen independent nonlinear algebraic equations with sixteen unknown coefficients consisting of zeroth and first order only, which are shown below.

Equation (7) renders

$$\begin{eqnarray}&&\begin{array}{l}{\rho }_{0}({v}_{r0}-{v}_{z1})t+\left[{\rho }_{0}({v}_{r2}-{v}_{z3})+{\rho }_{2}({v}_{r0}-{v}_{z1})\right]\\ \times \frac{{t}^{3}}{3}=-\frac{\dot{M}}{4\pi }.\end{array}\end{eqnarray} \tag{ 29 }$$

Equation (12) yields

$$\begin{eqnarray}&&{B}_{z1}=\frac{{B}_{r0}}{4},\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&{B}_{z3}=\frac{3}{4}{B}_{r2}.\end{eqnarray} \tag{ 31 }$$

Equation (9) produces

$$\begin{eqnarray}&&-\frac{1}{2}{v}_{r0}^{2}-{v}_{\varphi 0}^{2}+GM-\frac{5}{2}{c}_{s0}^{2}-\frac{1}{16\pi }\frac{{B}_{\varphi 0}^{2}}{{\rho }_{0}}=0,\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&\begin{array}{l}-\frac{1}{2}{\rho }_{2}{v}_{r0}^{2}-3{\rho }_{0}{v}_{r0}{v}_{r2}-{\rho }_{2}{v}_{\varphi 0}^{2}-2{\rho }_{0}{v}_{\varphi 0}{v}_{\varphi 2}\\ +2{\rho }_{0}{v}_{z1}{v}_{r2}+GM({\rho }_{2}-\frac{3}{2}{\rho }_{0})-\frac{9}{2}({\rho }_{2}{c}_{s0}^{2}+2{\rho }_{0}{c}_{s0}{c}_{s2})\\ +\frac{1}{4\pi }(-\frac{5}{2}{B}_{\varphi 0}{B}_{\varphi 2}-\frac{9}{4}{B}_{z1}^{2}-2{B}_{z1}{B}_{r2})=0.\end{array}\end{eqnarray} \tag{ 33 }$$

Equation (10) provides

$$\begin{eqnarray}&&\frac{1}{2}{v}_{r0}{v}_{\varphi 0}+\frac{3}{4}{\alpha }_{r\varphi }t{v}_{\varphi 0}{c}_{s0}+\frac{1}{16\pi }\frac{{B}_{r0}{B}_{\varphi 0}}{{\rho }_{0}}=0,\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&\begin{array}{l}\frac{1}{2}({\rho }_{2}{v}_{r0}{v}_{\varphi 0}+{\rho }_{0}{v}_{\varphi 0}{v}_{r2})-\frac{3}{2}{\rho }_{0}{v}_{r0}{v}_{\varphi 2}\\ +2{\rho }_{0}{v}_{z1}{v}_{\varphi 2}+\frac{1}{2}{\alpha }_{r\varphi }t\left(\frac{3}{2}{\rho }_{2}{c}_{s0}{v}_{\varphi 0}\right.\\ \left. +\frac{3}{2}{\rho }_{0}{c}_{s2}{v}_{\varphi 0}+\frac{7}{2}{\rho }_{0}{c}_{s0}{v}_{\varphi 2}\right)\\ +\frac{1}{16\pi }({B}_{\varphi 0}{B}_{r2}+7{B}_{r0}{B}_{\varphi 0})=0. \end{array}\end{eqnarray} \tag{ 35 }$$

Equation (11) gives

$$\begin{eqnarray}&&\begin{array}{l}-\frac{3}{2}{\rho }_{2}{v}_{r0}{v}_{z1}+{\rho }_{0}{v}_{z1}^{2}-{\rho }_{0}GM\\ +2({\rho }_{2}{c}_{s0}^{2}+4{\rho }_{2}{c}_{01}{c}_{s2}\\ +\frac{1}{4\pi }(2{B}_{r0}{B}_{r2}+2{B}_{\varphi 0}{B}_{\varphi 2}\\ +\frac{9}{4}{B}_{r0}{B}_{z1})=0,\end{array}\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&\begin{array}{l}-\frac{3}{2}{\rho }_{2}{v}_{r0}{v}_{z1}-\frac{3}{2}{\rho }_{0}{v}_{r2}{v}_{z1}-\frac{7}{2}{\rho }_{0}{v}_{r0}{v}_{z3}\\ +{\rho }_{2}{v}_{z1}^{2}+4{\rho }_{0}{v}_{z1}{v}_{z3}\\ +GM({\rho }_{2}-\frac{3}{2}{\rho }_{0})+4({\rho }_{2}{c}_{s0}{c}_{s2}+{\rho }_{0}{c}_{s2}^{2})\\ +\frac{1}{8\pi }(4{B}_{r2}^{2}+4{B}_{\varphi 2}^{2})-\frac{1}{4\pi }\\ \left(-\frac{9}{4}{B}_{r2}{B}_{z1}-\frac{17}{4}{B}_{r0}{B}_{z3}\right)=0.\end{array}\end{eqnarray} \tag{ 37 }$$

Equation (14) generates

$$\begin{eqnarray}&&({v}_{r0}{B}_{z1}-{v}_{z1}{B}_{r0})+{\alpha }_{r\varphi }^{ {\mathscr R} }t\\ &&\left(2{c}_{s0}{B}_{r2}+\frac{9}{4}{c}_{s0}{B}_{z1}\right)=0.\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}&&({v}_{r2}{B}_{z1}+{v}_{r0}{B}_{z3})-({v}_{z3}{B}_{r0}+{v}_{z1}{B}_{r2})+{\alpha }_{r\varphi }^{ {\mathscr R} }t\\ &&\left(2{c}_{s2}{B}_{r2}+\frac{9}{4}{c}_{s2}{B}_{z1}+\frac{17}{4}{c}_{s0}{B}_{z3}\right)=0.\end{eqnarray} \tag{ 39 }$$

Equation (15) yields

$$\begin{eqnarray}&&\begin{array}{l}({v}_{\varphi 0}{B}_{z1}-{v}_{z1}{B}_{\varphi 0})+\frac{7}{4}({v}_{r0}{B}_{\varphi 0}-{v}_{\varphi 0}{B}_{r0})\\ +{\alpha }_{r\varphi }^{ {\mathscr R} }t\left(2{c}_{s0}{B}_{\varphi 2}+\frac{7}{16}{c}_{s0}{B}_{\varphi 0}\right)=0,\end{array}\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}&&\begin{array}{l}3({v}_{\varphi 0}{B}_{z3}-{v}_{\varphi 2}{B}_{z1})-3({v}_{z3}{B}_{\varphi 0}-{v}_{z1}{B}_{\varphi 2})+\\ \frac{7}{4}[({v}_{r0}{B}_{\varphi 2}+{v}_{r2}{B}_{\varphi 0})-({v}_{\varphi 2}{B}_{r0}+{v}_{\varphi 0}{B}_{r2})]\\ +{\alpha }_{r\varphi }^{ {\mathscr R} }t[6{c}_{s2}{B}_{\varphi 2}+\frac{15}{16}({c}_{s2}{B}_{\varphi 0}+9{c}_{s0}{B}_{\varphi 2})] \\ \,\,\,\,\,=0.\end{array}\end{eqnarray} \tag{ 41 }$$

Equation (19) using Equations (20) and (21) renders

$$\begin{eqnarray}&&\begin{array}{l}\frac{3}{2}(\beta -1){v}_{r0}{c}_{s0}=(1-f){\alpha }_{r\varphi }^{ {\mathscr R} }\\ t\left(\frac{9}{4}{v}_{\varphi 0}^{2}+\frac{1}{16}\frac{{B}_{\varphi 0}^{2}}{4\pi {\rho }_{0}}\right),\end{array}\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&\begin{array}{l}\frac{3}{2}(\beta -1){\rho }_{0}{v}_{r2}{c}_{s0}^{2}+(9\beta -15){\rho }_{0}{v}_{r0}{c}_{s0}{c}_{s2}\\ +\frac{1+3\beta }{2}{\rho }_{2}{v}_{r0}{c}_{s0}^{2}+6(2-\beta ){\rho }_{0}{v}_{z1}{c}_{s0}{c}_{s2}\\ -2{\rho }_{2}{v}_{z1}{c}_{s0}^{2}=(1-f){\alpha }_{r\varphi }^{ {\mathcal R} }t\left[\left(\frac{9}{4}{\rho }_{2}{c}_{s0}{v}_{\varphi 0}^{2}\right.\right.\\ \left.+\frac{9}{4}{\rho }_{0}{c}_{s2}{v}_{\varphi 0}^{2}+\frac{21}{2}{\rho }_{0}{c}_{s0}{v}_{\varphi 0}{v}_{\varphi 2}\right)\\ +\frac{1}{4\pi }\left(4{c}_{s0}{B}_{\varphi 2}^{2}+4{c}_{s0}{B}_{r2}^{2}+\frac{81}{16}{c}_{s0}{B}_{z1}^{2}\right.\\ +9{c}_{s0}{B}_{r2}{B}_{z1}+\frac{1}{16}{c}_{s2}{B}_{\varphi 0}^{2}\\ +\left.\left.\frac{9}{8}{c}_{s0}{B}_{\varphi 0}{B}_{\varphi 2}\right)\right].\end{array}\end{eqnarray} \tag{ 43 }$$

In the above equations t = h(r)/r. Using Equations (29), (30), (31), (33), (34), (35), (36), (38), (39), (40) and (42), and after very complicated and tedious algebra, we systematically determine the values of all the first order coefficients in terms of zeroth order coefficients of hydrodynamic variables. Substituting them in Equations (32), (37), (41) and (43), we successfully reduce 15 equations to four nonlinear algebraic equations comprised of zeroth order coefficients v_z1, v_φ, B_r0 and B_φ0 only. In an accretion flow, when there is no net vertical flux, the vertical height of the accretion geometry is calculated from hydrostatic pressure balance, assuming the pressure and density at the outer accretion surface to be zero.

However, this physical condition ceases to exist when the outflow is incorporated in the system, and the scale-height of the coupled accretion-induced outflow region becomes difficult to ascertain. In these circumstances, an accretion-outflow coupled surface can be treated as a photospheric height, delineating between the accretion-outflow surface and the transition region, leading to the outflow decoupling from the inflow. As a first approximation, we treat t = h(r)/r as a parameter to get physically plausible solutions. In reality, however, the physical conditions to launch outflow would self-consistently determine the vertical height of the inflow-outflow surface, from where the outflow decouples from the accretion region.

One can then solve the stated four nonlinear equations through an iterative Newton-Raphson method with an appropriate initial guess. In this way we can obtain the values of zeroth and first order coefficients of the corresponding dynamical variables. The values of the coefficients are scaled by putting G = M = 1, and $\dot{M}$ in units of ${\dot{M}}_{{\rm{Edd}}}$.

In the next section, we evaluate them for appropriate choices of parameters $\dot{M}$, β and f.

In Section 1 we have analyzed the necessity to have a sub-Keplerian advective accretion regime to eject strong outflows and jets from the accretion region in the vicinity of BHs. The corresponding accretion regime can be possibly envisaged if the mass accretion rate or net mass flow rate is considerably sub-Eddington ($\dot{M}\,\lesssim {10}^{-2}\,{\dot{M}}_{{\rm{Edd}}}$). The ratio of gas to the total pressure β and the radiative cooling factor f directly depends on $\dot{M}$. Looking meticulously into the relevant Equations (29–43), we notice that the continuity equation is written in the integral form through Equation (29), unlike the other hydrodynamical equations. This is being deliberately done to preserve information on $\dot{M}$ in the flow, as $\dot{M}$ is the most fundamental parameter which determines the nature of the BH accretion paradigm. However, this constrains the number of dynamical equations. For our case, we have 16 unknown coefficients but 15 equations. To resolve this, we proceed in the following way.

$\dot{M}$ carries the information of density of the flow which means that if $\dot{M}$ is known, in principle, density also is known. Equation (29) reveals that density in the accretion-outflow region is a function of two unknown coefficients ρ₀ and ρ₂. ρ₀ is the signature of inflow whereas ρ₂ is that of outflow. If the outflow is discarded, ρ₂ loses its significance and $\dot{M}$ then becomes the usual mass accretion rate from where ρ₀ can be easily calculated. If r_j represents the outer radial boundary of the accretion-outflow coupled region, then at r ≥ r_j the net mass flow $\dot{M}$ is equivalent to the mass accretion rate of the flow. At r ≥ r_j, density of the accretion flow is just a function of p₀, which is then given by (also see Narayan & Yi 1994)

$$\begin{eqnarray}&&{\rho }_{0}=\frac{{\left[5+2\frac{1-\beta }{1-f}\frac{{\alpha }_{r\varphi }}{{\alpha }_{r\varphi }^{ {\mathscr R} }}\right]}^{3/2}}{12\pi \sqrt{2}}\frac{\dot{M}}{{\alpha }_{r\varphi }}.\end{eqnarray} \tag{ 44 }$$

${\alpha }_{r\varphi }^{ {\mathscr R} }$ arises due to a statistical averaging approach. Note that the value of ρ₀ computed in Equation (44) is not exactly the same as that of ρ₀ in Equation (29), where it is coupled to ρ₂. However, we presume that the density in the accretion flow does not change abruptly due to emanation of the outflow and jet, as only a small fraction of matter is ejected through the outflow. With this presumption, we calculate ρ₀ in Equation (29) using Equation (44), and supply its value to the other nonlinear equations for further computation.

One can notice that the parameter ${\alpha }_{r\varphi }^{ {\mathscr R} }$ is related to α_rφ. It would be quite convenient if one can estimate the value of ${\alpha }_{r\varphi }^{ {\mathscr R} }$ in relation to α_rφ, at least approximately. Pessah et al. (2006b) showed that during the late time of the exponential growth of instability, the ratio of Maxwell stress to Reynolds stress becomes

$$\begin{eqnarray}&&\frac{\bar{{t}_{r\varphi }^{ {\mathscr M} }}}{\bar{{t}_{r\varphi }^{ {\mathscr R} }}}=\frac{4-q}{q},\end{eqnarray} \tag{ 45 }$$

where 0 < q < 2. q is related to angular velocity through the relation Ω ∼ r^−q. q = 3/2 signifies a Keplerian flow whereas q > 3/2 implies a sub-Keplerian or an advective accretion flow. Although this relationship corresponds to late times of the linear phase of instability, even in the saturated state of instability with fully developed MHD turbulence the qualitative nature of this relationship remains similar; for 1.6 ≲ q ≲ 1.9, the magnitudes of the ratios of Maxwell to Reynolds stresses in the turbulent saturated state nearly approach the corresponding magnitudes obtained during the late times of the linear phase of instability (Pessah et al. 2006b). For the case of a strongly advective flow we are interested in, it would be reasonable to use the relation given in Equation (45) to get an approximate estimate of the value of ${\alpha }_{r\varphi }^{ {\mathscr R} }$ in relation to α_rφ, as

$$\begin{eqnarray}&&\frac{{\alpha }_{r\varphi }^{ {\mathscr R} }}{{\alpha }_{r\varphi }}=\frac{q}{4}.\end{eqnarray} \tag{ 46 }$$

The qualitative feature of the dynamical coefficients for an accretion-induced outflow which are physically plausible should satisfy the following properties of the flow variables: ${\bar{v}}_{r}\to$ negative, ${\bar{v}}_{z}\to$ positive and ${\bar{B}}_{\varphi 0}\to$ negative. All the other variables should have a positive value. The positivity of ${\bar{B}}_{r}$ and ${\bar{B}}_{z}$ is related to open magnetic field lines threading the accretion flow, across which the accreting matter diffuses, and then gets accelerated outwards along the poloidal field lines by extracting angular momentum from the flow. The term $-{r}^{2}{({\bar{B}}_{\varphi }{\bar{B}}_{z})}_{h}$ [last term in Eq. (10)], which is a magnetic torque on the accretion flow, is attributed to the transport of angular momentum from the h(r) surfaces of the accretion-outflow region outwards. This term should be positive in order to launch a jet, which removes angular momentum from the flow and decouples from the accretion region. This premise makes an obvious choice for ${\bar{B}}_{\varphi }$ to have a negative value. It should be remembered that the negative values of many quantities do not necessarily mean that their magnitude is negative, but it represents the direction of their flow.

Next we compute the value of the coefficients of all dynamical variables for three relevant choices of $\dot{M}$ with appropriate values of β and f, conducive to forming outflows and jets.

4.1. Case 1. For $\dot{M}={10}^{-4}\,{\dot{M}}_{{\rm{Edd}}}$

The particular choice of $\dot{M}$ corresponds to RIAF, which is linked observationally to the low-hard state of BHXRBs and LERGs/LLAGNs. This type of flow is significantly gas pressure dominated and strongly advective. The flow is considerably geometrically thick, optically thin and radiatively inefficient. We choose appropriate values of β and f corresponding to this $\dot{M}$ to get physically plausible solutions. As this flow is highly sub- Keplerian and strongly advective, an appropriate choice of value of q ∼ 1.9 has been taken. In Tables 1, 2 and 3, we present the computed values of dynamical coefficients for β = 0.95 and f = 0.1 with suitable choices of α and t. For β ∼ 0.95, Γ is ∼ 1.635. It is interestingly found that for q ≪ 1.9 and in the range of 1.9 ≪ q < 2, we never identified any physically valid solutions. However, for the much stronger advective paradigm ($\dot{M}\ll \% \lt \lt {10}^{-4}{\dot{M}}_{{\rm{Edd}}}$), one would still obtain physically correct solutions for 1.9 ≪ q < 2.

Table 1.  $\dot{M}={10}^{-4}\,{\dot{M}}_{{\rm{Edd}}}$, β = 0.95, f = 0.1, α = 0.3, t = 0.1, q = 1.9

ρ0	vr0	vφ0	vz1	cs0	Br0	Bφ0	Bz1
7.4864e–5	−0.0899	0.3692	0.7069	0.5867	0.0036	−0.0050	1.3e–3
ρ2	vr2	vφ2	vz3	cs2	Br2	Bφ2	Bz3
−3.4889e–4	−0.2421	0.3471	−33.1161	−0.3551	0.0523	−0.2176	0.0793

Table 2.  $\dot{M}={10}^{-4}\,{\dot{M}}_{{\rm{Edd}}}$, β = 0.95, f = 0.1, α = 0.5, t = 0.1, q = 1.9

ρ0	vr0	vφ0	vz1	cs0	Br0	Bφ0	Bz1
4.4918e–5	−0.1494	0.3694	0.3859	0.5891	0.0051	−0.0041	4.25e–3
ρ2	vr2	vφ2	vz3	cs2	Br2	Bφ2	Bz3
−1.1174e–4	−0.0925	0.3211	−7.0486	−0.2300	0.0257	−0.0603	0.0193

Table 3.  $\dot{M}={10}^{-4}\,{\dot{M}}_{{\rm{Edd}}}$, β = 0.95, f = 0.1, α = 0.5, t = 0.2, q = 1.9

ρ0	vr0	vφ0	vz1	cs0	Br0	Bφ0	Bz1
4.4918e–5	−0.2876	0.3636	0.6211	0.5977	0.0149	−0.0061	5.2e–3
ρ2	vr2	vφ2	vz3	cs2	Br2	Bφ2	Bz3
−2.1111e–4	−0.7225	0.6940	−9.7460	−0.3859	0.0598	−0.0724	0.0448

We will analyze the family of analytical solutions later, however, we note that the values of the coefficients are in conformity with desired physically valid solutions (mentioned earlier). We notice that they are obtained only at high a (≳0.3) and at a reduced vertical scale-height (t≳0.1) of the accretion-outflow coupled region. Although t is small, it corresponds to a geometrically thick accretion flow (t≳0.1, t measures the degree of flow thickness). We have elucidated the necessity of a geometrically thick accretion flow to eject outflow and jet in previous paragraphs. However, corresponding to α ∼ 0.3, we never found any physically acceptable solution for t > 0.1. For higher α ∼ 0.5, however, we obtain a solution at a maximum value of t ∼ 0.2. This infers that with the increase of α, plausible physical solutions can be obtained with a thicker accretion geometry. Nevertheless, very high α (> 0.5) accretion flow might not be realistic in nature, and hence we restrict our study to a maximum plausible value of α ∼ 0.5.

Obtaining physically valid solutions of an accretion-induced outflow at a reduced geometrical thickness (scale-height) of the accretion region, as compared to a non-magnetized accretion flow (without outflow) like advection dominated accretion flows (ADAFs) for a similar accretion paradigm, is owing to the fact that the magnetic stresses in the flow have a tendency to squeeze or compress the region by acting oppositely to the thermal pressure gradient, consequently reducing the scale-height of the accretion region. This squeezing effect has been discussed by several other authors (e.g., Soria et al. 1997; Casse & Keppens 2004) in the context of magnetized accretion flow. This can be noticed from the vertical momentum balance equation (Eq. (11)). As we intend to see the effect of magnetic field on the geometrical thickness of the accretion region, for simplicity we ignore the outflow in Equation (11). Further, we found that with the increase in z, the magnitude of all components of the magnetic field increases, however the density and thermal pressure decrease with increase in z. As we ignore the outflow, for simplicity of our calculation, we consider density and thermal pressure at accretion flow scale-height (h) to be negligible as compared to theirequatorial values; this will not alter the qualitative natureof our argument. Expanding the terms in Equation (11), and restricting up to h²/r², Equation (11) will reduce tothe magnetohydrostatic equilibrium equation, given by

$$\begin{eqnarray*}\begin{array}{ll}{\bar{\rho }}_{{\rm{eq}}}\displaystyle \frac{{h}^{2}}{{r}^{3}} & \sim {\bar{P}}_{{\rm{eq}}}-\displaystyle \frac{1}{8\pi }{r}^{-5/4}\\ & \times ({B}_{\varphi 0}{B}_{\varphi 2}+{B}_{r0}{B}_{r2}+\displaystyle \frac{9}{4}{B}_{r0}{B}_{z1})\displaystyle \frac{{h}^{2}}{{r}^{2}},\end{array}\end{eqnarray*}$$

where ${\bar{\rho }}_{{\rm{eq}}}$ and ${\bar{P}}_{{\rm{eq}}}$ are the density and thermal pressure at equatorial plane respectively. The scale-height of the accretion flow would then approximately be given by

$$\begin{eqnarray}&&h\sim \sqrt{\displaystyle \frac{{c}_{{\rm{s}}}^{2}{r}^{3}}{1+\displaystyle \frac{1}{8\pi {\rho }_{0}}({B}_{\varphi 0}{B}_{\varphi 2}+{B}_{r0}{B}_{r2}+\displaystyle \frac{9}{4}{B}_{r0}{B}_{z1})}}.\end{eqnarray} \tag{ 47 }$$

If we neglect the magnetic components, Equation (47) is then the usual hydrostatic scale-height of the accretion flow. Due to the presence of the magnetic field, the scale-height is now approximately reduced by a factor $\sqrt{1+\frac{1}{8\pi {\rho }_{0}}({B}_{\varphi 0}{B}_{\varphi 2}+{B}_{r0}{B}_{r2}+\frac{9}{4}{B}_{r0}{B}_{z1})}$.

As the region becomes more compressed, the thermal content of gas increases and the excess thermal pressure gradient will help in lifting plasma vertically outwards.

4.2. Case 2. For $\dot{M}={10}^{-3}\ {\dot{M}}_{\rm{Edd}}$

4.3. Case 3. For $\dot{M}={10}^{-2}\ {\dot{M}}_{\rm{Edd}}$

In this section, we analyze the family of solutions for advective flows in the accretion-outflow coupled region with $\dot{M}\lesssim {10}^{-2}\ {\dot{M}}_{\rm{Edd}}$. Although the mean flow variables vary in both r and z, we do not display three dimensional figures as they are very obscure and difficult to interpret. As outflow and jet effuse out from the inner region of the accretion flow, we restrict our analysis up to 50 Schwarzschild radii (r_g) within which we presume that the accretion and outflow are coupled, where r_g = 2GM/c². Also it has been stated by Kumar & Chattopadhyay (2013) that VLBI observations of M87 (Junor et al. 1999) have shown that the jet originates from the vicinity (∼ 50r_g) of a BH/compact object. However, we do not expect any outflow of the accreting plasma in the close vicinity of a BH, as the accretion flow remains highly bounded in the close vicinity of a gravitationally starved BH. So, we restrict our study up to 5 r_g in the inward radial direction, a quite reasonable choice for the inner radius of the accretion-outflow coupled region. Both r and z coordinates in the figures are expressed in units of r_g. We express $\dot{M}$ in our entire analysis in units of ${\dot{M}}_{\rm{Edd}}$. The dynamical solutions are shown in the following figures with G = M = 1. All the flow variables in the figures represent mean quantities.

Figures 1 and 2 display the variation of vertically averaged flow variables as functions of radial coordinate r in different accretion paradigms. Variation of ${\bar{v}}_{z}$ and ${\bar{B}}_{z}$ are shown only along coupled accretion-outflow surface h(r) as they are odd functions in z. Figure 1 indicates that with the decrease in $\dot{M}$, the poloidal components of the velocity (${\bar{v}}_{r},{\bar{v}}_{z}$) and sound speed consistently increase. Conversely, the magnitude of ${\bar{v}}_{\varphi }$ increases with the flow becoming less advective and more centrifugally dominated. It is interestingly found from the tables in Section 4 that the value of the large-scale poloidal magnetic field enhances with an increase in the geometrical thickness of the accretion flow. As one moves from a strongly advective regime to a moderately advective regime, there is a sharp fall in the value of poloidal component of the magnetic field. This is owing to the fact that the geometrical thickness of the accretion region corresponding to moderately advective accretion flow is much less as compared to the case in the strongly advective accretion paradigm, as shown in Figure 2(a) and (b). This demonstrates the dominating influence of vertical thickness of the accretion flow structure on the poloidal component of the magnetic field. In contrast, the toroidal component of the magnetic field ${\bar{B}}_{\varphi }$ always increases with the flow becoming less advective and more centrifugally/rotationally dominated (Fig. 2(c)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 3 shows variation of the poloidal component of velocity and magnetic field with α for different accretion paradigms along the radial distance r. ${\bar{v}}_{r}$ and ${\bar{B}}_{r}$ are vertically averaged quantities, whereas ${\bar{v}}_{z}$ and ${\bar{B}}_{z}$ are along the coupled accretion-outflow surface h(r). We choose the value of corresponding t for different α to be maximum as illustrated in Section 4. This is because the system has a greater tendency to relax itself to the maximum possible height available to render a physically plausible solution of the coupled accretion-outflow, as the geometrically thicker accretion flow is more conducive to propel plasma vertically outwards from the accretion region. We find that with a small increase in α from 0.3 to 0.5, the value of the poloidal component of velocity and the magnetic field increases for both accretion paradigms. We do not show the variation of other flow variables with α as their dependence on α is insignificant for a particular $\dot{M}$, which can be verified from the tables in Section 4.

Zoom In Zoom Out Reset image size

In Figures 4 and 5 we present the variation of the few dynamical variables in z (which are very relevant to outflow) at any arbitrary location along r, for ${10}^{-4}\ {\dot{M}}_{\rm{Edd}}\lesssim \dot{M}\lesssim {10}^{-2}\ {\dot{M}}_{\rm{Edd}}$. Figure 4(a) shows the dependence of ${\bar{v}}_{z}$ with z for different $\dot{M}$ for $\alpha =0.3$, similar to that in Figure 1. It is found that initially ${\bar{v}}_{z}$ increases rapidly in z, however there is a sudden deceleration of ${\bar{v}}_{z}$ as the flow approaches the coupled accretion-induced outflow surface. This is due to the fact that the inward vertical component of gravitational force (F_Gz) dominates near the coupled accretion-outflow surface. Also at low $\dot{M}$ corresponding to moderately advective accretion flow, the increase of ${\bar{v}}_{z}$ in the vertical direction is much steeper. The truncation of the curves at a particular z represents the corresponding vertical thickness t. Figure 4(b) shows the variation of ${\bar{v}}_{z}$ in z with α for two different accretion paradigms similar to that in Figure 3. The nature of the curves is similar to that in Figure 4(a). Nevertheless, with the increase in α, the value of ${\bar{v}}_{z}$ gets enhanced, and there is a steeper increase of ${\bar{v}}_{z}$ in z. The nature of the variation in the toroidal component of the magnetic field along z for different $\dot{M}$ is shown in Figure 4(c). It is seen that with the increase in $\dot{M}$, ${\bar{B}}_{\varphi }$ increases at a much faster rate along z. Variation of ${\bar{B}}_{\varphi }$ with α is insignificant (see tables in Sect. 4) and hence is not graphically displayed. In Figure 4(d) and 4(e), we show the variation of density in z corresponding to strongly advective and moderately advective accretion regimes, respectively. Both these figures indicate that $\bar{\rho }$ decreases with the increase in α for all $\dot{M}$. Also, there is a steeper fall of density in z with the increase in $\dot{M}$. Figure 5 depicts the variation of toroidal velocity ${\bar{v}}_{\varphi }$ in z. With the increase in $\dot{M}(\dot{M}\gt {10}^{-3}\ {\dot{M}}_{\rm{Edd}})$ as the flow becomes less advective and more centrifugally dominated, ${\bar{v}}_{\varphi }$ decreases in z for all relevant values of α and f. This feature of ${\bar{v}}_{\varphi }$ has already been remarked on in Section 4. In Figure 5(c) and 5(d), we show them for $\dot{M}\sim {10}^{-2}\ {\dot{M}}_{\rm{Edd}}$. We also find that with an increase in α the profile of ${\bar{v}}_{\varphi }$ in z becomes steeper.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

If the system has a large residual toroidal velocity (centrifugally dominated), it is possible that the angular momentum loss in the vertical direction due to the magnetic torque will be considerably high, owing to which the degree of angular momentum loss proportionately increases as the matter flows vertically outwards.

In Figure 6(b) and 6(a) we depict the profile of differential magnetic torque ($-{r}^{2}{\bar{B}}_{\varphi }{\bar{B}}_{z}$) along z, as well as differential magnetic torque acting on the coupled accretion-outflow surface $(-{r}^{2}{\bar{B}}_{\varphi h}{\bar{B}}_{zh})$ along r, respectively, corresponding to different accretion paradigms. This term represents the magnetic extraction of angular momentum by the outflowing matter, which is called magnetic braking. The curves show that with an increase in the value of $\dot{M}$ as well as with an increase in α (for a specific value of $\dot{M}$), the value of differential magnetic torque increases. As $\dot{M}$ increases, with the accretion flow becoming less advective and more centrifugally dominated, the extraction of angular momentum by the outflowing plasma is greatly enhanced. Consequently, the gas gets centrifugally accelerated and would leave the accretion region by removing angular momentum from the accreting matter. This mechanism predominantly determines the outward flow of matter for accretion flow with moderate advection. If the system has a large residual toroidal velocity (centrifugally dominated), it is possible that the angular momentum loss in the vertical direction due to the magnetic torque will be considerably high, owing to which the degree of angular momentum loss proportionately increases as the matter flows vertically outwards. If the system has a large residual toroidal velocity (centrifugally dominated), it is possible that angular momentum %loss in the vertical direction due to magnetic torque will be exceptionally high, owing to which the degree of angular momentum loss proportionately increases as the matter flow vertically outwards. As a consequence, there will be an eventual decrease of toroidal velocity ${\bar{v}}_{\varphi }$ in the vertical direction within the accretion-outflow coupled region, as seen in Figure 5(c) and (d), corresponding to accretion flow with $\dot{M}\sim{10}^{-2}\,{\dot{M}}_{{\rm{Edd}}}$. On the other hand, if the accretion flow is predominantly gas pressure dominated as in a strongly advective regime (with $\dot{M}\lesssim {10}^{-3}\,{\dot{M}}_{{\rm{Edd}}}$), the gas pressure gradient would play a more contributory role to lift the plasma vertically outwards with the help of magnetic forces, and the effective contribution of magnetocentrifugal acceleration to control the dynamics of outflowing matter gets curtailed as compared to that in a more centrifugally dominated accreting system.

Zoom In Zoom Out Reset image size

Figure 6(d) and 6(c) shows the variation of poloidal component of the magnetic field ${\bar{B}}_{{\rm{P}}}[=\sqrt{}({\bar{B}}_{r}^{2}+{\bar{B}}_{z}^{2})$] along z and at coupled accretion-outflow surface $({\bar{B}}_{{\rm{Ph}}})$ along r, respectively. The nature of the curves indicates that as high α renders the accretion-induced outflow to a greater geometrical thickness (commented earlier), the large-scale poloidal field gets strongly augmented with an increase in α for a specific $\dot{M}$ due to the dominating influence of vertical thickness on ${\bar{B}}_{{\rm{P}}}$. Thus with the increase in turbulent viscosity parameter α, the value of differential magnetic torque responsible for centrifugal acceleration of the outflowing plasma, as well as the large-scale poloidal magnetic field ${\bar{B}}_{{\rm{P}}}$, gets strongly augmented; consequently enhancing the transport of vertical flux outwards. In the moderately advective accretion paradigm with more centrifugal domination, the effective contribution to launch and eject matter vertically outwards from the accretion region arises mainly from the magnetocentrifugal acceleration. A small increase in the turbulent viscosity parameter α from 0.3 to 0.5 intensifies the process of extraction of the angular momentum due to magnetic torque. The eventual result is the enhanced transport of outward vertical flux with an increase in effective angular momentum transport in the z direction. This causes ${v}_{\varphi 2}$ to become negative at a lower f and higher β as compared to that for α = 0.3, as stated in the last two lines of Section 4.3.

In Figure 7(a) we show the variation of the ratio of ${\bar{v}}_{z}$ and poloidal Alfvén velocity ${\bar{v}}_{{\rm{AP}}}[={\bar{B}}_{{\rm{P}}}/\sqrt{}(4\pi \bar{\rho })]$ in z for different $\dot{M}$ and α. With an increase in $\dot{M}$ as the system becomes less advective and more centrifugally/rotationally dominated, there is an increase in the ratio of ${\bar{v}}_{z}/{\bar{v}}_{{\rm{AP}}}$. In contrast, with an increase in α for a specific $\dot{M}$ there is a sharp fall in the value of the above ratio. Figure 7(b) shows the variation of the ratio of ${\bar{v}}_{z}$ and net Alfvén velocity ${\bar{v}}_{AT}[=\bar{B}/\sqrt{}(4\pi \bar{\rho })]$ in z similar to that in Figure 7(a). Figure 7(c) depicts the profile of plasma ${\beta }_{{\rm{P}}}[={\bar{B}}^{2}/(8\pi \bar{\rho }{\bar{c}}_{{\rm{s}}}^{2})]$ in z for different $\dot{M}$ and α. We find that ${\beta }_{{\rm{P}}}$ always increases steadily in the vertical direction, and its value becomes strongly augmented with a small increase in α, however, it always remains mostly below equipartition for all relevant $\dot{M}$ and α. Moreover it is found that with the increase in $\dot{M}$, in general, ${\beta }_{{\rm{P}}}$ decreases.

Zoom In Zoom Out Reset image size

In Figures 8 and 9 we make a comparison of our magnetized accretion-outflow solutions with those of ADAF type solutions in Narayan & Yi (1994). In Figures 8 and 9, we compare the radial profiles for radial velocity, orbital velocity, sound speed and density, obtained for our magnetized accretion-outflow solutions with those of ADAF type solutions, corresponding to different $\dot{M}$ and different α. We found that the magnitude of radial velocity along r that we obtained in our accretion-induced outflow is less as compared to that obtained in self-similar ADAF. On the other hand, the magnitude of orbital velocity along r obtained in our case is higher as compared to that obtained in ADAF. Also there is a marginal decrease in the magnitude of sound speed or equivalently the temperature of the gas along r obtained in our accretion-induced outflow as compared to that obtained in ADAF. However, the magnitude of density of the gas along r in our accretion-induced outflow is found to be quite similar to that obtained in the case of ADAF.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Observationally it is found that the low/hard state of BHXRBs, which are supposed to be powered by geometrically thick strongly advective sub-Eddington (presumably with $\dot{M}\lesssim {10}^{-3}\,{\dot{M}}_{{\rm{Edd}}}$) and consequently quasispherical and RIAFs, emanates strong outflows and relativistic jets. Outflows and jets are not observed in high/soft states of BHXRBs which are powered by a geometrically thin and optically thick standard Keplerian accretion disk. The physics of the origin and launching of outflows/jets in Galactic BH systems (also called microquasars) is supposed to be similar with that corresponding to SMBHs in AGNs, as AGNs may be seen as scaled up BHXRBs (McHardy et al. 2006; Körding et al. 2006). A geometrically thick advective accretion flow having a substantial amount of gas pressure, with strong advection, is more conducive to effuse and accelerate plasma in the vertical direction out of the inner accretion region. Although we do not aspire to explore the exact mechanism of launching and ejection of jets, however, it is generally conceived that the origin, launching and ejection of outflow and jet from the accretion flow is an MHD process. In the present work we mainly focus on accretion powered hydromagnetic outflows.

Although the distinctive cause of the origin and launching of accretion powered hydromagnetic outflows is still inconclusive, however, it is certain that the dynamics of the outflowing matter should be intrinsically coupled to the accretion dynamics through the fundamental laws of conservation of matter, momentum and energy within the coupled accretion-induced outflow region, and should not be treated as dissimilar objects. Conservation laws are the most valuable foundation in physics, and play a significant role in understanding astrophysical outflows and jets. This is because the physical dynamics of the coupled inflow and outflow are essentially governed by the laws of conservation. The nature of the dynamical solutions in the accretion-outflow coupled region should then reflect upon the physical conditions/criteria to eject outflows. For the theoretical analysis of the accretion-outflow coupling, one needs to be very thoughtful about proper modeling of the system, which essentially needs to solve a complete set of MHD conservation equations in a 2.5-dimensional viscous, resistive and advective paradigm. In Section 2 we have endeavored to describe a robust form of accretion-outflow coupled MHD set of equations in the viscous, resistive and advective paradigm, upholding all the conservation laws in the 2.5-dimensional mean field MHD regime. The dynamical flow variables are then represented through their appropriate mean values, describing a mean flow in the disk, involving a time average with time scale much larger than the time scales of turbulent fluctuations. The mean flow variables vary both in (r, z). Statistical averaging gives rise to the emergence of various turbulent correlation terms, where we restrict our study to first order turbulent correlation. Note that the turbulent magnetic diffusivity and turbulent viscous terms in the induction and energy conservation equations arise only from the kinetic part of the turbulent stress tensor through Reynolds stress.

In this work, we have assumed isotropic turbulence and also neglected the contribution of other turbulent stress tensors apart from rφ, which may be dynamically important. The contribution of the rφ component would be dynamically more dominant as it is responsible for the radial transport of angular momentum outwards. Vertical transport of angular momentum occurs mainly through large-scale magnetic stresses. Nonetheless, in the future we would like to examine the possibility of their inclusion, as well as investigate the nature of the flow with anisotropic turbulence.

The inflow and outflow are governed by eight coupled integro-partial differential MHD equations in the cylindrical geometry. Limited observational inputs put constraints on the boundary conditions as well as the scaling relation between accretion and outflow. Owing to the fact that it is beyond the scope to have complete global numerical solutions of the said coupled partial differential MHD equations, this motivated us to invoke a necessary and proper quasi-analytical method to solve them. Ever since the work of Narayan & Yi (1994), use of power law self-similarity in studying the accretion flow dynamics, especially ADAFs, to explain the nature of LLAGNs has become widely popular. A realistic strongly advective accretion flow preserves self-similarity reasonably well, within an appreciable region of the flow (Narayan & Yi 1994; Narayan et al. 1997), and has been widely used to explain observational features in LLAGNs (see Ho 2008 for a review). Previously, self-similar methods had been indeed used to study outflow from the accretion disk on many occasions (see the references in Sect. 1). It is being found from many studies that self-similarity holds approximately well in the context of outflows (Fendt 2006; Pudritz et al. 2007; McKinney & Narayan 2007b,a; Stute et al. 2014). Keeping the essence of power law self-similarity, we sought a generalized n^th degree polynomial expansion for all the mean flow variables in two dimensions, and solve the complete set of coupled integro-partial differential MHD conservation equations within the accretion-outflow coupled region self-consistently, in the 2.5-dimensional viscous, resistive advective paradigm, where we have restricted up to order [h(r)/r]². The flow of matter in the accretion region is being considered to have reflection symmetry about the equatorial plane. Consequently the vertical component of the magnetic field B_z has odd symmetry about the equatorial plane, which is required from the divergence condition of the magnetic field. The odd symmetry configuration of magnetic field has been used previously on other occasions in the context of outflows/jets from the accretion flows (e.g., Lovelace et al. 1987; Samadi & Abbassi 2016). It has been argued by Samadi & Abbassi (2016) that odd symmetry may perhaps be more realistic due to the fact that the magnetic field is due to dynamo processes, and in the disk, the fastest growing dynamo field mode has odd symmetry (Brandenburg & von Rekowski 2007). We comment on this aspect in more details in paragraphs 7 and 8 of this section. Nonetheless, in the future, we would like to pursue a similar kind of study with an even symmetry configuration of magnetic field.

In Sections 4 and 5, we have analyzed the nature and behavior of our MHD solutions in the advective accretion paradigm with $\mathop{M}\limits^{.}\precsim {10}^{-2}\,{\mathop{M}\limits^{.}}_{\mathrm{Edd}}$. Although we have not intended to explore the physical mechanism of outflow/jet launching in the present study, the quasi-stationary dynamical solutions of the accretion-induced outflow carry information about the physical conditions/criteria to propel matter vertically outwards out of the accretion-outflow region. We have mainly focused within the accretion-outflow coupled region where the flow is essentially bounded. We obtain solutions at a reduced vertical thickness irrespective of the nature of the accretion paradigm which we have focused on. The magnetic field tends to compress or squeeze the accretion region by counterbalancing the thermal pressure gradient. We found that the large-scale poloidal component of the magnetic field is enhanced with the increase in geometrical thickness of the accretion flow, consistently. With the increase in $\mathop{M}\limits^{.}$ as ${10}^{-4}\,{\mathop{M}\limits^{.}}_{\mathrm{Edd}}\ll \mathop{M}\limits^{.}\precsim {10}^{-2}\,{\mathop{M}\limits^{.}}_{\mathrm{Edd}}$, and with the accretion flow becoming less advective and more centrifugally dominated with lesser geometrical thickness, i.e., for a flow with moderate advection, there is a sharp fall in the value of the poloidal component of the magnetic field, however with a strong enhancement in the value of the toroidal component of the magnetic field and consequently the differential magnetic torque ($-{r}^{2}{\bar{B}}_{\varphi }{\bar{B}}_{z}$).

This term (differential magnetic torque) is responsible for the magnetic extraction of angular momentum to magnetocentrifugally accelerate the outflowing plasma out of the radial accretion flow, and this predominantly determines the outward flow of matter in a moderately advective accretion paradigm which is more centrifugally dominated. However, with the decrease in $\mathop{M}\limits^{.}$ as the flow becomes strongly advective ($\mathop{M}\limits^{.}\precsim {10}^{-3}\,{\mathop{M}\limits^{.}}_{\mathrm{Edd}}$) and geometrically more thick with strong gas pressure and inefficient cooling, despite a decrease in the value of $-{r}^{2}{\bar{B}}_{\varphi }{\bar{B}}_{z}$, a consistent increase in v_z occurs, indicating that the gas pressure gradient might play a more contributory role to lift the plasma vertically outwards with the help of magnetic forces. The plasma in the accretion flow can be lifted outwards and ejected, only if some physical process can overcome the effect of inward vertical force due to central gravity. The dynamical behavior of the solutions indicates that in the advective paradigm both magnetocentrifugal acceleration and thermal pressure gradient, along with magnetic forces, will help in lifting and accelerating the plasma vertically outwards, and the gas material will diffuse outwards across magnetic field lines. However, the effective contribution of either magnetocentrifugal acceleration or thermal pressure gradient to lift the plasma vertically outwards depends on the degree to which the flow is advective. In fact, with the increase in mass accretion rate as the flow tends to become less advective and more centrifugally dominated with lesser geometrical thickness, in general, the efficacy of the disk to eject outflows diminishes.

In Paper II (in preparation), we have quantitatively demonstrated this aspect with the increase in $\mathop{M}\limits^{.}$ from ${10}^{-4}\,{\mathop{M}\limits^{.}}_{\mathrm{Edd}}$ to ${10}^{-2}\,{\mathop{M}\limits^{.}}_{\mathrm{Edd}}$; the accretion flow with $\mathop{M}\limits^{.}\sim {10}^{-2}\,{\mathop{M}\limits^{.}}_{\mathrm{Edd}}$ is least conducive to ejecting outflows.

We obtain dynamical solutions in the accretion-outflow coupled region only at high turbulent diffusive parameter α(≳0.3). The accretion-induced outflow solutions have a profound dependence on turbulent diffusive parameter α. It is being interestingly found from our solutions that the enhancement in α renders the accretion-induced outflow region to attain a greater geometrical thickness. Consequently, the poloidal component of magnetic field ${\bar{B}}_{{\rm{P}}}$, as well as the differential magnetic torque ($-{r}^{2}{\bar{B}}_{\varphi }{\bar{B}}_{z}$), gets strongly augmented, enhancing the transport of vertical flux outwards. Also the plasma beta ${\beta }_{{\rm{P}}}$ increases steadily in the vertical direction, and its value gets strongly augmented with a small increase in α, however, it always remains mostly below equipartition within the accretion-outflow coupled region.

Although we expect the accretion flow to have a large α owing to the advective nature of the flow, however, the values of α that we have obtained in our solution may have been slightly overestimated. Nonetheless, it is indeed being found from previous works (Narayan & Yi 1994, 1995; Narayan et al. 1997; Yuan et al. 2008) that strongly advective accretion flows, in general, are mostly plausible for large values of α(α ≳ 0.1). Gu & Lu (2000) also showed that the transition from an outer geometrically thin Keplerian disk to an advection-dominated accretion flow is possible for α > 0.5. McKinney & Narayan (2007a), in their general relativistic MHD simulation of the disk-outflow model, found a large turbulent viscosity parameter in the accretion disk in the vicinity of the BH. Further, King et al. (2007) suggested a typical range of α ~ (0.1 − 0.4) from observational evidence.

One of the important approximations we have used in our study is to treat the scale-height of the accretion-outflow coupled region as a parameter. In reality the physical conditions to launch outflow would consistently determine the vertical height of the inflow-outflow surface, from where the outflow decouples from the accretion region. Moreover, in our study we have neglected the effect of spin of the BH in the accretion dynamics and its subsequent impact on the outflow, as the self-similar technique can only be used in the Newtonian approximation. This restricts us from using this method to investigate the physical behavior of the system in the extreme vicinity of the BH, where general relativistic effects are indispensable. Although power law self-similarity is an analytical approximation, and the quantitative feature of the solutions may have been either overestimated or undervalued, the dynamical solutions show consistent and predictable behavior, and exhibit many physical insights on the nature of the accretion-induced outflow, as well as reflect upon the relevant physical conditions that propel and eject plasma out of the accretion flow. Power law self-similarity, thus, seems to be a reasonably good approximation, within the accretion-outflow coupled region.

A very important determinant to launch outflows/jets is the existence of a large-scale magnetic field in the disk. One of the most likely mechanisms for the origin of this large-scale field is through the internally/locally generated self-sustaining MHD dynamo within the disk, possibly involving a kind of magneto rotational instability (MRI) or in short, through 'MRI supported disk dynamo,' from a seed magnetic field (Balbus & Hawley 1998; Lesur & Ogilvie 2008). Recently, Salvesen et al. (2016) performed a local simulation of a vertically stratified accretion disk supported by an MRI dynamo, in the ideal MHD paradigm, using 'zero net vertical magnetic flux' and with the incorporation of a near-realistic vertical boundary condition (outflowing) to capture the physically correct effect of magnetic buoyancy in the disk. The authors found that this local simulation with zero net vertical magnetic flux renders the disk to be unstable against magnetic buoyancy with a considerable expulsion of toroidal magnetic field from the disk, and the magnetic field expelled by the magnetic buoyancy cannot be effectively replenished by the 'MRI-dynamo' generated toroidal magnetic field. This makes the disk settle into a weakly magnetized state with the magnetic field becoming dynamically less important for launching strong outflows and jets. It then appears that a disk with zero net vertical magnetic flux will be less susceptible to launching outflows and jets. This result seems to be particularly sensitive in the context to of our present study, as in a purely axisymmetric case (which we considered here), the odd symmetry configuration of magnetic field that we used will render zero net vertical magnetic flux. Salvesen et al. (2016), however, pointed out that the local approximation they used in their study is more appropriate in the context of a geometrically thin accretion disk, where the time scale for the buoyant field to escape vertically is less compared to the time scale of viscous radial advection of toroidal field. In case of geometrically thick advective accretion flow (with ${10}^{-4}\,{\mathop{M}\limits^{.}}_{\mathrm{Edd}}\precsim \mathop{M}\limits^{.}\precsim {10}^{-2}\,{\mathop{M}\limits^{.}}_{\mathrm{Edd}}$) which we have focused upon in our present study, nonetheless, buoyant expulsion of magnetic field vertically would presumably take place on a relatively larger time scale. Many previous works on the effect of magnetic buoyancy on accretion disk exist in the literature (e.g., Torkelsson 1993; Chakrabarti & D'Silva 1994; D'Silva & Chakrabarti 1994; Rozyczka et al. 1996; Ziegler 2001; Johansen & Levin 2008), with a few of them showing plausible suppression of buoyant escape of magnetic flux from the disk (Torkelsson 1993; D'Silva & Chakrabarti 1994; Ziegler 2001). On the other hand, MRI-dynamo activity is found to be highly sensitive to non-ideal (resistive) MHD effects (Lesur & Ogilvie 2008; Riols et al. 2015). A global MHD simulation in the context of a vertically stratified geometrically thick advective flow, likely supported by MRI-dynamo in the resistive paradigm with background zero net vertical magnetic flux and with appropriate vertical boundary condition, thus, would provide a more definitive understanding of the effect of magnetic buoyancy on disk magnetization, and consequently would shed more light on the nature of magnetic field configuration in the disk, as well as on the magnetic field strength, required for ejection processes.

Notwithstanding, in reality, the advective accretion flow is not expected to be purely axisymmetric (Mukhopadhyay & Chatterjee 2015), and hence, even with odd symmetry configuration of magnetic field, the advective flow would then be expected to retain, at least some, net vertical magnetic flux. Although the present study has been pursued in a purely axisymmetric framework (owing to mathematical simplicity), slight non-axisymmetry may not alter the qualitative nature of the solutions. A detailed non-axisymmetric analysis of such a system, however, is beyond the scope of the present work. There may be, however, viable alternative routes which can give rise to a dynamically important magnetic field in the disk: either through other types of internally generated MHD dynamo activities (e.g., Tout & Pringle 1992, 1996; Johansen & Levin 2008) or through inward advection of external field by the accretion flow from large radii; nevertheless, these require further investigations.

Nonetheless, more extensive study is required to understand the definitive criteria or condition in launching accretion powered outflows and jets. BH spin is a very important aspect that needs to be incorporated in the conservation equations, at least through the use of pseudo-Newtonian potentials (e.g., Ghosh & Mukhopadhyay 2007; Ghosh et al. 2014, 2016), while understanding accretion powered outflow dynamics or the correlated dynamics of accretion and outflow. It is found from the work of Bhattacharya et al. (2010) that the spin of the BH significantly influences accretion powered outflows/jets; for a rapidly rotating BH, the outflow power increases by ~ two orders of magnitude. Moreover, it is also essential that to have a complete and more realistic understanding of the dynamics of the accretion-outflow coupled region, a global numerical solution of such a system in an advective paradigm should be performed. Also, explicit inclusion of cooling/radiative processes is required for the completeness of energy conservation. The relevant dynamical solutions at the accretion-outflow coupled surface would then necessarily act as boundary conditions at the base of the jet. A more definitive understanding of the criteria to launch accretion powered outflows and jets, thus, requires a complete 2.5-dimensional viscous, resistive, advective global MHD numerical solution with the inclusion of BH spin, which is left for future work. In a subsequent work (Paper II, in preparation), we will investigate in detail the energetics of magnetized accretion-induced outflows and study the spectral behavior of accretion powered sources.

The author is thankful to the anonymous reviewer for insightful comments that helped us to improve the manuscript.

The integro-differential continuity equation, Equation (7), after substitution of the flow variables in the n^th polynomial order, is given by

$$\begin{eqnarray} &&\begin{array}{l}\displaystyle \sum _{n=0}^{\infty }\left[\frac{1}{2n+1}\displaystyle \sum _{m=0}^{n}{\rho }_{2(n-m)}{v}_{r2m}\\ +\frac{1}{e+c-2n+2}\displaystyle \sum _{m=0}^{n}{\rho }_{2(n-m)}{v}_{z(2m+1)}\right]\\ {\left(\frac{h}{r}\right)}^{2n+1}=-\frac{\mathop{M}\limits^{.}}{4\pi }.\end{array}\end{eqnarray} \tag{ A.1 }$$

Similarly the polynomial expansion of all the other height-integrated MHD equations are done, however we do not furnish all of them here as the structures of the equations are huge. As an example we show it for the radial momentum balance expression, Equation (9).

$$\begin{eqnarray}&&\begin{array}{l}\left[\displaystyle \sum _{n=0}^{\infty }\displaystyle \sum _{m=0}^{n}\displaystyle \sum _{l=0}^{m}(a-2l){\rho }_{2(n-m)}{v}_{r2(m-l)}{v}_{r2l}\\ -\displaystyle \sum _{n=0}^{\infty }\displaystyle \sum _{m=0}^{n}\displaystyle \sum _{l=0}^{m}{\rho }_{2(n-m)}{v}_{\varphi 2(m-l)}{v}_{\varphi 2l}\\ +\displaystyle \sum _{n=1}^{\infty }\displaystyle \sum _{m=1}^{n}\displaystyle \sum _{l=1}^{m}2l{\rho }_{2(n-m)}{v}_{z[2(m-l)+1]}{v}_{r2l}\\ +GM\displaystyle \sum _{n=0}^{\infty }\displaystyle \sum _{m=0}^{n}\left(\begin{array}{c}-3/2\\ m\end{array}\right){\rho }_{2}(n-m)\\ +\displaystyle \sum _{n=0}^{\infty }\displaystyle \sum _{m=0}^{n}\displaystyle \sum _{l=0}^{m}(e-2n+2d){\rho }_{2(n-m)}\\ {c}_{s2(m-l)}{c}_{s2l}+\frac{1}{4\pi }[\displaystyle \sum _{n=0}^{\infty }\displaystyle \sum _{m=0}^{n}(j-2m+1)\\ {B}_{\varphi 2(n-m)}{B}_{\varphi 2m}\\ +\displaystyle \sum _{n=1}^{\infty }\underset{m=1}{\overset{n}{\displaystyle \sum }}[k-2(m-1)]\,{B}_{z[2(n-m)+1]}\\ {B}_{z[2(m-1)+1]}\\ -\displaystyle \sum _{n=1}^{\infty }\displaystyle \sum _{m=1}^{n}2m{B}_{z[2(n-m)+1]}{B}_{r2m}]\right]\\ \frac{1}{2n+1}{\left(\frac{h}{r}\right)}^{2n}=0.\end{array}\end{eqnarray} \tag{ A.2 }$$