A&A 421, 1-8 (2004)
S. K. Chakrabarti1,2 - K. Acharyya2 - D. Molteni3
1 - S. N. Bose National Centre for Basic Sciences, Salt Lake, Kolkata 700098, India
2 - Centre for Space Physics, Chalantika 43, Garia Station Rd., Kolkata, 700084, India
3 - Dipartimento di Fisica e Tecnologie Relative, Viale delle Scienze, 90128 Palermo, Italy
Received 16 October 2003 / Accepted 16 January 2004
We present the results of several numerical simulations of two dimensional axi-symmetric accretion flows around black holes using Smoothed Particle Hydrodynamics (SPH) in the presence of cooling effects. We consider both stellar black holes and super-massive black holes. We observe that due to both radial and vertical oscillation of shock waves in the accretion flow, the luminosity and average thermal energy content of the inner disk exhibit very interesting behaviour. When power density spectra are taken, quasi-periodic variabilities are seen at a few Hz and also occasionally at hundreds of Hz for stellar black holes. For super-massive black holes, the time scale of the oscillations ranges from hours to weeks. The power density spectra have a flat top behavior with average rms amplitude of a few percent and a broken power-law behavior. The break frequency is generally found to be close to the frequency where the shock oscillates.
Key words: methods: numerical - black hole physics - hydrodynamics - instabilities - shock waves
X-rays emitted by galactic black hole candidates are known to exhibit quasi-periodic variations (QPVs). Today, there are almost a dozen confirmed stellar mass black hole candidates for which QPVs are regularly observed, some at normal frequencies (0.1-10 Hz) and some with frequencies of hundreds of Hz. While it is conceptually easier to explain such variations when there is a hard surface, as on a neutron star (e.g., Shaham 1987), in black hole candidates they are difficult to understand as there are no surfaces from which a perturbation could be reflected back so as to produce a sustained oscillation of significant amplitude.
In a very important work, Langer et al. (1981) first suggested that in the presence of cooling effects standing shocks on a white dwarf surface can oscillate and in fact the oscillation may also subside when cooling is increased (Chanmugam et al. 1985). While a black hole does not have a hard surface like neutron stars or white dwarfs, the centrifugal force in accretion matter close to a black hole can become strong enough to actually slow down matter and form centrifugal-pressure-supported axisymmetric standing shocks (Chakrabarti 1989) depending on specific energy and angular momentum (this region will be referred to as CENBOL, centrifugal pressure dominated boundary layer). In the particle dynamics picture, the appearance of a discontinuity can be argued this way: the centrifugal force grows more rapidly (for constant or almost constant specific angular momentum ) than the gravitational force and they become comparable at a distance of , where M is the mass of the central compact object. In fluid dynamics, because of the addition of thermal pressure and/or radiation pressure, such a balance takes place farther out where a shock may develop. The almost constancy of angular momentum comes into being primarily because, as has been shown by extensive work done earlier (Chakrabarti 1996a), the viscosity timescale is much longer than the infall time scale r3/2/(GM)1/2, being the kinematic viscosity. It was shown by Molteni et al. (1996, hereafter referred to as MSC96) and Ryu et al. (1997, hereafter RCM 1997) that these shocks may also oscillate in the presence of cooling. These oscillations are explained to be due to a resonance between the dominant cooling time scale and the infall time scale, and the cooling could be due to thermal/non-thermal radiative effects (as in MSC96), or to dynamical cooling due to outflows (as in RCM97) or both. This shock forms in low-angular-momentum (in the sense mentioned), quasi-spherical flows which may or may not surround a Keplerian disk located in the equatorial plane. In case the latter disk is present and the soft photons emitted are intercepted, inverse-Comptonized and re-emitted by the post-shock region, as in the two-component advective flow model (Chakrabarti & Titarchuk 1995; hereafter CT95), then, as the shock oscillates, the intensity of hard X-rays is also modulated at the same frequency.
If our shock oscillation solution is correct, then a large number of simple predictions could be made: (a) the harder radiation coming from the hot, post-shock region would primarily participate in QPVs while the softer radiation would have less power. This was verified by Chakrabarti & Manickam (2000, hereafter CM00) and Rao et al. (2000). (b) The stronger a shock is, the stronger is the thermally driven outflow rate. Thus, the outflow from the post-shock region depends on the shock strength, i.e., the spectral states (Chakrabarti 1999). This has also been tested by various observations (e.g., Dhawan et al. 2000). (c) In case of outflows, the electron density in the post-shock region is reduced and the Comptonized spectrum should be softened. Chakrabarti et al. (2003) showed that indeed the presence/absence of outflows does change the spectral slope in the hard X-ray region. These observations vindicate that the Comptonized photon, and the major part of the outflows both come out of the post-shock region. Thus, a thorough study of the behavior of CENBOL, especially when it is time dependent, is needed.
In the present paper, we carry out numerical simulation of thick,
inviscid, rotating axisymmetric flows, in the presence of
cooling effects for a large region of the parameter space.
The hot electrons in the CENBOL
reprocess soft photons and generally enhance their
energy before these photons escape. The degree by which the
energy is transferred depends on the so-called
Compton y parameter (Rybicki & Lightman 1979),
Unlike MSC96, where symmetry around the equatorial plane was assumed (and thus the simulation was carried out only in one quadrant), we consider both halves above and below the equatorial plane. In this way, allowance is made not only for radial oscillations of the shocks but also for vertical oscillations. Earlier, in a simulation of a non-shock flow it was observed that the outgoing winds could interact with the inflow and the interface might undulate (Molteni et al. 2001). So we expect that shocks should also have vertical oscillations in presence of winds. In the present simulations, we observe that the shocks do oscillate both vertically and horizontally for a range of inflow parameters, and the PDS of the light curve obtained by our numerical simulation does have characteristics similar to the PDS of the observed light curves of the black hole candidates. The PDS of a numerically simulated light curve has a flat top up to a break frequency, and an rms amplitude in the range of a few percent depending on the accretion rate of the flow. Considering that an axisymmetric shock fundamentally separates a flow into a fast moving (lower density) supersonic and a slow moving (higher density) sub-sonic region, it should not be surprising that its oscillation frequency will generally lie at or near the break frequency in the power law. Of course, shocks are not infinitesimally thin. Neither is it true that they exhibit only a single mode of oscillation. Thus broad and often multiple peaks are expected. This is precisely what we see. What is more, since the inner sonic point also separates the flow into two parts, namely, one with the disk-like behavior and the one with free-fall behavior, there is often a weak peak in the PDS at some hundreds of Hz, perhaps corresponding to an oscillation of this inner region. We shall show these as well.
In the next section, we present the basic setup of our numerical simulation. In Sect. 3 we present some results of the simulation for galactic and extra-galactic black holes and analyze them to demonstrate shock oscillations. From the PDS of the "simulated'' light curves, we show that there are sharp peaks at certain frequencies only. The PDS seems to have the familiar flat top shape with a break. In Sect. 4, we discuss possible cause of the behaviour of the shocks. Finally, in Sect. 5 we draw our conclusions.
We are interested in simulating the behavior of an axisymmetric, inviscid flow in the presence of cooling effects. Unlike an ordinary star, a black hole is capable of accreting matter even in the absence of viscosity when the pressure is large enough. However, if the specific angular momentum is less than the marginally stable value ( 3.67 GM1/c for a Schwarzschild black hole of mass M1), even a cool, non-viscid gas can enter a black hole without any problem. We chose the specific angular momentum of the flow to be lower compared to this value everywhere. This makes , the Keplerian value, everywhere in the flow. We call such a flow a "sub-Keplerian'' flow. The motivation for studying such a flow stems from the fact that accretion processes into a black hole are necessarily transonic, and the flow has to be sub-Keplerian near the horizon. Furthermore, a flow can have a substantial amount of sub-Keplerian matter itself because of accretion of winds from the companion. The wind of velocity , having "zero'' angular momentum with respect to the mass-lossing star of mass M2 and radius R2 has an angular momentum of , where cm is the radius of the accretion disk. For typical values: 1033 gm, cm and we get, in units of 2GM1/c, with . Hence, our choice of 1.75 which is averaged over the entire flow (thus could be much smaller than ) is justified. Recent observations do indicate the presence of both Keplerian and sub-Keplerian components (Smith et al. 2002) in half a dozen black hole components.
It was shown by Chakrabarti (1989, 1996a) that in the large region of the parameter space spanned by the specific angular momentum and energy a stable standing shock may form in accretion whenever the Rankine-Hugoniot conditions are satisfied. Here, the flow first passes through the outer saddle-type sonic point, then through a standing shock and finally through the inner saddle-type sonic point. It was subsequently shown by MSC96 and Ryu et al. (1997) that even when the parameters of the injected flow are outside this region oscillating shocks may form. In the present paper we use Smoothed Particle Hydrodynamics (SPH) to do the simulations as in MSC96, but we use both halves of the flow to inject matter. The code has already been tested for its accuracy against theoretical solutions already (Chakrabarti & Molteni 1993, MSC96) and we do not repeat this here. We just want to recall that the "pseudo''-particles have been chosen to be "toroidal'' in shape and they preserve angular momentum very accurately (Molteni et al. 1996). In MSC96, suggestion was made that radial oscillations of the shock may take place when the infall time in the post-shock region is comparable to the cooling time scale. Presently, we concentrate on the vast region of the parameter space where both the radial and vertical oscillations are present.
Table 1: Inputs and extracted parameters for the model runs.
In the context of galactic black holes one class of QPVs have been observed that are commonly known as Quasi-Periodic Oscillations or QPOs. Our solutions for QPVs are so generic that we believe that QPVs should be observed even for super-massive black holes. Actually there are indications that QPOs in supermassive black holes may have been observed (Halpern & Marshall 1996). Supermassive black holes are believed to be fueled by low angular momentum matter generated from winds of nearby stars. So, a similar consideration as that for wind accreting galactic black holes will apply. As far as Comptonization is concerned, we use the following simple procedure: (a) we use an accretion rate and run the code with bremsstrahlung cooling alone. (b) We get the temperature of the electrons and the optical depth of the "mean'' size of the CENBOL region ( in Table 1). To get the corresponding physical parameters (as if the Compton cooling were included to begin with) (c) we assume that the injected soft photons have energy 1 kev for the stellar mass black holes and 0.01 kev for the massive black holes and using that (d) we compute the enhancement factor from DLC91 for uniform spherical geometry for this injected photon and the electron temperature (assumed to be lower by a factor of than the proton temperature (e.g., Rees 1984)) obtained in each run. (e) We reduce the accretion rate by this to get and the optical depth by . For self-consistency, we ensure that thus obtained is the same as that used while computing from DLC91. We also compute the Compton y parameter to check if it has an acceptable value. This procedure re-formulates the problem of inclusion of the exact cooling processes into the case where a simplified cooling effect may be used. Since we are interested in the integrated energy loss due to cooling and not the spectral properties, this procedure should give reasonable answer.
Table 1 lists the parameters for the Model runs we report in this paper. There are basically two groups of input. In Group A, we consider the black hole to be super-massive ( ) and in Group B, we consider a stellar mass black hole ( ). Different runs are characterized by the injected matter density (gm/cc) at the outer boundary of the numerical grid given in the second column. In all the model runs, we choose the following parameters: the index in the polytropic relation (P is the isotropic pressure and is the gas density) is 5/3, the outer boundary , the specific angular momentum (in units of 2GM/c), injected radial velocity , the sound velocity a=0.04 and the vertical height h=15. The vertical height at the outer boundary is so chosen that the flow remains in hydrostatic equilibrium in the vertical direction at the point of injection . We measure all the distances in units of the Schwarzschild radius of the black hole , all the velocities in units of the velocity of light c and all the masses in units of , the mass of the black hole. c and G are the velocity of light and the gravitational constant respectively. Note that the marginally stable (lowest possible angular momentum with a stable Keplerian orbit) angular momentum is 1.83 in these units. In the third column we give the enhancement factor that we obtained using DLC91 for our run and in the fourth column we provide in units of , the Eddington rate, the effective accretion rate in the presence of Comptonization. In the fifth column we give the Compton y parameters. We assume soft injected photon of energy keV for galactic black holes, and 0.01 keV for supermassive black holes in order to calculate . In the sixth column, we give the average observed location of the shock. In the seventh column, we give the average number of "pseudo''-particles present in the disk in each run. In the eight column we present the QPV(QPO) frequency in Hz as obtained from PDS. In the ninth column, the inverse of the infall timescale corrected for a constant compression ratio R (4 for a strong shock) is provided as an indication of the expected oscillation time scale. In the tenth column, the ratio of the quantities in the eighth and ninth column is given. In the eleventh column, the rms amplitude in percentage is given, we shall discuss in the next section. As in MSC96, we choose the Paczynski-Wiita (1980) pseudo-Newtonian potential to describe the space-time around the Schwarzschild black hole. We believe that the result would remain unchanged if a true Schwarzschild geometry were used, since the transonic properties do not change very much. However, if a Kerr geometry were used, we expect the locations of the inner sonic points and shocks to move inward (Chakrabarti 1996b) giving rise to higher QPO frequencies.
One of the important conditions of a standing shock
is that the thermal pressure (P) plus ram pressure (
continuous across the shock in the steady state, i.e.,
In Group B Runs, within our parameter range, occasionally there were stronger shocks which oscillated radially, and the interaction between the outflow and the inflow produces vertical oscillations as well. Strong convection and turbulences were seen in the post-shock region which contributed to pushing the shock upstream. As a result, the shock location was not seen to change monotonically with accretion rate. The location is clearly determined by several effects: centrifugal force, thermal pressure, ram pressure, turbulent pressure and the cooling rate. The nature of the dependence has been discussed in earlier work in detail (Chakrabarti 1989; Chakrabarti & Molteni 1995; Molteni et al. 1994; MSC96) and we do not repeat this here. We will just recall that the shock location is independent of the location of the numerical boundary when the fundamental quantities such as specific angular momentum, accretion rate and specific energy at a given point (say, at the inner sonic point) are kept fixed. Cases B.1 and B.2 have very small Compton parameters and negligible Comptonization. Cases B.3 and B.4 are presented for academic purposes, just to indicate the effects of excess cooling. The Compton parameters are very large in these two cases.
For reference, we may add that the light crossing times of the horizon for the two classes of black hole are s and 10-4 s respectively. Since a steady shock (in the Schwarzschild geometry) may form anywhere between to , the QPO time-period, if assumed to be related to the infall time , may be close to 105 to 108 s for super-massive black holes, and close to 10-2 s to 10 s for stellar mass black holes. Since we chose the outer boundary of the simulation to be at , the infall time is roughly in units of . In order to trust the results of our simulations, we ran the code for a duration several hundred times of this timescale (typically, or more). For a black hole of , this time corresponds to only 5-6 s of real time.
First, we show examples of shock formation which exhibit radial and/or vertical oscillations. Figures 1a-c shows the locations of the SPH particles (dots) along with the velocity vectors (arrows) of every fifth particle for clarity. This case corresponds to case A.2. The matter distribution is shown at three different times (in units of ) illustrating the vertical and radial motions of the shock which oscillates between 11 to 15 Schwarzschild radii, mostly staying close to . Some matter can be seen bouncing back from the centrifugal barrier (boundary between the empty region along the vertical axis and the disk/jet matter) near the axis and forming giant vortices which interact with the inflow. These vortices push the post-shock region ( ) alternately into the upper and the lower halves.
|Figure 1: Snapshots of simulations of accretion disks around a black hole by Smoothed Particle Hydrodynamics. The dots are particle locations and arrows are drawn for every fifth particle for clarity. Time (in units of ) is marked in each box. Note the vertical as well as radial oscillation of the accretion shock wave located at .|
|Open with DEXTER|
In Figs. 2a-d we show the variation of the light curve in runs A.1-A.4 (marked by a-d in the figure). The light curves are the variations of the luminosity of bremsstrahlung radiation (in arbitrary units) for all the matter within r=50 flow as a function of time (in seconds). As the density is increased, the average thermal energy gradually goes down due to the presence of enhanced cooling and thus the luminosity also goes down. The average location of the shock decreases (MSC96, Table 1). The average number of pseudo-particles (Table 1) as well as the variation of the number of these particles also go down. As a result, with the increase in cooling loss, the light curves are found to be less noisy, and the amplitude of fluctuation is found to be lower.
|Figure 2: a)- d): Total bremsstrahlung luminosity of the accretion flow (in arbitrary units) as a function of time (in seconds). The injected densities are a) 10-14 gm/s; b) gm/s; c) gm/s and d) gm/s respectively. The disk becomes cooler with increasing accretion rate due to bremsstrahlung energy loss.|
|Open with DEXTER|
|Figure 3: a)- d): Power Density Spectra (PDS) of the four cases shown in Fig. 2. Quasi Periodic Variation frequencies (denoted by ) can be seen with timescales of hours to weeks. QPV peaks (marked with arrows) are located near break frequencies with flat top behavior at low frequency and power-law behavior at high frequency. See Table 1 for other properties.|
|Open with DEXTER|
In Figs. 3a-d we show the results of the Fast Fourier Transform of these light curves using FTOOLS provided by NASA, the same software package that is used to analyze observational results. We find the remarkable result of familiar power density spectra (PDS) with a QPV near the break frequency. The QPV frequencies, indicated by the arrows, are presented in Table 1. We find that, perhaps due to the presence of both the vertical and radial oscillations, there are two QPVs; the frequency of the stronger one gradually increases with the cooling, while that of the weaker one gradually decreases. The occurrence of QPV frequencies at or near the break frequency of the PDS can be understood by the following: according to our solution, QPVs occur due to oscillations of the shock waves that separate two "phases'' of the flow - the pre-shock flow is supersonic and weakly radiating, while the post-shock flow is subsonic and strongly radiating. The former region of the flow produces the "flat-top'' region in the PDS, while the latter region produces a PDS with a slope -2. The rms amplitude of the flat top region up to the break frequency is shown in Table 1. It has generally similar value at around 11.2% except in (d) where it is only around 3.3%. These rms values are similar to the observed results from black hole candidates as well. In Figs. 3a and b, where the cooling is weaker, the QPVs are less prominent and we had to search the locally significant peaks of the PDS for them. As the cooling becomes stronger, the peak becomes stronger.
In Figs. 4a-d, we present similar results as in Fig. 2 for stellar mass black holes (runs B.1-B.4 of Table 1). Here, too, the disk becomes cooler with increasing density and the light curve changes its character according to whether the shocks are strong or not. In the beginning of the simulation some transient behaviour is seen till the disk settles down within a second in real time. In Fig. 4c, the parameters are such that a strong shock forms and it heats up the disk while for the parameters in Fig. 4d a very weak shock forms. Shock oscillation produces a large amplitude noise in the B.3 case (Fig. 4c). In Figs. 5a-d, the PDS of these four runs are shown. The transient region has been excluded while doing the FFT. As in Figs. 3a-d, here too QPVs are observed. For small black holes QPVs occur at a few Hz and the QPV is located at or near the break frequency. Thus the characteristics are similar to those of QPOs observed from black hole candidates. Radial and vertical oscillations often produce multiple peaks in this case as well. Other details can be seen in Table 1. In Fig. 5c the peak is not very strong, perhaps due to the large amplitude noise that is produced in this case (See, Fig. 4c).
|Figure 4: a)- d): Total bremsstrahlung luminosity of the accretion flow (in arbitrary units) as a function of time (in seconds). The injected densities are a) 10-10 gm/s; b) gm/s; c) gm/s and d) gm/s respectively. The disk becomes cooler with increasing accretion rate due to bremsstrahlung energy loss.|
|Open with DEXTER|
Occasionally, peaks at a very high frequency could be seen, but these peaks are found to be transient. If the entire light curve is broken into smaller pieces, weak peaks at frequencies 100-300 Hz can be seen, though they do not persist. One example of this is shown in Fig. 6a where the PDS of the average thermal energy up to T=36 000 is plotted for the case B.4. In Figs. 6b and c we show the flow pattern at two times which are separated by only 18 units (0.0018 s). The shock located at shows two distinctly different shapes at these two times, bent up and bent down. Oscillation of a density enhancement at would have a time period of s, where R, the compression ratio, is 2 for a weak shock. The corresponding frequency is 333 Hz. The weak peak at Hz (Fig. 6a) may be due to this oscillation. When data for much longer time is analyzed, this weak peak disappears completely.
|Figure 5: a)- d): Power Density Spectra (PDS) of the four cases shown in Fig. 4. Quasi Periodic Variation frequencies (denoted by ) can be seen in the range of 2-10 Hz and also close to several hundred Hz. QPV peaks (marked with arrows) are located near break frequencies with flat top behavior at low frequency and power-law behavior at high frequency. See Table 1 for other properties.|
|Open with DEXTER|
|Figure 6: a)- c): a) Power density spectrum of the average thermal energy for case B.4 with time up to 3.6 s. Snapshots of simulations at closely separated times (marked) are shown in b) and c). Only the inner 10 Schwarzschild radii are shown. The dots are particle locations and arrows are drawn for every fifth particle for clarity. Time (in units of ) is marked in each box. The oscillation of this region causes QPV at 300 Hz (marked with an arrow in a)).|
|Open with DEXTER|
A shock remains steady, or "standing'' if Eq. (2) along with the specific energy and angular momentum flux remain continuous across it. These conditions are known as the Rankine-Hugoniot conditions (Landau & Lifshitz 1959). However, in the presence of an efficient energy extraction mechanism from the disk, whether through thermal/non-thermal radiation or through higher entropy flows along the shocked winds and outflows, the shock should move closer to the black hole as the pressure in the post-shock region (P+) drops. This is what is seen (see, Table 1), except in case B.3 where the turbulence was too strong to push the shock outward.
What about the oscillation time scale? We note that the cooling time scale (where is the specific energy and a "dot'' indicates its time derivative) competes with the infall time scale . MSC96 argued that when these two time scales agree roughly, i.e., roughly when the resonance takes place, the shocks may start oscillating. The infall time is not easy to compute, however. CM00 argued that the infall velocity need not be free-fall velocity as the post-shock flow is reduced at the shock by at least a factor of R, the compression ratio. Furthermore, the velocity is highly modified due to turbulent pressure and angular motion and could be very slowly varying in between the shock and the inner sonic point (CM00; Chakrabarti et al., in preparation). A third and perhaps the most important parameter in predicting the time scale is the rate of outflow from the post-shock region and the extent to which the total mass of the disk oscillates. These two cooling mechanisms, namely, through radiation loss and through outflow, are acting in opposite directions - enhanced thermal cooling reduces the pressure gradient force, thereby reducing the wind rate. The QPV timescale is expected to be close to . In Table 1, is presented. This ratio seems to be closer to 0.5 (justifying a near resonance condition) except when the cooling is very strong (runs A.4 and B.4) when the ratio drastically goes down. In these latter types of runs the shock moves very close to the black hole and becomes very weak. The QPV could be due to the domination of the interaction of the winds with the inflow. The two QPV frequencies often seen could be due to the two types of cooling.
In this paper, we presented the discovery of Quasi-Periodic Variations of emitted radiation through extensive time-dependent numerical simulations. In earlier studies, such as in MSC96 and CM00, it was suggested that shock oscillations might contribute to the observed QPOs. In the present work, we studied shock oscillations in low-angular-momentum axisymmetric flows around stellar and super-massive black holes for a range of accretion rates. Unlike MSC96, we used fully two dimensional inflow (upper and lower halves) and demonstrated that the shocks participate in vertical and radial oscillations. Using PDS of light curves obtained from time dependent numerical simulations of sub-Keplerian flows, we demonstrated that the QPVs could be the cause of the observed QPOs. We showed that the computed light curve produces power density spectra with characteristic QPVs located at or near the break frequency which separates the flat-top type and power-law type PDSs. For stellar black holes, we see QPVs at reasonable frequencies (2-20 Hz). Very often high frequency QPVs (at around 100-300 Hz) could be observed. We demonstrated that this could be due to oscillation at the very inner edge where the flow becomes supersonic. For Active Galaxies and Quasars, QPVs at very low frequencies ( 10-7-10-6 Hz) are observed. Their properties are similar to those of QPOs observed in these candidates.
Based on our simulations we propose that the QPOs observed in black hole candidates are genuinely related to the accretion shock oscillations. Variation of QPO frequencies from time to time in the same object are thus manifestations of the variation of the Keplerian and sub-Keplerian accretion rates. In our work, we assumed a simple Paczynski-Wiita potential to study the shock properties around Schwarzschild black holes. Based on our experience with transonic flows, we are confident that none of the properties would be affected if the simulations were carried out in true Schwarzschild geometry. However, if a rotating black hole is used, the shocks will form much closer to the black hole (Chakrabarti 1996b), resulting in higher QPV frequencies.
S.K.C. thanks the support of Indian Space Research Organization for a RESPOND project. Discussions with Prof. A. R. Rao (TIFR) and Mr. A. Nandi (SNBNCBS) are also acknowledged.