Effects of resistivity on standing shocks in low angular momentum flows around black holes

Black hole accretion is the most efficient process which can address the issue of power generated in the neighbourhood of a black hole. Historically, the study of black hole accretion has been based on two extreme cases of accretion process: radiatively inefficient flow called Bondi flow (Bondi (1952); Michel (1972)) and radiatively efficient one called Keplerian disk (Shakura & Sunyaev (1973); Novikov & Thorne (1973)). Both suffer from certain limitations. Spherical Bondi flow with zero angular momentum is quite fast and cannot explain the high luminosities associated with observational signatures around the black hole. However, in reality, accretion flow is supposed to have some amount of angular momentum associated with it. On the other hand, cold, thin, Keplerian disk cannot explain the issue of change of spectral states and associated temporal variabilities and it is not applicable for the close region around the black hole as pressure gradient and advective radial velocity terms are ignored.

Accretion flow onto the black hole is supposed to be supersonic at the event horizon and subsonic at a large distance as the accretion flow approaches the speed of light at the horizon with sound speed being of lesser value. So the flow with angular momentum must pass through at least one sonic point before plunging onto the black hole and should be advective (Liang & Thompson (1980)). In case of accretion onto a star, even a small angular momentum will stop the matter fall onto its surface because of the infinite potential barrier associated with the Newtonian potential. Whereas in the case of a black hole, gravity always wins over centrifugal force because of higher-order terms (Chakrabarti (1993)). Not only that, for given values of specific energy and specific angular momentum of accretion flow around the black hole, multiple sonic points may also exist with the possibility of standing shocks (Fukue (1987); Chakrabarti (1989)). There have been several interesting works which explored the solutions with standing shocks in pseudo-Newtonian potential (Paczyńsky & Wiita (1980)) taking into account various prescriptions for alpha parameter (Chakrabarti (1996); Becker, Das & Le (2008); Kumar & Chattopadhyay (2013)). General relativistic solutions for inviscid (Das et al. (2015)) and viscous disks (Chattopadhyay & Kumar (2016)) with standing shocks have also been reported. Even in the presence of magnetic field, formation of standing shocks in accretion flows have been explored (Takahashi et al. (2006); Fukumura et al. (2007)).

In recent time the model which has wide recognition is the advection dominated accretion flow (ADAF) (Narayan & Yi (1994); Narayan, Kato & Honma (1997)) solution which takes care of inner boundary condition around the black hole however, has only one sonic point close to the black hole. It should be noted that advective flow with multiple sonic points may not necessarily be ADAF-type especially when standing shock exists in the accretion flow (Chakrabarti (1996)). Overall, ADAF solutions occupy a small region of parameter space for given specific energy and specific angular momentum (Lu, Gu, Yuan (1999); Kumar & Chattopadhyay (2013, 2014)). The need of a sub-Keplerian component was presented in addition to Keplerian one, the sub-Keplerian component can undergo shock transition and form a hot, puffed up region like corona (Chakrabarti & Titarchuk (1995)). The properties of post-shock region formed from the natural course of flow dynamics can address issues like state transitions (Mandal & Chakrabarti (2010)), origin of hard power-tail and low frequency QPOs (Chakrabarti, Mondal & Debnath (2015)) and also the origin of outflows (Das et al. (2001); Singh & Chakrabarti (2011); Aktar, Das & Nandi (2015)).

In the last 25 years, there have been a significant amount of simulation works dedicated to explore the formation of standing shock in low angular momentum sub-Keplerian advective flows around black holes. Using smoothed particle hydrodynamics (SPH) simulations, stable standing shocks were shown to form in one-dimensional (Chakrabarti & Molteni (1993)) and two-dimensional setups (Molteni, Lanzafame & Chakrabarti (1994)) as predicted by semi-analytical solutions of inviscid flows (Chakrabarti (1989)). For the first time, the dependence of standing shock stability on values of viscosity parameters was shown by SPH simulations as well (Chakrabarti & Molteni (1995)). The origin of outflows from the post-shock region in accretion disks were shown in simulations using SPH (Molteni, Ryu & Chakrabarti (1996)), Eulerian total variation diminishing (TVD) (Ryu, Chakrabarti & Molteni (1997); Okuda, Teresi & Molteni (2007)) and Lagrangian TVD (Lee et al. (2016)). In the presence of cooling, the post-shock region may oscillate as the cooling time scale becomes comparable to free-fall time scale and can be responsible for quasi-periodic oscillations(QPOs) in case of stellar mass as well as supermassive black holes (Molteni, Sponholz & Chakrabarti (1996); Okuda, Teresi & Molteni (2007)). Besides the case of inviscid flow, viscosity can also induce shock oscillations and give rise to QPOs (Lanzafame, Molteni & Chakrabarti (1998); Chakrabarti, Acharrya & Molteni (2004); Lanzafame et al. (2008); Lee, Ryu & Chattopadhyay (2011); Das et al. (2014); Lee et al. (2016)). There have been some works regarding stability or instability of the shock and shocks seem to be stable against axisymmetric (Nakayama (1992, 1994); Nobuta & Hanawa (1994); Le et al. (2016)) as well as non-axisymmetric perturbations (Molteni, Toth & Kunetsov (1999); Gu & Foglizzo (2003); Gu & Lu (2006)). Recently it has established through numerical simulations that advective flow can be segregated into two components, Keplerian as well as sub-Keplerian, in the presence of viscous heating and cooling processes (Giri & Chakrabarti (2013); Giri, Garain & Chakrabarti (2015); Roy & Chakrabarti (2017)). All the above-mentioned simulation works addressed the accretion flow behaviour around a non-rotating black hole using pseudo-Newtonian potential. Recently general relativistic high-resolution shock-capturing simulation code was used to study the scenario in the Schwarzschild (Kim et al. (2017)) and Kerr (Kim et al. (2019)) space-time which further established the formation of standing shock in hydrodynamic flow around non-rotating as well as rotating black hole in full general relativistic treatment. However, till now there has been only one work taking into account different magnitudes of magnetic field strength in such flows in the presence of standing shocks (Okuda et al. (2019)). The long term evolution properties were investigated and long term flares in connection with Sgr A* could be explained.

Simulation works dealing with advective flows usually take into account two types of setup : torus equilibrium solution (e.g., (Stone & Pringle (2001); McKinney, Gammie (2002)) and Bondi flow along with arbitrary choice of specific angular momentum (Proga & Begelman (2003)). Our study involves a third and different kind of set up where we take initial conditions based on exact solutions of hydrodynamic equations (Chakrabarti (1989)). We are dealing with a big gap in parameter space which lie between the regime of high angular momentum torus and zero angular momentum Bondi flow. To make our study simpler, we are dealing with inviscid flows having constant and small specific angular momentum value which are lower than that of Keplerian value of the specific angular momentum for an inner-most stable orbit. Such low angular momentum flows are likely to be present in binary systems accreting winds from companion beside Roche lobe overflow as well as active galactic nuclei where winds from stellar clusters collide and lose angular momentum before getting accreted onto central black holes (Chakrabarti & Titarchuk (1995); Smith et al. (2001, 2002); Moscibrodzka et al. (2006)). The objective of our work is to study the effects of resistivity with the varying magnitude on the formation and stability of standing shocks in low angular momentum accretion flows around the black holes which have not been explored before.

Section 2 shows details of 1.5 dimensional (1.5D) theoretical solutions which have been used for the simulation set up. In section 3, basic equations solved by simulation code are presented. Besides, the details of computational domain, initial, and boundary conditions are described in section 3. Section 4 contains details of numerical results followed by section 5 where we present a summary and discussion of our work.

We consider a semi-analytical approach of solving the standard conservation equations under hydrodynamic (HD) framework. The calculations are done in cylindrical co-ordinates with co-ordinates $R$ and $z$. Axi-symmetry is assumed for the angular $\phi$ co-ordinate. For simplicity, we further assume that the flow velocity along the vertical direction is zero and therefore only integrate along the radial co-ordinate assuming vertically averaged dynamical quantities.

We define the scale radius as Schwarzschild radius $r_{g}$= $2GM/c^{2}$, with $M$ being the mass of the central compact object, G as the gravitational constant and $c$ is the speed of light. The matter that is accreted onto the central compact object has radial velocity given by $u_{\rm R}$, specific angular momentum $L$, and total specific energy $\cal E$. As the semi-analytical calculations are carried out using non-dimensional quantities, we define the following -

For studying ideal, inviscid flow onto a compact object, we deal with mass conservation equation,

where $\rho$ is the density and $h$ being the half-thickness of flow. Energy conservation gives us a relation of specific energy of the flow or Bernoulli constant,

Here $c_{s}$ is adiabatic sound speed in units of $c$ and $\Gamma$ is the adiabatic index. $c_{s}=\sqrt{\frac{\Gamma p}{\rho c^{2}}}$. $p$ is the thermal pressure and $\Phi$ is the non-dimensional gravitational potential given by $-1/2(x-1)$ for non-rotating black hole (Paczyńsky & Wiita (1980)) where $x=r_{sp}/r_{g}$ and $r_{sp}$ being the spherical radius.

Using vertical equilibrium condition, we evaluate the radial dependence for $h$

For given values of $\epsilon$ and $\lambda$, we solve equations (2) and (3) and look for transonic conditions. Differentiating the equations (2) and (3), we obtain

where $G=\lambda^{2}/2r^{2}-1/2(r-1)$ and $f=2r^{3/2}(r-1)$
(Chakrabarti (1989)). At critical points, the vanishing of left-hand side gives radial velocity,

and the vanishing of right hand side gives sound speed,

The subscript $crit$ and $K$ represent quantities at the critical points and Keplerian orbits respectively. In the case we obtain multiple critical points, we also check whether shock conditions are satisfied or not along the accretion flow. If shock conditions are satisfied, then there is a possibility of standing axisymmetric shock otherwise not. The shock location can be determined using an invariant quantity, $C$, across the shock which is given by

Here, $\mathcal{M}=v_{\rm R}/c_{s}$ is the Mach number of the accretion flow. The subscripts $-$ and $+$ represent quantities in the pre-shock and post-shock region respectively. Further details of the semi-analytical approach can be found in (Chakrabarti (1989)).

Fig. 1 shows the variation of Mach number, $\mathcal{M}$, of flow with radial distance from black hole obtained from exact theoretical solution solving conservation equations. The transonic flow passes through the outer critical point “a” and continues its journey towards the black hole. The flow chooses to undergo shock transition along “bc”, becomes subsonic then again accelerates towards the black hole and passes through inner critical point “d” to become supersonic before entering the black hole horizon. The shocked flow is preferrable in nature as the entropy generation is relatively higher compared to no shock flow.

To get the characteristic features of the flow, we examine the time evolution of shock location $R_{s}$ on the equator and total luminosity $L$ of the flow. The luminosity $L$ is calculated as follows, assuming that the gas is optically thin

where $q_{ff}$ represents the free-free emission rate per unit volume and integration is performed over all the whole computational domain. In order to have some estimate of how much matter is lost as an outflow from the system, mass outflow rate in the z-direction is given by,

where $v_{z}$ is the vertical velocity.

In the present study, the time variability of the flow will correlate with the magnetorotational instability (MRI). We check whether the flow is subject to the MRI and whether we can resolve the fastest growing MRI mode or not. The stringent diagnostics of spatial resolution for the MRI instability has been examined in 3D magnetized flow. Therefore, its application to our 2D magnetized flow may be limited. The critical wavelength of the instability mode is given by $\lambda_{\rm c}=2\pi v_{\rm A}/\sqrt{3}\Omega$, where $v_{\rm A}$ and $\Omega$ are the Alfven velocity and the angular velocity (Hawley & Balbus (1991); Balbus & Hawley (1998)). A criterion value $Q_{\rm x}$ of the MRI resolution is defined by

where $\Delta x$ is the mesh sizes $\Delta R$ and $\Delta z$ in the radial and vertical directions, respectively. When $Q_{\rm x}\gg$ 1, the flow is unstable against the MRI instability, otherwise, the flow is stable. Fig. 3 shows 2D contours of radial MRI-criterion $Q_{r}$ for cases of $\eta=10^{-6}$ and 1.0 at times $t$= $7\times 10^{6}$ and $8.7\times 10^{6}$ seconds, respectively. The analyses show $Q_{r}\gg$ 1 in most regions for both cases. So, both flows are unstable to the MRI. The contours are asymmetric to the equator in the former but symmetric in the latter.

We studied the effect of resistivity on standing shock in the magnetized flow around a black hole. The flow parameters of specific energy, $\epsilon$ = 1.98 $\times$ $10^{-6}$ and specific angular momentum, $\lambda=1.35$, with $\Gamma$ = 1.6 have been considered to address the flow behaviour around Sgr A*. For flows with lower resistivity $\eta=10^{-6}$ and $0.01$, the luminosity and the shock location on the equator vary quasi-periodically. These quasi-periodic oscillations are attributed to the interactive result between the outer oscillating standing shock and the inner weak shocks occurring at the innermost hot blob. The luminosity varies maximumly by a factor of ten around the average $L$ $\sim 3.0\times 10^{34}$ erg s^-1. The mass outflow rate is very large as a few tens of percent of the input accretion rate. The MHD turbulence seems to play important roles in the outward transport of not only angular momentum but also accreting gas. The power density spectra of luminosity variation show the peak frequencies which correspond to the periods of 5 $\times 10^{5}$s (5.8 days), 1.4 $\times 10^{5}$s (1.6 days) and 5$\times 10^{4}$s (0.6 day) , respectively. While for cases with higher resistivity $\eta=0.1$ and 1.0 the flow becomes steady and symmetric to the equator. Variable features of the luminosity disappear here and the steady standing shock is formed more outward compared with the hydrodynamical flow. The mass outflow rate is also high as $\sim$ a few tens $\%$ of the input gas. The high resistivity considerably suppresses the magnetic activity such as the MHD turbulence and tends to form the magnetized flow to be stable and symmetric to the equator. The low angular momentum magnetized flow model with low resistivity has possibility for the explanations of the high mass outflow rate $\sim$ as 10$\%$ of the Bondi accretion rate $\sim$ 1 $\times$ $10^{-5}$ $M_{\odot}yr^{-1}$ as suggested by Chandra observations (Wang et al. (2013)) and of the long-term flares with $\sim$ one per day and $\sim$ 5 – 10 days of Sgr A* in the latest observations by Chandra, Swift, and XMM-Newton monitoring of Sgr A*.