Iron Opacity Bump Changes the Stability and Structure of Accretion Disks in Active Galactic Nuclei

1.5 MB·

by  Yan-Fei Jiang

Draft version January 27, 2016
PreprinttypesetusingLATEXstyleemulateapjv.08/22/09
IRON OPACITY BUMP CHANGES THE STABILITY AND STRUCTURE OF ACCRETION DISKS IN ACTIVE
GALACTIC NUCLEI
Yan-Fei Jiang(姜燕飞)1,4, Shane W. Davis2 & James M. Stone3
1Harvard-SmithsonianCenterforAstrophysics,60GardenStreet,Cambridge,MA02138,USA
2DepartmentofAstronomy,UniversityofVirginia,P.O.Box400325,Charlottesville,VA22904-4325,USAand
3DepartmentofAstrophysicalSciences,PrincetonUniversity,Princeton,NJ08544,USA
Draft version January 27, 2016
ABSTRACT
6 AccretiondisksaroundsupermassiveblackholeshaveregionswheretheRosselandmeanopacitycan
1 bemuchlargerthantheelectronscatteringopacityprimarilyduetothelargenumberofbound-bound
0 transitions in iron. We study the effects of this iron opacity “bump” on the thermal stability and
2
verticalstructureofradiationpressuredominatedaccretiondisks,utilizingthreedimensionalradiation

n magneto-hydrodynamic simulations in the local shearing box approximation. The simulations self-
a consistently calculate the heating due to MHD turbulence caused by magneto-rotational instability
J and radiative cooling by using the radiative transfer module based on a variable Eddington tensor
5 in Athena. For a 5×108 solar mass black hole with ∼ 3% of the Eddington luminosity, a model
2 including the iron opacity bump maintains its structure for more than 10 thermal times without
showing significant signs of thermal runaway. In contrast, if only electron scattering and free-free

] opacity are included as in the standard thin disk model, the disk collapses on the thermal time scale.
E
The difference is caused by a combination of 1) an anti-correlation between the total optical depth
H andthetemperature,and2)enhancedverticaladvectiveenergytransport. Theseresultssuggestthat
. the iron opacity bump may have a strong impact on the stability and structure of AGN accretion
h
disks, and may contribute to a dependence of AGN properties on metallicity. Since this opacity is
p
relevant primarily in UV emitting regions of the flow, it may help to explain discrepancies between
-
o observation and theory that are unique to AGNs.
r Subject headings: accretion disks − (magnetohydrodynamics:) MHD − methods: numerical − radia-
t
s tive transfer
a
[
1. INTRODUCTION Stella & Rosner 1984; Merloni 2003). The first 3D local
1
shearing box simulations with turbulence driven by the
v The radiation from active galactic nuclei (AGNs) is
magneto-rotational instability (MRI, Balbus & Hawley
6 widely believed to be produced from the gas accreted
1991) and optically thick radiative cooling seemed con-
3 by supermassive black holes through their accretion
sistent with thermal stability (Turner 2004; Hirose et al.
8 disks (e.g., Malkan 1983). The standard α disk model
2009b), but more recent simulations performed in larger
6 (Shakura & Sunyaev 1973) is commonly used to model
domains and using more sophisticated radiative transfer
0 the structures of accretion disks in most AGNs, mainly
methodshavefoundthattheradiationdominatedaccre-
. because the temperature predicted by the α disk model
1 tiondisksstillshowsignsofthermalrunawayafterafew
is roughly consistent with the big blue bump in AGN
0 thermal time scales (Jiang et al. 2013a). However, no
spectra(Elvisetal.1994;Koratkar&Blaes1999). How-
6 clear evidence of rapid variability that might be consis-
ever,therearediscrepancieswhendetailedpredictionsof
1
tent with the thermal instability has yet been observed
: the α disk model are compared with AGN observations
v (see e.g. Koratkar & Blaes 1999). in most AGNs. Although there has been speculations
Xi From the theoretical point of view, the main puzzle is that recently discovered “change-look” AGNs (LaMassa
et al. 2015; Ruan et al. 2015) may be caused by these
theoutcomeoftheinflow(“viscous”)andthermalinsta-
r instabilities, it is difficult to understand why the insta-
a bilities in α disks that are radiation pressure dominated,
bilities would uniquely manifest themselves only in this
which should be the case for the inner region of AGN
specific subset of AGNs.
disks (Lightman & Eardley 1974; Shakura & Sunyaev
Observationally, many properties of AGNs cannot be
1976). The thermal instability, which should grow faster
explained with the standard α disk model (Koratkar &
thantheinflowinstability,existsbecausethecoolingrate
Q− in this regime is proportional to the midplane pres- Blaes 1999), in spite of the issue of thermal instability.
sureP whiletheheatingrateQ+ ∝P2 (Piran1978). TheSpectrumEnergyDistribution(SED)ofmostAGNs
z,0 z,0 showsaturnoveralwaysaround1000˚A,almostindepen-
When the pressure is perturbed to be smaller (larger)
than equilibrium value, Q+ decreases (increases) faster dent of the black hole mass (Zheng et al. 1997; Bonning
compared with the change of Q− and the disk will con- et al. 2007; Davis et al. 2007; Laor & Davis 2014, and
references therein). The predicted Lyman edge is also
tinue to cool down (heat up) and never return to the
not observed (Shull et al. 2012). The inferred size of the
originalequilibriumstate. Theexistenceoftheradiation
AGN accretion disks based on micro-lensing measure-
dominated thermal instability has been questioned (e.g.
ments (e.g. Morgan et al. 2010; Blackburne et al. 2011)
4EinsteinFellow or the lags between the variations in different continu-
2 Y.-F. Jiang et al.
ous bands (Edelson et al. 2015) is systematically larger 2. SIMULATIONSETUP
than the half-light radius of the α disk model. A soft X-
WesolvethesamesetofradiationMHDequationsun-
ray excess is also ubiquitously observed in bright AGNs
derthelocalshearingboxapproximationasequation(2)
(Crummy et al. 2006; Done et al. 2012; Done 2014), but
of Jiang et al. (2013a) given by. The simulation box is
is not easily explained by the α disk model if it results
located at a fiducial radius r = 20 Schwarzschild radii
from continuum emission. from a M = 5 × 108M 0black hole. Keplerian ro-
The order-of-magnitude temperature and density in BH (cid:12)
tation is assumed for the background flow with the or-
the AGN accretion disks can actually be roughly es- bital frequency Ω = 1.60×10−6 s−1 and shear parame-
timated based on a few assumptions independent of
ter q = −dlnΩ/dlnr = 3/2. We solve these equations
any accretion disk model. If we assume the accreted
withtheGodunovradiationMHDcodeasdescribedand
gas is optically thick with typical Eddington luminosity
testedinJiangetal.(2012)andDavisetal.(2012),with
L = 1.5 × 1046M /(108M ) erg/s, typical emis-
Edd BH (cid:12) the improvements given by Jiang et al. (2013b). Under
sion size to be the Schwarzschild radius r = 3.0 ×
s the shearing box approximation, heating is generated by
1013M /(108M ) cm, then the effective temperature
BH (cid:12) the dissipation from MRI turbulence with the energy ul-
is T ∼ 3.9×105(cid:0)M /108M (cid:1)−1/4 K. Density of the timately coming from work done by the shearing period
BH (cid:12)
gas depends on the assumed inflow velocity, which is boundary condition, while cooling is dominated by the
likelysubsonicwithrespecttotheradiationsoundspeed photons leaving from the top and bottom of the simula-
(Pringle 1981). As an order-of-magnitude estimate, we tion box.
take it to be the gas isothermal sound speed with the For the opacities describing the interactions between
effective temperature estimated above. With the typ- the radiation and gas, we calculate the total Rosseland
ical mass accretion rate M˙Edd = 10LEdd/c2 = 1.6 × mean opacity κt based on the opacity table from the
1026M /(108M ) g/s and inflow radius r , the typi- Modules for Experiments in Stellar Evolution (MESA,
BH (cid:12) s
cal density is ρ ∼ 2×10−9 g/cm3 for the 108M black Paxton et al. 2011) for solar metallicity in order to cap-
(cid:12) turethedependenciesoftheopacityontemperatureand
hole. The density and temperature in the α disk model
density correctly. This opacity table is originally from
are consistent with this simple estimate. However, this
the OPAL opacity project (Iglesias & Rogers 1996) and
temperature and density regimes are very similar as in
has been successfully used to study the stability and
the envelope of massive stars (Paxton et al. 2011; Jiang
structures of massive star envelopes (Jiang et al. 2015).
et al. 2015). With non-negligible metallicity, the OPAL
The opacity as a function of temperature and density is
opacity project (Iglesias & Rogers 1996) shows that the
showninFigure2ofJiangetal.(2015),whichalsocovers
Rosseland mean opacity can be significantly enhanced
the relevant parameter space for AGN disks. The most
compared with the electron scattering value due to met-
important feature is the opacity bump caused by a large
als. Inparticular,theironsproducethewell-knownopac-
itybumparoundthetemperature1.8×105 K.Theopac- number of bound-bound transitions, most due to iron
atoms,fortemperaturesaround1.8×105 K.Thiscanbe
ity drops rapidly with increasing or decreasing tempera-
a factor of ∼ 4 larger than the electron scattering opac-
ture and it only depends on the density weakly (Figure
ity for solar metallicity. In order to split this Rosseland
2 of Jiang et al. 2015). As significant metals are indeed
mean opacity into scattering and absorption terms (as
observed in the broad line regions of AGNs (Hamann
only the absorption opacity enters the energy exchange
& Ferland 1993; Dhanda et al. 2007), the iron opacity
terms in the radiation moment equations), we simply
bumpshouldexistintheAGNdisks. However,thestan-
dard α disk model only includes the electron scattering adopttheelectronscatteringopacityκes =0.34cm2 g−1
andfree-freeopacities. FollowingJiangetal.(2013a),we andsubtractthisfromtheRosselandmeanopacity. The
willcalculatetheturbulencefromMRIwithoutanyαas- Planck mean absorption opacity is taken to be the same
sumption and radiative transfer including the iron opac- as Rosseland mean absorption opacity for simplicity, al-
ity bump self-consistently based on 3D radiation MHD thoughitislikelytobeanunderestimatesincetheRosse-
simulations. Because optical depth is critical for the ra- landmeantendstomoreheavilyweightfrequencieswith
diativecoolingofthedisk,wewillstudyhowthethermal lowopacity. Forcomparison,wealsoperformsimulations
stabilityandstructuresoftheAGNdiskswillbeaffected with the “standard” thin disk model opacity, which only
by the iron opacity bump. Note that for accretion disks includes the electron scattering opacity κes, Plank-mean
in X-ray binaries around stellar mass black holes, the free-free absorption opacity κ =3.7×1053(cid:0)ρ9/E7(cid:1)1/2
aP g
disk midplane temperature in the inner region is too hot cm2 g−1 and Rosseland-mean free-free absorption opac-
for iron opacity to play an important role, but it may ity κ = 1.0×1052(cid:0)ρ9/E7(cid:1)1/2 cm2 g−1. Here E =
playaroleintheouterdiskdynamics. Anotherexample aF g g
where opacity effects significantly change the dynamics Pg/(γ −1) is the gas internal energy with gas pressure
of the disk is the hydrogen ionization instability in ac- Pg, density ρ and adiabatic index γ = 5/3. This combi-
cretion disks around white dwarfs (Lasota 2001; Hirose nation of opacities has been used in most previous simu-
et al. 2014). lations, whichfocussedonhotterdisksmoreappropriate
Thispaperisorganizedasfollows. Wedescribehowwe to∼10M(cid:12)blackholeX-raybinaries(Hiroseetal.2009b;
setup the simulations in § 2, and the initial and bound- Jiang et al. 2013a).
aryconditionsweusein§3. Ourprimaryresultsarede-
scribed in § 4, while § 5 compares the simulation results 3. INITIALANDBOUNDARYCONDITIONS
with the α disk model and discusses the implications for We construct the initial vertical profiles of the disk
AGN observations. based on hydrostatic equilibrium and diffusion equation
as in Hirose et al. (2009b) and Jiang et al. (2013a) but
3
withthetheironopacitybumpincludedself-consistently.
Table 1
We assume the dissipation profile dF /dz ∝ ρκ /τ0.5,
r,z t SimulationParametersoftherunOPALR20
where τ is the optical depth from the nearest surface of
the disk and F is the vertical component of the ra- Ω/s−1 1.60×10−6
r,z
diation flux. If the total optical depth from the disk r0 /cm 2.97×1015
midplane to the surface is τ , by symmetry, the radia- ρ0/gcm−3 1.00×10−8
tion flux as a function of τ0within the photosphere is T0 /K 2.00×105
F = F (cid:0)τ0.5−τ0.5(cid:1)/(cid:0)τ0.5−1(cid:1). Here we choose P0 /dyncm−2 2.77×105
r,z max 0 0 Hs /cm 2.81×1013
the initial maximum radiation flux Fmax = 7.92×1011 Σ/gcm−2 4.34×105
erg s−1 cm−2, and Fr,z is fixed to this value in the re- τ0 2.31×105
gion where τ < 1. We choose the midplane tempera- Fmax/ergs−1 cm−2 6.36×1012
ture 2.4 × 105 K and integrate vertically according to Note: TheparametersΩ,r0andΣarefixedforthesimulation,
the diffusion equation dEr/dτ = 3Fr,z/c, where c is the while τ0 and Fmax are time averaged properties between 60
speed of light. The total optical depth τ is chosen such and 125 orbits. The density ρ0, pressure P0, temperature T0
√ 0 and scale height Hs are the fiducial units we use to describe
that at τ = 1, Er = 3cFr,z. Initially, gas tempera- thesimulation. Theyareprettyclosetothemidplanedensity,
ture T is set to be the same as the radiation temper- temperatureandcharacteristicdiskscaleheight.
ature Tr ≡ (Er/ar)0.25, where the radiation constant 4.1. Simulation History
a =7.57×1015ergcm−3K−4. Theinitialdensityprofile
r
isconstructedbasedontheequationofhydrostaticequi-
librium dP/dz =κ F /c−Ω2z and dτ =−ρκ dz. The
t r,z t
midplane density is adjusted to be 10−8 g cm−3 such
that the initial total optical depth is τ = 9.5 × 105.
0
All the quantities above the photosphere are fixed to
be the same values as at τ = 1. Notice that in the
initial condition, only the surface density Σ is the con-
served quantity, while all the other quantities such as
ρ,T,τ ,F will adjust self-consistently during the sim-
0 r,z
ulation. The magnetic field is initialized in the same
way as in Jiang et al. (2013a) with the initial ratio be-
tween gas pressure and magnetic pressure to be 12 at
z =0. The boundary conditions are also the same as in
Jiang et al. (2013a). For the short characteristic mod-
ule we use to calculate the variable Eddington tensor
(VET), we use 80 angles per cell to capture the angular
distribution of the radiation field. Sizes of the simula-
tion box are all fixed to be L = 0.87H , L = 3.48H
x s y s
and L = 6.96H , where H is the length unit listed
z s s
in Table 1. The length unit is chosen based the total
radiation flux F we get from the simulation (Table
max
1) as H = κ F /(cΩ2). For the typical density ρ
s es max 0
and temperature T given in table 1, it is related to the
0
gas pressure scale height H ≡c /Ω and radiation pres-
g g
sure scale height H ≡ c /Ω as H = 8.57H = 2.25H , Figure 1. Top: historiesofthevolumeaveragedradiationenergy
r r s g r
where c is the isothermal sound speed for temperature densityEr (redline),gasinternalenergyEg (blackline)andmag-
g (cid:112) neticenergydensityEB (blueline)fortherunOPALR20withiron
T and radiation sound speed c = a T4/(3ρ ). We opacity bump. Middle: history of the volume averaged Maxwell
0 r r 0 0
use 64×128×512 grids for x,y,z directions so that we stress−BxBy (blackline)andReynoldsstressρvxδvy (blueline).
Bottom: historyoftheαparameter,whichistheratiobetweenthe
have roughly 32 grids per radiation pressure scale height
sumofthetotalMaxwellandReynoldsstressandthesumofthe
Hr. Following the convention in Athena (Stone et al. total radiation, gas and magnetic pressure. The vertical dashed
2008), unit of the magnetic field is chosen so that mag- lineseparatesthefirst60orbitswhenEddingtonapproximationis
netic permeability is one. adopted and the time with VET turned on. Units for the energy
densityandstressareP0 asgiveninTable1.
4. RESULTS
The initial evolutions of the disks are very similar for 4.1.1. The Run OPALR20 with Iron Opacity Bump
thecaseswithorwithouttheironopacity. Thediskcools For OPALR20, we first ran the simulation by setting
down slightly while MRI is still growing from the lami- the Eddington tensor f = 1/3I and disk lasted for 60
narinitialconditionduringthefirstfeworbits. However, orbits. Then the simulation is restarted with the short
once heating is generated by vigorous MHD turbulence characteristic module turned on to calculate the VET
from MRI, the disk undergoes quite different evolution self-consistently for another 75 orbits. Histories of the
histories for different opacities. We label the run with volume averaged energy densities for the whole simula-
iron opacity bump as OPALR20, while three runs with tion duration are shown in the top panel of Figure 1.
just electron scattering and free-free opacities for com- Although we have run the simulation for more than 10
parison as ESR20a, ESR20b and ESR20c. thermal time scales and E is more than 60 times larger
r
4 Y.-F. Jiang et al.
Jiang et al. (2013a), where Q+ and Q− diverge from
3 each other when the disk undergoes thermal runaway.
2 100 Because the disk is in thermal equilibrium, we also can-
1 10-1 not measure the dependence of Q+ and Q− on Pz,0 as
/Hs 0 10-2ρ we did in Jiang et al. (2013a).
z
1
10-3
2
3 10-4
3 2.5
2 2.0
1.5
1 1.0
/Hs 0 00..05 By
z 1 0.5
1.0
2 1.5
2.0
3 2.5
60 70 80 90 100 110 120 130
t/Orbit
Figure 2. Space-timediagramofthedensityρ(toppanel,inunit
o√f ρ0) and azimuthal magnetic field By (bottom panel, in unit of
2P0)fortherunOPALR20.
than E , the radiation energy density E , gas internal
g r
energy E and magnetic energy density E do not show
g b
any thermal runaway behavior as shown in Jiang et al.
(2013a). The radiation energy density E varies only by
r
afactorof2whileE isalmostaconstant. Magneticen-
g
ergy density E shows a larger variation amplitude and
b
it can change by a factor of 5 after the initial 60 or-
bits. HistoriesofthevolumeaveragedMaxwellstressand
Reynolds stress are shown in the middle panel of Figure Figure 3. Top: change of the heating Q+ and cooling Q− rates
1, while history of the equivalent α parameter is shown per unit area as a function of the midplane total pressure Pz,0
in the bottom panel of Figure 1. Here α is calculated as for the run OPALR20. Bottom: change of the total optical depth
the ratio between the time and volume averaged stress τ0 and flux weighted optical depth τ¯ as a function of Pz,0. The
midplanepressureisinunitofP0 whiletheunitsforQ+ andQ−
andtotalpressure,whichis0.025afterthefirst60orbits. areP0HsΩ.
The average ratio between the total Maxwell stress and
Reynolds stress from the MRI turbulence is 4.33 while
the average ratio between the total Maxwell stress and
4.1.2. Three Runs Without the Iron Opacity Bump
magneticpressureis0.25. Thesestatisticalpropertiesare
consistentwithpreviouslocalshearingboxorglobalsim- Toconfirmthatthedifferentbehaviorsweseebetween
ulationsofMRIturbulence,eitherwithisothermalequa- the run OPALR20 and the simulations shown in Jiang
tion of state or self-consistent radiative transfer (Turner et al. (2013a) are indeed caused by the iron opacity
etal.2003;Guanetal.2009;Hawleyetal.2011;Sorathia bump, we have done three comparison runs by only in-
et al. 2012; Jiang et al. 2013b). cluding the electron scattering and free-free opacities as
The space-time diagram of the density ρ and toroidal in Jiang et al. (2013a). For the run ESR20a, we use the
magneticfieldB forthissimulationafterthefirst60or- same surface density as in OPALR20. Because electron
y
bitsareshowninFigure2,wherethewell-knownbutter- scatteringopacityissmallerthantheironopacity,theto-
fly diagram is clearly observed. The turbulent magnetic talopticaldepthτ isonly34%ofthevalueinOPALR20.
0
fieldgeneratedbyMRIpeaksaroundz ≈±H . Thebut- For the run ESR20b, we increase the surface density by
s
terfly pattern can be attributed to regions of strong B a factor of 2 so that the total optical depth is closer to
y
buoyantly rising away from the midplane. This pattern the value in OPALR20. We do not increase the surface
of B reverses roughly every 10 orbits. The buoyantly density to match the τ in OPALR20, because then the
y 0
risingmagneticfieldprovidesenhancedpressuresupport surface density is larger than the maximum surface den-
and causes the density scale height near the surface to sityallowedbythethindiskmodelwithα=0.02−0.03.
rise and fall following the same pattern. These buoyant Initial conditions for the two runs are constructed in the
motions also affect the energy transport inside the disk, same way as described in Section 3. To test the effect of
and are discuss further in section 4.2. the initial condition, for the run ESR20c, we restart the
To confirm that the disk is in roughly thermal equi- simulation OPALR20 at 60 orbit by changing the opac-
librium, we compare the total heating Q+ and cooling ity to be electron scattering and free-free opacities while
Q− rates according to equations (2) and (3) in Jiang keeping all the other quantities unchanged. In this way,
et al. (2013a), which are shown in the top panel of Fig- ESR20c has exactly the same turbulence as OPALR20 to
ure 3. Indeed, Q+ and Q− track each other very well start with. All the other parameters of the three runs,
and midplane pressure only varies in a very small dy- such as box size and resolution, are the same as the run
namic range. This is quite different from Figure 2 of OPALR20.
5
Figure 4. Histories of the volume averaged energy densities Er Figure 6. ChangeoftheheatingQ+,coolingQ− ratesperunit
(redlines),Eg (blacklines)andEB (bluelines). Fromthetopto area(top)andtotalopticaldepthτ0aswellasfluxweightedoptical
bottom,theyareforthethreerunsESR20a,ESR20bandESR20c, depthτ¯(bottom)asafunctionofthemidplanetotalpressurePz,0
which only include the electron scattering and free-free opacities. fortherunESR20c. TheunitsarethesameasinFigure3.
UnitsoftheenergydensitiesareP0.
declines while the disks collapse. They do not reach any
radiation pressure dominated thermal equilibrium state
as the run OPALR20. Instead, they behave in a very
similar way as the simulation RSVET shown in Jiang
et al. (2013a). For the run ESR20c where heating from
the MRI turbulence exists from the beginning, the ini-
tialradiationenergydensityofthediskisprettycloseto
value as predicted by the radiation pressure dominated
standard thin disk solution with the same surface den-
sity and an equivalent α = 0.02. However, with only
electron scattering and free-free opacities, the disk does
not adjust itself to reach a radiation pressure dominated
equilibriumstate. Instead,E dropsbyoneorderofmag-
r
nitudecontinuouslywithin30orbits. TheMaxwellstress
alsodecreaseswhilethediskcollapses. Thedependences
ofheatingQ+ andcoolingQ− rateonthemidplanepres-
sure for the run ESR20c are shown in the top panel of
Figure6,whichareverysimilartothesimulationRML-
VET reported by Jiang et al. (2013a). The heating rate
doeshaveastrongersensitivitytothemidplanepressure
compared with the cooling rate when the thermal run-
awayhappens. ComparedwithOPALR20,thethreeruns
confirm that the different behaviors are indeed caused
by the iron opacity bump, as the opacity law is the only
difference between them.
Figure 5. HistoriesofthevolumeaveragedMaxwellstressforthe
three runs ESR20a, ESR20b and ESR20c. Units of the stress are
P0. 4.2. Vertical Structure of the Disk in the Run OPALR20
Inordertoinvestigatethereasonswhytheironopacity
Histories of the volume averaged energy densities and bump can make the radiation pressure dominated disks
Maxwell stress of the three runs are shown in Figure 4 lastmuchlonger,wefirststudythetimeaveragedvertical
and 5. For all the three cases, the disks continue to cool structures of disk in the run OPALR20. We compute
down and collapse within ∼ 4−6 thermal time scales. timeaveragesstartingat60orbitssothatonlytheVET
For ESR20a and ESR20b where we start the simulations portion of the runs is included.
from the laminar state, they collapse more quickly be- The horizontally averaged vertical profiles of ρ, T , T
g r
cause thereis noheating atthe beginning. The Maxwell andκ areshowninFigure7. Themidplanetemperature
t
stress reaches the peak within the initial ∼6 orbits and islargerthanT sothatthepeakoftheFeopacitybump
0
6 Y.-F. Jiang et al.
the radiation acceleration a . The time averaged verti-
r
cal profiles of these accelerations are shown in the top
panel of Figure 8, which shows that equation (1) is in-
deedsatisfiedverywell. Theverticalprofilesofradiation
pressure P , gas pressure P and magnetic pressure P
r g b
areshowninthebottompanelofFigure8. Nearthedisk
midplane when |z| (cid:46) 1.8H , the disk is radiation pres-
s
suredominated, which isalso consistentwith thefact a
r
is much larger than a and a in this region. Beyond
gas B
that, magnetic pressure becomes dominant. The mag-
netic pressure first increases with distance near the disk
midplaneandpeaksaround±0.7H . Thenitdropswith
s
height. This feature is commonly observed in previous
MRIsimulations(e.g.,Miller&Stone2000;Hiroseetal.
2006; Blaes et al. 2011; Jiang et al. 2014b). The ratio
between radiation pressure and gas pressure varies from
∼15 to more than 100.
Figure 7. Timeandhorizontallyaveragedverticalprofilesofden-
sityρ(toppanel),gas(Tg)andradiation(Tr)temperature(middle
panel), as well as the Rosseland mean opacity κt (bottom panel)
fortherunOPALR20.
occurs off the midplane. The opacity κ peaks around
t
≈0.5H anddropsbothwhenthetemperaturedecreases
s
towards the photosphere and when it increases towards
the midplane. At the peak, κ is more than three times
t
the value of the electron scattering opacity for the solar
metallicity we adopt. The iron opacity has a very weak
dependence on density (Figure 2 of Jiang et al. 2015).
The rapid drop around ±H is primarily because of the
s
drop in temperature.
Because local dynamic time scale 1/Ω is much shorter
than the thermal time scale, we expect the hydrostatic
equilibrium to be maintained very well. This means the
timeaveragedverticalaccelerationsduetovariousforces
should be roughly balanced as
Figure 8. Top: timeandhorizontallyaveragedverticalprofilesof
ag =ar+agas+aB, (1) accelerations due to gas agas, radiation ar and magnetic field aB
fortherunOPALR20. Thesumofthethreeaccelerationsisasum
where ag = Ω2z is the vertical component of the gravi- whileag istheverticalcomponentofthegravitationalacceleration
tational acceleration due to the central black hole. The duetothecentralblackhole. Bottom: averagedverticalprofilesof
accelerations due to radiation (a ), gas (a ) and mag- gaspressurePg,radiationpressurePr andmagneticpressurePB.
netic field (a ) are calculated asr gas Units for the accelerations are Ω2Hs while the pressure units are
B P0.
κ F
a = t r,z0,
r c
4.2.1. Change of the Optical Depth with Temperature
∂P
a =− g,
gas ρ∂z Becauseofthesensitivedependenceoftheironopacity
bumpwithtemperature,itwillonlyenhancetheopacity
a =− 1 ∂ (cid:0)B2+B2(cid:1) significantly within a narrow temperature range, which
B 2ρ∂z x y corresponds to a certain vertical range of the disk. The
1(cid:18) ∂B ∂B (cid:19) contribution of the iron opacity bump to the total opti-
+ρ Bx ∂xz +By ∂yz . (2) caldepthisdominatedbyτp =ρpκpHp,whereρp andκp
arethedensityandopacityattheopacitypeakwhileH
p
Here B and B are the radial and vertical components is proportional to the disk scale height. For a constant
x z
of the magnetic field and a includes both the accelera- surfacedensity,thiscanberewrittenasτ =ρ κ Σ/ρ ,
B p p p z,0
tionsduetomagneticpressuregradientandverticalcom- whereρ isthedensityatdiskmidplane. Whenthedisk
z,0
ponent of magnetic tension (Blaes et al. 2011). Notice becomeshotterandmidplanetemperatureincreases,the
that only the diffusive radiation flux F contributes to iron opacity peak will move to larger height and thus
r,z0
7
ρ /ρ decrease. When the disk becomes colder and
p z,0
midplane temperature drops, the iron opacity peak will
move towards the midplane and ρ /ρ increases. Be-
p z,0
cause κ only depends on density weakly, as long as the
p
whole iron opacity bump is included in the disk and it is
thedominantopacity, thismeansthetotalopticaldepth
willdecreasewhenthediskbecomeshotter,andincrease
with the disk becomes colder.
To confirm this, we plot the total optical depth 2τ
0
as a function of the midplane pressure (dominated by
P ) for the run OPALR20 after the initial 60 orbits at
r
the bottom panel of Figure 3, which shows a very clear
anti-correlationbetweenτ andP . Asimplepowerlaw
0 z,0
fitting shows that τ ∝ P−0.89 in this simulation. This
0 z,0
is completely different from the electron scattering dom-
inated case as shown in the bottom panel of Figure 6 for
the run ESR20c, where τ is almost independent of P
0 z,0
as expected. This is critical for the thermal instability,
becauseinthestandardthindiskmodel,thecoolingrate
Q− ∝ P /τ . With this close to linear anti-correlation
z,0 0
between τ and P , the sensitivity of the cooling rate
0 z,0
to midplanepressure is enhanced. This isone important
reason why the thermal instability is suppressed and the
disk can last much longer. Figure 9. Top: timeandhorizontallyaveragedverticalprofilesof
One crude way to understand the anti-correlation be- labframeradiationfluxFr,z,diffusiveradiationfluxFr,z0,advec-
tweenτ andP quantitivelyisbyassumingaconstant tiveradiationfluxFadv andPoyntingfluxFB. Bottom: averaged
0 z,0 vertical profiles of local dissipation rate. The solid black line is
ratio between radiation pressure and gas pressure. This
theworkdonebythestressateachheightqΩ(−BxBy+ρvxδvy).
isexpectedwithefficientconvectionintheradiationpres- The solid red, blue and green lines are dFr,z/dz, dFr,z0/dz and
sure dominated regime so that the radiation entropy per dFB/dz,respectively. ThedashedblacklineisthesumofdFr,z/dz
unit mass, which is proportional to P /P , is roughly a and dFB/dz. The dashed red line is the critical dissipation rate
constant. Thisisclearlynotthecaseforrthgewholediskas cΩ2/κes,whilethedashedbluelineiscΩ2/κt. Unitsfortheenergy
fluxanddissipationarecgP0 andcgP0/Hs.
this ratio increases from ∼15 near the disk midplane to
more than 100 near the surface. However, in the region Another important component of energy flux is the
when |z| (cid:46) 0.5H , P /P does not vary too much and Poynting flux
s r g
this is not a very bad assumption. Once we adopt this
F =B2v −(B·v)B . (5)
assumption, we have ρ /ρ = (T /T )3, where T B z z
p z,0 p z,0 z,0
is the midplane temperature. Because the iron opacity The mechanical energy flux and gas enthalpy flux are
peak only exists for a roughly fixed temperature T and small and we will not show them here. Time averaged
p
P ∝ T4 in the radiation pressure dominated regime, vertical profiles of these energy fluxes are shown in the
z,0 z,0
the total optical depth contributed by the iron opacity toppanelofFigure9. Thelocaldissipationratesateach
peak can be roughly estimated as τp ∝ Pz−,00.75, which is haneidghqt−a=redsFim/pdlyz.qTr−he=wdoFrkr,zd/odnze, bdyFrt,hz0e/sdtzr,esdsFaatdve/adchz
actuallyprettyclosetowhatwegetfromthesimulation. B B
height is q+ = qΩ(−B B +ρv δv ). These local heat-
x y x y
ing and cooling rates are shown in the bottom panel of
4.2.2. Vertical Energy Transport
Figure 9. In the thermal equilibrium state, local energy
AsshowninthetoppanelofFigure3,duringthether- conservation requires
mal equilibrium state, the total heating rate Q+ is bal-
anced by the total cooling rate Q−, where Q− is dom- q+ =qr−+qB−. (6)
inated by radiative cooling. In the standard thin disk
Timeaveragedverticalprofilesoftheleftandrighthand
model, only the diffusive radiation flux is considered for
sides are shown as the solid and dashed black lines in
theverticalradiativecoolingeverywhere. Thisisnotnec-
the bottom panel of Figure 9, which shows that they do
essarythecaseasphotonscanalsobeadvectedoutwith
agree very well.
the buoyantly rising fluid elements (Blaes et al. 2011;
Around z = ±H , there is a big difference between
Jiang et al. 2014a). The total energy flux carried by the s
F and F , because there is a significant advection
photons is the lab frame radiation flux F , which is the r,z r,z0
r,z flux F , which is comparable to F . The nature of
sum of the diffusive radiation flux F and radiation adv r,z0
r,z0 the advection flux will be discussed in the next section.
enthalpy flux
The importance of the advection flux was first pointed
outbyHiroseetal.(2009b)andfurtherexploredinlater
F =F +v E + v·P | . (3)
r,z r,z0 z r r z work (Blaes et al. 2011; Jiang et al. 2013a). Significant
Here v is the fluid velocity and P is the radiation pres- advective transport is also observed in the global radia-
r
sure tensor. The net advection flux is tion MHD simulation of super-Eddington accretion disk
byJiangetal.(2014a). Theadvectionfluxdropsaround
F =v E . (4) ±2H as energy is increasingly transported by the dif-
adv z r s
8 Y.-F. Jiang et al.
fusive flux. Because the advection flux transports en- neticpressure. Thefluctuationscanbequantifiedbythe
ergy from the region near the midplane of the disk to standard deviations, while the anti-correlation between
the part closer to the surface in spite of the large opti- density and magnetic pressure can be quantified by the
cal depth, only the dissipation rate dF /dz follows the cross correlation coefficient
r,z
heatingrateq+ andtheverticaldistributionofdFr,z0/dz (cid:16) (cid:17)
ismuchbroader. Thepeakoftheheatingrateq+isoffset (cid:104)(ρ−ρ) B2−B2 (cid:105)
from the midplane and consistent with the peak of the σρ,B2 = σ σ , (7)
magneticpressureP showninthebottompanelofFig- ρ B2
B
ure8. ThePoyntingfluxalsopeaksatthesamelocation
where (cid:104)·(cid:105) means horizontal average of the quantity and
as q+ and drops quickly with height.
ρandB2 arethehorizontallyaveragedρandB2 ateach
height. The standard deviations of ρ and B2 are σ and
ρ
4.2.3. Nature of the Advection Flux
σ . The top panel of Figure 11 shows the vertical pro-
B2
Theadvectionfluxplaysanimportantroleforthesta- filesofthestandarddeviationsofρandPt ≡PB+Pr+Pg
bilityandstructureofthedisk. Particularlywiththeiron scaledwiththehorizontallyaveragedρandPt. Thetime
opacitybump,itissignificantlyenhancedwithrespectto averaged vertical profile of σρ,B2 is shown at the bottom
theelectronscatteringcaseforthereasonswewillexplain panel of Figure 11. It is clear that density fluctuations
now. Theadvectionfluxisassociatedwiththebuoyantly are much larger than the fluctuations of total pressure.
rising pattern of the butterfly diagram shown in Figure The relative fluctuation σρ/ρ can reach 30% at ±Hs.
2. Around z = ±H , the local optical depth per typical There is a strong anti-correlation between density and
s
scaleheightHs isτ =ρκtHs =3.77×104,whichislarger magnetic pressure fluctuations as σρ,B2 is negative for
thanthecriticalopticaldepthτ ≡c/c =5.69×103 de- |z|(cid:46)1.5H .
c g s
fined in Jiang et al. (2015). According to the criterion The total flux F we get from the simulation only
max
ofefficientradiationpressuredominatedconvectionstud- corresponds to ∼2%−3% Eddington accretion rate de-
iedbyJiangetal.(2015),thephotondiffusiontimescale fined based on the electron scattering opacity. But we
is much longer than any local dynamic time scale when already observe significant advection flux in this simula-
τ >τ sothatphotonscanbetrappedandadvectedwith tionwithironopacity. ForthecomparisonrunsESR20a,
c
the buoyantly rising fluid elements. Around ±2H , the ESR20b and ESR20c, advection flux is negligible every-
s
local optical depth τ drops to 7.61×102, which is much where. Fortheelectronscatteringdominatedsimulations
smaller than τ . The diffusion time becomes short and presentedinJiangetal.(2013a), theadvectionfluxonly
c
the photons cannot be trapped anymore. The photons became comparable to the local diffusive flux when the
are released and transported out by the diffusive flux. total radiation flux was more than 10% Eddington lu-
Thefluidelementsmusteventuallyfallback(onaverage) minosity. Therefore, it seems that we are seeing an en-
since we do not observe a net outward mass flux when hancement of the advection flux in the simulation with
average over dynamo cycles. However, the falling fluid the additional Fe opacity.
elements are (on average) colder compared with rising One reason for the enhanced advective flux is that the
elements because radiative energy is lost due to photon higher opacity requires it. The amount of energy that
diffusion near the surface. Thus, there is a net vertical can be transported by the diffusive flux is actually con-
advective energy flux. Since the advection flux is much strained by the hydrostatic equilibrium, because in the
less sensitive to the optical depth, the net radiation flux radiation pressure dominated regime F = cΩ2z/κ
r,z0 t
at the surface of the disk can becomes much larger for and dF /dz = cΩ2/κ . However, the local heating
r,z0 t
a given Pz,0 than the value ∼ cPz,0/τ0 predicted by the rate from the MRI turbulence q+ is not limited by the
diffusion equation. opacity. AsshowninthebottompanelofFigure9,q+ is
Blaes et al. (2011) provide a detailed study on why muchlargerthandF /dz. Inordertomaintainather-
r,z0
the fluid elements rise buoyantly. Buoyant motion is not mal equilibrium state, the excess energy must be carried
drivenbythestandardconvectiveinstabilitybecausethe out by advection. Notice that dF /dz is significantly
r,z0
entropy profile is stable even when radiation entropy is below the red dashed curved corresponding to cΩ2/κ
included. Instead,itisthenonlinearoutcomeoftheanti- es
in Figure 9 because κ > κ . Comparison with Fig-
correlation between density and magnetic field fluctua- t es
ure 2 of Blaes et al. (2011) shows that dF /dz must
tions in the MRI turbulence. In the radiation pressure r,z0
instead follow the cΩ2/κ curve in electron scattering
dominated regime, large density and magnetic field fluc- es
dominated simulations, at least for the inner few scale
tuations can be produced by MRI turbulence (Turner
heights. Since dF /dz, and therefore F , must be
et al. 2003; Jiang et al. 2013b). If we consider horizon- r,z0 r,z0
higher in this electron scattering dominated case, there
tal variations at fixed height, low density regions have
is less need for advective flux to transport energy.
larger magnetic pressure to maintain the pressure bal-
ItalsoseemstobethecasethatthesimulationswithFe
ance horizontally. Therefore, low density, highly magne-
opacityarelessstabletomagneticbuoyancyinstabilities.
tized regions rise buoyantly (and vice versa). We have
We can considier this by comparing the hydrodynamic
confirmedthepresenceofthiseffectintherunOPALR20.
andmagnetohydrodynamicBrunt-V¨ais¨al¨afrequenciesas
Figure10showsthehorizontalsliceatz =H ofdensity,
s done by Blaes et al. (2011). Following equations (27),
magneticpressure,radiationpressureandMaxwellstress
(28) and (29) in Blaes et al. (2011), the square of the
at 100 orbit. The relatively lower density regions have
hydrodynamic Brunt-V¨ais¨al¨a frequency is
stronger magnetic pressure, and also relatively lower ra-
diation pressure to maintain the horizontal pressure bal- (cid:18) (cid:19)
1 dln(P +P ) dlnρ
ance. TheMaxwellstressiswellcorrelatedwiththemag- N2 =a r g − , (8)
g Γ dz dz
1
9
Figure 10. Horizontalsliceofdensity,magneticpressure,radiationpressureandMaxwellstressatz=Hs fortherunOPALR20at100
orbit.
where Γ is the first adiabatic index defined in equation V¨ais¨al¨a frequencies based on the horizontally averaged
1
(15) of Jiang et al. (2015). Negative N2 will indicate vertical profiles at each snapshot and then time aver-
that it is unstable to convection. In the regime when age them. Figure 12 shows the averaged vertical profiles
the diffusion time scale is longer than the local dynamic of N2, N2 and N2 . Consistent with Blaes et al.
mag mag,r
timescale, themagneticBrunt-V¨ais¨al¨afrequencycanbe (2011), N2 ispositivesothatthediskisstabletohydro-
defined as dynamic convection. However, N2 first increases with
Nm2ag =N2+ agcv2A2 dldnzB, (9) himeiugmht∼wit0h.0in3Ω±20.a7rHousnadnd±t1h.e2n−dr1o.p5Hsasn.dTrehaischmesinaimmuinm-
t does not exist in the electron scattering case shown in
wherev istheAlfv´envelocityandc2 ≡Γ (P +P )/ρ. Blaes et al. (2011) and it is caused by the iron opacity
A t 1 r g
When the diffusion time scale is much shorter than the peakaround ±0.7H (see thebottom panel ofFigure 7).
s
local dynamic time scale so that diffusion is rapid, the Becausemagneticpressurealsopeaksnear±0.7H and
s
relevant magnetic Brunt-V¨ais¨al¨a frequency is decreases with height after the peak, it makes the mag-
netic Brunt-V¨ais¨al¨a frequency N2 to become negative
a v2 dlnB mag
N2 = g A . (10) around ±Hs as shown in Figure 12, although magnetic
mag,r c2+v2 dz pressureisstillmuchsmallerthantheradiationpressure
t A
there. This means the region that is unstable to Parker
UndulatoryParkerinstabilitywillhappenwhenN2 or
mag instability is much deeper compared with the electron
N2 becomenegativeintheappropriateregime. Time scattering opacity dominated case. Above ±1.5H , dif-
mag,r s
averaged vertical profiles of N2, N2 and N2 are fusion is rapid and the relevant magnetic Brunt-V¨ais¨al¨a
mag mag,r
shown in Figure 12. Here we first calculate the Brunt- frequency is N2 , which is also negative. This mag-
mag,r
10 Y.-F. Jiang et al.
5. DISCUSSIONSANDCONCLUSIONS
By including the iron opacity bump in 3D local shear-
ing box radiation MHD simulation of AGN disks, we
show that the radiation pressure dominated accretion
disks can survive many thermal time scales without
showing significant thermal runaways. In contrast, if we
just change the opacity to be electron scattering plus
free-free as in the standard α disk model, the disk col-
lapses quickly. The iron opacity bump can make the ra-
diation pressure dominated accretion disks more stable
becauseitcausesthetotalopticaldepthtoanti-correlate
with the midplane pressure, and it enhances the advec-
tive cooling. Both effects increases the cooling rate and
make the cooling rate change as fast as the heating rate.
Theseresults havea widerangeofimplications forAGN
observations.
5.1. Comparisons with the α disk model
Neitheroftwoeffectswehaveidentifiedareincludedin
thestandardthindiskmodel. Althoughweonlyhaveso-
lutions for one set of parameters here, we expect the ad-
vection flux will also increase with increasing total pres-
sure, because the amount of energy that can be trans-
Figure 11. Top: time averaged vertical profiles of the standard
ported by the diffusive flux is limited by the hydrostatic
idzeovnitaatliloynsavoefrdaegnedsitvyaρluaesndfotrottahleprruesnsuOrePAPLtRs2c0a.ledBwoittthomth:ethimore- equilibrium. The excess heating will be transported by
averagedverticalprofileofthecrosscorrelationbetweenρandB2. the advection flux to maintain thermal equilibrium.
Steady state vertical structure of a radiation pressure
dominateddiskisconstrainedbytheequationsofhydro-
static equilibrium
dP
r =−ρΩ2z, (11)
dz
and thermal equilibrium
αP H Ω=F . (12)
z,0 ρ max
Hereweadoptthesameαansatzasinthestandardthin
disk model that the vertical integrated stress is related
to the midplane pressure P 2 and density scale height
z,0
(cid:115)
1 P
H ≡ z,0, (13)
ρ Ω ρ
z,0
where ρ is the midplane density. From the diffusion
z,0
equation
dP ρκ F
r =− t rz,0, (14)
dz c
we can define the flux weighted optical depth τ¯ as
1 (cid:90) Lz/2
aFnigdumreag1n2e.tohTyidmreodayvneraamgiecdBvreurtnitc-aVla¨pirso¨afil¨alesfreoqfutehnechieysdrNod2,ynNamm2aigc τ¯≡ Fmax 0 ρκtFr,z0dz. (15)
andNm2ag,r fortherunOPALR20. Then the cooling rate per unit area Q− = F =
max
cP /τ¯ and P ≈P in the radiation pressure dom-
rz,0 rz,0 z,0
netic pressure dominated region is Parker unstable, as inated regime. Because F = 0 at the disk midplane
rz,0
pointed out by Blaes et al. (2011). Because photons and where maximum density is located and F has con-
max
gas are tightly coupled around ±Hs (τ >τc), the Parker tributions from both the diffusive and advective compo-
instability enhances the buoyancy and theadvection flux nents, τ¯ is smaller than τ . Advection flux also changes
0
as shown in Figure 9. Above ±1.5H because of rapid
s
diffusion, even though N2 and N2 are negative, they 2 Notethatαdefinedthiswayistheratiooftheverticallyinte-
mag,r
do not cause any significant advection flux as in the in- gratedstresstotheproductofmidplanepressureandthedensity
scale, which can differ from the ratio of the vertically integrated
efficient convection case studied by Jiang et al. (2015).
stresstotheverticallyintegratedpressureshowninFigure1.