Issue 
A&A
Volume 657, January 2022



Article Number  A32  
Number of page(s)  21  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/202140682  
Published online  24 December 2021 
Algorithms and radiation dynamics for the vicinity of black holes
I. Methods and codes
Department of Astrophysics, Astronomy and Mechanics, Faculty of Physics, University of Athens, Panepistimiopolis Zografos, Athens 15784, Greece
email: leelamk@phys.uoa.gr
Received:
28
February
2021
Accepted:
8
October
2021
We examine radiation and its effects on accretion disks orbiting astrophysical black holes. These disks are thermally radiating and can be geometrically and optically thin or thick. In this first paper of the series, we discuss the physics and the formulation required for this study. Subsequently, we construct and solve the relativistic radiative transfer equation, or find suitable solutions where that is not possible. We continue by presenting some of the accretion disks we considered for this work. We then describe the families of codes developed in order to study particle trajectories in strong gravity, calculate radiation forces exerted onto the disk material, and generate observation pictures of black hole systems at infinity. Furthermore, we also examine the veracity and accuracy of our work. Finally, we investigate how we can further use our results to estimate the black hole spin and the motion of disk material subjected to these radiation forces.
Key words: accretion, accretion disks / black hole physics / radiative transfer / relativistic processes
© ESO 2021
1. Introduction
In 2017, one half of the Nobel Prize in Physics was awarded to Barry C. Barish and Kip S. Thorne for the first observation of gravitational waves from merging black holes. In 2019, the Event Horizon Telescope captured the first image of the shadow of the black hole at the core of galaxy M87. And even more recently, one half of the 2020 Nobel Prize in Physics was awarded to Roger Penrose for the discovery that black hole formation is a robust prediction of the general theory of relativity. In the past, black holes may have been considered exotic and mysterious; however, today they are of fundamental importance for astrophysics and are objects of strong ongoing scientific research and investigation. One aspect of that research involves the study of the dynamical effects of radiation in the immediate environment surrounding an astrophysical black hole. One should keep in mind that most known astrophysical black holes are accompanied by accretion disks of a temperature high enough to produce Xrays close to Eddington luminosities in stellar black holes. In fact, before the detection of gravitational waves, this was the only way to detect astrophysical black holes.
In the astrophysical setting of a black hole surrounded by an extended (noncentral) source of radiation (hot accretion disk, corona, and jet), the dynamical effect of radiation, namely radiation pressure, must be taken into account in detail when examining the stability and the evolution of the surrounding matter itself. We can understand this with an analogy from an astrophysically “milder” environment closer to us, namely the Solar System. The dynamical effects of solar radiation were first described by J. H. Poynting (Poynting 1903). Years later, H. P. Robertson (Robertson 1937) properly explained what is now known as the Poynting–Robertson effect and cleared up misconceptions that had puzzled many great scientists, such as J. Larmor for years. This effect is also called a “drag” because, even though the solar photons are emitted radially outward, they cause the dust grains in orbit absorbing them to slowly brake and infall onto the Sun, the Earth, or any other massive object that is close enough. This force is of relativistic origin and is proportional to the absorbing object’s azimuthal velocity. A simple nonrelativistic approximative calculation of the radiation force on electrons orbiting a central source of luminosity L is:
where V_{ϕ} is the target’s azimuthal velocity, c the speed of light, σ_{T} the Thomson cross section, and r the target’s radial distance from the central object.
The goal of this paper is to present the methodology we used and the codes written in order to study the dynamical effects from an extended source of radiation in general relativity, that is way beyond Eq. (1). The first question we would like to look into is what the magnitude of the radiation created by the hot accretion disk is and whether the effects caused by it could indeed be as negligible as they are often considered to be. We would also like to investigate what the effects of these radiation forces can be on the disk itself and what they could mean for the kinematics, the evolution, and the stability of the disk. We can ponder, for example, if this radiation could perhaps be one of the processes that trigger or regulate the accretion processes. Additionally, we explore the possibly distinct effects that different types of accretion disks can have due to their dissimilar geometrical characteristics, temperature, density gradients, etc. Finally, in this work we present the codes written for our study and show some examples of their results. These include geodesics and trajectories, radiation forces, photograph pictures of black holes surrounded by accretion disks, black hole spin estimation, and radiation induced accretion of material. The results will be presented and discussed in more details in the second part of this work.
Some of the first who studied radiation effects within General Relativity in environments relevant to the ones we examine in the present work are Abramowicz et al. (1990), Miller & Lamb (1993, 1996) and Lamb & Miller (1995). All the aforementioned papers, however, unlike our work, consider central sources of photons, which in some cases are also rotating. A different approach to relativistic radiation effects was followed by Bini et al. (2009, 2011, 2015). The work of Sądowski (2009; 2011; 2016) and Sądowski et al. (2016) that examine a variety of accretion disks and the effects of radiation and magnetic fields, along with their feedbacks to these systems were also studied and are worth mentioning. The studies of Fuerst & Wu (2004, 2007), Fuerst (2006) and Younsi et al. (2012) were also looked into thoroughly. In addition to the main accretion disk environment, we also chose to calculate the radiation force acting on material in outflow regions. We thus looked into original and fundamental accretion, jet and Blandford–Znajek process studies, such as Bondi (1952), Blandford & Znajek (1977), Komissarov (1999, 2001), Livio et al. (1999), Lee et al. (2000), Vlahakis & Königl (2004), McKinney (2005), Komissarov et al. (2007), Penna et al. (2013a), up to more recent general relativistic magnetohydrodynamic (GRMHD) simulation researches on the subjects, such as Nakamura et al. (2018), Parfrey et al. (2019), Park et al. (2019), Yuan et al. (2019), Event Horizon Telescope Collaboration IV (2019), Chatterjee et al. (2019), Mahlmann et al. (2020), Konoplya et al. (2021), Komissarov & Porth (2021).
One other important astrophysical effect related to the complex radiation field around a stellar black hole is the Cosmic Battery (Contopoulos & Kazanas 1998). This phenomenon operates in all environments including accretion disks, yet its impact is expected to be much more prominent when the central object is more compact, affected by the induced spacetime curvature (see Contopoulos & Kazanas 1998, Table 1). What happens is that the radiation emitted by the hot accretion disk is absorbed by the material itself, exerting on it the aforementioned radiation force. This force, however, is primarily acting upon the material electrons since f_{e}/f_{p} ∼ (m_{p}/m_{e})^{2}, where f_{p} and f_{e} the force on a proton and an electron and m_{p}, m_{e} their masses respectively. This results in the electrons moving with a different speed than the protons and hence in the generation of a ring current. This current leads consequently to the generation of a poloidal magnetic field that has notable consequences in the structure, equilibrium and evolution of the entire system and its possible outflows. This model was strongly criticized by BisnovatyiKogan et al. (2002), and was subsequently revisited in Contopoulos et al. (2006). Later on, applications were looked into, where notable effects were examined in relevant environments. Additionally, other topics were studied, such as the interaction of the Poynting – Robertson effect and the Cosmic Battery in Xray binaries (Kylafis et al. 2012) and the repositioning of the inner edge of the accretion disk due to radiation (Contopoulos & Papadopoulos 2012).
The primary objective of our research was to study the intensity and effects of radiation in stellar black hole environments and assemble information about systems where the aforementioned Cosmic Battery model could play a mentionable part or have noticeable impacts. Since, nevertheless, the research on supermassive black holes (SMBC) is much more extended and evolving in a much faster pace, particularly the past few years, we also considered that examining studies about the physics of SMBC systems would also be notably constructive, if not necessary. We have thus studied and compared, where possible, our work with the researches of Broderick & Loeb (2005, 2006a,b), Noble et al. (2007, 2011), Mościbrodzka et al. (2009, 2014, 2018) and Davelaar et al. (2018), all of which employ GRMHD.
The present work is a continuation of previous research in Koutsantoniou (2014) and Koutsantoniou & Contopoulos (2014), where we first investigated environments with noncentral radiation sources orbiting the central compact object. So far, our work has been the development of ray tracing codes along with custom postprocessing algorithms that allow us to redesign and improve the quality and speed of the codes without externally developed (i.e., “black box”) components. This was deemed necessary due to the size, duration and complexity of the subject.
The current paper presents and describes a leap forward regarding the quality and effectiveness of our codes as first presented in Koutsantoniou & Contopoulos (2014). Including a new process, we were able to increase the code’s resolution by more than a hundredfold, with only doubling the execution time. This allowed us to run a vastly increased number of simulations with finite optical depths at various heights, along and above the equatorial plane, inside and outside the accretion disk, compared to our previous work which involved only targets at the disk’s innermost stable circular orbit, hereafter ISCO. Finally, we should mention that the fully ray tracing profile of our codes allows us to look at these objects from very close (see also Davelaar et al. 2018). This permits us to study the systems and the incurring radiation forces, as well as obtain images of the black hole and the surrounding accretion disk for radii ranging from 1.25M for a rotating black hole, up to 25M or more, where M is the central black hole mass. This is important because as we see in the results, the closer we travel from the inner edge of the disk toward the black hole, the greater the radiation forces are and the faster their magnitude increases.
In this work, we present in Sect. 2 the mathematical formulation necessary to set up and use the Kerr metric and the locally nonrotating frames, and the methods to study particle trajectories and radiation effects. In Sect. 3 we present the various models of disks used in our work and the accompanying physics. In Sect. 4, we describe the five different families of codes written for our studies. Finally, in Sect. 5 we summarize our work and our codes, along with their possible extensions. In addition, we review the significance of this approach and mention the results we discuss in the second part of this work.
2. Mathematical formulation
We assume that the immediate environment around a rotating and accreting black hole, hereafter BH, can be adequately described using the Kerr metric. This suggests that the spacetime is determined by the central compact object that is axisymmetric, uncharged and possibly rotating. We also assume that the presence and motion of test particles does not affect the spacetime form or the stress–energy tensor. We hereafter use the geometrized unit system in which c = G = 1. We hence measure distances in units of gravitational radii r_{g} = GM/c^{2} = M. We also assume the Einstein notation for summation over double indices. Lastly, we denote spacetime components by Greek indices and space components by Latin indices.
2.1. The Kerr metric
The BH and the spacetime it creates, can be fully described using its mass M and spin parameter a. The Kerr metric in Boyer – Lindquist (hereafter BL) coordinates (t, ϕ, r, θ) is given by:
where:
with:
and the spacetime angular velocity is given by:
see Bardeen (1970) and Bardeen et al. (1972).
From the metric (2), we can determine the various characteristic surfaces present around a rotating BH. The event horizon arises from one of the poles of the g_{rr} component and is found at the outermost root of the equation Δ = 0:
The event horizon is thus a sphere of radius r_{evh} = 2M for a nonrotating Schwarzschild BH and r_{evh} = M for a maximally rotating one. The second characteristic surface is the static limit that constitutes the outer boundary surface of the ergosphere and can be found at the point where the g_{tt} component changes sign:
Studying the total energy of circular equatorial orbits (see e.g., Bardeen et al. 1972), we can see that there is a limiting case that describes particle orbits of infinite energy per unit rest mass. This is none other than the photon orbit, the innermost circular particle orbit. The radius r_{ph} of this photon ring is given by:
where the upper sign refers to direct and the lower sign to retrograde orbits. For a Schwarzschild BH with a = 0, the photon ring radius is r_{ph} = 3M, while for a maximally rotating BH with a = M, we have that r_{ph} = M for the direct and r_{ph} = 4M for the retrograde photon orbit.
Finally, another noteworthy set of trajectories are the equatorial circular orbits for massive particles and in particular the ISCO, whose radius is given by:
where:
where again the upper sign refers to direct and the lower sign to retrograde orbits. The ISCO starts from a value of r_{ISCO} = 6M for a = 0 and for a = M, reaches r_{ISCO} = M for a direct orbit. Before moving on, we remark here that the coincidence of the aforementioned characteristic surfaces for a maximally rotating BH at a radius r = M is deceptive and only an artifact of the BL coordinate system. These surfaces and orbits remain separate and distinct for a = M, as they differ in radial proper distance (e.g., Bardeen et al. 1972; Chandrasekhar 1983).
2.2. Locally nonrotating frames
The Kerr spacetime is stationary and axisymmetric but the central object rotation introduces complexity in both the physics of the problem and the mathematics required. First of all, the nondiagonality of the metric introduces cumbersome algebraic calculations when rising or lowering indices. In addition, physical difficulties arise when examining locations within the static limit and throughout the ergosphere. This is due to the fact that there cannot be static BL observers at these points, since the t basis vector becomes spacelike.
In order to simplify the calculations and remove various formulation problems inside the ergosphere, we choose to introduce a new set of observers and work in that new frame. The best choice of observers is one where said observers rotate with the spacetime geometry at the point we wish to study. We thus define the locally nonrotating frame (LNRF) or the zero angular momentum observer (ZAMO) at the point in question and describe the requested quantities using their projection on the Minkowskian orthonormal frame of the local observer. If required, we can then easily switch the calculated quantities from the LNRF into the BL frame. We denote quantities calculated in the LNRF by using hats over the component indices (e.g., ) and quantities calculated in the BL frame by unhatted indices (e.g., u^{α}). The transformation tensor and between the two frames has nonzero components:
Vectors u^{α} and , and tensors T^{αβ} and are transformed following the equations:
2.3. Particle trajectories
In order to study particle trajectories in Kerr spacetime, it is necessary to make full use of the particle’s integrals of motion and the respective conserved quantities. Let us assume a particle with rest mass m and fourmomentum p = (p_{t},p_{ϕ},p_{r},p_{θ}) in geodesic motion around a rotating uncharged BH. This particle has got four conserved quantities: the particle rest mass m, the total energy E = −p_{t}, the angular momentum component parallel to the rotation and symmetry axis L = p_{ϕ} and the Carter constant Q = p_{θ}^{2} + cos^{2}θ[a^{2}(μ^{2}−p_{t}^{2}) + p_{ϕ}^{2}/sin^{2}θ] (Carter 1968). This constant could perhaps be simply explained as a measure of how much a trajectory deviates from the equatorial plane. A particle in geodesic motion that starts in the equatorial plane and has Q = 0 will remain there indefinitely and a particle moving outside the equatorial plane with Q > 0 will at some point cross it. Let us note nonetheless, that the magnitude of Q does not relate linearly to the deviation from the equatorial plane motion.
The equations describing the particle motion are:
where the effective potentials are given by:
and λ is an affine parameter for massless particles and λ = τ/μ for massive particles, with τ the particle’s proper time (Bardeen et al. 1972; Wilkins 1972).
The above form of the equations is compact and elegant but hides various problems that appear when one attempts to solve them. The system appears problematic during numerical integration, since the square roots in the latter two Eq. (13) cause the quick accumulation of errors near the turning points. There are various solutions, such as reparameterization, in order to deal with this issue. In our study, we choose to work with the Hamiltonian and transform the above system accordingly. The new system of equations is then as follows:
We have therefore transformed the initial four equations of motion into a new system of eight differential equations. The new forms are smooth and do not have poles or other problems throughout their range and can be directly integrated. Let us note that the fifth and sixth of the above equations describe two of the motion’s conserved quantities, the conservation of energy and zmomentum respectively.
Finally, we define the coordinate angular velocity for a circular equatorial orbit as:
where u = (u^{t},u^{ϕ},u^{r},u^{θ}) is the fourvelocity of a particle.
2.4. Radiation and equations of motion
Our main goal in this subsection is to calculate the effects of radiation on the target particle dynamics. Thus, we begin from the formula that relates the target particle position with the acceleration a^{α}:
where x^{α} are the particle position components, τ the proper time and the Christoffel symbols or connection coefficients (Mueller & Grave 2009).
The acceleration in turn, can be given by the relativistic equation of motion:
where m is the rest mass of the target particle and f^{α} are the nongravitational fourforce components. In the present study this is the radiation fourforce on the target particle motion.
In order to calculate f^{α}, we need to know the flux of the radiation that generates it, as:
where σ is the particle cross section for the momentum transfer. The radiation flux fourvector F^{α} can in turn be calculated using the target particle covariant fourvelocity u_{μ} and the radiation stress – energy tensor T^{αβ} using the formula:
where is the projection tensor:
In order afterwards to acquire the BL radiation stressenergy tensor T^{αβ}, it is necessary to find the LNRF stressenergy tensor and then use Eq. (12):
where the are given by Eq. (11).
In order to calculate now, we make use of the formula:
where I(r,θ,ã,b̃) is the frequency integrated specific intensity of the radiation, the solid angle element with ã and b̃ the related appropriate local angles and a unit spacelike vector (Fig. 1). Simple calculations can give the vector components as:
Fig. 1.
Local sky around the target particle. Left: forward in time (the photon moves toward the target). Right: moving backwards in time (the photon moves away from the target). For the incoming photon, we define angles ã and b̃ similar to the typical polar angle θ and azimuthal angle ϕ of spherical coordinate systems. The disk in stripes at the top left represents the BH event horizon. 
Intensity I is, as expected, a function of the particle’s position in space, since different locations receive different amounts of radiation. From this, we have already excluded the ϕ coordinate due to the spacetime axisymmetry. Additionally, I also depends on the angles ã and b̃, since different amounts of radiation are received in different orientations of the local sky.
Finally, the frequency integrated specific intensity I is calculated by integrating the specific intensity I_{ν} across all the contributing frequencies:
This can be applied to any distribution such as a blackbody, thermal radiation, a singleenergy beam of light, or frequencyindependent radiation. In an environment where the radiation is emitted by just a surface layer of the source, the calculation of I is relatively easy and it is done with the method described in Sect. 2.6. On the contrary, in environments where the radiation is emitted by multiple layers or various objects, that method cannot be used. In order to calculate the specific intensity I_{ν} there, it is required to investigate the radiation transfer process and solve the appropriate equation. This will be addressed in the following subsection.
2.5. Radiative transfer
In this subsection we look into disks with finite optical depth, where photons are emitted by the material throughout their entire volume. The radiation is then regulated by its passage through the disk’s absorbing and emitting material. If this material is dense enough, then the ray reaching a point deep inside the disk accumulates a high enough optical depth. This means that the ray cannot have originated from outside the local disk material, but only from inside it. If, on the contrary, the material is of low density or of small quantity, then the ray will travel outside the local material. Later on, the ray can either reenter the emitting matter further along its path or escape to infinity. In these cases, the incoming ray will only bring in low intensity radiation. 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 a thermalized material of number density n. This material consists of particles that act as radiation absorbers. We define the absorption coefficient a_{ν}(cm^{−1}) at frequency ν as:
where σ_{ν}(cm^{2}) is the absorbing area cross section at a particular frequency. Assuming an initial specific intensity I_{ν(}erg cm^{−2} s^{−1} ster^{−1} Hz^{−1}) at frequency ν, the presence of the material’s radiation absorbing particles for a propagation length ds, will cause a decrease in this specific intensity of a propagating light ray given by:
Things are simpler for the emission coefficient j_{ν}(erg cm^{−3} s^{−1} ster^{−1} Hz^{−1}). When the light ray propagates for distance ds, it transverses emitting material of volume dV and its specific intensity increases as:
The radiative transfer equation combines the above two processes and describes the resulting effects on the light ray’s specific intensity as:
In order to express the solution of the above equation more elegantly, we introduce the concept of the optical depth τ_{ν} at frequency ν that is defined as:
We can then calculate the optical depth by integrating the above along the path of the light ray:
where s_{0} is an arbitrarily selected initial point of the scale. We can subsequently restate the RTE as:
Integrating this gives the solution to the radiative transfer equation:
The above magnitudes are not in most cases Lorentz invariant and thus cannot be used in the general solution of the various problems, unless restated in such a form (Misner et al. 1973). We begin by considering the phase space number density:
where N is the number of particles under examination and 𝒱 the phase space volume they occupy. By taking into account Liouville’s theorem in curved spacetime, we have that:
where again, λ is an affine parameter for massless particles. By combining the above with the conservation of particle number along the world line of the bundle, we obtain that:
which is the collisionless Boltzmann kinetic equation and hence 𝔑 is Lorentz invariant. The phase space volume is:
and hence:
Here, h is Planck’s constant. Since the specific intensity is defined as:
we can see that:
and therefore the Lorentz invariant specific intensity ℐ_{ν} is:
The optical depth, used to count photon fractions, is a scalar quantity and is thus invariant:
In order to find the Lorentz invariant absorption coefficient, we use Fig. 2. The tube width d = d ′ is the same in both the lab and the matter rest frame, since it is perpendicular to the direction of motion. Likewise, the xcomponent of the photon momentum k_{x} = k_{x}^{′}, remains unchanged. This subsequently means that k sin θ = k ′ sin θ ′ and thus ν sin θ = ν ′sin θ ′. From the Lorentz invariance of the optical depth τ_{ν} = a_{ν}s, we then have:
Fig. 2.
Left: emitting matter in the lab frame K, moving with velocity along the vertical axis in a tube of width d. A photon of frequency ν crosses the tube, traveling at an angle θ from the vertical axis. Right: emitting matter in its rest frame K′ in the tube. The photon with frequency ν ′ now appears to cross the tube at angle θ ′. 
Finally, for the emission coefficient we utilize Eqs. (32), (41) and (43) and conclude that:
The Lorentz invariant form of the RTE (32) will therefore be:
Since for the optical depth, it is dτ_{ν} = a_{ν}ds, Eq. (45) along with (43) and (44) gives:
In order to improve this, we also implicate the path length variation ds/dλ. By using the projection tensor h^{αβ} = g^{αβ} + u^{α}u^{β}, we have the photon velocity (υ ′)^{α} in the fluid frame K′ as:
where u^{α} is the fluid fourvelocity. By the above, we obtain that:
and for the frequency ratio, it is:
(see also Younsi et al. 2012). Quantities with an accent, such as ν ′ above, are henceforth measured in the local rest frame. Combining these, we have that:
Equation (46) combined with (43), (44) and (50) gives the differential form of the invariant RTE, the general relativistic radiative transfer equation (GRRTE) equation as:
Integration of the above gives the solution for the Lorentz invariant specific intensity:
The optical depth can be calculated as:
and Eq. (52) can then be rewritten as:
2.6. Intensity of single emission source radiation
In this subsection we describe the way to estimate the radiation 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, henceforth AD, 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 studied in the literature: Abramowicz et al. (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 nonrotating and rotating masses and radiating sources. We have also studied this subject in the previous work Koutsantoniou & Contopoulos (2014) examining fewer examples of totally opaque disks and a single observer position for each modelBH spin set. Here, apart from expanding into semiopaque disks discussed later on, we expand our analysis to more disk models and instead of having a single observer at the ISCO of each modelspin set, we fill the entire region of the system with a large amount of observers in different locations.
As we saw previously in Eq. (41), the Lorentz invariant specific intensity is I_{ν} = I_{ν}/ν^{3} and thus for the frequency integrated specific intensity it is:
for any two random points P_{1} and P_{2}. From this, we have that for the emitted and the received frequency integrated specific intensity of a photon, I_{em} and I_{rec} respectively, it is:
We note here that the frequency fraction ν_{rec}/ν_{em} that appears above does not depend on the frequencies involved, but only on the spacetime and the photon’s emission angle. This frequency fraction includes the effects of three different phenomena caused mainly by the spacetime properties. Firstly, it includes the effects of gravitational time dilation, which appear both in Schwarzschild and Kerr spacetimes. It also includes the frame dragging frequency shift due to the spacetime’s differential rotation that appears only in a Kerr spacetime. Finally, it includes the Doppler shift caused by the motion of the source’s emitting surface. This can exist in both Schwarzschild and Kerr spacetimes.
The gravitational time dilation causes a frequency shift for the received frequency ν_{rec} given by:
where ν_{em} the emitted frequency.
Assuming the photon source moves azimuthally with negligible radial and poloidal velocity components, the Doppler shift due to the emitting surface motion then introduces a change in frequency:
where the source threevelocity here, γ = (1−V^{2})^{−1/2 } the emitting material Lorentz factor and ψ the angle between the emitting matter velocity and the photon emission direction (Fig. 3). Let us mark that both the γ factor and the ψ angle are measured in the ZAMO frame at the point of emission.
Fig. 3.
Emission of a photon from the hot AD. The photon is emitted along n̂ from an element of matter moving with velocity . The photon is emitted at an angle ψ relative to the accreting material motion. 
Concluding the necessary transformations, we have the factor required for the implementation of the frame dragging effects, which is:
where k_{a} are the photon covariant fourmomentum components, which are also conserved quantities. The ratio k_{ϕ}/k_{t} depends only on the direction of the photon emission and from the previous statement is also a conserved quantity. Combining the above, we have for the received frequency that:
and therefore for the frequency integrated specific intensity:
In Fig. 4, we can see a breakdown of the process described above, where for visual simplicity we assume that the emission of photons is done by a central object.
Fig. 4.
Schematic of the frequency changes. A Doppler shift is used to move from the frame comoving with the emitting surface (left) to the LNRF at the radius of the surface (middle), taking into account the emission source rotation. Then we move to the receiving particle frame (right), 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. 
3. Accretion tori
In this section we present the accretion tori we used for our codes. Some of the tori are optically thick while others are stratified and semiopaque. In general, we examine tori of assorted geometrical shapes and diverse density profiles. This way we can better cover for example the various stages of matter infall into the BH and the different stages of disk evolution in Xray binaries. For the physics, evolution, magnetic fields role and various other phenomena present in the intriguing Xray binary environments, one could consult Esin et al. (1997), Verbunt (1999), Bildsten & Rutledge (2001), Fender (2002), Haggard et al. (2004), Müller (2004), Done et al. (2007), MeyerHofmeister et al. (2009), van Haaften et al. (2012), Heinke et al. (2013).
A broad and thorough study of compact objects, accretion disks and the many physical phenomena observable in such environments was presented by Blandford et al. (2002). In order for us to design and choose acceptable disk models with a good balance between model quality and computational time, we looked into Shakura & Sunyaev (1973), Cunningham (1975, 1976), Abramowicz et al. (1978, 1988, 1996), Kozlowski et al. (1978), Narayan & Yi (1994, 1995), Lasota (1999), Igumenshchev et al. (2003), Narayan et al. (2003, 2012), Narayan & McClintock (2008), Sądowski (2009), Noble et al. (2011), Penna et al. (2013b), Fuerst (2006) and references therein, Sadowski (2011) and references therein. We have thus constructed two different AD groups, one with optically thick disks and one with semiopaque or translucent disks. The main repercussions of the disk’s thermal radiation, along with its impact onto the disk’s geometry, vertical height and optical thickness were described and studied in Thorne & Price (1975), Inoue & Hoshi (1987), Takahashi et al. (1995), Beloborodov (1998, 1999, 2001), BisnovatyiKogan (2001), Abramowicz & Fragile (2013). For the finer points that differentiate the optically and geometrically thin and thick disks and segregate the categories, one could look for example into Artemova et al. (1996), Quataert (2001) and Dubus (2003).
We note here, that the ADs we considered for this work are just a sample for the study of the most commonly considered models. Many other frequently used AD models, such as Novikov–Thorne disks (Novikov & Thorne 1973) for example, could also be implemented and studied using our codes, should the need arise. The reason why we chose to consider more simplistic perhaps disk models than the latest GRMHD researches (e.g., Chatterjee et al. 2019; Mahlmann et al. 2020), is because of the number of code executions this work required. This results in any model improvement, such as a more realistic, a nonaxisymmetric or a timevariant AD, greatly increasing the total execution time.
The first group of disks we examine is used to represent tori that increase their density abruptly and very close to their outer surface. They are mostly rotationally supported and totally optically thick. This practically means that there is no reason to solve the GRRTE for these tori and instead another methods of calculation must be used. We use these models to describe physical tori that are either cold or compact, or both. An interesting case they could also be used to describe, is transient stages of the Xray binary systems (see e.g., Tauris & van den Heuvel 2006). During the quiescent stages of these systems, their ADs tend to be cooler and at a larger distance from the compact object (e.g., Esin et al. 1997; Narayan & McClintock 2008). They subsequently remain in a similar state for an indefinite amount of time. At some point later on, they start increasing their temperature and swelling up while reducing their density and density gradient. From that point on, we can no longer describe them using these opaque models and must instead employ semiopaque tori models. Additional and detailed information about the evolution and stages of Xray binaries can be found in Tauris et al. (2000), Podsiadlowski et al. (2002) and Chen & Podsiadlowski (2016).
The second group of tori, the semiopaque ones, describes more common and familiar perhaps cases of ADs. There is a measurable density and temperature gradient. For these cases, we must use a ray tracing process and solve the GRRTE along the photon trajectory. This way, we calculate how much radiation is produced by the hot material in every step and how much of this is absorbed away by it. Depending on the direction and angle of motion of the traveling photon, it can at times be absorbed by the disk material and, at other times, it can traverse part of the disk without it being absorbed. This means that at some points the disk is optically thick and at other times optically thin, hence the name semiopaque. In these cases, we observe effects such as transparency and limb darkening (Fig. 19, bottom left and right).
All the tori we mention here can rotate in various ways. In our program executions we assume that the disk material can rotate circularly (u^{r} = 0) with the typical coordinate angular velocity Ω (Eq. (16)) or with more detailed profiles such as the one in Eq. (75) below. Additionally, we examine cases where the disk material follows inspiral motion profiles (u^{r} ≠ 0) attempting to mimic the SANE (Standard And Normal Evolution) and MAD (Magnetically Arrested Disk) models (Narayan et al. 2003, 2012; Penna et al. 2013b). In each execution run, our codes give results for all of the aforementioned different velocity profiles for the AD material.
We should also mention here that although in many of the AD simulations and studies it is assumed that the material is in a stationary condition, geometrically at least, in reality it is far from that. There are increased amounts of turbulence, instabilities (see e.g., Tchekhovskoy et al. 2011; McKinney et al. 2012; Narayan et al. 2012) and other phenomena taking place, often in smaller scales, that are at times ignored. One such phenomenon is the flow and diffusion of angular momentum throughout the different disk sectors. The main cause of this diffusion is considered to be the material viscosity and is treated in various ways. One of the best known and more frequently used methods is the αviscosity approach (Shakura & Sunyaev 1973; Abramowicz et al. 1988). This method, nevertheless, cannot be applied to all disk models, such as non αdisks, where other solutions must be found.
Another important phenomenon that is known and generally mentioned in such works but often ultimately ignored, is the existence of magnetic fields. The presence of magnetic fields usually gives rise to very important phenomena, such as the Blandford–Znajek process (Blandford & Znajek 1977; Livio et al. 1999; Komissarov 2001; McKinney 2005; Komissarov et al. 2007; Penna et al. 2013a), that affect the structure of the disk and its stability, and determine its evolution. In theory, the BlandordZnajek process is, along with the Penrose process (Penrose & Floyd 1971), one of the two most promising phenomena responsible for launching astrophysical jets. The presence of substantial poloidal magnetic field makes the extraction of spin energy and angular momentum from a rotating BH possible. This happens due to the escape of angular momentum from the rotating magnetosphere inside the ergosphere. This mechanism can accurately describe the formation and ejection of jets from spinning SMBHs and is also considered to play a pivotal part in gammaray bursts.
An additional noteworthy effect brought on by the existence of magnetic fields and very important for the AD dynamics, is the generation of magnetorotational instabilities (MRI). This strong fluid instability occurs when a conductive AD is situated in a magnetic field and is rotating differentially, with its inner regions rotating faster than its outer regions. The freely moving charges of the material are subjected to the Lorentz force, due to the presence of the magnetic field. Any fluid element deviating even to a small extent from circular motion, has its trajectory further destabilized by a force increasing proportionally to the displacement from the circular orbit. This causes the disk to become unstable and consequently turbulent. The MRI can therefore have important consequences, particularly on the distribution and flow of angular momentum throughout the AD and its diffusion toward the outer layers and components (see Balbus & Hawley 1991, 1992, 1998; Hawley & Balbus 1991, 1992; Hawley et al. 1995; Krolik 1999a,b; Turner et al. 2003; Pessah et al. 2007). Other phenomena the MRI is also expected to influence is the formation of active galactic nuclei (Krolik 1999a), the production of Xrays in compact object systems (Blaes 2004), as well as gammaray bursts (Wheeler 2004).
3.1. Optically thick accretion tori
In this subsection we describe the models we used for opaque tori. We built some of these tori by assuming simplistic disk cross section shapes, such as polygons. These tori can be viewed either as toy models or as initial condition “snapshot” states. One could then go on to study the evolution of these tori taking into account the presence of radiation effects. Some of the other models we considered are more complex. We built those models selfconsistently by assuming that the material of the disk is supported and kept in place by its rotation.
Optically thick ADs are generally expected to be geometrically thin (e.g., Shakura & Sunyaev 1973), even though that is not always a canon, as described in the aforementioned accretion disk studies. This is caused by the “inefficiency” of the radiation: since the disk material is opaque, the radiation transmitted by its hot components cannot reach other, more distant parts of the disk. This results in each local material component to have a significantly lesser “inflating” radiation pressure element than a “deflating” gravitational force element. The result is an AD of smaller geometrical thickness and a much larger pressure gradient, specifically close to its outer surface. We, however, investigate both the cases of geometrically thin and thick opaque ADs.
A matter of particular importance is the surface temperature distribution we assume for these tori. For our calculations, we considered two separate cases, an isothermal disk and a disk whose temperature follows T ∝ r^{−3/4} (Shakura & Sunyaev 1973). The first case, albeit unnatural, is the simplest possible one could imagine and is thus perhaps easier to understand and effortlessly anticipate certain results.
The second case is to assume that the disk temperature distribution is caused by the material accreted onto the central compact object. If we assume that the object’s luminosity is equal to the Eddington luminosity L_{Edd} = 4πGℳm_{p}c/σ_{T}, then we have an Eddington accretion rate:
where G is the gravitational constant, ℳ the disk mass, and m_{p} the mass of the proton. Assuming a large enough amount of scatterings, the disk material can be adequately described by blackbody radiation. Then, its temperature will be given by:
where σ_{SB} is the Stefan – Boltzmann constant (see Longair 2011).
We notice here, nevertheless, the problem that arises if we simplistically hypothesize the above. If we assume that the accretion luminosity is equal to the Eddington luminosity, then the disk cannot be geometrically thin. This is because, as the accretion luminosity increases, the radiation pressure exerted onto the material keeps getting larger and finally comparable to local gravitational forces. As this happens, the disk keeps inflating by gaining height and width and thus gradually turning into a geometrically thick and optically thin torus (e.g., Thorne & Price 1975). The easiest way to bypass such problems is to assume that the accreting object radiates only a fraction ϵ of the Eddington luminosity:
and thus, it is:
Finally, after having picked any of the disk temperature profiles, we can have its effect on the emitted photon from:
In our work, we have considered so far six different models for optically thick accretion tori. Specified by their given names, we have the models band, disk, slab, wedge, torus and opaque rotationally supported torus:
(a) Band (toy, snapshot model): a cylindrical surface of halfheight h (from its highest point to the equatorial plane) at the distance of the respective ISCO for the selected spin parameter (Fig. 5). We can freely choose the halfheight h without restrictions and we can also adjust the cylinder radius. We use this mode when we wish to study for example the radiation effects on only the innermost surface of an AD.
Fig. 5.
ADs “Band” (left) and “Disk” (right) constituted by twodimensional surfaces (Sect. 3.1a and b respectively). In the center of each image is in black the BH event horizon and around it in gray, its ergosphere. We see the AD in the yellow (hotter) and red tones (colder). 
(b) Disk (toy, snapshot model): an infinitesimally thin disk at the equatorial plane (Fig. 5). Its inner radius is at the distance of the ISCO and its outer radius at any multiple of this distance. The innermost and outermost radius of the disk can easily be adjusted.
(c) Slab (toy, snapshot model): a disk of halfheight h from the radius of the ISCO to a distance of three times the ISCO radius. The cross section of the disk is a rectangle (Fig. 6). The halfheight h can be freely chosen without restrictions. The innermost and outermost radius of the AD can also be adjusted.
(d) Wedge (toy, snapshot model): a disk whose cross section is an isosceles trapezoid and is centered above and below the equatorial plane. Its inner radius is equal to the radius of the ISCO and its outer radius to three times that. We construct the disk in such a way that an angle with its vertex at the BH (the origin) extending outward, reaches the ISCO cylinder intersecting a ring of halfheight h. The angle sides continue extending outward in the same direction until crossing the outer edge of the disk (Fig. 6). The disk halfheight h, as well as its inner and outer radius can easily be modified.
(e) Torus (toy, snapshot model): a disk with a circular cross section. The center of the circle is at coordinates (r,θ) = (2r_{ISCO},π/2) and its radius is equal to the ISCO radius. The disk inner edge is therefore at r = r_{ISCO} and the outer edge at r = 3r_{ISCO}. The cross section center and radius of the disk can be adjusted at will (Fig. 6).
(f) Opaque rotationally supported torus (ORST – selfconsistent model): a rotationally supported torus. This disk is one of the more complex cases considered for optically thick disk examples. The disk we consider here is stationary and axisymmetric and has its rotation axis aligned with the rotation axis of the BH. In our work we have assumed that the two angular velocity vectors are collinear, but it is simple to consider the opposite case in order to study retrograde disks. We then assume that the disk acceleration is what creates this setup and specifies its shape. The acceleration along the particle trajectory is given by:
where are the Christoffel symbols which we can calculate from the metric (2), using the formula:
(more details can be found in ChoquetBruhat et al. 1977). The torus we have assumed here, like many of the tori in works of this type, has negligible radial and poloidal velocity components, so (67) can be slightly simplified. We then get for the acceleration components that:
The first two of the above equations are in accordance with our initial assumption that the disk is stationary and axisymmetric respectively. The last two equations respectively are what can give the surface of constant acceleration, the isobaric surfaces through:
Using the above equations, this gives:
This leads to the pair of differential equations:
where:
and Ω = u^{ϕ}/u^{t} is, as before, the material’s angular velocity. The last information necessary to solve the above and have the resulting torus is its inner edge at the equator. We find this by solving the equation for marginal stability orbits:
In order to solve Eq. (74), we must define the angular velocity function. We use here the angular velocity profile proposed and explained in Fuerst & Wu (2004, 2007) and Younsi et al. (2012):
where r_{K} is the equatorial plane radius at which the material moves with Keplerian velocity. Here, the parameter n_{p} corresponds to pressure forces and is responsible for the geometry of the torus determining its thickness. Tori solutions for various r_{K} and n_{p} values are shown in Fig. 7, while the selected tori used in our simulations are displayed in Fig. 6.
Fig. 7.
ORST cross sections (Sect. 3.1f). In (a) we can see the effects caused by the change of the radius of Keplerian rotation speed r_{K}. In (b) we see the tori shapes and sizes for various values of the parameter n_{p}. 
The last thing remaining for ADs of this kind is to calculate the frequency integrated specific intensity. This is done by using the method described in Sect. 2.6 and then applying Sect. 2.4 to determine the radiation flux or force and the ensuing acceleration caused by the disk’s hot material.
3.2. Semiopaque tori
In this subsection we refer to the semiopaque and transparent disk models considered in our work. Again as before, some of the models are more simplistic than others (toy or snapshot models) and some are based on specific physical conditions and are more complex (selfconsistent models).
For our research, we considered up to this point five disk models that could fit in the semiopaque or transparent category. Specified by the name of the considered models, we have the following cases:
(a) No torus: there is no AD around the central BH. Photon trajectories continue until they either cross the event horizon or escape the system by crossing an adjustable outer radius boundary. We use this mode when we wish to study just geodesics for example.
(b) Semiopaque pressure supported polish doughnut (PS PD – selfconsistent model): a stationary and axisymmetric radiation pressure supported polish doughnut. We construct this accretion torus following Abramowicz et al. (1978) and Kozlowski et al. (1978). As in previous tori, we assume here that the material has no significant radial or poloidal velocity components. The material thus has fourvelocity u^{α} = (u^{t},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,Γ = 4/3 and , with ℏ the reduced Planck constant, μ the mean molecular weight and β the ratio of gas pressure to total pressure. We also assume that the disk rotation follows (75). The AD is then described by:
where ξ(r,θ) = ln[Γ−1+Γκn(r,θ)^{Γ − 1}] is a function of the disk number density n(r,θ) and the acceleration components are a_{α}(r,θ) = u_{α}; β(r,θ)u(r,θ)^{β}. The maximum number density point lies in the equatorial plane at r_{center} = r_{K} and the number density there is n_{center} = 10^{18} cm^{−3}. Cross sections of the tori number density used are shown in Fig. 8.
Fig. 8.
Polish doughnut number density cross sections (Sect. 3.2b,c). The center of all tori lies at (12,0) and has a number density of 10^{18} cm^{−3}. The top row shows the disk for BH spin parameters a = 0 and a = 0.5, while the bottom row for a = 0.9 and a = 0.998. 
(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) Semiopaque LFM torus (toy, snapshot model): a stationary and axisymmetric semiopaque torus of circular cross section. The cross section center lies on the equatorial plane at r_{center} = 2r_{ISCO} and the torus cross section has a radius r_{torus} = r_{ISCO}. The torus thus stretches from an inner radius of r_{inner} = r_{ISCO} to an outer radius r_{outer} = 3r_{ISCO} and has a maximum height h_{torus} = r_{torus} = r_{ISCO}. The center number density is n_{center} = 10^{18} cm^{−3} and decreases to zero moving toward the torus surface. For this torus the product a_{ν} ⋅ 2r_{torus} is ∼1 − 5 throughout the cross section. Images of the cross section number density are shown in Fig. 9.
Fig. 9.
LFM model number density cross section for any spin parameter (Sect. 3.2d,e). The center of the tori lies at (2r_{ISCO},0) and has a number density of 10^{18} cm^{−3}. 
(e) Translucent LFM torus (toy, snapshot model): a translucent LFM torus of circular cross section. It is the same as the previous model, but without its material absorbing any of the photons crossing it.
The disk models we discussed above are responsible for giving us most importantly the material’s number density n(r,θ). From that, we can then obtain other useful quantities for the matter, one of which is the material temperature. Following standard procedures, we have here:
where ρ the (volumetric mass) density.
Continuing on, we can obtain firstly the necessary material’s absorption coefficient from Eq. (26) as:
where σ_{ν} is the absorption cross section best chosen for the processes under study. Further on, we find the emission coefficient of the material by using the thermal emission and blackbody radiation properties we have for the assumed disk. The thermal emission assumption dictates that the emission coefficient will be given by:
where B_{ν}(T) the Planck function and T the corresponding temperature. In order to procure the Planck function, we make use of the blackbody attributes of the material and have:
where k the Boltzmann constant. Combining then the above equations, we have for j_{ν} that:
We additionally remark here that in the above equations, we could also add shaping functions to modify the emission and the absorption of the material. We could this way study other disk models or perhaps different physical properties. A more general form of the above functions could thus be:
where C_{abs} and C_{em} are absorption and emission coefficients respectively and f a shaping function of the number density n(r,θ), the material temperature T[n(r,θ)] and the photon energy function E[n(r,θ)]. We can see one such case for example in Younsi et al. (2012).
After concluding the calculations described above, we have the resulting radiation intensity I(r,θ,ã,b̃). We can then apply the method of Sect. 2.4, obtaining the stress – energy tensor, the flux and the force of the radiation.
4. Algorithms and codes
In this section we present the codes we developed from 2012 to 2021, used in our work and kept improved and cross checked since. We explain their capabilities and show some of their results. We also note here that all of the codes used and presented, were designed in order to be executed with extremely limited computational resources. This has important consequences on the design, speed and effectiveness required from the codes.
4.1. Code Omega
Code Omega was the first we created in our work and it is a central part of all following codes, as its main purpose is to calculate photon trajectories. It works for a Schwarzschild and a Kerr spacetime, but it can be easily modified in order to work in other spacetime models as well, for example Kerr – Newman, Reissner – Nordström, Friedmann – Lemaître – Robertson – Walker etc.
The code studies the tori models referred to and explained in Sects. 3.1 and 3.2 and depicted in the accompanying Figs. 5–9. It is also possible to add new AD models to the code, such as disks resulting from combinations of the previously mentioned models. Examples of such tasks are shown in Fig. 10. Additionally to the above, Omega has an option “Star”, where instead of a noncentral accretion disk, it studies the radiation field produced by a central spherical star of isotropic and uniform radiation. We have mostly used this option so far in order to crosscheck and validate the veracity of our results in comparison to other studies, as we see further on. This code feature can also be modified as well, in order to study objects more complex, such as nonspherical or nonuniformly emitting stars.
Fig. 10.
ADs created from combinations of previously referred models (from top to bottom): wedge & slab, polish doughnut & torus, torus & wedge and polish doughnut & spheres. 
Omega solves the particle trajectory equations mentioned in Sect. 2.3 for a photon and finds the trajectory and the point of origin of said photon. This point of origin could be on the hot AD, the BH event horizon, or a location outside and far away from the system. If all that the backward photon trajectory intersects with is the event horizon or the system exterior, then no radiation or energy is carried to the AD material and the target particle.
Depending on the environment in which we use it, Omega can have different outputs. In its original form, which is visual, the program has an interface that allows the user to select primary properties for the environment, such as the disk model, the BH spin parameter and the disk height. Also, the user can select important options for a trajectory, including its maximum length, its point of origin and angle of emission and the two emission angles ã and b̃. Finally, there are some additional visualization options that include the choice of frame size of the visual box and the depiction of obscured parts of the outer disk. The code’s dynamic output picture shows the BH event horizon and its ergosphere, the AD and the requested photon trajectory. The particle trajectory is drawn in different styles and colors for escaping photons, photons infalling in the BH and photon trajectories starting from the AD. In addition, some trajectory information are displayed in the picture, including the photon energy and angular momentum, the trajectory’s Carter constant and a confirmation of the photon momentum magnitude conservation. All the above can be seen in Fig. 11.
Fig. 11.
Control interface of the Omega code, as described in Sect. 4.1. In light green is the requested photon trajectory and with the cyan point we denote its point of origin. 
The latest Omega versions are the most compact, since they now are functions and parts of more complex codes, outputting thus only numerical data and no visual information. Omega runs for a single photon trajectory per execution and returns key information for this trajectory. Firstly, it reveals the existence or absence of incoming radiation in the requested direction. Additionally, it gives the precise photon trajectory numerically as well as the coordinates of the emission source.
In order to verify the validity of the calculated particle trajectories, we subjected Omega to various tests. Firstly, we verified that Omega can successfully predict elementary yet fundamental particle trajectories, such as the notable orbits presented in Bardeen et al. (1972) and discussed in Sect. 2.1. In Fig. 12 we see that we indeed obtain such orbits, like the photon ring (massless particles) and the ISCO (massive particle). Additionally, we confirm that if the BH mass is reduced to zero, the spacetime is no longer a curved Schwarzschild or Kerr one. Instead, it degenerates into a flat Minkowskian form, causing the particle orbits to turn into straight lines.
Fig. 12.
Tests for the validity of Omega. We confirm that the code successfully produces fundamental particle orbits, such as the ones mentioned in Sect. 2.1. In each picture we have the BH situated at the axes origin and with gray is its event horizon. In (a) we have in red the photon ring produced numerically by the code and on top of it, in the black dashed circle, its theoretically expected position. In (b) we see in blue the ISCO produced numerically and on top of it, in the black dotted circle, its theoretically expected position. In (c) we set the BH mass M = 0 (with light gray showing where its horizon would have been) and examine the trajectory of a photon. With any BH mass M > 0, the photon would have to infall and cross the event horizon, but with no BH mass present, the spacetime is Minkowski and the photon travels in a straight line. 
Continuing on, we compare our code results to the work of Bini et al. (2015) concerning photon trajectories and we note agreement. In Fig. 13 we see for example free photon trajectories produced by Omega. The photons are launched perpendicularly and upward from various points in the equatorial plane of a Schwarzschild BH (we compare to Bini et al. 2015, Fig. 1).
Fig. 13.
Free photon trajectories in the poloidal x − z plane. Photons are emitted perpendicularly upward from various points of the equatorial plane of the system. The closest to the BH emission points are (3, 0) and (−3, 0) and taking into account the angle of emission, these photons should and indeed do follow the photon ring trajectory (compare to Bini et al. 2015, Fig. 1). 
Finally, we examine massive particle trajectories and compare Omega results to the work of Levin & PerezGiz (2008) and Levin & Grossman (2009), with which we are also in accord. In Fig. 14, for example, we crosscheck our code results with what is shown in Levin & PerezGiz (2008) and we again notice agreement.
Fig. 14.
Omega code successfully creating massive particle orbits. These trajectories are in agreement with the corresponding orbits shown in Levin & PerezGiz (2008). 
4.2. Code Infinity
The Infinity code is now an extended and fully automatic program that can run complete simulations for a large stream of different points and situations. These runs do not require to be for the same object or of the same environment. It automatically creates output files and saves numerical data, images and information possibly useful for other runs or aggregated results. It also has features to prevent the loss of data and computational time during power or network outages.
In the Infinity code, we choose a point of interest anywhere in the system under study and the program sets an observer there. It then scans the entire sky around this observer, solving the GRRTE and returning the radiation reaching the observer from each of the angles of the local sky. It then calculates the radiation stress – energy tensor and creates various images and sky maps of the radiation. It continues on to run a procedure that increases the code resolution by a significant number of times (see Fig. 15).
Fig. 15.
Two instances of the first version of the resolution enhancement process (625 times here). In the top row of each part is the original image created by the AD observation. On the bottom is the resulting image after the resolution enhancement. These functions are now rewritten and incorporated into Infinity to enhance results by a number of times. 
The program continues on to calculate the radiation flux and the force fourvectors applied to various observers. This includes observers at rest in the local frame rotating with ω = ω(r,θ) (Eq. 5) as seen from infinity, observers in circular orbits moving with Ω (Eqs. (16) and (75)) depending on the selected disk model and observers on accreting material approximating the SANE and MAD models. In addition to the above, we can also have observers in two possible outflow regions, “jet regions” if you will, close to the system rotation axis and above a certain height (Fig. 16). Let us note here that, even though we have made our choice on certain velocity profiles for the various disk models, this can easily be adjusted to suit the needs of other disk models with different material velocities.
Fig. 16.
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 AD 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 BH’s ergosphere and has thus a radius r = 2M. The actual region is considered to exist above a certain height, further away from the BH. In this region, we can have matter at rest in the local rest frame and matter flowing outward with a certain velocity. We may also consider a subregion there (see e.g., Asada et al. 2016; Park et al. 2019), a narrower cone or cylinder with radius r = M with stronger outflow and faster velocity or even with velocity of the opposite direction. 
The program then concludes by outputting its results and creating various save files. These include a Mollweide map picture (Fig. 17) for the radiation in the specific observer location, pictures of the important matrices of the simulation and tables of the assorted profiles’ force components and the radiation stress – energy tensor matrix. Additionally to the figures, videos of flights around and if possible through the disk for all the models we studied, can be found on youtube.com, under the name of this work’s creator “Leela Elpida Koutsantoniou”^{1}.
Fig. 17.
Infinity code: Mollweide sky maps of the frequency integrated specific intensity for tori with downward angular velocity vector. The color scale for each image is displayed below it, from minimum to maximum. Top: wedge AD model (opaque, Sect. 3.1d), bottom: LFM AD model (semiopaque, Sect. 3.2d). 
In order to certify the legitimacy of the code, we crosschecked our results with other works. We start by assuming that the radiation source is a central nonrotating star that emits photons isotropically from each point of its surface. We then assume that the star can have various radii and set observers at assorted distances from its surface. We then use Infinity and observe the system, noting the apparent size of the central object and the ensuing radiation stress–energy tensor. Finally, we look into our results and compare them to their expected values, given by analytical formulae in Abramowicz et al. (1990). As we can see in Table 1, we are in good agreement with the theoretically expected values. If, however, we desire so, we can have even better estimations by increasing the code resolution.
Comparison of Infinity results with formulae in Schwarzschild spacetime.
Keeping then in mind that in general the radiation source is not expected to be static but instead rotating, we expand the above check and look into environments around rotating stars. We thus look into rotating stars in Kerr spacetime in Table 2 and verify that Infinity does give good approximations of the central object dimensions, as given using formulae in Lamb & Miller (1995) and Miller & Lamb (1996). Finally, in Table 3, we show the diminution of calculation errors with increasing code resolution. Even though a resolution of 0.4 degrees proves to be satisfactory to produce results of good quality, if the need arises for results of greater precision, Infinity can simply be executed with a higher resolution.
Comparing Infinity with formulae in Kerr spacetime.
Infinity code convergence with resolution.
4.3. Code Elysium
Elysium code has as its main purpose to design and create a recording screen a specifically userselected distance away from the BH and the AD system. Depending on the selected program resolution, the screen has the corresponding amount of “pixels” and from each of those, a light ray is emitted perpendicularly to the screen and moves toward the disk and the BH (Fig. 18). Depending on what this ray will meet along its path, it returns information about its origin and the radiation received. We can see some examples in Fig. 19. The pictures of the second row, as well as more that will be shown further on, agree with some of the assorted pictures shown in Younsi et al. (2012).
Fig. 18.
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 solid green lines are light rays that intersect the AD and thus carry radiation. The dashed purple lines end up in the event horizon and the dotted orange ones at infinity, both without crossing any part of the AD and thus attributing nothing to the radiation total. 
Fig. 19.
Elysium results: pictures showing radiation specific intensity for four different AD models, as mentioned in Sect. 3: torus (opaque, Sect. 3.1e), wedge (opaque, Sect. 3.1d), translucent (semiopaque, Sect. 3.2c) and PD (semiopaque, Sect. 3.2b). Radiation intensity color legends shown on the top from minimum to maximum. The second row pictures are similar to those shown in Younsi et al. (2012). 
We remark here, that Elysium appears similar to the aforementioned Infinity code but is decidedly different, since it is practically its complement. Elysium therefore has equal quality and quantity of capabilities, options and result information as Infinity. The difference of the two codes is literally, as well as figuratively, a point of view. Infinity starts at the single end point of ray trajectories and integrates the equations going backwards in time. This way, it finds if the intersecting light pathway can possibly traverse a light source at any point inside the system in question. Elysium, on the other hand, has multiple starting points, the “pixels” of the screen, and integrates the equations of motion to see if an observer sitting at the pixel’s location can see any part of the disk. There are various advantages and disadvantages in both methods and depending on the type of information required each time, we can choose the most appropriate and fast code of the two.
The reason why we deemed it worthy to present and discuss this imaging code is because this code turns out to be much more than a simple imaging mechanism. The information we can gather from its various processes and results can deliver important facts about the system under study. These include information about the central BH spin (see also Sect. 4.4), the system mass and its distribution, the angular momentum allocation and flow and other such useful data.
Elysium is the best choice to make when for example we want to study or compare to results of BH observations created using GRMHD simulations, or perhaps to construct the anticipated results these will give. One such example where we could utilize Elysium is to study bright, hot material orbiting in close proximity to the central BH, such as the cases described and discussed in Broderick & Loeb (2005, 2006a). We could also use it to construct the expected finite resolution disk images of an AD with various inclinations, similarly to what is presented in Noble et al. (2007). Additionally, we could use it to investigate the radiative status of accretion flows, such as in Mościbrodzka et al. (2009), the observational appearance of radiatively inefficient accretion flows (RIAFs) and jet outflows such as in Mościbrodzka et al. (2014), or even the expected images of the SMBC in the Galactic center and its outflow regions, examined in Mościbrodzka et al. (2018). Finally, we could use it to examine possible ways to test the legitimacy of the Kerr metric or even General Relativity in the vicinity of compact objects, such as described and examined in Broderick & Loeb (2006b).
4.4. Code Tranquillity
The main purpose of Tranquillity is to find the AD inclination and then use it to make an assessment of the central BH spin. This is done by using information given by the incoming radiation emitted from the AD and traveling closer and further away from the BH and its gravitational well.
In the code’s execution, an observer at infinity looks toward a BH and the AD swirling around it. The observer records in very high resolution and measures the relative position of the AD and its first Einstein ring^{2}, the gravitationally lensed higherorder image, produced by photons circling around the BH. For this Einstein ring we henceforth use the designation “echo ring”^{3}, as this term is more intuitive, making the phenomenon easier to perceive and visualize. If the BH is nonrotating, the shapes of the AD and the echo will be concentric. If the BH is however rotating, there will be a divergence of the echo’s center, tied to the BH spin. This is more or less expected since the photons forming the echo have to travel much closer to the BH and the perturbed spacetime around it. Therefore, the higher the BH spin is, the higher this divergence will be.
The first step of Tranquillity is to determine the inclination of the AD plane compared to the line of sight. Then, the code calculates the aforementioned echo divergence compared to the AD image. Finally, by using the inclination and divergence results, we can have an estimation of the central BH spin.
The calculations for the AD inclination were rather successful, since the average declination was below 0.8 degrees for the 240 cases of different inclinations and spin parameters examined (see Table 4). Nonetheless, about 1.5% of the cases examined, gave inclination errors above average. The estimation errors were about five degrees or less, and were caused by specific conditions of the setup that had a BH first echo ring appear in very particular and peculiar locations. The model selected for the AD of the object in question, does not appear to play any significant part in the inclination assessment so far.
Tranquillity code inclination estimation errors.
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 on the model examined, different disk models will follow slightly different divergence plots, directly affected by the disk’s inner edge radius (see Abramowicz et al. 2010). We clearly state here however, that our code does not use or rely at all on the inner edge of the disk, but only on the appearance of the entire disk as a whole. The disk model thus, does not prove to be an insurmountable problem, since it only causes a small recalibration to the divergence plot, as we see in the second part of this work.
4.5. Code Burning Arrow
The Burning Arrow code has as a main purpose to study the BH massive particles orbit degradation due to the hot disk radiation.
In order to study the particle motion, the code must solve the general relativistic equations of motion, equivalent to the Classical Newton’s Laws of Motion. Starting from the first law, in absence of general relativistic forces, the particle in question will follow a geodesic through spacetime. This geodesic obeys the equation:
where τ here is the proper time, u^{α} = dx^{α}/dτ is the fourvelocity and the Christoffel symbols, with commas used for partial derivatives. As in Newton’s equation, the zero term on the right stands for the absence of acceleration. Solving the above, gives the various geodesic solutions that describe among others, special circular orbits such as the ISCO, the photon sphere etc. The previous equation can also be rewritten as:
by using the chain rule in order to have derivatives by the coordinate time x^{0} = t. This, nevertheless, has proved to be more prone to error accumulation in our work and so we choose to work with the proper time equations instead.
The above equation gives later on, rise to the equivalent of Newton’s second law of motion. It also gives us a way to calculate the accelerations present in a problem, such as those generated by the radiation forces in our case. As in the Infinity code, we consider different cases of velocity behaviors, such as the SANE and MAD models as mentioned before, in addition to the typical Bardeen et al. (1972) angular velocity profile.
In the beginning of the execution and depending on the study requested, the code looks up into the Infinity code result files and reads the appropriate ones, relevant to each case. It then uses these results to generate proper functions that give the radiation stress – energy tensor at every point of spacetime used in the problem at hand. This gives the fouracceleration due to radiation at each point from the extension of the aforementioned equation as:
Since the radiation acceleration is known in spacetime from the previous results, we can then solve these eight differential equations and get the fourvelocity and fourpositions describing the sought particle trajectory.
The Burning Arrow output for degenerating orbits can be seen in Fig. 20. Without the presence of radiation, the material electrons would continue moving in circular orbits, shown in the figure with the dashed circles. The presence of radiation, nevertheless, causes the electron orbits to destabilize. The electrons then, either fall into the central object or leave the AD system (gray annulus) and escape to infinity. We show here some of the results of the radiation effects on three different velocity profiles. In the following results part of this work, we examine and discuss thoroughly the noteworthy and interesting repercussions of these interactions.
Fig. 20.
Burning Arrow results: degradation of equatorial orbits around a nonrotating BH due to the semiopaque (LFM, Sect. 3.2d) disk’s thermal radiation. With solid black is the Bardeen et al. (1972) velocity profile, with dashed cyan and with dotted red an approximation of a SANE and a MAD model velocity profile respectively (Narayan et al. 2003, 2012; Penna et al. 2013b). 
5. Discussion
In this work we presented the physics and the setup required in order to examine the radiation field produced by the hot accretion disk orbiting an astrophysical black hole. We subsequently described the five families of codes written to study these objects and investigate some of their facets. These codes study photon as well as massive particle trajectories and examine the dynamics of the black hole and accretion disk systems from close up and from infinity. Additionally, the codes create inclination  spin templates for assorted disk models in an attempt to assess the black hole’s spin from observational data. In more, the codes study particle trajectories under the influence of relativistic effects and the presence of the disk’s thermal radiation.
We looked into optically thin and thick, and geometrically thin and thick accretion disks. Also, we chose to apply certain radial and azimuthal velocity and temperature profiles that follow and agree with the majority of studies in environments such as the ones we study here. Various alterations and reassessments can be made nonetheless, in order to better suit other situations, other disk models, prograde or retrograde material motion, spacetime environments etc.
In total, for this task we ran more than 8 thousand simulations and we studied more than 4 billion photon trajectories. In order to examine and confirm the validity of our results, we corroborated that our codes successfully calculate analytically described or assorted expected trajectories for massless and massive particle, such as those explained in Bardeen et al. (1972) and Levin & PerezGiz (2008). Additionally, we showed that we correctly obtain trajectories of free photons moving in strong gravity environments, such as those appearing in Bini et al. (2015). In more, we adapted our codes to look toward central sources of photons, namely stars, instead of accretion tori in order to compare and crosscheck our results with other, relevant studies including Abramowicz et al. (1990), Miller & Lamb (1993, 1996), Lamb & Miller (1995), with which we found agreement for observational and physical quantities. Finally, we saw in brief here and will see in more detail further on, that our codes create observationlike images of various accretion disks at infinity which are in agreement with other similar works, such as Fuerst & Wu (2004, 2007), Younsi et al. (2012).
The radiation field in such environments proves to be very complex for many reasons, some of which are the general relativistic effects present, and also the fact that the radiation source is no longer central, but instead extended in a large azimuthal and poloidal volume of space. In the second part of his work (Koutsantoniou, in prep.), we look into the results of our codes and discuss photon trajectories and bundles, disk images and induced radiation forces, photographs of systems from infinity, their usage to estimate a black hole spin and particle orbit degeneration due to radiation. We then assess how we can use this information to draw conclusions about the physics and phenomena of these systems, their dynamical condition, their equilibrium and perhaps their evolution.
Direct link: youtube.com/channel/UCJ4v8rSg390gt9kQfVtfx5A
An Einstein, Chwolson or echo ring is created when light from the hot disk passes close to the BH and reaches an observer. Due to the BHAD setup we study here, the relative locations of these objects causes the gravitational lensing taking place in the system to display a full ring of light, created by the accretion disk’s upper and lower parts.
We refer to this ring as an echo due to its similarities with the light echo phenomenon, that is analogous to the sound echo, but with light instead of sound. The light echo feature, however, is more often connected to nova and supernova events (e.g., V838 Monocerotis, 2002 outburst), than with BHs (e.g., V404 Cygni, 2015 outburst).
Acknowledgments
This work was in part supported by the General Secretariat for Research and Technology of Greece and the European Social Fund in the framework of Action “Excellence”. Part of this work was performed at the Research Center for Astronomy and Applied Mathematics of the Academy of Athens. The algorithms discussed in this work can be found at the following addresses: Omega: https://gitlab.com/leelamichaels/Omega.git, Infinity: https://gitlab.com/leelamichaels/Infinity.git, Elysium: https://gitlab.com/leelamichaels/Elysium.git, Tranquillity: https://gitlab.com/leelamichaels/Tranquillity.git, Burning Arrow: https://gitlab.com/leelamichaels/Burning_Arrow.git
References
 Abramowicz, M. A., & Fragile, P. C. 2013, Liv. Rev. Rel., 16, 1 [Google Scholar]
 Abramowicz, M., Jaroszynski, M., & Sikora, M. 1978, A&A, 63, 221 [NASA ADS] [Google Scholar]
 Abramowicz, M. A., Czerny, B., Lasota, J. P., & Szuszkiewicz, E. 1988, ApJ, 332, 646 [Google Scholar]
 Abramowicz, M. A., Ellis, G. F. R., & Lanza, A. 1990, ApJ, 361, 470 [NASA ADS] [CrossRef] [Google Scholar]
 Abramowicz, M. A., Chen, X. M., Granath, M., & Lasota, J. P. 1996, ApJ, 471, 762 [NASA ADS] [CrossRef] [Google Scholar]
 Abramowicz, M. A., Jaroszyński, M., Kato, S., et al. 2010, A&A, 521, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Artemova, I. V., BisnovatyiKogan, G. S., Bjoernsson, G., & Novikov, I. D. 1996, ApJ, 456, 119 [NASA ADS] [CrossRef] [Google Scholar]
 Asada, K., Nakamura, M., & Pu, H.Y. 2016, ApJ, 833, 56 [NASA ADS] [CrossRef] [Google Scholar]
 Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [Google Scholar]
 Balbus, S. A., & Hawley, J. F. 1992, ApJ, 400, 610 [Google Scholar]
 Balbus, S. A., & Hawley, J. F. 1998, Rev. Mod. Phys., 70, 1 [Google Scholar]
 Bardeen, J. M. 1970, ApJ, 162, 71 [Google Scholar]
 Bardeen, J. M., Press, W. H., & Teukolsky, S. A. 1972, ApJ, 178, 347 [Google Scholar]
 Beloborodov, A. M. 1998, MNRAS, 297, 739 [NASA ADS] [CrossRef] [Google Scholar]
 Beloborodov, A. M. 1999, ASP Conf. Ser., 161, 295 [Google Scholar]
 Beloborodov, A. M. 2001, Adv. Space Res., 28, 411 [NASA ADS] [CrossRef] [Google Scholar]
 Bildsten, L., & Rutledge, R. E. 2001, in The Neutron Star  Black Hole Connection, eds. C. Kouveliotou, J. Ventura, & E. van den Heuvel, 567, 245 [NASA ADS] [CrossRef] [Google Scholar]
 Bini, D., Jantzen, R. T., & Stella, L. 2009, CQG, 26, 055009 [NASA ADS] [CrossRef] [Google Scholar]
 Bini, D., Geralico, A., Jantzen, R. T., Semerák, O., & Stella, L. 2011, CQG, 28, 035008 [NASA ADS] [CrossRef] [Google Scholar]
 Bini, D., Geralico, A., Jantzen, R. T., & Semerák, O. 2015, MNRAS, 446, 2317 [CrossRef] [Google Scholar]
 BisnovatyiKogan, G. S. 2001, Discrete Dyn. Nat. Soc., 6 [Google Scholar]
 BisnovatyiKogan, G. S., Lovelace, R. V. E., & Belinski, V. A. 2002, ApJ, 580, 380 [NASA ADS] [CrossRef] [Google Scholar]
 Blaes, O. M. 2004, in Accretion Discs, Jets and High Energy Phenomena in Astrophysics, eds. V. Beskin, G. Henri, F. Menard, et al., 78, 137 [NASA ADS] [Google Scholar]
 Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433 [Google Scholar]
 Blandford, R., Agol, E., Broderick, A., et al. 2002, in Astrophysical Spectropolarimetry, eds. J. TrujilloBueno, F. MorenoInsertis, & F. Sánchez, 177 [Google Scholar]
 Bondi, H. 1952, MNRAS, 112, 195 [Google Scholar]
 Broderick, A. E., & Loeb, A. 2005, MNRAS, 363, 353 [NASA ADS] [CrossRef] [Google Scholar]
 Broderick, A. E., & Loeb, A. 2006a, MNRAS, 367, 905 [Google Scholar]
 Broderick, A. E., & Loeb, A. 2006b, J. Phys. Conf. Ser., 54, 448 [NASA ADS] [CrossRef] [Google Scholar]
 Carter, B. 1968, Phys. Rev., 174, 1559 [Google Scholar]
 Chandrasekhar, S. 1983, The Mathematical Theory of Black Holes [Google Scholar]
 Chatterjee, K., Liska, M., Tchekhovskoy, A., & Markoff, S. B. 2019, MNRAS, 490, 2200 [Google Scholar]
 Chen, W.C., & Podsiadlowski, P. 2016, ApJ, 830, 131 [CrossRef] [Google Scholar]
 ChoquetBruhat, Y., DeWittMorette, C., & DillardBleik, M. 1977, Analysis Manifolds and Physics [Google Scholar]
 Contopoulos, I., & Kazanas, D. 1998, ApJ, 508, 859 [Google Scholar]
 Contopoulos, I., & Papadopoulos, D. B. 2012, MNRAS, 425, 147 [NASA ADS] [CrossRef] [Google Scholar]
 Contopoulos, I., Kazanas, D., & Christodoulou, D. M. 2006, ApJ, 652, 1451 [NASA ADS] [CrossRef] [Google Scholar]
 Cunningham, C. T. 1975, ApJ, 202, 788 [NASA ADS] [CrossRef] [Google Scholar]
 Cunningham, C. 1976, ApJ, 208, 534 [NASA ADS] [CrossRef] [Google Scholar]
 Davelaar, J., Bronzwaer, T., Kok, D., et al. 2018, Comput. Astrophys. Cosmol., 5, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Done, C., Gierliński, M., & Kubota, A. 2007, A&Arv, 15, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Dubus, G. 2003, EAS Pub. Ser., 7, 283 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Esin, A. A., McClintock, J. E., & Narayan, R. 1997, ApJ, 489, 865 [NASA ADS] [CrossRef] [Google Scholar]
 Event Horizon Telescope Collaboration IV, 2019, ApJ, 875, L4 [CrossRef] [Google Scholar]
 Fender, R. 2002, in Relativistic Outflows from Xray Binaries (‘Microquasars’), eds. A. W. Guthmann, M. Georganopoulos, A. Marcowith, & K. Manolakou, 589, 101 [Google Scholar]
 Fuerst, S. V. 2006, PhD Thesis, Mullard Space Science Laboratory, University College London, Holmbury St. Mary, Dorking, Surrey RH5 6NT, UK [Google Scholar]
 Fuerst, S. V., & Wu, K. 2004, A&A, 424, 733 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fuerst, S. V., & Wu, K. 2007, A&A, 474, 55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Haggard, D., Cool, A. M., Anderson, J., et al. 2004, ApJ, 613, 512 [NASA ADS] [CrossRef] [Google Scholar]
 Hawley, J. F., & Balbus, S. A. 1991, ApJ, 376, 223 [NASA ADS] [CrossRef] [Google Scholar]
 Hawley, J. F., & Balbus, S. A. 1992, ApJ, 400, 595 [NASA ADS] [CrossRef] [Google Scholar]
 Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742 [Google Scholar]
 Heinke, C. O., Ivanova, N., Engel, M. C., et al. 2013, ApJ, 768, 184 [NASA ADS] [CrossRef] [Google Scholar]
 Igumenshchev, I. V., Narayan, R., & Abramowicz, M. A. 2003, ApJ, 592, 1042 [Google Scholar]
 Inoue, H., & Hoshi, R. 1987, ApJ, 322, 320 [NASA ADS] [CrossRef] [Google Scholar]
 Komissarov, S. S. 1999, MNRAS, 308, 1069 [Google Scholar]
 Komissarov, S. S. 2001, MNRAS, 326, L41 [Google Scholar]
 Komissarov, S., & Porth, O. 2021, New Astron. Rev., 92, 101610 [Google Scholar]
 Komissarov, S. S., Barkov, M. V., Vlahakis, N., & Königl, A. 2007, MNRAS, 380, 51 [Google Scholar]
 Konoplya, R. A., Kunz, J., & Zhidenko, A. 2021, JCAP, 12, 002 [Google Scholar]
 Koutsantoniou, L. E. 2014, Master’s Thesis, Department of Astrophysics, Astronomy and Mechanics, Faculty of Physics, University of Athens, Panepistimiopolis Zografos, Athens 15784, Greece [Google Scholar]
 Koutsantoniou, L. E., & Contopoulos, I. 2014, ApJ, 794, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Kozlowski, M., Jaroszynski, M., & Abramowicz, M. A. 1978, A&A, 63, 209 [Google Scholar]
 Krolik, J. H. 1999a, Active Galactic Nuclei: from the Central Black hole to the Galactic Environment [Google Scholar]
 Krolik, J. H. 1999b, ASP Conf. Ser., 161, 315 [NASA ADS] [Google Scholar]
 Kylafis, N. D., Contopoulos, I., Kazanas, D., & Christodoulou, D. M. 2012, A&A, 538, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lamb, F. K., & Miller, M. C. 1995, ApJ, 439, 828 [NASA ADS] [CrossRef] [Google Scholar]
 Lasota, J. P. 1999, Phys. Rep., 311, 247 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, H. K., Wijers, R. A. M. J., & Brown, G. E. 2000, Phys. Rep., 325, 83 [NASA ADS] [CrossRef] [Google Scholar]
 Levin, J., & Grossman, R. 2009, Phys. Rev. D, 79, 043016 [CrossRef] [Google Scholar]
 Levin, J., & PerezGiz, G. 2008, Phys. Rev. D, 77, 103005 [CrossRef] [Google Scholar]
 Livio, M., Ogilvie, G. I., & Pringle, J. E. 1999, ApJ, 512, 100 [NASA ADS] [CrossRef] [Google Scholar]
 Longair, M. S. 2011, High Energy Astrophysics [Google Scholar]
 Mahlmann, J. F., Levinson, A., & Aloy, M. A. 2020, MNRAS, 494, 4203 [Google Scholar]
 McKinney, J. C. 2005, ApJ, 630, L5 [NASA ADS] [CrossRef] [Google Scholar]
 McKinney, J. C., Tchekhovskoy, A., & Blandford, R. D. 2012, MNRAS, 423, 3083 [Google Scholar]
 MeyerHofmeister, E., Liu, B. F., & Meyer, F. 2009, A&A, 508, 329 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Miller, M. C., & Lamb, F. K. 1993, ApJ, 413, L43 [NASA ADS] [CrossRef] [Google Scholar]
 Miller, M. C., & Lamb, F. K. 1996, ApJ, 470, 1033 [CrossRef] [Google Scholar]
 Misner, C. W., Thorne, K. S., & Wheeler, J. A. 1973, Gravitation (San Francisco: W.H. Freeman and Co) [Google Scholar]
 Müller, A. 2004, Ph.D. Thesis: Black hole astrophysics: Magnetohydrodynamics on the Kerr geometry, Landessternwarte Heidelberg, Germany [Google Scholar]
 Mościbrodzka, M., Gammie, C. F., Dolence, J. C., Shiokawa, H., & Leung, P. K. 2009, ApJ, 706, 497 [Google Scholar]
 Mościbrodzka, M., Falcke, H., Shiokawa, H., & Gammie, C. F. 2014, A&A, 570, A7 [Google Scholar]
 Mościbrodzka, M., Falcke, H., & Noble, S. 2018, in Fourteenth Marcel Grossmann Meeting  MG14, eds. M. Bianchi, R. T. Jansen, & R. Ruffini, 3519 [Google Scholar]
 Mueller, T., & Grave, F. 2009, ArXiv eprints [arXiv:0904.4184] [Google Scholar]
 Nakamura, M., Asada, K., Hada, K., et al. 2018, ApJ, 868, 146 [Google Scholar]
 Narayan, R., & Yi, I. 1994, ApJ, 428, L13 [Google Scholar]
 Narayan, R., & Yi, I. 1995, ApJ, 444, 231 [Google Scholar]
 Narayan, R., & McClintock, J. E. 2008, New Astron. Rev., 51, 733 [CrossRef] [Google Scholar]
 Narayan, R., Igumenshchev, I. V., & Abramowicz, M. A. 2003, PASJ, 55, L69 [NASA ADS] [Google Scholar]
 Narayan, R., Sądowski, A., Penna, R. F., & Kulkarni, A. K. 2012, MNRAS, 426, 3241 [Google Scholar]
 Noble, S. C., Leung, P. K., Gammie, C. F., & Book, L. G. 2007, CQG, 24, S259 [NASA ADS] [CrossRef] [Google Scholar]
 Noble, S. C., Krolik, J. H., Schnittman, J. D., & Hawley, J. F. 2011, ApJ, 743, 115 [NASA ADS] [CrossRef] [Google Scholar]
 Novikov, I. D., & Thorne, K. S. 1973, Black Holes (Les Astres Occlus), 343 [Google Scholar]
 Parfrey, K., Philippov, A., & Cerutti, B. 2019, Phys. Rev. Let., 122, 035101 [NASA ADS] [CrossRef] [Google Scholar]
 Park, J., Hada, K., Kino, M., et al. 2019, ApJ, 871, 257 [Google Scholar]
 Penna, R. F., Narayan, R., & Sądowski, A. 2013a, MNRAS, 436, 3741 [Google Scholar]
 Penna, R. F., Sądowski, A., Kulkarni, A. K., & Narayan, R. 2013b, MNRAS, 428, 2255 [NASA ADS] [CrossRef] [Google Scholar]
 Penrose, R., & Floyd, R. M. 1971, Nat. Phys. Sci., 229, 177 [NASA ADS] [CrossRef] [Google Scholar]
 Pessah, M. E., Chan, C.K., & Psaltis, D. 2007, ApJ, 668, L51 [Google Scholar]
 Podsiadlowski, P., Rappaport, S., & Pfahl, E. D. 2002, ApJ, 565, 1107 [NASA ADS] [CrossRef] [Google Scholar]
 Poynting, J. H. 1903, MNRAS, 64, A1 [Google Scholar]
 Quataert, E. 2001, ASP Conf. Ser., 224, 71 [NASA ADS] [Google Scholar]
 Robertson, H. P. 1937, MNRAS, 97, 423 [NASA ADS] [Google Scholar]
 Rybicki, G. B., & Lightman, A. P. 1986, Radiative Processes in Astrophysics (WileyVCH) [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 500, 33 [NASA ADS] [Google Scholar]
 Sądowski, A. 2009, ApJS, 183, 171 [CrossRef] [Google Scholar]
 Sadowski, A. 2011, Ph.D. Thesis: Slim accretion disks around black holes, Nicolaus Copernicus Astronomical Center, Warsaw, Poland, [arXiv:1108.0396] [Google Scholar]
 Sądowski, A. 2016, MNRAS, 459, 4397 [CrossRef] [Google Scholar]
 Sądowski, A., Lasota, J.P., Abramowicz, M. A., & Narayan, R. 2016, MNRAS, 456, 3915 [CrossRef] [Google Scholar]
 Takahashi, A., Fukue, J., Sanbuichi, K., & Umemura, M. 1995, PASJ, 47, 425 [NASA ADS] [Google Scholar]
 Tauris, T. M., & van den Heuvel, E. P. J. 2006, Formation and evolution of compact stellar Xray sources, 39, 623 [NASA ADS] [CrossRef] [Google Scholar]
 Tauris, T. M., van den Heuvel, E. P. J., & Savonije, G. J. 2000, ApJ, 530, L93 [NASA ADS] [CrossRef] [Google Scholar]
 Tchekhovskoy, A., Narayan, R., & McKinney, J. C. 2011, MNRAS, 418, L79 [NASA ADS] [CrossRef] [Google Scholar]
 Thorne, K. S., & Price, R. H. 1975, ApJ, 195, L101 [NASA ADS] [CrossRef] [Google Scholar]
 Turner, N. J., Stone, J. M., Krolik, J. H., & Sano, T. 2003, ApJ, 593, 992 [NASA ADS] [CrossRef] [Google Scholar]
 van Haaften, L. M., Nelemans, G., Voss, R., Wood, M. A., & Kuijpers, J. 2012, A&A, 537, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Verbunt, F. 1999, ASP Conf. Ser., 160, 21 [NASA ADS] [Google Scholar]
 Vlahakis, N., & Königl, A. 2004, ApJ, 605, 656 [Google Scholar]
 Wheeler, C. J. 2004, Adv. Space Res., 34, 2744 [NASA ADS] [CrossRef] [Google Scholar]
 Wilkins, D. C. 1972, Phys. Rev. D, 5, 814 [CrossRef] [Google Scholar]
 Younsi, Z., Wu, K., & Fuerst, S. V. 2012, A&A, 545, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Yuan, Y., Blandford, R. D., & Wilkins, D. R. 2019, MNRAS, 484, 4920 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1.
Local sky around the target particle. Left: forward in time (the photon moves toward the target). Right: moving backwards in time (the photon moves away from the target). For the incoming photon, we define angles ã and b̃ similar to the typical polar angle θ and azimuthal angle ϕ of spherical coordinate systems. The disk in stripes at the top left represents the BH event horizon. 

In the text 
Fig. 2.
Left: emitting matter in the lab frame K, moving with velocity along the vertical axis in a tube of width d. A photon of frequency ν crosses the tube, traveling at an angle θ from the vertical axis. Right: emitting matter in its rest frame K′ in the tube. The photon with frequency ν ′ now appears to cross the tube at angle θ ′. 

In the text 
Fig. 3.
Emission of a photon from the hot AD. The photon is emitted along n̂ from an element of matter moving with velocity . The photon is emitted at an angle ψ relative to the accreting material motion. 

In the text 
Fig. 4.
Schematic of the frequency changes. A Doppler shift is used to move from the frame comoving with the emitting surface (left) to the LNRF at the radius of the surface (middle), taking into account the emission source rotation. Then we move to the receiving particle frame (right), 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. 

In the text 
Fig. 5.
ADs “Band” (left) and “Disk” (right) constituted by twodimensional surfaces (Sect. 3.1a and b respectively). In the center of each image is in black the BH event horizon and around it in gray, its ergosphere. We see the AD in the yellow (hotter) and red tones (colder). 

In the text 
Fig. 6.
Opaque ADs cross sections for the spins studied (Sect. 3.1c–f). 

In the text 
Fig. 7.
ORST cross sections (Sect. 3.1f). In (a) we can see the effects caused by the change of the radius of Keplerian rotation speed r_{K}. In (b) we see the tori shapes and sizes for various values of the parameter n_{p}. 

In the text 
Fig. 8.
Polish doughnut number density cross sections (Sect. 3.2b,c). The center of all tori lies at (12,0) and has a number density of 10^{18} cm^{−3}. The top row shows the disk for BH spin parameters a = 0 and a = 0.5, while the bottom row for a = 0.9 and a = 0.998. 

In the text 
Fig. 9.
LFM model number density cross section for any spin parameter (Sect. 3.2d,e). The center of the tori lies at (2r_{ISCO},0) and has a number density of 10^{18} cm^{−3}. 

In the text 
Fig. 10.
ADs created from combinations of previously referred models (from top to bottom): wedge & slab, polish doughnut & torus, torus & wedge and polish doughnut & spheres. 

In the text 
Fig. 11.
Control interface of the Omega code, as described in Sect. 4.1. In light green is the requested photon trajectory and with the cyan point we denote its point of origin. 

In the text 
Fig. 12.
Tests for the validity of Omega. We confirm that the code successfully produces fundamental particle orbits, such as the ones mentioned in Sect. 2.1. In each picture we have the BH situated at the axes origin and with gray is its event horizon. In (a) we have in red the photon ring produced numerically by the code and on top of it, in the black dashed circle, its theoretically expected position. In (b) we see in blue the ISCO produced numerically and on top of it, in the black dotted circle, its theoretically expected position. In (c) we set the BH mass M = 0 (with light gray showing where its horizon would have been) and examine the trajectory of a photon. With any BH mass M > 0, the photon would have to infall and cross the event horizon, but with no BH mass present, the spacetime is Minkowski and the photon travels in a straight line. 

In the text 
Fig. 13.
Free photon trajectories in the poloidal x − z plane. Photons are emitted perpendicularly upward from various points of the equatorial plane of the system. The closest to the BH emission points are (3, 0) and (−3, 0) and taking into account the angle of emission, these photons should and indeed do follow the photon ring trajectory (compare to Bini et al. 2015, Fig. 1). 

In the text 
Fig. 14.
Omega code successfully creating massive particle orbits. These trajectories are in agreement with the corresponding orbits shown in Levin & PerezGiz (2008). 

In the text 
Fig. 15.
Two instances of the first version of the resolution enhancement process (625 times here). In the top row of each part is the original image created by the AD observation. On the bottom is the resulting image after the resolution enhancement. These functions are now rewritten and incorporated into Infinity to enhance results by a number of times. 

In the text 
Fig. 16.
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 AD 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 BH’s ergosphere and has thus a radius r = 2M. The actual region is considered to exist above a certain height, further away from the BH. In this region, we can have matter at rest in the local rest frame and matter flowing outward with a certain velocity. We may also consider a subregion there (see e.g., Asada et al. 2016; Park et al. 2019), a narrower cone or cylinder with radius r = M with stronger outflow and faster velocity or even with velocity of the opposite direction. 

In the text 
Fig. 17.
Infinity code: Mollweide sky maps of the frequency integrated specific intensity for tori with downward angular velocity vector. The color scale for each image is displayed below it, from minimum to maximum. Top: wedge AD model (opaque, Sect. 3.1d), bottom: LFM AD model (semiopaque, Sect. 3.2d). 

In the text 
Fig. 18.
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 solid green lines are light rays that intersect the AD and thus carry radiation. The dashed purple lines end up in the event horizon and the dotted orange ones at infinity, both without crossing any part of the AD and thus attributing nothing to the radiation total. 

In the text 
Fig. 19.
Elysium results: pictures showing radiation specific intensity for four different AD models, as mentioned in Sect. 3: torus (opaque, Sect. 3.1e), wedge (opaque, Sect. 3.1d), translucent (semiopaque, Sect. 3.2c) and PD (semiopaque, Sect. 3.2b). Radiation intensity color legends shown on the top from minimum to maximum. The second row pictures are similar to those shown in Younsi et al. (2012). 

In the text 
Fig. 20.
Burning Arrow results: degradation of equatorial orbits around a nonrotating BH due to the semiopaque (LFM, Sect. 3.2d) disk’s thermal radiation. With solid black is the Bardeen et al. (1972) velocity profile, with dashed cyan and with dotted red an approximation of a SANE and a MAD model velocity profile respectively (Narayan et al. 2003, 2012; Penna et al. 2013b). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.