Evolution of bimodal accretion flows

We consider time-dependent models for rotating accretion flows onto 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 expected for the associated total luminosity. In particular, rapid oscillations of the transition radius are present on a timescale 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 QPO phenomenon.


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 scale. This so-called SSD model cannot explain, however, the hard spectral component observed in most compact accreting systems, e.g., in X-ray binaries or in 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 et al. 1998) where the heat, that 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 BHXB's (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 this line, with the transition radius as a free parameter and with the ADAF and SSD parts treated as separate entities, allows to fit the spectra of persistent X-ray sources (Narayan et al. 1996;Manmoto et al. 1997;Esin et al. 2001). Furthermore, this model does also ⋆ E-mail: J.Gracia@lsw.uni-heidelberg.de explain 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 do 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 . 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 et al. 2000). Honma (1996) and  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.

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 generalises the stationary transition model of  to the non-stationary case.
Assuming vertical hydrostatic equilibrium, the vertical pressure scale height is H = cs/ΩK, with cs the isothermal sound speed and ΩK 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 where R is the cylindrical radius and l is the specific angular momentum. Viscous transport is modeled 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, ν = αcsH, which results in µ = Σν = ΣαcsH. 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 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., with κT an effective diffusion coefficient. Following (Honma 1996; we parametrise κT in analogy with the α-viscosity prescription as κT = αTcsH. In general, αT α is expected, since 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 + Pg∇·v. The radiative cooling rate is Q − rad = 2cρHκP (aRT 4 − U/2H), with κP the Planck-mean opacity for free-free absorption.
Finally, the radiation energy equation is where Q + src = Q − rad . Sinks for the radiative energy density originate from radiative losses through the upper and lower disk surfaces Q − esc and from radial diffusion Q − diff = ∇·F rad . Following , we write 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 = Pg + P rad . We apply Pg = (γ − 1)E, with γ the ratio of specific heats, together with P rad = U/3. We discretise 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 timescales. We avoid this problem by treating these terms implicitly. Implicit methods for finite volume discretisation of our equations result in a large sparse matrix, where only the neighbourhood of the diagonal is populated. The coupling between individual equations gives raise 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 letter we restrict attention to the model parameters (ṁi, α, αT), whereṁi 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 = 10M⊙ and γ = 1.5 in all models. As initial data we take a self-similar ADAF profile of Narayan & Yi (1994) with the advection efficiency parameter f a step-function, i.e., with f dropping from unity (initial ADAF) to zero (initial SSD) at an arbitrarily chosen initial transition radius. At the outer boundary, located at 2000Rg -where Rg = GM/c 2 is the gravitational radius -we fix the accretion rate toṁi 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.2 Rg, we apply free outflow boundary conditions.

A FIDUCIAL SOLUTION
The transition radius Rtr 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 Rg. In Fig. 1 we plot Figure 1. Radial profile of the contributions to the gas energy balance (eq. 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.
the contributions to the energy balance for a fiducial model with parameters (ṁi, α, αT) = (0.021, 0.4, 0.4). The transition in the initial conditions was located at Rtr = 55Rg. These parameters were chosen to facilitate comparison with other published results.
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 et al. 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 .

EVOLUTION OF THE TRANSITION RADIUS
The transition radius Rtr derived from our solutions is a function of time. We find that Rtr is always located within the transition region. After a short initial relaxation phase, we generally observe a quasi-periodic time variation of Rtr, superimposed on a much slower outward directed drift. This behaviour is shown in Fig. 2(a). The mass flow rate through the transition region is not constant in time, but modulated with the same quasiperiodic pattern as Rtr. The transition region tends to abruptly deposit cold gas into the ADAF. 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, Rtr grows from initially 55 Rg to finally 123 Rg at t = 92 P tr K . The associated drift velocity decreases with increasing Rtr 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 = 92 P 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 . 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.

The oscillation cycle
The period of an oscillation cycle ∆tosc increases with time. In Fig. 2(a) we plot the innermost location and the outermost location which the transition radius takes during the complete evolution of our fiducial model. The outer location is treated as the generic transition radius and fitted by an analytical function Rtr(t). 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 ∆tosc.
As shown in Fig. 2(b), a strong correlation is found between ∆tosc and the generic transition radius Rtr(t). Furthermore, ∆tosc is very close to (within a factor of order unity) the Keplerian orbital period at the generic transition radius, PK(Rtr(t)), independent of model parameters (ṁi, α, αT). However, the width of the transition region is a function of the adiabatic index γ  and therefore we expect it to modify the proportionality factor between ∆tosc and PK(Rtr(t)).
A necessary condition for a narrow transition region is super-Keplerian rotation within the transition region, i.e. l > lK (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  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 κ. will preferably grow or decay Kato 2001, for a review). If κ 2 is negative the flow is dynamically unstable. If κ 2 is positive the perturbation will oscillate with frequency ω 2 near κ 2 max .  further argue that the oscillations will be trapped within the transition region, since they are reflected at the transition region boundaries.
In our solutions we find an upper limit κ 2 max /Ω 2 K ≈ 2, i.e., the oscillation period ∆tosc/P tr K has a lower limit 1/ √ 2 of order unity, consistent with our numerical results and with our proposition that the oscillations have a physical nature. On the other hand, we note that κ 2 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 ∆tosc 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.
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, 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 contin-uation 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 Xray flux is expected to originate. Hence the evolution of our models is strongly affected by the existence of oscillations in the transition region.

Radial drift or the strong ADAF principle
The outward drift of the transition region which we see in our simulations does also require 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 SSDlike 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 Rtr(t) does always drift outward, thereby increasing the total entropy of the entire solution. The drift is stopped once Rtr 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 radia- Figure 3. Dynamics of the oscillations in the ADAF-SSD transition region. On each panel radius R = (2Rg · · · 2000Rg) increases from left to right on a logarithmic scale, while time t = (0 · · · 0.92 P tr K (62 Rg)) runs from bottom to top on a linear scale. The individual panels show logarithm of temperature log T (left), and logarithm of vertically integrated radial momentum density log(|Σv R |) (right). T is measured in units of (γ − 1)c 2 µmp/k B and Σv R in units ofṀ Edd /(2πRg). The outer black contours indicate the zone of super-Keplerian rotation, i.e. the transition region. The inner one, the transition radius Rtr.
tive 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).

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 frequency 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 outward into the SSD. These oscillations give rise to rapid periodic modulations of the ADAF and the transition region. We estimated the total disk luminosity by integrating the radiative cooling rate over the disk surface and again find a periodic modulation with a variability of a few percent. This could possibly explain the QPOs observed in X-Ray transient systems (e.g. Nowak & Lehr 1998, and references therein).  estimated the oscillation frequency in bimodal disc transitions to a small fraction ∼ 10 −3 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 1-dimensional modeling. The resolution of this issue does require at least 2-dimensional simulations in the parameter range that allows for a transition in the 1-dimensional stationary case. This region of the parameter space, however, has not been addressed in a 2-dimensional simulation by Hujeirat & Camenzind (2000). The gap between vertical disc evaporation models (Meyer & Meyer-Hofmeister 1994;Meyer et al. 2000) and the present 1-dimensional radial transition physics therefore remains to be closed in a future work.