Algorithms and Ray Tracing For Tori Around Black Holes: Leela Elpida Koutsantoniou

Astronomy & Astrophysics manuscript no.
output c
ESO 2020
August 31, 2020
Algorithms and ray tracing for tori around black holes

Leela Elpida Koutsantoniou1,2,?
1
Research Center for Astronomy and Applied Mathematics, Academy of Athens, Athens 11527, Greece
2
Department of Astrophysics, Astronomy and Mechanics, Faculty of Physics, University of Athens
Received ? ?, ?; accepted ? ?, ?
ABSTRACT
We present our study and results about accretion disks around stellar black holes. We build and solve the relativistic radiative transfer
equations and examine particle trajectories in the environment or moving with the disk. These particles have various velocity profiles
and absorb thermal radiation from the hot disk itself. We consider six different disk models, four opaque and two semi-opaque and
various velocity profiles, depending on the target’s location in each system. We then describe and explain the five families of codes
written for our newer expanded work and present various images and results for forces, black hole images, spin estimation and particle
trajectories that also include radiation contributions.
Key words. accretion, accretion disks – black hole physics – radiative transfer
1. Introduction Miller 1995. All these works consider a central source of pho-
tons, that is in some cases rotating. In more recent years, it is
Black holes have been for the past years, one of the most talked worth mentioning many years of work by Bini et al. (2009, 2011,
about subjects of Astrophysics, along with gravitational waves. 2015), who have given interesting results.
Pictures of these objects, their environments and their progeni-
Our work initially started by studying Contopoulos &
tors have been studied and discussed in length and efforts were
Kazanas 1998 and considering the Cosmic Battery model. This
made to further study these far away phenomena. Something
study states that radiation is not only causing various effects
however, that is often neglected and sidelined, could prove to be
on matter, but is in cases even enough to generate magnetic
of great importance to all these objects and their surroundings:
fields from zero, due to the Poynting–Robertson drag. This work
radiation.
was later on revisited and extended in Contopoulos, Kazanas,
Radiation, along with its pressure and effects appear negli-
& Christodoulou 2006 and Contopoulos & Papadopoulos 2012.
gible in settings so dense and energetic as black holes, but one
Further on, it was used to examine the effects of radiation onto
should keep in mind that objects like these, are accompanied by
X-ray binaries and the evolution of the hardness-intensity dia-
accretion disks. These accretion disks are usually of a tempera-
gram in Kylafis et al. 2012. Our work continued on to study the
ture high enough to produce X-rays and large numbers of pho-
effects of radiation in the much more complicated environments
tons of this energy cannot be deemed unimportant or negligible
of hot radiating accretion disks orbiting BH in Koutsantoniou
at any case.
(2014) and then Koutsantoniou & Contopoulos 2014. We also
Noticeable effects of the Solar radiation, far inferior than a
examined the possibility and the capability of the radiation cre-
black hole’s X-ray radiation, were first described by J. H. Poynt-
ating magnetic fields from zero.
ing (Poynting 1903). Years later, H. P. Robertson (Robertson
1937) properly explained what is now known as the Poynting– Last but not least, we studied the work of Sadowski
˛ et al.
Robertson effect and cleared misconceptions that had puzzled 2016 and Sadowski
˛ 2016 that examine a variety of accretion
many great scientists, such as J. Larmor for years. This effect is disks and the effects of radiation and magnetic fields, along with
also called a drag because, even though the Solar photons are their feedbacks to these systems.
emitted radially outwards, they cause the dust grains in orbit ab- The current paper is a continuation of our previous works
sorbing them, to slowly brake and infall onto the Sun, the Earth but involves and also describes a great leap forward in both the
or any other close enough massive object. This force is a purely codes quality and effectiveness. The new integrations allow us to
relativistic effect and is proportional to the absorbing object’s az- use several different machines running codes and automatically
imuthal velocity. A simple non relativistic formula can calculate saving and sharing the results. Including a new process, there is
this force to good accuracy as the possibility to increase the codes resolution by 25 − 625 times
or in cases even more, by only doubling the execution time and
φ Vφ LσT Vφ that in only certain cases. This allows us to now run a vastly in-
fPR = − frad
r
=− , (1)
c 4πr2 c2 creased number of simulations compared to our previous work
and instead of having a sole target at each disk’s innermost sta-
where Vφ is the target’s azimuthal velocity.
ble circular orbit, hereafter ISCO, we spread absorbing targets
Moving towards our initial goal of black hole radiation, we
throughout the various systems, inside the disk, outside and if
must first look at important early work such as Abramowicz,
possible throughout it. We also greatly extend by leaving the
Ellis, & Lanza 1990 and Miller & Lamb 1993, 1996, Lamb &
equatorial plane and spreading the targets to all the sphere’s vol-
?
email:[email protected] ume from close to the vertical z-axis, all the way down to the
Article number, page 1 of 40
A&A proofs: manuscript no. output
equator. Finally, it is worth mentioning that in contrast to the rotating one. The second characteristic surface is the static limit
vast majority of BH studies, our fully ray tracing portion of our that constitutes the outer boundary surface of the ergosphere and
codes allows us to look very close to these objects. Even though can be found at the point where the gtt component changes sign:
this possibility makes ray tracing codes rather demanding, it al- √
lows us to get images of the BH and the accretion disk for radii rergo = M + M 2 − a2 cos2 θ. (7)
ranging from to 1.25M for rotating BH, up to 25M or more.
In this work, we start in §2 by providing the mathematical Finally, a noteworthy set of radii are the equatorial circular orbits
formulation necessary in order to set up and use the Kerr metric for massive particles and in particular the ISCO, that is given by
and the locally nonrotating frames and then the methods to study
particle trajectories and radiation effects. In §3 we present the h p i
various models of disks used in our work and the accompanying rIS CO = M 3 + Z2 ∓ (3 − Z1 ) (3 + Z1 + 2Z2 ) , (8)
physics. In 4, we present and describe the five different families
of codes written for our studies and show many new results from where
these studies. Finally, in 5 we summarize our conclusions about !1/3 "
a2
#
the results and discuss the possibilities of future extensions. a 1/3 a 1/3
Z1 = 1 + 1 − 2 1+ + 1−
M M M
2. Mathematical Formulation r
3a2
We assume that the immediate environment around a rotating Z2 = + Z12 . (9)
and accreting black hole, hereafter BH, can be adequately well M2
described using the Kerr metric. This suggests that the space- and the upper signs refer to direct and the lower signs to retro-
time and its form are determined by the central compact object grade orbits. The ISCO starts from a value of rIS CO = 6M for
that is axisymmetric, uncharged and possibly rotating. We also a = 0 and for a = M, reaches rIS CO = M for a direct orbit.
assume that the test particle’s presence and motion do not affect
the spacetime form or the stress-energy tensor. We will hereafter
use the geometrized unit system in which c = G = 1 and assume 2.2. Locally nonrotating frames
the Einstein notation for summation over double indices. Space-
time components will be denoted by Greek indices and space The Kerr spacetime is stationary and axisymmetric but the cen-
components by Latin indices. tral object rotation introduces complexity in both the physics of
the problem and the mathematics required. First of all, the non-
diagonality of the metric introduces cumbersome algebraic cal-
2.1. The Kerr Metric culations when rising or lowering indices. In addition, physical
difficulties arise when examining locations within the static limit
The black hole and the spacetime it creates, can be fully de- and throughout the ergosphere. This is due to the fact that, there
scribed using its mass M and its spin parameter a. The Kerr met- cannot be static Boyer – Lindquist (hereafter BL) observers at
ric in Boyer – Lindquist (t, φ, r, θ) coordinates is then given by these points, since the t basis vector becomes spacelike.
ds2 = gαβ dxα dxβ = In order to simplify the calculations and remove the prob-
lems in the physics at points inside the ergosphere, we choose
= −e2ν dt2 + e2ψ (dφ − ωdt)2 + e2µ1 dr2 + e2µ2 dθ2 , (2) to introduce a new set or observers and work in that new frame.
where The best choice of observers is one where said observers rotate
with the spacetime geometry at the point we wish to study. We
Σ∆ 2ψ Asin2 θ 2µ1 Σ thus define the locally nonrotating frame (LNRF) or the zero
e2ν = , e = , e = , e2µ2 = Σ , (3)
A Σ ∆ angular momentum observer (ZAMO) at the point in question
with and describe the requested quantities using their projection on
the Minkowskian orthonormal frame of the local observer. If re-
∆ = r2 − 2Mr + a2 , quired, we can then easily switch the calculated quantities from
the LNRF into the BL frame. Quantities calculated in the LNRF
Σ = r2 + a2 cos2 θ , are denoted by using hats over the component indices, e.g. uα̂ ,
2 and quantities calculated in the BL frame are denoted by unhat-
A = r2 + a2 − a2 ∆sin2 θ , (4)
ted indices, e.g. uα . The transformation tensor Eαµ̂ and Eµ̂α be-
and the spacetime angular velocity given by tween the two frames has nonzero components
gφ t 2Mra
ω=− = , (5) ettˆ = eν , eφ̂t = −ωeψ , eφ̂φ = eψ , er̂r = eµ1 , eθ̂θ = eµ2 ,
gφ φ A
see Bardeen (1970) and Bardeen, Press, & Teukolsky (1972).
From the metric (2), one can determine the various charac- ettˆ = e−ν , eφtˆ = ωe−ν , eφφ̂ = e−ψ , err̂ = e−µ1 , eθθ̂ = e−µ2 . (10)
teristic surfaces present around a rotating black hole. The event
horizon arises from the pole of the grr component and is found Vectors uα and uµ̂ , and tensors T αβ and T µ̂ν̂ are transformed fol-
at the outermost root of the equation ∆ = 0: lowing the equations
√
revh = M + M 2 − a2 . (6) ua = eαµ̂ uµ̂ , uµ̂ = eµ̂α uα
The event horizon is thus a sphere of radius revh = 2M for a non-
rotating Schwarzschild black hole and revh = M for a maximally T αβ = eαµ̂ eβν̂ T µ̂ν̂ , T µ̂ν̂ = eµ̂α eν̂β T αβ . (11)
Leela Elpida Koutsantoniou : Algorithms and ray tracing for tori around black holes
2.3. Particle trajectories d pφ

= 0,
dλ
In order to study particle trajectories in Kerr spacetime, it is nec-
∂ ∆ ∂ 1 ∂ Vr + ∆Vθ
! ! !
essary to make full use of the integrals of the particle motion d pr
= −pr 2 − pθ 2 + ,
and the respective conserved quantities. Let us assume a parti- dλ ∂r 2Σ ∂r 2Σ ∂r 2Σ∆
cle with rest mass µ and four-momentum p = pt , pφ , pr , pθ in
∂ ∆ ∂ 1 ∂ Vr + ∆Vθ
! ! !
d pθ
geodesic motion around a rotating uncharged black hole. This = −pr 2 − pθ 2 + . (14)
particle has got four conserved quantities: the particle rest mass dλ ∂θ 2Σ ∂θ 2Σ ∂θ 2Σ∆
µ, the total energy E = −pt , the angular momentum component We have thus transformed the initial four equations of motion
parallel to the rotation and symmetry
axis L = pφ and the Carter into a new system of eight differential equations. The new forms
p2
constant Q = pθ + cos θ a µ − pt 2 + sinφ2 θ (Carter 1968).
2 2 2 2 are smooth and do not have poles or other problems throughout
their range and can be directly integrated. Note that the fifth and
This constant could perhaps be simply explained as a measure sixth of the above set of equations, describe two of the motion’s
of how much a trajectory deviates from the equatorial plane. A conserved quantities, the conservation of energy and angular mo-
particle in geodesic motion that starts in the equatorial plane and mentum respectively.
has Q = 0 will remain there indefinitely and a particle moving Finally, we also define here the coordinate angular velocity
outside the equatorial plane with Q > 0 will at some point cross for a circular equatorial orbit as
it.
The equations describing the particle motion are (Bardeen dφ uφ M 1/2
et al. 1972; Wilkins 1972): Ω= = t = 3/2 (15)
dt u r + aM 1/2
dt T
Σ = −a aEsin2 θ − L + r2 + a2 , where u = ut , uφ , ur , uθ is the four-velocity of a particle.

dλ ∆
!
dφ L T
Σ = − aE − +a , 2.4. Radiation and equations of motion
dλ sin2 θ ∆
Here in this paragraph, we discuss the quantities we wish to cal-
dr p
Σ = ± Vr , culate and seek the appropriate series of equations that can allow
dλ us to estimate them.
dθ Our main goal is to calculate the effects of radiation on the
Σ = ± Vθ ,
p
(12) target particle dynamics. Thus, we begin from the formula that
dλ
relates the target particle position with the acceleration
where the effective potentials are given by
d 2 xα dxµ dxν
+ Γαµν = aα ,

T = r2 + a2 E − aL, 2
(16)
dτ dτ dτ
where xα are the particle position components, τ the proper
h i
Vr = T 2 − ∆ µ2 r2 + (L − aE)2 + Q ,
time and Γαµν the Christoffel symbols or connection coefficients
L2
" #
(Mueller & Grave 2009).
Vθ = Q − cos2 θ a2 µ2 − E 2 + , (13) The acceleration aα in turn, can be given by the relativistic
sin2 θ
equation of motion
and λ is an affine parameter for massless particles and λ = τ/µ
for massive particles, with τ the particle’s proper time. fa
The above form of the equations is compact and elegant but aα = , (17)
m
hides various problems that appear when one attempts to solve
them. The system appears to be highly problematic during nu- where m is the rest mass of the target particle and f α are the 4-
merical integration, since the presence of the square roots in the force components. Notice that since this is the relativistic equa-
latter two equations cause the quick accumulation of errors near tion of motion, this f α only includes nongravitational forces.
the turning points. There are various solutions, such as reparame- Since we wish to investigate the effects of radiation on the par-
terization, in order to deal with this issue. In our work, we choose ticle motion, this 4-force is for us specifically the radiation 4-
to work with the Hamiltonian and we transform the above system force.
accordingly. The new system of equations is then as follows: In order to calculate the 4-force f α , one needs to know the
flux of the radiation that generates it, as
dt 1 ∂
=− (Vr + ∆Vθ ) , f α = σF α , (18)
dλ 2Σ∆ ∂pt
dφ 1 ∂ where σ is the particle cross section for the momentum transfer.
=− (Vr + ∆Vθ ) , The radiation flux 4-vector F α can in turn be calculated us-
dλ 2Σ∆ ∂pφ
ing the target particle covariant 4-velocity uµ and the radiation
dr ∆ stress–energy tensor T αβ using the formula
= pr ,
dλ Σ
F α = hαβ T µβ uµ , (19)
dθ 1
= pθ ,
dλ Σ where hαβ is the projection tensor
d pt
= 0, haβ = δαβ + uα uβ . (20)
dλ
2.5. Radiative transfer

In this paragraph we look into disks with finite optical depth,
where photons are emitted by the material throughout its entire
volume. The radiation is then regulated by its passage through
the disk’s absorbing and emitting material, then escapes the disk
and travels throughout the entire system. We thus attempt here to
find a way to calculate the specific intensity Iν of this radiation.
For this reason we look into the radiative transfer equation (see
Rybicki & Lightman 1986), hereafter RTE, and the necessary
changes required to obtain it in a Lorentz invariant form.
We begin by assuming we have a thermalized material of
number density n. This material consists of particles that will
act as radiation
absorbers. We define the absorption coefficient
aν cm at frequency ν as
−1
Fig. 1: Local sky around the target particle. For the incoming
photon, we define angles ã and b̃ similar to the typical polar an- aν = nσν , (25)
gle θ and azimuthal angle φ of spherical coordinate systems. The
disk in stripes at the top left represents the black hole event hori- where σν cm2 is the absorbing area cross section at
zon. said
frequency. Assuming an initial specific intensity
Iν ergcm−2 s−1 ster−1 Hz−1 at frequency ν, the presence of
In order then to acquire the BL radiation stress–energy tensor the material’s radiation absorbing particles for a propagation
T αβ , it is in general necessary to find the LNRF stress–energy length ds, will cause a decrease in this specific intensity of a
tensor T µ̂ν̂ and then use eq. (11): propagating light ray given by
T αβ = eαµ̂ eβν̂ T µ̂ν̂ , (21) dIν = −aν Iν ds. (26)
eαµ̂
where the are given by eq. (10). Things are simpler for the emission coefficient
In order to calculate now T µ̂ν̂ , we make use of the formula

jν ergcm−3 s−1 ster−1 Hz−1 . When the light ray propagates
for distance ds, it transverses emitting material of volume dV
Z
T µ̂ν̂ = I r, θ, ã, b̃ nµ̂ nν̂ dΩ̃,

(22) and has its specific intensity increased as
where dΩ̃ = sin ã dã db̃ is the solid angle element with ã and dIν = jν ds. (27)
µ̂ µ̂
appropriate local angles (fig. 1), n = p /p and
tˆ
b̃ the related
The radiative transfer equation combines the above two pro-
I r, θ, ã, b̃ the frequency integrated specific intensity of the ra-
cesses and describes the resulting effects on the light ray’s spe-
diation. Simple calculations can give the nµ̂ vector components cific intensity as
as
dIν
nφ̂ = sin ã sin b̃, nr̂ = cos ã, nθ̂ = sin ã cos b̃. (23) = −aν Iν + jν . (28)
ds
Intensity I is, as expected, considered to be a function of the
particle’s position in space, since different positions receive dif- In order to express the solution of the above equation more el-
ferent amounts of radiation. From this, we have already excluded egantly, we introduce the concept of the optical depth τν at fre-
the φ coordinate due to the spacetime axisymmetry. In addition, quency ν that is defined as
I also naturally depends upon the angles ã and b̃, since different dτν = aν ds. (29)
amounts of radiation are received when looking towards differ-
ent parts of the local sky. The optical depth can then be calculated by integrating the above
Finally, the frequency integrated specific intensity I can, as along the path of the light ray
its name suggests, be calculated by integrating the specific inten- Z s
sity Iν across all the contributing frequencies
τν (s) = aν s0 ds0 ,

Z (30)
s0
I = Iν dν. (24)
ν where s0 is the arbitrarily selected initial point of the scale. The
This can be applied for anything, such as a blackbody dis- radiative transfer equation can then be restated as
tribution, thermal radiation, a single-energy beam of light or
dIν jν
frequency-independent radiation. In an environment where the = −Iν + . (31)
radiation is emitted by just a surface layer of the source, the cal- dτν aν
culation of I is relatively easy and it is done in the method de-
Integrating this gives the solution to the radiative transfer equa-
scribed in §2.6. However, in environments where the radiation is
tion
emitted by multiple layers or various objects, the method of said Z s
paragraph cannot be used. In order to calculate there the spe- 0
Iν (s) = Iν (s0 ) e +
−τν
jν s0 e−[τν (s)−τν (s )] ds0 .

cific intensity Iν , it is required to investigate the radiative transfer (32)
s0
processes and solve the appropriate equation. This is generally a
rather lengthy subject, both in theory and especially in applica- The above magnitudes are not in most cases Lorentz invari-
tions and will be separately addressed in the following §2.5. ant and thus cannot be used in the general solution of the various
and the matter rest frame, since it is perpendicular to the direc-

tion of motion. Likewise, the x-component of the photon mo-
mentum k x = k x 0 , remains unchanged. This however, means that
k sin θ = k0 sin θ0 and thus ν sin θ = ν0 sin θ0 . From the Lorentz
invariance of the optical depth τν = aν s, we then have
aν = νaν = ν0 aν 0 = Lorentz invariant. (42)
Finally, for the emission coefficient we utilize eqs. (31), (40)
and (42) and conclude that
jν jν 0
jν = = 0 2 = Lorentz invariant. (43)
ν 2 ν
Fig. 2: (left) Emitting matter in the lab frame K, moving with The Lorentz invariant form of the RTE (31) will thus be
velocity →−u along the y axis in a tube of width d. A photon of dIν jν
frequency ν crosses the tube, travelling at an angle θ from the y = −Iν + . (44)
dτν aν
axis. (right) Emitting matter in its rest frame K 0 in the tube. The
photon with frequency ν0 now appears to cross at angle θ0 . Since for the optical depth, it is dτν = aν ds, eq. (44) along with
(42) and (43) gives
dIν jν
problems, unless restated in such a form (Misner et al. 1973). We = −aν Iν + 3 . (45)
begin by considering the phase space number density ds ν
In order to improve this, we also need a way to implicate the
N path length variation ds/dλ. By using the projection tensor hαβ =
N= , (33)
V gαβ + uα uβ , we have the photon velocity (υ0 )α in the fluid frame
where N is the number of particles under examination and V the K 0 as
phase space volume they occupy. By taking into account Liou- υ0 α = hαβ kβ = kα + kβ uβ uα ,

(46)
ville’s theorem in curved spacetime, we have that
where uα is the fluid 4-velocity. By the above, we obtain that
dV
=0 (34) ds
q
dλ = − υ0 α obs = − gαβ (υ0 )α (υ0 )β = −kβ uβ obs

(47)
and by combining the above with the conservation of particle dλ obs
number along the world line of the bundle, we obtain that and for the frequency ratio, it is
kβ uβ obs

dN ν
= 0, (35) = , (48)
dλ ν0 kα uα |λ
which is the collisionless Boltzmann kinetic equation and hence (Younsi et al. 2012). Quantities with an accent such as ν0 above,
N is Lorentz invariant. The phase space volume is are henceforth considered measured in the local rest frame.
dV = d3 xd3 p = dAdth3 ν2 dνdΩ (36) Combining these, we have that
ds ν
and hence = −kα uα |λ . (49)
dλ ν0
dN
N= 3 2 . (37) Eq. (45) combined with (42), (43) and (49) gives the differential
h ν dAdtdνdΩ
form of the invariant RTE equation as
Since the specific intensity is defined as
jν 0
!
dIν
hνdN = − kα uα |λ −aν 0 Iν + 0 3 . (50)
Iν = , (38) dλ ν
dAdtdνdΩ
Integration of the above gives the solution for the Lorentz invari-
we can see that ant specific intensity
1 Iν Rλ
N= (39) Iν (λ) = Iν (λ0 ) e λ0 ν
a 0 (ζ) k uα | dζ
α ζ
−
h4 ν 3 Z λ 0 Rλ
jν (ξ) ξ aν 0 (ζ) kα uα |ζ dζ
and therefore the Lorentz invariant specific intensity is − e kα uα |ξ dξ (51)
λ0 ν03
Iν
Iν = = Lorentz invariant. (40) The optical depth can be calculated as
ν3
Z λ
The optical depth, used to count photon fractions, is a scalar τν (λ) = − aν 0 (ζ) kα uα |ζ dζ (52)
quantity and is thus invariant λ0
τ = Lorentz invariant. (41) and equation (51) can then be rewritten as

Z λ 0
In order to find the Lorentz invariant absorption coefficient, jν (ξ) −[τν (λ)−τν (ξ)]
we use fig. 2. The tube width d = d0 is the same in both the lab Iν (λ) = Iν (λ0 ) e −τν (λ)
− e kα uα |ξ dξ. (53)
λ0 ν03
2.6. Intensity of single emission source radiation

In this paragraph we will describe the way to estimate the radi-
ation received by a target, when said radiation is emitted by a
single emission source. This means that the photons are emitted
by a skin surface of the accretion disk and do not traverse any
of its material. This happens in the case where the disk is totally
optically thick. Various parts of this procedure have been stud-
ied in the literature: Abramowicz, Ellis, & Lanza (1990) studied
radiation emitted by a central nonrotating star in Schwarzschild
spacetime. Miller & Lamb (1996) also studied the environment
around emitting stars and expanded this work by examining non-
rotating and rotating masses and radiating sources. We have also
studied this subject in the previous work Koutsantoniou & Con-
topoulos (2014) examining fewer examples of totally opaque
disks and a single observer position for each model-BH spin set.
Here, apart from expanding into semi-opaque disks discussed
later on, we expand our analysis to more disk models and instead Fig. 3: Emission of a photon from the hot accretion disk. The
of having a single observer at the ISCO of each model-spin set, photon is emitted along n̂ from an element of matter moving
we fill the entire region of the system with a large amount of with velocity υφ̂ . The photon is emitted at an angle ψ relative to
observers in various positions. accreting material motion.
As we saw previously in eq. (40), the Lorentz invariant spe-
cific intensity is Iν = Iν /ν3 and thus for the frequency integrated
specific intensity it is Concluding the necessary transformations, we have the fac-
tor required for the implementation of the frame dragging ef-
I1 I2 fects, which is
= 4 (54)
ν1 ν2
4
1 + ωrec pφ /pt
νrec = νem , (58)
for two random points. From this, we have that for the emitted 1 + ωem pφ /pt
frequency Iem and the received frequency Irec of a photon, it is
!4 where pa are the photon covariant four-momentum components,
νrec which are also conserved quantities and thus do not need more
Irec = Iem . (55) subscript specifications. The ratio pφ /pt depends only on the di-
νem
rection of the photon emission and from the previous statement
Note here that the frequency fraction νrec /νem that appears is also a conserved quantity. Combining the above, we have for
above does not depend on the frequencies involved, but only the received frequency that
on the spacetime and the photon’s emission angle. This fre-
gtt,em 1/2 1 + ωrec pφ /pt
!
quency fraction includes the effects of three different phenomena 1
caused mainly by the spacetime properties. Firstly, it includes νrec = νem (59)
1 + ωem pφ /pt γ 1 − υφ̂ cos ψ

gtt,rec
the effects of gravitational time dilation, which appear both in
Schwarzschild and Kerr spacetimes. It also includes the frame and for the frequency integrated specific intensity that
dragging frequency shift due to the spacetime’s differential rota-
tion that appears only in a Kerr spacetime. Finally, it includes the gtt,em 2 1 + ωrec pφ /pt 4
! !
1
Doppler shift caused by the motion of the source’s emitting sur- Irec = Iem . (60)
face. This can exist in both Schwarzschild and Kerr spacetimes. gtt,rec 1 + ωem pφ /pt γ4 1 − υφ̂ cos ψ4
The gravitational time dilation causes a frequency shift for
the received frequency νrec given by In fig. 4, one can see a breakdown of the process described
above, where for visual simplicity we assume that the emission
gtt,em 1/2
!
of photons is done by a central object.
νrec = νem , (56)
gtt,rec
where νem the emitted frequency. 3. Accretion tori
Assuming the photon source moves azimuthally with negli-
In this section we present the accretion tori we used for our
gible radial and poloidal velocity components, the Doppler shift
codes. Some of the tori are optically thick while others are strat-
due to the emitting surface motion then introduces a change in
ified and semi-opaque.
frequency
The first group, are used to represent tori that increase their
1 density abruptly and very close to their outer surface. They are
νrec = νem , (57) mostly rotationally supported and totally optically thick. This
γ 1 − υφ̂ cos ψ practically means that there is no reason to solve the RTE for
−1/2 these tori and instead another method of calculation must be
where υ = υφ̂ the source velocity, γ = 1 − υ2

the emitting used. These tori could be used to describe physical tori that are
material Lorentz factor and ψ the angle between the emitting either cold or compact or both. An interesting case they could
matter velocity and the photon emission direction (fig. 3). Notice also be used to describe, is transient stages of the x-ray binary
that both the γ factor and the ψ angle are measured in the ZAMO systems. During the quiescent stages of these systems, their ac-
frame at the point of emission. cretion disks tend to be cooler and at a larger distance from the
Fig. 4: Schematic of the frequency changes. A Doppler shift is used to move from the frame comoving with the emitting surface
to the LNRF at the radius of the surface, taking into account the emission source rotation. Then we move to the receiving particle,
accounting for two more changes in frequency, the gravitational time dilation due to the change of radial distance from the source
to the target and the different effects of frame dragging, again because of the change in radial distance.
compact object. They remain in a similar state for an indefinite the Eddington luminosity LEdd = 4πGMm p c/σT , then we have
amount of time and at some point later on, they start rising their an Eddington accretion rate
temperature and swelling up while reducing their density and
density gradient. From that point on, they can no longer be de- 4πGMm p
ṀEdd = , (61)
scribed by these opaque models but require usage of the semi- cσT
opaque tori models.
The latter group of tori, the semi-opaque ones, describes where G is the gravitational constant, M the disk mass, m p the
more common and familiar perhaps cases of accretion disks. mass of the proton, c the speed of light and σT the Thomson
There is a measurable density and temperature gradient. For cross section. Assuming a large enough amount of scatterings,
these cases, one must use a ray tracing process and solve the the disk material can be adequately well described by blackbody
RTE along the photon trajectory, calculating how much radia- radiation and then for the temperature it is
tion is produced by the hot material in every step and how much !1/4
of this is absorbed away by it. Depending on the direction and 3G ṀM
angle of motion of the travelling photon, it can be at times ab- T= r−3/4 , (62)
8πσ
sorbed by the disk material and at other times, it can traverse part
of the disk without it being absorbed. This means that at some where σ the Stefan – Boltzmann constant (see Longair 2011).
points the disk is optically thick and at other ones, optically thin, Notice here however, the problem that arises if we simplis-
and hence the name semi-opaque. In these cases, effects such as tically hypothesize the above. If we assume that the accretion
limb darkening and transparency can be observed. luminosity of the object is equal to the Eddington luminosity,
the disk then cannot be geometrically thin. This is because, as
the accretion luminosity increases, the radiation pressure exerted
3.1. Optically thick accretion tori onto the material keeps getting larger and finally comparable to
In this paragraph we will describe the models we used for opaque local gravitational forces. As this happens, the disk keeps inflat-
tori. Some of these tori are built by assuming simplistic disk ing by gaining height and width and thus gradually turning into
cross section shapes, such as polygons. Others are more com- a geometrically thick and optically thin torus. The easiest way
plex and are built self-consistently by assuming that the material to bypass such problems, is to simply assume that the accreting
of the disk is supported and kept in place by its rotation. object only has a fraction of the Eddington luminosity
Optically thick accretion disks are in general geometrically
Ṁacc = ṀEdd (63)
thin. This is caused by the “inefficiency” of the radiation. Since
the material of the disk is opaque, the radiation transmitted by and thus, it is
its hot components cannot reach other, further away parts of the
disk. This results in each local material component to have a !1/4
3G ṀM
significantly lesser “inflating” radiation pressure element than a T= r−3/4 . (64)
“deflating” gravitational force element. The result is an accretion 8πσ
disk of smaller geometrical thickness and a much larger pressure
gradient, specifically close to its outer surface. Finally, after having picked any of the disk temperature pro-
A matter of particular importance is the surface temperature files, we can have its effect on the emitted photon from
distribution we assume for these tori. For our calculations, we σS
considered two separate cases, an isothermal disk and a disk Iem = T (rem )4 , (65)
π
whose temperature follows T ∝ r−3/4 . The first case, albeit un-
natural, is the simplest possible, one could imagine and is thus where σS the Stefan – Boltzmann constant.
perhaps easier to understand and conceivably anticipate certain We should also mention here that although in the vast ma-
results. jority of accretion disk simulations and studies we assume so,
The second case is to assume that the disk temperature dis- the disk material in reality is in a far from stationary condi-
tribution is caused by the material accreted onto the central com- tion. There are increased amounts of turbulence, instabilities and
pact object. If we assume that the object’s luminosity is equal to other phenomena taking place, often in miniature scales, that are
frequently ignored. One such phenomenon is the flow and dif- where Γαβρ are the Christoffel symbols that can be calculated
fusion of angular momentum throughout the various disk sec- from the metric (2), using the formula
tors. The main cause of this diffusion is considered to the ma-
terial viscosity and is treated in various ways, one of the best 1 αµ
Γακλ =

g gµκ,λ + gµλ,κ − gκλ,µ . (67)
known and more frequently used ones being the α-viscosity 2
method (Shakura & Sunyaev 1973; Abramowicz et al. 1988).
The torus we have assumed here, as the majority of tori in
This method however, cannot unfortunately be applied to all disk
works of this type, has negligible radial and poloidal velocity
models, such as in non α-disks, where other solutions must be
components, so (66) can be slightly simplified. We then get
found.
for the acceleration components that
Another important phenomenon that is known and generally
mentioned in such works but ultimately ignored, is the existence at ≡ 0 (68)
of magnetic fields. Magnetic field effects usually give rise to very
important effects that affect the structure of the disk and deter- aφ ≡ 0
mine its evolution. One such important effect is for example the ∆
"
Σ − 2r2 t 2 #
2 φ 2
uφ

existence of magnetorotational instabilities (MRI) that can have ar = − M u − asin θu + rsin2
θ
important effects, particularly on the distribution and flux of an- Σ Σ2
gular momentum throughout the torus and its diffusion towards sin θ cos θ 2Mr h t 2
( i2 2 )
θ
the outer layers and components (see Balbus & Hawley 1991). a =− au − r + a u + ∆ uφ
2 φ
.
Σ Σ2
In our work, we have considered so far six different models
for optically thick accretion tori. Specified by their given names, The first two of the above equations are in accordance with
we have the models band, disk, slab, wedge, torus and opaque our initial assumption that the disk is stationary and axisym-
rotationally supported torus: metric respectively. The last two equations are what can give
the surface of constant acceleration, the isobaric surfaces us-
(a) Band: a ring that consists of a vertical surface of half-height h ing
(from its highest point to the equatorial plane) at the distance
of the respective ISCO for the selected spin parameter. The aα uα = 0. (69)
half-height h can be freely chosen without restrictions. The
ring radius can easily be modified. Using the above equations, this gives
(b) Disk: an infinitesimally thin disk at the equatorial plane. Its
Σ r r
inner radius is at the distance of the ISCO and its outer radius ar ur + aθ uθ = a u + Σaθ uθ = 0. (70)
at twice this distance. The innermost and outermost radius of ∆
the disk can easily be adjusted. This leads to the pair of differential equations
(c) Slab: a disk of half-height h from the radius of the ISCO to
a distance of three times the ISCO radius. The cross section dr Y
of the disk is a rectangle (fig. 5). The half-height h can be = √ ,
dξ Y 2 + ∆X 2
freely chosen without restrictions. The innermost and outer-
most radius of the accretion disk can easily be adjusted. dθ X
(d) Wedge: a disk whose cross section is an isosceles trapezoid =−√ , (71)
and is centered over and under the equatorial plane. Its inner
dξ Y + ∆X 2
2
radius is equal to the radius of the ISCO and its outer radius where
to three times that. The disk is constructed in such a way that Σ − 2r2 −1 2
an angle with its vertex at the location of the BH (the origin) X=M Ω − asin2 θ + rsin2 θ
extending outwards, reaches the ISCO cylinder intersecting Σ 2
a ring of half-height h. The angle sides continue extending "

2Mr −1
#
2 2

outwards in the same direction until crossing the outer edge Y = sin θ cos θ aΩ − r − a + ∆
2
(72)
Σ2
of the disk (fig. 5). The disk half-height h, as well as its inner
and outer radius can easily be modified. and Ω = uφ /ut is as before the material’s angular velocity.
(e) Torus: a disk with a circular cross section. The center of the The last information necessary to solve the above and have
circle is at coordinates (r, θ) = (2rIS CO , π/2 ) and its radius is the resulting torus, is its inner edge at the equator. This is
equal to the ISCO radius. The disk inner edge is therefore at given by solving the equation for marginal stability orbits
r = rIS CO and the outer edge at r = 3rIS CO . The cross section
Σ − 2r2 2 a2 Mrsin2 θ
" 2 !#
center and radius of the disk can be adjusted at will (fig. 5). r
(f) Opaque rotationally supported torus (ORST): a rotationally 2aMsin θ4
− r +a + 2
Ω3 +
Σ Σ2 Σ
supported torus. This disk is one of the more complex cases
considered for optically thick disk examples. The disk we
 ! 
 Σ − 2r2  6Mr r + a
2 2
 
2Mr 

  
+ + Σ +
 
consider here is stationary and axisymmetric and has its ro- M 3∆ −  1 − r

·
Σ Σ Σ
2

 
 
tation axis aligned with the rotation axis of the black hole.


   


In our work we have assumed that the two angular velocity
vectors are collinear, but it is simple to consider the oppo- Σ − 2r2 6aM 2 rsin2 θ
· sin2 θΩ2 − Ω+
site case in order to study retrograde disks. We then assume Σ2 Σ
that the disk acceleration is what creates this setup and spec-
∂Ω Σ − 2r2
!
2Mr
ifies its shape. The acceleration along the particle trajectory +∆sin θΩ 2
−M 1− = 0. (73)
is given by ∂r Σ2 Σ
Duα In order to solve eq. (73), we must define the angular velocity
aα = = uα ;β uβ = uα ,β uβ + Γαβρ uβ uρ , (66) function. We use here the angular velocity profile proposed
dτ
and explained in (Fuerst & Wu 2004, 2007) and Younsi et al.

(2012)
√ r n p
M K
Ω ($) = Ω (r sin θ) = √ , (74)
(r sin θ) + a M
3/2 r sin θ
where rK is the equatorial plane radius at which the material
moves with Keplerian velocity and the parameter n p corre-
sponds to pressure forces and is responsible for the geometry
of the torus determining its thickness. Tori solutions for var-
ious rK and n p values are shown in fig. 6, while the selected
tori used in our simulations are displayed in fig. 5.
The last thing remaining for accretion disks of this kind, is to
calculate the frequency integrated specific intensity. This is done
by using the method described in §2.6 and then applying §2.4 to
determine the radiation flux or force and the ensuing acceleration
caused by the disk’s hot material.
3.2. Semi-opaque tori

In this paragraph we refer to the semi-opaque and transparent
disk models considered in our work. Again as before, some of
the models are more simplistic than others and some are based
on specific physical conditions and are more complex.
For our research, we considered up to this point five disk
models that could fit in the semi-opaque or transparent category.
Specified by the name of the considered models, we have the
following cases:
(a) No torus: there is no accretion disk around the central black
hole. Photon trajectories continue until they either cross the
event horizon of the black hole or escape the system entirely
by crossing an adjustable outer radius boundary.
(b) Semi-opaque pressure supported polish doughnuts (PS PD):
a stationary and axisymmetric pressure supported pol-
ish doughnut. An accretion torus constructed following
Abramowicz et al. 1978 and Kozlowski et al. 1978. As in
previous tori, we assume here that the material has no sig-
nificant radial or poloidal velocity
components and thus has
α t φ

a 4-velocity u = u , u , 0, 0 . We follow Younsi et al.
2012 and assume that the torus has a polytropic equation
of state P = κnΓ ,where n the material number density and
h i4 1/3
κ = ~c 45 (1 − β)/ π2 µm p β , with ~ the Planck con-
stant, c the speed of light, µ the mean molecular weight, m p
the proton mass and β the ratio of gas pressure to total pres-
sure. The torus is then described by
∂r ξ (r, θ) = −ar (r, θ)
∂θ ξ (r, θ) = −aθ (r, θ) , (75)
where ξ (r, θ) = ln Γ − 1 + Γκn(r, θ)Γ−1 is a function of the
h i
disk number density n (r, θ) and the acceleration components

are aα (r, θ) = uα;β (r, θ) u(r, θ)β . Cross sections of the torus
used are shown in fig. 7.
(c) Translucent PS PD: a translucent pressure supported polish
doughnut. It is the same as the above torus, but displays no
absorption of photons by the disk material.
(d) Semi-opaque LFM torus: a stationary and axisymmetric
semi-opaque torus of circular cross section. The cross sec-
tion center lies on the equatorial plane at rcenter = 2rIS CO Fig. 5: Opaque accretion disks cross sections for the spins stud-
and the torus cross section has a radius rtorus = rIS CO . The ied.
torus thus stretches from an inner radius of rinner = rIS CO
The disk models discussed above are responsible for giving

us most importantly the material’s number density n (r, θ). From
that, one can then obtain other useful quantities for the matter,
one of which is the material temperature. Following regular pro-
cedures, we can have here
#1/3
~c 45 (1 − β)
"
T (r, θ) = ρ(r, θ)1/3 , (76)
k π2 µm p β
where ρ the (volumetric mass) density.

Continuing on, we can obtain firstly the necessary material’s
absorption coefficient from eq. (25) as
aν (r, θ) = σν n (r, θ) , (77)
where σn is the absorption cross section best chosen for the pro-
cesses under study. Moving forward on, we can find the corre-
sponding emission coefficient of the material by making use of
the thermal emission and blackbody radiation properties we have
for our assumed disk. The thermal emission assumption gives
that the emission coefficient is given by
jν (r, θ) = aν (r, θ) Bν (T ) , (78)
where Bν (T ) the Planck function and T the corresponding tem-

perature. In order to procure the Planck function, we make use
of the blackbody attributes of the material and have
2hν3 /c2
Bν (T ) = , (79)
exp (hν/kT ) − 1
where k the Boltzmann constant. Combining then the above
equations, we have for jν that
2hν3 /c2
jν (r, θ) = σν n (r, θ) . (80)
exp [hν/kT (r, θ)] − 1
Notice here however, that in the above equations, one could also
add shaping functions to modify the emission and absorption
of the material to study other disk models or perhaps different
physical properties. A more general form of the above functions
could thus be
aν (r, θ) = Cabs σν f (n (r, θ) , T [n (r, θ)] , E [n (r, θ)]) ,
jν (r, θ) = Cem f (n (r, θ) , T [n (r, θ)] , E [n (r, θ)]) . (81)

One case of such choices can be seen for example in Younsi et al.
2012.
Fig. 6: ORST cross sections. In (a) we can see the effects caused After concluding the calculations
described above, we have
by the change of the radius of Keplerian rotation speed rK . In the resulting radiation intensity I r, θ, ã, b̃ . We can then apply
(b) we see the tori shapes and sizes for various values of the the method of §2.4 and obtain the stress – energy tensor and the
parameter n p . flux and force of the radiation.
to an outer radius router = 3rIS CO and has a maximum 4. Algorithms and codes
height htorus = rtorus = rIS CO . The center number density is
In this section we will present the codes we created and used in
ncenter = 1018 cm−3 and decreases to zero moving towards the
our work and explain their capabilities. We will also show many
torus surface. For this torus the product aν · 2rtorus is ∼ 1 − 5
results of the various studies we have performed in assorted se-
throughout the cross section. Images of the cross section are
tups of disk models, spin parameters and other. Before moving
shown in fig. 8.
along, the reader should be made aware of the fact that all of
(e) Translucent LFM torus: a translucent LFM torus of circular the codes used in our work and presented here, were designed
cross section. It is the same as the previous model, but within order to be executed in extremely limited computational re-
out its material absorbing any of the photons crossing it. sources. This has important consequences on the design, speed
and effectiveness required from the codes.
Fig. 7: Polish doughnut number density cross sections. The center of all tori lies at (12, 0) and has a number density of 1018 cm−3 .
The top row shows the disk for black hole spin parameters a = 0 and a = 0.5, while the bottom row for a = 0.9 and a = 0.998.
Fig. 8: LFM model number density cross section for any spin
parameter. The center of the tori lies at (2rIS CO , 0) and has a
number density of 1018 cm−3 .
4.1. Code Omega
The code named Omega was the first created in our work and
its first version was written in 2013. It is so far, one of the most
important parts of all following codes.
The code’s main goal is to find photon trajectories. It
works for a Schwarzschild and a Kerr spacetime, but it can
be easily modified in order to work in other spacetime mod-
els as well, e.g. Kerr – Newman, Reissner – Nordström, Fried-
mann–Lemaître–Robertson–Walker etc.
The code studies various tori models which were referred Fig. 9: Accretion disks created from combinations of previously
and explained in §3.1 and §3.2 and some were also depicted in referred models: wedge & slab, torus & wedge, polish doughnut
figs. 5, 7 and 8. It is also possible to add new accretion disk & torus and polish doughnut & spheres.
models to the code. Adding disks resulting from a combination
of the previously mentioned models, is a quick and easy matter.
Examples of such tasks are shown in fig. 9. or energy is carried to the accretion disk material and the target
Omega solves the particle trajectory equations mentioned in particle.
§2.3 for a photon and finds the trajectory and the point of origin Depending on the environment in which it is used, Omega
of said photon. This point of origin could be on the hot accretion code can have different outputs. In its original form, which is vi-
disk, the black hole event horizon or a location outside and far sual, the program has an interface that allows the user to select
away from the system. If all that the photon trajectory intersects primary properties for the environment, such as the disk model,
with is the event horizon or the system exterior, then no radiation the black hole spin parameter and the disk height. Also, one can
select important options for a trajectory, including its maximum server from each of the angles of the local sky. It then calcu-
length, its point of origin and angle of emission and the two lates the stress – energy tensor for the radiation and creates var-
emission angles ã and b̃. Finally, there are some additional op- ious images and sky maps of the radiation. It continues on to
tical options that include the choice of frame size of the visual run a procedure that increases the code resolution by a signif-
box and the depiction of obscured parts of the outer disk. The icant number of times (see fig.12). In our runs, we choose to
code’s dynamic output picture shows the black hole event hori- increase the resolution by five times but it is easily possible to
zon and its ergosphere, the accretion disk and the requested pho- change that to a higher or smaller number. It is also possible to
ton trajectory. The photon trajectory is drawn in different styles choose to run more than one applications of the process. This
and colors for escaping particles, particles infalling in the black was a rather hard procedure to create but was necessary, because
hole and particle trajectories starting from the accretion disk. In as mentioned in the beginning of this section, the available com-
addition, some trajectory information are displayed on the pic- putational resources are severely restricted and thus running the
ture, including the photon energy and angular momentum, the desired quality and amount of codes is not possible since that
trajectory’s Carter constant and the initial “TARGET” momen- would require many months or some years.
tum magnitude. Finally, at the bottom appear the requested ex- The program continues on to calculate the radiation flux and
periment “SOURCE” information, the received photon source force four-vectors applied on various observers. This includes
coordinates and the momentum magnitude at emission. All the observers at rest in the local frame, rotating with ω = ω (r, θ) (eq.
above can be seen in fig. 10. 5) as seen from infinity, observers in circular orbits moving with
Another version of the Omega code is available, where in- Ω (eqs. 15 and 74) depending on the selected disk model, ob-
stead of a single trajectory, a bundle of photon orbits are dis- servers on accreting material attempting to mimic the SANE and
played. Again, the photon orbit representation depends on its MAD models (Narayan et al. 2012; Penna et al. 2013; Narayan
origin, as described above. Since more than one orbits are dis- et al. 2003), and finally observers in two possible outflow re-
played in a single picture, no additional numerical information gions close to the system rotation axis and above a certain height
for the trajectories are displayed in this case. In addition to the (fig. 13). Let us note here that, even though we have made our
previous version’s options, in this one the user also has the extra choice on certain velocity profiles for the various disk models,
choice of how many photon orbits will be displayed. An example this can easily be adjusted to suit the needs of other disk models
of the code’s output image can be seen in fig.11 with different material velocities. Notice also, that such a change
The last Omega code version is the most plain, yet by far the would only affect the final step of the program, that only takes
most important one of all. It is not a code anymore but a func- a few minutes to run, and does not require for all the code to be
tion which solves the photon equations of motion and outputs re-executed from the beginning, something that would require
only numerical data and no visual information. It runs for a sin- hours or a day.
gle photon trajectory per execution and returns key information The program then concludes by outputting its results and cre-
for this trajectory. Firstly, it gives the existence or absence of in- ating various save files. Firstly, it incorporates its minimum and
coming radiation in the requested direction. It then also gives the maximum radiation values to the relevant files for each model.
precise photon trajectory numerically as well as the coordinates Then, it writes on the selected model’s four and three-force val-
of the emission source. The reason Omega is considered so im- ues files for the various observer velocities. It then goes on to
portant is because it is called and used by all almost all of the create the numerical raw data file for the simulation and a Moll-
subsequent codes we have created so far and is thus one of the weide map picture for the radiation of the specific observer lo-
most important parts of our entire work. cation. Finally, it creates a graphical output file that includes
Before moving on to the next family of codes, we should short but important information on the simulation it is describ-
mention here that even though Omega is built in order to study ing, along with several pictures of the important matrices for the
photon trajectories, it is fairly easy to modify it in order to make simulation. It also includes tables for the assorted profiles’ force
it follow massive particle trajectories. components and the radiation stress – energy tensor matrix.
In Table 1, we can see the amount of Infinity codes ran for
each model and spin parameter. In fig. 14, we can see some Moll-
4.2. Code Infinity
weide projection sky maps of the incoming radiation in the same
The Infinity code is the next step after Omega. Its first version models mentioned in Table 1. Also, videos of flights around and
was written in 2013 and it has had many improving versions if possible through the disk for all the models we studied, can
since. It started off as a code that would run and give the radiation be found on youtube.com, under the name of this work’s creator
information about a single target at a single and specific point in "Leela Elpida Koutsantoniou".
space and is now an extended and fully automatic program that Continuing on, we can see figs. (19) – (24) and examine them
can run entire simulations for a large stream of different points in detail. These pictures show the 4-force components sign and
and situations that do not even require being of the same kind. magnitude at several points outside, near and inside (if possi-
It automatically creates various files, saving numerical data, im- ble) the various disk models mentioned in Table 1 and fig. 14.
ages and information possibly useful for other runs or aggregated Each line of these pictures refers to target particles belonging to
results. It also has another interesting ability worth mentioning one of the aforementioned velocity groups: “halo” includes par-
here. Infinity code can be stopped at any time during execution ticles moving only with ω due to the spacetime rotation, “disk”
and then restarted later on when the system is again available includes particles in circular orbits with speed Ω, “SANE” and
and the calculations simply continue on from where they were “MAD” refer to particles in inspiral orbits, attempting to mimic
left off. This is very useful in many instances, including cases of these two accretion situations and “outflow” refers to rotating
power outage. target particles that also have a radial velocity component point-
In the Infinity code, the user chooses a point of interest any- ing outwards.
where in the system under study and the program sets an ob- The first four lines of each picture are poloidal plots: the hor-
server there. It then scans the entire sky around said observer, izontal axis measures the cylindrical distance $ from the central
solving the RTE and returning the radiation reaching the ob- axis and the vertical measures the height z from the equatorial

Fig. 10: The control interface of the Omega code, as described in §4.1
izontal axis is z and measures the height from the equatorial

plane and the vertical axis is $, measuring the cylindrical dis-
tance from the system rotation axis.
Also, each of the pictures consists of two parts: part A shows
the force components measured in the BL frame and part B mea-
sured in the ZAMO frame. Even though the plot color schemes
change for the various disk models, the reader should keep in
mind that in all cases, point zero is always assigned to white
color. In addition, warmer colors (hues of red, pink, orange etc.)
are used to show positive force values and cooler colors (blue,
green, cyan etc.) to show negative force values. The columns in
each plot group show in order the distribution of the t-, φ-, r- and
θ- force components.
Some of the results presented by the plots are expected, while
others are not. Notice that the divergence of the BL halo f t plot
Fig. 11: Output image of the photon bundle edition of the Omega from zero gives us a measurement of the numerical and compu-
code. The user here has selected the depiction of 12 photon tra- tational errors present in our calculations. These errors appear to
jectories. The solid (blue) lines in the center are trajectories of be around 15 orders of magnitude bellow the forces we are at-
photons emitted from the accretion disk. The dashed (purple) tempting to calculate, so they are well within acceptance limits.
lines in the left are trajectories that would originate from the Looking at the first column of a picture group, we can see the
event horizon and the dotted (red) lines in the right are photons effects of the radiation on the energy absorption amount and rate
that should have been emitted from somewhere outside the sys- by the target particles. The images provided by the simulations
tem. are in qualitative agreement with the expected results, since we
can see the energy transfer getting stronger, the closer we get
to the central object and thus the faster we rotate and absorb
plane, which is also a symmetry plane for the accretion disks. energy as well. The second column shows the azimuthal compo-
The fifth and last line of each picture shows plots similar to the nent of the radiation force. Notice that this appears to be nega-
previous ones, but these plots are rotated 90 degrees: the hor- tive in the vast majority of cases and changes only for areas very
Table 1: The number of Infinity code executions for the various disk models and spin parameters. The first four columns are some
of the opaque disks we described earlier. The slab is discussed in §3.1 (c), wedge in (d), torus in (e) and ORST in (f). The remaining
columns are two of the semi-opaque disks we mentioned in §3.2. PD can be found in (b) and LFM in (d).
—————– Disk model —————–

slab wedge torus ORST LFM PD
s 0 168 170 163 184 496 426 Total
p 0.5 182 177 176 188 478 456 number
i 0.9 156 158 148 179 424 460 of
n 0.998 - - - 191 414 460 runs
Sums 506 505 487 742 1812 1802 5854
close to the ergosphere, or for cases with large black hole spins Table 2: Declinations of the Tranquillity code estimation of the
(a ≥ 0.9). This is in agreement with the results presented in our disk inclinations for various black hole spins (columns) and in-
previous work Koutsantoniou & Contopoulos 2014. The third clinations (lines). The spin is shown in the top row and the in-
column shows the radial 4-force recorded for each disk model clination, set by hand for the disk, is measured in degrees and is
and it is also in qualitative agreement with the expected results. displayed on the first column. The values in the matrix are the
When the target is further out the disk, it receives more radiation errors in the Tranquillity code execution and are also measured
coming from its interior local sky hemisphere than its exterior in degrees. Two of the three worst cases of erroneous results,
one, and is thus pushed further outwards. For a target further in- caused by the echo ring, are shown here, but are still rather small.
side the disk, the opposite occurs: the outer hemisphere receives
more radiation than the inner one, and the particle is pushed to- 0.2 0.7 0.9
wards the black hole. The fourth column shows the θ- force in 0◦ 0 0 0
the poloidal plane and its effects on the system particles. No- 10◦ 0.396 0.191 0.071
tice that as expected, this force disappears in the equatorial plane 20◦ 1.247 0.683 0.554
and becomes negligible, the closer we get to the vertical rotation 30◦ 1.276 0.815 0.677
axis. In almost all the cases this force is negative, which means 40◦ 1.471 0.887 0.695
it pushes the material towards the rotation axis. In the very few 50◦ 1.200 0.917 0.667
cases where the force sign is positive, note that these points are 60◦ 1.249 0.776 0.555
in very close proximity to the ergosphere and can only see a 70◦ 0.987 3.874 2.002
small part of the sky around them. Hence, the received radiation 80◦ 0.534 0.348 0.370
is governed by the Einstein ring of the opposite side of the disk, 90◦ 0 0 0
across the black hole.
Some other important notes for these pictures are also the
following. First of all, even though the forces may seem small hole (fig. 15). Depending on what this ray will meet along its
or negligible, we should always keep in mind that these forces path, it returns information about its origin and the radiation re-
∼ 10−13 dyn act upon the electrons and therefore a mass ∼ 10−27 . ceived.
This means that the acceleration due to radiation is ∼ 1014 cm/s2 This code is similar to the aforementioned Infinity code but
and it is primarily acting upon the material electrons since is practically its complement. Elysium therefore has equal qual-
fe / f p ∼ (m p /me )2 . Also, the reader should be aware that the ity and quantity of capabilities, options and result information
large red and blue/green dots in the plots act in the same manner as Infinity (fig. 16). The difference of the two codes is literally,
as the diffuse colors and represent about half of the runs. These as well as figuratively, a point of view. Infinity starts at the sin-
dots are there only to help us distinguish between positive and gle end point of ray trajectories and integrates the equations go-
negative value forces when the colors are perhaps too faint to ing backwards in time in order to find if the intersecting light
see, e.g. close to zero, where the force colors approach white. pathway can possibly traverse a light source at any point inside
Finally, looking closely at these pictures, one can notice that the system in question. Elysium has multiple starting points, the
the t-force zero curves are not located at the same places as the “pixels” of the screen, and integrates in a normal way, going
φ-, r- and θ-force zero curves. This was not an anticipated event forward in time to see if an observer sitting at the pixel’s loca-
and it is worth mentioning explicitly since it gives rise to per- tion can see any part of the disk. There are various advantages
haps unexpected 3-force components close and through the disk. and disadvantages in both methods and depending on the type of
(High resolution images of these plots can be found in the elec- information required each time, we can choose the most appro-
tronic form of this work and can be zoomed in and studied better priate and fast code of the two.
and easier.)
4.4. Code Tranquillity
4.3. Code Elysium
Tranquillity code was written in 2019. Its main purpose was to
The Elysium code was written in 2014 and kept up to date and find the quickest way possible to have an estimation of the ac-
improved since. Its main goal is to design and create a record- cretion disk inclination and to see if an assessment of the central
ing screen a specifically user selected distance away from the black hole spin is possible.
black hole and accretion disk system. Depending on the selected The calculations for the inclination were rather successful,
program resolution, the screen has the corresponding amount of since the average declination was below 0.8 degrees for the 120
“pixels” and from each of those, a light ray is emitted perpen- cases of different inclinations and spin parameters examined (2).
dicularly to the screen and moves towards the disk and the black However, about 3% of the cases examined, gave inclination er-
Fig. 13: The various regions of a system we consider in a run.

“Halo” only includes observers at rest in the local rest frame
and rotating with ω. “Outer disk” is the main accretion disk of
the system. Material there can move in four different ways: at
rest in the local rest frame rotating with ω, purely azimuthally
with angular velocity Ω, or by slowly infalling in a manner like
SANE and MAD. “Inner disk” is the region of infalling material
and thus includes only matter at rest in the local rest frame and
infalling material mimicking SANE and MAD. The “Outflow”
region is part of a cylinder centered at the system’s rotation axis
that includes the black hole’s ergosphere and has thus a radius
r = 2M. The actual region is considered to exist above a certain
height, further away from the black hole. In this region, we can
have matter at rest in the local rest frame and matter flowing
outwards with a certain velocity. We may also consider a sub-
region there, a narrower cylinder with radius r = M with stronger
outflow and faster velocity.
rors above average. The estimation errors were about five de-
grees or less, and were caused by specific conditions of the setup
that had a black hole first echo ring appear in very particular and
peculiar locations. The model selected for the accretion disk of
the object in question, does not appear to play any significant
part in the inclination assessment so far.
In fig. 17, we can see the divergence results given by the ex-
ecution of 120 runs of the Tranquillity code for an accretion disk
of varying inclinations, ranging from 0◦ up to 90◦ , compared to
the line of sight. An angle of 90◦ , means the disk is seen edge
on. The code was ran for inclination angles every 5◦ , plus runs
for an angle of 89◦ . There is a clearly visible trend for the diver-
gence evolution of the same black hole spin across the various
inclinations.
An important note here is that the disk model adopted plays
a very important part in the divergence calculations. Since the
inner edge of the disk examined can be at very different radii,
depending of the model examined, different disk models will
follow different divergence plots, directly affected by the disk’s
inner edge radius (see Abramowicz et al. 2010). This however,
Fig. 12: Three instances of the first version of the resolution en- does not prove to be an insurmountable problem, since it only
hancement process. The program was initially written in order to seems to cause a recalibration to the divergence plot.
improve the appearance of the output images. However, so much Finally, keeping in mind the very good quality of the afore-
attention and detail was put in it, that it proved to be good enough mentioned Tranquillity inclination results, one could use obser-
to improve the simulation results themselves. The program func- vational data for an estimation of the inner edge of the accretion
tions were thus rewritten in order to work in the appropriate way disk, and by moving vertically in the appropriate divergence plot,
to enhance the run results by a number of times. can have an assessment of the central black hole spin parameter.
Fig. 14: Mollweide projection sky maps of the frequency integrated specific intensity. Each image shows the run results for a
different disk model. Below each image, the model’s color scale is displayed from minimum to maximum.
4.5. Code Burning Arrow commas used for partial derivatives. As in Newton’s equation,
the zero term on the right stands for the absence of acceleration.
The Burning Arrow code was first written in late 2019 and has Solving the above, gives the various geodesic solutions that de-
had some improved versions since. Its main purpose is to study scribe among others, special circular orbits such as the ISCO, the
the black hole massive particle orbit degradation due to the hot photon sphere etc. The above equation can also be rewritten as
disk radiation.
In order to study the particle motion, the code must solve the d 2 xµ α β α β µ
µ dx dx t dx dx dx
general relativistic equations of motion, equivalent to the Classi- + Γ αβ − Γαβ = 0, (83)
cal Newton’s Laws of Motion. Starting from the first law, in ab- dt2 dt dt dt dt dt
sence of general relativistic forces, the particle in question will by using the chain rule in order to have derivations by the coordi-
follow a geodesic through spacetime. This geodesic obeys the nate time x0 = t. This however has proved to be more prominent
equation to error accumulation in our work and so we choose to work with
the proper time equations instead.
d 2 xµ
+ Γµαβ uα uβ = 0, (82) The above equation gives later on, rise to the equivalent of
dτ2 Newton’s second law of motion and a way to calculate the ac-
α celerations present in a problem, such as those generated by the
where τ is the proper time, uα = dx dτ is the 4-velocity and radiation forces in our case. As in the Infinity code, we con-
µ 1 µν

Γαβ = 2 g gαν,β + gβν,α − gαβ,ν the Christoffel symbols, with sider different cases of velocity behaviors, such as the SANE
Fig. 15: Images that show the workings of the Elysium code. On the left is a run for a 7 × 7 pixels screen with j = 0.5, at an
inclination 45◦ . On the right is the center pixels column of a 100 × 100 screen with j = 0 at inclination 85◦ . The continuous lines
are light rays that intersect the accretion disk and thus carry radiation, while the dashed ones end up in the event horizon and the
dotted ones at infinity, attributing nothing to the radiation total.
and MAD models as mentioned before, in addition to the typical adding accretion and radial movement to the material. The aris-
Bardeen et al. 1972 angular velocity profile. ing radiation field close to the BH and the disk is studied from
In the beginning of the execution and depending on the study close up and in detail and new information is gathered about
requested, the code looks up into the Infinity code result files and such environments. The radiation forces, acting firstly upon the
reads the appropriate ones, relevant to each case. It then uses material electrons and then upon the protons (or ions), prove to
these results to generate proper functions that give the radiation be larger and stronger than expected and create accelerations far
stress – energy tensor at every point of spacetime used in the from negligible. This means that radiation effects should not and
problem at hand. This gives the 4-acceleration due to radiation cannot be ignored in calculations for such objects or similar se-
aµrad at each point from the extension of the aforementioned equa- tups.
tion as A matter worth perhaps looking further into, could be the
effects of these significant radiation forces in the dynamics and
d 2 xµ the stability of the disk itself. The modifications brought upon
aµrad = + Γµαβ uα uβ . (84)
dτ2 these systems could be noteworthy, since long considered con-
cepts such as the ISCO, must take radiation effects and perhaps
Since the radiation acceleration is known in spacetime from the corrections into account. Notice also that, since the codes writ-
previous results, one can solve these 8 differential equations and ten are adjustable, various other new accretion, temperature or
get the 4-velocity and 4-positions describing the sought particle radiation profiles could be studied.
trajectory. Another matter worth examining, is the radiation effects and
Finally, the last step for the code is to draw the evolution of possible consequences on outflowing material close to the ob-
the various velocity profile orbits. As an example, we show in fig. ject axis. Is the radiation perhaps of notable importance on the
18 the degradation of equatorial orbits, due to radiation from a very early stages of material outward movement? The disk radi-
semi-opaque LFM torus around a non-rotating black hole. These ation is also pushing the material closer to the rotation axis but
emitted massive particles start from various radii from the axes does it offer any significant contribution to the initial, albeit very
origin and thus encounter different environments of material and small, collimation of the outflow material? Finally, something
radiation density each. that could be considered, is the calculated radiation’s effects on
further, more complex subjects, such as the disk evolution in X-
ray binaries and the observed hardness-intensity diagrams.
5. Conclusions & Discussion
Our study of accretion disk models and radiative transfer in all References
optically thin and thick, and geometrically thin and thick disks Abramowicz, M., Jaroszynski, M., & Sikora, M. 1978, A&A, 63, 221
gave a substantial quantity and quality of information. Some Abramowicz, M. A., Czerny, B., Lasota, J. P., & Szuszkiewicz, E. 1988, ApJ,
of the results are new and in cases unexpected, especially after 332, 646

Abramowicz, M. A., Ellis, G. F. R., & Lanza, A. 1990, ApJ, 361, 470
Abramowicz, M. A., Jaroszyński, M., Kato, S., et al. 2010, A&A, 521, A15
Balbus, S. A. & Hawley, J. F. 1991, ApJ, 376, 214
Bardeen, J. M. 1970, ApJ, 162, 71
Bardeen, J. M., Press, W. H., & Teukolsky, S. A. 1972, ApJ, 178, 347
Bini, D., Geralico, A., Jantzen, R. T., & Semerák, O. 2015, MNRAS, 446, 2317
Bini, D., Geralico, A., Jantzen, R. T., Semerák, O., & Stella, L. 2011, Classical
and Quantum Gravity, 28, 035008
Bini, D., Jantzen, R. T., & Stella, L. 2009, Classical and Quantum Gravity, 26,
055009
Carter, B. 1968, Physical Review, 174, 1559
Contopoulos, I. & Kazanas, D. 1998, ApJ, 508, 859
Contopoulos, I., Kazanas, D., & Christodoulou, D. M. 2006, ApJ, 652, 1451
Contopoulos, I. & Papadopoulos, D. B. 2012, MNRAS, 425, 147
Fuerst, S. V. & Wu, K. 2004, A&A, 424, 733
Fuerst, S. V. & Wu, K. 2007, A&A, 474, 55
Koutsantoniou, L. E. 2014, Master’s thesis, Department of Astrophysics, As-
tronomy and Mechanics, Faculty of Physics, University of Athens, Panepis-
timiopolis Zografos, Athens 15784, Greece
Koutsantoniou, L. E. & Contopoulos, I. 2014, ApJ, 794, 27
Kozlowski, M., Jaroszynski, M., & Abramowicz, M. A. 1978, A&A, 63, 209
Kylafis, N. D., Contopoulos, I., Kazanas, D., & Christodoulou, D. M. 2012,
A&A, 538, A5
Lamb, F. K. & Miller, M. C. 1995, ApJ, 439, 828
Longair, M. S. 2011, High Energy Astrophysics
Miller, M. C. & Lamb, F. K. 1993, ApJ, 413, L43
Miller, M. C. & Lamb, F. K. 1996, ApJ, 470, 1033
Misner, C. W., Thorne, K. S., & Wheeler, J. A. 1973, Gravitation
Mueller, T. & Grave, F. 2009, arXiv e-prints, arXiv:0904.4184
Narayan, R., Igumenshchev, I. V., & Abramowicz, M. A. 2003, PASJ, 55, L69
Narayan, R., Sadowski,
˛ A., Penna, R. F., & Kulkarni, A. K. 2012, MNRAS, 426,
3241
Penna, R. F., Sadowski,
˛ A., Kulkarni, A. K., & Narayan, R. 2013, MNRAS, 428,
2255
Poynting, J. H. 1903, MNRAS, 64, A1
Robertson, H. P. 1937, MNRAS, 97, 423
Rybicki, G. B. & Lightman, A. P. 1986, Radiative Processes in Astrophysics
(Wiley-VCH)
Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 500, 33
Sadowski,
˛ A. 2016, MNRAS, 459, 4397
Sadowski,
˛ A., Lasota, J.-P., Abramowicz, M. A., & Narayan, R. 2016, MNRAS,
456, 3915
Wilkins, D. C. 1972, Phys. Rev. D, 5, 814
Fig. 16: Photograph pictures taken by running Elysium for the Younsi, Z., Wu, K., & Fuerst, S. V. 2012, A&A, 545, A13
ORST. In the left column, there are photographs taken from a
screen with 45◦ and the right with 85◦ inclination. Each line has
a different black hole spin parameter. Moving from top to bottom
these are j = 0, j = 0.5 and j = 0.998. The reader could compare
this picture with the similar Younsi et al. 2012, figure 3.
Fig. 17: The Tranquillity code collective results, where the cir-
cles represent execution results. There is a clear trend for the
divergence evolution for every black hole spin. In this plot, we
also include the three worst code results, depicting them with the
dashed colored lines.

Fig. 18: Degradation of equatorial orbits around a non-rotating BH due to the semi-opaque (LFM) disk’s thermal radiation. With
black is the Bardeen et al. (1972) velocity profile, with cyan and with orange an approximation of a SANE and a MAD model
velocity profile respectively (Narayan et al. 2012; Penna et al. 2013; Narayan et al. 2003).

Fig. 19: (a). Forces for the slab opaque disk model §3.1(c) for BH spin a = 0. BL is on the top and ZAMO on the bottom.

Fig. 19: (b). Forces for the slab opaque disk model §3.1(c) for BH spin a = 0.5. BL is on the top and ZAMO on the bottom.

Fig. 19: (c). Forces for the slab opaque disk model §3.1(c) for BH spin a = 0.9. BL is on the top and ZAMO on the bottom.

Fig. 20: (a). Forces for the wedge opaque disk model §3.1(d) for BH spin a = 0. BL is on the top and ZAMO on the bottom.

Fig. 20: (b). Forces for the wedge opaque disk model §3.1(d) for BH spin a = 0.5. BL is on the top and ZAMO on the bottom.

Fig. 20: (c). Forces for the wedge opaque disk model §3.1(d) for BH spin a = 0.9. BL is on the top and ZAMO on the bottom.

Fig. 21: (a). Forces for the torus opaque disk model §3.1(e) for BH spin a = 0. BL is on the top and ZAMO on the bottom.

Fig. 21: (b). Forces for the torus opaque disk model §3.1(e) for BH spin a = 0.5. BL is on the top and ZAMO on the bottom.

Fig. 21: (c). Forces for the torus opaque disk model §3.1(e) for BH spin a = 0.9. BL is on the top and ZAMO on the bottom.

Fig. 22: (a). Forces for the ORST disk model §3.1(f) for BH spin a = 0. BL is on the top and ZAMO on the bottom.

Fig. 22: (b). Forces for the ORST disk model §3.1(f) for BH spin a = 0.5. BL is on the top and ZAMO on the bottom.

Fig. 22: (c). Forces for the ORST disk model §3.1(f) for BH spin a = 0.9. BL is on the top and ZAMO on the bottom.

Fig. 22: (d). Forces for the ORST disk model §3.1(f) for BH spin a = 0.998. BL is on the top and ZAMO on the bottom.

Fig. 23: (a). Forces for the LFM semi-opaque disk model §3.2(d) for BH spin a = 0. BL is on the top and ZAMO on the bottom.

Fig. 23: (b). Forces for the LFM semi-opaque disk model §3.2(d) for BH spin a = 0.5. BL is on the top and ZAMO on the bottom.

Fig. 23: (c). Forces for the LFM semi-opaque disk model §3.2(d) for BH spin a = 0.9. BL is on the top and ZAMO on the bottom.

Fig. 23: (d). Forces for the LFM semi-opaque disk model §3.2(d) for BH spin a = 0.998. BL is on the top and ZAMO on the bottom.

Fig. 24: (a). Forces for the PD semi-opaque disk model §3.2(b) for BH spin a = 0. BL is on the top and ZAMO on the bottom.

Fig. 24: (b). Forces for the PD semi-opaque disk model §3.2(b) for BH spin a = 0.5. BL is on the top and ZAMO on the bottom.

Fig. 24: (c). Forces for the PD semi-opaque disk model §3.2(b) for BH spin a = 0.9. BL is on the top and ZAMO on the bottom.

Fig. 24: (d). Forces for the PD semi-opaque disk model §3.2(b) for BH spin a = 0.998. BL is on the top and ZAMO on the bottom.

Algorithms and Ray Tracing For Tori Around Black Holes: Leela Elpida Koutsantoniou

Uploaded by

Copyright:

Available Formats

Algorithms and Ray Tracing For Tori Around Black Holes: Leela Elpida Koutsantoniou

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Algorithms and Ray Tracing For Tori Around Black Holes: Leela Elpida Koutsantoniou

Uploaded by

Copyright:

Available Formats

Astronomy & Astrophysics manuscript no.

Algorithms and ray tracing for tori around black holes

2.3. Particle trajectories d pφ

2.5. Radiative transfer

and the matter rest frame, since it is perpendicular to the direc-

τ = Lorentz invariant. (41) and equation (51) can then be rewritten as

2.6. Intensity of single emission source radiation

a ring of half-height h. The angle sides continue extending "

and explained in (Fuerst & Wu 2004, 2007) and Younsi et al.

3.2. Semi-opaque tori

disk number density n (r, θ) and the acceleration components

The disk models discussed above are responsible for giving

where ρ the (volumetric mass) density.

aν (r, θ) = σν n (r, θ) , (77)

jν (r, θ) = aν (r, θ) Bν (T ) , (78)

where Bν (T ) the Planck function and T the corresponding tem-

aν (r, θ) = Cabs σν f (n (r, θ) , T [n (r, θ)] , E [n (r, θ)]) ,

jν (r, θ) = Cem f (n (r, θ) , T [n (r, θ)] , E [n (r, θ)]) . (81)

4.1. Code Omega

Article number, page 12 of 40

izontal axis is z and measures the height from the equatorial

—————– Disk model —————–

Fig. 13: The various regions of a system we consider in a run.

Article number, page 17 of 40

Article number, page 18 of 40

Article number, page 19 of 40

Article number, page 20 of 40

Article number, page 21 of 40

Article number, page 22 of 40

Article number, page 23 of 40

Article number, page 24 of 40

Article number, page 25 of 40

Article number, page 26 of 40

Article number, page 27 of 40

Article number, page 28 of 40

Article number, page 29 of 40

Article number, page 30 of 40

Article number, page 31 of 40

Article number, page 32 of 40

Article number, page 33 of 40

Article number, page 34 of 40

Article number, page 35 of 40

Article number, page 36 of 40

Article number, page 37 of 40

Article number, page 38 of 40

Article number, page 39 of 40

Article number, page 40 of 40

You might also like