Magnetized Accretion Flows: Effects of Gas Pressure
We study how axisymmetric magnetohydrodynamical (MHD) accretion flows depend on adiabatic index in the polytropic equation of state. This work is an extension of Moscibrodzka & Proga, where we investigated the dependence of 2-D Bondi-like accretion flows in the hydrodynamical (HD) limit. Our main goal is to study if simulations for various can give us insights into to the problem of various modes of accretion observed in several types of accretion systems such as black hole binaries (BHB), active galactic nuclei (AGN), and gamma-ray bursts (GRBs). We find that for , the fast rotating flow forms a thick torus that is supported by rotation and gas pressure. As shown before for , such a torus produces a strong, persistent bipolar outflow that can significantly reduce the polar funnel accretion of a slowly rotating flow. For low , close to 1, the torus is thin and is supported by rotation. The thin torus produces an unsteady outflow which is too weak to propagate throughout the polar funnel inflow. Compared to their HD counterparts, the MHD simulations show that the magnetized torus can produce an outflow and does not exhibit regular oscillations. Generally, our simulations demonstrate how the torus thickness affects the outflow production. They also support the notion that the geometrical thickness of the torus correlates with the power of the torus outflow. Our results, applied to observations, suggest that the torus ability to radiatively cool and become thin can correspond to a suppression of a jet as observed in the BHB during a transition from a hard/low to soft/high spectral state and a transition from a quiescent to hard/low state in AGN.
keywords:accretion, accretion discs – black hole physics – X-ray: binaries – galaxies: nuclei – methods: numerical –MHD
Many astrophysical objects powered by accretion on to a compact objects like e.g. active galactic nuclei – AGN, black hole binaries – BHB, gamma-ray bursts – GRBs, show signatures of outflows, either in a form of wind or jet, or both. Large scale, collimated jets are being detected from radio to -ray wavelengths in AGN and microquasars (e.g. Königl 2006, Markoff 2006). Recently an increasing radio emission from the core of some BHB was associated with small scale outflow/jet-like structures (Fender et al., 2004). The observations of BHB suggest that the power of the jet-like outflow may be related to the evolutionary stage of the object.
The BHB exhibit several X-ray spectral states, which can change in time along a hysteresis loop on the luminosity-hardness diagrams. BHB (as well as galactic nuclei) spend most of their time in the quiescent/off mode during which their X-ray luminosity is low (Gallo et al., 2003). The discovery of 17 out of 20 known BHB was possible only during an outburst (X-ray novae) accompanying the transition from the quiescent/off to a low/hard spectral state after which the X-ray brightness increases three or more orders in magnitude (Remillard & McClintock, 2006). In the low/hard spectral state the non-thermal power-law component dominates in the X-ray energy band. After the transition into the so-called high/soft state the spectrum is dominated by the thermal multi-blackbody emission coming from a cold Keplerian disc (Shakura & Syunyaev 1973). The source returns to the quiescent state through the soft-state, which closes the hysteresis cycle. It is also well established that spectral states of BHB can change on relatively short time-scales from days and months to years (Remillard & McClintock, 2006) and that the X-ray lumonisty of the BHB does not necessarily correlate with its spectral state (Done et al., 2007).
The spectral state changes probably reflect the changes in the accretion rate/pattern and the formation/suppression of an outflow. For example, in the low/hard state, the BHB radio luminosity is strongly correlated with the X-ray luminosity (, Gallo et al. 2003). This correlation does not hold after the transition to the high/soft state, when the radio emission becomes suppressed (Corbel et al. 2001, Gallo et al. 2003). The X-ray variability and radio activity strongly support a model, where in the low/hard state, a standard thin accretion disc is truncated and does not reach the black hole horizon. During low/hard state, the inner disc that used to be during the high/soft state is replaced by a hot thick flow, which probably acts as a lunching site for a steady jet with the Lorentz factor, , typically less than 2. (e.g. Heinz & Sunyaev 2003; Meier 2005; Remillard & McClintock 2006; Fender et al. 2004 and references therein).
The transition from low/hard to high/soft spectral state occurs via an intermediate state, the so-called very high state. The very high state is characterized by strong black body and substantial power-law emission. In the very high state, the radio emission exhibits flaring. The flaring can be due to a fast non-stationary optically thin outflow with high Lorentz factor, 2 (Fender et al., 2004).
A radio-X-ray correlation, similar to that seen in BHB, is also expected to work in much larger scale systems like AGN. Because of much longer evolutionary time-scales in AGN, in comparison to BHB, the problem may be studied statistically only. The expectation of the correlation is supported by the fact that the radio-loud quasars have higher X-ray fluxes than radio-quiet quasars (Zdziarski et al., 1995). In this scenario, the AGN counterparts of the ‘quiescent’ state are Low Luminosity AGN (LL AGN), because they lack the ’big blue bump’ in their spectra, are weak in X-rays and are loud in radio in comparison to AGN (e.g. Ho 2003; Kawakatu et al. 2009). Our own galactic center - Sgr A* - is considered as an extreme example of LL AGN. The models of the quiescent/off phase of AGN and BHB, are usually based on the two assumptions: either most of the accretion energy is transferred into the outflow kinetic energy which results in a reduction of the mass-accretion rate, , and lower luminosity, or the accretion is highly inefficient, or both (Ichimaru 1977, Narayan et al. 1998,Yuan et al. 2003).
From the theoretical point of view, it is still unclear what induces the accretion phase transition in BHB, AGN, and GRBs. Many authors considered 1-D, vertically integrated, models where the cold disc evaporates, transitions into a hot phase and then re-condensates (e.g. Liu et al. 1999, Meyer et al. 2000, Różańska & Czerny 2000, Liu et al. 2002, Meyer et al. 2007). In 1-D models, the mass-accretion rate is the key parameter responsible for the disc behaviour. These models include such physical processes as mass transfer between hot and cold phase, thermal conduction, and radiative transport. However, the 1-D models do not include an outflow formation and do not yet attempt to reproduce observed correlation between X-rays and radio luminosity probably caused by outflow, during low/hard state. As underlined in Abramowicz et al. (2000), at least a 2-D treatment of the problem is needed to correctly model the accretion state transition process.
In this paper, we present results of 2-D axisymmetric magnetohydrodynamical (MHD) simulations of a slowly rotating accretion flow on to a black hole (BH). We investigate the time evolution and properties of inflows/outflows for various values of an adiabatic index, . Our work is an extension of work of Mościbrodzka & Proga (2008) (hereafter MP08) who considered models for various but in the hydrodynamical (HD) limit (see also Proga & Begelman 2003a, hereafter PB03a, for HD simulations with , work of Proga & Begelman 2003b, hereafter PB03b, who considered MHD models only for , and finally work of Proga et al. (2003) who considered a MHD collapsar model for GRBs where the inner accretion flow is radiation pressure dominated and is effectively ).
In MP08, we showed that in the HD limit, many properties of the accretion flow strongly depend on . Most relevant to this work is that the thickness of the torus decreases with decreasing whereas in units of the Bondi rate, increases with decreasing . Here, we want to check how affects magnetized flows. For example, is the efficiency of producing outflows sensitive to the thickness of the torus? This question is related to the observations of BHBs and AGN. One could argue that magnetically driven outflows should be the strongest from geometrically thick discs/torii with significant pressure support and the weakest, if any, from geometrical thin discs with little pressure support (e.g. Livio et al. 1999; Meier 2001).
To our best knowledge, there are no direct simulations that would confirm the above expectation mainly because previous simulations considered thick, pressure/rotation supported torii. What is clear is that most of these previous global MHD simulations of accretion flows show outflows (e.g. Uchida & Shibata 1985; Stone & Norman 1994; Stone & Pringle 2001; Hawley & Balbus 2002; PB03b; Proga et al. 2003;De Villiers et al. 2005; Mizuno et al. 2004; Kato et al. 2004; McKinney & Gammie 2004; Proga 2005, Fragile & Meier 2009). However, it is unclear under what conditions (e.g., for what torus thickness) the torus stops to produce an outflow.
Our very simplistic approach to control the torus thickness is dictated by a complexity of the relevant physics. A full treatment of an accretion disc, especially geometrically thin one, requires modeling energy dissipation, radiative heating and cooling, and radiative transfer. In our approach, for = 5/3 and 4/3 the torus will be hot and supported against gravity by gas pressure and rotation. As such it will be geometrically thick. On the other hand, for the torus will be supported only by rotation and therefore it will be Keplerian and geometrically thin (or at least considerably thinner than for ). Therefore, our study of effects of will allow us to gain some insights into the problem of the outflow formation for various physical conditions and different phases of activity like e.g. soft and hard spectral states in BHB and AGN activity modes.
Here we will not be able to immediately compare our results with observations in any great detail (as it was possible in Mościbrodzka et al. 2007 where we computed and compared radiative properties predicted by PB03b’s MHD simulations with observations of Sgr A*). However, our approach enable us to produce results that can be easily and consistently compared with a great body of data from complementary numerical simulations that were preformed within one framework and with the same numerical tools (i.e.,PB03a; PB03b; Proga et al. 2003; Proga 2005, MP08).
To calculate the structure and evolution of an accreting flow, we solve the equations of MHD
where , P, , , and is the mass density, gas pressure, fluid velocity vector, internal energy density, gravitational potential, and magnetic field vector, respectively. Lagrangian derivative is defined as: . We relate and through the equation of state , where is an adiabatic index.
We perform simulations using the pseudo-Newtonian potential introduced by Paczyński & Wiita (1980) (here after PW)
where is the Schwarzschild radius. This potential approximates general relativistic effects in the inner regions, for a non-rotating BH. In particular, the PW potential reproduces the last stable circular orbit at as well as the marginally bound orbit at .
2.2 Initial conditions and boundary conditions
Generally, the setup and numerical methods used in our simulations are as in PB03b. Here, we only briefly summarize the initial and boundary conditions and spell out the differences between our and PB03b simulations.
Our calculations are performed in spherical polar coordinates . We assume axial symmetry about the rotational axis of the accretion flow ( and ). For the initial conditions of the fluid variables we follow PB03b and adopt a Bondi accretion flow. In particular, we adopt and computed using the Bernoulli function and mass accretion rate for spherically symmetric Bondi accretion with the PW potential. We set and mass of the central black hole which make our results directly applicable to the galactic center (e.g. it will be useful in our future calculations of synthetic spectra as in Mościbrodzka et al. 2007). However, our results could be rescaled and applied to other systems. We specify the sound speed at infinity, through a normalized Schwarzschild radius (, where is the Bondi radius. The normalized BH radius characterizes the gas temperature. We consider models with .
We specify the initial conditions by adopting a non-zero specific angular momentum l at the outer boundary. We consider a case in which the angular momentum at the outer radius depends on the polar angle via
This formula was used by PB03a and MP08. We note that PB03b use a step function so that this is one of the differences between the setup of our and PB03b’s simulations. However, as discussed by PB03b and shown below the details of the l distribution are not very important. We express the angular momentum on the equator as
where is the ‘circularization radius’ on the equator in units of for the Newtonian potential (i.e., at ).
To generate the magnetic field, we use vector potential, , i.e., . The initial magnetic field configuration is a pure vertical magnetic field. The zero divergence vertical field is defined by the potential =(, , ). We scale the magnitude magnetic field by plasma parameter at the inner boundary, .
where is the initial pressure at the inner radius.
Our standard computational domain is defined to occupy the radial range and the angular range . The domain is discretized into zones with 140 zones in the direction and 100 zones in the direction. We fix zone size ratios, , and for . With such defined grid we are able to resolve the fastest growing mode of the magnetorotaitonal instability (MRI) only in a case of thick torus (=5/3). See § 3.6.
The boundary conditions are specified as follows. At the poles, (i.e., and ), we apply an axis-of-symmetry boundary condition. At both the inner and outer radial boundaries, we apply an outflow boundary condition for all dynamical variables. As in PB03b, to represent steady conditions at the outer radial boundary, during the evolution of each model we continue to apply the constraints that in the last zone in the radial direction, , (where ), and the density is set to at all times. Note that we allow to float. For the magnetic field, at the outer boundary, we apply an outflow condition. At the inner boundary, we follow Stone & Pringle (2001) and use negative stress condition (i.e. we enforce at ). To reduce the problem of high Alfvénic velocities in regions of low density, we set a lower limit to the density on the grid as and enforce it at all times in the simulations.
3 Results of the simulations
The parameters used for the simulations are summarized in Table 1. The Table columns (1)-(8) show respectively the number of a given run; the BH radius and the Bondi radius ratio (); the circularization radius compared to the Bondi radius (); the adiabatic index (); the initial magnetic plasma parameter () at ; the final time () at which we stopped each simulation in units of the Keplerian orbital time at ( for ); the same final time but in units of the Keplerian orbital time at , and finally the time-averaged mass accretion rate on to the BH in units of the Bondi accretion rate () over 1000 at the end of each run. The shorter duration of our simulations (compared to the duration of PB03b’s simulations) is caused by the high density gradients which appear in our runs for low , and which lead to significantly smaller time steps.
3.1 Velocity field and sonic surfaces evolution
In our simulations, the evolution of the poloidal velocity field initially proceeds in a similar manner as in the corresponding HD simulations (see e.g. panels 3 and 4 in fig.2 in PM08 for each ). Namely, after a short episode of the radial inflow, the matter with relatively high-l forms a torus in the equatorial plane around the BH (high-l is defined here as and low-l as , where is a critical specific angular momentum). For =5/3, 4/3, and 1.2, initially when the torus forms, a shock wave is produced. The shock wave forms because some of the gas forming a torus turns around and moves outward along the equator due to the centrifugal force and gas pressure or tries to move closer to the poles and accrete. For 1.01, the gas is the most gravitationally bound compared to other ’s. As such it accumulates near BH and does not flow out along the equator contrary to the higher cases.
In its early phase of the evolution, the torus is not able to accrete on to a BH, independently of , because the MHD turbulence has not yet fully developed to transoport the angular momentum. The matter in the polar regions does not feel any centrifugal barier, thus it flows radially toward the BH. In our MHD simulations and the HD simulations studied in PM08, the sonic surfaces topology evolves in a similar manner. It changes from a circular shape into two lobes where the lobes are elongated along the poles. As we reported in PM08, the evolution of the sonic surface is the slowest for =1.01, due to a lack on an equatorial outflow.
After a short time the magnetic field effects start to dominate the accretion flow evolution. The torus starts to accrete toward the BH. The biggest and expected difference between HD and MHD simulations is the development of bipolar outflows in MHD simulations regardless of . PB03b’s simulations well illustrate the general evolution of the magnetic fields. The polar outflows form at small radii () close to the poles and expand to the outer boundary of the computational domain along the symmetry axis. We find that the time at which the outflows form depends on : the lower , the later the outflow forms. Although initially HD and MHD flows evolve similarly, at the end of the simulations, the geometry and character of the accretion flows in the HD and MHD simulations, are significantly different. At the end of the simulation, the MHD flows are composed of a bipolar outflow propagating in the polar regions, and the inflow near the equator. The equatorial region is turbulent and strongly asymmetric for all ’s.
3.2 Mass accretion rate evolution
Fig. 1 shows the as a function of time, for each studied as listed in Table 1. The time is given in units of the orbital time at the inner boundary of the computational domain. A pair arrows in each panel mark the moments at which the polar outflow forms, in the Northern (N) and Southern (S) directions. The moment of the formation of the large scale polar outflow is defined as a moment when the matter stops accreting through the polar regions and the radial velocity sign changes sign from negative to positive. We determine this by inspecting the velocity maps and checking that the outflow velocity becomes larger than escape velocity at the outer boundary.
The upper-left panel in Fig. 1 shows for run 1. Initially, decreases from 1 to 0.3 and then remains constant until 2500. The value for of 0.3 is consistent with that obtained in the HD counterpart (see MP08 and PB03b), and it corresponds to the matter accreting in the polar funels in the initial HD phase of the simulation. However, later around t= 3700, suddenly decreases to a minimum value (below 0.001), and oscillates between the minimum value and 0.03. We note that the moment of the S outflow formation does not correspond exactly to the moment of the reduction to its minimum value. In general (the details will be discussed in the next sections) the mass accretion rate reduction is caused by the mentioned bipolar outflow, which pushes the radially infalling material away from the polar funels. Such a vacuuming of the polar regions makes the torus the only source of gas available for accretion. During later phases of the evolution (i.e., ), shows narrow spikes and dips similar to those found by PB03b, corresponding to the episodes of enhanced and suppressed torus accretion, respectively.
The upper-right panel in Fig. 1 shows the evolution for =4/3 (run 2). At the beginning but after the torus formation, 0.2, which is consistent with mass accretion rate found previously in the HD simulations. At this stage, the gas is being accreted from the polar regions. Later when an outflow forms at , suddenly decreases. For , the evolution is characterized by narrow spikes, dips (similarly to run 1), and additionally broader features with higher amplitudes. The mass accretion rate in the spikes can reach a level of about 0.02 and corresponds to the torus accretion. During the dips, the has a minimum value below 0.001, when the torus stops accreting. The rates seen in the two broad features (around t=5500 and t=7000) reach the level of 0.1. At even later times (), is almost constant around 0.01. The broad features have been found before by PB03b, and correspond to the non-rotating gas accretion that reenter the BH vicinity in the equatorial plane. Near the end of the simulation, the minima (dips) reappear.
For run 3 (the lower-left panel in Fig. 1), the evolution of the mass accretion rate differs significantly from runs 1 and 2. In particular, there are no sudden drops in . Instead, evolves quite smoothly, despite the formation of bipolar outflows. The curve is variable, and it shows two significant minima around t=7000 and at the end of the simulation. The mean level of is around 0.1, which is rather high compared to those of runs 1 and 2.
The lower-right panel in Fig. 1, presents for the model with =1.01 (run 4). Here, shows a strong time variability although the time averaged rate remains high compared to other models. There is one significant minimum in at the final stages of the simulation. The panel shows many periods of increased accretion, during which is higher than 1!
Interestingly, for =1.2 and 1.01, in spite of the presence of the bipolar outflows, does not decrease as sharply as for =5/3 and 4/3. It also does not show the characteristic features such as spikes and dips which were also found by PB03b. The different character of the mass accretion rate evolution in runs: 1, 2 and 3, 4, reflects two distinct outflow formation mechanisms.
3.3 Inflow/outflow rates
To investigate the mechanisms responsible for the polar outflow formation, we examine inflow/outflow/total fluxes of various physical quantities as functions of the distance from BH. The fluxes are averaged over the whole range and over 1000 before the end of each run.
We compute the fluxes of: mass, internal energy, kinetic energy, magnetic energy and total energy. The mass flux at a given radius, is given by:
where , is the density, and is the radial velocity. To calculate the mass inflow rate at a given radius, we include only contributions with , whereas to calculate the mass outflow rate , we include only contributions with . To compute the net (total) mass flow flux , we include all contributions regardless of the sign ( ). To compute the internal, kinetic, magnetic and total energy fluxes, we use the same approach as for the mass fluxes with small modification. Namely, in Eq. 9 is substituted by e, , and for the internal, kinetic, magnetic and total energy, respectively. We normalize all energy fluxes with .
In Fig. 2 and Fig. 3, we only show the total mass flux and total energy fluxes, respectively. We present it as a ) function so that the positive values will indicate the domination of an outflow and the negative - an inflow. The solid, dotted, dashed and dot-dashed lines correspond to runs 1, 2, 3, and 4, respectively (as marked in the panels). Below, we will also the main conclusion which will be based our analysis of all studied fluxes not just those presented in the figures.
The total mass flux (Fig. 2) as a function of radius is not constant. This indicates that overall the flows did not reach a steady state. However, in run 1 and 2, the total mass flux is constant in the inner parts of the computational domain only up to few . For models 3 and 4, the flux is nearly constant over a wider radial range, i.e., up to .
Our analysis of the various fluxes showed that in run 1, the total internal energy flux changes sign from negative (domination of the inflow flux) to positive (domination of the outflow flux) very close to the BH (). The similar behaviour characterizes run 2 and 3. Run 4 differs from other cases: its total internal energy flux is negative for all radii.
There are a few general trends in the behaviour of all fluxes in all runs. The strongest contribution to the total energy inflow flux (upper panel in Fig. 3) in the outer region of the computational domain (which we define here as the region for ), is from the internal energy, regardless of the value of , whereas the total energy outflow flux (middle panel in Fig. 3) is dominated by kinetic energy (with one exception for for which the internal energy outflow flux is even greater than kinetic one but only up to ).
In the inner part of the flows (for ), where the most of the outflow is produced, the total energy inflow flux is dominated by the kinetic energy for all runs.
The total energy outflow flux in the close vicinity to the BH (i.e. for ) is dominated by magnetic energy in runs 1, 2, and 3. In run 4 (), the total energy outflow flux is dominated by the kinetic and internal energy – the contributions from both components are comparable. This suggests that the centrifugal forces may be responsible for an outflow in run 4, rather than magnetic forces (runs 1-3). We are elaborating this point below.
3.4 Formation of a stationary outflow
In this section, we investigate how the large scale bipolar outflow is formed. We will do it in the context of previous work, namely PB03b, Proga et al. (2003) and Proga (2005). We will concentrate on identifying the flow components and checking if they are by nature the same ones as those found in previous simulations. Therefore we first briefly summarize the results from PB03b.
As mentioned in the introduction, PB03b considered models similar to ours but only for . They found that the nature of a weakly magnetized, slightly rotating accretion flow is controlled by magnetic fields when . As in the HD inviscid case, studied in PB03a (see also MP08), the magnetized material with forms an equatorial torus. However, in the MHD case, the torus accretes on to BH because of magnetocentrifugal instability (MRI; Balbus & Hawley 1991) and magnetic braking. The magnetic field in the torus amplifies and produces both a magnetic corona and an outflow. The latter two can be strong enough to prevent accretion of the low-l material through the polar regions, which was the source of accretion in the HD inviscid case.
Very similar time evolution and flow behaviour was found in simulations of the MHD collapsar model for gamma-ray bursts (Proga et al., 2003) which were similar to those of PB03b’s simulations, except for much more sophisticated micro-physics. In particular, Proga et al.’s considered an equation of state that included contributions from gas, photons and electrons. The most relevant aspect of Proga et al.’s simulations for us here is that the inner torus is radiation pressure dominated and the effective adiabatic index of the inner torus is . A general conclusion from PB03b and Proga et al. 2003 is that for the torus formation and evolution, including the development of a corona and outflow are qualitatively similar. In particular, the torus outflow is strong enough to suppress polar funnel accretion and to propagate out to a large distance from the BH. This result is perhaps unsurprising given that the torus is geometrical thick for .
The previous simulations showed that slightly rotating accretion flows are very variable and complex. In particular, the inner radius of the torus and the strength of the outflow change with time. For these reasons, only a later, much more detailed analysis of PB03b and Proga et al. (2003) simulations showed an additional flow component: a very collimated outflow driven from an infalling gas with (Proga, 2005). To articulate the basic physics that occurs in production of this collimated outflow, Proga (2005) performed also some new simulations. He considered gas with so that after an initial transient, the flow in the HD case accretes directly on to a BH without forming a torus. However, in the MHD case, even with a very weak initial magnetic field, the flow settles into a configuration with four components: (1) an equatorial inflow, (2) a bipolar outflow, (3) polar funnel outflow, and (4) polar funnel inflow. The second flow component of the MHD flow represents a simple yet robust example of a well-organized inflow/outflow solution to the problem of MHD jet formation.
In each of our present simulation, the accretion flow forms a torus and outflow although some low-l () inflow may occur for most of the time. We confirm that an outflow may form in two manners. First type of an outflow is driven by centrifugal and magnetic forces as it was already described in Proga (2005). The second type of outflow is a wind driven from a torus by the toroidal magnetic field gradient as described in PB03b (2003b and references therein). In the first case, the outflow is collimated and removes only a small fraction of the low-l matter from the polar funnel and does not affect the torus evolution, i.e. the torus is separated from the collimated outflow by an accreting stream of low-l matter. In the second case, the outflow is a broad wind from the torus. This wind can remove all low-l matter from the polar funnel (by pushing it sideways). The two types of outflows contribute the time variability of described in the previous section and of the flow structure detailed here.
In Figs. 4 and 5, we show the two-dimensional structures of the density and velocity fields along with the contours of the angular momentum of the gas. We also show the maps of the magnetic field parameter , including total magnetic filed (), at the times of the formation of an outflow on N and S side.
For =5/3 (run 1), the high-l matter (, dotted lines on the angular momentum contours) does not reach the BH and a torus does not form yet ( here, we define the torus as a high-l matter accreting in the equatorial region). Instead in the inner flow, the low- gas and the high-l gas are mixed. Such a state of accretion was not found in the previous simulations because here we assume a continuous distribution of l with at , whereas in the previous simulations of PB03b a step-function was assumed.
In run 1, the initial outflow is driven by the centrifugal and magnetic forces along the rotation axis as found by Proga (2005) (but not in PB03b). As we mentioned in the previous subsection does not decrease significantly during the formation of the initial outflow (Fig. 1 left-upper panel), but it does decrease slightly some time after the outflow is formed.
The polar outflow is highly collimated by the infalling ambient gas. The outflow that forms on the opposite side of the torus later during the simulation is very similar to that just described. At the moment of the outflow formation, the torus is not developed yet. The outflow is marked by the solid contours in the angular momentum panel in Fig. 4 (upper middle panel). As one can see this outflow is made of the slowly rotating matter near the symmetry axis.
In the later phases of the simulations (Fig. 5), after the bipolar outflow expansion, the high-l gas forms a well developed torus that extends inward to the last stable orbit and accretes on the BH. Subsequently, smaller scale (within 20 ) magnetic torus corona and outflow develop. In this phase, the torus outflow starts to push aside the low-l matter accreting from the polar regions. For the time up to t 6000, the high- and the low-l gases mix with each other, near the BH. Later, when the low-l matter is significantly pushed aside by the fully developed magnetic torus corona and outflow, the accretion flow becomes very similar to that in PB03b (see also Sect. 3.6).
For =4/3 (run 2, second row of panels in Figs. 4 and 5), the outflow is formed in the similar way as found by PB03b i.e. it is a collimated wind from a torus magnetic corona. In this case, the torus is more bound than in the case of =5/3. Therefore the high-l matter does not mix with the low-l matter near the equator close to the BH (the upper panels of fig. 8 in MP08 for the HD run 1). The magnetic torus corona and outflow form relatively fast. The corona in this case extends up to . Because the torus is geometrically thick, and the low-l matter is pushed aside very efficiently from the polar region, the polar inflow is completely suppressed ( see Fig. 5 middle panel corresponding to ). The moment of the sudden decrease of (Fig. 1) corresponds to the time of the outflow formation.
Next, we describe what happens for =1.01 ( run 4, bottom panels in Figs. 4 and 5). Already in the HD simulations (runs G, H, I, J in fig.8 in MP08), one may notice that because of much lower gas pressure, the torus is geometrically thinner compared to the cases with higher . Here, a strongly turbulent magnetic corona forms above the disc. The corona is geometrically thick, but it does not expand into the polar regions as in the previous cases. As in run 1, initially the outflow is caused by the centrifugal and magnetic forces on both sides (as in Proga 2005). However, this collimated outflow does not completely prevent the low-l matter inflow toward the BH from the polar regions. Therefore, remains always high (close the Bondi rate). In the lower panels in Figs. 4 and 5, the dotted contours in the angular momentum panels indicate that the high-l matter is separated from the outflow by the low-l accreting matter marked as the solid contours. The low- corona is separated from the low- collimated outflow with the high- region of the accreting matter. As we mentioned in the previous section, the moment of the outflow formation cannot be clearly identified in the mass accretion rate curve (Fig. 1). The persistently accreting low-l matter separates the polar outflow and magnetic corona for most of the simulation time. In this case, the accretion flow consists of four parts (from the equator to the pole): (1) geometrically thin torus in the equatorial plane, (2) highly turbulent magnetic corona, (3) low-l matter inflow through polar funnels, and (4) polar outflow.
In the intermediate case of =1.2 (run 3, third row of panels in Figs. 4 and 5), although the magnetic corona forms, the polar outflow forms in both manners as described for runs 1 and 2. In Fig. 4, one notices the corona formation on the N and S sides. The corona extends beyond , which is clearly seen in Fig. 4. First, on the S pole, the expanding corona causes the outflow. The outflow on the opposite pole, is caused by the magneto-centrifugal mechanism described in Proga (2005). The mass accretion rate, just as in the case of =1.01, does not experience a sudden reduction during its evolution because the outflow is weak so that the low-l matter persistently accretes and mainly contributes to the total mass accretion rate. Nevertheless, the mass accretion rate for =1.2 is slightly lower than that for =1.01.
In our models, the collimation of the jet-like outflow depends on the index and it is related to the mechanism of its formation rather than the torus geometrical thickness. For and (S side), the jet-like outflow at the final time of the simulation is not well collimated, because it is driven from the torus corona. To investigate the problem of collimation, one should perform the simulation over a longer time period to allow the flow to reach the quasistationary state at larger radii. For and (N side) the outflow is more collimated in comparison to the coronal jet-like outflow because it is formed by the low-l matter and collimated by the accreting low-l matter.
3.5 Episodic coronal outflows for low
For =1.01 (as well as for =1.2), remain high for most of the simulations time. However, at the late phases of the simulations, we observe a few significant minima in the mass accretion rate curves (lower panels in Fig.1). These minima are caused by small scale coronal outflows which temporary suppresses accretion of the low-l matter in the polar funel.
Fig. 6 presents the sequence of maps of the logarithmic density, angular momentum contours, and plasma parameter , before, during, and after an episode of an outflow in run 4. To better show the structure of the inner flow, we present only a half of the computational domain above the equator within 500 . Before the outflow, the magnetic torus corona extents to 100 (the two top rows of panels). During an outflow, the corona expands up to 300 (the third and fourth rows of panels). Later, the matter flows back toward the BH. Before and after the outflow, the torus corona and polar outflow are separated by a stream of the low-l matter. After the outburst, increases back to its mean high level (the bottom row of panels). The whole cycle lasts about 3500 .
3.6 The structure of the torii
In this section, we return to our investigation of the time evolution of for various . Here, we identify accretion states that correspond to various episodes of torus accretion. Figs. 7 shows the segments of the evolution for runs 1, 2, 3, and 4, which contain the typical ‘accretion states’ that appeared for a given model. The characteristic times, are marked with arrows denoted as S1, S2 (in run 1 and 2), S3 (only in run 2), S4 and S5 (in run 3 and 4).
We consider the flow to be unstable (MRI) when the following relation holds (Balbus & Hawley, 1998):
where is the Alfvén speed, is the angular velocity, is the distance from the BH and h is the half of the disk thickness. We define h as the height at which the density of the matter decreases by e along the vertical direction in comparison to the density at the equator. For run 1, 2, and 3, the disk height can be approximated with its radius h=r in Eq. 10. For run 4, we adopt the numerical value of ( roughly corresponds to 1/3 r). In the vertical direction, the torus typically contains and 24 grid points, for =5/3, 4/3, 1.2, and 1.01, respectively. In the radial direction, there are and equaly spaced in logarithmic scale zones within first 20 and 100 , respectively. To find, whether the fastest growing mode of the MRI is resolved with our computational mesh, we compute its critical wavelength as defined by eq. 113 in Balbus & Hawley (1998). We checked the ratio between and (the grid spacing) as a function of position on the grid. If the ratio is larger than one, the grid resolves torus for MRI. In regions where the ratio is smaller than one the grid does not resolve the critical wavelength .
We have analysed the flow properties at each accretion state as PB03b did. For =5/3, we found that the accretion state S1 corresponding to the narrow peak in Fig. 7 (panel for =5/3) appears when the high-l torus accretes on to the BH. We recognize S1 as ‘B’ state found in PB03b. During S1 the torus accretes because of the magneto-rotational instability (MRI). We checked that the current grid does resolve the fastest mode of the MRI. within . The accretion at the peak is dominated by the accretion from torus, ( reaches 0.06) and typically stays at this level for less than . S1 appears in the mass accretion rate curve many times (i.e. all spikes correspond to the torus accretion).
The accretion from torus is subsequently accompanied by the polar funnel inflow (‘leaking’) of the matter that was left in the polar regions after the outflow. The density in the polar funnel is relatively low in comparison to the torus density (or in comparison to the low-l matter removed from the polar funnel by the outflow). Its value is artificially set by the assumed parameter (see Sect. 2.2). The matter from the polar funnel contributes to the mass accretion rate and sets its mean level around 0.02 (Fig. 7).
During state S2 in run 1, is very low and it corresponds to a state when the torus is truncated from the BH. The truncation is due to magnetic effects (ie., a ‘magnetic barrier’) as in ‘C’ state described in PB03b (see also Proga & Zhang 2006 and fig. 1 there). We find that the poloidal magnetic field is larger than the toroidal component for during this state. The magnetic field lines are stretched for . The MRI does not operate in the inner . The matter from the polar funnel also flows outward during state S2. Therefore, drops much below the mean value of 0.02. The duration of S2 state is comparable to the duration of S1 and all moments with correspond to the same S2 state.
In our run 1, we do not see an accretion state corresponding to the ‘D’ state found in PB03b. In the ‘D’ state, should increase above 0.1 level, because of the slowly rotating matter previously pushed away by the coronal outflow, reentering the BH vicinity. We attribute lack of the ’D’ state in run 1 to a relatively short duration of our simulations. We performed our simulation for much shorter time (Table 1) in comparison to PB03b (2.5 ). In our simulations, the low-l matter is far away () from the center at the end of the simulation.
For =4/3, S1 state corresponds to the episode of torus accretion, but has a much lower amplitude than narrow peaks found for =5/3 ( S1 is representative of all narrow peaks with amplitude of in Fig. 7, panel for =4/3). We find that MRI causes the torus accretion, but our grid resolution was not adequate to capture the fastest growing mode of MRI, within . Thus, the mass accretion rate from the torus is lower than it would be if the fastest mode of MRI were resolved. In S1 state, the torus accretion is accompanied by accretion of the matter with density of from the polar region. This polar funnel accretion sets the mean accretion rate at the level of 0.01, similar to the mean accretion rate obtained for =5/3.
During S2 state (Fig. 7, panel for =4/3), the torus is pushed away from the BH by the magnetic barrier (as in S2 state for =5/3 case).
The mentioned ‘D’ state appears in run 2 despite of a short duration of the simulation. S3 state (Fig. 7, panel for =4/3) is characterized by the the sudden increase of up to the level of and much longer duration compared to state marked as S1. In Fig. 7, panel for =4/3, one can see only one state S3, but in Fig.1 one can find two broad features corresponding to state S3. During state S3, the low-l matter accretes on to the BH near the equator. This matter has been pushed away from the BH by the expanding magnetic corona during the bipolar outflow formation but not as far as in run 1. Therefore the weakly rotating matter may reenter the region close to the BH in a much shorter time-scale.
In the intermediate case of =1.2 (run 3), the mass accretion rate remains relatively high during the whole simulation (Fig. 7, panel for =1.2) due to almost persistent low-l polar radial inflow. Therefore, contrary to run 1 and 2, it is difficult to distinguish any characteristic patterns in the mass accretion rate curve. S4 marks peaking at 0.2, when the torus accretes along with the low-l matter. However, in this case the grid does not resolve the fastest growing mode of the MRI for . The torus accretion is dominated by the slower MRI modes. The main contribution to comes from a polar inflow of the low-l matter. The mass accretion rate is slightly higher than average because the matter from polar regions simultaneously accretes from both, N and S sides (S4). State S4 appears only a few times during run 3. For most of the time, the matter in the polar regions accretes from one side only. S4 is similar to the S3 in run 2, but during S3 low-l matter approaches the BH in the equatorial plane rather than from the poles, like in S4.
In S5 state, the mass accretion rate is dominated mainly by the torus, because almost all matter is removed from the polar funnel due to the episodic coronal outburst described in § 3.5. At that time becomes lower than 0.01. In run 3, the high-l torus never undergoes the truncation because the poloidal magnetic field is very weak and does not form a barrier for accretion. For S5 appears two times at the end of the simulation. Note that S5 is similar to state S1 in run 1 and 2, but here the torus accretes inefficiently, therefore S5 appears as a dip (not as a spike like in e.g. run 1) above the mean level of accretion.
For =1.01 (run 4), remains relatively high (Fig. 7, panel for =1.01) for most of the duration of the simulation. The low-l matter accretes quite persistently (as in run 3) along the torus surface (S4), on its N or S sides, or both. Although MRI operates in the torus, again the accretion fastest mode is not resolved with the simulation. Thus, the torus accretion rate is underestimated. In run 4, the torus never experiences a truncation. The low-l matter inflow is only occasionally stopped by the turbulent corona but such episodes are short since the outflow is weak (S5).
4 Discussion and Conclusions
In this article we present results from the axisymmetric 2-D MHD simulations of a low angular momentum accretion flow on to the BH. Our work is an extension of MP08, which in turn was initiated by PB03a and PB03b (see also Proga et al. 2003). Here, we concentrate on investigation of effects of gas pressure on the flow structure and evolution. In particular, we study how the flow properties change as a function of the adiabatic index, . We also compare our MHD runs to their HD counterparts presented in MP08.
We find that in the early phase of the evolution, the MHD runs behave very similar to those in MP08. Namely, the matter with forms a torus near the equator. As expected, the torus thickness decreasing with decreasing . The similarity in the MHD and HD initial evolution, is a consequence of our assumption of a initially very weak magnetic field. At later times, however, the MHD torii do not persistently accumulate matter near the BH as their HD counterparts do but rather accrete on the BH due to MRI and magnetic braking. In addition, the MHD torii are far more turbulent and also slightly geometrically thinner than their HD counterparts. However, most importantly contrary to the HD torii, the MHD torii produce bipolar outflows.
In MP08, we noticed that for , the torus is excited by the inflowing material and shows small amplitude quasi-periodic oscillations in the direction perpendicular to the equator. In the MHD runs, we do not observe clear torus oscillations in the power spectrum of the mass accretion rate or in any other way . This difference is due to the fact that in the MHD simulations the high-l torus is more turbulent and does not have a well defined outer boundary what would separate it clearly from the inflow (such a boundary is very clear in the HD simulations - see the bottom row of panels fig. 9 in MP08).
One of the properties of the innermost torus that is common for the HD and MHD simulations is the cylindrical distribution of the angular velocity in torus for all .
The bipolar outflow and angular momentum transport by magnetic effects transforms the structure of the flow from a polar-inflow/equatorial outflow (or equatorial accumulation) in the HD models into a polar-outflow/equatorial accretion in the MHD models. However, the torus outflow is strong only in models with =5/3, 4/3 and 1.2. For =1.01, the outflow is very weak and cannot propagate throughout the inflow.
In the MHD runs with a strong torus outflow, the total mass accretion rate of the gas on to the BH is reduced because is dominated by the polar funnel accretion of the low-l gas which the torus outflow suppresses. The mass accretion rate is also more variable in the MHD simulations than in the HD ones. This variability is caused mainly by mixing of the low- and high-l gas. As we showed in MP08, the degree of the mixing increases with the distances from the center for =4/3 and 1.2. The mixing causes episodes of an enhanced polar funnel accretion in the HD models (bursts). In the MHD models, we observe similar bursts not only for =4/3 and 1.2 but also for =5/3. In both HD and MHD models, can increase by a factor of few. For =1.2 in the MHD run, we do not observe very strong flares seen in the HD counterpart.
Our new simulations for =5/3 are consistent with the previous results obtained by PB03b. We found however, some minor differences between PB03b and our simulations. These differences occur during the initial phase of the evolution which we attribute to the differences in the initial conditions described in § 2. In addition, we do not observe here the episodes of enhanced accretion of the low-l matter (call by PB03b the accretion state D). This difference is due to a shorter duration of our simulations compared to PB03b’s simulations (see the previous section).
The torus outflow is not the only outflow we find in our simulations, i.e., we see also a highly collimated outflow driven by the centrifugal and magnetic forces. As found by (Proga, 2005), this outflow appears regardless of the high-l torus formation and is produced from the low-l inflow (see Proga et al. 2003). The ‘outflow-from-inflow’ does not significantly influence the torus evolution. If the matter in the polar region were not rotating (e.g. a distribution of l were a step function at the outer boundary) this outflow would not form as it was the case in PB03b’s simulations.
The most important new result in the present simulations is that for a geometrically thin disc model () the torus outflow is strongly suppressed. Generally, the structure of the torus and its magnetic corona depend on . For =5/3, gas pressure is significant and the flow is hot and the torus corona extends far from the equator. For =1.01, the magnetic corona above the disc is highly turbulent, and it is small and gravitationally bound. The magnetic field in the corona is less tangled in the flows for than the flow with =1.01.
As mentioned above, the torus outflow changes the mass accretion rate. For =5/3, dropped suddenly about one order of magnitude after the outflow develops. This drop in the accretion rate is due to depleting of the polar funnel. We observe the same behaviour of in the run with =4.3. For =1.2, the outflow is significantly weaker and is often one-sided. For =1.01, the torus outflow is very weak and is stalled by the inflow. The mass accretion rate remains high for the most of the simulation time in comparison to PB03b runs and our runs for =5/3 and 4/3.
In runs 1 and 2, we observe the accretion state ’C’ which occurs when the torus is truncated by the magnetic barrier. This state was first identified by PB03b. However, in the runs 3 and 4 (=1.2 and 1.01) we observe the stationary torus accretion (called by PB03b as state B) with the state C never occurring. This indicates that for low , the role of the magnetic field is reduced not only in producing a torus outflow but also in changing the torus accretion. We note that the reduced role of magnetic fileds may be also due to resolution of our models. The MRI cannot be resolved with points on mesh for , because the disk thickness decreases with (half of the disk thickness for and 1.2 and for ). We plan perform simulations with higher resolution in the near future.
The direct comparison of the simulations with the observations of BHB and AGN spectra are beyond the scope of this work. Therefore, we can only make some remarks on the observational consequences of our models. For example, the strong torus outflow for and the suppression of the outflow =1.01 may be relevant to the observed transition from the hard/low state to the soft/high state. The episodic outburst found in models with =1.2 and 1.01, may correspond with the flares seen in a very high states in BHB. To relate our results with observations, we plan to carry out calculations of broad-band spectra predicted by our models following Mościbrodzka et al. (2007).
Our models are axisymmetric, which may affect the results. In axisymmetry we cannot model the toroidal field MRI and a poloidal magnetic field can dissipate according to the antidynamo theorem (Moffatt, 1978). We note that even in the HD limit, the axisymmetric models differ from non-axisymmetric ones. For example, Janiuk et al. (2008) found that in fully 3-D simulations the inner torus processes. It would be important to check if the torus will also precess in the MHD limit and what consequences it would have on the outflow formation and the polar-funnel accretion. Recent work on outflows driven by radiation from a precessing disc show that the flow geometry can be dramatically changes although the mass accretion rate can be similar to that in the non-precessing disc case (Kurosawa & Proga, 2008). We plan to address these issues in our future work.
Future simulations should also include radiative cooling and heating processes. Recently, Fragile & Meier (2009) pointed out that the radiative cooling may be of the same order as artificial cooling caused by lack of the treatment of kinetic and magnetic dissipative processes in non-conservative codes (as one used by us). Our models also do not resolve the geometrically thin discs. Currently, global simulations of the magnetized geometrically thin disc is beyond computational abilities.
We acknowledge support provided by the Chandra awards TM7-8008X (M.M. and D.P.) and TM8-9004X (D.P.) issued by the Chandra X-Ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS 8-39073.
- Abramowicz et al. (2000) Abramowicz M. A., Björnsson G., Igumenshchev I. V., 2000, PASJ, 52, 295
- Balbus & Hawley (1991) Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214
- Balbus & Hawley (1998) Balbus S. A., Hawley J. F., 1998, Reviews of Modern Physics, 70, 1
- Corbel et al. (2001) Corbel S., Kaaret P., Jain R. K., Bailyn C. D., Fender R. P., Tomsick J. A., Kalemci E., McIntyre V., Campbell-Wilson D., Miller J. M., McCollough M. L., 2001, ApJ, 554, 43
- De Villiers et al. (2005) De Villiers J.-P., Hawley J. F., Krolik J. H., Hirose S., 2005, ApJ, 620, 878
- Done et al. (2007) Done C., Gierliński M., Kubota A., 2007, A&AR, 15, 1
- Fender et al. (2004) Fender R. P., Belloni T. M., Gallo E., 2004, MNRAS, 355, 1105
- Fragile & Meier (2009) Fragile P. C., Meier D. L., 2009, ApJ, 693, 771
- Gallo et al. (2003) Gallo E., Fender R. P., Pooley G. G., 2003, MNRAS, 344, 60
- Hawley & Balbus (2002) Hawley J. F., Balbus S. A., 2002, ApJ, 573, 738
- Heinz & Sunyaev (2003) Heinz S., Sunyaev R. A., 2003, MNRAS, 343, L59
- Ho (2003) Ho L. C., 2003, in Collin S., Combes F., Shlosman I., eds, Active Galactic Nuclei: From Central Engine to Host Galaxy Vol. 290 of Astronomical Society of the Pacific Conference Series, The Central Engines of Low-Luminosity AGNs. p. 379
- Ichimaru (1977) Ichimaru S., 1977, ApJ, 214, 840
- Janiuk et al. (2008) Janiuk A., Proga D., Kurosawa R., 2008, ApJ, 681, 58
- Kato et al. (2004) Kato Y., Mineshige S., Shibata K., 2004, ApJ, 605, 307
- Kawakatu et al. (2009) Kawakatu N., Nagao T., Woo J.-H., 2009, ApJ, 693, 1686
- Königl (2006) Königl A., 2006, Memorie della Societa Astronomica Italiana, 77, 598
- Kurosawa & Proga (2008) Kurosawa R., Proga D., 2008, ApJ, 674, 97
- Liu et al. (2002) Liu B. F., Mineshige S., Meyer F., Meyer-Hofmeister E., Kawaguchi T., 2002, ApJ, 575, 117
- Liu et al. (1999) Liu B. F., Yuan W., Meyer F., Meyer-Hofmeister E., Xie G. Z., 1999, ApJ, 527, L17
- Livio et al. (1999) Livio M., Ogilvie G. I., Pringle J. E., 1999, ApJ, 512, 100
- Markoff (2006) Markoff S., 2006, in Kannappan S. J., Redfield S., Kessler-Silacci J. E., Landriau M., Drory N., eds, New Horizons in Astronomy: Frank N. Bash Symposium Vol. 352 of Astronomical Society of the Pacific Conference Series, Accretion and Jets in Microquasars and Active Galactic Nuclei. p. 129
- McKinney & Gammie (2004) McKinney J. C., Gammie C. F., 2004, ApJ, 611, 977
- Meier (2001) Meier D. L., 2001, ApJ, 548, L9
- Meier (2005) Meier D. L., 2005, Ap&SS, 300, 55
- Meyer et al. (2000) Meyer F., Liu B. F., Meyer-Hofmeister E., 2000, A&A, 361, 175
- Meyer et al. (2007) Meyer F., Liu B. F., Meyer-Hofmeister E., 2007, A&A, 463, 1
- Mizuno et al. (2004) Mizuno Y., Yamada S., Koide S., Shibata K., 2004, ApJ, 615, 389
- Moffatt (1978) Moffatt H. K., 1978, Magnetic field generation in electrically conducting fluids. Cambridge, England, Cambridge University Press, 1978. 353 p.
- Mościbrodzka & Proga (2008) Mościbrodzka M., Proga D., 2008, ApJ, 679, 626(MP08)
- Mościbrodzka et al. (2007) Mościbrodzka M., Proga D., Czerny B., Siemiginowska A., 2007, A&A, 474, 1
- Narayan et al. (1998) Narayan R., Mahadevan R., Grindlay J. E., Popham R. G., Gammie C., 1998, ApJ, 492, 554
- Paczyński & Wiita (1980) Paczyński B., Wiita P. J., 1980, A&A, 88, 23
- Proga (2005) Proga D., 2005, ApJ, 629, 397
- Proga & Begelman (2003a) Proga D., Begelman M., 2003a, ApJ, 582, 69
- Proga & Begelman (2003b) Proga D., Begelman M., 2003b, ApJ, 592, 767(PB03b)
- Proga et al. (2003) Proga D., MacFadyen A. I., Armitage P. J., Begelman M. C., 2003, ApJ, 599, L5
- Proga & Zhang (2006) Proga D., Zhang B., 2006, MNRAS, 370, L61
- Remillard & McClintock (2006) Remillard R. A., McClintock J. E., 2006, ARA&A, 44, 49
- Różańska & Czerny (2000) Różańska A., Czerny B., 2000, A&A, 360, 1170
- Shakura & Syunyaev (1973) Shakura N. I., Syunyaev R. A., 1973, A&A, 24, 337
- Stone & Norman (1992a) Stone J. M., Norman M. L., 1992a, ApJS, 80, 753
- Stone & Norman (1992b) Stone J. M., Norman M. L., 1992b, ApJS, 80, 791
- Stone & Norman (1994) Stone J. M., Norman M. L., 1994, ApJ, 433, 746
- Stone & Pringle (2001) Stone J. M., Pringle J. E., 2001, MNRAS, 322, 461
- Uchida & Shibata (1985) Uchida Y., Shibata K., 1985, PASJ, 37, 515
- Yuan et al. (2003) Yuan F., Quataert E., Narayan R., 2003, ApJ, 598, 301
- Zdziarski et al. (1995) Zdziarski A. A., Johnson W. N., Done C., Smith D., McNaron-Brown K., 1995, ApJ, 438, L63