{'fraction_non_alphanumeric': 0.03734691722972973, 'fraction_numerical': 0.01272170608108108, 'mean_word_length': 4.6416021441334125, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Transonic behaviour and stability analysis of quasi-viscous black hole accretion

Analytical studies of black hole accretion usually presumes the stability of the stationary transonic configuration. Various authors in the past several decades demonstrated the validity of such an assumption for inviscid hydrodynamic flow. Inviscid approximation is a reasonable approach for low angular momentum advection dominated flow in connection to certain supermassive black holes at the centres of the galaxies (including our own) fed from a number of stellar donors. Introduction of a weak viscosity, as a first order linear correction involving the viscosity parameter, however, may sometimes provide a more detail understanding of the observed black hole spectra. The transonic behaviour of the stationary solutions have been studied for the aforementioned quasi-viscous accretion for all possible geometric configurations of axisymmetric flow. For a sufficiently low range of the viscosity parameter, transonic solutions containing one or three critical points have been found for allowed ranges in the astrophysical parameters under the post-Newtonian pseudo-Schwarzschild scheme. With the introduction of such viscosity parameter, the only feasible critical points are of saddle and spiral types in contrary to the inviscid case where centre type points were formed instead of spiral ones. Introduction of linear perturbations on stationary flow solutions and their time evolution in both standing and radially propagating wave forms have been examined (completely analytically) in detail. Our analysis shows that similar kind of secular instability exists in all the considered disk models at large distance in the asymptotic limit, however, the model itself is valid only within a certain length scale and the disks sustain within that length scale only for at least a considerable time scale.


Introduction
Viscous transport of angular momentum is a necessary phenomena to allow the matter to fall in through the event horizon for a rotation dominated axially symmetric accretion flow onto astrophysical black holes [1,2]. A proper viscosity prescription is imperative for the study of such flow having the Keplerian distribution of angular momentum.
Forty two years after the discovery of the standard accretion disk model [3,4], exact modeling of viscous transonic black hole accretion (including appropriate heating and cooling mechanism), however, is still quite an arduous task. This is true even for flow dynamics analyzed from a purely Newtonian framework -let alone for the general relativistic accretion disks. Over the past several years, various authors have pointed out that for low angular momentum advection dominated accretion (which can model the flow structure onto Sgr A*, for instance), large radial velocity of matter close to the black hole may imply that the infall time scale is much smaller compared to the corresponding viscous time scale [5]. This further indicates that the low angular momentum advective accretion may practically be considered to possess almost constant angular momentum. For such sub-Keplerian hot accretion disk, the assumption of inviscid flow is surely not unjustified from an astrophysical point of view. Such practically inviscid accretion flow may manifest multi-transonic behavior where steady standing shock waves may be developed [6][7][8][9][10][11][12][13][14][15][16][17][18][19].
Far away from the accretor, especially where the flow is subsonic, inflow time scale may become comparable to the viscous time scale. Introduction of viscous transport of angular momentum may be of considerable relevance to probe the accretion dynamics in such region even for low angular momentum accretion flow to assist the inwards drift of the accreting matter. However, the exact prescription for viscosity in an accretion disk is still a matter of much debate [1,20]. It has really been demonstrated that for low angular momentum, sub-Keplerian accretion disks, viscosity can be implemented as a perturbative effect, up to a certain length scale, which involves a first order correction in the standard α-viscosity parameter [3] about the zeroth order inviscid solution, introduced as a 'quasiviscous' disk model in [21,22]. This method reduces the problem of solving a second-order nonlinear differential equation, Navier-Stokes equation [23], to that of an effective first-order differential equation. Such quasi-viscous prescription helps to capture the important properties of the axisymmetric accretion on a comparatively large length scales, without compromising on the fundamentally simple and elegant features of the inviscid flow model [24]. One obtains the stationary integral flow solutions without solving the energy conservation equations -to avoid the complex formulation of the corresponding heating and cooling mechanism.
Existing works on the quasi-viscous flow, however, concentrates on one particular geometrical configuration of flow -axially symmetric accretion in hydrostatic equilibrium along the vertical direction. Such a flow configuration is usually associated with thin disk models. The quasi-viscous formalism, however, requires that the flow possess with ab-initio non-zero value of the radial flow velocity along the equatorial plane. Such a flow configuration with reasonably low angular momentum is better described by a quasi-spherical flow rather than a thin disk structure with its radius dependent flow thickness. A conical wedge shaped geometrical configuration of the rotating flow serves as the appropriate description of the accretion profile under consideration. Such a geometry is be a better alternative of the flow described by a standard thin disk structure.
Results obtained so far in the literature for quasi-viscous disk accretion for flow in vertical equilibrium are to be compared with the calculations performed using a more general theoretical set up capable of studying various features of hot sub-Keplerian accretion disk in all three different flow geometries considered in the literature, i.e., flow with constant thickness, conical wedge shaped flow, and flow in hydrostatic equilibrium along the vertical direction (see, e.g., section 4 of ref. [25]). It is to be noted that the corresponding expressions for the flow thickness in any of the flow configurations in the existing literature are derived from a set of idealised assumptions. In reality, the rigorous derivation of the flow thickness may be accomplished by using the framework of non-LTE radiative transfer (e.g., [26,27]) or, by using the Grad-Shafranov equations for the MHD-related aspects of the flow (e.g., [28][29][30]).
In the present work, we provide a generalized formalism capable of studying the quasi-viscous flow in all possible geometric configuration of the accreting flow structure under the influence of all available post-Newtonian pseudo-Schwarzschild black hole potentials proposed in the literature. We find the multi-critical properties of such flow and demonstrate how to categorise the nature of the corresponding critical points -i.e., how to identify which critical point is saddle type and which critical point is spiral type by following a dynamical systems approach [31]. However, a global understanding of the flow topologies is obtained upon performing a complete numerical investigation of the nonlinear stationary equations describing the dependence of the velocity on the radial distance. While this is a general practice, it is still possible to semi-quantitatively capture the essential behaviour of the transonic solution without resorting to numerical techniques through an alternative means reported in the recent literature on accretion [17,[32][33][34]. Equipped with the mathematical formalism of a general dynamical systems approach, these studies have shown that one can derive a clear analytical conception of some of the global features of the flow by analysing the local properties of the critical points.
Most of the aforementioned works dealing with multi-critical flow essentially address the stationary transonic flow solutions corresponding to the steady state accretion. Transient phenomena are not quite unusual in the large scale astrophysical set up. Also, recently some of the authors [21,22] demonstrated that the introduction of a very small amount of viscosity may give rise to secular instabilities as well. It is thus essential to ensure that the information collected from the stationary integral flow solutions are reliable by establishing that the steady states are stable -at least for the astrophysically relevant time scales. Such a task can be accomplished by perturbing the set of equations describing the steady state quasi-viscous accretion disk and thereby demonstrating that such perturbation does not diverge with time. We perform a linear stability analysis of the stationary quasi-viscous flow solutions to analyse such perturbation.
We observe that for all three flow geometries considered here, a certain length scale can be defined in terms of the αviscosity parameter as well as the density and the temperature of the ambient medium -beyond which the quasi-viscous prescription does not hold good in the sense that the first order correction becomes comparable with its corresponding unperturbed value of the constant specific angular momentum, λ 0 . Numerical estimate for such length scale varies with different geometric configurations of flow. We generalised these studies to all the disk models available in the literature.

Basic Formalism
The accreting material is assumed to possess low angular momentum and is considered to be an advective adiabatic fluid with a weak viscous effect driven by a gravitational field described in terms of a generalised pseudo-Schwarzschild potential, Φ(r).
For an inviscid flow, the system is governed by the equation of continuity and the Euler equation. For a viscous flow one needs to incorporate the angular momentum conservation equation as well to include the effect of viscosity. Due to axisymmetry of the problem and by averaging all the relevant physical quantities over local height, all the governing equations are formulated and solved on the equatorial plane itself.
For a fluid accreting in the form of a thin disk like structure onto some accretor, the continuity equation is obtained as, Where, Σ = ρH is the local surface density of the disk, ρ being the local height averaged volume density, H being half of the local disk height. v is the radial drift velocity of the fluid along the equatorial plane, where r is the radial position on the equatorial plane. The earlier works [21,22] under this quasiviscous approach in accretion disk problem dealt with the disk model where the disk height was determined at each value of r, from the hydrostatic equilibrium condition along the lateral or vertical direction. Apart from that, in the literature, the conical disk model is used where the quasi-spherical flow through a certain solid angle, maintaining azimuthal symmetry and the symmetry with respect to the equatorial plane, is considered. Lastly the simplest and primitive model of a disk with a constant thickness or height has been considered here.
The dependence of H on radial distance r varies in different disk models. These are [35], Here, c s (r) = ∂P/∂ρ = γP /ρ is the local polytropic sound speed in the flow medium with equation of state as P = Kρ γ [36]. Now these can be generalised into a single parametric form, such that, The conservation equation of specific angular momentum in the flow is given by [1], where,ω is the local angular velocity of the flow, G = 2πνΣr 3 ∂ω ∂r is the torque and ν = αc s H in accordance to Shakura & Sunyaev prescription for kinematic viscosity [3]. From the above equation it is clear that even in stationary condition i.e., under the condition ∂ ∂t (r 2ω ) = 0, the effective specific angular momentum r 2ω (=λ eff , say) is not a constant, rather is dependent on r. It would however be independent of r (a constant λ 0 ), for inviscid flow, i.e., when α = 0. Hence, in the weak viscosity limit λ eff may be expressed as a position dependent perturbative correction to λ 0 proportional to α, the Shakura & Sunyaev (1973) dimensionless viscosity parameter as, such that αF << 1.
Therefore, using ρ = (γK) −n c 2n s and n = 1 γ − 1 [36], H and f can be expressed as, Here, f is inherently negative as v is negative for accretion flows. Now, equation (4) may similarly be written as, where f 1 = ρc s H 2 r 3 and . Now, putting the λ eff from equation (5) in equation (10) we get, by assuming ∂ ∂t (λ eff ) ≈ 0 and neglecting all higher orders of α and retaining only the first order, one may introduce the 'quasi-viscous' approximation. Under the same approximation the radial Euler equation becomes, The equation (12) thus becomes the effective dynamical equation in our quasi-viscous model.

Stationary Solutions
For ∂ ∂t = 0, the equations (7) and (12) take the form It is clear that from equation (13), the first integral, f (=ρvrH), is a constant; while −4πf is equal to the mass accretion rate (Ṁ ). However, due to presence of the α dependent term, the second equation (14) cannot provide a closed algebraic form for the first integral which otherwise would represent total specific energy for an inviscid flow.
Also under stationary conditions, the expression of F , using the condition (13) in equation (11), is found to be Hence replacing F by F s in equation (14) and using 1 ρ Again, by using 1 ρ , equation (13) reduces to, Stationary flow equations (16) and (17) can be combined to obtain the velocity gradient as, dv dr =

Critical point analysis
Critical point condition (all corresponding quantities have been labelled with subscript 'c') is given by the condition as N = D = 0. Therefore, we have, Equation (19a) is solved for the critical point radius, r c , for any constant values of γ and α with F sc (c sc ) (using equations (15) and (8)) and with a given choice for the two flow parameters, i.e., λ 0 and the entropy accretion rate (Ṁ); whileṀ is related to the conserved mass accretion rate (Ṁ = −4πf ) through c sc as expressed in [37], using equation (19b). Hence one gets, . Again using equations (15) and (3), it is obtained that, . It may be noted that for the models 'C' and 'H', ǫ = 0 and for the model 'V, nǫ = 1/2 and h = √ K, so that for all three disk models the c sc , H c and F sc and hence location of the critical point (r c ), do not involve K.
Once we determine the r c values, the critical values (or values at fixed points) for flow velocity, sound speed, and local height can subsequently be determined at each evaluated critical radius. Now to find the radial velocity gradient at the critical point one has to take limit of dv/dr as r tending to the critical point and making use of L'Hospital rule one may find a quadratic equation for dv/dr| c as [18], where, with the definitions, Hence the solution we have, Using equation (24) in equation (17) one can obtain the value of dc s dr at the critical point. One may check that the condition, α = 0, implies∆ = 0, which gives back the inviscid results. After knowing r c , v c , c sc , dv/dr| c and dc s /dr| c , one may use standard numerical algorithm such as Runge-Kutta method on the equations (18) and (17) to draw the phase portrait explicitly, but the instead of doing that we follow a more elegant method (see refs. [17,32]) to extract the qualitative features of it using the techniques of dynamical systems described in the following sections.

Categorization of critical points
In what follows we provide a detailed description of the procedure to be followed to understand the nature of the critical points and the topological behaviour of the phase orbits around them. To accomplish such task equation (18) is parametrised by some arbitrary parameter τ (similar approach in invisid cases as followed in [17,22,[32][33][34]38]) and presented as a set of coupled autonomous first-order differential equations as, The system of above equations are then linearised about the fixed points (critical points) as a next step. By applying linear perturbations about the critical point quantities as v 2 = v 2 c + δv 2 , r = r c + δr, c 2 s = c 2 sc + δc 2 s and F s = F sc + δF s , it is possible to derive a set of two autonomous first-order linear differential equations in the δv 2 − δr plane as, using the interrelations, The coefficients have the following forms, After decoupling the equations (27) & (28) with suitable linear transformation and assuming solutions of type ∼ e Ωτ for the transformed dependent variables, the eigenvalues Ω signifying the growth rates may turn out to be, solution of which occur in the form, Here, P , Q and ∆ are defined as, The nature of the possible critical points can be predicted by investigating the form of Ω [31]. If Q < 0, the critical point has to be a saddle. For Q > 0, if P = 0 the point is a centre. The previous two varieties are the only ones found in an inviscid flow. But if P = 0 the critical point must become a node and in addition if ∆ 2 < 0 the point will turn out to be a spiral (node). This analysis therefore delivers the qualitative features of the phase portrait without explicitly obtaining the stationary integral solutions.
For our present work, under the condition α = 0, only saddle and spiral type critical points are obtained. In this paper, we denote the two real eigenvalue solutions by Ω 1,2 (relevant for saddle type critical points) and denote Ω re by P/2 and Ω 2 im by ∆ 2 /4 (relevant for centre or spiral type critical points).

Overall Methodology
The overall plan of this work is as follows. We have already constructed the set of linear energy momentum conservation equations including the presence of the quasi viscous effects, as well as the mass conservation equation. From the modified Euler-like equation and the continuity equation, we found out the equation (19a) for locating the critical points, velocity gradient in eq. (24) and sound speed gradient, conditions for categorising the types of the critical points in a generalised parametrised form in the previous section. In the following sections, based on the generalised conditions, we find the critical points for various geometric configurations of axially symmetric accretion flow. Using the corresponding critical point conditions, the possibility of having multi-transonic flow is discussed, and it is pointed out that stationary shocks may form as a consequence of having multiple critical points. The topological nature of the phase orbits for the stationary transonic integral solutions is discussed. Once we acquire the details of the stationary properties of the flow profile, we get into the stability analysis of such accretion configurations. To accomplish such task, we linear-perturb the stationary flow solutions to check whether such perturbation diverges, and conclude on the nature of the stability of the flow based on such observation.
For inviscid polytropic flow, one integrates the Euler and the continuity equation to obtain two global first integrals of motion -the conserved specific flow energy (the Bernoulli constant) and the mass accretion rate, respectively, see, e.g., [1,15] for further details. The term 'global' signifies that such first integrals remain conserved along the entire streamline in the steady state, i.e., along the entire integral flow solution (phase orbit) characterised by a set of specified initial boundary conditions. For quasi-viscous accretion, however, the modified Euler-like equation (12) becomes an integro-differential equation, and hence can not be integrated analytically. No such complications arise while obtaining the integral solution of the continuity equation. One thus obtains a globally conserved mass accretion rate but the impossibility of obtaining the integral solution of the Euler like equation does not allow to have any specific energy term for the quasi viscous flow which is conserved along the streamline. Such energy like variable can only be defined, as is explained in detail in the subsequent sections.
Presence of non-analytically integrable term in the Euler like equation alters the solution scheme for obtaining the flow profile. For inviscid flow, total specific energy was globally conserved and hence the overall flow behaviour could be parameterized by [E, λ, γ] where E is the globally conserved specific energy of the flow, λ is constant specific angular momentum and γ is the fixed value of the adiabatic index since the atomicity of the gas has been considered to be invariant. For quasiviscous flow, neither E nor λ will be considered as an appropriate boundary condition. Since we do not incorporate any heating or cooling mechanism in our work (viscosity has been introduced as a perturbative effect), γ can still be used as the initial boundary condition for characterising the flow. Instead of [E, λ], we will be using Ṁ , λ 0 to characterise the integral solutions.

Multi-transonic behaviour
Equation (19a) may be solved for a fixed set of Ṁ , λ 0 , γ to obtain the location of the critical point. It is not necessary that the solution will be unique. One may have at most at most three real roots for the equation (19a) -i.e., at most three critical points -may form, all of them being situated outside the black hole event horizon. In subsequent sections, we explicitly demonstrate that the critical points formed for our our quasi-viscous model can be of saddle type or of spiral type. A real physical transonic solution can pass through a saddle type critical point only. For a multi-critical flow, usually two saddle type critical points embrace a spiral type critical point.
If the stationary integral solutions passing through two different saddle type points can be connected by a stationary shock, true multi-transonic flow profiles are formed. All these issues have been demonstrated analytically in subsequent sections.

'Accretion' type and 'Wind' type
For astro-physically relevant set of Ṁ , λ 0 , γ , we obtain the integral solutions passing either through one or through three critical points. Solutions containing three critical points can be categorised into two distinct groups. If we define  [15] )passing through the inner saddle type critical point may launch albeit discontinuously through a standing shock, onto the wind type solution passing through the outer sonic point. Hence one obtains multi-transonic wind solutions with shocks in wind instead of accretion branch. In this case, the conventional ingoing accretion solution will be mono-transonic, in spite of obtaining three critical pints by solving the equation (19a) for a fixed set of initial boundary condition characterised by Ṁ , λ 0 , γ .

Expressions for different flow models:
In Table 1, we display the the expressions for different variables (already introduced in the previous sections) corresponding to three disk geometries. The space gradient of flow velocity for various models is evaluated from equation (18) as, where, β 2 = 2 γ + 1 . The space gradient of sound speed is expressed, from equation (17) by, The real root(s) of the following equations, evaluated from equations (19a), (15), will be the location(s) of the critical point(s) in the corresponding flow models, The flow velocity and the sound speed are related at critical points and the relations can be evaluated from equation (19b), for different flow geometries. For the 'C' and 'H' case, sonic and critical points bear identical meaning, whereas, for 'V' they differ by a constant factor 2/(γ + 1).
The slopes of the transonic flow lines at a critical point (saddle) are evaluated by L'Hospitals rule as obtained in eqs. (21), where the constants A , B and C for different flow models are obtained as, where we write, The nature of critical points for each flow geometry is evaluated by determining the coefficients (A, B, C, D) of the equations (27)- (28) and for various flow geometries those are found to be,

Qualitative understanding of the stationary flow configuration
Results demonstrated in this section are obtained for the Paczynski and Wiita (1980) potential [39], where r is the radial coordinate scaled by the Schwarzschild radius. Similar results may be obtained for any other specific pseudo-Schwarzschild black hole potential.
With the help of the mathematical procedure outlined in section 3.1, the critical points have been obtained by solving equation (19a) for r c (scaled by Schwarzschild radius) in terms of the system parameters, i.e., λ 0 andṀ for a given choice of γ and α. We look for physically allowed single critical points and three critical points (identified as outer-middle-inner critical points) sets with variation in α values for all the three flow geometries.
A subset of the λ 0 -Ṁ parameter space which corresponds to three critical point solutions satisfying the condition, that the local flow energy E 0,outer > E 0,inner as mentioned in section 4, where may support multicritical incoming solution [40]. The other option i.e. E 0,outer < E 0,inner cannot support the multicritical incoming solution in spite of having 3 critical points on the phase plot, rather it would correspond to the single critical point incoming solution through the inner one. Accordingly, in the parameter space the subregions are labelled by A and A 1 respectively. In figure (1), we have shown the single ('O' for outer or 'I' for inner one) and wedge shaped multitransonic i.e. three critical point solutions in the parameter space ofṀ and λ 0 , a feature common to all the three different disk models (V, C & H) considering three values of the α -varying between 0 and 0.1 for γ = 4 3 . However, the span of the [Ṁ, λ 0 ] for multi critical solutions varies between disk models. It covers maximum span for theṀ parameter for the 'H' case and minimum for the 'V' case. Further, for any of the disk geometries, the parameter space corresponding to three critical point solution is classified into A and A 1 region based on the prescription mentioned above. It should be noted that a point on the boundary line between the subregions A and A 1 indicates equal value for the local flow energy for the outer and inner critical points. The multicritical accretion region seems to contract with increasing α values. In particular, with higher value of α, these regions turn to monotransonic 'O' regions in the lower λ 0 side. Next, we have followed the procedure described in section 3.2 to assess the nature of the critical points (to visualise the flow topology) by calculating the Ω values for each critical point obtained by using equations (31a) -(33c). In figures (2 -4), we have presented the variation in eigenvalues against λ 0 for three critical points appearing for the range of λ 0 shown there (in both the regions A and A 1 ) for all the disk models (V, C & H) for 0 ≤ α ≤ 0.1 with γ = 4 3 andṀ=2 × 10 −5 units. As expected, we obtain two equal and opposite real roots for Ω (i.e., Ω 1 and Ω 2 ) signifying saddles as inner and outer critical points and purely imaginary roots (i.e., Ω re = 0 and Ω 2 im < 0) indicating centres as middle critical point with α = 0 (inviscid flow) for all the disk models. However, non-zero values of α produce unequal real values, with opposite signs, of Ω 1 and Ω 2 for the saddles and Ω re = 0 alongwith Ω 2 im < 0 indicating the switching over from the centre type middle critical points to spirals, according to the prescription already mentioned in the section 3.2.  In the usual inviscid formalism, the stationary solutions are flows passing through the saddles or closed paths about a center [17,[32][33][34]. However, the slightest presence of viscous perturbation alters the centers to spirals . Mathematically, α = 0 implies P = 0 and vice versa. For any value of P , we obtain Q < 0 indicating a saddle and Q > 0 with ∆ < 0 implying a centre (for P = 0) or a spiral (for P = 0) . Therefore, all the outer or inner critical points have been found to be saddle and no simple node (with ∆ > 0 and Q > 0) have been obtained for our accreting system. Nodal points, however, cannot be ruled out completely as [38] have shown. Therefore, the resultant picture that emerges as      points, one type located after another.
Once the nature of all the physically relevant critical points have been understood, a complete qualitative picture of the flow solutions passing through these points (if they are saddle points), or in the neighbourhood of these points (if they are centre/ spiral points), can be constructed, along with an impression of the direction that these solutions can have in the phase portrait of the flow. A physical multitransonic flow can result connecting the inner and outer saddle points via formation of standing shocks [16,[40][41][42][43][44]. Under shock formation, the multicritical accretion solution implies that the flow connects infinity to event horizon passing through both the outer and inner critical points via shock formation which allows a discontinuous jump from the supersonic branch of the seperatrix for outer critical point to the subsonic part of the seperatrix of the inner saddle point for a physically allowed accretion process.

Introduction of Time Dependent Perturbation
So far we have analysed the characteristics of the stationary flows in this quasi-viscous scheme without bothering about the sustainability of such a flow. But for such a flow to exist, one has to investigate the stability of such a flow from the time dependent flow equations. In this section thus we attempt to study the linear stability analysis of the quasi-viscous accretion disk about its stationary solutions. Let us introduce small time dependent perturbations (marked by a~sign at the top of the respective variables) on the stationary solution of the flow variables (denoted by suffix 0) at each radial distance r in order to track the time evolution of the perturbation and make predictions about the stability of the background stationary configuration, determined by it. The perturbed quantities appear as, v(r, t) = v 0 (r) +ṽ(r, t), (40a) ρ(r, t) = ρ 0 (r) +ρ(r, t), (40b) f (r, t) = f 0 (r) +f (r, t), (40c) g 1 (r, t) = g 10 (r) +g 1 (r, t), (40d) g 2 (r, t) = g 20 (r) +g 2 (r, t), (40e) F (r, t) = F 0 (r) +F (r, t). (40f) To investigate the stability of the stationary solution, one has to identify whether these perturbing terms grow with time or not. Substituting these expressions given in equation (40) into equation (12) one gets, keeping up to the linear perturbing terms, and the equation (7) gives, Again, Now from (13), Hence from equation (42), From the relation f (ρ, v, r) = ρvrH we get, and hence, -where ∂f /∂v| 0 = ρ 0 rH 0 = f 0 /v 0 and ∂f /∂ρ| 0 = v 0 rH 0 = f 0 /ρ 0 . From the equations (43) and (44) we get Finally after taking another ∂/∂t on the last equation and using equations (43), (44) and (45), one gets an inhomogeneous wave equation, Calculation ofF and ∂F ∂t F can be decomposed into two parts such as F = F s + F N , where F s is already defined in equation (15) and being the additional contribution under no-stationary condition. Now after perturbation, the perturbed quantity F s = F s0 +F s0 and F N =F N ; because at steady state F N 0 = 0 From equation (15) we getF From equation (44) it follows that Substituting equations (50) and (45) in equation (49) From equation (47) we can write Now combining equations (15), (51) and (52),

Solution of the Wave Equation
Now choosing a trial wave solution,f one gets may explore the behavior of the spatial part. Hence substituting the above in equation (46), the equation becomes, The solution of this equation may be of two types. One is of standing wave type, wheref takes the form A ω (r, t) exp(−iω r t) where A ω (r, t) is real and ω may be complex with a real part ω r . The other is of travelling wave type wheref takes the form A ω (r, t) exp [i(s r (r) − ωt)], both A ω (r, t) and ω being real though s r (r) may be the imaginary part of a spatially dependent complex exponent is(r) in the solution.

Standing Wave Analysis
The stationary transonic solution connects the outer boundary at infinity to inner boundary at event horizon pass through the critical point (saddle type). Complying the respective boundary conditions, there may always be a rapid time variation or fluctuations at different points, which may be captured by a standing wave solution.
Upon integrating the above equation (55) so that the surface integral terms, evaluated at the boundaries, vanish we get finally something like such that where Evidently the exp(αξ(r)) enters into A ω (r, t) and the solution diverges making the situation unstable if ξ(r) diverges to positive infinity at any particular length scale.

Travelling Wave Analysis
Disturbance may propagate too, which may be of travelling wave type. Rearrangement of equation (55) gives and Taking the spatial part of the perturbation as where the function s(r) itself is represented as a power series of the form The integral term in equation (60) can, through some suitable algebraic substitutions, be recast as where S itself is given by another power series as To a leading order Substituting for g ω in equation (60) and collecting the same powers of ω gives us coefficients of ω 2 ,ω and ω 0 Making coefficient of which leads to, Similarly making coefficient of which gives us k 0 = − 1 2 ln f 0 g 20 (r)c s0 √ 1 + εg 10 (ρ) ± αλ 2 0 c 2 s0 f 0 With large enough value of ω, characteristic to any noise as perturbation, these leading order terms determine the nature of the solution.

Stability under perturbations
For standing wave analysis the behavior of ξ(r) determines how rapidly the perturbation grows. Now after putting the specific forms of g 1 , g 2 and ǫ, for the standing wave solutions one gets for H model, ξ(r) ∼ r −2 while for C, ξ(r) ∼ r 1 and for V, ξ(r) ∼ r 5/2 ; which readily implies that at large radial distance for C & V the perturbation rapidly grows with time making the disk unstable but for H it diminishes with radial distance making the disk practically stable for large amount of time. Thus in the case of standing wave solution the asymptotic behaviour of ξ(r) shows the disk may become highly unstable if it spans large enough radial distances; though there is a question of time scale that ought to be within astrophysically realisable timescale. If ξ(r) diverges at some r, the disk eventually gets destabilized at that distance for whatever small value of α under the present quasi-viscous scheme. But actually this disk analysis is valid only under the condition αF << 1. As, it is evident that ξ(r) ∼ḡ 3 (r)/r 2 cannot diverge within this limit, irrespective of the flow geometry (the minus sign comes because v and f , both are negative for accretion).
In the case of traveling wave too, it is easy to see that k −1 ∼ r and for α = 0, k 0 ∼ ln(r). But if α = 0, k 0 ∼ αḡ(r) which may make is(r) diverge for r → ∞. But again αF << 1 puts restriction on that. The rest of the terms, with higher powers of 1/ωwhere ω is generally large enough (≫ 1), are not expected to make the whole thing diverge. So here too within the region where the present scheme of viscosity remains valid, there is no imminent threat to the stability of the disk.
Here this analysis is valid only under the condition when the disk is formed within αF s << 1 and only under this condition the disk is stable. As, the validity of the model is restricted within the region where, g 2 (r) << f −2ρ 2ǫ+1 h 2 c s which may be estimated by plugging in corresponding values of the quantities in the right hand side.

Concluding Remarks
In the Sakura-Sunayev model, viscosity is introduced in the entire set of equations governing the disk dynamics to understand the flow profile. But, it takes a lot of mathematical difficulties to solve the equations explicitly. On the other hand in this present quasi-viscous approach, viscosity is introduced as a perturbative correction to the constant specific angular momentum in the radial Euler equation for inviscid flow, which saves a lot of mathematical complexities. Our quasi-viscous approach, however, as shown is valid within a limited region of the disk depending upon the value of viscosity parameter, α. But, within this region we compared the qualitative nature of the transonic accretion flow topologies over three different disk geometries. The formulation of the problem and the stability analysis was done in a single parametrized form for all the three geometries. This analysis too has been done here without explicit numerical integration by using the methodology of dynamical system analysis. The stationary solutions conform with the earlier studies using Sakura-Sunayev scheme [40] and quasi viscous scheme [22]-both of which dealt with vertical hydrostatic equilibrium disk model only.
In some earlier works [21,22] the respective authors identified the existence of a secular instability under linear perturbation about the stationary solutions using the time dependent equations raising doubt about the sustainability of such a stationary flow. Here too, we identified the same sort of secular instability -at least in another disk geometry, along with the previously identified oneasymptotically at a large length scale. But, we have shown that such divergence in the time evolution of the perturbation cannot take place within the length scale in which the model itself is applicable. That conclusion is valid both in stationary and travelling wave type perturbations on the stationary flow background. Hence we conclude that the stationary flow sustains at least for a considerable period of time for sufficiently low value of viscosity parameter (0 < α ≤ 0.1) for which the present approach stands.