text:
# Extending the Shakura-Sunyaev approach to a strongly magnetized accretion disc model

## 1 Introduction

The well known and most widely used model of the accretion disc was proposed and elaborated by Shakura (1972) and Shakura & Sunyaev (1973). In this model the disc is vertically supported by the thermal pressure (Shakura & Sunyaev, 1973). Turbulent viscosity is invoked in the Shakura–Sunyaev model to explain the angular momentum transfer required by the accretion flow. As originally pointed out in Lynden-Bell (1969) and Shakura & Sunyaev (1973) a magnetic field can also contribute to the angular momentum transport. A robust mechanism of the excitation of magnetohydrodynamical (MHD) turbulence was shown to operate in accretion discs due to the magneto-rotational (MRI) instability (Balbus & Hawley, 1998). The growth of the MRI leads to the excitation of turbulent magnetic fields and self-sustained MHD turbulence. The contribution of Maxwell stresses to the transport of angular momentum is usually larger than Reynolds stresses. However, the magnetic energy observed in many numerical experiments was smaller than the thermal energy of the gas in the disc (Brandenburg, 1998). Simulations of the non-linear stage of MRI are typically local simulations in a shearing box of an initially uniform small part of the disc. Attempts to expand the computational domain to include a wider area of radii and azimuthal angle (Hawley & Krolik, 2001; Hawley, 2001; Armitage, Reynolds & Chiang, 2001) are underway. However, even before the recent focus on the MRI Shibata, Tajima & Matsumoto (1990) observed the formation of transient low $`\beta `$ state in a shearing box simulations of the non-linear Parker instability in an accretion disc.

Vertical stratification has been considered in the shearing box approximation (Brandenburg et al., 1995; Miller & Stone, 2000). In particular, Miller & Stone (2000) investigated discs with initial Gaussian density profiles supported by thermal pressure. The initial seed magnetic field grows and starts to contribute to the vertical pressure balance. The computational domain extends over enough vertical scale heights to enable Miller & Stone (2000) to simulate the development of a magnetically dominated corona above the disc surface. In the case of an initial axial magnetic field, Miller & Stone (2000) observed that the saturated magnetic pressure dominates thermal pressure not only in the corona but everywhere in the disc. As a consequence, the thickness of the disc increases until it reaches the axial boundaries of the computational box. The formation of low $`\beta `$ filaments in magnetized tori was also observed in global MHD simulations by Machida, Hayashi & Matsumoto (2000). Although further global MHD simulations of vertically stratified accretion discs are needed, this numerical evidence suggests that magnetically dominated thin discs may exist.

Previously, analytic models of thin accretion discs with angular momentum transfer due to magnetic stresses were considered by Eardley & Lightman (1975) and Field & Rogers (1993a, b). Both these works included magnetic loops with size of the order of the disc thickness. In Eardley & Lightman (1975), the magnetic loops were confined to the disc. Loop stretching by differential rotation was balanced by reconnection. The reconnection speed was a fraction of the Alfvén speed. Radial magnetic flux was considered as a free function of the radius. Vertical equilibrium and heat transfer were treated as in Shakura & Sunyaev (1973), with the addition of the magnetic pressure in the vertical support. No self-consistent magnetically dominated solutions were found in model of Eardley & Lightman (1975).

In contrast, dominance of the magnetic pressure over the thermal and radiation pressure was postulated from the beginning by Field & Rogers (1993a, b) and verified at the end of their work. These authors assumed that the ordered magnetic field in the disc, amplified by differential rotation, emerges as loops above the surface of the disc due to Parker instability. Because the radial magnetic field in the disc has an intially sectorial structure, the loops above the disc come to close contact and reconnect. All dissipation of magnetic field occurs in the corona in the model of Field & Rogers (1993a, b). Such a corona was assumed to be consisting of electrons and some fraction of positrons and no outflow from the disc is present. Electrons and positrons are accelerated to relativistic energies at the reconnection sites in the disc corona and subsequently emit synchrotron and inverse Compton photons. Because reconnection was assumed to occur at loop tops, Field & Rogers (1993a, b) found that up to 70 per cent of the energy released in reconnection events in the corona will be deposited back to the surface of the disc in the form of relativistic particles and radiation. Only thin surface of optically thick disc is heated and cools by the thermal emission, which is the primary source of soft photons for the inverse Compton scattering by relativistic particles in the corona.

Since the characteristic velocity of rise of the loops of the buoyant magnetic field is of the order of the Alfvén speed, it takes about the time of a Keplerian revolution for the loop of the magnetic field to rise (e.g., Beloborodov 1999). This is also about the characteristic dissipation time of the magnetic field in shocks inside the disc (see section 2). The model we explore here differs from that of Field & Rogers (1993a, b) in that the dissipation of the magnetic energy occurs essentially inside the disc and the heat produced is transported to the disc surface and radiated away. Observations of hard X-ray flux indicate the presence of hot coronae where a significant fraction of the total accretion power is dissipated. For example, the X-ray band carries a significant fraction of the total luminosity of Seyfert nuclei: the flux in the 1–10 keV band is about $`1/6`$ of the total flux from IR to X-rays, and the flux in 1–500 keV band is about 30–40 per cent of the total energy output (Mushotzky, Done & Pounds, 1993). Another example is the low/hard state of galactic black hole sources, where the borad band spectrum is completely dominated by a hard X-ray power law, rolling over at energies of $`150\text{keV}`$ (Nowak, 1995; Done, 2002). Also, in the so called very high state, some of galactic black hole X-ray sources show both thermal and non-thermal (power law) components, with the ratio of non-thermal to total luminosity of 20–40 per cent (Nowak, 1995). Reconnection events and particle acceleration should also happen in rarefied strongly magnetized corona of the disc in our model and could cause observed X-ray flaring events. However, we do not consider the coronal dynamics here, and instead just focus on the structure and the emission spectrum of the disc itself.

Models of magnetized accretion discs with externally imposed large scale vertical magnetic field and anomalous magnetic field diffusion due to enhanced turbulent diffusion have also been considered (Shalybkov & Rüdiger, 2000; Campbell, 2000; Ogilvie & Livio, 2001). The magnetic field in these models was strong enough to be dynamically important. But those models are limited to the subsonic turbulence in the disc and the viscosity and magnetic diffusivity are due to hydrodynamic turbulence. Angular momentum transport in those models are due to the large scale global magnetic fields. Both small scale and large scale magnetic fields should be present in real accretion discs. Here we consider the possibility that the magnetic field has dominant small scale component, that is magnetic field inside the disc consists mostly of loops with size less than or comparable to the thickness of the disc.

We consider vertically integrated equations describing the radial structure of the magnetically dominated turbulent accretion disc and provide the solutions for the radial dependences of the averaged quantities in section 2. In section 3 we analyse the conditions for a magnetically supported disc to be self-consistent. In section 4 we calculate thermal emission spectra of magnetically supported disc taking into account scattering by free electrons.

## 2 Radial Disc Structure

Here we neglect effects of general relativity and do not consider the behaviour of the material closer than the radius of the innermost circular stable orbit $`r_\mathrm{s}`$. We assume a non-rotating black hole with $`r_\mathrm{s}=3r_\mathrm{g}`$, where the gravitational radius of the black hole of mass $`M=M_8\times 10^8\mathrm{M}_{\mathrm{}}`$ is $`r_\mathrm{g}=2GM/c^2=3\times 10^{13}\text{cm}M_8`$. We assume that accretion occurs in the form of a geometrically thin accretion disc and verify this assumption in section 3. We consider a disc of half-thickness $`H`$, surface density $`\mathrm{\Sigma }`$, averaged over the disc thickness volume density $`\rho =\mathrm{\Sigma }/(2H)`$, accretion rate $`\dot{M}`$, and radial inflow velocity $`v_r`$, $`v_r>0`$ for the accretion. We take $`\mathrm{\Omega }_\mathrm{K}=(GM/r^3)^{1/2}=7.2\times 10^4\text{s}^1M_8^1(r_\mathrm{g}/r)^{3/2}`$ to be the angular Keplerian frequency at the radial distance $`r`$ from the black hole. Then, equation of mass conservation reads

$$\dot{M}=2\pi r\mathrm{\Sigma }v_r\text{.}$$
(1)
In the stationary state $`\dot{M}`$ does not depend on time. Equation of angular momentum conservation is (e.g., Shapiro & Teukolsky (1983))

$$f_\varphi (r)2H2\pi r^2=\dot{M}\left(\sqrt{GMr}\zeta \sqrt{GMr_\mathrm{s}}\right)\text{,}$$
(2)
where $`f_\varphi `$ is the tangential stress at a radius $`r`$, which acts on the inner part of the disc, factor $`\zeta `$ accounts for the unknown torque acting on the disc at the inner edge $`r_\mathrm{s}`$. In the standard Shakura–Sunyaev model, $`\zeta =1`$, which corresponds to a zero torque inner boundary. The tangential stress at the innermost stable orbit, $`f_\varphi (r_\mathrm{s})`$ is related to the parameter $`\zeta `$ as follows

$$f_\varphi (r_\mathrm{s})2H2\pi r_\mathrm{s}^2=(1\zeta )\dot{M}\sqrt{GMr_\mathrm{s}}\text{.}$$
(3)
The stress $`f_\varphi (r_\mathrm{s})>0`$, if the torque spins up the gas at $`r_\mathrm{s}`$, and $`f_\varphi (r_\mathrm{s})<0`$, if the torque retards the rotation of the gas. Since the magnetic and turbulent pressures exceed the thermal pressure, by assumption, we neglect thermal pressure contribution to the vertical pressure balance. Assuming equipartition between turbulent and magnetic pressure in supersonic MHD turbulence (Stone, Ostriker & Gammie, 1998; Stone, 1999; Miller & Stone, 2000) it is easy to show that the equation of vertical equilibrium is solved to give approximately the result

$$H=\frac{v_\mathrm{A}}{\mathrm{\Omega }_\mathrm{K}}=\frac{1}{\mathrm{\Omega }_\mathrm{K}}\sqrt{\frac{HB^2}{2\pi \mathrm{\Sigma }}}\text{,}$$
(4)
where $`v_\mathrm{A}={\displaystyle \frac{B}{(4\pi \rho )^{1/2}}}`$ is the average Alfvén velocity in the disc. We use the ’$`\alpha `$–prescription’ for the magnetic viscosity in the disc, i.e. taking the tangential stress to be proportional to the sum of the magnetic pressure, $`{\displaystyle \frac{B^2}{8\pi }}`$, and the turbulent pressure, $`\rho v_\mathrm{t}^2/2`$. We do not consider global transport of angular momentum and energy across the disc due from large scale magnetic fields. We assume the magnetic field to be sufficiently tangled on scales of order of the thickness of the disc, such that large scale field is small compared to this tangled field produced by turbulence. In equipartition, $`{\displaystyle \frac{B^2}{8\pi }}\rho v_\mathrm{t}^2/2`$ and $`\alpha `$–prescription becomes

$$f_\varphi =\alpha \frac{B^2}{4\pi }\text{.}$$
(5)
The dissipation of the magnetic and kinetic energy causes heat input in the disc which is balanced by the heat losses due to radiation. If the cooling is efficient enough such that the time of radial advection of the heat due to the accretion flow is much longer than a few Keplerian periods, the heating and cooling balance, and establish an equilibrium kinetic temperature in the disc. However, we assume that this temperature is insufficient for the associated thermal pressure to contribute significant vertical support. The total radiated energy from the unit surface of the disc will be the same as in the standard Shakura–Sunyaev model. This energy flux is independent of the viscosity mechanism assumed, but depends on the inner torque boundary condition (3) (see also Gammie (2000) for the effects of varying the torque at $`r_\mathrm{s}`$). Thus, the system of equations (1), (2), (4) and (5) decouples from the energy balance equation.

Equation (4) can be solved to give $`H`$ as

$$H=\frac{B^2}{2\pi \mathrm{\Sigma }\mathrm{\Omega }_\mathrm{K}^2}\text{.}$$
(6)
With the prescription for the viscous stress, equation (5), the angular momentum conservation equation (2) becomes

$$\alpha \frac{B^4}{2\pi \mathrm{\Sigma }}=\dot{M}\left(1\zeta \sqrt{\frac{r_\mathrm{s}}{r}}\right)\mathrm{\Omega }_\mathrm{K}^3\text{.}$$
(7)
It is remarkable that the value of $`\alpha `$ can be constrained in our model. In the framework of the local approach all the work done by non-gravitational forces on a patch of the disc is reduced to the work done by viscous stress $`f_\varphi `$. Physically, this work results from the action of turbulent and Maxwell stresses. Heating occurs from dissipation of supersonic turbulence. The rate of such heating can be expressed through the kinematic viscosity coefficient $`\nu `$ in the usual way (heating per unit volume)

$$\rho T\frac{ds}{dt}=\frac{f_\varphi ^2}{\nu \rho }\text{,}$$
(8)
where $`T`$ is the average temperature of gas inside the disc, $`s`$ is the entropy per unit mass (e.g., see in Shapiro & Teukolsky (1983)). Viscous stress, $`f_\varphi `$, for Keplerian shear is

$$f_\varphi =\frac{3}{2}\nu \rho \mathrm{\Omega }_\mathrm{K}\text{.}$$
(9)
Comparing this to the alternative expression (5) for $`f_\varphi `$ we see that $`\alpha `$ and $`\nu `$ are related by

$$\nu =\alpha \frac{2}{3\mathrm{\Omega }_\mathrm{K}}v_\mathrm{A}^2\text{.}$$
(10)
Using expressions (9) and (10) to substitute in equation (8) one can obtain the rate of heating per unit area of the disc $`Q`$ as follows

$$Q=2H\rho T\frac{ds}{dt}=3\alpha \rho v_\mathrm{A}^3\text{.}$$
(11)
One can also rewrite $`Q`$ as

$$Q=\frac{3}{4\pi }\frac{GM\dot{M}}{r^3}\left(1\zeta \sqrt{\frac{r_\mathrm{s}}{r}}\right)$$
(12)
the same expression familiar in the standard disc model (e.g., see Shapiro & Teukolsky (1983)). The total energy density of the turbulent pulsations and magnetic fields is $`{\displaystyle \frac{B^2}{8\pi }}+{\displaystyle \frac{1}{2}}\rho v_\mathrm{t}^2=\rho v_\mathrm{A}^2`$ under the assumption of equipartition. In the stationary state, this turbulent energy dissipates with the rate $`Q`$. Therefore, the characteristic time of the dissipation of the turbulence is

$$\tau _\mathrm{A}=\frac{2H\rho v_\mathrm{A}^2}{Q}=\frac{2}{3\mathrm{\Omega }_\mathrm{K}\alpha }\text{.}$$
(13)
On the other hand, dissipation of the supersonic turbulence occurs on the time-scale $`\tau _\mathrm{t}`$ of the flow crossing the largest flow coherence size $`l_\mathrm{t}`$ of the turbulence. Direct dissipation at the shock fronts dominates the turbulent cascade of energy down to the microscopic resistive scale and leads to the enhanced rate of the dissipation. The question of the dissipation of supersonic MHD turbulence has been studied in connection with the interstellar turbulence, which is observed to be highly supersonic. Direct numerical simulations of both steady-state driven and freely decaying MHD turbulence (Stone et al. 1998; Stone 1999; Ostriker, Stone & Gammie 2001) all confirm this picture. Even if initially the motion is set up to be incompressible in the numerical experiments, shocks develop rapidly due to the non-linear conversion of the waves and the dissipation becomes dominated by the dissipation on shocks. This dissipation time is $`\tau _\mathrm{t}=l_\mathrm{t}/v_\mathrm{t}`$. The coefficient of turbulent viscosity $`\nu `$ is also related to the largest flow coherence size and turbulent velocity as

$$\nu =\frac{1}{3}l_\mathrm{t}v_\mathrm{t}\text{.}$$

We already know that $`v_\mathrm{t}=v_\mathrm{A}`$ from equipartition of magnetic and kinetic energies in turbulence. Also, $`\nu `$ is expressed through $`\alpha `$ by relation (10). This allows to estimate the largest flow coherence size of the turbulence

$$l_\mathrm{t}=2\alpha H\text{.}$$
(14)
Consequently, the flow crossing time of the largest coherence size becomes $`\tau _\mathrm{t}={\displaystyle \frac{2\alpha }{\mathrm{\Omega }_\mathrm{K}}}`$. It should be that $`\tau _\mathrm{t}=\tau _\mathrm{A}`$ with $`\tau _\mathrm{A}`$ given by expression (13). This is possible only when $`\alpha =1/\sqrt{3}`$. Due to approximate nature of all calculations which lead us to this value for $`\alpha `$, it is only meaningful to state that $`\alpha `$ should be of order of $`1`$. We assume $`\alpha =1`$ for all further estimates. For $`\alpha 1`$ the largest flow coherence size of the turbulence becomes of order of the thickness of the disc $`l_\mathrm{t}2H`$, and the turbulent viscosity coefficient takes its largest possible value $`\nu Hv_\mathrm{A}`$, still compatible with the local viscous description of the disc. The dissipation time-scale for the magnetic turbulence is $`\tau _\mathrm{A}1/\mathrm{\Omega }_\mathrm{K}`$. It is very probable that such large scale turbulence will lead to the buoyant rising of the magnetic field loops into the corona, subsequent shearing by the differential rotation and reconnection of the loops. This can result in the formation of the hot corona or acceleration of particles to relativistic energies (Field & Rogers, 1993a, b). The formation of a magnetized corona and the emission spectrum from the corona are important, however, here we focus on the disc.

The free parameters are $`\dot{M}`$ and $`\zeta `$. Also, we need to specify one more physical condition, since the dependence of $`B`$ on radius in equation (7) is undetermined. Such a condition should come from the modelling of supersonic turbulence. Lacking a detailed model, we assume that the radial dependence of the vertically averaged magnetic field in the disc is the power law

$$B=B_{10}\left(\frac{r}{10r_\mathrm{g}}\right)^\delta \text{,}$$
(15)
where $`B_{10}`$ is the strength of the magnetic field at $`10r_\mathrm{g}`$, $`\delta >0`$ is some constant. Accretion rate $`\dot{M}`$ can be related to the total luminosity of the disc $`L`$ and the radiated fraction $`ϵ`$ of the rest mass accretion flux $`\dot{M}c^2`$. At the luminosity of an AGN

$$L=1.3\times 10^{46}l_\mathrm{E}M_8\text{erg s}^1\text{,}$$
(16)
the mass flux is

$`\dot{M}=0.23\mathrm{M}_{\mathrm{}}\text{yr}^1\left({\displaystyle \frac{l_\mathrm{E}}{ϵ}}\right)M_8=`$
$`1.4\times 10^{25}\text{g s}^1\left({\displaystyle \frac{l_\mathrm{E}}{ϵ}}\right)M_8\text{.}`$ (17)
Here we denote the ratio of disc luminosity to the Eddington luminosity by $`l_\mathrm{E}=L/L_{\mathrm{edd}}`$. The value of $`ϵ`$ is determined by the inner boundary condition at $`r=r_\mathrm{s}`$. Typically, $`ϵ0.1`$. Using such parametrization and the expression (15) for the magnetic field, one can obtain the following radial profiles of $`H`$, $`\mathrm{\Sigma }`$ and $`\rho `$ from equations (6) and (7)

$`{\displaystyle \frac{H}{r_\mathrm{g}}}={\displaystyle \frac{\mathrm{\Omega }_\mathrm{K}\dot{M}𝒢}{\alpha B^2r_\mathrm{g}}}=2.1\times 10^1{\displaystyle \frac{l_\mathrm{E}}{2ϵ}}\left({\displaystyle \frac{B_{10}}{10^4\text{G}}}\right)^2\times `$
$`M_8^1𝒢\left({\displaystyle \frac{r}{10r_\mathrm{g}}}\right)^{2\delta 3/2}\text{,}`$ (18)
$`\mathrm{\Sigma }={\displaystyle \frac{\alpha B^4}{2\pi \mathrm{\Omega }_\mathrm{K}^3\dot{M}𝒢}}=5.1\times 10^3\text{g cm}^2\left({\displaystyle \frac{l_\mathrm{E}}{2ϵ}}\right)^1\times `$
$`\left({\displaystyle \frac{B_{10}}{10^4\text{G}}}\right)^4M_8^2𝒢^1\left({\displaystyle \frac{r}{10r_\mathrm{g}}}\right)^{9/24\delta }\text{,}`$ (19)
$`\rho ={\displaystyle \frac{\alpha ^2B^6}{4\pi \mathrm{\Omega }_\mathrm{K}^4\dot{M}^2𝒢^2}}=4\times 10^{10}\text{g cm}^3\left({\displaystyle \frac{l_\mathrm{E}}{2ϵ}}\right)^2\times `$
$`\left({\displaystyle \frac{B_{10}}{10^4\text{G}}}\right)^6𝒢^2M_8^2\left({\displaystyle \frac{r}{10r_\mathrm{g}}}\right)^{66\delta }\text{,}`$ (20)
where $`𝒢=1\zeta \sqrt{{\displaystyle \frac{r_\mathrm{s}}{r}}}`$. We see that the disc becomes geometrically thicker and less dense when magnetic field decreases: weaker magnetic field leads to weaker tangential stress (equation (5)); $`H`$ increases in order to accomodate constant angular momentum flux for the same $`\dot{M}`$ such that $`HB_{10}^2`$ (equation (2)); $`\mathrm{\Sigma }`$ and $`\rho `$ should decrease strongly, $`\mathrm{\Sigma }B_{10}^4`$ and $`\rho B_{10}^6`$, in order to ensure vertical equilibrium with larger $`H`$ and less pressure support from $`B^2`$ (equation (6)); radial inflow velocity increases as $`v_rB_{10}^4`$ to ensure constant mass flux.

Let us now summarize the similarities and differences between our model and the standard Shakura-Sunyaev model. If we replace the thermal pressure in the standard model by the sum of the magnetic and turbulent pressures, the equations for mass conservation (1), angular momentum conservation (2), the viscosity prescription (5) and vertical pressure support (4) are the same as in the standard model. The pressure in the standard model is determined by the rate of the cooling of the disc, while the $`\alpha `$ coefficient can be an arbitrary function of $`r`$, $`\alpha (r)<1`$. In our model we have the magnetic pressure unspecified in its radial dependence as soon as it exceeds the thermal pressure, but $`\alpha 1`$ for all $`r`$. The latter results from the much faster dissipation of supersonic turbulence than subsonic turbulence assumed in the standard model. Both our model and the standard model have only one undetermined function of radius, ($`\alpha (r)`$ in the standard model and $`B(r)`$ in our model). The determination of this free function would eventually come from detailed modelling of the MHD turbulence.

## 3 Estimates of the Disc Parameters

Now let us use the solution for the disc structure provided by equations (6), (7) and (15) and obtain constraints on free parameters of the model, such that our model of thin magnetized accretion disc is self-consistent. Using equations (7) and (6) to substitute for $`\dot{M}`$ in equation (1) we can express the radial inflow velocity as

$$\frac{v_r}{v_\mathrm{K}}=\frac{\alpha }{𝒢}\frac{H^2}{r^2}\text{.}$$
(21)
The factor $`𝒢`$ vanishes at $`r=r_\mathrm{s}`$ for $`\zeta =1`$ and so the surface density $`\mathrm{\Sigma }`$ of the disc diverges near $`r=r_\mathrm{s}`$. The same divergence also occurs in standard disc model (Shakura & Sunyaev, 1973). In reality, of course, viscous torque does not vanish at $`r=r_\mathrm{s}`$, $`\zeta 1`$ but close to $`1`$ and there is no divergence of $`\mathrm{\Sigma }`$ at $`r=r_\mathrm{s}`$. Determining $`\zeta `$ would require a full general relativistic treatment of the accretion flow close to the black hole and is beyond the immediate scope of this work. Only parts of the disc close to $`r_\mathrm{s}`$ are sensitive to the exact value of $`\zeta `$. For $`r>r_\mathrm{s}`$, the disc structure is approximated well with the formulae in section 2 for $`\zeta =1`$. Provided that $`r>r_\mathrm{s}`$ and $`\alpha 1`$ one can see that the radial inflow velocity for a thin disc ($`Hr`$) is always a small fraction of the Alfvén velocity which, in turn, is a small fraction of the Keplerian velocity. Therefore, the dissipation time-scale $`\tau _\mathrm{A}`$ is always much shorter than the radial inflow time-scale. In the stationary case this means that the advective terms in the energy balance equation can be always neglected. Energy balance becomes local: the rate of gas heating $`Q`$ should be balanced by the cooling due to radiative processes. Now we consider the physics of radiative cooling which determines the disc temperature.

A necessary condition for the existence of a magnetically dominated disc is that the vertical escape time for radiation must be shorter than $`\tau _\mathrm{A}`$, so that the energy density of radiation, $`\epsilon _\mathrm{r}=aT^4`$, remains smaller than the energy density of the magnetic field $`\epsilon _\mathrm{B}={\displaystyle \frac{B^2}{8\pi }}`$. For optically thin, geometrically thin discs this condition is always satisfied since the inverse of escape time $`c/H\mathrm{\Omega }_\mathrm{K}`$. As we will see, Thomson scattering is the dominant source of opacity in most cases of optically thick discs. The average time it takes for a photon to escape out of optically thick disc with optical depth $`\tau 1`$ is $`\tau H/c`$. For Thomson scattering $`\tau =\tau _\mathrm{c}=Hn_\mathrm{e}\sigma _\mathrm{T}=\kappa _\mathrm{T}\mathrm{\Sigma }/2`$, where $`n_\mathrm{e}`$ is the number density of free electrons, $`\sigma _\mathrm{T}`$ is the Thomson cross section and $`\kappa _\mathrm{T}=0.4\text{cm}^2\text{g}^1`$ is the Thomson scattering opacity. For simplicity we assume the composition of the disc to be completely ionized hydrogen. Then, $`n_\mathrm{e}`$ is equal to the number density of protons in the disc, $`n`$. The necessary condition now becomes

$$\frac{H^2}{c}n\sigma _\mathrm{T}<\frac{1}{\mathrm{\Omega }_\mathrm{K}}\text{.}$$
(22)
Using equations (6) and (7) one can rewrite the condition (22) as

$$\frac{l_\mathrm{E}}{2ϵ}\frac{r_\mathrm{g}}{H}𝒢^2<1\text{.}$$
(23)
We express this and all the subsequent conditions in terms of free parameters of the model: $`B_{10}`$, $`\delta `$, $`l_\mathrm{E}/ϵ`$, and $`M_8`$. Using the expression (18) for $`H`$, the necessary condition (23) becomes

$$4.7\times \left(\frac{B_{10}}{10^4\text{G}}\right)^2M_8𝒢\left(\frac{r}{10r_\mathrm{g}}\right)^{3/22\delta }1\text{.}$$
(24)
The condition for the disc to be optically thick for Thomson scattering $`\tau _\mathrm{c}1`$ becomes

$`2.0\times 10^3\left({\displaystyle \frac{l_\mathrm{E}}{2ϵ}}\right)^1\left({\displaystyle \frac{B_{10}}{10^4\text{G}}}\right)^4\times `$
$`M_8^2𝒢^1\left({\displaystyle \frac{r}{10r_\mathrm{g}}}\right)^{9/24\delta }1\text{.}`$ (25)
In an optically thick disc, radiation is transported by turbulent motions and radiative diffusion. The characteristic time-scale for the turbulent transport of radiation to the surface of the disc is $`1/\mathrm{\Omega }_\mathrm{K}`$. When condition (22) is satisfied, the diffusion of radiation dominates the advection due to turbulent motions. Thus, we can neglect turbulent convective transport of radiation for any values of parameters, whenever steady state magnetically dominated model of the accretion disc is considered. In the steady state, the radiation flux from each surface of the disc must be equal to $`Q/2`$ with the dissipation rate $`Q`$ given by expressions (11) or (12). The effective temperature of the escaping radiation flux is determined by this dissipation rate $`Q`$ as

$`T_{\mathrm{eff}}=\left({\displaystyle \frac{2Q}{ac}}\right)^{1/4}=7.5\times 10^4\text{K}\left({\displaystyle \frac{l_\mathrm{E}}{2ϵ}}\right)^{1/4}\times `$
$`\left({\displaystyle \frac{r}{10r_\mathrm{g}}}\right)^{3/4}𝒢^{1/4}M_8^{1/4}\text{,}`$ (26)
and is the same in our model as in the standard disc model. True absorption of photons in free-free transitions also occur in the disc. For an approximate estimate of the radiative conditions in the disc we consider Rosseland mean of the free-free absorption opacity

$$\overline{\kappa }_{\mathrm{ff}}=6.4\times 10^{22}\text{cm}^2\text{g}^1\rho T^{7/2}\text{,}$$
(27)
where $`\rho `$ is expressed in $`\text{g cm}^3`$ and $`T`$ is expressed in $`{}_{}{}^{\mathrm{°}}\text{K}`$. When the effective optical thickness of the disc $`2\overline{\tau }_{}=\mathrm{\Sigma }\sqrt{(\kappa _\mathrm{T}+\overline{\kappa }_{\mathrm{ff}})\overline{\kappa }_{\mathrm{ff}}}`$ is large, then local thermal equilibrium is established in the disc and the radiative flux is described by diffusion approximation. Only in the thin surface layer at a distance from the surface less then the thermalization optical depth $`\overline{\tau }_{}1`$, the spectrum of radiation deviates substantially from that of the black body (see section 4 below). The solution of the diffusive radiation transport gives the usual result relating the temperature at the midplane of the disc, $`T_{\mathrm{mpd}}`$, with the effective temperature at the surface (Shakura & Sunyaev, 1973; Krolik, 1999)

$`T_{\mathrm{mpd}}T_{\mathrm{eff}}\left({\displaystyle \frac{3\kappa _\mathrm{T}\mathrm{\Sigma }}{16}}\right)^{1/4}=3.3\times 10^5\text{K}\times `$
$`\left({\displaystyle \frac{B_{10}}{10^4\text{G}}}\right)M_8^{1/4}\left({\displaystyle \frac{r}{10r_\mathrm{g}}}\right)^{3/8\delta }\text{.}`$ (28)
The averaged temperature in the disc is close to $`T_{\mathrm{mpd}}`$ with the actual profile being determined by the details of the vertical dependence of the dissipation rate in shocks. As in the case of the inner radiation dominated part of the standard disc, $`T_{\mathrm{mpd}}`$ does not depend on the accretion rate. However, it is directly proportional to the value of the magnetic field and has radial dependence governed by the $`\delta `$.

The dominance of the magnetic and turbulent energy compared to the energy density of radiation is expressed as $`aT_{\mathrm{mpd}}^4\rho v_\mathrm{A}^2`$. One can substitute here for $`T_{\mathrm{mpd}}`$ from equations (28) and (26). Heating rate $`Q`$ is given by equation (11). After using equations (1820) to manipulate with $`\mathrm{\Sigma }`$, $`B`$ and $`H`$, one can reduce the condition of magnetic pressure dominance over the radiation pressure to be

$$\frac{9\alpha }{4c}\sigma _\mathrm{T}nH^2\frac{1}{\mathrm{\Omega }_\mathrm{K}}\text{.}$$
(29)
This differs only by a factor of order unity from the necessary condition for the escape of radiation, equation (22). The assertion that the condition (24) implies smallness of radiation pressure is thus confirmed by direct calculation. The next condition we need to impose is that the gas pressure is small compared to the magnetic pressure and turbulent stresses, which is equivalent to the statement that the turbulence is highly supersonic. For a thermalized radiation field, this condition is $`nkT_{\mathrm{mpd}}\rho v_\mathrm{A}^2`$. Using expressions (28) and (26) for the temperature and substituting for the density and $`H`$ from expressions (1820) one can express the condition $`nkT_{\mathrm{mpd}}\rho v_\mathrm{A}^2`$ as

$$\left(\frac{9k^4\sigma _\mathrm{T}}{32\pi ^2acm_\mathrm{p}^5}\right)^{1/2}\alpha ^{9/2}\mathrm{\Omega }_\mathrm{K}^{17/2}B^{10}\dot{M}^4𝒢^41\text{.}$$
(30)
Substituting for $`\mathrm{\Omega }_\mathrm{K}`$, $`B`$, and $`\dot{M}`$ in (30) and taking the $`1/5`$ power we obtain

$`13\times M_8^{9/10}\left({\displaystyle \frac{r}{10r_\mathrm{g}}}\right)^{51/20+2\delta }𝒢^{4/5}\times `$
$`\left({\displaystyle \frac{l_\mathrm{E}}{2ϵ}}\right)^{4/5}\left({\displaystyle \frac{B_{10}}{10^4\text{G}}}\right)^21\text{.}`$ (31)
Because the free-free opacity $`\overline{\kappa }_{\mathrm{ff}}`$ strongly depends on temperature and density (equation (27)), the vertical density and temperature distributions are needed to evaluate $`\overline{\tau }_{}`$. Simulations in Miller & Stone (2000) show a sharp density drop off by two orders of magnitude at the surfaces of the disc (see figure 11 in Miller & Stone (2000)). Inside the slab bounded by this density drop off the magnetic field and kinetic energy are in approximate equipartition. For the purposes of calculating effective optical depth for absorption we approximate the density profile between $`z=H`$ and $`z=H`$ inside the disc as constant and assume zero density at the disc surfaces at $`z=\pm H`$. Modest variations of magnetic field and density across the disc height observed in numerical simulations support this and also suggest the assumption of a uniform turbulence dissipation rate across the thickness. With these approximations, the solution for diffusive radiation transport in the vertical direction $`z`$ can be written (Shakura & Sunyaev, 1973; Krolik, 1999)

$$T^4=T_{\mathrm{mpd}}^4\left[\frac{8}{3\tau _\mathrm{c}}+\left(1\frac{z}{H}\right)\right]\text{.}$$
(32)
The Eddington approximation near the surface $`z=H`$ of the disc was used to obtain this temperature distribution, and $`T_{\mathrm{eff}}=\sqrt[4]{2}T|_{z=H}`$ according to the solution (32). This solution is not valid near the surface of the disc in the region dominated by Compton scattering but gives $`T`$ in the bulk of the disc if $`\overline{\tau }_{}1`$. Using the expression (32) for $`T`$ one has for the optical depth of the disc for free-free absorption

$`\overline{\tau }_{\mathrm{ff}}={\displaystyle _0^H}\rho \kappa _{\mathrm{ff}}dz=6.4\times 10^{22}\rho ^2T_{\mathrm{mpd}}^{7/2}H\times `$
$`{\displaystyle _0^1}\left(1\xi +{\displaystyle \frac{8}{3\tau _\mathrm{c}}}\right)^{7/8}𝑑\xi \text{,}`$ (33)
where $`\xi =z/H`$, $`\rho `$ and $`T`$ are expressed in $`\text{g cm}^3`$ and $`{}_{}{}^{\mathrm{°}}\text{K}`$. The value of the integral over $`\xi `$ determines how much the actual value $`\overline{\tau }_{\mathrm{ff}}`$ is larger than the value of $`\overline{\tau }_{\mathrm{ff}}`$ obtained if one assumes $`T=T_{\mathrm{mpd}}`$ everywhere in the disc. The $`\xi `$ integral in equation (33) is calculated to be

$`{\displaystyle _0^1}\left(1\xi +{\displaystyle \frac{8}{3\tau _\mathrm{c}}}\right)^{7/8}𝑑\xi =`$
$`8\left[\left(1+{\displaystyle \frac{8}{3\tau _\mathrm{c}}}\right)^{1/8}\left({\displaystyle \frac{8}{3\tau _\mathrm{c}}}\right)^{1/8}\right]\text{.}`$
Because $`\tau _\mathrm{c}1`$ the number in square parenthesis is close to $`1`$, so we can omit it and obtain the final result for $`\overline{\tau }_{\mathrm{ff}}`$

$$\overline{\tau }_{\mathrm{ff}}=5\times 10^{23}\rho ^2T_{\mathrm{mpd}}^{7/2}H\text{.}$$
(34)
When one substitutes for $`T_{\mathrm{mpd}}`$, $`H`$, and $`\rho `$ their expressions through the magnetic field and the accretion rate, the expression (34) becomes

$`\overline{\tau }_{\mathrm{ff}}=8.6\times 10^2M_8^{25/8}\left({\displaystyle \frac{r}{10r_\mathrm{g}}}\right)^{147/1613\delta /2}𝒢^3\times `$
$`\left({\displaystyle \frac{l_\mathrm{E}}{2ϵ}}\right)^3\left({\displaystyle \frac{B_{10}}{10^4\text{G}}}\right)^{13/2}\text{.}`$ (35)
Generally, for $`\delta 1`$, $`M_81`$, $`l_\mathrm{E}/ϵ0.1`$, $`B_{10}10^4\text{G}`$, and $`r10r_\mathrm{g}`$, $`\overline{\tau }_{\mathrm{ff}}1`$. However, because of steep dependence of $`\overline{\tau }_{\mathrm{ff}}`$ on $`l_\mathrm{E}`$, $`M_8`$, and $`B_{10}`$, the value of $`\overline{\tau }_{\mathrm{ff}}`$ can become large for lower accretion rates, more massive black holes, and stronger magnetic fields. The ratio of $`\overline{\tau }_{\mathrm{ff}}`$ to $`\tau _\mathrm{c}`$ is

$`{\displaystyle \frac{\overline{\tau }_{\mathrm{ff}}}{\tau _\mathrm{c}}}=4.2\times 10^5M_8^{9/8}\left({\displaystyle \frac{r}{10r_\mathrm{g}}}\right)^{75/165\delta /2}𝒢^2\times `$
$`\left({\displaystyle \frac{l_\mathrm{E}}{2ϵ}}\right)^2\left({\displaystyle \frac{B_{10}}{10^4\text{G}}}\right)^{5/2}\text{.}`$ (36)
We see that for the typical values of the parameters above the ratio $`\overline{\tau }_{\mathrm{ff}}/\tau _\mathrm{c}10^2`$ but free-free optical depth can exceed Compton scattering optical depth for smaller values of $`l_\mathrm{E}`$ and at larger radii. When $`\overline{\tau }_{\mathrm{ff}}\tau _\mathrm{c}`$ the effective optical thickness of the disc becomes $`\overline{\tau }_{}=\sqrt{\overline{\tau }_{\mathrm{ff}}\tau _\mathrm{c}}`$ or

$`\overline{\tau }_{}=13\times M_8^{41/16}\left({\displaystyle \frac{r}{10r_\mathrm{g}}}\right)^{219/3221\delta /4}𝒢^2\times `$
$`\left({\displaystyle \frac{l_\mathrm{E}}{2ϵ}}\right)^2\left({\displaystyle \frac{B_{10}}{10^4\text{G}}}\right)^{21/4}\text{.}`$ (37)
When $`\overline{\tau }_{\mathrm{ff}}\tau _\mathrm{c}`$ the effective optical thickness is equal to $`\overline{\tau }_{\mathrm{ff}}`$.

Finally, the ratio of the disc semi-thickness $`H`$ to the radius $`r`$ using (18) is

$$\frac{H}{r}=2.1\times 10^2\frac{l_\mathrm{E}}{2ϵ}\left(\frac{B_{10}}{10^4\text{G}}\right)^2M_8^1𝒢\left(\frac{r}{10r_\mathrm{g}}\right)^{2\delta 5/2}\text{.}$$
(38)
Now we summarize conditions when our model is valid:

1. The necessary condition (24), which is also the condition for the dominance of the magnetic pressure over the radiation pressure.
2. The condition (31) for the dominance of the magnetic pressure over the thermal pressure of gas in thermalized optically thick disc.
3. The condition for the disc to be optically thick. This is either condition (25) $`\tau _\mathrm{c}1`$ or the condition that $`\overline{\tau }_{\mathrm{ff}}`$ given by expression (35) is greater than $`1`$.
4. The condition that the gas and radiation inside the disc are in local thermal equilibrium, $`\overline{\tau }_{}1`$, where $`\overline{\tau }_{}=\overline{\tau }_{\mathrm{ff}}`$ if $`\overline{\tau }_{\mathrm{ff}}/\tau _\mathrm{c}>1`$, and $`\overline{\tau }_{}`$ is given by expression (37) if $`\overline{\tau }_{\mathrm{ff}}/\tau _\mathrm{c}<1`$.
5. $`H/r1`$.

We varied the parameters $`l_\mathrm{E}/ϵ`$, $`\delta `$, $`M_8`$ to obtain the allowed region of our disc model in the $`B_{10}`$, $`r/r_\mathrm{g}`$ plane, using the above five conditions. These plots are shown in Figs. 1 – 2.

Depending on the power $`\delta `$ in the dependence of magnetic field $`Br^\delta `$, optically thick magnetically dominated accretion discs can exist only at a limited interval of radii. For $`\delta 1`$ (as on our plots for $`\delta =5/4`$) a thin magnetically dominated disc (shadowed regions on the plots) is possible for $`5r_\mathrm{g}<r<1000r_\mathrm{g}`$. The window for the strength of the magnetic field is not very wide: about one order of magnitude or less. This window is narrower for low masses of the central black hole and is wider for higher masses. The value $`l_\mathrm{E}/(2\alpha ϵ)=10`$ corresponds to about the Eddington accretion rate for $`ϵ0.1`$, and because of $`\alpha 1`$ in our model (see section 2). Higher values of $`l_\mathrm{E}/(2\alpha ϵ)`$ correspond to higher accretion rates. Allowed values of the magnetic field are higher for higher accretion rates. The magnetic fields in the discs around higher mass black holes are smaller than in the discs around lower masses black holes as is temperature of the disc ($`T_{\mathrm{mpd}}`$) and the surface radiation flux. For large luminosities ($`l_\mathrm{E}/(2\alpha ϵ)1`$) the inner disc cannot be optically thick for true absorption but can be optically thick to free electron scattering. Comptonisation becomes significant for in the inner regions at high luminosities (see section 4). We leave consideration of Comptonised regimes for future work.

Physically, the limitations on the magnetic field strength can be understood as follows: suppose one decreases $`B_{10}`$ while keeping $`\dot{M}`$ constant. Then, $`HB_{10}^2`$ is increasing; $`\mathrm{\Sigma }B_{10}^4`$, and $`\rho B_{10}^6`$, both decreasing (equations (1820)). Scattering opacity through the disc $`\tau _\mathrm{c}\mathrm{\Sigma }`$ strongly decreases, so the heat is transported to the surface faster and $`T_{\mathrm{mpd}}B_{10}`$ decreases (equation (28)); thermal and radiation pressures decrease as $`P_{\mathrm{th}}B_{10}^7`$ and $`P_{\mathrm{rad}}B_{10}^4`$ respectively; plasma parameter $`\beta `$ decreases, so the disc becomes more magnetically dominated; $`\overline{\tau }_{\mathrm{ff}}`$ and $`\overline{\tau }_{}`$ both decrease since their decrease due to lower $`\mathrm{\Sigma }`$ overcomes the increase of the absorptive opacity from the drop of the temperature. Therefore, there exists only a limited interval of $`B_{10}`$ such that $`\beta <1`$ still the disc is optically thick to true absorption.

## 4 Radiation Spectra of Optically Thick Magnetically Dominated Disc

### 4.1 Calculation of spectra

Free-free, bound-free and cyclotron emission could contribute to the radiation spectrum. In Appendix A we show that, because of the self-absorption in the dense disc, the total flux of cyclotron emission from the disc surface is negligibly small compared to the total radiated power $`Q`$. This power is entirely due to free-free and bound-free radiative transitions. Cyclotron and synchrotron emission can be important in the rarefied and strongly magnetized disc corona (Ikhsanov 1989; Field & Rogers 1993a, Di Matteo et al. 1997a), but our focus here is on the disc.

We perform a simplified calculation of the emergent spectrum. We assume local thermodynamic equilibrium and do not consider effects of the temperature change with depth. This is justified when the spectrum is formed in thin layer near the disc surface. For simplicity we do not include the bound-free contribution to the opacity. Free-free opacities are the dominant source of thermal absorption for $`T10^5\text{K}`$, so our simplified spectrum is most relevant for smaller masses of the central black hole, for which the inner disc is hotter. Our goal here is to capture the effect of the magnetic field on the shape of the spectrum. We consider only optically thick disc models with both $`\tau _\mathrm{c}1`$ and $`\tau _{}1`$. The electron scattering opacity $`\kappa _\mathrm{T}`$ does not depend on frequency in the non-relativistic limit, whereas free-free absorption opacity $`\kappa _{\mathrm{ff}}(\nu )`$ is a function of frequency:

$$\kappa _{\mathrm{ff}}(\nu )=1.5\times 10^{25}\text{cm}^2\text{g}^1\rho T^{7/2}\frac{1e^x}{x^3}g(x,T)\text{,}$$
(39)
where we denote $`x={\displaystyle \frac{h\nu }{kT}}`$. Thus, the parameters of the disc model $`B_{10}`$ and $`\delta `$ will affect the spectrum by means of $`\kappa _{\mathrm{ff}}`$ dependence on $`\rho `$ and $`T/T_{\mathrm{eff}}`$. Gaunt factor $`g(x,T)`$ is slowly varying function of $`x`$ and $`T`$ between approximately $`0.3`$ and $`5`$ in wide range of frequencies and temperatures (Rybicki & Lightman, 1979). Moreover, in the wide temperature interval $`10^3\text{K}<T<10^8\text{K}`$ Gaunt factor is between $`0.5`$ and $`2`$ for frequencies $`0.1<x<10`$ near the maximum of thermal emission. It is quite reasonable to set $`g(x,T)=1`$, which we do in all further calculations. Note, that $`\kappa _{\mathrm{ff}}`$ behaves as $`1/x^2`$ for $`x1`$, so free-free absorption is always more important at lower frequencies than electron scattering.

The energy transfer due to repeated scatterings (Comptonisation process) is characterized by Compton $`y(\nu )`$ parameter (Rybicki & Lightman, 1979)

$$y=\frac{4kT}{m_\mathrm{e}c^2}\text{Max}(\tau _{\mathrm{es}},\tau _{\mathrm{es}}^2)\text{,}$$
(40)
where the optical depth for Thomson scattering $`\tau _{\mathrm{es}}(\nu )`$ must be measured from the effective absorption optical thickness $`\tau _{}(\nu )1`$ and is given by (formula \[7.42\] in Rybicki & Lightman (1979))

$$\tau _{\mathrm{es}}(\nu )\left(\frac{\kappa _\mathrm{T}/\kappa _{\mathrm{ff}}(\nu )}{1+\kappa _{\mathrm{ff}}(\nu )/\kappa _\mathrm{T}}\right)^{1/2}\text{.}$$
(41)
If $`y(\nu )1`$, photons do not change their initial frequency $`\nu `$ in the process of repeated scatterings before they escape the surface of the disc. While if $`y1`$ the Comptonisation effects become significant. We calculated values of $`y`$ using expressions (40), (41) and (39). One can see that $`y`$ is always monotonically rising with frequency, therefore, the Comptonisation effect at higher frequencies is always more significant than at lower frequencies. On the other hand, there is very little radiation at $`h\nu 5kT`$ due to exponential cut off in thermal spectra. It turns out that in most cases $`y(\nu )1`$ for optically thick disc models and for $`h\nu <5kT`$. Exceptions are the cases of high accretion rates $`l_\mathrm{E}/ϵ1`$. In those cases, inner parts of the accretion disc ($`r<30r_\mathrm{g}`$) can have $`y1`$ and Comptonisation will influence the highest frequencies of the disc spectrum. Ignoring these exceptions, we did not take into account Comptonisation in the following calculations and assume coherent scattering. This assumption has been checked a posteriori for self-consistency.

In the case of coherent scattering the approximate expression of the radiative flux per unit surface of an optically thick medium is given by (Rybicki & Lightman, 1979)

$$F_\nu \frac{4\pi h\nu ^3/c^2}{\left(e^{h\nu /kT}1\right)\left(1+\sqrt{1+\kappa _\mathrm{T}/\kappa _{\mathrm{ff}}(\nu ,T)}\right)}\text{.}$$
(42)
$`F_\nu `$ approaches black body spectrum $`\pi B_\nu (T)`$ in the limit $`\kappa _{\mathrm{ff}}\kappa _\mathrm{T}`$ and becomes modified black body spectrum

$$F_\nu =2\pi B_\nu (T)\sqrt{\kappa _{\mathrm{ff}}/\kappa _\mathrm{T}}$$
(43)
in the limit $`\kappa _{\mathrm{ff}}\kappa _\mathrm{T}`$. A part of the disc can emit black body spectrum at lower frequencies $`\nu \nu _0`$ and modified black body spectrum at higher frequencies $`\nu \nu _0`$, where $`\nu _0=\nu _0(\rho ,T)`$ is defined such that $`\kappa _{\mathrm{ff}}(\nu _0)=\kappa _\mathrm{T}`$. However, if frequency averaged free-free opacity is larger than Thomson opacity, $`\overline{\kappa }_{\mathrm{ff}}(T,\rho )>\kappa _\mathrm{T}`$, then $`h\nu _0>kT`$ and almost all the radiation is emitted as a black body spectrum. In the opposite limit, $`\overline{\kappa }_{\mathrm{ff}}(T,\rho )\kappa _\mathrm{T}`$, one has $`h\nu _0kT`$ and the spectrum is mostly modified black body, transiting to black body only at lowest frequencies, $`\nu <\nu _0`$. Dotted line on Figs. 12 separates regions with $`\overline{\kappa }_{\mathrm{ff}}>\kappa _\mathrm{T}`$ and $`\overline{\kappa }_{\mathrm{ff}}<\kappa _\mathrm{T}`$. We see that the optically thick disc has $`\overline{\kappa }_{\mathrm{ff}}\kappa _\mathrm{T}`$ and emits a modified black body spectrum in the inner parts but may become absorption dominated in the cooler outer parts, where black body spectrum will be emitted. At lower accretion rates, $`l_\mathrm{E}/(2ϵ)<10^4`$, all the surface of the disc will emit black body spectrum (with different $`T`$ at different radii, of course). The locations of the regions with mostly black body and modified black body spectra over the disc radii are fairly insensitive to the black hole mass $`M`$. The surface temperature $`T`$ in formula (42) should be determined by equating the total emitted flux, the integral of $`F_\nu `$ over all frequencies, to the half of the heat production rate $`Q`$ per unit disc surface (half is to account for two surfaces of the disc). Introducing variable $`x={\displaystyle \frac{h\nu }{kT}}`$ one can write this energy balance condition as

$$\frac{Q}{2}=4\pi \frac{k^4T^4}{h^3c^2}_0^{\mathrm{}}\frac{x^3dx}{(e^x1)\left(1+\sqrt{1+\kappa _\mathrm{T}/\kappa _{\mathrm{ff}}}\right)}\text{.}$$
(44)
Further introducing Stefan–Boltzmann constant $`\sigma ={\displaystyle \frac{2\pi ^5k^4}{15c^2h^3}}`$ into right hand side of equation (44) and effective temperature $`Q/2=\sigma T_{\mathrm{eff}}^4`$ (equation (26)) into left hand side of equation (44) we can transform equation (44) to

$$\frac{T^4}{T_{\mathrm{eff}}^4}=\frac{\pi ^4/15}{{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{2}{1+\sqrt{1+\kappa _\mathrm{T}/\kappa _{\mathrm{ff}}}}}{\displaystyle \frac{x^3dx}{e^x1}}}\text{.}$$
(45)
Since $`{\displaystyle \frac{\pi ^4}{15}}={\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{x^3dx}{e^x1}}`$, $`T`$ is always larger than $`T_{\mathrm{eff}}`$. This is in accordance with general thermodynamic argument that the black body is the most efficient emitter of all. Equation (45) together with expression (39) for $`\kappa _{\mathrm{ff}}`$, expression (20) for $`\rho `$, expression (26) for $`T_{\mathrm{eff}}`$, and $`\kappa _\mathrm{T}=0.4\text{cm}^2\text{g}^1`$ has been solved numerically to determine $`T(r)`$. After one knows $`T=T(r)`$, it is possible to integrate $`F_\nu (r)`$ (42) over the disc surface to obtain the spectral distribution of the total energy emitted by the disc

$$E_\nu =2_{r_{\mathrm{in}}}^{r_{\mathrm{out}}}2\pi rF_\nu 𝑑r\text{,}$$
(46)
where the $`2`$ accounts for the two surfaces of the disc.

When a significant interval of radii exists where the emitted spectrum is a modified black body, e.g. $`h\nu _0<kT`$, it is possible to get an approximate analytic expression for $`E_\nu `$. For $`\nu >\nu _0`$ we use expression (43) for $`F_\nu `$, which becomes

$$F_\nu =2.6\times 10^3\frac{\text{erg}}{\text{s}\text{cm}^2\text{Hz}}T^{5/4}\rho ^{1/2}x^{3/2}\frac{e^x}{\sqrt{1e^x}}\text{.}$$
(47)
Integrating expression (47) over frequencies by integrating from $`0<x<\mathrm{}`$, we obtain

$$\frac{Q}{2}=8.2\times 10^7\frac{\text{erg}}{\text{s}\text{cm}^2}T^{9/4}\rho ^{1/2}\text{,}$$

which can be solved for the temperature using expression (12) for $`Q`$. Let us denote the temperature found in this way by $`T_\mathrm{s}`$. The expression for $`T_\mathrm{s}`$ is

$`T_\mathrm{s}=2.1\times 10^5\text{K}\left({\displaystyle \frac{l_\mathrm{E}}{2ϵ}}\right)^{8/9}M_8^{8/9}\left({\displaystyle \frac{B_{10}}{10^4\text{G}}}\right)^{4/3}\times `$
$`\left({\displaystyle \frac{r}{10r_\mathrm{g}}}\right)^{8/3+4\delta /3}𝒢^{8/9}\text{.}`$ (48)
Now we use $`T_\mathrm{s}`$ to substitute in the expression (47) for $`F_\nu `$, and then to expression (46) to obtain $`E_\nu `$. It is convenient to switch from the integration in $`r`$ to the integration in $`x`$ in equation (46). We do so by writing $`dx={\displaystyle \frac{h\nu }{kT_\mathrm{s}^2}}{\displaystyle \frac{dT_\mathrm{s}}{dr}}dr`$ and expressing $`x`$ through $`r`$ for a given $`\nu `$ from equation (48). This procedure can be done analytically if one puts $`𝒢=1`$, that is our analytical approximation do not describe spectrum emerging close to the inner edge of the disc, where $`𝒢`$ deviates from $`1`$ significantly. Carrying out calculations one obtains

$`E_\nu =4.2\times 10^{44}\times 10^{\frac{11.73}{2\delta }}{\displaystyle \frac{\text{erg}}{\text{s}\text{Hz}}}\left({\displaystyle \frac{l_\mathrm{E}}{2ϵ}}\right)^{\frac{43\delta }{3(2\delta )}}\times `$
$`M_8^{\frac{83\delta }{3(2\delta )}}\left({\displaystyle \frac{B_{10}}{10^4\text{G}}}\right)^{\frac{1}{2\delta }}{\displaystyle \frac{1}{(2\delta )}}\nu ^{\frac{4\delta 5}{4(2\delta )}}\times `$ (49)
$`{\displaystyle _{x_{\mathrm{in}}}^{x_{\mathrm{out}}}}x^{\frac{3(2\delta 3)}{4(\delta 2)}}{\displaystyle \frac{e^x}{\sqrt{1e^x}}}𝑑x\text{,}`$
where $`x_{\mathrm{in}}={\displaystyle \frac{h\nu }{kT_\mathrm{s}(r_{\mathrm{in}})}}`$ and $`x_{\mathrm{out}}={\displaystyle \frac{h\nu }{kT_\mathrm{s}(r_{\mathrm{out}})}}`$. We take for typical estimates $`r_{\mathrm{in}}=10r_\mathrm{g}`$ and $`r_{\mathrm{out}}=1000r_\mathrm{g}`$ as the inner and outer edges of the disc. Although inner parts of the disc contribute significantly to the total emitted power and determine the most energetic part of the spectrum, the calculation of the spectrum there must be performed by taking into account factor $`𝒢`$, not to mention relativistic effects. The outer extension of the disc at $`1000r_\mathrm{g}`$ is somewhat arbitrary, but the disc beyond $`1000r_\mathrm{g}`$ is too cool to be described by our simple radiation model and, in the case of AGNs, even the existence of the disc for $`r1000r_\mathrm{g}`$ is questionable because of the instability to the gravitational fragmentation. The value of the integral in formula (49) decreases exponentially for $`x_{\mathrm{in}}>1`$. This corresponds to an exponential cutoff in the spectrum for $`h\nu >kT(r_{\mathrm{in}})`$. If $`x_{\mathrm{in}}1`$ but $`x_{\mathrm{out}}1`$, then the value of the integral is almost independent on $`\nu `$ and is a slowly varying function of $`\delta `$. We see that $`E_\nu \nu ^{(4\delta 5)/(84\delta )}`$ in this case. Thus, $`E_\nu `$ is rising for $`\delta >5/4`$ and declining for $`\delta <5/4`$. If both $`x_{\mathrm{in}}1`$ and $`x_{\mathrm{out}}1`$, then it is possible to see from the expansion of the integral in expression (49) that $`E_\nu \nu `$. At frequencies below $`h\nu =kT(r_{\mathrm{out}})`$ the whole disc surface would contribute with the low frequency tails of modified black body spectra, which scale as $`F_\nu \nu `$ (equation (47)). Therefore, it is easy to understand the scaling $`E_\nu \nu `$ for $`h\nu <kT(r_{\mathrm{out}})`$. However, we do not see the latter spectral index in calculated spectra because $`\kappa _{\mathrm{ff}}`$ becomes comparable to $`\kappa _\mathrm{T}`$ already at the frequencies where $`x_{\mathrm{out}}>1`$. The $`E_\nu \nu ^{(4\delta 5)/(84\delta )}`$ law extends down to the frequency at which $`x_{\mathrm{in}}=x_{0\mathrm{s}}(r_{\mathrm{in}})`$, where $`x_{0\mathrm{s}}={\displaystyle \frac{h\nu _0}{kT_\mathrm{s}(r)}}`$ is found by equating $`\kappa _{\mathrm{ff}}=\kappa _\mathrm{T}`$. For $`x_{0\mathrm{s}}1`$ one obtains using expression (39) for $`\kappa _{\mathrm{ff}}`$ together with expression (48) for $`T_\mathrm{s}`$ and expression (20) for $`\rho `$

$`x_{0\mathrm{s}}1.5\times 10^1\left({\displaystyle \frac{l_\mathrm{E}}{2ϵ}}\right)^{46/27}M_8^{46/27}\left({\displaystyle \frac{B_{10}}{10^4\text{G}}}\right)^{32/9}\times `$
$`\left({\displaystyle \frac{r}{10r_\mathrm{g}}}\right)^{46/932\delta /9}𝒢^{46/27}\text{.}`$ (50)
For $`x<x_{0\mathrm{s}}`$ spectrum gradually transits to the sum of a local black body with $`T=T_{\mathrm{eff}}`$, which has $`E_\nu \nu ^{1/3}`$ (Shakura & Sunyaev, 1973; Shapiro & Teukolsky, 1983; Krolik, 1999). Finally, for $`h\nu kT(r_{\mathrm{out}})`$ the spectrum is the sum of $`\nu ^2`$ low frequency black bodies of $`\nu <\nu _0`$. When the outer part of the disc, where $`\overline{\kappa }_{\mathrm{ff}}>\kappa _\mathrm{T}`$, is sufficiently truncated then the $`E_\nu \nu ^{1/3}`$ part of the spectrum may be absent and the spectrum will transit directly from $`E_\nu \nu ^{(4\delta 5)/(84\delta )}`$ power to $`E_\nu \nu ^2`$ power.

In summary, we see that magnetically dominated accretion discs have power law spectra with the spectral index depending on the radial distribution of magnetic field strength such that, $`E_\nu \nu ^{(4\delta 5)/(84\delta )}`$. This contrasts the standard weakly magnetized $`\alpha `$-disc which shows a declining modified black body formed from the inner radiation dominated disc with $`E_\nu \nu ^{2/5}`$.

### 4.2 Modified black body spectrum in a standard disc

As a side remark we note that the value for the spectral index $`\gamma ={\displaystyle \frac{\nu d\mathrm{ln}E_\nu }{d\nu }}`$ close to $`0`$ found for the latter regime of accretion disc by Shakura & Sunyaev (1973) (text on page 349 after equation \[3.11\] of that work) is different from ours $`\gamma =2/5`$. It is easy to follow the exact prescription of Shakura & Sunyaev (1973), namely calculate integral \[3.10\] in their work for spectrum \[3.2\] and temperature \[3.7\]. As a result we obtain $`\gamma =2/5`$ rather than $`\gamma =0.04`$ given in Shakura & Sunyaev (1973). We need to point out this discrepancy because it is widely stated in many textbooks on accretion discs (Shapiro & Teukolsky, 1983; Krolik, 1999) with the reference to Shakura & Sunyaev (1973) that high luminosity accretion discs has almost flat plateau in its spectrum before the exponential cut off corresponding to $`kT(r_\mathrm{s})`$. However, the flat spectrum $`E_\nu \nu ^{1/29}`$ is produced by part (b) of the standard disc model, where gas pressure dominates over radiation pressure. The spectral index $`\gamma =1/29`$ is close to the $`\gamma =0.04`$ given in Shakura & Sunyaev (1973) but the radial dependence of the surface temperature in zone (b) is $`T_\mathrm{s}r^{29/30}`$ rather than $`T_\mathrm{s}r^{5/3}`$ given by their formula \[3.7\]. Thus, the standard $`\alpha `$-disc possessing both (b) and (a) zones should have spectrum steepening from plateau to $`\nu ^{2/5}`$ and then exponentially cutting off at the temperature of the inner edge. Because the intervals of $`r`$, where approximate analytic expressions for emitted spectrum are valid, do not typically exceed one order of magnitude (the same is true for our disc model as well), one does not see “pure” extended power laws when calculating spectra numerically by using general expressions (42) and (44). For example, Wandel & Petrosian (1988) found $`\gamma =2/5`$ slope in the narrow interval between 1000 Å and 1450 Å by numerically integrating disc spectra.

### 4.3 Results of spectrum calculations

We present results of the simplified analytical integration of the spectrum using equation (47) as well as more exact numerical integration using equation (42), solving for $`T`$ from equation (45) and integrating equation (46). Function $`𝒢(r)`$ was kept in numerical calculations, so the results are applicable to the innermost parts of the disc, where the most of energy is radiated. For a given $`M_8`$ and $`l_\mathrm{E}/ϵ`$, an optically thick magnetically dominated discs exist within $`5r_\mathrm{g}<r<1000r_\mathrm{g}`$ only for $`\delta `$ in the interval of about $`1`$ to $`1.4`$. In Figs. 311 we illustrate models for the four choices of parameter sets:

1. $`M_8=1`$, $`\frac{l_\mathrm{E}}{2ϵ}=0.1`$, $`\delta =5/4`$, $`B_{10}=3\times 10^3\text{G}`$,
2. $`M_8=1`$, $`\frac{l_\mathrm{E}}{2ϵ}=0.1`$, $`\delta =1`$, $`B_{10}=3\times 10^3\text{G}`$,
3. $`M_8=1`$, $`\frac{l_\mathrm{E}}{2ϵ}=0.1`$, $`\delta =1.4`$, $`B_{10}=5\times 10^3\text{G}`$,
4. $`M_8=1`$, $`\frac{l_\mathrm{E}}{2ϵ}=10^3`$, $`\delta =5/4`$, $`B_{10}=7\times 10^2\text{G}`$.

The dependencies of $`H`$ on $`r`$ given by equation (18) for four parameter sets are plotted in Fig. 3. For comparison we also plot the half-thickness $`H`$ in the Shakura–Sunyaev model of the disc with the same accretion rate (parametrized by $`l_\mathrm{E}/ϵ`$) and the same mass of the central object. We use approximate analytic expressions for the parameters of the disc ($`H`$, $`\rho `$, $`\mathrm{\Sigma }`$, $`T_{\mathrm{mpd}}`$) in the radiation dominated zone (a) of the Shakura–Sunyaev model and thermal pressure dominated zone (b) (Shakura & Sunyaev, 1973). The magnetically dominated disc is thicker than the standard disc. For higher accretion rates, the standard disc has a concave shape due to the transition from inner zone (a) to intermediate zone (b), which allows illumination of the outer parts of the disc by the inner hotter and brighter parts. Magnetically dominated disc has convex shape, which exclude such illumination.

In Fig. 4 we plot the dependencies of the column thickness through the disc $`\mathrm{\Sigma }`$ (equation (19)) on the radius and also compare to Shakura–Sunyaev standard disc. The magnetically dominated disc is much less massive than the standard disc. Both $`\mathrm{\Sigma }`$ and $`\rho `$ are smaller for magnetically dominated discs, and only in the inner $`10r_\mathrm{g}`$ are the densities comparable.

The dependencies of mid-plane temperature $`T_{\mathrm{mpd}}`$ on radius given by equation (28) are shown in Fig. 5. On the same figure we also plot $`T_{\mathrm{mpd}}`$ in the Shakura–Sunyaev model and $`T_{\mathrm{eff}}`$ given by equation (26), which is the same for magnetically dominated and standard discs. Because of the lower column density of the magnetically dominated disc, $`T_{\mathrm{mpd}}`$ is less than for standard $`\alpha `$-discs.

The dependencies of magnetic plus turbulent pressure, $`B^2/(4\pi )`$, radiation pressure in the disc mid-plane $`P_{\mathrm{rad}}=aT_{\mathrm{mpd}}^4/3`$, and thermal pressure in the disc mid-plane $`P_{\mathrm{th}}=nkT_{\mathrm{mpd}}`$ are presented in Fig. 6. We see that the assumption of magnetic pressure dominance is well satisfied for our models except in the innermost regions, $`r10r_\mathrm{g}`$, for higher accretion rates $`l_\mathrm{E}/(2ϵ)=0.1`$, where radiation pressure becomes comparable to the magnetic pressure. The latter fact limits the existence of magnetically dominated regime in the innermost parts of accretion discs for higher luminosities. Plasma parameter $`\beta `$ defined as $`\beta =8\pi (P_{\mathrm{rad}}+P_{\mathrm{th}})/B^2`$ decreases with radius and varies from $`1`$ to $`10^2`$ in our models.

In the limit $`\beta =1`$ magnetic pressure is comparable to the largest of radiation or thermal pressures and our strongly magnetized disc model transforms into Shakura–Sunyaev model with $`\alpha =1`$. If $`\delta =3/4`$ in our model, then the radial scaling of the magnetic and turbulent pressures, $`B^2/(4\pi )`$, is the same as that of the radiation pressure inside the disc, $`aT_{\mathrm{mpd}}^4`$, aside from the factor $`𝒢`$. If $`\delta =51/40`$, then the radial scaling of $`B^2/(4\pi )`$ is the same as that of the thermal pressure inside the disc, $`nkT_{\mathrm{mpd}}`$. Therefore, by choosing $`\delta =3/4`$ and adjusting the magnitude of $`B_{10}`$ one can construct the model with approximately constant $`\beta `$ in the zone where radiation pressure exceeds thermal pressure. By choosing $`\delta =51/40`$ one can construct constant $`\beta `$ model in the zone where thermal pressure dominates radiation pressure. We illustrate this in Fig. 7, where we show the dependencies of pressures, $`\rho `$, $`H`$, and $`\mathrm{\Sigma }`$ on $`r`$ for our model with $`\delta =3/4`$, $`l_\mathrm{E}/(2ϵ)=0.5`$, and $`M_8=1`$, and for the Shakura–Sunyaev model with $`\alpha =1`$ for the same accretion rate $`l_\mathrm{E}/(2ϵ)`$ and $`M_8`$. The transition from zone (a) to zone (b) in this Shakura–Sunyaev model occurs at $`r_{\mathrm{ab}}=360r_\mathrm{g}`$. The breaks on the curves corresponding to the Shakura–Sunyaev model occur at $`r=r_{\mathrm{ab}}`$ in Fig 7. We adjusted $`B_{10}`$ such that the magnetic pressure will be in equipartition with the radiation pressure in our model. Then, as it is seen from the top-left plot in Fig. 7, thermal pressure is less than magnetic and radiation pressures for $`r`$ less than some $`r_\mathrm{c}`$ and exceeds magnetic and radiation pressures for $`r>r_\mathrm{c}`$, so our model is not applicable for $`r>r_\mathrm{c}`$. The subsequent three plots show that $`r_\mathrm{c}r_{\mathrm{ab}}`$. Two right plots and bottom-left plot in Fig. 7 show that $`\rho `$, $`H`$, and $`\mathrm{\Sigma }`$ in our model for $`r<r_\mathrm{c}`$ are very close to $`\rho `$, $`H`$, and $`\mathrm{\Sigma }`$ in Shakura–Sunyaev $`\alpha =1`$ model in the radiation pressure dominated zone $`r<r_{\mathrm{ab}}`$. A similar conclusion holds for the transition of our model with $`\delta =51/40`$ to a Shakura–Sunyaev zone (b) model for $`r>r_{\mathrm{ab}}`$ and $`\beta =1`$. The radiation spectra of our model in the limit $`\beta =1`$ also approach that of the Shakura–Sunyaev model, as shown by direct numerical calculations. The power law modified black body spectrum $`E_\nu \nu ^{(4\delta 5)/(84\delta )}`$ derived in section 4.1 becomes $`E_\nu \nu ^{2/5}`$ for $`\delta =3/4`$ and $`E_\nu \nu ^{1/29}`$ for $`\delta =51/40`$, which is coincident with the modified black body power laws for the Shakura–Sunyaev zone (a) and zone (b) spectra (see section 4.2).

The dependencies of optical depths through the half disc thickness on radius are shown in Fig. 8 for four parameter sets. The three curves plotted are: $`\overline{\tau }_{\mathrm{ff}}`$ given by equation (35), $`\tau _\mathrm{c}=\kappa _\mathrm{T}\mathrm{\Sigma }/2`$, and the effective optical thickness $`\overline{\tau }_{}={\displaystyle \frac{\mathrm{\Sigma }}{2}}\sqrt{(\kappa _\mathrm{T}+\overline{\kappa }_{\mathrm{ff}})\overline{\kappa }_{\mathrm{ff}}}`$. For $`\delta =5/4`$, the effective optical thickness $`\overline{\tau }_{}`$ is almost constant throughout the disc, but when $`\delta `$ deviates from $`5/4`$, $`\overline{\tau }_{}`$ starts to approach $`1`$ either at the inner or at the outer edge of the disc and so our model breaks down at those radii. With the decrease of the accretion rate, the disc becomes cooler and denser so the absorbing opacity rises and becomes larger than the scattering opacity in the outer parts of the disc.

In Fig. 9 we show: $`T_{\mathrm{eff}}(r)`$ given by equation (26), $`T_{\mathrm{mpd}}(r)`$ given by equation (28), $`T_\mathrm{s}(r)`$ given by equation (48), and $`T(r)`$ by solving equation (45) numerically. Note that $`T_{\mathrm{eff}}`$ and surface temperature $`T`$ are always smaller than the $`T_{\mathrm{mpd}}`$ for an optically thick disc. For low accretion rate, $`l_\mathrm{E}/(2ϵ)=10^3`$, $`\overline{\kappa }_{\mathrm{ff}}>\kappa _\mathrm{T}`$ and $`TT_{\mathrm{eff}}`$. In this case, $`T_\mathrm{s}`$ is ill defined and values of $`T_\mathrm{s}<T_{\mathrm{eff}}`$ are unphysical on the plot for $`l_\mathrm{E}/(2ϵ)=10^3`$ and also for $`r>100r_\mathrm{g}`$ on the plot for $`\delta =1`$, $`l_\mathrm{E}/(2ϵ)=0.1`$. The temperature $`T_\mathrm{s}`$ becomes a good approximation for $`T`$ when $`T_\mathrm{s}>T_{\mathrm{eff}}`$ (scattering dominates over absorption in the surface layer). Unlike the values and slope of $`T_{\mathrm{mpd}}(r)`$, which substantially increases with increasing $`\delta `$, the value of $`T`$ is less sensitive to $`\delta `$: only inner parts of the disc becomes slightly hotter for larger values of $`\delta `$. Both $`T`$ and $`T_{\mathrm{eff}}`$ are changed significantly when the accretion rate or mass $`M`$ are changed.

Fig. 10 shows the results of calculating Comptonisation $`y`$ parameter according to equations (40) and (41). We conservatively set $`x=5`$ for the calculation of $`y`$. Then, $`y`$ is the function of radius alone. On the same figure we also show $`y`$ in the regime of modified black body spectrum, using $`T_\mathrm{s}`$ and writing the simplified version of equation (40) as $`y={\displaystyle \frac{4kT_\mathrm{s}}{m_\mathrm{e}c^2}}{\displaystyle \frac{\kappa _\mathrm{T}}{\kappa _{\mathrm{ff}}}}`$. We see that Comptonisation is not important for our models even in the inner disc.

Energy spectra $`E_\nu `$ are presented in Fig. 11. We normalized frequency to the characteristic frequency of an effective black body from the inner disc, namely, we plot versus $`\nu /\nu _{\mathrm{eff}}`$, where

$$h\nu _{\mathrm{eff}}=kT_{\mathrm{eff}}(10r_\mathrm{g})=6\text{eV}\times \left(\frac{l_\mathrm{E}}{2ϵ}\right)^{1/4}M_8^{1/4}\text{,}$$

so the spectra plotted cover the range from EUV to infrared for $`M=10^8\mathrm{M}_{\mathrm{}}`$. We checked that the total thermal energy emitted from the disc between $`r_{\mathrm{in}}`$ and $`r_{\mathrm{out}}`$ calculated as an integral of the spectrum over frequencies is equal to the surface integral of the dissipation $`Q`$ (expression (12)):

$`E={\displaystyle _0^{\mathrm{}}}E_\nu 𝑑\nu ={\displaystyle _{r_{\mathrm{in}}}^{r_{\mathrm{out}}}}2\pi rQ𝑑r=`$
$`9.5\times 10^{45}\text{erg s}^1{\displaystyle \frac{l_\mathrm{E}}{ϵ}}M_8{\displaystyle \frac{r_\mathrm{g}}{r_{\mathrm{in}}}}\times `$ (51)
$`\left[1{\displaystyle \frac{r_{\mathrm{in}}}{r_{\mathrm{out}}}}{\displaystyle \frac{2}{3}}\zeta \left(\sqrt{{\displaystyle \frac{r_\mathrm{s}}{r_{\mathrm{in}}}}}{\displaystyle \frac{r_{\mathrm{in}}}{r_{\mathrm{out}}}}\sqrt{{\displaystyle \frac{r_\mathrm{s}}{r_{\mathrm{out}}}}}\right)\right]\text{.}`$
All spectra shown in Fig. 11 were computed by integrating from $`r_{\mathrm{in}}=3.1r_\mathrm{g}`$ to $`r_{\mathrm{out}}=1000r_\mathrm{g}`$ with the factor $`𝒢(r)`$ taken into account. We also show the spectrum calculated by using the approximate analytic expression (49) for a modified black body in its regime of validity ($`x_{\mathrm{in}}>x_{0\mathrm{s}}(r_{\mathrm{in}})`$). Because the analytic expression was obtained by setting $`𝒢=1`$ in equation (48) for $`T_\mathrm{s}`$, it overestimates the temperature in the inner parts of the disc by a factor of a few and does not describe the high energy part of the spectra correctly. The lower frequency at which the sum of modified black bodies is still a good approximation, increases with the overall increase of absorption in the disc. For parameter set 4 above, with $`l_\mathrm{E}/(2ϵ)=10^3`$, pure modified black body cannot be used at all, so the corresponding panel in Fig. 11 does not show a second curve. The spectrum for $`\delta =5/4`$ (solid line) shows flat plateau extending by more than an order of magnitude in frequency in accordance with analytical result. Although the declining top part of the spectrum for $`\delta =1`$ and rising top part for $`\delta =1.4`$ are apparent, the interval of frequencies, where a modified black body approximation works, becomes small and blends with the $`\nu ^{1/3}`$ spectrum of the sum of local black bodies. Thus no dependence on $`\delta `$ is evident here. The top right plot for low luminosity $`l_\mathrm{E}/(2ϵ)=10^3`$ is indistinguishable from the sum of the local black body spectra. All spectra behave like $`\nu ^2`$ for low frequencies. In Fig. 11 we also show spectra of a Shakura–Sunyaev $`\alpha `$-disc with the same $`l_\mathrm{E}/ϵ`$ and $`M_8`$ parameters to compare with plots of a magnetically dominated disc with viscosity parameter $`\alpha =0.1`$. These spectra were calculated in the same way we calculated the spectra of the magnetically dominated disc: first we found the surface temperature $`T(r)`$ by solving equation (45) with the $`\rho (r)`$ profile taken from standard $`\alpha `$-disc model, and then integrated equation (46) with $`F_\nu `$ given by expression (42). For low accretion rates of order of $`10^2`$ of Eddington accretion rate and smaller, the spectra of Shakura–Sunyaev disc are very close to the sum of local black bodies with temperatures $`T_{\mathrm{eff}}(r)`$ (Wandel & Petrosian, 1988). In general the spectra of our magnetically dominated discs are close to the spectra of Shakura–Sunyaev discs, so it seems to be difficult to distinguish between them observationally. This means that sources with discs previously thought to be thermally supported could actually be magnetically supported.

## 5 Discussion and Conclusions

We have found self-consistent solutions for thin, magnetically supported turbulent accretion discs assuming the tangential stress $`f_\varphi =\alpha (r){\displaystyle \frac{B^2}{4\pi }}`$. When compared to the standard $`\alpha `$-disc models (Shakura & Sunyaev, 1973) magnetically dominated discs have lower surface and volume densities at the same accretion rate. This is due to the more efficient angular momentum transport by supersonic turbulence and strong magnetic fields than the subsonic thermal turbulence of the standard model. As a result, magnetically dominated discs are lighter and are not subject to self-gravity instability. In the limit of plasma $`\beta =1`$, magnetic pressure is comparable to the largest of radiation or thermal pressures and our strongly magnetized disc model transforms into the Shakura–Sunyaev model with $`\alpha =1`$.

When we derived the disc structure, we made no explicit distinction between turbulent and magnetic pressure support and angular momentum transfer. As such, our model would be valid in any situation in which the magnetic and turbulent kinetic energies are comparable to, or greater than the thermal energy density. The assumption that the kinetic and magnetic energies are nearly comparable is natural because turbulence should result in the amplification of small scale magnetic fields in highly conducting medium due to dynamo action. Typically, in a sheared system, the magnetic energy can be even slightly larger than turbulent kinetic energy since the magnetic energy gains from the additional shear. We find that the thermal spectrum from the surface of the magnetically dominated disc in the optically thick regime is close to the spectrum of the standard Shakura–Sunyaev disc.

The issue arises as to how the magnetic field could reach sonic or supersonic energy densities. To obtain sonic turbulence and produce a $`\beta =1`$ disc, the MRI might be sufficient. To obtain a $`\beta <1`$ supersonic turbulence may require something else. One possibility in AGN appeals to the high density of stars in the central stellar cluster surrounding AGN accretion discs. Passages of stars through the disc might be an external source of supersonic turbulence analogous to the supernovae explosions being the source of supersonic turbulence in the Galaxy. Stars pass through the disc with the velocities of order of Keplerian velocity, which is much larger than the sound speed in the disc. We consider the support of turbulence by star-disc collisions in Appendix B and find that statistically speaking, star-disc collisions are unlikely to provide enough energy to sustain supersonic turbulence in most AGN accretion discs, however the possibility remains that a small number out of a large population could become magnetically dominated.

Indeed whether a disc could ever really attain a magnetically dominated state is important to understand. The present answer from simulations is not encouraging, but not completely ruled out. Further global MHD simulations of turbulence in vertically stratified accretion discs with realistic physical boundary conditions are needed along with more interpretation and analysis. Magnetic helicity conservation for example, has not been fully analyzed in global accretion disc simulations to date, and yet the large scale magnetic helicity can act as a sink for magnetic energy since magnetic helicity inverse cascades.

As an intermediate step in assessing the viability of low $`\beta `$ discs, it may be interesting to assess whether they are stable. One can take, as an initial condition, the stationary model of the magnetically dominated accretion disc given by expressions (15), (1820) with the initial magnetic field satisfying all constraints of our model and falling into the shadowed regions on plots in Figs. 12.

One point of note is that magnetically dominated discs may be helpful (though perhaps not essential, if large scale magnetic fields can be produced (Blackman, 2002; Blackman & Pariev, 2003)) in explaining AGN sources in which 40% of the bolometric luminosity comes from hot X-ray coronae. If the non-thermal component in galactic black hole sources is attributed to the magnetized corona above the disc (e.g., (Di Matteo et al.1999), also Beloborodov (1999) discusses possible alternatives), then magnetically dominated discs can naturally explain large fractions, up to 80 % (, Di Matteo et al.1999), of the accretion power being transported into coronae by magnetic field buoyancy (although $`\beta 1`$ disc solutions are also possible, (Merloni, 2003)). Though coronae can form in systems with high $`\beta `$ interiors, the percentage of the dissipation that goes on in the interior vs. the coronae could be $`\beta `$ dependent.

The main purpose of our study was simply to explore the consequences of making a magnetically dominated analogy to Shakura and Sunyaev, and filling in the parmeter regime which they did not consider. In the same way that we cannot provide proof that a disc can be magnetically dominated, they did not present proof that a disc must be turbulent, but investigated the consequences of their assumption. We also realize that the naive $`\alpha `$ disc formalism itself can be questioned and its ultimate validity in capturing the real physics is limited. Nevertheless it still has an appeal of simplicity.

Finally, we emphasize that our model does not describe dissipation in the corona and interaction of the corona with the disc. Further work would be necessary to address relativistic particle acceleration and emission, illumination of the disc surface by X-rays produced in the corona and subsequent heating of top layers of the disc, and emergence of magnetized outflows.

## Acknowledgments

V.I. Shishov is thanked for the remark improving presentation of the results. Two of authors (VP and EB) want to acknowledge their stay at the Institute of Theoretical Physics in Santa Barbara, where part of this work was done. SB acknowledges his stay at the University of Rochester when this work was initiated. This research was supported in part by the National Science Foundation Under Grant No. PHY99-07949. VP and EB were also supported by DOE grant DE-FG02-00ER54600.

## Appendix A On cyclotron emission

Since the characteristic temperature inside the disc ($`T_{\mathrm{mpd}}`$ given by equation (28)) is non-relativistic, cyclotron emission of an electron occurs at frequencies close to multiples of the gyrofrequency $`\omega =s\omega _B`$, where $`s=1,2,3,\mathrm{}`$, and

$$\omega _B=\frac{eB}{m_\mathrm{e}c}1.7\times 10^{11}\text{s}^1\left(\frac{B_{10}}{10^4\text{G}}\right)\left(\frac{r}{10r_\mathrm{g}}\right)^\delta \text{.}$$
(52)
However, because the magnitude of the magnetic field varies across the disc, the resulting emission blends many discrete gyrolines. The typical range of cyclotron emission is cm radio waves for AGNs and submillimetre to infrared for stellar mass black holes. Characteristic plasma parameters in our disc for the case of supermassive black holes are similar to those encountered in solar chromosphere, where plasma effects are important for the generation and propagation of radio waves (Zheleznyakov, 1970). One should expect that collective plasma effects will influence the cyclotron radiative process at such low frequencies. The plasma frequency is

$`\omega _\mathrm{L}=\sqrt{{\displaystyle \frac{4\pi e^2n}{m_\mathrm{e}}}}=8.8\times 10^{11}\text{s}^1\left({\displaystyle \frac{l_\mathrm{E}}{2ϵ}}\right)^1\left({\displaystyle \frac{B_{10}}{10^4\text{G}}}\right)^3\times `$
$`𝒢^1M_8\left({\displaystyle \frac{r}{10r_\mathrm{g}}}\right)^{33\delta }`$ (53)
and the ratio of cyclotron to plasma frequencies is

$`{\displaystyle \frac{\omega _B}{\omega _\mathrm{L}}}=2\times 10^1\left({\displaystyle \frac{r}{10r_\mathrm{g}}}\right)^{3+2\delta }𝒢\times `$
$`\left({\displaystyle \frac{B_{10}}{10^4\text{G}}}\right)^2\left({\displaystyle \frac{l_\mathrm{E}}{2ϵ}}\right)M_8^1\text{.}`$ (54)
We see that typically $`\omega _B\omega _\mathrm{L}`$, so that cyclotron emitted radiation can not propagate for some disc parameters. However, even if $`\omega _B\omega _\mathrm{L}`$, the plasma affects cyclotron radiation. As summarized in Zhelezniakov (1996) collective effects suppress the emission on the first harmonic, $`s=1`$, such that it becomes of order of the emission on the second harmonic, $`s=2`$. The emissivity on higher harmonics, $`s>2`$, is smaller by the factor $`(kT/m_\mathrm{e}c^2)^{s2}`$ than on the second and first harmonics. This occurs only at high enough plasma densities $`n{\displaystyle \frac{B^2}{4\pi }}{\displaystyle \frac{3kT}{m_\mathrm{e}^2c^4}}`$, which translates into $`c^2v_\mathrm{A}c_\mathrm{s}m_\mathrm{p}/m_\mathrm{e}`$, where $`c_\mathrm{s}3kT/m_\mathrm{p}`$. In vacuum, note that the emissivity on the first harmonic is $`(m_\mathrm{e}c^2/kT)`$ times larger than the emissivity on second harmonic. The latter condition is narrowly satisfied for small radii of our optically thick disc models and the larger the $`r`$ the better it is satisfied.

Cyclotron self-absorption also occurs in narrow lines centred on multiples of $`\omega _B`$. At some frequency $`\omega `$, emission and absorption occurs only in spatially narrow resonant layers inside the disc, where the magnetic field strength matches the frequency, i.e. $`s\omega _B(B)\omega `$ is small. The width of emission and absorption frequency intervals is determined mainly by the thermal Doppler shifts $`\mathrm{\Delta }\omega /\omega \sqrt{kT/m_\mathrm{e}c^2}`$ (Zhelezniakov, 1996). The width of such resonant layers can be estimated as $`H\sqrt{kT/m_\mathrm{e}c^2}`$. Zhelezniakov (1996) (chapter 6) gives the expression for the optical thickness through such gyroresonance layers on the second harmonic $`\omega =2\omega _B`$, which for our disc is

$`\tau _{\mathrm{cyc}}{\displaystyle \frac{\omega }{c}}{\displaystyle \frac{\omega _\mathrm{L}^2}{\omega ^2}}{\displaystyle \frac{kT}{m_\mathrm{e}c^2}}H=1.8\times 10^{12}\left({\displaystyle \frac{B_{10}}{10^4\text{G}}}\right)^4\times `$
$`\left({\displaystyle \frac{r}{10r_\mathrm{g}}}\right)^{39/84\delta }𝒢^1\left({\displaystyle \frac{l_\mathrm{E}}{2ϵ}}\right)^1M_8^{9/4}\text{.}`$ (55)
For all our models $`\tau _{\mathrm{cyc}}`$ is always very large and the emission is always strongly self-absorbed. Cyclotron photons are not subject to Compton scattering by free electrons, since the wavelength of the emission is always larger than Debeye radius in plasma, so electrostatic shielding of charges prevents them from scattering. Under such circumstances the resulting cyclotron flux from each gyroresonance layer is that of the black body with the local plasma temperature in the gyroresonance layer. Due to the overlapping of all layers, the resulting spectrum is a black body spectrum of width $`2\omega _B`$. Since $`\mathrm{}\omega _BkT`$, the total flux of cyclotron emission from the disc surface is negligibly small.

## Appendix B Star-disc collisions as possible source of turbulence

When a star passes through a disc, it creates strong cylindrical shock propagating in the surrounding gas in the disc. The aftershock gas is heated to temperatures exceeding the equilibrium temperature in the accretion disc. As the shock weakens, this heating decreases until at some distance from the impact point the incremental heating becomes comparable to the equilibrium heat content. The scale substantially affected by a star passage is $`xR_{}v_\mathrm{K}/c_\mathrm{s}`$, much larger than the radius of the star $`R_{}`$. The shock front can become unstable and turbulence can occur in the aftershock gas. The heated gas becomes buoyant, rises above the disc and falls back because of gravity. Fall-back occurs with supersonic velocities and can further excite turbulence. Turbulence will derive energy from both heating by star passages and shear of the flow. The energy, which can be derived from shear, is equal to $`Q`$ given by expression (12). It is possible that star-disc collisions might mainly be a trigger for the available shear energy to be converted into supersonic turbulence, and additional energy deposited into the disc by star-disc collisions is negligible. However, it seems unreasonable that the star-disc collisions can influence the global structure of the accretion disc unless the energy deposition from them is some fraction of the energy necessary to sustain turbulence level $`Q`$ in the disc.

The energy deposition rate by stars per unit surface of the disc is

$$Q_{}\pi R_{}^2\mathrm{\Sigma }v_{}^2\frac{1}{2}n_{}v_{}\text{,}$$
(56)
where $`v_{}v_\mathrm{K}`$ is the typical velocity of stars at radial distance $`r`$, $`R_{}R_{\mathrm{}}`$ is the average radius of stars, $`n_{}=n_{}(r)`$ is the number density of stars, $`\mathrm{\Sigma }`$ is the surface density of the disc given by expression (19) in our model. Only accretion discs orbiting supermassive black holes $`M10^6\mathrm{M}_{\mathrm{}}`$ can be influenced by star-disc collisions.

The resolution of observations is only enough to estimate the number density of stars at about 1 pc for M32 and M31 and about 10 pc for nearest ellipticals. In line with these observations we assume a star density $`n(\text{1 pc})10^410^6\mathrm{M}_{\mathrm{}}\text{pc}^3`$ at 1 pc distance from the central massive black hole (Lauer et al., 1995). To estimate $`n_{}`$ for $`r10^3r_\mathrm{g}`$ we need to rely on the theory of central star cluster evolution. The gravitational potential inside the central 1 pc will be always dominated by the black hole. Bahcall & Wolf (1976) showed that, if the evolution of a star cluster is dominated by relaxation, the effect of a central Newtonian point mass on an isotropic cluster would be to create a density profile $`nr^{7/4}`$. However, for small radii ($`0.11\text{pc}`$) the physical collisions of stars dominate two-body relaxations. Also, regions near the black hole will be devoid of stars due to tidal disruption and capture by the black hole. Numerical simulations of the evolution of the central star cluster, taking into account star-star collisions, star-star gravitational interactions, tidal disruptions and relativistic effects were recently performed by Rauch (1999). Rauch (1999) showed that star-star collisions lead to the formation of a plateau in stellar density for small $`r`$ because of the large rates of destruction by collisions. We adopt the results of model 4 from Rauch (1999) as our fiducial model. This model was calculated for all stars having initially one solar mass. The collisional evolution is close to a stationary state, when the combined losses of stars due to collisions, ejection, tidal disruptions and capture by the black hole are balanced by the replenishment of stars as a result of two-body relaxation in the outer region with a $`nr^{7/4}`$ density profile. Taking into account the order of magnitude uncertainty in the the observed star density at 1 pc, the fact that model 4 has not quite reached a stationary state can be accepted for order of magnitude estimates. For $`M=10^8M_8\mathrm{M}_{\mathrm{}}`$ we approximate the density profile of model 4 as

$`n=n_5\times 10^5{\displaystyle \frac{\mathrm{M}_{\mathrm{}}}{\text{pc}^3}}\left({\displaystyle \frac{r}{1\text{pc}}}\right)^{7/4}\text{for}r>10^2\text{pc}\text{,}`$
$`n=n_5\times 3\times 10^8{\displaystyle \frac{\mathrm{M}_{\mathrm{}}}{\text{pc}^3}}\text{for}5r_\mathrm{t}<r<10^2\text{pc}\text{,}`$ (57)
$`n\text{decreases}\text{for}r<5r_\mathrm{t}\text{,}`$
where $`r_\mathrm{t}=2.1\times 10^4\text{pc}\times M_8^{1/3}=20r_\mathrm{g}M_8^{2/3}`$ is the tidal disruption radius for a solar mass star, and $`n_5={\displaystyle \frac{n(1\text{pc})}{10^5\mathrm{M}_{\mathrm{}}/\text{pc}^3}}`$. The region $`r<r_\mathrm{t}`$ is completely devoid of stars.

We see that star-disc collisions cannot excite turbulence and strong magnetic fields in the very inner part of the accretion disc, for $`r<r_\mathrm{t}`$, and such excitation should be weak for $`r_\mathrm{t}<r<5r_\mathrm{t}`$. The relative width of the star depleted region, $`5r_\mathrm{t}/r_\mathrm{g}`$, decreases with increasing $`M`$. For $`M2\times 10^9\mathrm{M}_{\mathrm{}}`$ $`r_\mathrm{t}=3r_\mathrm{g}`$ and star-disc collisions happen all over the disc. For $`M<3\times 10^6\mathrm{M}_{\mathrm{}}`$ $`5r_\mathrm{t}>10^3r_\mathrm{g}`$ and for $`M<3\times 10^5\mathrm{M}_{\mathrm{}}`$ $`r_\mathrm{t}>10^3r_\mathrm{g}`$ and star-disc collisions are unimportant for the structure of the accretion disc. Let us assume that it should be $`Q_{}=fQ`$ where the fraction $`f`$ is less than unity but not much less than unity. Further, we use expression (12) for $`Q`$, expression (19) for $`\mathrm{\Sigma }`$, and the value for $`n`$ in the constant density core of the star cluster, second raw in expression (57), to substitute into $`Q_{}=fQ`$. Since the relation $`Q_{}=fQ`$ should be satisfied for all values of $`r`$, the value of $`\delta `$ is determined and turns out to be $`\delta =3/2`$. Solving the rest of the equation for the magnitude of magnetic field $`B_{10}`$ at $`\delta =3/2`$ we obtain

$$\frac{B_{10}}{10^4\text{G}}=6.5\times 10^4\left(\frac{l_\mathrm{E}}{2ϵ}\right)^{1/2}M_8^{3/4}n_5^{1/4}f^{1/4}\text{.}$$
(58)
This value of $`B`$ required by energy input from star-disc collisions should fall into the range of constraints for the $`B_{10}`$ listed in the end of section 3.

We explored all feasible range of parameters $`M_8`$, $`l_\mathrm{E}/ϵ`$, $`f>10^3`$, $`n_5<10^3`$ and found that the magnetic field calculated from expression (58) is always too strong to fall in the allowed range of parameters discussed at the end of section 3. In particular, the constraint that magnetic and turbulent pressure dominate thermal and radiation pressure is violated. The minimum number density of stars necessary to satisfy this constraint at the most favourable values of other parameters still plausible for some AGN ($`M_8=40`$, $`l_\mathrm{E}/ϵ=10^{10}`$, $`f=10^3`$), turns out to be $`n_510^6`$. Such a high number density of stars would imply total mass in stars of order of $`10^{11}\mathrm{M}_{\mathrm{}}`$ inside the central parsec from the central black hole. This mass exceeds observational and theoretical limitations.