Bondi-Hoyle Accretion around the Non-rotating Black Hole in 4D Einstein-Gauss-Bonnet Gravity

In this paper, the numerical investigation of a Bondi-Hoyle accretion around a non-rotating black hole in a novel four dimensional Einstein-Gauss-Bonnet gravity is investigated by solving the general relativistic hydrodynamical equations using the high resolution shock capturing scheme. For this purpose, the accreated matter from the wind-accreating X-ray binaries falls towards the black hole from the far upstream side of the domain, suphersonically. We study the effects of Gauss-Bonnet coupling constant $\alpha$ in 4D EGB gravity on the accreated matter and shock cones created in the downstream region in detail. The required time having the shock cone in downstream region is getting smaller for alpha>0 while it is increasing for alpha<0. It is found that increases in alpha leads violent oscillations inside the shock cone and increases the accretion efficiency. The violent oscillations would cause increase in the energy flux, temperature, and spectrum of X- rays. So the quasi-periodic oscillations (QPOs) are naturally produced inside the shock cone when -5 \leq alpha \leq 0.8. It is also confirmed that EGB black hole solution converges to the Schwarzschild one in general relativity when alpha \rightarrow 0. Besides, the negative coupling constants also give reasonable physical solutions and increase of alpha in negative directions suppresses the possible oscillation observed in the shock cone.


INTRODUCTION
Einstein's general theory of relativity has successfully passed many observational tests during the last decades (Akiyama1 et. al. 2019;Akiyama2 et. al. 2019;Akiyama3 et. al. 2019;Abbott1 et. al. 2016;Abbott2 et. al. 2016). The existence of super-massive black hole at the center of galaxy M87 was one of the most impressive discoveries. The observed M87 black hole shadow indicates that the observed data well consistent with the prediction of the general theory of relativity. The observed gravitational waves created during the coalescing of two stellar mass black holes made by LIGO detectors was another discovery to see the ripples of the space time predicted by Einstein's general theory of relativity. The observed electromagnetic spectrum emitted from an accretion disk around the black holes could allow us to define many properties of the black holes, such as mass, spin, etc. (Frank et. al. 2002;Yuan & Narayan 2014;Nampalliwar & Bambi 2019). All recent observations indicate that testing the general theory of relativity and its alternatives in a week and a strong gravitational regions are hot topics in this century. On this end, in order to extract more detailed information about the strong gravitational regimes around the black holes, we study the alternative theory to the Einstein's theory of the relativity called four dimensional Einstein-Gauss-Bonnet (EGB) gravity which was first formulated by Glavan & Lin (2020).
A novel four-dimensional EGB gravity theory was proposed in Glavan & Lin (2020) where the Gauss-Bonnet coupling constant ⋆ E-mail: orhan.donmez@aum.edu.kw (OD) is α → α/(D − 4) in the limit D → 4 and bypasses the Lovelocks theorem. As a consequence of this limit, the Gauss-Bonnet coupling constant gives non-trivial contributions to the gravitational dynamics. It is also shown that the theory is free from Ostrogradsky instability and preserves the number of degrees of freedom. And the static spherically symmetric black hole, which has two horizons instead of one compared to the Schwarzschild black hole, was discovered (Cheng et. al. 2020). The horizon is rescaled with the Gauss-Bonnet coupling constant α and the causal structure of the black hole radically is altered with a repulsive effect. Although this 4D EGB theory is currently under debate (Julio et. al. 2020;Gurses et. al. 2020;Ai 2020), it is worthy and meaningful to study the spherical symmetric black hole solution of 4D EGB gravity.
The newly discovered four-dimensional static and spherically symmetric Gauss-Bonnet black hole was used to reveal many interesting features of some astrophysical problems. Lots of work have been done the different astrophysical phenomena the viability of the solution. The properties of the black hole and its shadow (Konoplya & Zinhailo 2020; Roy & Chakrabarti 2020), the innermost stable circular orbit for massive particles and photons Zhang et. al. 2020), generating and radiating black hole (Ghosh & Maharaj 2020;, the gravitational lensing by black holes (Jin et. al. 2020;Islam et. al. 2020), observational constraints on the Gauss-Bonnet constant (Feng et. al. 2020;Clifton et. al. 2020) and Greybody factor and power spectra of the Hawking radiation (Konoplya & Zinhailo 2020;Zhang et. al. 2020) were studied with all details.
The wind accretion scenario onto the black hole is one of the important physical phenomena to explain soft and hard X− rays observed by different X− ray telescopes. The wind accretion process can be explained by using the Bondi-Hoyle model (Edgar 2004). A black hole moving inside the gas cloud captures mass which causes either acceleration or deceleration of it. As a result of this the shock cone or bow shock would be formed around the black hole. This is the one of hot topics studied by different researchers using the analytical and the numerical techniques in the last few decades. Analytic representation of the wind accretion for a perfect fluid onto a Schwarzschild black hole has been found by Tejeda & Aguayo-Ortiz (2019) and numerical treatments from Newtonian and relativistic perspective were extensively studied by Donmez et. al. (2011);Donmez (2012);Cruz-Osorio & Lora-Clavijo (2016);Cruz-Osorio et. al. (2017).
In the present paper, building onto the previous findings, we would like to study the creation of the shock cones and their dynamical evolution in case of Bondi-Hoyle accretion in the background of the 4D EGB spherically symmetric black hole in the strong gravitational region, numerically. Bondi-Hoyle accretion is an important mechanism to explain Quasi-Periodic Oscillations (QPOs) observed by X − ray telescope. The companion star in the X−ray binary system provide the fast wind which causes to accreate towards to the black hole (Orosz et. al. 2008). The accelerating single black holes due to the other black holes could have Bondi-Hoyle accretion (Lora et. al. 2013). Studying the QPOs on accretion disk in vicinity of a black hole is important to explore the physical parameters of the black hole such as its spin and mass. QPO arises in the inner accretion disk of the black hole binary and could be created in terms of an oscillating, precessing hot flow in the truncated-disk geometry due to the strong shocks. Therefore, it would be a further step to know whether the Gauss-Bonnet coupling constant α can play a dominant role in the oscillation of the shock cone created during the Bondi-Hoyle accretion. To reveal all these details, we numerically model the Bondi-Hoyle accretion in the vicinity of 4D EGB black hole. We compute the shock cone structures for different values of Gauss-Bonnet constant α and compare them with standard general relativistic solution, Schwarzschild black hole.
The plan of the paper is as follows: In Section 2, we briefly give summary of the recently proposed non-rotating black hole solution in 4D EGB gravity and definition of the horizon of the black hole. In section 3, we describe the conserved form of the general relativistic hydrodynamical equations, lapse function, and shift vectors in 4D EGB black hole necessary in our numerical simulations. In Section 4, we describe the theoretical framework of the Bondi-Hoyle accretion for pressureless gas, and initial and boundary conditions used in our numerical simulation. The numerical results and discussion are also given to show dependencies of the disk dynamics, creation of shock cones, accretion rates, and mode power to Gauss-Bonnet coupling constant α. In Section 5, we discuss the implications from our numerical results and draw the direction of future work. Unless specified, we use geometrized units throughout the paper for the speed of light and gravitational constant, G = c = 1.

NON-ROTATING BLACK HOLE SOLUTION OF 4D EGB GRAVITY
The theory of the static and spherically symmetric solution of the gravity can be driven by summing Einstein-Hilbert action and higher order Lovelock invariants with the vanishing bare cosmological constant Glavan & Lin (2020).
where MP , R, α, and G are the reduced Planck mass, the Ricci scalar of the space-time, the Gauss-Bonnet coupling constant, and the Gauss-Bonnet invariant, respectively (Glavan & Lin 2020;Cheng et. al. 2020). The theory of this black hole was already established for D 5 (Boulware & S. Deser 2020). But the Gauss-Bonnet term G does not contribute to the 4D gravitational dynamics. The reason is that the equation of a motion contains term D − 4 and it will disappear while D = 4. After rescaling the coupling constant α → α D−4 , the factor D − 4 would be removed. So when the Lovelock theorem is bypassed, it gives nontrivial dynamics, and the static and spherically symmetric black hole solution was discovered (Glavan & Lin 2020).
The static and spherically symmetric black hole solution in four dimensional EGB gravity has the following form where Here M is the mass of the black hole. The coupling constant parameter plays a dominant role defining the black hole horizon. By solving f (r) = 0, we can reach the following two horizons for non-rotating black hole in 4D EGB gravity, where, as seen in Fig.1, the Gauss-Bonnet coupling constant is restricted to α < 1. If α > 0, one has two horizons. One degenerate horizon is seen when α = 0. Otherwise, α < 0, the black hole has only one horizon. It was belied that static and spherically symmetric ansatz in Eq.2 might not give the real solution (Glavan & Lin 2020) for a short radial distance r 3 < −8αM . As shown in Fig.1, the black hole in 4D EGB gravity exists when α < 0. Therefore we would like to numerically model Bondi-Hoyle accretion for possible values of Gauss-Bonnet coupling constant either 0 < α < 1 or α < 0 and to compare them with the theory and Schwarzschild solution. As it was also noted in reference , the Inner Stable Circular Orbit (ISCO) in the novel 4D Gauss-Bonnet gravity can be greater or less than the one in Schwarzschild black hole solution depending on the coupling constant α.  where T µν and J µ are the stress-energy tensor and the matter current density, respectively. The latin indices run from 1 to 3 and Greek indices run from 0 to 3. In order to define the characteristic waves of the GRH equation, we simplify them neglecting the magnetic field and any type of the viscosity. Therefore the stress-energy tensor for a perfect fluid is, where specific enthalpy h = 1+ǫ+ P ρ . ρ, P , u µ , and g µν are the restmass density, ideal gas equation of state (pressure), four-velocity, and four metric, respectively. The equation of state is P = (Γ−1)ρǫ. Here, Γ is adiabatic index which defines the compressibility of the fluid.
The form of GRH equations given in Eq.5 is not suitable to use in High Resolution Shock Capturing scheme (HRSC). Hence, GRH equations are written in conservation form using 3 + 1 formalism and we end up with the following first-order, flux-conservative hyperbolic system, where U , F i , and S are conserved quantities, fluxes, and sources, respectively. Detailed formulations of these three vectors are given as the functions of lapse functionα, shift vectors β i , the determinant of the three-metric γ, four-velocity u µ , three velocity vi, Lorentz factor W =αu 0 = (1 − γijv i v j ) −1/2 , rest-mass density ρ, the components of the momentum vector Si, pressure P , enthalpy h, and the 4D Christoffel symbol Γ α µν . The spectral decomposition of the Jacobian matrix of the system ∂ F i /∂ U is needed to use HRSC scheme and they are given in Donmez (2004) with the full details.
The static and spherically symmetric solution of the black hole in 4D Gauss-Bonnet gravity metric given in Eq.2 is used to define the source term appearing in the right hand side of Eq.7. The fourmetric, three-metric, their inverses, the 4D Christoffel symbol, and Lorentz factor can be easily driven by having straightforward calculations. The lapse function for the EGB black hole isα and the shift vectors are, 4 BONDI-HOYLE ACCRETION ONTO 4D EGB BLACK HOLE

Theoretical Framework
A homogeneous supersonic flow that has a relative velocity v∞ and density ρ∞ at infinity is deflected by point mass M due to the warped space-time (Hoyle & Lyttleton 1939). The Bondi-Hoyle accretion condition for particles is defined with an impact parameter ζ.
As it is seen in Fig.2, if the pressure effect is negligible, the trajectory of the particle can be defined by conventional orbit. And then accretion condition is, where ζHL is commonly known as accretion radius or the stagnation point. The suggested accreated matter could be computed by considering the trapped material inside the gravitational potential. The accretion rate suggested by Hoyle & Lyttleton (1939) is, The material falling towards the black hole that does not encounter to θ = 0 axis would not accreate and it would escape from the system as seen in dashed (blue) lines in Fig.2. The theory given above does not exactly describe the motion of the gas flow accreated onto the black hole in the presence of the gas pressure in the strong gravitation region but it helps us define what the physical parameters are at infinity and how the accretion mechanism works.
Bondi-Hoyle accretion is studied in the vicinity of the Schwarzschild and Kerr black hole and called spherically symmetric accretion. During the accretion, the Bondi radius would be created and the accreated flow outside the radius is subsonic, and the disk mass-density is almost uniform. But the gas inside the radius becomes supersonic and it will be accreated towards to the black hole (Donmez et. al. 2011;Donmez 2012). The strong gravity focuses the material behind the black hole (down-stream region) and the material would be accreated. In this paper, we will numerically model the Bondi-Hoyle accretion in the vicinity of the static and spherically symmetric black hole in 4D Gauss-Bonnet gravity to reveal the effect of Gauss-Bonnet coupling constant, not only on the shock cone structures but also on the oscillation properties and QPOs of those shock cones.

Initial and Boundary Conditions
Numerical simulation of the Bondi-Hoyle accretion around static and spherically symmetric non-rotating black hole defined in 4D EGB gravity gives us an opportunity to understand how the dynamics of the accretion disk and shock cone would be effected by the Gauss-Bonnet coupling constant α. For this purpose, we numerically solve the GRH equations on the equatorial plane (r, φ), i.e. θ = π/2 (see Eq.7) using the HRSC scheme based on approximate Riemann solver, further described in Donmez (2004).
As explained in section 3, we need to know pressure to evolve GRH equation depending on given initial values. The pressure is evaluated by using the perfect fluid equation of states P = (Γ−1)ρǫ with adiabatic index Γ = 4/3. In order to perform the numerical simulation of the Bondi-Hoyle accretion on the equatorial plane, the gas is injected from an outer boundary of the upstream region with the following velocities, V∞ is called the asymptotic velocity of the gas at infinity. The radial and angular velocities of the gas injected from the outer boundary are also given in terms of the components of the three-metric. The Eq.12 guarantees that V 2 = ViV i = V 2 ∞ is valid everywhere along the computational domain.
In order to model the Bondi-Hoyle accretion numerically, a homogeneous supersonic flow is injected from the upstream region [π/2, 3π/2] at the location of the outer boundary using the same analytic prescriptions given in Eq.12. The computational domain is extended from rmin, reported in Table 1 for different models, to outer boundary rmax = 100M which is fixed for all models. φ goes from 0 to 2π. The sound speed and the rest-mass density are chosen as cs,∞ = 0.1 and ρ∞ = 1, respectively at the boundary. And then the gas pressure is computed using the expression p = c 2 s,∞ ρ(Γ−1)/[Γ(Γ−1)−c 2 s,∞ Γ] accordingly (Donmez et. al. 2011). The asymptotic velocity is chosen as V∞ = 0.3 with the parameters mentioned above and given in Table 1.
The initial injected values of a wind flow are used at upstream region at the outer boundary. The continuation is needed to avoid matter reflected back towards the black hole at the outer boundary of the downstream region (from 3π/2 to π/2) along the radial coordinate. It is treated simply by using the zeroth-order extrapolation for all primitive variables. The inner boundary of computational domain along the radial distance close to the black hole horizon is Table 1. The physical parameters and times: Gauss-Bonnet coupling constant α, the inner boundary location of computational domain r in , time to need to create a fully formed shock cone t critial (required time for the steady state), and total simulation time t total . M is the mass of the black hole. implemented by using outflow boundary condition, handled simply copying first values to ghost zones. Lastly, the periodic boundary condition is adopted along the φ direction.

Numerical Results and Discussion
In order to reveal the effects of Gauss-Bonnet coupling constant α on the accreated material and the creation of shock cone, we first focus on the injected matter from the upstream region of the computational domain. Later, we wait some time to shock cone reaching the steady state and it is called t critical . The critical times are varying from 1500M to 1920M , shown in Table.1. The evolution of these models are, at least, followed up to t = 16470M which is sufficient to extract oscillation properties of the shock cone after reaches the critical time. The critical time is getting smaller for α > 0 while it is increasing for α < 0. The morphology of the accreated shock cones formed in downstream region are given for the various values of Gauss-Bonnet coupling constant α, seen in Fig.3 where we report the logarithmic rest-mass density of the accreated matter. The snapshots refer to t ∼ 15000M , which is much later than the critical time needed to reach the steady state. The results due to different α show a qualitatively similar behavior, although quantitative differences appear in all models. The development of the shock cone and its strong shock locations, also seen in Fig.6, could cause a chaotic non-linear phenomena inside the cone. The shock cones are attached to the numerical horizon, rin given in Table 1 for different values of α. It is know that the certain amount of matter would pass across the cone and fall into the black hole.
Besides, understanding the dynamic evolution of an accretion disk due to the Bondi-Hoyle accretion and calculating the mass accretion rates allow us to find out more important features of the shock cone. The mass accretion rate is, where all the physical quantities are defined in Section 3. Here, we report the mass accretion rate at a fixed radial distance r = 6.1M for various values of α. The mass accretion onto the 4D EGB black hole is clearly sensitive to the α, as can be seen in Fig.4. It is clearly seen that increasing the value of α towards zero leads to more violent phenomena inside the shock cone and thus creates severe oscillations after the shock cone reaches the steady state. The  high amplitude oscillation is an indicator of instability fully developed during the evolution. We have also confidence that the inner boundary of the computational domain does not play any role on oscillation properties of the shock cone. The oscillating properties of the mass accretion rates are not the same although the inner boundaries for α = −12 and α = −7 are in the same locations, On the other hand, oscillation inside the shock cone is considerably dissipated by the larger values of the negative α.
After showing the behavior of the accretion rates in vicinity of 4D EGB black hole, we can now discuss the maximum value of oscillation strength using the mass accretion rate. Fig.5 represents the maximum oscillation amplitude (∆Φ) of the shock cone with various Gauss-Bonnet coupling constant α after the shock cone is fully formed. This is clear evidence that the significant oscillations happen when α is getting closer to zero as compared to high negative values. Together with the dependency of the accretion rate to α, we have also investigated the convergence of EGB gravity solution to the general relativistic one. As it can be seen in Fig.5, (∆Φ) around the 4D EGB black hole converges to Schwarzschild one when α → 0.
To have deeper understanding of Fig.3, we extract the information along the angular direction (φ) at a fixed radial coordinate r = 6.5M . The left panel of Fig.6 shows one-dimensional profiles of the rest-mass density for different values of Gauss-Bonnet coupling constant α and Schwarzschild solution. The strong shocks are created at the border of the shock cone. Therefore, a sharp transition is seen in density. The sharp transition location has a lowest density throughout the disk since along this location gas falls towards to the black hole along the streamlines, supersonically. The locations of streamlines seen in left side of the left panel of Fig.6 versus α are given in the right panel of the same figure. Increasing in α (going from −7 to 0.8) causes a slight exponential increase in the location of streamline, seen in the right panel of Fig.6. Hence the shock opening angle appeared in the downstream side of the computational domain expands along the angular direction.
In order to uncover the effects of Gauss-Bonnet coupling constant α on the instability created during evolution and especially inside the shock cone after it reaches the steady state, we numerically study the Fourier mode analyze to compute the saturation point and to analysis the characterization of the instability. The Fourier mode m = 1 is performed and the growth rates are obtained for the the rest-mass density power using the following equation (Donmez 2014), where rout and rin are outer and inner boundary of the computational domain, respectively. The real and imaginary parts are Re(wm(r)) = 2π 0 ρ(r, φ)cos(mφ)dφ and Im(wm(r)) = 2π 0 ρ(r, φ)sin(mφ)dφ. As seen in Fig.7, the instability mode grows in the beginning of the simulation until t ∼ 304M for α = −7 but the time is slightly getting greater when α converges to 0. And then we have witnessed decreasing in the mode for a short time scale (t ∼ 200M ). Later, the m = 1 mode creates an exponential growth again until t ∼ 1400M . The m = 1 mode power reaches a maximum for the α getting close to zero and 0 < α < 1, although the modes in all models saturate almost at the same time t ∼ 1400M . On the other hand, both α close to zero and Schwarzschild black hole power modes exhibit the same behavior during the evolution.
Fourier analysis of the computed data allows us to extract a lot of information about emerged phenomenology and their connections with Gauss-Bonnet coupling constant α. In Fig.8, for instance, we obtain the power spectra using the rest-mass accretion rate computed during the time evolution for different values of α and Schwarzschild black hole. We have also conducted an extra test to find the dependency of power spectra to radial position and found that the mode of oscillation does not depend on the radial position. It means that the modes are global eigenmodes of the oscillating shock cone. As it is seen in Fig.8, there are two genuine eigenmodes appearing in the power spectra and the rest are the results of nonlinear couplings of these genuine modes. For Bondi-Hoyle accretion around the Schwarzschild solution, f1 = 9.5Hz and f2 = 17Hz are genuine eigenmodes, while 20Hz, 27.7Hz, 32Hz, and 45.5Hz are the nonlinear couplings of those genuine eigenmodes with a 2Hz of error bar. These nonlinear couplings are 2f1, 3f1, 2f2, and 2f1 + f2, respectively. Similarly, for α = 0.8, f1 = 6Hz and f2 = 13.5Hz are genuine eigenmodes, while 21Hz, 24Hz, and 46.5Hz are the nonlinear couplings of those genuine eigenmodes with a 2Hz of error bar. These nonlinear couplings are f1 + f2, 4f1, and 5f1 + f2, respectively. For α = −0.8, f1 = 11Hz and f2 = 15.5Hz are genuine eigenmodes, while 36Hz, 42Hz, and 47.5Hz are the nonlinear couplings of those genuine eigenmodes with a 2Hz of error bar. These nonlinear couplings are 3f2 − f1, f1 + 2f2, and 3f2, respectively. For α = −5, f1 = 30Hz, f2 = 39.5Hz, and f3 = 63Hz are genuine eigenmodes, while 49Hz, 58.3Hz, 78Hz, and 89Hz are the nonlinear couplings of those genuine eigenmodes with a 2Hz of error bar. These nonlinear couplings are 2f2 − f1, 3f2 − 2f1, 2f2, and 3f2 − f1, respectively. This behavior is expected behavior of the nonlinear equations in the physical system in the limit of small oscillations (Landau & Lifshitz 1976). It is also important to note that the genuine modes and their nonlinear couplings show different behavior for α = −5. There are three main differences between α close to zero (Schwarzschild solution) and negative value of α when α = −5. First, there are two genuine modes for −0.8 α 0.8 while there are three for α = −5. Second, the frequency of the first genuine mode f1 is getting greater when α goes from 0.8 to −5. Third, the frequencies of these genuine eigenmodes are much higher in α = −5 than the other three cases. As it is expected, the amplitude of genuine modes are getting smaller when α is getting larger in negative direction.

Comparison of 4D EGB and Schwarzschild Black Holes
To illustrate the consistency of the results found from the static spherically symmetric black hole solution in a novel 4D EGB gravity with the Schwarzschild one, the Schwarzschild solution of the Bondi-Hoyle accretion using the same initial conditions is performed. Fig.9 shows the comparison of the accretion rates between EGB and Schwarzschild solutions using certain values of the Gauss-Bonnet coupling constant α. The accretion rates are plotted after the initial phase of relaxation of the shock cone. As it is also seen in the previous section, the computed accretion rate for the α = −0.0001 shows similar behavior (oscillation + amplitude) with the Schwarzschild solution. Obviously, tendency of α to reach zero would give the solution in the general relativity. But decreasing in the parameter α causes a significant change in the oscillation properties of the accreated disk. So that a smaller α would cause to form a less luminous, cooler, and less efficient shock cone around the black hole. This is in agreement with a thin accretion solution around the 4D EGB black hole (Cheng et. al. 2020).
Power spectrum data has a clear signature about the creation of shock cone instability for various values of Gauss-Bonnet coupling constant α and Schwarzschild black hole. Fig.10 represents how the maximum value of the m = 1 power mode (A mod (max)) changes with α and Schwarzschild black hole. A mod (max) increases with an increasing α and it goes to the Schwarzschild solution when α → 0 (e.g. see the star and the square symbols at α ∼ 0 in Fig.10). For the positive values of α, it is seen in Fig.10 that increasing in α produces higher values of the maximum power spectrum.

CONCLUSION
We have performed a systematic investigation of Gauss-Bonnet coupling constant α during the Bondi-Hoyle accretion onto a nonrotating black hole in 4D EGB gravity. The effect of α on the shock cone and its oscillation properties have been extensively studied when −12 α < 1. The numerical simulation of Bondi-accretion onto Schwarzschild black hole using the same initial condition is also performed to compare the general relativistic solution with 4D EGB gravity one. We investigate how the physical features of the shock cone and its oscillation properties depend on α and find the following outcomes summarized in Table 2.
(i) we have found that the Bondi-Hoyle accretion around 4D EGB black hole with a negative Gauss-Bonnet coupling constant produces physical solution even it is very close to the black hole horizon. The same was also confirmed by  for the geodesic motions of time-like and null particles. (ii) The numerical simulations reveal that moving matter from upstream side of the computational domain creates a steady state shock cone in a short time scale. The required time having the shock cone in downstream region is called critical time which decreases for α > 0 while increases for α < 0. (iii) The key indicator of having instability inside the shock cone is a high amplitude oscillation. The increase in the α that is getting close to zero, would lead to more violent phenomena inside the shock cone. Therefore the severe oscillations are created after the shock cone reaches the steady state. In addition, the oscillation inside the shock cone is considerably dissipated by the larger values of the negative α. (iv) The shock opening angle slightly depends on the variation of α.
Going from α = −7 to α = 0.8 would lead an exponential growth in the location of the shock cone.  (v) The power mode m = 1 reaches a saturation point almost at the same time for all models. The m = 1 mode power gets the maximum value when α is close to zero and 0 < α < 1. In addition, the power modes calculations for the Schwarzschild black hole and for α close to zero show the same behavior during the evolution.
In addition to recovering many important features of the shock cones for varying α, the novel features of the shock cones which traps the pressure modes inside the high density are also extracted by using Fourier transform. QPOs are computed for the Schwarzschild and different values of α, i.e. α = 0.8, α = −0.8, and α = −5. The global genuine modes and their nonlinear oscillation counter- Table 2. The effects of Gauss-Bonnet coupling constant α to the accretion mechanism, shock cones, and physical interpretation are summarized. D represents decreasing while I represents increasing. SS, SC, P D stand for Steady State, Shock Cone , and Power Mode, respectively.
There are three genuine modes found in α = −5 while it is two for α = 0.8 and α = −0.8 including the Schwarzschild one. It is also important to note that the higher absolute frequencies are found in case of α = −5. As a result, having three genuine modes need to be confirmed by doing more accurate simulations in the negative α direction. Finally, we have carried out a comparison of 4D EGB solution with Schwarzschild one. Obviously, it is seen from our numerical simulations that the tendency of α to reach zero would yield solution in general relativity.