Cyclic and heteroclinic flows near general static spherically symmetric black holes

We investigate the Michel-type accretion onto a static spherically symmetric black hole. Using a Hamiltonian dynamical approach, we show that the standard method employed for tackling the accretion problem has masked some properties of the fluid flow. We determine new analytical solutions that are neither transonic nor supersonic as the fluid approaches the horizon(s); rather, they remain subsonic for all values of the radial coordinate. Moreover, the three velocity vanishes and the pressure diverges on the horizon(s), resulting in a flowout of the fluid under the effect of its own pressure. This is in favor of an earlier prediction that pressure-dominant regions form near the horizon. This result does not depend on the form of the metric and it applies to a neighborhood of any horizon where the time coordinate is timelike. For anti-de Sitter-like $\text{f}(R)$ black holes we discuss the stability of the critical flow and determine separatrix heteroclinic orbits. For de Sitter-like $\text{f}(R)$ black holes, we construct polytropic cyclic, non-homoclinic, physical flows connecting the two horizons. These flows become non-relativistic for Hamiltonian values higher than the critical value allowing for a good estimate of the proper period of the flow.


I. INTRODUCTION
General relativity is one of the most well tested theories in physics, however, there seem to be indications that it might be modified at sufficiently large scales (as well as small scales). The most important indication of the modification of general relativity comes from the observations made on the Supernova type Ia (SN Ia) and Cosmic Microwave Background (CMB) radiation [1][2][3]. These observations indicate that our universe is undergoing accelerated expansion. This could be explained by dark energy, and the vacuum energy in quantum field theories could have been used as a proposal for dark energy [4,5]. However, the problem with this proposal is that the vacuum energy in quantum field theory is much more than the dark energy required to explain the present rate of expansion of the universe. There seem to be serious limitations on modifying quantum field theories such that the vacuum energy is reduced to fit the amount of dark energy in the universe. In fact, it has been argued that such modifications will lead to a violation of the weak equivalence principle [6,7].
The action for general relativity has also been modified to explain the accelerated expansion of the universe, and currently f(R) gravity is one of the most well studied modifications of general relativity [8][9][10][11][12][13]. This is because the f(R) gravity theories are known to produce an accelerated expansion of the universe [14][15][16]. Furthermore, if a cosmological constant exists, it will not have any measurable effect for most astrophysical phenomena [17,18]. However, the f(R) gravity theories can have astrophysical consequences. In fact, astrophysical consequences have also been used to constraint certain type of f(R) gravity models [19,20]. So, it becomes both interesting and important to study astrophysical phenomena using f(R) gravity. Several methods for the static spherically symmetric solutions in f(R) gravity are studied in Refs [21,22]. Regular black holes in f(R) gravity are studied in Refs. [23][24][25]. Myung discussed the stability of f(R) black holes [26]. Further, there are many applications of f(R) gravity, e.g. gravity waves, brane models, effective equation approach, LHC test etc. [27][28][29] An important astrophysical effect of black holes is that they tend to accrete matter, and such accretion on a black hole have been thoroughly studied [30][31][32][33]. As the first studies of the accretion around a black hole were done by Bondi in the Newtonian framework [34], this effect is now known by the name of the Michel-type accretion. In his work, Bondi studied the hydrodynamics of polytropic flow, and demonstrated that settling and transonic solutions exist for the gas accreting onto compact objects. The relativistic versions the Michel-type accretion have also been studied using the steady state spherically symmetric flow of a test gas around a black hole [35,36]. It may be noted that the luminosity spectra and the effect of an interstellar magnetic field in ionized gases [37], the effect of radiative processes [37][38][39], and the the effect of rotation [40] on accreting processes have also been studied. Recently, the Michel-type accretion of perfect fluids for a black hole in the presence of a cosmological constant has also been studied [41][42][43]. Jamil and collaborators studied the effects of phantom energy accretion onto static spherically symmetric black holes and the primordial black holes and found the masses of black holes to decrease and vanishing near the Big Rip [44][45][46][47]. The accretion on topologically charged black holes of the f(R) theories and the Einstein-Maxwell-Gauss-Bonnet black hole has also been investigated by focusing on both inward and outward flows from the accretion disk [48,49]. Using the fact that data from the high-mass X-ray binary Cygnus X-1 has been used to constrain the values of the parameters for the f(R) gravity theories [50], in this paper, we will rather analyze some other aspects of the Michel-type accretion for a black hole in a theory of f(R) gravity.
The order of the paper is as follows. In Sec. II we discuss the general equations for spherical accretion including conservation laws for any static metric. We particularly show that the pressure of the perfect fluid for such spherically symmetric flows is, up to a sign, the Legendre transform of the energy density. This leads to a nice differential equation allowing the determination of the energy density, enthalpy, or pressure knowing one of the equations of state. In Sec. III, without restricting ourselves to a specific static black hole, we study the accretion phenomenon using the Hamiltonian dynamical system in the plane (r, v) where r is the radial coordinate and v is the three-dimensional speed of the fluid. We discuss sonic and non-sonic critical points for ordinary fluids as well as for non-ordinary matter. In Sec. IV we write down the metric for static spherically symmetric black hole in a particular model of f(R) gravity [51] and discuss some of its properties. In Sec. V we study the isothermal fluid and various subcases. There we provide examples of new solutions among which critical flows and purely subsonic flows with vanishing speed and divergent pressure on the horizon as well as separatrix heteroclinic orbits by restricting the analysis to an f(R) anti-de Sitter-like black hole. We also determine solutions that are purely supersonic and solution with transonic flows. We discuss the stability of some of these flows. In Sec. VI we apply the results of our Hamiltonian dynamical analysis to polytropic fluids. In Sec. VII we again consider the accretion of a polytropic fluid onto an f(R) black hole solution where the function f(R) is modeled by (a) Hu-Sawicki [52] and (b) Starobinsky [53] formulas. The last Section contains the conclusion and discussions of the above derivations.
Throughout the paper we have used the common relativistic notations. The chosen metric signature is (−, +, +, +) and the geometric units G = c = 1.

II. GENERAL EQUATIONS FOR SPHERICAL ACCRETION
In this section, in Sec. III, and in the first part of each of Sec. V and Sec. VI we consider any static spherically symmetric metric of the form without specifying the form of the metric coefficient f . Our results will apply to any black hole of that form and to any horizon in a neighborhood of which the time coordinate is timelike. In the second part of each of Sec. V and Sec. VI we consider some applications to an f(R) anti-de Sitter-like, to Schwarzschild, and to an f(R) de Sitter-like black holes.
In this section, we define the governing equations for spherical accretion. Here, we are considering the gas as a perfect fluid. We analyze the accretion rate and flow of a perfect fluid in f(R) gravity. For this, we define the two basic laws of accretion i.e. particle conservation and energy conservation. We assume that the fluid is simple containing a single particle species; the fluid could be made of different particle species with low reactions rates or no reactions at all. Let n be the baryon number density in the fluid rest frame and be the intrinsic four velocity of the fluid where τ is the proper time. We define the particle flux or current density by J µ = nu µ . From the law of particle conservation, there will be no change in the number of particles i.e. neither particles are created nor destroyed. In other words, we say that for this system, the divergence of current density is conserved where ∇ µ is the covariant derivative. On the other hand, the stress-energy (SET) for a perfect fluid is given by where e denotes the energy density and p is the pressure. The Michel-type accretion is steady state and spherically symmetric [41][42][43], so all the physical quantities (n, e, p, u µ ) and others that will be introduced later are functions of the radial coordinate r only. Furthermore, we assume that the fluid is radially flowing in the equatorial plane (θ = π/2), therefore u θ = 0 and u φ = 0. For ease of notation we set u r = u. Using the normalization condition u µ u µ = −1 and (1), we obtain, On the equatorial plane (θ = π/2), the continuity equation (3) yields or, upon integrating, where C 1 is a constant of integration. This shows that, in a unit of proper time, the particle flux πr 2 nu through a sphere a radius r remains constant for all r. The thermodynamics of simple fluids is described by the two equations [61] dp = n(dh − Tds), de = hdn + nTds, (8) where T is the temperature, s is the specific entropy (entropy per particle), and is the specific enthalpy (enthalpy per particle) 1 .
A theorem in relativistic hydrodynamics [61,62] states that the scalar hu µ ξ µ is conserved along the trajectories of the fluid: where ξ µ is a Killing vector of spacetime generator of symmetry. In the special case we are considering in this work ξ µ = (1, 0, 0, 0) is timelike yielding where C 2 is a constant of integration. This equation can be derived directly upon evaluating where we have used T µ ν = nhu µ u ν + (nh − e)δ µ ν . Since, the flow is stationary, any time derivative vanishes (∇ t (nh − e) ≡ 0), hence the result. If the fluid had a uniform pressure, that is, if the fluid were not subject to acceleration, the specific enthalpy h reduces to the particle mass m and Eq. (10) reduces to mu µ ξ µ = cst along the fluidlines. This is the well know energy conservation law which stems from the fact that the fluid flow is in this case geodesic. Now, if the pressure throughout the fluid is not uniform, acceleration develops through the fluid and the fluid flow becomes non-geodesic; the energy conservation equation mu µ ξ µ = cst, which is no longer valid, generalizes to its inertial equivalent [61] hu µ ξ µ = cst as expressed in Eqs. (10) and (11).
It is well known that a perfect fluid (4) is adiabatic; that is, the specific entropy is conserved along the evolution lines of the fluid (u µ ∇ µ s = 0). This is easily established using the conservation of the SET, Eq. (3), and the second equation in (8). First, rewrite T µν as 1 If m is the baryonic mass, then ρ = mn is the mass density. Now, if h = h/m and s = s/m denote the enthalpy and entropy per unit mass, respectively, then ρh = nh and ρs = ns. In terms of (h, s, ρ), Eqs. (8) and (9) take the forms dp = n(dh − Tds), de = hdρ + ρTds, and h = (e + p)/ρ. nhu µ u ν + (nh − e)g µν , then project the conservation formula of the SET onto u µ In the special case we are considering in this work where the fluid motion is radial, stationary (no dependence on time), and it conserves the spherical symmetry of the black hole, the latter equation reduces to ∂ r s = 0 everywhere, that is, s ≡ const.. Thus, the motion of the fluid is isentropic and equations (8) reduce to dp = ndh, de = hdn.
Equations (7) and (11) are the main equations that we will use to analyze the flow of a perfect fluid in the background of f(R) black hole.
Another formula that will turn useful in the subsequent sections is the barotropic equation. Notice that the canonical form of the equation of state (EOS) of a simple fluid is e = e(n, s) [62]. Since s is constant, this reduces to the barotropic form From the second equation (14) we have h = de/dn yielding where the prime denotes differentiation with respect to n. Now, the first equation (14) yields which we integrate by parts to derive Here we identify, up to a sign, the Legendre transform of the energy density F. This conclusion is purely thermodynamic and it does not depend on the symmetric properties of the flow (presence of a timelike Killing vector and spherical symmetric flow); rather, it is valid for any isentropic flow (s constant everywhere). The conclusion states that the pressure is the negative of the Legendre transform of the energy density and that an EOS of the form p = G(n) is not independent of an EOS e = F(n). The relationship between F and G can be derived upon integrating the first differential equation In a locally inertial frame, the three-dimensional speed of sound a is given by a 2 = (∂p/∂e) s [63]. Since the entropy s is constant, this reduces to a 2 = dp/de. Using (14), we derive a useful formula needed for the remaining sections Using (16), this reduces to Another useful formula is the three-velocity of a fluid element v as measured by a locally static observer. Since the motion is radial in the plane θ = π/2, we have dθ = dφ = 0 and the metric (1) implies the decomposition in the standard special relativistic way [64,65] as seen by a locally static observer. The latter measures proper distances and proper times by dℓ = dr/ f and dτ 0 = f dt corresponding to radial dr and time dt changes, respectively, and measures the three-velocity v of the fluid element by where we have used u r = u = dr/dτ, u t = dt/dτ, u t = − f u t , and (5). This implies and (7) becomes In relativistic hydrodynamics one usually derives the above formulas on considering the woldlines of a fluid element and that of a locally static observer. If u and u 0 are the respective four-velocities, we have [62,66] where U is the relative four-velocity, that is, the velocity of the observer attached to the fluid element relative to the locally static observer with the property u 0 · U = 0, where the dot represents the scalar product with respect to the metric (1). Γ is the Lorentz factor Γ ≡ −u 0 · u = dτ 0 /dτ [62,66]. In the case of radial motion in the θ = π/2 plane, we have u = (u t , u, 0, 0) = u t ∂ t + u∂ r , Here u t and u = u r are as defined in (2) and V r = dr/dτ 0 = f v. Since ∂ r is not a unit four vector, rather it is v, and not V r , the three velocity that the locally static physical observer, who uses the orthonormal basis (∂ t / f , f ∂ r , ∂ θ /r, ∂ φ /r), measures. Squaring (26) we obtain where we have used U · U = g rr V r V r = v 2 in the last expression. The expressions (24) are rederived from (26), (27), and (28). All the above expressions remain valid for an observer outside the horizon, more precisely, for an observer where the time coordinate is timelike. We define the value v h of v on the horizon(s) r h as the limit of the continuous three velocity field v(r) as r approaches r h from within the region where the time coordinate is timelike ( f > 0):

III. HAMILTONIAN SYSTEMS
We have derived two integrals of motion (C 1 , C 2 ) given in (7) and (11). Either of these integrals, or any combination of them, can be used as a Hamiltonian for the fluid flow. The simplest Hamiltonian system has one degree of freedom, in which case the Hamiltonian H is a two-variable function (x, y). Let H be the square of the lhs of (11): Now, we need to fix the two dynamical variables (x, y) on which H depends and the time variablet of the Hamiltonian dynamical system. There are different ways to fix the dynamical variables; one may choose (x, y) to be (r, u) [43], (r, v 2 ) [43], (r, n) [67], (r, h), or even (r, p). The time variablet for the dynamical system is any variable on which H (30) does not depend explicitly so that the dynamical system is autonomous. In Sec. II we have seen that, under the symmetry requirements of the problem, h is an explicit function of the baryon number density n only; this applies to the pressure p too. So, if (x, y) are chosen to be (r, h) (resp. (r, p)), the Hamiltonian (30) takes the form where we have used (7) (resp. H = h(p) 2 f (r) + C 2 1 r 4 n(p) 2 ). This conclusion does not extend to other dynamical variables, that is, if one chooses (x, y) to be, say, (r, v), it is not true to assume h = h(r) or h = h(v), for, by (7) and (24), n is a function of (r, v) and so is h. With h = h(r, v), the Hamiltonian (30) of the dynamical system reads where we have used (24) to eliminate u 2 from (30). We have thus fixed the dynamical variable to be (r, v). No use has been made of (7) to derive (32); use of it will be made in the derivation of the critical points (CPs), particularly, of the sonic points. From now on, partial derivatives will be denoted as

A. Sonic points
In the remaining part of this section, we assume that the parametric Hamiltonian of the dynamical system is given by (32). In this section we use (32) to derive the CPs of the dynamical system and derive them in the Appendix B VIII C using (31).
With H given by (32), the dynamical system readṡ (here the dot denotes thet derivative). In (33) it is understood that r is kept constant when performing the partial differentiation with respect to v in H ,v and that v is kept constant when performing the partial differentiation with respect to r in H ,r . We will keep using this simple notation in the subsequent steps of this section. The CPs of the dynamical system are the points (r c , v c ) where the rhs's in (33) are zero. Evaluating the rhs's we find The rightmost formula in (20) yields ; (37) and if v is kept constant we have the equation r 2 n f = const. which upon differentiating with respect to r we obtain Finally, the system (33) readṡ Let us assume that h is never zero and finite (the same applies to n). The rhs's vanish if where f c = f (r)| r=c and f c,r c = f ,r | r=c . The second equation expresses the speed of sound at the CP, a 2 c , in terms of r c which will allow to determine r c once the EOS a 2 = dp/de [or e = F(n)] is known. The remaining needed ingredient is a simplified expression for n/n c . If we write the constant C 2 1 in (25) as where we have used (41). Using this in (25) we obtain As we shall see in the subsequent sections, there will be two types of fluid flow approaching the horizon, in the one type the speed v vanishes and in the other one the speed approaches that of light in such a way that the ratio (1 − v 2 )/ f may remain finite. In the former type of motion, the number density n diverges on the horizon independently of the expression of f . An expression for u 2 c is derived upon substituting (41) into (24), then making use of (42) Another sonic CP is the point corresponding to f c = 0 and a 2 c = 1. But the roots of f c = 0 may coincide with the horizons r h of the black hole. This implies that the fluid becomes ultra-stiff as it approaches the horizon where r c = r h (the fluid is not necessarily ultra-stiff for all r). This conclusion does not apply to f(R) gravity only; rather, to any static spherically symmetric metric of the form (1). To the best of our knowledge, this result has not been announced elsewhere. Now, by (25), since f c = 0 we must necessarily have v 2 c = 1. This point, however, may fail to behave as a CP in the mathematical sense, for the rhs's of (39) and (40) may become undetermined or may have nonzero values there. This point, (r = r h , v = 1), may behave as a focus point as we shall see in the next section.

B. Non-sonic critical points
From (39), we see that f c = 0 and f c,r c = 0 may lead to a non-sonic CP. However, this CP would be a double root of f = 0, which is out of the scope of this paper where we only consider non-extremal black holes.
Another obvious CP, which lies within the scope of f(R) gravity, corresponds to h(r c ) = 0 (39) and (40). This is not possible for ordinary matter but is the case for non-ordinary matter with negative pressure. When this is the case, h may vanish at some point with no special constraint on v 2 and a 2 . This means that for non-ordinary fluids, the flow may not become transonic at all. We will not pursue this discussion here, for it is out of the scope of this work. In the next section, however, we will pursue this discussion for ordinary matter where it is generally admitted that "the flow must be supersonic at the horizon, though it is necessarily subsonic at a large distance" [68]. We will explicitly show, through physical solutions, the existence of subsonic flow for all values of the radial coordinate. Moreover, the speed of the flow vanishes as the fluid approaches the horizon, so the flow does not necessary become supersonic nor transonic near the horizon [69,70]. Our conclusion remains true even for the Schwarzschild black hole. We believe that the use of standard methods for tackling the accretion problems has masked many features of them.
The conclusions made in this section, concerning the sonic CP [from (39) to (45)], do not apply to f(R) gravity only, for we have not fixed the form of the metric coefficient f yet; they apply to any static metric with g tt = −1/g rr and g θθ = r 2 .
Applications are given in the following sections where we consider three models of f(R) gravity.

IV. BLACK HOLE IN f(R) GRAVITY
Recently, an interesting model of f(R) gravity has been proposed [51], and the motion of test particles around a black hole in this theory has been investigated. The Lagrangian for this model of f(R) theory is given by [51], where Λ is the cosmological constant, R c is a constant of integration 2 , and α, d are free parameters of this the-2 R c is merely a constant of integration which is used to balance the dimensions of R. Its value, which "is not sensitive to the SNIa data" [54], is not known by any physical theory and can only be determined using astronomical constraints as suggested by Safari and Rahvar [54].
ory. The limit that is relevant for astrophysical scale corresponds to R ≫ Λ and d 2 (6α 2 ) −1 R ≫ 2α. In this limit, we obtain f(R) This limit constrains the accelerating expansion [54]. It is useful to introduce a parameter β = α/d in terms of which both limits of the theory can be studied [51]. In this theory, the metric for a spherically symmetric black hole with mass M takes the form, If Λ = 0, (47) reduces to a special case of Kiselev black hole [57,58] and if β = 0, (47) reduces to Schwarzschildde-Sitter or Schwarzschild-anti-de-Sitter black hole.
The present model of f(R) can explain the flat rotation curve of galaxies, consistent with solar system tests and also explains the pioneer anomaly/acceleration. For details concerning the motivation for this particular model of f(R) theory, we refer the reader to the original work by Saffari and Rahvar [54]. Of course the present analysis can also be done for other f(R) black holes such as Eq. 32 of Ref [55] and will be reported elsewhere. However due to the generality of our work, further analysis will be trivial as was the case with f(T) gravity black holes [56].
It is well-known that f(R) theory has a representation equivalent to a particular class of scalar-tensor (ST) theories namely, the Brans-Dicke (BD) theory i.e. a scalar field being non-minimally coupled to gravity or curvature with vanishing kinetic term of the scalar field. This description holds for both metric and palatini f(R) theories [71,72]. Furthermore, the no-hair theorem for black holes in a general ST theory suggests that the Schwarzschild solution is the only asymptotically flat, exterior, vacuum, static and spherically symmetric solution to ST theory [73]. However, it does not rule out the existence of non-asymptotically flat ST black holes without hair. For instance, the Reissner-Nordström Anti-de Sitter kind of topological black holes are derived in BD-Maxwell ST theory [74]. In the same context, we study a non-asymptotically flat f(R) black hole.
The roots of f = 0, or equivalently, the roots of P = 0, where P ≡ 3r f = −Λr 3 + 3βr 2 + 3r − 6M is a polynomial of degree 3, determine all possible horizons of (47). If Λ > 0, the equation P = 0 has always some negative root, which we ignore because of the physical singularity at r = 0, and it may have two positive roots or a double positive root depending on the values of its coefficients. These two positive roots, if any, determine the event and cosmological horizons. In this case, the fluid flow would be confined in the space region enclosed by the two horizons. If there are no positive roots, the metric coefficient g tt is positive for all r > 0; this case is not interesting.
We will be interested in the cases where the positive roots of P = 0 are single. Assuming Λ < 0 (anti-de Sitter-like black hole) and β ≥ 0, then if β 2 > −Λ, P = 0 has either two negative roots and one positive root or one double negative root and one positive root; if 0 ≤ β 2 ≤ −Λ, P = 0 has one single positive root. On converting the polynomial P(r) into the Weierstrass polynomial w(z) ≡ 4z 3 − g 2 z − g 3 by the transformation r = z + β/Λ, we can parameterize the roots of P = 0 based on the parametrization of the roots of w(z) as given in the Appendix A VIII [60]. The horizon is given by if P = 0 has at least two real roots; if P = 0 has only one real root. Here g 2 and g 3 are defined by and ∆ and the angle 0 ≤ η ≤ π are defined as in Eqs. (A.2) and (A.4), respectively. Now, assuming Λ > 0 (de Sitter-like black hole) and β ≥ 0, P = 0 has always one negative root and will have two positive roots, corresponding to an event horizon r eh and a cosmological horizon r ch > r eh if 2(β 2 + Λ)r + > 6MΛ − β where r + is the positive root of P ′ (r) = 0. When this the case, the roots are given where g 2 and g 3 are defined by (50). ∆ and the angle 0 ≤ η ≤ π are defined as in Eqs. (A.2) and (A.4), respectively. To have a common notation with the case Λ < 0, we will for short denote r eh and r ch by r h . The scalar invariants R, R µν R µν , and R µνσρ R µνσρ are given by which reduce to the Schwarzschild values I 1 = I 2 = 0 and I 3 = 48M 2 /r 6 if β = Λ = 0. Clearly r = 0 is the curvature singularity, which is not removable.

V. ISOTHERMAL TEST FLUIDS
Isothermal flow is often referred to the fluid flowing at a constant temperature. In other words, we can say that the sound speed of the accretion flow remains constant throughout the accretion process. This ensures that the sound speed of accretion flow at any radii is always equivalent to the sound speed at sonic point [75]. Here our system is adiabatic, so it is more likely that the flow of our fluid is isothermal in nature. Therefore, in this section we find the general solution to the isothermal equation of state of the form p = ke, that is of the form p = kF(n) (15) with G(n) = kF(n) (19). Here k is the state parameter constrained by (0 < k ≤ 1) [42]. Generally, the adiabatic sound speed is defined as a 2 = dp/de. So by comparing the adiabatic sound speed to the equation of state, we find a 2 = k.
The differential equation (19) reads yielding where we have chosen the constant of integration 3 so that (9) and (16)    where the dot denotes differentiation with respect to the new timet.
For a subsequent physical discussion we need an expression for the pressure. With p = ke, we obtain upon substituting (44) Since the Hamiltonian (59) remains constant on a solution curve, if the latter approaches the horizon (any horizon) from within the region where t is timelike, f approaches 0, and so the speed v must either approach 1 or 0 so that the Hamiltonian retains the same constant value (otherwise, the Hamiltonian would always assume a 0 value on the horizon regardless its constant value elsewhere). In former case (v → 1), the pressure (60) may remain finite in a neighborhood of the horizon. In the latter case (v → 0), the pressure diverges as the solution curve approaches the horizon. This is a very general conclusion which holds for any metric coefficient f and any horizon of the black hole. If the latter is of de Sitter type (Λ > 0), a pressure-dominant region may form near both the event and cosmological horizons. This is in favor of a proposal that a pressuredominant region would form near the horizon [76].
If f (r) = 0 has a single root as r approaches r h (corresponding to an event, a cosmological, or any horizon in a neighborhood of which t is timelike), which is our case, then, in the latter case (v → 0), as the curve approaches the horizon f ∼ (r − r h ) and v 2k ∼ f 1−k , thus v 2 ∼ (r − r h ) (1−k)/k . Using this in (60) we see that the pressure diverges, as the curve approaches the horizon, as If r h is a double root of f = 0, we obtain Before we proceed further let us see what the constraints on k to have a physical flow are. Along a solution curve, the Hamiltonian of the dynamical system (59) is constant [where the constant is proportional to C 2 (11)]. A global flow solution that extends to spatial infinity corresponds to where (α > 0, v 1 , |v ∞ | ≤ 1) are constants. Inserting this in the Hamiltonian (59) reduces to : (b): Using the metric (47), each case splits into two subcases as follows.
Thus, for ordinary fluids we deduce and for non-ordinary fluids (−1 ≤ k < 0) we deduce On comparing the leading terms in the expansion (62), we see that the fluid flow for ordinary matter is faster at spatial infinity than it is for non-ordinary matter.
In the following we will analyze the behavior of the fluid by taking different cases for the state parameter k. For instance, we have k = 1 (ultra-stiff fluid), k = 1/2 (ultra-relativistic fluid), k = 1/3 (radiation fluid) and k = 1/4 (sub-relativistic fluid). For the case of the metric (47), Eq. (42) reduces to and we keep in mind that a 2 = k in (59). The system (59) and (73)  Ultra-stiff fluids are those fluids in which isotropic pressure and energy density are equal. For instance, the usual equation of state for the ultra-stiff fluids is p = ke i.e. the value of state parameter is defined as k = 1. This reduces (42) or (73) to f c = 0, thus r c = r h (48,49). The Hamiltonian (59) reduces to Since the Hamiltonian in Eq.(74) is a constant, one immediately obtains 4 v ∼ 1/r 2 .
It is clear from (74) that the point (r, v 2 ) = (r h , 1) is not a CP of the dynamical system, as was noticed in the previous section. Notice that H no longer depends on f ; thus, this expression, and the following conclusions, are valid for any metric of the form (1). From (74) we see that, for physical flows (|v| < 1), the lower value of H is H min = 1/r 4 h : H > H min . As shown in Fig. 1, physical flows are represented by the curves sandwiched by the two black curves, which are contour 4 For the cases k = 1 and k = 1/2 we have expressed explicitly v as a function of r as in Eqs. (75) and (83); it is possible to do the same for the other cases k = 1/3 and k = 1/4 [see Eqs. (89) and (92)] but the expressions of v(r) would be cumbersome. That's why we preferred a numerical analysis in this section. It is worth mentioning that the expressions (75) and (83) may be derived from the metric and the conservation laws using the classical approach for accretion [35].   (Fig. 1).
plots of H(r, v) = H min . The upper curves where v > 0 correspond to fluid outflow or particle emission and the lower curves where v < 0 correspond to fluid accretion. If H 0 > H min is the value of the Hamiltonian on a solution curve, then in the (r, v) plane the curve is the plot v = ±1/( √ H 0 r 2 ). Using this we can evaluate all the other quantities, for instance (44) becomes for any solution curve H 0 > H min = r −4 h , and n n c for the solution curve through (r, v 2 ) = (r h , 1) (H 0 = H min ), which all depend on f . A contour plot of H (74), depicted in Fig. 1, shows two type of motion: (a) purely subsonic accretion (black, magenta, or blue curves where v < 0) or subsonic flowout (black, magenta, or blue curves where v > 0) for H > H min = r −4 h , and (b) purely supersonic accretion or flowout (along the red and green curves) for H < H min = r −4 h . The flow in (b), along the green and red curves, is however unphysical, for the speed of the flow exceeds that of light on some portions of the curves. A brief elaboration is given in Table I. B. Solution for ultra-relativistic fluid (k k k = 1/2) Ultra-relativistic fluids are those fluids whose isotropic pressure is less than the energy density. In this case, the equation of state is defined as p = e 2 yielding k = 1/2. Using this expression in (73) reduces to This polynomial has always one and only one positive root if Λ < 0 and β ≥ 0. Converting this polynomial into the Weierstrass one w(z) by the transformation r c = z + 3β/(2Λ), the CP r c is given either by (see Appendix A) if Q = 0 has at least two real roots or by (80) if Q = 0 has only one real root. Here g 2 and g 3 are defined by and ∆ and the angle 0 ≤ η ≤ π are defined as in In the limit β → 0, we recover the Schwarzschild antide Sitter spacetime and Eq. (79) reduces to The Hamiltonian (59) takes the simple form It is clear from this expression that the point (r, v 2 ) = (r h , 1) is not a CP of the dynamical system. For some given value of H = H 0 , Eq. (82) can be solved for v 2 . We find where g(r) ≡ f /(H 0 r 4 ). The plot in Fig. 2 , 0), where the speed vanishes, then a subsonic flowout until (r c , v c ), followed by a supersonic flowout, or (4) (lower plot) a subsonic accretion followed by a supersonic accretion which ends inside the horizon. In the upper plot, we have a supersonic outflow followed by a subsonic motion. The summary of this is given in Table II. The fluid flow in Type (3) from (r c , −v c ) to (r c , v c ) describes a heteroclinic orbit that passes through two different saddle CPs: (r c , −v c ) and (r c , v c ). It is easy to show that the solution curve from (r c , −v c ) to (r c , v c ) reaches (r c , v c ) ast → −∞, and the curve from (r c , v c ) to (r c , −v c ) reaches (r c , −v c ) ast → +∞; we can change the signs of these two limits upon performing the transformation t → −t and H → −H.
The flowout of the fluid, which starts at the horizon, is caused by the high pressure of the fluid, which diverges there (61): The fluid under effects of its own pressure flows back to spatial infinity.
It is clear from Fig. 2 that, after watching the subsonic branches of the blue and magenta solution curves, there is no way to support the claim, recalled at the end of Sec. III, that "the flow must be supersonic at the horizon" [68]. For these new solutions the speed of the fluid increases during the accretion from 0, according to the analysis made from (62) to (72), to some value below v c where dv/dr = 0, then decreases to 0 at the horizon, and the process is reversed during the flowout. It is easy to show, using (83), that the point where the speed is maximum is r c , as shown in Fig. 2. Thus, the flow does not necessary become supersonic nor transonic near the horizon [69,70]. This conclusion does not depend on the presence of a negative cosmo-   (Fig. 2).
logical nor on a nonvanishing constant β: such solutions exist even for a Schwarzschild black hole, as the subsonic branches of the blue and magenta solution curves in Fig. 3 show.
Curiously enough, such solutions were never discussed in the literature. This is probably due to the fact that the pioneering works on this subject did not employ the Hamiltonian dynamical system approach to tackle the problem. These new solutions are related to the instability and fine tuning problems in dynamical systems. To see that consider the asymptotic behavior of (82). Since f ∼ −(Λ/3)r 2 as r → ∞ and since H remains constant on a solution curve, we must have v ∼ v 1 r −1 (v 1 < 0 during accretion), which agrees with (62) and (66). Asymptotically, Eq. (82) reads which is used to determine the value of |v 1 | by Notice that as |v 1 | increases, H ∞ decreases. Now consider the lower plot of Fig. 2 where the subscript "b" is for black. If one decreases the value of the asymptotic speed, that is, the value of |v 1 | by ǫ: |v 1 | → |v 1b | − ǫ, as is the case of the subsonic magenta curve of Fig. 2, then H ∞ increases by a corresponding amount: This small perturbation in the value of |v 1 | leads the flow to completely change course, by deviating from the black critical curve, and to undergo a purely subsonic motion along the subsonic magenta curve. Conversely, a small increase in the value of the asymptotic speed (of the coefficient |v 1 |) would lead the flow to follow the red curve adjacent to the black critical curve. Thus, the black critical curve is certainly unstable and in practical situations it would not be easy to fix the value of |v 1 |, which is an average value for the pressure is not zero, by fine tuning it to have a critical motion, that is, a motion that becomes supersonic beyond the CP and reaches the speed of light as the fluid approaches the horizon.
This stability issue is related to the character of the CPs (r c , −v c ) and (r c , v c ) that are saddle points of the Hamiltonian function. As is well known saddle points of the Hamiltonian function are also saddle points of the Hamiltonian dynamical system. Further analysis of stability requires linearization of the dynamical system and/or use of Lyapunov's theorems [77][78][79] and their variants [80].
Another type of instability is the flowout that starts in the vicinity of the horizon (r = r h + 0 + , v = 0 + ) under the effect of a divergent pressure. This flowout is unstable, for it may follow a subsonic path (the magenta or blue curves) or a critical path (the black curve) through the CP (r c , v c ) and becomes supersonic with a speed approaching that of light. From a cosmological point of view, this point (r = r h , v = 0) looks like an attractor where solution curves converge and a repeller from where the curves diverge [80].
The motion along the rightmost branches of the green and red curves is unphysical. Along the leftmost branches of these curves, we have an accretion starting from the leftmost point of the branch until the horizon where the speed vanishes and the pressure diverges, followed by a flowout back to the same starting point. To realize such a flow one needs to have a sink and source at the leftmost point of these branches.

C. Solution for radiation fluid (k k k = 1/3)
Radiation fluid is the fluid which absorbs the radiation emitted by the black hole. It is the most interesting case in astrophysics. Here, the value of state parameter k = 1/3. Eq. (73) leads us to which is solved by The Hamiltonian (59) takes the simple form It is clear from this expression that the point (r, v 2 ) = (r h , 1) is not a CP of the dynamical system. Eq. (89) can be solved for v 2 , and a contour plot of it can be depicted, which reveals the same characteristics of the plot shown in Fig. 2; We observe the same types of motion as in the case k = 1/2. Sub-relativistic fluids are those fluids whose energy density exceeds their isotropic pressure. Taking the value of the state parameter k = 1/4, Eq. (73) leads to This polynomial has either two distinct positive roots or a double positive root if Λ < 0 and β ≥ 0. Converting this polynomial into the Weierstrass one w(z) by the transformation r c = z − β/(2Λ), the two CPs r c1 < r c2 are given by (see Appendix A) where g 2 and g 3 are defined by It is clear from this expression that the point (r, v 2 ) = (r h , 1) is not a CP of the dynamical system. A contour plot of H (92) is depicted in Fig. 4 in the (r, v) plane. There are two saddle points (r c1 , v c ) and (r c1 , −v c ) and two centers (r c2 , v c ) and (r c2 , −v c ). Let (r rm , v c ) and (r rm , −v c ) be the rightmost points of the upper and lower plots, respectively. If we assume that dv/dr remains continuous as the fluid crosses the saddle CPs, the accretion motion starts from the rightmost point (r rm , −v c ) on the black curve in the lower plot. If the motion is subsonic it proceeds along the upper branch in the lower plot, goes through the CP (r c1 , −v c ), then crosses the horizon. Otherwise, if the motion is supersonic it proceeds along the lower branch in the lower plot, goes again through the CP (r c1 , −v c ) until v vanishes as the fluid approaches the horizon [this is obvious from (92) where v vanishes whenever f does too], then the fluid goes again through the CP (r c1 , v c ) and follows the upper branch of the upper plot undergoing a supersonic motion until the rightmost point of the upper plot (r rm , v c ). First, by similar arguments as those given in the case k = 1/2, it can be shown that such motion is unstable.
Secondly, the motion may become periodic but it is too hard to achieve that by (a) fine tuning the speed of the fluid at (r rm , −v c ) and (b) realizing a source at (r rm , −v c ) and a sink at (r rm , v c ).
The fluid flow along the branch of the curve from (r c , −v c ) to (r c , v c ) describes a heteroclinic orbit that passes through two different saddle CPs: (r c , −v c ) and (r c , v c ). It is easy to show that as the flow approaches, from within the heteroclinic orbit, one or the other saddle CP the dynamical-system's timet goes to ±∞.
Here again the flowout of the fluid, which starts at the horizon, is caused by the high pressure of the fluid, which diverges there (61).
As we have done in the case k = 1/2, we consider the fluid flow where r decreases but v > 0 or r increases but v < 0 as unphysical since the fluid is taken as a test matter and we have neglected its backreaction on the metric of the black hole. As far as a fluid element is taken as a test particle, such a motion is not possible in the background of the black hole metric. This is why a flow along a closed path in Fig. 4, or "homoclinic" as some authors call it, is unphysical. We do not know if homoclinic orbits exist in a more realistic model where the backreaction of the fluid is taken into consideration.
For the clarity of the plot, Fig. 4 has been plotted for unphysical parameters M = 1, β = 0.5, and Λ = −0.075; for astrophysical values of the the parameters (Λ → 0 − ), the difference r c2 − r c1 becomes so large to be represented on a sheet of paper. The constraint that two CPs exist is to have two positive roots for the polynomial in (90): N(r) = Λr 3 + 3β 2 r 2 + 6r − 21M. With Λ < 0 and β > 0, the polynomial has a local minimum (at some negative value of r) and a local maximum at The heteroclinic orbit exists if N(r c ) = 0 has two positive CPs; that is, if N(r s ) > 0 yielding generalizing the expression derived in Ref. [81]. This should be read as a constraint on β. In the limit Λ → 0 − , this reduces to and the expressions of the two positive CPs and the horizon read It is easy to show that r c1 > r h .
In the astrophysical limit Λ → 0 − we find, for general values of k, the following constraints on β In this limit, the CPs are expressed as while the expression of r h (96) is independent of k.

VI. POLYTROPIC TEST FLUIDS
A very interesting approach to describe the motion of fluid is by constructing its models. The prototype of such model is Chaplygin gas. The Chaplygin gas model leads to very interesting results. Some of them are discussed in Ref [82][83][84][85][86]. There are many variations of the Chaplygin gas model have been proposed in the literature. One of them is the modified Chaplygin gas model [87,88]. In astrophysics, the modified Chaplygin gas is the most general exotic fluid. Its equation of state is: where A and B are constants and (0 < α < 1). If we put A = 0, B = −k and α = −γ, we get the polytropic equation of state i.e. p = G(n) = Kn γ , where K and γ are constants. For ordinary matter, one generally works with the constraint γ > 1. In this work, we only observe the constraint γ = 1. Inserting p = G(n) = Kn γ in the differential equation (19) yields The solution provides the energy density e = F by where a constant of integration has been identified with the baryonic mass m. This yields (16) The three-dimensional speed of sound is found from (21) by On comparing (102) and (103) we see that similar to an expression for h derived for the accretion onto a black hole in a string cloud background [59]. Using (44) in (102), we obtain where Inserting (105) where m 2 has been absorbed into a re-definition of (t, H). A couple of remarks concerning the fluid flow onto an anti-de Sitter-like f(R) black hole are in order. For ordinary matter K > 0 and f c,r c > 0 (since we are interested in the cases where For the case (a) the sum of the terms inside the square parentheses in (107) is positive while the coefficient . So, the Hamiltonian too diverges. Since the latter has to remain constant on a solution curve, we conclude that there are no global solutions in this case (solutions that extend to spatial infinity). This conclusion remains true even if Λ = 0 provided β = 0. If Λ = 0 and β = 0 (the Schwarzschild metric), the global solutions do not exist if |v ∞ | = 1 (62) and exist otherwise provided 0 < α ≤ 2 if |v ∞ | = 0 or 0 < α if 0 < |v ∞ | < 1.
For the case (b), since Y < 0, we can make it such that in order to have global solutions. For instance, if we restrict ourselves to v having an expansion in powers of 1/r with a vanishing three-dimensional speed at spatial then, on observing (108) we find α = 3, δ ≥ 4, and This is another, rather much harder, fine tuning problem. Here Y depends on n c , so is v 1 : Unless v 2 1 is the rhs of (110), there will be no global solutions to this case too.
For non-ordinary matter, since K < 0, the above two cases are reversed, that is, for γ > 1 it is possible to have global solutions, again with a fine tuning problem, while for γ < 1 (γ = 0) there are non global solutions.
In the following we provide two curve solutions for an anti-de Sitter-like f(R) black hole in the cases γ > 1 (non-global solution) and γ < 1 (global solution) and a curve solution for a de Sitter-like f(R) black hole in the case γ > 1. First, using (44) Since at the CPs we have a 2 c = v 2 c (41), we replace a 2 in (111) and in (42) by v 2 c and solve the system (111) and (42) to find the CPs (r c , v c ). We rewrite these latter equations after making the substitution a 2 c = v 2 c as Here we keep using f to show the general character of these equations. Inserting (113) into (112) we can first solve numerically for r c then get v c from (113). Since the signs of both sides of (112) must be the same, we conclude that, for γ < 1, v 2 c > γ − 1 (which is always satisfied) and that, for γ > 1, v 2 c < γ − 1. Notice that the solution curves do not cross the r axis at points where v = 0 and r = r h , for otherwise the Hamiltonian (107) would diverge there. We recall that r h is the unique horizon of an anti de Sitter-like f(R) black hole or it represents either the event horizon r eh or the cosmological horizon r ch of a de Sitter-like f(R) black hole. The curves may cross the r axis at the unique point r = r h in the vicinity of which v behaves as if f = 0 has a single root at r h . We see that only solutions with 1 < γ < 2 may cross the r axis. Here H(r h , 0) is the value of the Hamiltonian on the solution curve, which is the limit of H(r, v) as (r, v) → (r h , 0). This can be evaluated at any other point on the curve. The pressure p = Kn γ diverges at the horizon as     For an electromagnetic source, The trace of (118) yields reducing (118) to where G µν is the Einstein tensor. On comparing (120) with the field equations of general relativity, we see that R 0 /4 plays the role of an effective cosmological constant and T µν /[1 + f ′ (R 0 )] is an effective SET. If the vector potential A µ = (−Q/r, 0, 0, 0), we obtain the spherically symmetric solution given by (1) with 5 For n > 0, Eq. (119) has always the root R 0 = 0. Notice that the model (122) has been introduced in order to keep |f ′ (R 0 )| ≪ 1, which ensures stability. Hence, we rule out the case 0 < n < 1 which would yield |f ′ (R 0 )| → ∞ as R 0 → 0. For n ≥ 1, the root R 0 = 0 reduces (121) to Reissner-Nordström black hole.
From now on we take n = 2. Since we want that one of the other roots of (119) be R 0 = q 1 M 2 , we substitute (123) and (124) into (119) to obtain yielding With the numerical values in (123) and (124), the four values of c 1 and c 2 are all negative and one should keep those values that ensure |f ′ (R 0 )| ≪ 1 With f (r) given by (121), the rhs of (113) reads For the plots of Fig. 7, we used Eqs. (112) and (128) to find the critical points. The graphs show that accretion is insensitive to the values of the constants (c 1 , c 2 ) and to the value of f ′ (R 0 ) whose effect is to modify the value of the charge in (121).

VIII. CONCLUSION
We have developed a Hamiltonian dynamic system for tackling a variety of problems ranging from accretions, matter jets, particle emissions to cosmological and astrophysical applications whenever conservation laws apply. There are several choices for the dynamical variables arguments of the Hamiltonian. The advantage of using the three velocity is that this entity is bounded (by −1 and 1) and it does not diverge in contrast with the pressure and the baryon number density, and other densities, which may diverge on the horizons. Throughout the paper we kept using the metric coefficient f (r) to emphasize the general character of the derived mathematical expressions. Since the scope of the model of accretion is fairly wide and applies to all static spherically symmetric solutions (asymptotically flat or else), the present analysis can also be done for other f(R) black holes as well as f(T) black holes [56]. Due to the generality of our work, further analysis will be trivial.
Our general results that applies to all metrics of the form (1) and to all perfect fluids, independently of the form of the EOS, are as follows. The Michel-type accretion of a perfect fluid is characterized by • The thermodynamic state functions are determined upon integrating a first order differential equation; • If the three velocity vanishes on the horizon(s), the particle number density n diverges there independently of the expression of f and of that of the EOS. Since the specific enthalpy h is never zero for ordinary matter, this implies that the sum e + p diverges there at least as fast as n; • The fluid may become ultra-stiff as it approaches the horizon(s).
By applying the Hamiltonian dynamic system to f(R) gravity we have performed a detailed analysis of the Michel-type accretion onto a static spherically symmetric black hole in f(R) gravity. Not every model of f(R) theory can predict black holes unless the function f(R) satisfies certain viability conditions such as f ′ (R) > 0 and f ′′ (R) > 0, and asymptotically de Sitter phase at present time (see further details in [93]).
To understand the nature of the f(R) black hole and to distinguish it from the known General Relativity black holes, it is worthwhile to study their astrophysical features such as the accretion of various kinds of fluids and their dynamics near them. Using the isothermal and polytropic equations of state, we showed that the standard method employed for tackling the accretion problem has masked some important properties of the fluid flow.
Accretion of isothermal perfect fluids is is characterized by • Existence of subsonic flows for all values of the radial coordinate. These solutions represent neither transonic nor supersonic flows as the fluid approaches the horizon; • Existence of solutions with vanishing three velocity as the fluid approaches the horizon. As v → 0, the fluid cumulates near the horizon resulting in a divergent pressure which pushes the fluid backward (flowout or a wind of the fluid under the effect of its own divergent pressure). These solutions, as the one depicted in Fig. 3, exist even in the case of a Schwarzschild black hole; • If the CP is a saddle point, the critical solution curve divides the (r, v) plane into regions where the flow is physical in some of them (corresponding to higher values of the Hamiltonian) and unphysical in the others (corresponding to lower values of the Hamiltonian); • The existence of separatrix heteroclinic orbits is subject to no constraint. We have checked this conclusion for the f(R) model of Ref. [51] and for Schwarzschild black hole and this should apply to all black holes; • For the f(R) model of Ref. [51], the existence of two CPs (one saddle and one center), with a possibly periodic flow inside a finite region of space, constraints the values of β not to exceed some lower limit; • Instability of the critical flow.
The polytropic test fluid has nearly no global solutions for the f(R) model of Ref. [51] unless one can deal with the fine tuning problem consisting in fixing the speed at spatial infinity in terms of the number density. Among the solutions we derived for the polytropic test fluid no saddle CP occurs. Moreover, the subsonic flow appears to be almost non-relativistic. This features appear quite different from the General Relativity black holes [90].
de Sitter-like f(R) black holes are characterized by the presence of closed, but non-homoclinic orbits, joining the event horizon to the cosmological horizon. Such cyclic curves are maintained by the high pressure present in the vicinity of the two horizons and do not require the presence of source-sink system for their realization. For γ > 3/2, the proper period of the cyclic flow converges to a finite value and has a logarithmically divergent limit for γ = 3/2. Comparison of the solutions (Figs. 6 and 7) show that the accretion is insensitive to the f(R) model. For the other sonic point, f c = 0 and a 2 c = 1, the rhs of (B.5) is manifestly zero. The rhs of (B.6) is also zero by (25) and (41). The latter provides the value of f c,r c as the limit r c → r f and a 2 c → 1.