Evolution of bimodal accretion flows

Abstract

We consider time-dependent models for rotating accretion flows on to black holes, where a transition takes place from an outer cooling-dominated disc to a radiatively inefficient flow in the inner region. In order to allow for a transition of this type we solve the energy equation, both for the gas and for the radiation field, including a radiative cooling flux and a turbulent convective heat flux directed along the negative entropy gradient. The transition region is found to be highly variable, and a corresponding variation is expected for the associated total luminosity. In particular, rapid oscillations of the transition radius are present on a time-scale comparable with the local Keplerian rotation time. These oscillations are accompanied by a quasi-periodic modulation of the mass flux at the outer edge of the advection-dominated accretion flow (ADAF). We speculate about the relevance for the high-frequency quasi-periodic oscillation phenomenon.

accretion, accretion discs, black hole physics, convection, hydrodynamics, instabilities

1 Introduction

The standard model for an optically thick, Keplerian accretion disc (Shakura & Sunyaev 1973) explains the continuous emission observed from a variety of accretion sources both on the galactic and extragalactic scales. This so-called SSD model cannot, however, explain the hard spectral component observed in most compact accreting systems, e.g. in X-ray binaries or in active galactic nuclei (AGNs). This component is generally fitted by a power law and ascribed to a thin, radiatively inefficient plasma. Such conditions are met in an advection-dominated accretion flow (ADAF; for a review see Narayan, Mahadevan & Quataert 1998) where the heat, which is liberated locally by a disc annulus due to viscous dissipation, is advected with the radial bulk flow and eventually disappears into the central black hole, without being emitted towards the observer.

The presence of reflection components and the Fe Kα line observed in X-ray spectra of black hole X-ray binaries (BHXBs) (Ebisawa et al. 1996) strongly suggest that cold gas from a disc and hot gas must coexist in a small volume. A possible configuration is the radial transition from an outer SSD to an inner ADAF. A bimodal disc model along these lines, with the transition radius as a free parameter and with the ADAF and SSD parts treated as separate entities, allows one to fit the spectra of persistent X-ray sources (Narayan, McClintock & Yi 1996; Manmoto, Mineshige & Kusunose 1997; Esin et al. 2001). Furthermore, this model also explains the different spectral states observed in some soft X-ray transients by varying the location of the transition radius (Esin et al. 1998). A number of these sources also show high-frequency quasi-periodic oscillations (QPO) in their Fourier-power spectra (Nowak & Lehr 1998). These examples demonstrate the necessity for a time-dependent treatment.

In a consistent theory, the disc transition radius is obtained naturally as a function of flow parameters and boundary conditions. For the stationary limit, this issue was addressed by Manmoto & Kato (2000). If the transition region between the inner ADAF and the outer SSD is sufficiently narrow (i.e. δR/R≪ 1), specific connection conditions can be derived which in turn significantly constrain the parameter space available for global solutions. The most restrictive condition is that the ADAF must support a minimum outward energy flux into the SSD, which is then liberated as radiation. The nature of this energy flux is still a matter of debate, in particular within the one-dimensional approximation (Abramowicz, Björnsson & Igumenshchev 2000). Honma (1996) and Manmoto & Kato (2000) specifically considered energy transport by turbulent convection, and conclude that the transition radius for stationary flows is then constrained to a radial domain with boundaries depending on the model parameters.

2 The Disc Model

We investigate time-dependent transition models by solving the full Navier-Stokes equations for a vertically integrated axisymmetric disc. In the energy equation for the gas we include dissipation of a turbulent convective energy flux, and in addition we solve a radiative transfer equation in the diffusion approximation. Our disc model thus generalizes the stationary transition model of Manmoto & Kato (2000) to the non-stationary case.

Assuming vertical hydrostatic equilibrium, the vertical pressure scale height is H=c_s/Ω_K, where c_s is the isothermal sound speed and Ω_K is the Keplerian angular velocity associated with the pseudo-Newtonian gravitational potential for the space-time of a non-rotating black hole with mass M (Paczynsky & Wiita 1980). As a result of the vertical averaging procedure (e.g. Narayan & Yi 1995a), the thermodynamic variables which appear in the model equations are all vertically integrated quantities, such as surface density Σ= 2Hρ, total pressure P, internal energy density E and radiation energy density U.

The continuity equation and the momentum equations in radial and azimuthal direction are

(1)

(2)

(3)

where R is the cylindrical radius and l is the specific angular momentum. Viscous transport is modelled as a Navier-Stokes stress, i.e. the vertically integrated stress tensor τ= 2μσ, with σ the shear tensor and μ the vertically integrated dynamical viscosity. We adopt the α-viscosity law of Shakura & Sunyaev (1973) as a model for anomalous turbulent viscosity, ν=αc_sH, which results in μ =Σν=Σαc_sH.

The thermal disc structure is governed by the magnitude of the various heating and cooling terms which contribute to the energy equation for both the gas and radiation field. The gas energy equation is written as

(4)

The vertically integrated viscous heating rate is given by Q⁺_vsc= (τ·∇) · v, with τ as given above.

As the entropy gradient in ADAFs decreases outward, ADAFs are convectively unstable (Honma 1996). Thus we additionally consider a bulk energy transport due to a convective energy flux F_trb with the associated net heating rate given by Q⁺_trb=−∇·F_trb. A natural assumption is to set F_trb proportional to the local entropy gradient, i.e.

(5)

where κ_T is an effective diffusion coefficient. Following Honma (1996), Manmoto et al. (2000) and Manmoto & Kato (2000) we parametrize κ_T in analogy with the α-viscosity prescription as κ_T=α_Tc_sH. In general, α_T≤α is expected, because convective turbulence contributes to the viscosity, but not all sources of viscosity produce a bulk energy transport (Ryu & Goodman 1992). The rate at which heat is advected with the flow is Q⁻_adv=ΣT(v · ∇)s=∇·Ev+P_g∇· v. The radiative cooling rate is Q⁻_rad= 2cρHκ_P (a_RT⁴−U/2H), where κ_P is the Planck mean opacity for free-free absorption.

Finally, the radiation energy equation is

(6)

where Q⁺_src=Q⁻_rad. Sinks for the radiative energy density originate from radiative losses through the upper and lower disc surfaces Q⁻_esc and from radial diffusion Q⁻_diff=∇·F_rad. Following Manmoto & Kato (2000), we write

(7)

where κ_R is the Rosseland mean opacity including electron scattering and free-free absorption. The equations are closed by an equation of state for the pressure, P=P_g+P_rad. We apply P_g= (γ− 1)E, with γ the ratio of specific heats, together with P_rad=U/3.

We discretize the equations on a non-uniform radial mesh, and use a time semi-implicit, operator-split scheme to advance the solution forward in time. The most restrictive constraint on the time-step in explicit schemes is due to radiative and viscous terms as a result of the associated time-scales. We avoid this problem by treating these terms implicitly. Implicit methods for finite volume discretization of our equations result in a large sparse matrix, where only the neighbourhood of the diagonal is populated. The coupling between individual equations gives rise to a few off-diagonal entries. If this coupling is weak, e.g. if the time increment is small, these off-diagonal terms may be neglected at low order and the resulting matrix can be inverted very efficiently. In fact, we first solve each equation independently, then update the dynamical variables and finally recompute the global defects with the full system. The procedure is now repeated until global convergence is reached. This scheme was originally proposed by Hujeirat & Rannacher (1998) and is naturally advantageous for moderately coupled systems. If the coupling is too strong, the efficiency in terms of achievable physical simulation time over computational CPU-time is very low, and may even drop below that of an explicit scheme.

In this paper, we restrict attention to the model parameters , where is the mass accretion rate of the initial model in Eddington units. The black hole mass and the adiabatic index are held fixed at M= 10 M_⊙ and γ= 1.5 in all models. As initial data we take a self-similar ADAF profile of Narayan & Yi (1994) where the advection efficiency parameter f is a step function, i.e. f drops from unity (initial ADAF) to zero (initial SSD) at an arbitrarily chosen initial transition radius. At the outer boundary, located at 2000R_g– where R_g=GM/c² is the gravitational radius — we fix the accretion rate to and the specific angular momentum to the local Keplerian value. Vanishing gradient boundary conditions are applied to the remaining variables. At the inner boundary at 2.2R_g we apply free outflow boundary conditions.

3 A Fiducial Solution

The transition radius R_tr is defined as the location where the vertical optical depth is equal to unity. Variables at the transition radius are denoted by the subscript ‘tr’, and the subscript ‘K’ refers to Keplerian motion. Time is measured in units of the Keplerian orbital period taken at the initial transition radius, P^tr_K, unless explicitly noted otherwise. Our unit of length is the gravitational radius R_g. In Fig. 1 we plot the contributions to the energy balance for a fiducial model with parameters . The transition in the initial conditions was located at R_tr= 55R_g. These parameters were chosen to facilitate comparison with other published results.

Figure 1.

Radial profile of the contributions to the gas energy balance (equation 4) for our fiducial solution. Note that Q⁺_trb does act as a cooling or heating rate at different radii, i.e. it changes sign in the transition region. The dots superimposed on the Q⁺_trb curve indicate the spatial resolution.

We generally identify three domains within the flow according to the magnitude of the different terms in the stationary energy equation (4), even though our solutions are intrinsically time-dependent and cannot be directly compared with steady-state profiles. With increasing radial distance from the central black hole we identify the inner ADAF where the optical depth is very low. Next we define the transition region as the zone where the rotation is super-Keplerian. Super-Keplerian rotation is present whenever the transition region is narrow (cf. Abramowicz, Igumenshchev & Lasota 1998). Between the outer edge of the transition region and the outer edge of the disc we observe the SSD where the optical depth is very large. The radial profiles are never stationary, as discussed below. Nevertheless, the time-averaged profiles are rather similar to the stationary solutions of Manmoto & Kato (2000).

4 Evolution of the Transition Radius

The transition radius R_tr derived from our solutions is a function of time. We find that R_tr is always located within the transition region. After a short initial relaxation phase, we generally observe a quasi-periodic time variation of R_tr, superimposed on a much slower outward directed drift. This behaviour is shown in Fig. 2(a).

Figure 2.

(a) Temporal evolution of the transition radius R_tr (dotted line). The plus signs + and crosses × indicate the innermost and outermost location of the transition radius. The latter is fitted by an analytical function (solid line). (b) Dependence of the oscillation cycle period Δt_osc on the transition radius for our fiducial model. In both panels time is measured in units of the Keplerian orbital period P^tr_K at the initial transition radius located at R_tr= 55R_g.

The mass flow rate through the transition region is not constant in time, but modulated with the same quasi-periodic pattern as R_tr. The transition region tends to deposit cold gas into the ADAF abruptly. Whenever this happens, the transition region grows at the expense of the ADAF.

On top of these rapid oscillations we observe a slow outward drift of the transition radius and consequently, the SSD shrinks at this very rate while the ADAF grows. For the fiducial model introduced above, R_tr grows from initially 55R_g to finally 123R_g at t= 92P^tr_K. The associated drift velocity decreases with increasing R_tr or, accordingly, with increasing time. A drifting behaviour of this type is perceived in all our simulations. However, the overall bimodal disc structure is never disrupted and the transition region is always narrow in the sense defined above.

At t= 92P^tr_K the time-step of our simulation dropped drastically to very small values, and the simulation was halted. This behaviour was seen in all our simulations (Gracia 2002). Interestingly, the values we measure for the transition radius at that time is nearly equal to the maximum transition radius predicted for a stationary bimodal disc model with the same parameters and boundary conditions (Manmoto & Kato 2000). Near this maximum transition radius the equations become very stiff, and thus the convergence of our numerical method decreases while the time-step freezes. The issue of the drift as such is further discussed in Section 4.2.

4.1 The oscillation cycle

The period of an oscillation cycle Δt_osc increases with time. In Fig. 2(a) we plot the innermost and outermost locations which the transition radius reaches during the complete evolution of our fiducial model. The outer location is treated as the generic transition radius and fitted by an analytical function . The innermost location was in turn considered as the perturbed transition radius in the course of the oscillation cycle and used to derive the oscillation period Δt_osc.

As shown in Fig. 2(b), a strong correlation is found between Δt_osc and the generic transition radius . Furthermore, Δt_osc is very close to (within a factor of the order of unity) the Keplerian orbital period at the generic transition radius, , independent of model parameters . However, the width of the transition region is a function of the adiabatic index γ (Manmoto & Kato 2000) and therefore we expect it to modify the proportionality factor between Δt_osc and .

A necessary condition for a narrow transition region is super-Keplerian rotation within the transition region, i.e. l > l_K (Abramowicz et al. 1998). In our calculations the specific angular momentum profile decreases outward and smoothly matches the Keplerian rotation profile of the SSD. However, if l decreases with increasing radius, ∂_Rl < 0, the flow is Rayleigh unstable. Consequently, a narrow ADAF-SSD transition region is the location of putative instability.

Outward decreasing specific angular momentum profiles have previously been obtained by e.g. Honma (1996) and Manmoto & Kato (2000), as seen from their figures. The specific angular momentum profile of Abramowicz et al. (1998) does not decrease outward. They used an adiabatic equation of state, which yields broader transition regions, i.e. l changes slower near the outer edge of the transition region and matches the SSD profile without actually decreasing.

Local stability analysis of a rotating flow predicts linear perturbations to oscillate with frequency ω that is bounded by the epicyclic frequency κ. Perturbations with ω²≈κ²_max will preferably grow or decay (Kato & Manmoto 2000; Kato 2001, for a review). If κ² is negative the flow is dynamically unstable. If κ² is positive the perturbation will oscillate with frequency ω² near κ²_max. Kato & Manmoto (2000) further argue that the oscillations will be trapped within the transition region, as they are reflected at the transition region boundaries.

In our solutions we find an upper limit of κ²_max/Ω²_K≈ 2, i.e. the oscillation period Δt_osc/P^tr_K has a lower limit of the order of unity, consistent with our numerical results and with our proposition that the oscillations have a physical nature. On the other hand, we note that κ²_max sensitively depends on the angular momentum gradient in the transition region, which is rather poorly resolved in our current models. Therefore, we cannot conclude from our simulations that Δt_osc is always bounded by the above limit.

We mentioned earlier that perturbations generated within the transition region are expected to be effectively trapped within the same. In reality, this simple picture does not survive the confrontation with numerical experiments entirely unharmed. Perturbations from within the transition region do back-react with the flow outside. Fig. 3 illustrates the dynamical evolution of the flow during an oscillation cycle at high temporal resolution.

Figure 3.

Dynamics of the oscillations in the ADAF-SSD transition region. On each panel radius R= (2R_g⋯ 2000 R_g) increases from left to right on a logarithmic scale, while time t=[0… 0.92 P^tr_K (62R_g)] runs from bottom to top on a linear scale. The individual panels show logarithm of temperature logT (left-hand panel), and logarithm of vertically integrated radial momentum density log(|Σv_R|) (right-hand panel). T is measured in units of (γ− 1) c²μm_p/k_B and Σv_R in units of . The outer black contours indicate the zone of super-Keplerian rotation, i.e. the transition region. The inner black contour indicates the transition radius R_tr. This figure is available in colour in the on-line version of the journal on Synergy.

The oscillation cycle appears to be initiated by a perturbation as seen, e.g. in the radial momentum flux density, approaching the inner boundary of the transition region. The perturbation leaks into the ADAF, and cold matter with high angular momentum detaches from the transition region and penetrates the ADAF as seen in Fig. 3. The ejected blob carries nearly constant angular momentum, such that the region of super-Keplerian rotation, i.e. the transition region, grows far into the ADAF. When the perturbation reaches the sonic radius, it is partly reflected as an outgoing wave and sweeps up hot plasma towards the transition. In continuation the ADAF relaxes slowly back into the unperturbed state. This causes high variability not only of the transition region, but also of the inner ADAF, where most of the X-ray flux is expected to originate. Hence the evolution of our models is strongly affected by the existence of oscillations in the transition region.

4.2 Radial drift or the strong ADAF principle

The outward drift of the transition region which we see in our simulations also requires interpretation. It is well known that, for a given accretion rate, the radial extent of the ADAF (if existent) is restricted by an upper limit R_τ (e.g. Chen et al. 1995; Abramowicz et al. 1995). At every radius R > R_τ, the optical depth exceeds unity and radiative cooling is efficient, such that the ADAF solution branch is no longer available. For certain R < R_τ, the ADAF and SSD-like solutions are both available in principle. From entropy considerations we expect that nature realizes the ADAF whenever possible, a proposition referred to as the strong ADAF principle (Narayan & Yi 1995b).

Our models do indeed seem to support this hypothesis. The transition radius does always drift outward, thereby increasing the total entropy of the entire solution. The drift is stopped once R_tr has reached the outermost location, R_trb, that is allowed from stationary analysis. Note that R_trb is smaller than R_τ because turbulent convective energy transport is taken into account in our models.

For a given set of parameters, the specific entropy of the ADAF is, at any radius of interest, larger than the specific entropy of the SSD. The particular form applied for turbulent convective energy transport (namely along the negative entropy gradient), implies that the transition region is always supplied with energy from the inner ADAF. As long as this energy surplus cannot be completely shed by radiative losses, the transition region is heated to nearly the local virial temperature and the flow switches to the ADAF mode. As a consequence, the SSD is slowly evaporated radially at its inner edge. This process does not stop until radiative energy shedding in the transition region completely dominates energy supply by convection — exactly the definition of the outermost allowed transition radius R_trb. Our models do indeed evolve in general within the parameter space of Manmoto & Kato (2000, see their figs 8, 9 and 10).

5 Concluding Remarks

We have modelled the time-evolution of bimodal accretion discs where an outer cooling-dominated disc matches an inner ADAF-type flow through a narrow transition region. Our models show an intrinsically time-dependent behaviour that can be understood from the dynamical instabilities predicted for stationary ADAF-SSD transitions. Specifically, oscillations with frequencies close to Ω_K are excited in the transition region. We associate these oscillations with a mutual local instability that originates from the radially decreasing specific angular momentum profile, a characteristic feature of any narrow-transition bimodal disc.

We observe that the waves excited in the transition region leak into the inner ADAF, but do not essentially propagate outwards into the SSD. These oscillations give rise to rapid periodic modulations of the ADAF and the transition region. We estimated the total disc luminosity by integrating the radiative cooling rate over the disc surface and again find a periodic modulation with a variability of a few per cent. This could possibly explain the QPOs observed in X-ray transient systems (e.g. Nowak & Lehr 1998, and references therein). Kato & Manmoto (2000) estimated the oscillation frequency in bimodal disc transitions to a small fraction ∼10⁻³ of Ω_K, albeit under the assumption of small frequency and using a rather crude estimate for the disc structure. Our numerical results show frequencies very close to Ω_K at the transition radius.

We furthermore observe that the initial transition radius does slowly drift to larger radii, i.e. even though a stationary transition is allowed for a wide range of radii, the disc appears to favour the outermost possible transition radius. This fits well into the strong ADAF principle hypothesis (Narayan & Yi 1995b). Within our modelling, convective energy transport is the key to make this drift possible, i.e. to evaporate the SSD radially until the transition reaches its outermost allowed location. But even if we neglect convective energy transfer as modelled in this paper, the large entropy gradient across the transition still persists. Significant transfer of heat from the ADAF to the SSD is a necessary condition for the existence of a stationary transition front. Whether turbulence is the relevant physical mechanism for this heat flux and whether the shape of the transition front is actually close to cylindrical cannot be addressed by one-dimensional modelling. The resolution of this issue requires at least two-dimensional simulations in the parameter range that allows for a transition in the one-dimensional stationary case. This region of the parameter space, however, has not been addressed in a two-dimensional simulation by Hujeirat & Camenzind (2000). The gap between vertical disc evaporation models (Meyer & Meyer-Hofmeister 1994; Meyer, Liu & Meyer-Hofmeister 2000) and the present one-dimensional radial transition physics therefore remains to be closed in a future work.

Acknowledgments

Part of this work was supported by the Deutsche Forschungsgemeinschaft (DFG) (Sonderforschungsbereich 439 and Sonderforschungsbereich 359).