Radiation-MHD simulations of blackhole accretion flows and outflows;~disk, outflow, and corona~

Ken OHSUGA (NAOJ)
Hiroyuki R. TAKAHASHI, Tomohisa KAWASHIMA (NAOJ),
Mariko NOMURA (Tsukuba Univ.), Shin Mineshige (Kyoto Univ.)

THREE MODELS
Abramowicz et al. 1995

ADAF/RIAF
Optically thin, Geometrically thick,and hot disk
LLAGN, Low-Hard state of BHB
Ichimaru 1977, Narayan, Yi 1994

Standard Disk
Optically thick, geometricallythick, cold disk
Seyfert, High-soft sate of BHB
Shakura & Sunyaev 1973

Slim Disk
Optically and geometricallythick, super-Eddington disk
NLS1? ULX?
Abramowicz et al. 1988

MAGNETO ROTATIONAL
INSTABILITY(c)Matsumoto(c)Hirose1. (a) illustrates the range of spectral shapes seen from the GBH. (b) compresses all the
spectral information into intrinsic colours for black holes (black), atoll systems (red points),
Z sources (cyan points) and the odd Z source Cir X–1 (blue points). The blue shaded area
indicates the section of the colour-colour diagram corresponding to the high/soft state spectra(Fig. 1a) which is only occupied by black holes.8)

DISK CORONA

Done & Gierlinski 2004

Disk spectra= soft-comp. + PL-comp.
PL-component is thought tobe produced by theCompton upscattering in
the disk corona.
However, disk-corona is not
yet well understood.

2. Sketched geometries corresponding to the spectral states as a function of Lh /Ls , together
with the eqpair spectra for GBH (blue). At high L/LEdd the identical accretion flow in atoll
systems has an additional component from an optically thick boundary layer (red).

OUTFLOW (UFO)+ $
4  >Tombesi et al. 2011

Blueshifted
absorption line

$4    KE9F 1  (i)40% of AGNs (ii)~0.1c (iii) Lkin~Ljet~0.1Lph
-    3    Observed UFO features imply that the powerful disk winds
!          
2-3
are ejected around the BHs (~10 Rs).


= &1  2 &M8*6@&5   1;&       =M5988598E DISK & OUTFLOW
図 6-2
Rad.
pressure
1/2 ページ = 左右 80 ミリ × 天地 70 ミリ(electron
scattering)
図 6-2
1/2 ページ = 左右 80 ミリ × 天地 70 ミリ

BH

ULX？

ブラックホール

図 6-2

L/LEdd

NLS1

slim
disk
スリム円盤

降
着
1/2
率 ページ = 左右 80 ミリ × 天地 70 ミリ

ブラックホール

スリム円盤

降
着
率
エディントン限界

Rad. pressure (line)

降着率ブラックホール
エディントン限界

スリム円盤
標準円盤

standard
disk
標準円盤
ライアフ

XRB(HS state)

Sy

QSO

エディントン限界

ライアフ
標準円盤

Mag. pressure
RIAF

ライアフ

XRB(LH state)

LLAGN

MBH

Shakura & Sunyaev 73; Ichimaru
77, Abramowics et al. 88;
Narayan & Yi 94

JET FROM RIAF
Kato et al. 2004

The magnetic field lines are
twisted and the toroidal
components are enhanced
around the black hole.
The jets are accelerated by
the magnetic pressure
force.

LINE DRIVEN WIND FROMSTANDARD DISK
Nomura, Ohsuga et al. submitted

The radiation force forspectral lines accelerates
the matter, driving the diskwind.
Line force is effective in thecase of AGNs. Because thedisk effectively emits UVphotons.
This model can explain theabsorption features ofUFOs (Nomura-san’s talk).

OUTFLOW FROM SUPER-EDD. DISK
Takahashi, Ohsuga in prep.

Radiation-dominated thick disk forms. Jets and disk wind are
powered by the radiation force for the electron scattering.

TODAY’S PLAN
•We will show that our simulations of the super-

Eddington flows can explain the basic properties of the ULXs.
(i)luminosity (ii)spectra(iii)time variation (iv)ULX bubbles•We also demonstrate that super-Eddington accretion

onto NSs is feasible. Pulsating ULX is super-Eddington source?

•We will introduce some important simulations (very high state,
BH spin, initial magnetic fields, end so on).

BASIC EQ. : RADIATION-MHD
+

t
( v)
+
t
e
+
t
E0
+
t
B
=
t

· ( v) = 0
·

BB
4

vv

· (ev) =

p

|B|2
p+
8

=

4
· v + 2 J2c

4+

B + c E0

· (E0 v) =· F0

v : P0 + 4

B

v

4
J
c

4
J=
c

B

B

RHD termsc

F0

MHD terms

c E0

BASIC EQ. : GR RADIATION-MHD
Continuity

Gaussʼs law

induction eq.

Energy-Mom.
of MHD
Energy-Mom.
of Radiation

We numerically solve above equations with using Kerr-Schildspace-time in order to avoid the singularity at the event horizon.

SUPER-EDDINGTON FLOW
flux force ~ dragRad.
Energydensity
~0.1c

~0.4c

Takahashi, Ohsuga 2014

The disk luminosity is a few
or several LEdd.
The radiation pressure
force drives jets and disk
winds from the radiationpressure-dominated disk.
The highly relativistic jets
cannot be reproduced,
because the radiation dragprevents acceleration.

03

02

20 02

20

20

10 01

10

10

0

01

R/r

APPARENT LUMINOSITY

0

10

20

51

41

31

21

40

0

50

10

20

30

40

0

50

10

20

30

r/RS

r/RS

r/RS

(g) Log(E0) [erg cm-3](h) Log[E0/(e+B2/8π)](i) Log(Tgas/Trad)Sre[ )0E(goL )g(

30

12

13

14

15

16

17

0

1

2

3

-1

•The radiative flux is mildly

80 08

80

Radiative
Flux

70 07

0

1

2

3

40

5

4

5

80

collimatedsince the disk
70
70 is opticallyandgeometricallythick.
60
60

60 06

•Thus, apparent luminosity is much

50 05

50

50

z/RS

Fig. 5.— Two-dimensional distribution
the gas temperature, (c) the plasma-β, (dcomponent, (f) the magnetic pitch, (g) th
energies, (i) and the ratio of the gas temp
and black arrows in panel (a) indicate thephotosphere, at which the optical thicknes
vector.
10 01

z/RS

z/RS

R/z

S

40
40
largerthantheEddington
here, at which the optical thickness measured from the upper boundary is unity. Theluminosity
arrow in panel (g) sho
k arrows in panel (a) indicate the velocity vectors whose magnitude exceed the escape velocity. The dashed l
exceptthe edge-on(e.g., over t
30 0of30radiationfor
3 the gas temperature to the
30 view(i) and the ratio
temperature.
All values
are time-averaged
ent, (f) the magnetic pitch, (g) the radiation energy, (h) the ratio of the radiation energy to the sum of t
22L
caseof of
Edd for ≲20° in thetemperature, (c)field, (e) the sam
20 the20 energies via the toroidal component
02 plasma-β, (d) the magnetic
20
5.— Two-dimensional distribution of the various quantities for Model A: (a) the density overlaid with theMdot~100LEdd/c2, Ldisk~3LEdd).
40 04

10

r/RS

10

r/RS

r/RS

r/RS

Super-Eddington
flows
can
10
20
30
40
50
10
20
30
40
50
03
02
01
10
20
30
40
5
r/R
R/r
r/R
explainr/Rthe large X-ray
10
10
10
5.—
noitTwo-dimensional
ubirtsid lanoisnem
distribution
id-owT —.5of2011.thegiF variousluminosity
quantities for Model
A:ULXs.
(a) the density overlaid with theof
Ohsuga,Mineshige
20
20 energies via the toroidal component
temperature,d( , β-amsalp(c)
e20
htthe
)c(plasma-β,,erutare pm(d)et thesag magneticeht
of field, (e) the sam
10

20

30

0
0

10
20

10

20

30

S

0
0

S

40

50

0
0

10

20

30

S

40

50

0
0

10

20

30

40

5

S

nt,
ht )(f)g( the,hctimagnetic
p citengam
pitch,eht (g)) f( ,thetnenradiation
opmoc energy, (h) the ratio of the radiation energy to the sum of t㻌ᬮ䛛䛟‽┦ᑐㄽ䝣䜯䞁䝛䝹 䞉 䜰

SPECTA
Without
ComptonWith only thermal Compton㻌ᬮ䛛䛟‽┦ᑐㄽⓗ㏿ᗘ䛾䝣䜯䞁䝛䝹 䞉 䜰䜴䝖䝣䝻䞊㻌ᬮ䛛䛟‽┦ᑐㄽⓗ㏿ᗘ䛾䝣䜯䞁䝛䝹 䞉 䜰䜴䝖䝣䝻䞊jet㻌ᬮ䛛䛟‽┦ᑐㄽⓗ㏿ᗘ䛾㻌෭䛯䛟㧗ᐦᗘ䛺䝣䜯䞁䝛䝹 䞉 䜰䜴䝖䝣䝻䞊㻌෭䛯䛟㧗ᐦᗘ䛺ప㏿䜰䜴䝖䝣䝻䞊

ప㏿䜰䜴䝖䝣䝻䞊

wind㻌෭䛯䛟㧗ᐦᗘ䛺ప㏿䜰䜴䝖䝣䝻䞊

㻌෭䛯䛟㧗ᐦᗘ䛺ప㏿䜰䜴䝖䝣䝻䞊

⾪ᧁἼຍ⇕㡿ᇦ BH
㻌㍽ᑕᅽඃໃ෇┙⾪ᧁἼຍ⇕㡿ᇦ

㻌㍽ᑕᅽඃໃ෇┙⾪ᧁἼຍ⇕㡿ᇦ
㻌㍽ᑕᅽඃໃ෇┙㻌㍽ᑕᅽඃໃ෇┙Trapping
photons

BH
BH

disk

BH

Thermal &Bulk Compton

Hard X-ray photons are generated by thermal and bulk comptonization
2of Comptonization on the(Kawashima-san’s
SEDs for a model
with
Ṁ
=
200L
/c
. T
E
Talk).

ys spectrum for a face-on observer (0◦ −10◦ ) without Compton scatter0 (red),ackbodye showsviewingor Ṁ ≈

Kawashima et al. 2012
Kawashima et al.
(data; Gladstone 2009)

Emergent
spectrais(12pp),sensitive
The AstrophysicalJournal, 752:18
2012 June 10 tothe mass
accretion
rate is higher,angle.
our calculations inthethisWhen
MBHandinclination
paper show that SEDs become softer, although our previous

NGC1313 X-2

Kawashima et al.

paper (Kawashima et al. 2009) expected that SEDs become
harder because of Compton upscattering by thermal electrons
in the jet. Kawashima et al. (2009) evaluated the photon indicesby calculating the Compton y-parameters in the region τeff ! 1.
The effective optical depth is calculated by integrating thefree–free absorption and the scattering coefficients from theouter boundary at r = 500rs to the inner boundary at r = 3rs
along the rays with polar angle θ = constant, for various θ .
(Here, we note that the computational domain for radiationhydrodynamic simulation in Kawashima et al. 2009 was smaller
than that in this work.) Since Kawashima et al. (2009) did
not take into account the downscattering by cool outflows,
they overestimated the spectral hardening when the accretion
rate is high. When the mass accretion rate is higher, a greaternumber of escaping photons tends to enter into the cool outflow
surrounding the jet and they are subsequently downscattered,because the funnel becomes narrower and jet becomes thicker
for Thomson scatterings. This is the reason why the SEDs
become softer when the mass accretion rate is higher, in
this work. Nevertheless, the conclusion in Kawashima et al.
(2009) that the spectral state changes from the slim disk one
to the Comptonizing outflows is valid in highly supercriticalaccretion flows because the SED of the photons escaping fromsuch highly supercritical flows is determined by the Compton
downscatterings in the outflow.

1

0.1

IC342 X-1
1

0.1

NGC5204 X-1
1

Figure 4. Dependence of SEDs on viewing angle i for the models with Ṁ ≈1000LE /c2 (upper panel) and Ṁ ≈ 200LE /c2 (lower panel).

4.3. Comparison with Observations of ULXs

Simulated spectra nicely
fit the observations.

Recent X-ray observatories (Suzaku, XMM-Newton, andChandra) provide us with precise spectra of ULXs below
10 keV. In Figure 7, we compare the SEDs obtained in this workand the XMM-Newton pn data of three ULXs: NGC 1313 X-2,
IC 342 X-1, and NGC 5204 X-1, which are kindly provided by
J. C. Gladstone. The SEDs of our results are normalized by the
X-ray luminosity shown in Figure 3, while the SEDs of ULXs

0.1

0.3

1

10

Figure 7. Comparison of SEDs obtained from our radiation hydrodynamic

CLUMPY OUTFLOW

Jet

Super-Eddington disk+ Jet

Takeuchi, Ohsuga, Mineshige 2013

accretion

Our simulations succeededin reporoducing theclumpy outflow.

ClumpyOutﬂow

Time-dependent, Clumpy
outflow with wide angle

ULX BUBBLE
Hashizume, Ohsuga et al. 2015

wide-angle

escvvr> sc
e
v
<
vr

BH

At the direction of 45°,
we find vr>vesc near theblack hole.
In 𝜃~45°-80°, outflow
velocity gradually
increases and exceeds
vesc at r ~1000-4000Rs.
Such a wide angle
outflow is accelerated byradiation force.

1.0
0.1

0°

Wide angle outflow
would produce theshock excited bubbles
around ULX.

Kinetic Power (vr>vesc)
Poutflow~3x1039erg/s(10Msun)

angle

90°

Kinetic power isconsistent withobservations(3-5x1039erg/s).

Super-Eddington model
would explain the shock
excited bubbles.

Hα + [ N II ]Cha

SUPER SOFT SOURCES & ULXULX Soft ULX ?

SSS ?

Surface of τeff=1.0 islocated around 103Rs.

The super-Eddington flows
exhibit the very soft spectra
for edge-on observer. Thus,
super soft sources might be
the off-set ULXs.

SUPER-EDDINGTONACCRETION ONTO NS
Publisher: NPG; Journal: Nature: Nature; Article Type: Physics letter

FIGURES

Bachetti et al. 2014

FPMA
FPMB

011

008-009

007

006

Is the central object NS ?
004

0.15

ObsID 002

0.2

0.1
Period (ms) - 1372.5

Pulsating ULX is discovered !

Chandra0.25

Is super-Eddington accretion
onto NSs feasible ?

1 b
0.5
0
−0.5
−1
30

Pulsed flux (%)

25
20

c

5
0

Δχ

Couts / s / module

0.3 a

−5
0

1
Phase

2

15
10

03 - 10 keV
03 - 30 keV
10 - 30 keV

5
0
56680

56690

56700
MJD

56710

56720

K. Ohsuga

Shell-shaped
structureappears.

NS

Ohsuga 2007

High conversion
eﬃciency (NS)[Vol. 59,

BH

Downloaded from http://pasj.oxford

ate onto the NS/BH and mass-output rate bynormalized by the critical rate, mP acc (circles)
unctions of the normalized mass-input rate.

Lkin > Lph (NS)Lkin < Lph (BH)
Downloaded from http://pasj.oxfordjournals.org/ by guest on November 20, 2014

K. Ohsuga[Vol. 59,

Fig. 4. Luminosity (circles) and kinetic energy output rate (squares)
as functions of normalized mass-accretion rate onto the central object,
MP acc c 2 =LE , for MP input c 2 =LE = 3 ! 102 ; 103 ; 3 ! 103 from left toONGOING SIMULATIONS
super-Eddington disk accretion
onto weakly magnetized NS(Takahashi & Ohsuga)super-Eddington pole accretion
onto strongly magnetized NS(Kawashima, Ohsuga, Mineshige)densityRad.
EnergydensityNSNSDISK TRUNCATION & DISK CORONA
Takahashi, Ohsuga, 2015-100

10 [Rs]Gas
Temperature>109K

In our GR radiation-MHD
simulations, standard disk is
truncated at ~ several Rs.
Gas temperature is very
high ,>109K, within thetruncation radius (disk corona).

~107K
MBH=10Msun
Mdot~LEdd/c2
spin parameter ~ 0.9

Disk structure for each mode
here and dashed curves indicat
DISK TRUNCATION & DISK CORONA
Overheated region sandwiches
the standard disk.
Compton upscattering in thecorona region would play animportant role for producing thehard X-ray photons.

Gas Temperature(r-z plane)

30

20

10

0

-10

-20

03

02

01

0

01-

02-

Very High State ?
Corona region(>109K)

Standard disk region(~107K)

edom hcae rof erutcurts ksiD

– 19 –HEATING OF CORONA & BH SPIN
1012
1011
1010

Gas Temperature(equatorial plane)
a0r4
Rotating BH (a=0.9)
a0r4-2
Non-rotating BH
a0r2
(a=0)

Gas temperature forrotating BH is about
10 time higher thanthat for non-rotating
BH.

a9r4

109
108
107
106 0

5

10

15

20

25

Fig. 5.— Radial profiles of < Tgas >ρ,θφ (solid) and < Trad >ρ,θφ (dashed).

30*Gas temperature might
decrease via Compton
cooling.

HEATING OF CORONA & BH SPIN
– 21 –

Outward Poynting Flux
a0r4

10

a9r4

-4.0

5

Outward Poynting Flux isenhanced via the rotationof BH.

-5.0

The matter around therotating BH is more
effectively heated up by the-6.0
dissipation of thetransported electromagnetic energy.

0-5

Non rotating BH

-10
-10

Rotating BH

-7.0

-5

0

5

10

ACCRETION, JET, & BH SPIN
Mass accretion rate islarger in retrograde
BHs than in prograde
BHs.
However, outflow
efficiency is higher inprograde BHs than in
retrograde BHs.

Time evolution of accretion rate and outflow efficiency

accretion rateaccretion rate

outflow efficiency
outflow efficiency& column).
McKinney
Figure 1. Initial conditions and time-evolution of our fiducial models, ATchekhovskoy
0.9f (left column) and A0.9f (rightSee Supporting2012
Information for movie[Panels a,f] A vertical slice through the initial conditions. Colour shows fluid frame rest-mass energy density (see colour bar) and thin black lines show leveof constancy of enclosed magnetic flux, h i, which represent field lines in the image plane. Magnetic field is axisymmetric, with B' = 0 everywhere anSimulations of Magnetized Advection Dominated Accretion

5

Figure 1. Poloidal plane of the grid used in the simulations, shown at two zoom levels.

INITIAL MAGNETIC FIELDS & OUTFLOW
B: 8 loops

B: single loop

Structure of themagnetic fields affectthe flow !
Figure 2. Initial configuration of the ADAF/SANE simulation. The top two panels show
the mid-plane
flux threading
the equatorial
Figure3. Similar density
to Fig. 2andbuttheformagnetic
the ADAF/MAD
simulation.
The main difference is that here the torus has a single loop of field centered at radius r = 300.
plane as a function of radius. Note the extended size of the initial torus, which is requiredfortheextremely
longdurationof
thissimulation.
also the around the BH, leading to the MAD state.
As a result, accretion causes magnetic flux of one sign toNote
accumulate
multiple oscillations in the magnetic flux, which prevents the accreting gas from reaching the MAD state. The lower two panels show the logarithm of thedensity ρ and the gas-to-magnetic pressure ratio β of the initial torus in the poloidal plane.

6Single loop B (right)→ Outflowing Motion

R. Narayan, A. Sadowski,
˛
R. F. Penna, A. K. Kulkarni

Multi loop B (left)→ Convective Motion

over the course of the simulation. The normalization of the magnetic field is adjusted such that the gas-to-magnetic pressure ratio,β, in the equatorial plane has a minimum value ∼ 100 for each of
the eight loops. Instead of using multiple poloidal loops, anotherway of setting up an ADAF/SANE simulation is to use a toroidal
initial field (e.g., Model A in Igumenshchev et al. 2003 and Model
A0.0BtN10 in McKinney et al. 2012).
The initial magnetic field of the ADAF/MAD run forms a single poloidal loop centered at r = 300 (Fig. 3). The gas accreted by
the BH in this simulation has the same orientation of the poloidal
magnetic field throughout the run, so the accretion flow is maintained in the MAD state. The minimum value of β in the initial
torus is ∼ 50.
The magnetic field construction is described in detail in
Penna et al. (2012).3
2.3

Preliminary Discussion of the Simulations

The two panels in Fig. 4 show snapshots from the end of the
ADAF/SANE and ADAF/MAD simulation. In each panel, theblack and white streaks and red arrows show velocity streamlines

Convective
motion3

Narayan et al. 2012

In the notation of Penna et al. (2012), the ADAF/SANE magnetic field
has rstart = 25M , rend = 550M , and λB = 2.5. The ADAF/MAD
magnetic field has rstart = 25M , rend = 810M , and λB = 25.

in the poloidal plane at azimuthal angle φ = 0, and the dashedlines correspond to one density scale height. The main difference
between the two simulations is that the SANE run exhibits moreturbulence compared to the MAD run.
Following Penna et al. (2010), we define the mass accretion
rate Ṁ , the accreted specific energy e, and the accreted specificangular momentum j, at radius r and time t, as follows:! !
Ṁ (r, t) = −ρur dAθφ ,
(1)
θ

φ

! !
Ė(r, t)1e(r, t) ==Ttr dAθφ ,(2)
Ṁ (r, t)Ṁ (r, t) θ φ
! !
˙ t)
J(r,
1
j(r, t) =Tφr dAθφ ,(3)=−Ṁ (r, t)Ṁ (r, t) θ φ√
where dAθφ = −gdθdφ is an area element in the θ-φ plane,ρ is rest mass density, uµ is the four-velocity, and Ttr and Tφr arecomponents of the stress-energy tensor describing the radial flux ofenergy and angular momentum, respectively,

Outflow

Ttr

=

Tφr

=

(ρ + Γu + b2 )ur ut − br bt ,(ρ + Γu + b2 )ur uφ − br bφ .

(4)
(5)

The quantity u is the internal energy of the gas, Γ is its adiabatic
index which is set to 5/3 in both simulations, and bµ is a four-vector
which describes the fluid frame magnetic field (see Gammie et al.

GR MHD SIM. OF TILTED DISK(c) C. Fragile

rotation axis (BH)

rotation axis (BH)

Most of simulations are assumed the rotation axis of
BH is aligned with that of the disk.
Simulations of tilted disk performed by Fragile exhibit
strong time variation !

SUMMARY•Our simulations of the super-Eddington flows can explain thebasic features of the ULXs.

(i)luminosity (ii)spectra (iii)time variation (iv)ULX bubbles•Calculations of the observed spectra (including metal lines) are
important feature work.

-Our hydrodynamics simulations + radiation transfer
calculations (Odaka-san’s code)-Multi-frequency radiation-MHD simulations would beattempted in super Computer K or post-K computer.

•Accretion flows onto BHs/NSs and launching outflows are important
targets for ASTRO-H