{'fraction_non_alphanumeric': 0.03838489774514945, 'fraction_numerical': 0.016081104701975178, 'mean_word_length': 4.658951533135509, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Beyond the standard accretion disc model: coupled magnetic disc--corona solutions with a physically motivated viscosity law

We present a systematic, analytical study of geometrically thin, optically thick accretion disc solutions for magnetized turbulent flows, with an alpha-like viscosity prescription. Under the only assumptions that (1) Magneto-Rotational instability (MRI) generates the turbulence that produces the anomalous viscosity needed for accretion to proceed, and that (2) the magnetic field amplified by the instability saturates due to buoyant vertical escape, we are able to self-consistently solve the disc structure equations including the fraction of power, f, that is carried off by vertical Poynting flux (and likely dissipated outside the optically thick disc). For low-viscosity discs, we obtain stable high-f solutions at low accretion rates, when gas pressure dominates, and unstable, low-f, radiation pressure dominated solutions at high accretion rates. For high viscosity discs, instead, a new thermally and viscously stable, radiation pressure dominated solution is found, characterized by f~1 and appearing only above a critical accretion rate (of the order of few tenths of the Eddington one). We discuss the regimes of validity of our assumptions, and the astrophysical relevance of our solutions. We conclude that our newly discovered thin disc solutions, possibly accompanied by powerful, magnetically dominated coronae and outflows, should be seriously considered as models for black holes accreting at super-Eddington rates.


INTRODUCTION
The so-called standard model for accretion discs (Shakura & Sunyaev 1973) describes an optically thick, geometrically thin solution to the accretion equations in which all our ignorance about the mechanisms of angular momentum transport and turbulent viscosity is side-stepped by making the (physically reasonable) Ansatz that the turbulent stress tensor scales with the local pressure, with a constant of proportionality that we will denote here αSS. Such an approach proved enormously productive, and despite the great development of numerical models in the last decade, standard disc modeling still remains at the very heart of accretion disc phenomenology and the basic link between theory and observations. The recent evolution of numerical studies, instead, has proved itself fundamental to shed light on the very nature of the αSS prescription, by elucidating the crucial role of MHD turbulence for the enhanced transport properties of accretion discs . In particular, non-radiative flows are ideal objects of study by means of global three dimensional simulations (see  and reference therein).
On the other hand, the observed spectrum of thin accretion discs around many compact objects suggests that a corona is also present (Nandra & Pounds 1994). Such coronae are most likely a by-product of the internal disc dynamics (as in the model originally proposed by Galeev, Rosner & Vaiana 1979) and local 3-D MHD simulations of stratified discs (Miller & Stome 2000), although adopting a thermodynamics which is most likely inappropriate for the problem at hand, may support this idea.
The main aim of the present work is to incorporate, in the simplest possible way, our current knowledge on the properties of turbulent MHD flows in the standard accretion disc theory. This will allow us, on the one hand, to extend the predictive power of the standard model by self-consistently including in the theory the magnetic properties of the flow and the generation of a magnetic corona by means of vertical Poynting flux; and on the other hand, to consider MHD turbulent discs in a regime where global, 3-D numerical simulations are currently unable to make predictions, i.e. for highly luminous, radiatively efficient thin discs.
Thus, we assume the accretion proceeds through a geometrically thin disc,à la Shakura & Sunyaev (1973), and we employ vertically integrated equations. Furthermore, we assume that magnetorotational instability (MRI; see , and reference therein) is the primary source of the turbulent viscosity and, consequently, that the turbulent magnetic stresses responsible for the angular momentum transport in the disc scale with magnetic pressure. Then, the accretion disc structure is fully described once the relationship between magnetic and disc pressure (given either by gas or radiation) is established.
The issue of the viscous scaling can be translated into the uncertainty about the mechanisms by, and the level at, which the disc magnetic field saturates. In the following, we make the working as-sumption (in a sense justified a priori by the need of explaining the observed presence of strong coronae in many accretion powered systems) that magnetic field is disposed of mainly by vertical buoyancy (Stella & Rosner 1984;Sakimoto & Coroniti 1989), and verify a posteriori its validity. Then we will show how under this condition, and taking into account the property that MRI grows less rapidly in radiation dominated discs, the scaling of magnetic pressure can be determined. This in turn specifies the closure relation that fully describes a coupled accretion disc-corona system.
Furthermore, we will show that by solving the structure equations a new solution is found, which has no analog in the standard Shakura-Sunyaev theory, and is relevant for highly viscous (and highly magnetized) flows accreting at a rate of the order of (or above) the Eddington one.

BASIC EQUATIONS
Numerical simulations of MHD turbulent discs clearly show that the non-linear outcome of the MRI is a fully three dimensional, anisotropic turbulence, that exhibits strongly correlated fluctuations in the azimuthal and radial components of velocity (Reynolds stresses) and magnetic field (Maxwell stress). In accordance with those results, and for the sake of simplicity, we will assume that the vertically averaged anomalous stress is dominated by Maxwell stresses, and therefore where Pmag = B 2 8π is the magnetic pressure, k0 is a factor of order unity that depends on the relative intensity of radial and azimuthal field components 1 .
This in turn implies that the local heating rate in a thin Keplerian disc is given by where ΩK Keplerian angular velocity, cs = Ptot/ρ is the isothermal sound speed.
The growth rate of MRI is influenced by the ratio of the gas to magnetic pressure, as demonstrated both by analytical studies of the linear regime of the instability (Blaes & Socrates 2001) and by numerical MHD simulations of radiation pressure dominated discs (Turner, Stone & Sano 2002): both found that, due to the substantial compressibility of MHD turbulence in radiation pressure dominated discs, where cg = Pgas/ρ is the gas sound speed. The magnetic field, with initial (highly subthermal) amplitude B0, rapidly grows to B ∼ B0e σt b , where H is the total disc scaleheight and is the mean buoyant rise-time. In order to estimate the saturation field, an estimate for the upward drift velocity vD caused by buoyancy of magnetic field is needed. In the case of uniform discs permeated by a net mean field, where the magnetic field can be modeled as a collection of flux tubes that retain their individuality over time scales much longer than dynamical time, it has been shown (Stella & Rosner 1984) that vD ≃ vA(cg/cs)(π/2CD) 1/2 (a/H) 1/2 . Here vA = |B|(4πρ) −1/2 is the Alfvén speed, a is a flux tube cross section and CD a drag coefficient, whose numerical value can be estimated to be ∼ 1-10. If the tube cross section a were independent on the ratio of gas to radiation pressure, such an expression for the upward drift velocity would lead to a final scaling of magnetic pressure with the gas pressure only, as already argued in Stella & Rosner (1984). However, the applicability of a flux tube analysis to the case of MRI turbulent flows is not justified. The work of Stone et al. (1996) clearly demonstrates that magnetic concentrations, although always present, do not persist as coherent structures. Moreover, as already mentioned, the turbulence in radiation pressure dominated discs is compressible. Then, in the vertical fluid momentum equation any vertical gradient of the magnetic pressure will not be balanced by either radiation (because radiation is diffusive) nor gas pressure (because azimuthal field energy density may exceed the true gas pressure, see below). Vertical motions are excited, with approximately (see e.g. Eq. 4.19 in Blaes 2002) which would imply (for a mainly azimuthal field) a direct proportionality between the vertical drift velocity and the Alfvèn speed. We therefore assume vD = bvA, with b constant of the order of unity. If we then put together Eqs.
(3), (4) and (6), we have that Pmag ln(B/B0) ∝ PtotPgas, where Ptot = Pgas + P rad is the sum of gas plus radiation pressure. Neglecting the logarithmic dependence on the initial field, this in turn suggests that, in MRI dominated turbulent flows, the field will saturate at an amplitude such that (Taam & Lin 1984;Burm 1985;Szuszkiewicz 1990;Merloni & Fabian 2002a): where we have introduced the constant where β = Ptot/Pmag is the usual plasma parameter. It is worth noting here the relationship between the constant α0 and the Shakura-Sunyaev coefficient αSS ≡ Pmag/Ptot (which is not constant in our theory). We have The condition that the turbulence be subsonic (vA < ∼ cs) implies αSS < ∼ 1/2 (or, equivalently, β > ∼ 2), and from (9) we see that α0 can still be larger than unity, when radiation pressure dominates in the disc.
The magnetic flux escaping in the vertical direction may dissipate a substantial fraction of the gravitational binding energy of the accreting gas outside the optically thick disc, with obvious deep implications for the spectrum of the emerging radiation (Haardt & Maraschi 1991). How much of this flux really escapes into the low-density environment above and below the disc (in the so-called corona), depends on the details of vertical transport of magnetic flux tubes and vertical stratification of the disc, which is ignored in our model. Nonetheless, we can safely assume that the fraction of power dissipated in the corona is determined by the ratio f of the vertical Poynting flux to the local heating rate Q+ (Svensson & Zdziarski 1994). The simplest way to estimate the vertical Poynting flux is to just assume that This translates into (see Eq. 2) In the following, we will assume that the numerical factor (of the order of unity) k1 = 3k0/2b = 1. Differences in its exact value will affect the numerical value of the critical viscosity parameter α0 (see below), but won't affect the nature of the solutions we find.
Thanks to the closure relation (11), we are now able to build a self-consistent accretion disc-corona solution. In order to do so, exactly as in the case of the standard Shakura-Sunyaev solution, we have to solve simultaneously four conservation equations (conservation of vertical momentum, mass, angular momentum and energy) and supplement them with the equations of state for the pressure and opacity as functions of density and temperature. The only modification with respect to the standard derivation is the reduction of the locally dissipated energy by the factor (1 − f ) (Svensson & Zdziarski 1994), with f given by Eq. (11), and the adoption of the viscosity law (7).
For simplicity, we present here the solutions in Newtonian approximation, as in Shakura & Sunyaev (1973), with an efficiency ǫ0 = 0.082 appropriate for non-rotating black holes (a full general relativistic analog of the equations below can be found in the appendix of Merloni & Fabian 2002b). We concentrate only on the innermost regions of the disc, where opacity is dominated by electron scattering (κes = 0.4 cm 2 g −1 ). Introducing dimensionless units for the central mass (m ≡ M/M⊙), accretion rate (ṁ ≡ ǫ0Ṁ c 2 /L Edd ) and distance from the center (r ≡ R/RS = Rc 2 /2GM ), we obtain, for the radiation pressure dominated solution, the following expressions for disc scaleheight, pressure (in dyne), density (in g cm −3 ), both calculated at the disc mid-plane, and temperature (in K): supplemented with the closure equation for f : Analogously, in the gas pressure dominated part of the disc, we have: (1 − f ) 9/10 = 4.7 × 10 2 (α0m) 1/10 r −21/20 (ṁJ(r)) 4/5 . (15) In the above formulae, the function J(r) = 1 − rin/r, with rin = 3 corresponding to the innermost stable circular orbit of a free particle around a non-rotating black hole, describe the notorque at the inner boundary condition.

STUDY OF THE SOLUTIONS
To find all the family of possible solutions, we proceed in the following way: for each set of the parameters (m,ṁ, r and α0) we find the roots of both Eq. (13) and (15) in the interval 0 < f < 1. We then substitute the obtained values of f into the appropriate expressions for the gas and radiation pressure and finally discard all the non-consistent solution (i.e. those found with the formula for gas pressure dominated discs that gives P rad > Pgas and vice versa).
The value of the viscosity parameter α0 = k1/2 = 0.5 is a critical point for the set of equations describing the coupled disc-corona system, in that it separates two qualitatively different regimes. In the next two sections we will examine these two regimes, beginning with the low viscosity case.

α0 < 1/2: an almost standard solution
It turns out that for α0 < 1/2 each set of parameters define a unique solution. Its properties have already been discussed (Merloni & Fabian 2002a), but we present them here in a more general framework. The fraction of disc power dissipated in the corona, f , tends to its maximum value, √ 2α0, when gas pressure dominates (low accretion rates), and decreases as the accretion rate increases and radiation pressure becomes more and more important. This is shown in Fig. 1, where, for m = 10 and α0 = 0.25, we plot the radial dependence of the coronal fraction for different values of the accretion rate.
As in the standard case, radiation pressure dominated part of the disc are thermally and viscously unstable, although the combined effects of the modified viscosity law (7) and the stabilizing effect of the corona, reduce the extent of the unstable region of the parameter space (Szuszkiewicz 1990).

α0 > 1/2: a new stable solution at high accretion rates
The nature of the solutions changes qualitatively if α0 > 1/2. Because in gas pressure dominated parts of the disc f ∼ √ 2α0, there can be no solutions there with f < 1. That is to say that the condition α0 < 1/2 should hold in gas pressure dominated discs, and these high viscosity solutions can only exist if radiation pressure dominates.
For each value of α0 > 1/2 the accretion rateṁ(r, f ) (Eq. 13) has a minimum for f ∈ (0, 1). Thus, solutions exists only for accretion rates larger that a critical value,ṁcr(r) = min f ∈(0,1) [ṁ(f, r)] ≡ṁ(r, fcr). This corresponds to a critical value, fcr, of the fraction of power vertically transported by Poynting flux. Forṁ(r) >ṁcr(r) two radiation pressure dominated solutions are found for each set of parameters, one with f < fcr and the other with f > fcr (see Fig. 2). The first is a "standard" (unstable) one, characterized by a low value of the fraction of power  transported vertically by Poynting flux, and is just the continuation of the radiation pressure dominated branch discussed in the previous section at high values of α0. The second one, which is discussed here for the first time, appears due to the feedback effect of the closure relation (11). It is a (mildly) radiation pressure dominated solution endowed with a powerful corona, in that f → 1 as the accretion rate increases.
The values ofṁcr(r) and fcr can be found by solving ∂ṁ(r, f )/∂f = 0, or, equivalently, the equation f 5 − 12α 2 0 f + 8α 2 0 = 0. For every r, re obtain fcr ≃ 0.7417 for α0 = 0.5, quickly converging towards fcr = 2/3 for increasing α0. The critical accretion rate, though, depends on r. The lowest possible value of the critical accretion rate, that can be simply denoteḋ mcr, appears at r ≃ 4.64 (which is the radius where the function J(r) 16/9 r −7/3 has a minimum). We have numerically computeḋ mcr as a function of α0. The relation can be approximated by the analytic expressioṅ which is accurate to better than 3% for α0 > 1 and to better than 10% for α0 > 0.6. The above formula also makes manifest the main differences that would be found if we chose supermassive black holes instead that stellar mass ones. Namely, the critical accretion rate separating gas from radiation pressure dominated solutions and that determining the existence of our new high viscosity, high-f solution would be reduced by a factor m −1/8 . The qualitative properties of the solutions would be left unchanged. Fig. 2 shows, for m = 10 and α0 = 0.75, the radial dependence of the coronal fraction for different values of the accretion rate. The two solutions appear forṁ > ∼ 0.19 and show an opposite trend as the accretion rate increases: while in the first case f decreases and the disc becomes less dense and more radiation pressure dominated, in the second case f increases, the disc becomes denser and less radiation pressure dominated.
It is worth studying the stability property of this newly discovered solution. Thermal stability requires (Piran 1978) We have (∂ ln Q−/∂ ln T ) = 4 and, from (12), Q+ ∝ṁ(1−f ) ∝ T 9/2 (1 − f ). Therefore, the stability criterion becomes: Because the numerator of the second term is always positive for α0 > 1/2, we have that for f > fcr (which is a pole for the second term in the above equation) indeed the criterion is satisfied, while the opposite is true for f < fcr. On the other hand, viscous stability requires that (∂Ṁ /∂Σ) > 0. We have Σ ∝ ρH ∝ṁ −1/9 (1 − f ) −1 , and the stability condition can be written as which is satisfied when thermal stability condition is also satisfied. The very important fact that our new discovered high-f solutions are both viscously and thermally stable may be an indication that, whenever allowed to do so, the system will chose those solutions over the low-f ones. The real physical state of a high viscosity disc at high accretion rates, though, is likely more complicated and inherently time dependent. Indeed, the very fact that for α0 > 0.5 no gas pressure dominated solutions exist, make us believe that in the outer part of the disc magnetic turbulence will self-regulate to keep the viscosity parameter below its critical value. Thus, the value of the viscosity parameter may not be constant throughout the disc, but more likely increases inwards, as also suggested by high-resolution numerical MHD simulations of the inner parts of non radiative discs (Reynolds & Armitage 2001; Krolik . For any given radial profile of the viscosity parameter α0, the radius at which it reaches its critical value will therefore play an important role in determining the physical properties of the disc. We may speculate that in a full time dependent solution both the viscosity parameter and the local accretion rate will adjust in order to avoid discontinuities in the flow there, but the exact nature of the final disc structure is impossible to predict from our stationary solutions.

DISCUSSION
Beside the thin disc approximation, which is straightforward to verify a posteriori for all our solutions but the low-f , high-ṁ ones (as in the standard SS case), our results are rooted in two further assumptions. The first is that the growth rate of MRI is given by (3) in radiation pressure dominated discs, the second is that buoyancy is the main mechanism of field saturation in turbulent magnetized discs. Let us examine these two issues one at a time.
As discussed in Blaes & Socrates (2001), the applicability of Eq. (3) depends on the capability of radiative diffusion to destroy radiation sound waves on the length scale of the MRI, and requires that the azimuthal magnetic field component dominates the vertical one, as should be the case in thin discs with Keplerian velocity profiles (see also the numerical simulations of initially random fields in . Let us define the Alfvén speed associated to the vertical field component, vAz = Bz(4πρ) −1/2 , and the (small) parameter δ ≡ vAz/vA. The diffusion wavenumber for Keplerian discs is (Blaes & Socrates 2001): Then, Eq. (3) should be valid if the characteristic wavenumber of the instability, ΩK/vAz is of the order of (or larger than) the diffusion wavenumber k diff . By using the expression (12) for the density, (11), and the definition of vAz, we obtain the following condition for the validity of our assumption: Noting that max[J(r) 4/9 r −21/36 ] ≃ 0.2, this translates into 2 1.5f m 1/18 1 ṁ α0 which is always satisfied for the low-f radiation pressure dominated solutions, while, for the f ∼ 1 solutions, it breaks down for very high accretion ratesṁ ≫ 0.55α0m −1/8 δ −9/4 −1 ≃ 1.6ṁcrδ −9/4 −1 . It is clear that in order for our viscosity scaling (7) to be applicable to the high-f solutions at super-Eddington rates, higher azimuthal field strengths (and low δ) are required. However, ifṁ grows too large, and Eq. (22) is not satisfied, then the magnetic stresses, and the magnetic pressure, will be proportional to the total pressure. The condition of subsonic turbulence would then imply again α0 < 1/2. The fraction of power transported vertically by Poynting flux, f will tend again to its maximum value ∼ √ 2α0, as in the gas-pressure dominated solutions. Therefore, it seems plausible to envisage that the high-f solutions could extend to arbitrarily large accretion rates, regardless to the viscosity law and the detailed properties of the MRI field, as long as α0 is sufficiently high in these regimes.
The second main assumption that was used in the derivation 2 We adopt here the notation A = 10 x Ax.
of our solutions was to consider buoyancy the main mechanism by which the disc gets rid of the magnetic field. In fact, two other competing mechanisms should be considered: turbulent diffusion and dissipation. Let us first compare the typical buoyancy timescale (t b ≃ H/2bvA) with the turbulent diffusion one (t d ≃ H 2 /ν). The kinematic viscosity coefficient ν can be recovered from the standard viscous form of the stress to obtain ν = 2τ Rφ /(3ρΩK). We then have Therefore, for any physical self-consistent, f < 1 solution, buoyancy is always a more efficient mechanism of field escape than field diffusion, unless k1 = 3k0/2b is substantially larger than unity.. On the other hand, the field may dissipate by cascading down the subsonic turbulence and leave its energy inside the disc. The typical timescale for this process is set by the coherence time of the turbulent flow, that we have assumed to be τt ≡ 1/σ ∼ vA/(cgΩK). This would give which implies that our working assumption is justified only for large viscosity solutions. In practice, if α0 is small, the actual maximal fraction of power dissipated into the corona will be further reduced by a factor 1 + 0.25/bα0. However, the high viscosity, high-f solutions discussed in section 3.2 will be little affected by dissipational losses of energy inside the disc.

CONCLUSIONS
We have presented a full analysis of thin accretion disc solutions under the assumptions that the magnetic field generated and amplified by MRI (which is at the origin of MHD turbulent stress in the disc) saturates mainly due to buoyant transport in the vertical direction. From these assumptions, a viscosity prescription is derived, in which magnetic pressure, and therefore turbulent stresses, scale proportionally to the geometric mean of gas and total (gas plus radiation) pressure. This prescription, in turn, allows us to uniquely determine, for every set of the parameters (m,ṁ, r and α0) the fraction of the accretion power f which is transported vertically by Poynting flux. For low values of the viscosity parameter the solution is unique, with f tending to its maximum value, √ 2α0, when gas pressure dominates (low accretion rates), and decreasing as the accretion rate increases and radiation pressure becomes more and more important. For 2α0 > 1 there are no solutions for gas pressure dominated regions, while in the radiation dominated part of the flow two solutions are possible for every value of the accretion rate. The first is a "standard" (unstable) one, with f ≪ 1, while the second has f → 1 (corona dominated) as the accretion rate increases and is thermally and viscously stable.
Two main results from the present study are important for the interpretation of observed spectral energy distributions of accretion-powered systems. The first is the property of low-α0 solutions to be more corona-dominated (i.e. to have the largest possible value of f ) at low accretion rates. This has been discussed in Merloni & Fabian (2002a) in the context of low-luminosity black hole. It is indeed an observational fact that the relative strength of the hard, optically thin, power-law-like component in their X-ray spectra increases with respect to the optically thick, quasi-thermal one asṁ decreases. The second is the existence of a new highf solution at high accretion rates for high values of the viscosity parameter.
As in the case of corona-dominated solutions at lowṁ, it is plausible to speculate that the large Poynting fluxes in the vertical direction that characterize these solutions are to be associated with the generation of powerful coronae and, possibly, powerful outflows/jets, too (Merloni & Fabian 2002a). Indeed, there are observational hints that black holes believed to accrete at Eddington or super-Eddington rates bear many similarities with their low-rate counterparts. This class of objects may include galactic black holes in the Very High state (see Done 2002, and references therein), whose X-ray spectra show a power-law component, likely of non-thermal origin (Poutanen & Coppi 1998;Gierliński et al. 1999), often as strong as the quasi-thermal one; broad line radio galaxies, where the weakness of X-ray reflection features has been interpreted as evidence for highly ionized discs surrounded by powerful coronae at high accretion rates (Ballantyne, Ross & Fabian 2002); powerful blazars, where it has been shown (Ghisellini & Celotti 2002) how the relativistic jet in these sources may dominate the total output power at very high (inferred)ṁ. Furthermore, some Ultra Luminous X-ray sources (ULX; Makishima et al. 2000) have shown power-law like spectra, analogous to those of black hole candidates in their very high state, at inferred accretion rates of a few times the Eddington rate for a 30M⊙ black hole (Kubota, Done & Makishima 2002).
In order to assess the relevance of our new solution in all the above cases, vertically integrated models should be abandoned, and the complicated problem of determining the disc structure in the vertical direction, together with the disc-corona boundary, tackled. This would allow us to establish the relationships between the (vertically integrated) parameter f , and the observationally determined fraction of power dissipated in the corona. In particular, it is only by carefully studying the disc-corona transition region that will be possible to assess how super-Eddington (if at all) our new high α0 solutions can be. In principle, if a high fraction of accretion power could be dissipated in the corona without disrupting the optically thick disc, there would be no limitation to the total luminosity that the source could emit, because the radiation is shining on the disc from outside. Such a possibility to obtain super-Eddington luminosities is physically different from the one recently proposed by Begelman (2001;, although both depend on the existence of magnetized radiation pressure dominated thin discs. In fact, the non-linear analysis of Gammie (1998) and Blaes & Socrates' (2001) photon bubble instability carried on in Begelman (2001) is applicable to the low-f radiation dominated solutions we have presented here, which is much more strongly radiation pressure dominated than the high-f ones, while the thinner discs of the high-f solutions have almost equal radiation and gas pressures, and therefore should not be too strongly inhomogenous.
Finally, we point out that the existence of multiple solutions for the same set of parameters (high α0 case), may give rise to interesting hysteresis effects during spectral transition caused by long term variation of the accretion rate (Maccarone & Coppi 2003). A full time-dependent analysis of our solution is needed in order to make definite predictions in this sense.