Oscillating shocks in the transonic viscous, variable Γ accretion flows around black holes

We investigate the time evolution of the transonic-viscous accretion flow around a non-rotating black hole. The input parameters used for the simulation are obtained from semi-analytical solutions. This code is based on the TVD routine and correctly handles the angular momentum transport due to viscosity. The thermodynamic properties of the flow are described by a variable adiabatic index equation of state. We regenerate the inviscid and viscous steady-state solutions, including shocks, using the simulation code and compare them with the semi-analytical solutions. The angular momentum piles up across a shock due to shock-jump conditions and viscous transport of angular momentum. This will push the shock-front outward and can result in shock oscillation or a complete destabilization of shock. We study how shocks behave in the presence of viscosity. As the viscosity parameter ( α ) crosses a critical value, the previously steady shock becomes time-dependent, eventually leading to oscillations. The value of this critical viscosity depends on the injection angular momentum ( λ ou ) and the specific energy ( ϵ ). We estimated the posteriori bremsstrahlung and synchrotron cooling, and the net radiative output also oscillates with the frequency of the shock. We also study the variation of frequency, amplitude, and mean position of oscillation with α . Considering a black hole with a mass of 10 M ⊙ , we observed that the power spectrum exhibits a prominent peak at the fundamental frequency of a few to about tens of Hz, accompanied by multiple harmonics. This characteristic is frequently observed in numerous accreting black hole candidates.


INTRODUCTION
Theory of accretion onto the compact objects started its journey with the seminal papers of (Hoyle & Lyttleton 1939) and (Bondi 1952).However, accretion physics gained popularity after the discovery of quasars and X-ray binaries.In 1964, Salpeter computed the luminosity from accreting matter onto black holes by using the Bondi accretion model, but the luminosity obtained was orders of magnitude lower.Bondi flow is the radial inflow model; therefore, the gravitational energy released increases the infall velocity, which makes the infall timescale of accreting matter much smaller than the radiation cooling timescale.Hence, the radiated luminosity is low.In order to increase the infall time scale compared to the radiative cooling timescale, matter must possess angular momentum.From the black hole inner boundary condition, the matter is most likely to possess sub-Keplerian angular momentum closer to the horizon but may possess higher angular momen-tum at a larger radius.Therefore, a general accretion disc requires viscosity to effectively remove angular momentum from the system.The foundational work on a viscous accretion disc was presented by Shakura & Sunyaev (1973).This model assumed a Keplerian angular momentum distribution, along with a geometrically thin and optically thick accretion disc.This disc is known as the Keplerian disc (KD).The spectra computed from the KD by summing up the black body emission from each of the annuli (each with different temperatures) matched the thermal part of the spectrum.As a result, KD is widely used to explain the observed optical continuum spectra of luminous active galactic nuclei (AGN) or the soft X-ray continuum spectra of the X-ray binaries.However, the KD did not adequately consider the pressure and advection terms, and it failed to satisfy the inner boundary condition near a black hole, except for the ad-hoc termination of the disc at r ≤ rms (where rms represents the marginally stable orbit or the Innermost Stable Circular Orbit, ISCO).Moreover, KD could not reproduce the non-thermal powerlaw part of the spectrum associated with objects harbouring black holes.Therefore, the search for a general accretion model that can explain the high-energy non-thermal radia-tion was investigated by various groups.The next important model is called the thick disc model, also sometimes called the 'Polish doughnut' model, was developed by Fishbone & Moncrief (1976); Paczynski (1980); Paczyńsky & Wiita (1980); Abramowicz et al. (1980) etc.In this model, the pressure gradient force along the radial as well as the transverse direction was considered, although the advection term was neglected.This accretion model was rotation-dominated and hot; therefore, it was puffed up in the form of a torus.Hence the name 'thick' accretion disc.Thick discs were shown to be susceptible to Papaloizou-Pringle (PP) instability (Papaloizou & Pringle 1984).However, it was also shown that for various configurations, thick discs are PP stable.Simultaneous with the development of these models, there were many efforts to develop models with advection terms present in the equations of motion.In 1980 Liang & Thompson studied transonic, rotating accretion solutions in Schwarzschild metric.This paper invoked a Bondi-type transonic solution but in the presence of rotation.They also showed that, unlike the Bondi flow, a rotating transonic solution may admit multiple sonic points.Subsequent studies (Fukue 1987;Chakrabarti 1989) have further demonstrated that two standing shocks can exist between the two physical, X-type critical points.The inner shock jumped onto an accelerating part of the subsonic branch while the outer shock jumped onto the decelerating part of the same subsonic branch.It was shown that the outer shock position out of the two available is stable, using analytical methods (Nakayama 1992) and also by using smooth particle hydrodynamic simulations (Chakrabarti & Molteni 1993).The instability of the inner shock was also studied in details in few other papers (Nakayama 1993(Nakayama , 1994;;Nobuta & Hanawa 1994;Nakayama 1996).Fukue (2019a,b) has recently investigated various types of shock in accretion disc.Fukue (2021) calculated shock conditions in the accretion disc by considering αp viscosity prescription.Chakrabarti & Molteni (1995) proposed a new form of viscous stress tensor that conserved the angular momentum across any shock jump.Based on this prescription, a number of studies on viscous flows were undertaken (Chakrabarti 1996;Gu & Lu 2004;Chakrabarti & Das 2004;Chattopadhyay & Das 2007;Das & Chattopadhyay 2008).However, many papers both in the analytical regime and numerical simulations, were also written using traditional viscosity prescription, where the viscous tensor was assumed to be proportional to the shear tensor (Lanzafame et al. 1998;Lu et al. 1999;Becker et al. 2008;Das et al. 2009;Lee et al. 2011;Kumar & Chattopadhyay 2013;Das et al. 2014;Lee et al. 2016).One of the popular versions of accretion flows in the advective regime is called ADAF or Advection Dominated Accretion Flow (Ichimaru 1977;Narayan & Yi 1994), which is a radiatively inefficient flow in which most of the dissipated gravitational energy is advected to the central black hole.The initial solution for ADAF was self-similar and completely subsonic.However, subsequent global solutions of ADAF, considering the influence of strong gravity, demonstrated that the flow might be self-similar far from the central object but has to be transonic near the horizon (Chen et al. 1997).Lu et al. (1999) actually made a thorough classification of viscous, transonic solutions in the advective regime and showed that ADAF solutions (monotonic, viscous accretion solutions with single sonic point closer to the horizon) are a part of the general advective solutions.
One of the major propositions for the advective tran-sonic disc is the identification of the post-shock region to be the major contributor for inverse-Comptonized hot photons (Chakrabarti & Titarchuk 1995), in other words, the illusive Compton cloud.Even more interesting was the fact that Molteni et al. (1996a) showed that the post-shock disc could eject a part of the accreting matter as bipolar jets.Moreover, Molteni et al. (1996b) showed resonant oscillation of the post-shock region, which might explain the lowfrequency quasi-periodic oscillation (QPO).The resonant oscillation was achieved by the enhanced radiative cooling rates in the post-shock region.On the other hand, a number of simulations showed that viscosity could also trigger shock oscillation (Lanzafame et al. 1998;Das et al. 2014;Lee et al. 2011;Giri et al. 2015;Lee et al. 2016).These simulations showed that beyond a critical viscosity, the post-shock disc may oscillate which would produce an oscillation with similar frequency of high energy photons emitted by the post-shock disc.Impressive as these results might be, all of these simulations were launched with super-sonic injection, and therefore, those were not proper outer boundary conditions.It was pointed out in steady-state analytical works that shock location may or may not shift outwards with the increasing viscosity parameter.It all depends if the average angular momentum of the post-shock disc (PSD) is higher or lower than the angular momentum at the outer boundary (Kumar & Chattopadhyay 2013).Moreover, these simulations used a fixed adiabatic index (Γ) equation of state (EoS) of the gas that constitutes the accretion disc.It may be noted that an accretion disc may extend from near the black hole horizon to a few thousand Schwarzschild radii, and the temperature of the flow over this entire spatial range may vary for two to three orders of magnitude.Over such a range of temperature distribution, the fixed Γ EoS is not tenable since, in principle, Γ is a function of temperature.It was also shown that Γ is a function of temperature as well as the composition of the gas (Chattopadhyay 2008;Chattopadhyay & Ryu 2009).
In the present paper, we study time dependant transonic, viscous accretion flow onto a non-rotating black hole, using the variable Γ EoS proposed by Chattopadhyay (2008); Chattopadhyay & Ryu (2009) and its acronym is CR EoS.The angular momentum transport will depend on the local temperature, viscosity parameter, and the Mach number.It may be noted that supersonic flow is less susceptible to viscous momentum transfer compared to a subsonic flow.Since a CR EoS computes the temperature correctly, therefore, studying the effect of viscosity on flows with CR EoS is very important.Further, if one wants to study only shocked solutions, then supersonic injection may make partial sense.However, if the interest encompasses all kinds of accretion solutions including Bondi type, shocked as well as ADAF type solution then one has to study solutions with sub-sonic injection.Even for shocked accretion flow, supersonic injections produce insignificant angular momentum transport in the pre-shock part of the accretion disc, while there is significant angular momentum transport in the post-shock disc, which would render an angular momentum imbalance in the computation setup of the disc.In the case of the subsonic outer boundary condition, the angular momentum transport in the pre-shock disc (precisely, in the subsonic region beyond the outer sonic point) is more than the supersonic part and therefore, the average angular momentum distribution on either side of the shock-front might be comparable and therefore may stabilise the shock oscillation.CR EoS has been implemented in simulation codes studying various problems (Chattopadhya et al. 2013;Joshi et al. 2022a,b;Joshi & Chattopadhyay 2023).However, since angular momentum is so important for accretion discs, we have replaced the azimuthal velocity by the specific angular momentum as one of the primitive variables and used the corresponding eigenstructure to write the TVD code.Using this new code, we would first like to test the code in the steady state inviscid and viscous rotating accretion flows by comparing the simulation results with the semianalytical solutions.In the process, we would like to check all possible accretion solutions in the advective regime, with subsonic outer boundary conditions.Our main intention is to study time-dependent shocked accretion solutions for a wide range of different conditions.We would like to see how viscosity may affect the accretion shock in the time-dependent regime.If the shock is oscillating, we would like to quantify the range of viscosity parameter which would admit shock oscillation.And if the shock oscillation due to viscosity can be associated with QPOs.
We have organised our paper as follows.We first describe the time-dependent conserved form of the equations of motion.Then, we present the equation of state (EoS), steady state semi-analytical viscous accretion flows solution in subsections 2.1 and 2.2.In section 2.3 we discuss our simulation code.In section 3 we present the verification of our code considering both the inviscid as well as viscous flows.We present our results from the simulations in section 4. Finally, we discuss our results in section 5 in terms of the observed properties of the accreting black holes.

ASSUMPTIONS AND EQUATIONS
This study focuses on examining the properties of viscous accretion flow around a non-rotating black hole, where the pseudo-Newtonian potential (Paczyńsky & Wiita 1980) is utilized to account for the strong gravitational field.We consider that the flow is conical or wedge flow.The geometry of the accretion flows is shown in Fig. 1, where θ0 is the co-latitude of the surface of the accretion disc.The central black spot represents the black hole (BH).The fluid conserved variables or the state variables are represented as qs, and the primitive variables are given as ws, Here, ρ is the fluid rest mass density, Mr = ρvr, M θ = ρv θ are the momentum density in the r and θ directions.M λ = ρλ is the angular momentum density, E = ρv 2 /2 + e represents the total energy density, which includes the internal energy and the kinetic energy.Here, p is the pressure, e is the energy density of the gas, vr, v θ and v ϕ are velocity components along r, θ and ϕ while λ = rv ϕ is the specific angular momentum.In addition The equations of motion in the conserved form for 1D spherical coordinate system are given as Where, the fluxes corresponding to the qs are given as and the source terms of the equations of motion are given as The set of equations represented by Eq. (2) consists of the continuity, three momentum equations, and the energy equation.The second row of S contains the gravitational force density, which is presented as GMBHρ/(r − rg) 2 , obtained from the Paczńsky-Wiita potential (Paczyńsky & Wiita 1980), which is given as Φ = −GMBH/(r − rg), where G is the gravitational constant and MBH is the black hole mass and the Schwarzschild radius is rg = 2GMBH/c 2 .Moreover, ρv 2 ϕ /r is the force due to the rotation.Additionally, the angular momentum equation includes a source term arising from viscous dissipation, denoted as S λ .The fifth row of S represents the work done due to the gravity, which is proportional to GMBHρvr/(r − rg) 2 , and the energy source term is SE.
Here, S λ and SE are given as, here, fc is cooling parameter, where 0 ≤ fc < 1. W rϕ is the r−ϕ component of viscous stress tensor.The local heat gained and lost by the flow are given by Q + = W 2 rϕ /η and Q − .Here We assume that the fc = 0 in our study, until mentioned otherwise.The viscous stress is given by, where, Ω is the angular velocity, ηv = ρν is the dynamic viscosity coefficient, ν = (αp)/(ρΩ k ) is the kinematic viscosity, α is the Shakura-Sunyaev viscosity parameter.Ω k is the local Keplerian angular velocity and is given by, An additional closure relation for the thermodynamic variables e, p, ρ is required to solve the equations of motion (Eq. 2) and is given by the CR EoS.

Equation of state
Taub (1948)  where, f is given as The fluid mass density, ρ in the above equation can be expressed as ρ = nimi = neme(2 − ξ + ξ/η), where η = me/mp, ξ = np/ne, and np, ne, mp, and me presented the proton number density, the electron number density, the proton rest mass, and electron rest mass, respectively.In this paper, we assume that the accretion flow consists solely of protons and electrons (e − − p + ).As a result, the value of ξ is 1.Additionally, Θ = p/ρc 2 serves as a measure of temperature and τ = 2 − ξ + ξ/η.The expression for specific enthalpy is also provided as, The polytropic index N is given as, And the adiabatic index is The sound speed is cs = Γ Θ.It can be inferred from equations ( 11) and ( 12) that the polytropic index and, consequently, the adiabatic index are dependent on Θ and ξ.Therefore, we don't have to provide them as independent variables.By examining equation ( 12), it becomes apparent that when the temperature is high (Θ ≫ 1), Γ approaches 4/3, and when the temperature is low (Θ ≪ 1), Γ approaches 5/3.

Steady state equations of motion
In the steady state ∂/∂t ≡ 0.Moreover, we assume the velocity vector to have the following components v ≡ (vr, 0, v ϕ ).
The equations of motion are solved under these conditions.The mass accretion rate is found by integrating the continuity equation and is given as where G = 4π cos θ0, θ0 is the co-latitude (see Fig. 1).Since accretion is a process where the radial velocity is directed inward so the accretion rate is negative (see, Lu et al. 1999).
Mass accretion rate Ṁ is a constant of motion.Another constant of motion of the flow is obtained by integrating the momentum balance equation, and the energy equation in the steady state is given by (see Molteni et al. 2001;Gu & Lu 2004;Becker et al. 2008;Kumar & Chattopadhyay 2013, 2014), It is the specific energy of the flow and is sometimes called as the generalized Bernoulli parameter.It remains constant throughout the flow even in the presence of viscous dissipation, and λ0 is the specific angular momentum at which the shear is zero presently, it is at the horizon and λ is the local specific angular momentum.The derivation of Eq. 14 is straightforward.The second and the third terms encompass centrifugal term and the viscous transport of angular momentum of the flow.For inviscid flow, α = 0 and λ = λ0, then we retrieve the expression of the Bernoulli parameter in the inviscid limit.Using the energy balance equation and the EoS equation given by equation 8, we can obtain the temperature gradient as Manipulating all the flow equations along with the EoS and using the geometric units we obtained, Where, The azimuthal component of the momentum equation in the steady state on integration is given by At large distances from the BH horizon, the matter is subsonic, while near the horizon, it is supersonic, therefore making it inherently transonic.As a result of this transonic behaviour, there are specific locations where D → 0 and Num → 0 in equation ( 16).Such a location is commonly referred to as the sonic point or critical point.Depending on flow parameters (ϵ, λ0), we can have single or multiple sonic points.Again for multiple sonic points, for the particular parameters, flows can have a shock transition in between the two sonic points.The initial point of integration is a point asymptotically close to the horizon (rin = 1.001 rg).At rin, we supply a guess value of Θin and compute the value of vrin from given values of constants of motion ϵ and λ0.We then integrate equations 16, 15 & 19 simultaneously to obtained the values of vr, Θ and λ.We iterate over the values of Θin for a given set of ϵ and λ0 to check for the solutions that pass through the sonic points (Becker et al. 2008;Kumar & Chattopadhyay 2014).The solutions may be smooth or may pass through two sonic points connected by a shock.The shock condition in an accretion disc where the viscous tensor is proportional to the shear tensor where given by several authors Becker et al. (2008); Kumar & Chattopadhyay (2013, 2014) is essentially conservation of mass flux, momentum flux, angular momentum flux, and the energy flux, The suffix '±' represents the post-shock and pre-shock region of the flow, respectively.On top of this, we assume shear less shock front (dΩ/dr|+ = dΩ/dr|−), which is a reasonable assumption for steady shock.The angular momentum jump condition across the shock is given by Here k = −vr+(λ+ − λ0)/Θ+.Equation 23 can be restructured as, Since, across a shock vr+ < vr− and Θ+ > Θ−, so λ+ > λ−.Jump condition for the velocity and temperature are given by Where Where The shock condition assuming hydrostatic equilibrium in the vertical direction has been computed and presented in Kumar & Chattopadhyay (2013, 2014).

Simulation code and numerical method
Our simulation code is implemented using the total variation diminishing(TVD) scheme, which was first introduced by (Harten 1983).This is a second-order accurate finite difference scheme and is Eulerian in nature.The TVD scheme has proven to be a robust approach to simulating astrophysical phenomena, as it can efficiently capture shocks.(2022a,b).Our code compute the spatial and temporal evolution of conserved quantities ρ, ρvr, ρv θ , ρλ and E using a Roe-type Riemann solver.Here vr and v θ represent the velocity of the flow in r and θ direction, respectively.In most Eulerian codes without viscosity, the azimuthal momentum density (ρv ϕ ) is treated as a conserved quantity instead of the angular momentum density (ρλ).However, in our code, we use ρλ instead of ρv ϕ to accurately compute the angular momentum such that it remains constant when there is no viscous dissipation.In a spherical coordinate system, the conserved variable(q) and the primitive variable(w) are shown in equation 1. Multiplying Eq.( 2) by r 2 , we obtain, Here q = r 2 q, and accordingly the fluxes are redefined as F = r 2 F r and source term as S = r 2 S. The eigenvalues of the above equation are given by, The right eigenvectors (R) obtained are Left eigenvectors (L) corresponding to right eigenvectors are obtained from the orthogonality condition, Where Ri and Lj are the component of the right and left eigenvectors respectively.The updating of state vector qn to qn+1 follows the procedure of (Harten 1983), ), ( 32) Here k = 1 − 5 goes for the five characteristic modes.The parameters ξ k implicitly control numerical viscosity, and they are defined for 0 ≤ ξ k ≤ 0.5.The updated fluid state variables q are obtained by incorporating the source terms with qs The upgraded primitive variables (w) are computed from the q.Recalling the Eq. 1, The relations between the first four components of state and primitive variables are easy The relation between w5 and q5 is not straightforward.From Eq. 1 we know E = ρv 2 /2 + e and combining the expression of EoS (Eq.8), we obtain, Here C is known at each of the computations cells from the updated values of E, ρ, and v 2 .Combining Eq.( 9) with Eq.( 40), we obtain the cubic equation of Θ as, where, Eq.( 41) has an analytic solution.We follow the standard methods described by Abramowitz & Stegun (1970) to solve the cubic equation and found that the CR EoS admits only one unique root.Defining b1 = 9K1/27; b2 = 12K2/27; and b3 = 16(1 − C)/(27ητ 2 ), then The number of roots of Eq. ( 41) will depend on the value of Ξ.In the case of CR EoS, we have Ξ < 0, and Q < 0, which means the cubic equations has three real but unequal roots.The physical root of Eq.( 41) is, Where, cos( θ) = T √ −Q 3 .Once we get Θ, we retrieve the pressure from the definition of Θ i.e. p = ρΘ.
The computational box is divided into equal radial grids.The incoming matter enters the computational box through the outer boundary.We set the density of the incoming gas, denoted as ρou, to a value of 1.This choice is made under the assumption of no cooling or self-gravity, allowing for the density to be scaled out of the calculations.At the outer boundary, we provide the radial velocity vr and the temperature Θ of the flow.These values are obtained from analytical results.To simulate the behaviour near the black hole horizon, we have introduced an absorbing inner boundary situated at 1.5rg.Within this boundary, all material is fully absorbed into the black hole.The background material in the simulation is characterized by a density of ρ bg = 10 −6 .Pressure at the outer boundary and the background material has been calculated from Θ.

CODE VERIFICATION
Two tests have been presented to demonstrate that the code can handle the transonic flows in the vicinity of black holes.We test the code with the steady-state analytic accretion solution around a black hole in spherical geometry, both for the inviscid flow and flows in the presence of viscosity.

Inviscid flow
In this section, we compare the theoretically obtained steadystate solutions with the results of the simulation code in the absence of viscosity.In the absence of viscosity, angular momentum is constant throughout the flow.For strong gravity, if the flow angular momentum λ < λms then accretion is possible.Here λms is the specific angular momentum of a particle rotating in ISCO or the marginally stable orbit.For the inviscid case, we choose λ < λms.The classification  of the ϵ − λ parameter space was first presented by Fukue (1987); Chakrabarti (1989).If the flow angular momentum is low, then the solution may pass through a single sonic point located far from the central object, and the solution will be Bondi type.However, for low λ but higher values of ϵ, the sonic points are located closer to the horizon.For flows with higher or intermediate λ, the flow may harbour multiple sonic points.With changing angular momentum, the spherical symmetry imposed by gravity is modified due to the presence of rotation, which causes the formation of multiple sonic points.If the rotation is strong enough, then it can trigger shock transitions in accretion flows.In Fig. 2 the domain which admits multiple sonic points or multiple critical points (abbreviated as MCP) in the ϵ − λ parameter space is marked with the red dashed curve.While the shock parameter space is marked by the black solid curve.On this parameter space, we marked four regions in the ϵ-λ parameter space as 'a' (ϵ, λ ≡ 1.001, 1.55), 'b' (ϵ, λ ≡ 1.001, 1.72), 'c' (ϵ, λ ≡ 1.001, 1.8), and 'd' (ϵ, λ ≡ 1.005, 1.8).In Fig. 3a-d This solution corresponds to higher λ and admits multiple sonic points.The flow becomes supersonic through the outer sonic point.The supersonic flow then 'feels' the rotational barrier due to higher values of λ.This barrier along with the thermal pressure, acts as a barrier to the supersonic flow and discontinuously slows down the flow in the form of a shock or r sh = 31.83(vertical jump) and jumps to the subsonic branch.The flow again becomes supersonic and dives into the BH after passing through the inner sonic point (rci = 2.617).Dotted curves present all possible solutions through the two physical sonic points available, but the physically correct accretion solution is the one traced by the simulation (red open circles).In simulations, the shock is captured in 6 cells.The analytically determined outer sonic point is at rco = 123.306,while in the simulations it is at 123.323, a fractional difference of only 0.0138%.The analytical inner sonic point is at rci = 2.617 while numerically it is at 2.667, a mere 1.9% in fractional change.The error is more in locating the inner sonic point however, it is indeed quite low.In vrou = −1.336× 10 −3 , Θou = 2.099 × 10 −3 , respectively.In both these cases, inner sonic points obtained with our simulation code matched the analytical result at around 2% accuracy.It may be noted that in the case of Fig. 3(c), we did not inject with the flow variables of α type topology.The main reason is that the α type topology is not global, but even if we choose to inject with flow values corresponding to the α type branch of Fig. 3(c), the injected matter jumps to the global solution (the black, solid curve) and closely follows it.

Same viscosity parameter different ϵ and λ0
In this section, we compare the theoretical steady-state, viscous solutions with the results from the simulation code and see how our code manages to tackle viscous dissipation.Viscosity redistributes and transports the angular momentum outwards and, in the process, heats up the flow.So, a code should be able to capture the angular momentum distribution accurately in order to understand the accretion disc dynamics properly.It may be noted that in the presence of viscosity, the generalized Bernoulli parameter ϵ (Eq.14) is a constant of motion, while λ0 is the constant of integration, therefore the parameter space for viscous accretion is marked in the ϵ − λ0.In Fig. 4, we mark the domain that admits multiple sonic points or the MCP region with a red dashed curve for the viscosity parameter α = 0.01.The shock parameter space is marked by the black solid curve.Due to the reduction of the angular momentum along the flow in the presence of viscosity, it causes a shift in the MCP and shock parameter space toward the lower end of the λ0 scale, resulting in a narrower range.The cusp of the MCP region for inviscid flow is around ϵ ∼ 1.009 (see, Fig. 2), while in Fig. 4, the cusp is at ϵ = 1.0052.Similarly, the shock region is also smaller.We marked four points in the ϵ-λ0 parameter space as a (ϵ, λ0 ≡ 1.0006, 1.The flow is for low angular momentum and due to viscosity, angular momentum transport will be such that the solution passes through an outer type single sonic point rco = 187.15(indicated by the thick blue dot).In Fig. 5b, we compare the analytical (black solid) with the simulation result (red open-circle) corresponding to the point 'b' of Fig. 4. The injection parameters are vrou = −2.091× 10 −2 , Θou = 6.324 × 10 −4 , λou = 1.7212.This solution corresponds to higher λ0 and admits multiple sonic points.The flow becomes supersonic through the outer sonic point (rco = 184.89).The supersonic flow then 'feels' the rotational barrier due to higher values of λ.This, in addition to the thermal pressure, acts as a barrier to the supersonic flow and discontinuously slows down the flow in the form of a stationary shock at the radial distance r sh = 21.18(vertical jump) and jumps to the subsonic branch.The flow again becomes supersonic and dives into the BH after passing through the inner sonic point (rci = 2.71).The dotted curve through rco presents the part of the supersonic branch through which the matter would have accreted if there was no shock.Again, the shock is resolved in 6 cells.Accretion shocks, if driven by the r − ϕ component of the viscous stress tensor, are thin shocks.However, if the viscosity originates from the θ−ϕ component or r−θ component of the stress tensor, then the shock front might be of a finite width (Zeldovich & Raizer 1966).Since r − ϕ component of the viscous stress tensor is dominant, the type of shock we study in this work is thin.In Figs. 5 c & d we compare analytical and simulation results for flow parameters corresponding to points 'c' and 'd' of Fig. 4. Since the flow in these two panels is for even higher angular momentum, the accretion flow velocity remains largely subsonic and eventually becomes supersonic to dive into the BH through the inner type sonic point rci.The injection parameters are, vrou = −1.419× 10 −2 , Θou = 6.709 × 10 −4 , λou = 1.794, and vrou = −2.317× 10 −3 , Θou = 1.598 × 10 −3 , λou = 2.431, respectively.
As we are examining viscous accretion flows, it is necessary to ensure that our code can handle the transport of angular momentum.In Fig. 6a-d, we have compared the theoretical λ distribution (solid, black) with the simulation (red, open circle) corresponding to the solutions presented in Fig. 5a-d.At the shock, the λ jumps to a higher value in the postshock flow, as is clear from Eq. 23.Here also, wherever |M | dips, which implies hotter and slower flow, Most interestingly, in Fig. 6b, the simulations show that both the post-shock branch through rci and the pre-shock branch through rco in the shocked solutions are solutions with the same λ0.In Kumar & Chattopadhyay (2013), we proposed that λ0 for the branches would be the same and obtained the solution, in this paper, we show that even numerical simulations confirm our proposal.The solution represented by Figure 5c is not a shock solution, it is a smooth but not a monotonically increasing or decreasing function.The angular momentum piles up in regions where the magnitude of the Mach number is low.It implies that, if the flow is slow (|M | < 1) and hot (Θ high), viscous transport of λ is significant (Fig. 6c).

Same ϵ and λ0 but different α
We would like to study how viscosity affects the flow variables for given constants of motion.Kumar & Chattopad-hyay (2013) analytically showed that by keeping ϵ and λ0 same if we increase viscosity parameter (α) solution topology changes.This amounts to studying accretion flows with the same inner boundary condition.Kumar & Chattopadhyay showed in their theoretical study that an accretion flow with α = 0, which is flowing through an outer sonic point and does not have shock, may undergo shock transition with the increase of α.On the other hand, solutions with shocks in inviscid limit will become smooth solutions with the increase of viscosity.All these depend on the angular momentum transport due to viscosity.In this paper, we show, through simulations, one case where a solution has no shock in the inviscid regime but generates shock as viscosity is increased, and with further increase in the viscosity parameter, the solution becomes shock-free, all the time keeping ϵ and λ0 same.Since ϵ and λ0 is same for all the solutions, but α is different for each of them.Therefore, the outer boundary conditions of all three cases presented in Figs.7 are different.
We compare analytical accretion solutions for flow parameters ϵ = 1.005 and λ0 = 1.67 with numerical simulation solutions.For the numerical simulation, we injected flow variables obtained from the analytical solution.The rou = 350 and the code had 3072 uniform grids.The injection parameters are mentioned in the caption of Fig. 7 and from the panels left to right, the injection speed, temperature, and angular momentum increase with increasing α.Black solid lines represent analytic accretion solutions, whereas red open circles show the simulated solutions.In Fig. 7a & d, we plot the Mach number M and the specific angular momentum λ as a function of r for inviscid (α = 0) accretion flow.In Fig. 7b & e we plot the accretion solution with α = 0.01, and the accretion solution for α = 0.03 are plotted in Fig. 7c & f.Values of injection parameters, i.e., velocity, temperature, and angular momentum, are taken from the analytic solution.In the lower panels (d, e, f), the λ distribution is shown, and λ is strictly constant for inviscid flow or α = 0.0, and for the viscous flow, the angular momentum is no more constant but the simulation code quite accurately computes it.There is a shock in the accretion disc at 18.8rg for α = 0.01 from the analytic result.But from the simulation, we get a shock at around 17.37rg.For this set of flow constants ϵ, λ0, we checked that analytically, accretion flow will harbour a steady shock if 0.008 ≤ α ≤ 0.0279.But in simulation, we find shock in the range 0.009 ≤ α ≤ 0.0285.It may be noted that we only compared the steady semi-analytical accretion solution for the same inner boundary condition.There is no shock in Fig. 7 (c & f) for α = 0.03.Since, at the outer boundary, the flow variables have very different values therefore it is not wise to conclude at α = 0.03, there is no shock.The angular momentum at the outer boundary i. e., λou is high, so the magnitude of infall velocity vr is low, and the flow remains subsonic till very close to the black hole.Therefore, there is no shock.

RESULTS: OSCILLATING SHOCKS
We study the effect that viscosity might have on a shocked accretion flow for a given outer boundary condition.We start with an inviscid solution and then go on increasing α.In the previous sections, we injected with a relatively shorter boundary, although the injection was subsonic i. e., M (rou) < 1.In order to study shock oscillation, we would need a larger computation region.It is clear that the accretion solution crucially depends on ϵ and λ.In Fig. 7 we fixed the inner boundary (ϵ & λ0) and changed α.Now we study how α affects accretion solutions with the same outer boundary.

Simulation set up
Viscosity transports angular momentum outward, but importantly it is dependent on Θ, apart from the shear.Therefore, it would be less effective in the supersonic region and more effective in the sub-sonic region.As a result, we see stronger angular momentum transport in the sub-sonic region.This implies there would be an accumulation of λ in the postshock region.Similarly, the dissipation of viscous heat in the post-shock disc is also greater compared to the pre-shock disc.
In the viscous transonic flow, a shock is stationary if the total pressure (including the ram pressure, gas pressure, and rotational barrier) remains conserved across the shock front.As we increase the viscosity parameter α, the angular momen-  tum will be redistributed, and therefore, beyond a certain critical value, it might destabilize the shock.
To properly study the effect of viscosity on accretion flow, we considered six models L1, L2, L3 and E1, E2, E3 (see, Table 1).We obtain analytical inviscid solutions corresponding to ϵ = 1.0001, but for λ = 1.76, , 1.77, 1.78, and called them models L1, L2 and L3, respectively.From the analytical solutions, we chose the rou, vrou, Θou, and λou.Keeping these values fixed as the outer boundary, we then increased the viscosity parameter α in the simulation code.The injection parameters of models E1, E2, and E3 on the other hand, is obtained for λ = 1.78 and ϵ = 1.00005, 1.0001, & 1.00015, respectively.From these analytical solutions the injection parameters rou, vrou, Θou, and λou are chosen (see Table 1).Again, keeping these injection parameters fixed, we change the viscosity parameter α.For all the models, we used 6826 uniform grids.The unit of time for models L1 -L3 and E1-E3 is tunit = 10 −2 (MBH/M⊙)s.
All these models are chosen such that the respective inviscid solutions harbours shock.For all the models, we start with the injection parameters with α = 0 till it reaches the steady state, and then we increase the α and study the changes induced by viscosity.Since all the models start with the inviscid flow, the angular momentum of the inviscid flow becomes the angular momentum at the outer boundary when viscosity is turned on.First, we study three cases where we kept ϵ same and studied flows with three different angular momentum, which we called L1 -L3 models.Then we kept the angular momentum the same but studied cases for three different ϵ and called them models E1 -E3.

Models L1 -L3
For the model L1, once the steady state is reached in the inviscid limit, we found a standing stable shock at around 40.78rg.Then the viscosity is turned on in the simulation.Viscosity is more effective in the subsonic, hot region of the flow.And post-shock disc, being hot and subsonic, will transport λ more efficiently, and there will be a pile-up of angular momentum in the post-shock disc.
This enhanced rotation of the disc near the shock front is likely to push the shock front outward.If the inviscid flow has a steady shock, then a small amount of viscosity will eventually produce a steady shock in another position.If the viscosity is increased further due to the piled up post-shock angular momentum, the shock position will overshoot the equilibrium position, and the shock starts to oscillate.The critical viscosity for the model L1 for the shock starts to oscillate is α = 0.025.Figure 8 a-f, is the comparison of the flow variables vr, M , ρ, p, Θ and λ, respectively, from model L1.The inviscid steady state solution (black, solid) is at code time t1 = 200 and the time-dependent viscous flow for α = 0.035 (red, solid) was obtained at t2 = 380.For the inviscid flow, the angular momentum is constant.However, for viscous flow, λ is much higher than the pre-shock flow, which causes the shock oscillation.The λ transport is much stronger in the post-shock region but in a small region near the shock front dλ/dr ∼ 0. This also causes a small spike in vr and M , while a tiny dip in ρ and Θ is also seen.
In Fig. 9, we have shown the time evolution of the shock position for different α for model L2. Figure 9a shows the shock position with time for α = 0.01, for which a stationary shock forms at around 61.02 rg.As viscosity is increased further, the stationary shock becomes unstable, and the shock front starts to oscillate.Time evolution of shock position with the same outer boundary condition but different viscosity parameters α = 0.02, α = 0.025 and α = 0.03 are shown in Fig. 9b, 9c and 9d, respectively.In this model, there is significant shock oscillation for α = 0.03.The shock oscillates with a mean position of about 113.49rg and has an amplitude of about 7.98rg.
For model L3, the amplitude and the mean shock position increase while frequency decreases with increasing α.In Fig. 10a, we compared the amplitude of the shock oscillation as a function of α for models L1 (green), L2 (red), and L3 (blue).Comparison of the corresponding oscillation frequency is presented in Fig. 10b and the mean position of oscillating shock is presented in Fig. 10c.The amplitude of the oscillation initially increases but tends to taper off with the increase of α.The amplitude for a given α increases with increasing λou.Frequency has an opposite trend, while the mean shock position follows the trend of the amplitude.In Fig. 10a, the oscillation starts at some critical viscosity parameter, say at α l marked by the black star, while ends with a black dot, representing an upper limit of viscosity αu beyond which there is no persistent oscillation.For model L2, α l = 0.016 and αu = 0.038.Therefore, the r sh vs t plot for α = 0.01 in Fig. 9a, shows no oscillation of the shock.If α > αu, then the shock oscillates but propagates outward and eventually leaves the computation domain, indicating that the force balance do not allow for a persistent oscillation.

Models E1 -E3
In the case of models E1-E3, the injection parameters are given in Table 1, and the computational domain has a 6827 uniform grid.Each grid has a size of 0.146rg, the same as model L1-L3.In other words, the resolution is kept exactly  .Snapshot of M and λ at times t 1 = 348 (black), t 2 = 351 (red), and t 3 = 588 (green) in code unit shown in (a) and (b).Time t 1 = 348 presents the inviscid steady state, t 2 = 351 present a dynamical state when a secondary shock (SS) forms, and in t 3 = 588 SS vanishes but primary shock oscillates for model E3.Viscosity parameter is α = 0.02. the same.The injected angular momentum of models E1-E3 is exactly the same as that of model L3, but from E1 to E3, we inject progressively hotter and slower matter.In Fig. 11a, 11b, 11c and 11d, we plot the shock location r sh with time of model E1, for α = 0.0075, α = 0.0125, α = 0.015 and α = 0.0175, respectively.The critical viscosity which triggers shock oscillation is α l = 0.0125.The oscillation amplitude and mean shock position increased with increasing α, but the frequency decreased.We zoomed the r sh vs t plot for α = 0.0175 (Fig. 11d) in Fig. 11e.We then compute the corresponding average acceleration fr in the immediate postshock flow at each time step and plotted them as a function of time.The average force is defined as fr = r * r sh frdr r * r sh dr Here r * is the position where the Mach number has its lowest magnitude or where d|M |/dr ∼ 0. And fr = dΦ/dr + λ 2 /r 3 − dp/(ρdr) is the net radial acceleration.The acceleration is averaged from the shock location to the region where the Mach number has the minimum value.The fr oscillates between positive and negative values, and it has a maximum value at the minimum shock location and a minimum value at the maximum r sh (Fig. 11f).And it is this oscillating average acceleration in the immediate post-shock region that causes the shock oscillation.The excess λ pile up in the post-shock disc due to the presence of viscosity is driving the shock front outward.While as the post shock region expands the thermal gradient force becomes weaker, and gravity manages to restore the shock to its original location, which the r sh overshoots and oscillation continues.However, beyond a certain value of α i. e., αu, the median shock location continues to move outward, and the shock becomes unstable.
In Fig. 12a, we plot and compare the amplitude of oscillation as a function of α for model E1 (red), model E2 (blue), and E3 (green).The black stars and dots indicate the two critical viscosity parameters α l and αu.In Fig. 12b, we plot the corresponding frequencies, and in Fig. 12c, we plot the variation of mean position with α.The character of the oscillations in models E1 and E2 are similar to the L3 model (i.e., amplitude/frequency increase/decrease with α).For model E3 (green), the amplitude starts to increase with α but eventually saturates.The frequency, on the other hand, increases with α, and it saturates too.The mean shock position monotonically increases for E1 and E2, but remains roughly the same for E3 with the increase of α.In Fig. 13a & b, we compare snapshots of M and λ at times t1 = 348 (black), t2 = 351 (red) and t3 = 588 (green) in code unit, for model E3 with α = 0.02.Time t1 = 348 represents the steady state, inviscid flow.At time t2 = 351 (red), a secondary shock is observed to form.At time t3 = 588 the secondary shock vanishes, but the primary shock continues to oscillate.It may be noted that the secondary shock is formed when the oscillating shock front (r sh ) approaches or exceeds 100rg.When the shock front recedes to a large distance from the horizon, then the temperature distribution of the post-shock disc varies significantly.Therefore, closer to the horizon, the rate at which the angular momentum will be transported and the rate at which λ would be transported near the shock front will be different and will cause one or two peaks in λ.This causes formation of secondary shocks (see, Lee et al. 2011Lee et al. , 2016)).The tendency to form a secondary shock was also seen in Fig. 8, where a local dip in λ corresponded with a local spike in vr and M , although the effect didn't develop into a secondary shock.In this figure, the secondary shock did form, but for a short time, and in later times, it washed away.This is quite different from the cases reported in Lee et al. (2011), where persistent secondary shock was reported.

Effect of shock oscillation
The post-shock is hotter and denser and so would emit highenergy photons.When the post-shock disc oscillates, the radiation emitted by it will also oscillate.Therefore, it is intriguing to investigate the time dependence of emission from such an oscillatory accretion disc.As radiative loss from the disc, we estimate the bremsstrahlung and synchrotron losses a posteriori from the disc as the representation of radiative processes.The emissivity resulting from bremsstrahlung (measured in erg cm −3 s −1 ) is described by Novikov & Thorne (1973) and given as, Te 1 + 4.4 × 10 −10 Te (43) and the emissivity due to the synchrotron (in erg cm −3 s −1 ) is given by Shapiro & Teukolsky (1983), At any time t, the total loss estimated posteriori would be Here, F(r sh ) represents the fitting function for the Comptonization parameter.This function is computed upon obtaining the accretion shock (similar to Kumar & Chattopadhyay 2014).The form of this analytic function is given by Where the values of parameters are A = 5.915, c1 = 100.527,c2 = 59.314, c3 = −6131787530.355,B = 1.972.For Qsyn, the magnetic field is assumed stochastic and is estimated assuming the magnetic pressure is 1% of the gas pressure.We followed Kumar et al. (2014) approach to estimate electron temperature.We assume the mass of the BH as MBH = 10M⊙.In Fig. 14, we plot the time series of shock oscillation (panel a), the time series of the total luminosity Lc (panel b) and the corresponding power density spectrum (panel c), for the model E1 and α = 0.02.Assuming an accretion rate of 0.1 ṀEdd we obtained the emission and then computed the power spectrum (PDS) for emission luminosity Lc We obtained a fundamental frequency of ν fund = 0.964 Hz and three harmonics.
In Fig. 15 a and b, we plot the time series of shock position r sh , the total power (Lc) for the model L1 and α = 0.03.The BH mass and accretion rate same as above.Interestingly, there appears to be a phase gap between the oscillation in emission and shock oscillation.To analyze the oscillation of radiation, we employ the methods described in Vaughan (2005).The PDS of the radiation is depicted in Fig. 15c.The PDS reveals that the quasi-periodic oscillation of radiation is characterized by the fundamental frequency ν fund = 7.49 Hz.The first significant harmonics are observed at 15.07 Hz (approximately 2ν fund ), 22.71 Hz (approximately 3ν fund ), and 30.04 Hz (approximately 4ν fund ) and so on.From Figs. 14 and 15, it is clear that the radiative efficiency of these flows is very low (∼ 10 −3 ).Realistic cooling may not qualitatively modify the results.We also included cooling in a very simplified manner and presented it in Appendix 5, which shows considering cooling may not dramatically change the solution.
It may be noted that the frequencies depicted in Figs. 10 & 12 are the fundamental frequencies of the shock oscillation.Although there is a phase shift between the oscillation of the radiation with that of the shock, the frequencies are the same.The peaks in emission are just after the shock minima.So, as the shock compresses the post-shock flow, making it hotter and denser the radiation peaks just after that.This has been seen before Lee et al. (2011Lee et al. ( , 2016)); Das et al. (2014).Table 2 lists all the models and various parameters related to shock oscillation.The frequencies mentioned are the fundamentals of oscillation.It shows that depending on the energy and angular momentum of the flow, the oscillations can range from frequencies that are less than Hz to a few Hz of frequencies.

SUMMARY AND DISCUSSION
In this study, we have conducted time-dependent numerical simulations to investigate the dynamics of viscous accretion flows surrounding black holes.For our investigations, we used a numerical code that implements the Eulerian TVD (Total Variation Diminishing) scheme along with a variable adiabatic index equation of state.Moreover, we have assumed a fully ionized electron-proton plasma.In addition, used the angular momentum density (ρ λ) as one of the state variables and the specific angular momentum (λ) as one of the primitive variables and then solved the equations of motion.As a result, this code ensures strict conservation of angular momentum in the absence of viscosity and ensures proper angular momentum transfer such that the angular momentum distribution is correct in the presence of viscosity.Therefore, with the use of CR EoS, we obtained a correct temperature distribution, and using the proper angular momentum equation instead of the equation for azimuthal speed, we obtained the correct angular momentum distribution within the framework of the Eulerian, upwind numerical scheme.Through rigorous testing, we demonstrated that the code accurately reproduces the angular momentum and temperature distribution similar to that in the analytical solution for viscous transonic flows (see Fig. 6).In order to be satisfied, we regenerated all types of accretion solutions with our simulation code, and as solutions reached a steady state, we compared them with the analytical steady-state solutions for α ≥ 0. These were the tests of steady-state solution where α is fixed but for different ϵ & λ0.We also regenerated the analytical accretion solutions for the same inner boundary condition by changing the outer boundary condition for a given α with our simulation code.And then finally, we compared timedependent accretion solutions by keeping the outer boundary the same but changing the α.We showed the accretions starting with no shock in the inviscid solution may harbour a shock due to the transfer of angular momentum.Of course, above critical viscosity, the shock will also disappear.The shock considered here is thin.In accretion discs, the r − ϕ component of the stress dominates over other components and that makes the shock front (if present) to be thin.In analytical/semi-analytical solutions, the shock front is an infinitesimally thin surface, but the numerical simulation code takes 5-6 cells to resolve the shock.Different dissipation processes may make the shock front spread over a small but finite width.
In our time-dependent study, we focused specifically on the behaviour of shock in the presence of viscosity.We considered six models to study them, three models varied the injection angular momentum but injection temperature and velocity is computed by keeping ϵ fixed, and three models varied the injection velocity and temperature according to varying ϵ values but at a given angular momentum at the outer boundary.In this paper, it has been demonstrated that for low values of viscosity parameter, the shock front in a disc tends to settle into a position, which might be different from the inviscid case.In the case of higher viscosity, the rate of angular momentum transfer increased, leading to a faster expansion of the shock front.As the shock front surpassed a potential equilibrium position, it initiated oscillations (Fig. 9, Fig. 11).This particular value of α can be denoted as critical viscosity.It is important to note that there is no specific value of the critical viscosity parameter but rather dependent on the initial outer boundary conditions ( vrou, Θou, λou).For instance, as the λou and Θou of the accretion flow increase for a given vrou, the value of critical viscosity decreases, but for increasing Θou for a given value of λou, it tends to have a higher critical viscosity parameter.This implies that the threshold at which the shock front undergoes oscillation is influenced by the specific characteristics of the flow, such as angular momentum and energy.In our simulations, we observed that the rate of outward angular momentum transport in the subsonic region is higher than in the supersonic region, and therefore, near the sink, the λ distribution tends to flatten out.The post-shock λ distribution is lower for steady-state, while in the case of oscillating shock, it is higher.It is the higher angular momentum transport that drives the shock oscillation.It has been shown in this paper that the viscosity transports angular momentum in the post shock disc more efficiently, causing the angular momentum to pile up.If the piling up angular momentum is such that along with the thermal gradient term, the average acceleration of the mass of post shock fluid pushes the shock front outward.As the shock front moves outward the thermal gradient force weakens, allowing gravity to bring it back, and then the shock front oscillates, generating what is called oscillating shocks.For model L3 i. e., higher λou, the mean shock position and the oscillation amplitude increase with increasing α.For models E1-E2, which imply hotter flow at the outer boundary, the mean shock position and its amplitude decrease with the viscosity parameter.But for E3 the mean shock position decreases and then becomes flat with increasing α.The power density spectra of the time series of the emission for models E1 and L1 show generally a C-type QPO with a higher Q factor and Lc is also low, but the QPO for E1 is broader compared to L1.Interestingly, the mean position of shock of model E1 with α = 0.02 is around few × 100rg, while for L1 (α = 0.03) it is ∼ 67rg.Oscillating a larger mass of fluid induces a degradation in the quality of oscillation.
In our simulations, we have neglected radiative cooling.However, we kept an option of including a simpler form of cooling by reducing the total viscous heat dissipation by a factor 0 ≤ fc < 1. Typically similar to what was done by Narayan & Yi (1994).But in the bulk of the paper, we considered fc = 0, i.e., no cooling, in order to study only the effect of viscosity.However, we showed in the appendix that even if cooling is considered to be 10% of the viscous heating, the solution is affected marginally.The median shock location would go inward marginally, and the oscillation frequency would increase marginally.However, the whole approach to cooling in this paper is very simplistic.Posteriori computation of bremsstrahlung and synchrotron emission from these accretion flows is also very low.However, radiative cooling is a tricky topic, and the complications involved are beyond the scope of the present paper (see for example Sarkar & Chattopadhyay 2019;Sarkar et al. 2020Sarkar et al. , 2023)).We would take up this study in all seriousness in the near future.
There are few examples of simulations of viscous accretion flow in the literature (Lanzafame et al. 1998;Lee et al. 2011;Das et al. 2014;Lee et al. 2016).These previous papers were all studied with the supersonic injections and fixed Γ EoS.Since viscosity is not very effective in the supersonic region, in all of these simulations, angular momentum distribution is almost constant or has a very smooth gradient in the preshock region.Viscosity is effective in the subsonic region, so in the subsonic post-shock region the angular momentum transport is very effective, resulting in a notable pile-up of angular momentum in the post-shock disc.Such a significant imbalance in angular momentum across the shock front causes large amplitude oscillation.On the contrary, in the present work, we injected the matter with subsonic velocities and used variable Γ EoS.Moreover, the injection parameters were chosen from analytical inviscid solutions, which harbours a steady shock.It is well known that an accretion shock solution will pass through the outer sonic point, then a shock, and then enter the black hole through the inner sonic point.The angular momentum transport between injection point rou and the outer sonic point is significant, and it is also similar in the post-shock region.So, the angular momentum transport in the post-shock and pre-shock disc is comparable.As a result, the oscillation amplitude is much less than in the previous simulations with supersonic injections.Hence, when the mean shock position is greater than 100rg, even then, persistent secondary shocks do not form as was reported by Lee et al. (2011Lee et al. ( , 2016)).We found a temporary secondary shock may form (model E3, α = 0.02), but it got washed away.It may also be noted that the effect of shock oscillation does not travel beyond the outer sonic point, and therefore, the shock oscillation does not affect the outer sonic point and the outer boundary.In other words, the supersonic region between the shock and the outer sonic point acoustically disconnects the shock oscillation from affecting the outer boundary.Therefore, as the shock oscillates, it is primarily the radiation from the shocked region that oscillates, and the radiation from the pre-shock region remains relatively unaffected by shock oscillation.
This is probably the first numerical simulation of viscous transonic accretion flows onto black holes using a variable Γ EoS.Moreover, the dependence of properties of timedependent accretion shocks on flow parameters is represented in Figs. 10 & 12, are also probably obtained for the first time.The frequency of the shock oscillation may decrease or increase, with the increase of viscosity parameter depending on the energy and angular momentum of the flow, and similarly, the mean position of the oscillating shock can increase or decrease with the increase of viscosity parameter.Finally, low-frequency C-type QPO around a stellar-mass black hole, ranging from frequencies less than a Hz to tens of Hz, can be easily explained by oscillating shock.

Figure 1 .
Figure 1.Schematic representation of accretion disc geometry.Here θ 0 is the co-latitude.
It has been widely used to study a range of astrophysical problems, including those described in Lee et al. (2011); Chattopadhyay et al. (2012); Raychaudhuri et al. (2021); Joshi et al.

Figure 2 .
Figure 2. The domain of multiple sonic points (dashed) and shock (solid) in ϵ − λ Parameter space in the inviscid accretion flows.
, we compare the theoretical solutions (absolute values of Mach number M = vr/cs with r) for the flow parameters corresponding to the points a, b, c, and d points of Fig.2.The simulation result is obtained by injecting the values of flow variables at an injection radius rou taken from the analytical solution at that point.Since for accretion vr < 0, so M is also negative, and that is why we plot absolute value of Mach number or |M |.In Fig.3a, we compare the theoretical accretion curve (black solid) with the simulation result (red open-circle) corresponding to point 'a' of Fig.2.The dotted line is an excretion-type flow that fulfills the outflow-type boundary condition.The injection radius for this simulation is also the outer boundary of simulation rou = 350, and injection flow variables at the rou are vrou = −1.722× 10 −2 , Θou = 7.869 × 10 −4 .The numerical simulation result follows the analytical curve closely and in fact, it becomes transonic around the analytical value of rco.The theoretically obtained sonic point is at rco = 126.414and is located at the crossing point of the dotted and solid black curves, while the sonic point obtained through numerical simulations is at 126.425, a discrepancy of only 0.0087%.In Fig.3b, we compare the analytical accretion (black, solid) with the simulation result (red open-circle) corresponding to the point 'b' of Fig. 2. The injection parameters are rou = 350, vrou = −1.706× 10 −2 , Θou = 7.871 × 10 −4 .

Figure 6 .
Figure 6.Comparison of λ vs r from numerical simulation (red open circle) and the analytic solution (black solid).The viscosity coefficient is α = 0.01.This figure corresponds to Fig. 5.The numerical solution is obtained using 4096 uniform grid cells.

Figure 12 .
Figure 12.Variation of amplitude (a) and frequency of the shock oscillation (b) with α assuming black hole mass to be M BH = 10M ⊙ .Flow parameter λou = 1.78.Red line is for model E1, blue line is for model E2, green line is for model E3.
Figure13.Snapshot of M and λ at times t 1 = 348 (black), t 2 = 351 (red), and t 3 = 588 (green) in code unit shown in (a) and (b).Time t 1 = 348 presents the inviscid steady state, t 2 = 351 present a dynamical state when a secondary shock (SS) forms, and in t 3 = 588 SS vanishes but primary shock oscillates for model E3.Viscosity parameter is α = 0.02.

Figure 14 .Figure 15 .
Figure 14.(a) Variation of shock position with time (in seconds), (b) Variation of total luminosity (ergs −1 ) with time, (c) Power density spectrum (PDS)of the emission for the Model E1 and α = 0.02.The BH mass is 10 M ⊙ .

Table 1 .
Details of injection parameters for different models.

Table 2 .
Properties of shock oscillation for different models with different values of α.