GR Simulations of Transonic Accretion Flows

# General Relativistic Numerical Simulation of sub-Keplerian Transonic Accretion Flows onto Black Holes: Schwarzschild Spacetime

Jinho Kim, Sudip K. Garain, Dinshaw S. Balsara and Sandip K. Chakrabarti
Department of Physics, University of Notre Dame, Notre Dame, IN 46556, USA
S. N. Bose National Center for Basic Sciences, Block JD, Sector III, Salt Lake, Kolkata, 700098, India
Indian Center for Space Physics, 43 Chalantika, Garia St. Rd., Kolkata, 700084, India
E-mail: sgarain@nd.edu
Accepted XXX. Received YYY; in original form ZZZ
###### Abstract

We study time evolution of sub-Keplerian transonic accretion flows onto black holes using a general relativistic numerical simulation code. We perform simulations in Schwarzschild spacetime. We first compare one-dimensional simulation results with theoretical results and validate the performance of our code. Next, we present results of axisymmetric, two-dimensional simulation of advective flows. We find that even in this case, for which no complete theoretical analysis is present in the literature, steady state shock formation is possible.

###### keywords:
accretion, accretion discs – black hole physics – hydrodynamics – shock waves – methods: numerical
pubyear: 2017pagerange: General Relativistic Numerical Simulation of sub-Keplerian Transonic Accretion Flows onto Black Holes: Schwarzschild SpacetimeGeneral Relativistic Numerical Simulation of sub-Keplerian Transonic Accretion Flows onto Black Holes: Schwarzschild Spacetime

## 1 Introduction

Independent of the source of matter supply, an accretion flow around a black hole is necessarily transonic, i.e., it must pass through one or more sonic point(s). An accretion flow on a gravitating star, whose specific angular momentum is everywhere or almost everywhere lower than that of the local Keplerian value may be termed as a sub-Keplerian flow. Theoretical calculations as well as numerical simulations of sub-Keplerian transonic accretion flows around black holes have shown that such flows can have standing shocks for appropriate choice of the flow parameters such as the specific energy and specific angular momentum of the accreting material (Chakrabarti 1989, 1990a, hereafter C89, C90, respectively). Paradoxical as it may sound, the flow decides to slow down just a few to few tens of gravitational radii away from the horizon, simply because the centrifugal barrier becomes very strong as compared to the gravity. However, the barrier is not unsurmountable and flow simply passes through a shock transition to satisfy boundary conditions, as is normal in many astrophysical circumstances. The shock can be simply standing or propagating away depending on the flow parameters and the dissipative and cooling processes present in the flow (\al@chakraba1989,chakraba1990; \al@chakraba1989,chakraba1990; Chakrabarti 1996a, and references therein). Using numerical simulations with Smoothed Particle Hydrodynamics (SPH), these shocks have been shown to be stable in both one and two dimensional simulations (Chakrabarti & Molteni, 1993; Molteni, Lanzafame, & Chakrabarti, 1994, hereafter CM93, MLC94, respectively). They were found to propagate away for high enough viscosity (Chakrabarti & Molteni, 1995). The post-shock region, which is known as CENtrifugal pressure supported BOundary Layer (or CENBOL), as in any other astrophysical flows, is found to extremely useful to explain detailed spectral properties of the black hole accretion disc quite satisfactorily (Chakrabarti & Titarchuk, 1995; Chakrabarti, 1997, hereafter CT95, C97, respectively). The so-called Two-Component Advective Flow (TCAF) model, proposed in CT95 and C97, describes the most general structure of an accretion disc which consists of two transonic components, namely, an optically thick (optical depth, ), geometrically thin Keplerian disc and an optically thin (), geometrically thick sub-Keplerian flow. Shock is formed in the sub-Keplerian component which has higher radial infalling velocity than the other component. In the post-shock region, i.e., inside the CENBOL, these two components mix up due to turbulence and heat (CT95) and form an optically slim (), geometrically thick disc and continues its journey towards the central object. The CENBOL itself becomes responsible for the non-thermal power-law component of the observed spectrum.

Recently, TCAF model has been included in HEASARC’s spectral analysis software package XSPEC (Arnaud, 1996) and is being used to study the spectral as well as timing properties of several black hole candidates (Debnath, Chakrabarti & Mondal, 2014; Mondal, Debnath & Chakrabarti, 2014; Debnath, Mondal & Chakrabarti, 2015; Chakrabarti, Mondal, & Debnath, 2015; Debnath et al., 2015; Jana et al.,, 2016; Chatterjee et al.,, 2016). It has been demonstrated that this model explains the observed data very well. In outbursting sources, the shocks were found to be propagating as well - steadily moving towards the black hole in the rising phase and moving away from it in the declining phase. Simulations of inviscid accretion flow with presence of a power-law cooling show that shock can oscillate about a certain mean location, particularly when there is resonance between the cooling and the infall time scales. This may be considered to be an explanation of the origin of low frequency quasi-periodic oscillations (LFQPOs; Molteni, Sponholz & Chakrabarti, 1996; Chakrabarti, Acharyya & Molteni, 2004). Such oscillations have been shown to be present and stable when real Compton cooling is present in the flow as well (Garain, Ghosh & Chakrabarti, 2012, 2014, hereafter, GGC12, GGC14, respectively). These intriguing properties of the CENBOL clearly demand conducting robust numerical experiments, which we set out to do in a series of papers.

Several numerical experiments are already present in the literature which study shock structures in the sub-Keplerian flow around black holes. For one and two-dimensional inviscid, adiabatic flows, it has been shown that a sub-Keplerian axisymmetric flow with and without shock is stable (CM93; MLC94; Molteni, Ryu, & Chakrabarti 1996, hereafter MRC96). New codes have been tested against such non-linear solutions as well (Toth, Keppens & Botchev, 1998). Numerical simulations of adiabatic viscous flow also demonstrate the stability of these standing shock waves (Giri & Chakrabarti, 2012; Lee et al., 2016). More recently, two-dimensional axisymmetric numerical simulations of viscous flow in presence of power-law and Compton cooling show that an advective flow actually splits into two components when appropriate viscosity parameters and cooling processes are chosen (Giri & Chakrabarti, 2013; Giri, Garain & Chakrabarti, 2015, hereafter, GC13, GGC15, respectively). As the sub-Keplerian flow advects towards the central black hole, angular momentum transport and condensation due to cooling on the equatorial plane help the flow to segregate into two distinct advective components, each being separately transonic. The component near the equatorial plane has been shown to have the properties very similar to a standard Keplerian disc (Shakura & Sunyaev, 1973). This component is surrounded by the sub-Keplerian component which will have the steady shock. Thus, these results show that a stable TCAF formation is indeed possible.

In case of magnetized accretion disc, magnetorotational instability (MRI; Balbus & Hawley, 1991; Hawley & Balbus, 1992) can make the disc turbulent and this turbulence may transport angular momentum outwards efficiently. However, it has been shown that turbulence triggered by MRI produce the value of the alpha viscosity parameter 0.01 (Brandenburg et al., 1995; Hawley et al., 1995, 1996; Arlt & Rüdiger, 2001; Smak, 1999; King, Pringle & Livio, 2007; Kotko & Lasota, 2012). In the literature, many authors from various groups published results of analytical viscous solutions as well as viscous hydro-dynamic simulations where it has been shown the shock is stable if the viscosity parameter is lower than a critical viscosity parameter (Chakrabarti, 1990b, 1996d; Chakrabarti & Molteni, 1995; Lanzafame, Molteni & Chakrabarti, 1998; Lanzafame et al., 2008; Lee, Ryu & Chattopadhyay, 2011; Das et al., 2014; Lee et al., 2016). Also, there have been many stability studies of shock (Nakayama, 1992, 1994; Nobuta & Hanawa, 1994; Gu & Foglizzo, 2003; Gu & Lu, 2006), but it was shown that even under non-axisymmetric perturbations, the shock tends to persist, albeit, as a deformed shock (Molteni, Tóth, & Kuznetsov, 1999).

Despite all these important developments, almost all the above mentioned simulations have been performed using the so-called pseudo-Newtonian potential (Paczyński & Wiita, 1980) which mimics the Schwarzschild spacetime. This potential retains the particle properties of the Schwarzschild geometry in the sense that the marginally bound and marginally stable orbits are located at exactly the same places as in GR. However, several properties in the strong gravity limit just outside of the horizon are not retained with precision. For instance, the energy released at the marginally stable orbit, or the velocity of matter on the horizon are different. Thus, though the results were satisfactory, it was difficult to judge if any of the effects observed were artefacts of the potential. Fortunately, even in general relativistic (GR) framework, in Kerr spacetime, such shocks have been shown to exist (C90; Chakrabarti 1996b, hereafter C96b). And it became clearer that the shocks in black hole accretion are indeed possible only because of the presence of the inner sonic point between the marginally bound and marginally stable orbits where strong gravity is important. Strong gravity forces the flow to have an inner sonic point. Most importantly in accretion flow configurations, the solution passing through the inner sonic point has higher entropy than that passing through the outer sonic point. So the flow generates entropy at the shock and then passes through the inner sonic point. The location of the sonic point is the closest indicator of a stable fluid (unlike the marginally stable orbit, which is the closest indicator for a particle trajectory). Thus there are all the more reasons to verify these results using a full relativistic framework. In this work, we perform a general relativistic simulation of the sub-Keplerian flow in Schwarzschild spacetime. To our knowledge, no numerical experiment has so far been performed which tests the possibility of a stable CENBOL formation in general relativity.

This paper is organized as follows. In Section 2, we present the analytical method to calculate the one-dimensional flow properties in Schwarzschild spacetime. In Section 3, we present the general relativistic equations which are solved numerically and the numerical procedure we use for doing this. In the next Section, we present the results for one and two dimensional simulations. Finally, we present our conclusions.

In this paper, we choose as the unit of distance, as unit of angular momentum, and as unit of time. In addition, we choose the geometric units ( is gravitational constant, is the mass of the black hole and is the unit of light). Thus , and angular momentum and time are measured in dimensionless units.

## 2 Analytical Solution

The general relativistic study of transonic flows has been done extensively by Chakrabarti (\al@chakraba1990,chakraba1996b; \al@chakraba1990,chakraba1996b). Therefore, we do not describe all the details here. However, for completeness, we only mention the important equations for Schwarzschild spacetime.

For this calculation, we use Boyer-Lindquist coordinates (). The line element in Schwarzschild spacetime is given as follows,

 ds2=gμνdxμdxν=−(1−2r)dt2+(1−2r)−1dr2+r2dθ2+r2sin2θdϕ2 (1)

We are interested in the flow close to the equatorial plane, so is assumed for analytical study.

In absence of viscosity and any heating or cooling, one can find the conserved specific energy as (C96b)

 ϵ=hut=11−na2ut, (2)

where, is the polytropic index, being the adiabatic index and is the enthalpy, being the sound speed. Also,

 ut=[1−2r(1−V2)(1−Ωl)]1/2. (3)

Here,

 Ω=uϕut=−lgttgϕϕ=lr2(1−2r), (4)

and is the specific angular momentum. Also,

 V=v(1−Ωl)1/2, (5)

where

 v=(−ururutut)1/2. (6)

The entropy accretion rate (\al@chakraba1989,chakraba1996b; \al@chakraba1989,chakraba1996b) is given by

 ˙μ=(a21−na2)nV(1−Ωl)1/2utr2 (7)

We follow the usual solution procedures use in transonic flows (\al@chakraba1989,chakraba1990; \al@chakraba1989,chakraba1990) to calculate and radial dependence of other required quantities. By differentiating equations (2) and (7) with respect to and eliminating terms involving , we find following expression as the gradient of :

 dVdr=V(1−V2)[1−2ra2+3a2−Ωl1−Ωl(r−3)]r(r−2)(a2−V2) (8)

We can readily see that for , we recover the similar expression for Bondi flow onto a black hole (Eq. 1.29 of C90). At the sonic point, both numerator and the denominator vanish and one obtains the so-called sonic point condition as

 Vc=ac,anda2c=12rc−3[1−Ωl1−Ωl(r−3)], (9)

where, is called the sonic radius.

To find a complete solution from the horizon to infinity, one needs to supply the specific energy and the specific angular momentum . If the supplied parameters allow the accretion solution to have a shock in the sense that Rankine Huguniot conditions are satisfied, then the shock location can be found by determining a constant, , which remains invariant across the shock (\al@chakraba1989,chakraba1990; \al@chakraba1989,chakraba1990). It is found that this invariant quantity is the same in Schwarzschild space-time and in pseudo-Newtonian potential (C90), as this is a local equation. Therefore, we use the following expression,

 C=[ΓM+(1/M)]22+(Γ−1)M2 (10)

to determine the shock location for this calculation (Chakrabarti & Das, 2001). Here, we use as the definition of Mach number.

## 3 Numerical Simulation Procedure

We use so called “Valencia formulation” to numerically solve the relativistic hydrodynamic equation (Banyuls et al., 1997). This formulation gives flux conservative form of the system of hydrodynamics equation in the framework of 3+1 formalism. It has been applied very successfully in computational fluid dynamics. In our coordinate system, the conservative variables () and primitive variables () are

 q=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝DSrSθSϕτ⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠≡⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝ρWρhW2vrρhW2vθρhW2vϕρhW2−P−D⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠,w=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝ρvrvθvϕP⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠. (11)

Here is the fluid rest mass density is the pressure and is the specific enthalpy. They are measured in the comoving frame of the fluid. is the fluid velocity measured by Eulerian observer. is the Lorentz factor and defined as . Here, are the spatial part of the metric components . The radial and angular velocity in the Eulerian frame can be expressed in terms of in Eq. (6) and in Eq. (4):

 vr=urW=(1−2r)12vvϕ=uϕW=(1−2r)−12Ω (12)

Assuming axisymmetry (), the hydrodynamical equations in the curved spacetime that described in Eq. (1) can be written as follows:

 ∂(√γq)∂t+∂(√−gfr)∂r+∂(√−gfθ)∂θ=√−gΣ, (13)

where

 fr = ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣DvrSrvr+PSθvrSϕvrτvr+Pvr⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦, fθ = ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣DvθSrvθSθvθ+PSϕvθτvθ+Pvθ⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦, Σ = ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣0−ρhW2r(1r−2(1+vrvr)−vθvθ−vϕvϕ)+2Prcotθ(ρhW2vϕvϕ+P)0−ρhW2vrr(r−2)⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦. (14)

Here and are the determinants of spatial and spacetime metric, respectively. From the Schwarzschild metric shown in Eq. (1), we have , and . For the spherically symmetric cases, . Eq. (14) consists of the continuity, three momentum and energy equations. We clearly see the conservation of total rest mass (baryon number) and angular momentum in the first and fourth rows of Eq. (14). Particularly, the terms of in the momentum equation contains the gravitational and the centrifugal forces. The gravitational force which is purely radial in the Schwarzschild metric is shown in the first term of the second row of . The second and third terms in the second row as well as the first term in the third row represents the centrifugal force by the rotation velocity. Since the centrifugal force, exerted by the , is not purely radial, it contributes to both and momentum equations. (c.f. The centrifugal force exerted by the is purely radial. Therefore, it only appears in the radial momentum equation.) Note that the last terms in the of the momentum equations are the additional terms in the spherical polar coordinates system. The fifth row of is the sink or source of the energy contributed by gravity.

We use the ideal gas equation of state which can be written in the following form:

 P=(Γ−1)ρe, (15)

where is the specific internal energy. The above equation of state provides the expression of specific enthalpy:

 h=1+ΓΓ−1Pρ. (16)

In this paper, we solve the above hydrodynamic equations using a numerical code developed by Kim et al. (2012). The details of the code can be found in Kim et al. (2012). The most useful property of this code is that it can be applied to any spacetime metric in any coordinate system. This code uses finite volume methods to ensure local conservation of the fluid in the computational grid. Therefore, the code can guarantee total mass and angular momentum conservations which appear in the first and fourth rows of Eq. (14). For the treatment of the discontinuous behaviors of the fluid such as shocks, rarefactions or contact discontinuities, the High Resolution Shock Capturing (HRSC) techniques are applied in the code. We use the third order slope limiter proposed by Shibata (2003) which is based on the minmod function. For the flux approximation, we use the HLL method (Harten et al., 1983). The HLL method has some dissipation but the results are very stable. For the time integration, we use the third order three stage Strong Stability-Preserving (SSP) Runge-Kutta method which is known as Shu-Osher method (Shu & Osher, 1988).

## 4 Results

We use inflow boundary condition at the outer boundary located at . The inner boundary is placed at and we use the extrapolated values of the primitive variables for the inner ghost cells. These inner ghost cells are located outside of the event horizon. In this paper, we present results of one-dimensional and two-dimensional simulations. For better resolution close to the central black hole, we logarithmically binned the radial direction in 300 zones for all the simulation results presented here. For two-dimensional simulations, in addition to above mentioned radial binning, we used 100 equi-spaced zones in polar direction.The innermost radial zone () has a size of and the outermost zone has a size of 1.28.

For these simulations, we need values of the primitive variables (see, Eq. (11)) at . For all the simulations, we set . The density of the incoming matter at is normalized to . In absence of self-gravity and heating or cooling, the density is scaled out and the simulation results remain valid for any accretion rate (MRC96; Giri at al., 2010, hereafter GCSR10). We evaluate by solving Eq. (8) for a given pair of conserved flow variables. Subsequently, the sound speed is evaluated using Eq. (2). Next, using Eqs. (5), (6) and (12), we evaluate and . Pressure of the incoming matter is evaluated from the sound speed using Eq. (16). These values are maintained in the ghost zones of the outer boundary of our computational domain. As an initial condition, we put floor values for the density to be (in normalized unit) and the corresponding floor value of pressure inside the computational domain. Floor value of the pressure is chosen such that the sound speed (or temperature) of the background matter is same as that of incoming matter (\al@mrc96,gc2013; \al@mrc96,gc2013). Thus, once we know , we evaluate the value of the pressure floor using Eq. (16) by substituting in this equation. Then, the value of pressure floor is . Initially, the velocity components in the floor grids are set to zero. Thus, initially, as the matter rushes towards the black hole, it fills the vacuum rapidly. After the matter reaches the inner boundary, the flow starts to feel the pressure and centrifugal force.

### 4.1 One-dimensional, spherically symmetric Bondi accretion

For a given , the analytical structure of a spherically symmetric Bondi accretion flow is completely determined ( for this flow). For the comparative study, we choose . This choice gives and at , and the sonic point at . Note that, just the energy was sufficient to determine all the quantities as the other conditions come from transonicity (Eq. (9); C90).

In Fig. 1, we compare the radial variations of Mach number, defined as . The black solid line represents the analytical result, whereas the red plus signs represent the numerical simulation result. Clearly, they are indistinguishable. Dynamical time, (defined as the time required for matter to reach inner boundary from in steady state) is found to be for this simulation and we ran the simulation for more than (). Also, we ran some more cases by varying the number of radial zones ranging from 150 to 600 and verified that the results remained converged. Note that the inner boundary places at for the lower resolution case in order to prevent the inner ghost cells from locating inside the event horizon.

### 4.2 One-dimensional accretion flow with non-zero angular momentum

When angular momentum is present, the flow structure changes significantly depending on its strength. As discussed earlier, depending on the values of and , the accretion flow may or may not have any shock. In this Section, we present the results of two simulations in one-dimension. One solution has a shock and the other does not. As discussed in C89 (see Fig. 3); C90 (see Chap. 3 and 6), the accretion solution having a shock first passes through the outer sonic point, makes a transition to the sub-sonic branch at the shock and then passes through the inner sonic point just before being supersonically accreted by the black hole. On the other hand, when and are such that it is not possible to have a shock, the flow passes only through the outer or inner sonic point before disappearing behind the horizon. However, such a flow still can slow down as it approaches the black hole because of the centrifugal force.

In Fig. 2(a), we present the radial variation of Mach number, , at the final steady state for the accretion flow which has a shock. As before, the black solid line represents the analytical result and the red plus signs represent the simulation result. and are chosen for this simulation. Analytical calculation gives and at . The outer sonic point, shock and inner sonic point are located at , and respectively. As can be seen, our 1D simulation has captured all the locations properly. For this case, dynamical time is found to be 1150 and the simulation was run till (). We have also computed and from the simulation result and verified that these two quantities remain constant along the flow.

Fig. 2(b) shows the radial variation of for a case where there is no shock in the accretion solution, and is expected to pass only through the outer sonic point. We choose and to compute the outer boundary values for this simulation. For this case, dynamical time is found to be 920 and the simulation was run till (). The outer sonic point is located at for these parameters. As can be seen from the Figure, analytical and numerical results are nearly inseparable and hence, we see that numerical simulation has captured this location very well. We may mention in passing that the simulations presented in the literature are primarily for supersonic injection to save computational time. However, our efficient method allows us to inject the flow at sub-sonic Mach number.

In order to verify that the shock location is not affected by the outer boundary, we rerun the case presented in Fig. 2(a) by moving to 200, 500, 1000 and 2000. All these simulations have been run using 300 zones in the radial direction. In Fig. 3, we present these simulation results. Green, blue, magenta and cyan (cross, star, open squrare and filled square) points represent the simulation results for respectively. In the inset, we show the zoomed in part around the shock location. Comparison of shock location for various shows that it does not get affected if we change . Slight mismatch for and may be due to the grid size variation at the shock location.

### 4.3 Two-dimensional simulations

Realistically, an accretion disc is three dimensional. However, assuming axisymmetry we can study the disc structure in two-dimensions. On the other hand, the theoretical formalism of transonic flows (Chakrabarti, 1996b), is developed for flows in vertical equilibrium which is quite thin. Indeed, the shock location and Mach number variation depend on the model assumption. Most interestingly, it was shown (GCSR10) that the pre-shock flow behaves as a conical flow, while the post-shock behaves as a flow in vertical equilibrium. So a direct comparison with theoretical result is not possible in a two dimensional simulation. On the other hand, in subsection 4.2, we showed that our simulation result matches with the theory very well. Thus the results obtained in two dimensions may be trusted.

For the results presented here, we do the simulation in coordinates. Simulation domain extends from to in radial direction and in polar direction. The incoming matter enters the simulation box at through one-tenth of polar zones starting from equatorial plane (\al@gcsr2010,ggc2012,ggc2014; \al@gcsr2010,ggc2012,ggc2014; \al@gcsr2010,ggc2012,ggc2014). For this simulation, we evaluate the outer boundary values using vertical equilibrium model. We choose and , and this choice gives and at . For these parameters, the one-dimensional analytical calculation predicts the shock to be located at . However, in a full 2D simulation, presence of turbulence due to centrifugal barrier is expected to shift the shock farther out. As the incoming matter hits the centrifugal barrier, some matter bounce back and interact with the incoming matter. Thus a turbulence is generated. This turbulent pressure seems to be comparable to the other pressure effects, such as thermal and ram pressure (\al@mlc94,mrc96; \al@mlc94,mrc96). This turbulent pressure shifts the location of the shock further out compared to the location that we calculate using one-dimensional analytical method. Interestingly this is precisely what we see.

In Fig. 4, we present the contours of normalized density inside the accretion disc at four different times. Associated color bar represents the values of the normalized density. Velocity vectors are over-plotted with density. The length of a vector is proportional to the magnitude of velocity at that location. The simulation is carried out till the time of 30000 and we measure the dynamical time to be . Thus, the simulation continued for at least 20 dynamical times without any significant time evolution and hence, we believe that the system has reached a steady state.

Fig. 4(a) shows the density and velocity vectors at time 158 (), a little after the simulation started. Fig. 4(b) shows the same at time 496 (), soon after matter reaches . It can be seen that a shock structure already started forming by this time. The shock front can be identified by the jump in density color contour close to the axis and also a sudden reduction of velocity vector lengths. The contours clearly look like those of a thick accretion disc as discussed in MLC94 also for a pseudo-Newtonian simulation. As we move up in the vertical direction from the equatorial plane, the shock front bends outward. This is explained in MRC96; Chakrabarti 1996c as due to the reduction of gravitational pull with height, while the centrifugal force which remains almost the same, pushes the shock outward. This shock front moves away radially from the black hole and stabilizes at a radial distance of . This can be seen in Fig. 4(c) which is a snapshot at nearly 1.5 dynamical times, i.e., 2250. Fig. 4(d) shows the snapshot at a final time of 30000 (). There appears to be practically no change in the shock location or the flow behaviour at this period of the simulation. We ran this simulation at resolution and found that the results presented here are converged. We also ran this simulation at resolution by placing at 2.5. Even for this lower resolution, the shock location and the overall structure of the disc is found to be very similar to the presented results. Velocity vectors in Figs. 4(c) and 4(d) also show that a strong outflow emerges from the post-shock region or CENBOL. In order to prove that it indeed leaves the system supersonically, we plot in Fig. 5 the contours of radial Mach number. As the color code would indicate, the inward flow became highly supersonic in the pre-shock region and highly subsonic in the post-shock region. It became supersonic again closer to the black hole as it moves to satisfy the boundary condition on the horizon. The outflow behaves exactly the opposite way. It starts subsonically from CENBOL surface and becomes supersonic by the time it reaches at about ten Schwarzschild radii. The raggedness of the contours on the funnel like vortex surface close to the axis is due to the artefact of finite sized grids whose size increase with radial distance from the black hole.

In Fig. 6, we present the time evolution of the shock location close to the equatorial plane. We dynamically compute the shock location by picking up the position where the Mach number changes from to . We have omitted the initial transient time from this plot when the floor grids were being filled by the incoming flow for the first time. Fig. 6 shows that the shock location has achieved a steady state, although it oscillates slightly about its mean location of , possibly due to the finite grid size at this distance.

## 5 Conclusions

It has been demonstrated that TCAF model explains spectral and temporal properties of the black hole accretion discs very well. The CENBOL region which is primarily responsible for emission of hard photons and the outflows, relies on the formation of shock in the sub-Keplerian component in this model. However, the formation of CENBOL was taken for granted only using simulations carried out in pseudo-Newtonian geometry although its genesis lies in properties of an advective flow in the strong gravity limit. In this paper, we performed numerical simulations in fully general relativistic framework in Schwarzschild spacetime to see whether a steady shock formation is still possible or not. We presented results of one-dimensional and two-dimensional simulations. For one-dimensional simulations, we tested the code by comparing our results with the steady state results obtained using analytical methods. We ran our code for sufficiently long time (over 20 dynamical time) the steady state was found to be achieved. Comparisons of radial Mach number variations were shown for spherically symmetric Bondi accretion flow as well as for axisymmetric accretion flow with and without shock. We found a very good match between the two methods even at moderate resolution and these validated the performance of our simulation code. This gave us confidence to delve into unchartered territory of running the code in two-dimensions, where, strictly speaking, a fully self-consistent theoretical result was missing. Here again, we found steady shock formation and the post-shock region can be easily identified to be the CENBOL region used by CT95. We found that a centrifugal force driven supersonic jet is formed from the surface of the CENBOL. We also plotted the time variation of the shock location close to the equatorial plane. We found that the location slightly oscillates about a mean value of , but overall structure is steady.

The results presented here were obtained from purely hydrodynamical simulations in Schwarzschild spacetime. Analytically, it has been shown in \al@chakraba1990,chakraba1996b; \al@chakraba1990,chakraba1996b that the presence of black hole rotation will affect the structure of an accretion disc. The parameter space of the shock formation is different for prograde and retrograde flows, and locations of the sonic points as well as shock can be closer or away depending on whether the so-called “spin-orbit” coupling term (arising out of the product of the spin vector of the black hole and the orbital angular momentum vector of the matter) is positive or negative. We will numerically study the effects of black hole rotation on the accretion disc and shock behaviour. Our code would be easily used to study the dragging of inertial effects also which is expected to affect the outflows and CENBOL behaviour. Furthermore, in presence of dissipations, such as viscosity and/or heating/cooling, the structure and the location of the CENBOL will change significantly as shown by previous simulations (\al@ggc2012,ggc2014,gc2013,ggc2015; \al@ggc2012,ggc2014,gc2013,ggc2015; \al@ggc2012,ggc2014,gc2013,ggc2015; \al@ggc2012,ggc2014,gc2013,ggc2015), performed using pseudo-Newtonian potential. We will explore the effects of these dissipations in future studies using our general relativistic numerical simulation codes and study how the segregation of the flow into two components may happen in general relativistic framework.

## Acknowledgements

DSB acknowledges support via NSF grants NSF-DMS-1361197, NSF-ACI-1533850, NSF-DMS-1622457. Several simulations were performed on a cluster at UND that is run by the Center for Research Computing. Computer support on NSF’s XSEDE and Blue Waters computing resources is also acknowledged.

## References

• Arlt & Rüdiger (2001) Arlt, R., & Rüdiger, G. 2001, A&A, 374, 1035
• Arnaud (1996) Arnaud K.A., 1996, ASP Conf. Ser., Astronomical Data Analysis Software and Systems V, ed. G.H. Jacoby & J. Barnes, 101, 17
• Balbus & Hawley (1991) Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214
• Banyuls et al. (1997) Banyuls F., Font, J. A., Ibáñez, J. M., Martí, J. M., & Miralles, J. A. 1997, ApJ, 476, 221
• Brandenburg et al. (1995) Brandenburg A., Nordlund A., Stein R. F., & Torkelsson U., 1995, ApJ, 446, 741
• Chakrabarti (1989) Chakrabarti S. K., 1989, ApJ, 347, 365 (C89)
• Chakrabarti (1990a) Chakrabarti S. K., 1990, Theory of Transonic Astrophysical Flows, World Scientific, Singapore (C90)
• Chakrabarti (1990b) Chakrabarti S. K., 1990, MNRAS, 243, 610
• Chakrabarti (1996a) Chakrabarti S. K., 1996a, ApJ, 464, 664 (C96a)
• Chakrabarti (1996b) Chakrabarti S. K., 1996b, MNRAS, 283, 325 (C96b)
• Chakrabarti (1996c) Chakrabarti S. K., 1996c, ApJ, 471, 237
• Chakrabarti (1996d) Chakrabarti S. K., 1996, Phys. Rep., 266, 229
• Chakrabarti (1997) Chakrabarti S. K., 1997, ApJ, 484, 313
• Chakrabarti & Molteni (1993) Chakrabarti S. K., Molteni D., 1993, ApJ, 417, 671
• Chakrabarti & Molteni (1995) Chakrabarti S. K. & Molteni D., 1995, MNRAS, 272, 80
• Chakrabarti & Titarchuk (1995) Chakrabarti S. K., & Titarchuk L. G., 1995, ApJ, 455, 623 (CT95)
• Chakrabarti & Das (2001) Chakrabarti S. K. & Das S., 2001, MNRAS, 327, 808
• Chakrabarti, Acharyya & Molteni (2004) Chakrabarti S. K., Acharyya K., Molteni D., 2004, A&A, 421, 1
• Chakrabarti, Mondal, & Debnath (2015) Chakrabarti S.K., Mondal S., & Debnath D., 2015, MNRAS, 452, 3451
• Chatterjee et al., (2016) Chatterjee D., Debnath D., Chakrabarti S. K., Mondal S., & Jana A., 2016, ApJ, 827, 88
• Das et al. (2014) Das S., Chattopadhyay I., Nandi A., Molteni D., 2014, MNRAS, 442, 251
• Debnath, Chakrabarti & Mondal (2014) Debnath D., Chakrabarti S.K., & Mondal S., 2014, MNRAS, 440, L121
• Debnath, Mondal & Chakrabarti (2015) Debnath D., Mondal S., & Chakrabarti S.K., 2015, MNRAS, 447, 1984
• Debnath et al. (2015) Debnath D., Molla A. A., Chakrabarti S .K., & Mondal S., 2015, ApJ, 803, 59
• Garain, Ghosh & Chakrabarti (2012) Garain S. K., Ghosh H., Chakrabarti S. K., 2012, ApJ, 758, 114
• Garain, Ghosh & Chakrabarti (2014) Garain S. K., Ghosh H., Chakrabarti S. K., 2013, MNRAS, 437, 1329
• Giri at al., (2010) Giri K, Chakrabarti S. K., Samanta M. M., & Ryu D., 2010, MNRAS, 403, 516
• Giri & Chakrabarti (2012) Giri K., Chakrabarti S. K., 2012, MNRAS, 421, 666
• Giri & Chakrabarti (2013) Giri K., Chakrabarti S. K., 2013, MNRAS, 430, 2836
• Giri, Garain & Chakrabarti (2015) Giri K., Garain S. K., & Chakrabarti S. K., 2015, MNRAS, 448, 3221
• Gu & Foglizzo (2003) Gu, W. M., Foglizzo, T. 2003, A&A, 409, 1
• Gu & Lu (2006) Gu, W. M., Lu, J. F. 2006, MNRAS, 365, 647
• Harten et al. (1983) Harten Amiram, Lax Peter D & van Leer, Bram  1983, SIAM Review, 25, 35
• Hawley & Balbus (1992) Hawley, J. F., & Balbus, S. A. 1992, ApJ, 400, 595
• Hawley et al. (1995) Hawley J. F., Gammie C. F., Balbus S. A., 1995, ApJ, 440, 742
• Hawley et al. (1996) Hawley J. F., Gammie C. F., Balbus S. A., 1996, ApJ, 464, 690
• Jana et al., (2016) Jana A., Debnath D., Chakrabarti S. K., Mondal S., & Molla, A. A., 2016, ApJ , 819, 107
• Kim et al. (2012) Kim J., Kim H. I., Choptuik M. W., & Lee H. M. 2012, MNRAS, 424, 830
• King, Pringle & Livio (2007) King, A. R., Pringle, J. E., & Livio, M. 2007, MNRAS, 376, 1740
• Kotko & Lasota (2012) Kotko, I., & Lasota, J.-P., 2012, A&A, 545, 9
• Lanzafame, Molteni & Chakrabarti (1998) Lanzafame, G., Molteni, D., & Chakrabarti, S. K. 1998, MNRAS, 299, 799
• Lanzafame et al. (2008) G. Lanzafame, P. Cassaro, F. Schilliró, V.Costa, G.Belvedere & R. A. Zappalá, 2008, A&A, 482, 473
• Lee, Ryu & Chattopadhyay (2011) Lee S. J., Ryu D., & Chattopadhyay I., 2011, ApJ, 728, 142
• Lee et al. (2016) Lee S, Chattopadhyay I, Kumar R, Hyung S., & Ryu D., ApJ, 831, 33
• van Leer (1977) van Leer, B. 1977, Journal of Computational Physics, 23, 276
• Molteni, Lanzafame, & Chakrabarti (1994) Molteni D., Lanzafame G., Chakrabarti S. K., 1994, ApJ, 425, 161 (MLC94)
• Molteni, Ryu, & Chakrabarti (1996) Molteni D., Ryu D., Chakrabarti S. K., 1996, ApJ, 470, 460 (MRC96)
• Molteni, Sponholz & Chakrabarti (1996) Molteni D., Sponholz H., Chakrabarti S. K., 1996, ApJ, 457, 805
• Molteni, Tóth, & Kuznetsov (1999) Molteni, D., Tóth, G., & Kuznetsov, O. A. 1999, ApJ, 516, 411
• Mondal, Debnath & Chakrabarti (2014) Mondal S., Debnath D., & Chakrabarti S. K., 2014, ApJ, 786, 4
• Nakayama (1992) Nakayama, K. 1992, MNRAS, 259, 259
• Nakayama (1994) Nakayama, K. 1994, MNRAS, 270, 871
• Nobuta & Hanawa (1994) Nobuta, K., & Hanawa, T. 1994, PASJ, 46, 257
• Paczyński & Wiita (1980) Paczyński B., Wiita P. J., 1980, A&A, 88, 23
• Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
• Shu & Osher (1988) Shu C.-W., & Osher S. 1988, Journal of Computational Physics, 77, 439
• Shibata (2003) Shibata M. 2003, Phys. Rev. D, 67, 024033
• Smak (1999) Smak, J. 1999, AcA, 49, 391
• Toth, Keppens & Botchev (1998) Toth G., Keppens R, & Botchev M. A., 1998, A&A, 332, 1159
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters