RAIKOU: A General Relativistic, Multi-wavelength Radiative Transfer Code

We present a general relativistic, ray-tracing radiative transfer code RAIKOU for multi-wavlength studies of spectra and images including the black hole shadows around Kerr black holes. Important radiative processes in hot plasmas around black holes, i.e., (cyclo-)synchrotron, bremsstrahlung emission/absorption and Compton/inverse-Compton scattering, are incorporated. The Maxwell-J\"uttner and single/broken power-law electron distribution functions are implemented to calculate the radiative transfer via both of the thermal and the nonthermal electrons. Two calculation algorithms are implemented for studies of both the images and broadband spectra. An observer-to-emitter algorithm, which inversely solve the radiative transfer equation from the observer screen to emitting plasmas, is suitable for efficient calculations of the images, e.g., the black hole shadows, and spectra without the Compton effects. On the other hand, an emitter-to-observer algorithm, by which photons are transported with a Monte-Carlo method including the effects of Compton/inverse-Compton scatterings, enables us to compute multi-wavelength spectra with their energy bands broadly ranging from radio to very-high-energy gamma-ray. The code is generally applicable to accretion flows around Kerr black holes with relativistic jets and winds/coronae with various mass accretion rate (i.e., radiatively inefficient accretion flows, super-Eddington accretion flows, and others). We demonstrate an application of the code to a radiatively innefficent accretion flow onto a supermassive black hole.


INTRODUCTION
Astrophysical black holes (BHs) are thought to be powered by the gravitational energy release via the accretion flows. The energy release in the accretion flows will result in the emission of an enormous amount of electromagnetic radiation (e.g., Shakura & Sunyaev 1973;Abramowicz et al. 1988;Narayan & Yi 1994) and the ejection of powerful outflows (i.e., relativistic jets and mildly relativistic winds), which are theoretically suggested to be accelerated by magnetic processes (Blandford & Znajek 1977;Blandford & Payne 1982;Lynden-Bell 1996), and radiative force (Sikora 1981). In order to study BH spacetimes and physics of these plasma dynamics, computations of general relativistic radiative transfer (GRRT) will be a powerful tool, because GRRT calculations of images and spectra will enable us to directly compare the theoretical models with the observed data.
One of the most important observable features of BHs is their shadows or silhouettes (Bardeen 1973;Luminet 1979;Falcke et al. 2000;Takahashi 2004;Broderick & Loeb 2006;Noble et al. 2007;Dexter et al. 2012;Mościbrodzka et al. 2016;Chael et al. 2018), because their images are formed as a consequence of the significant light-bending effect due to the extremely strong gravity of the BHs. Recently, the Event Horizon Tele-scope (EHT) detected the image of BH shadow in the elliptical galaxy M87* (Event Horizon Telescope Collaboration 2019a,b,c,d,e,f) and it provided us a strong and important evidence of the supermassive BHs (SMBHs), whose formation scenario is still controversial (Rees 1984;Volonteri 2012).
The diameter of the observed photon ring, which is a bright ring image formed by the photons passing through near the unstable photon orbit around the BHs, is ∼ 40µas (Event Horizon Telescope Collaboration 2019d). The BH mass can be estimated to be 6.5 × 10 9 M (Event Horizon Telescope Collaboration 2019f) because the diameter of the photon ring is roughly proportional to the BH mass, i.e., ∼10r g , where r g is the gravitational radius defined as r g ≡ GM/c 2 , G is the gravitaional constant, M is the BH mass, and c is the speed of light. However, the constrain of the BH spin is still grand challenge because the diameter of the photon ring depends weakly on the BH spin, in practice the diameter changes within ±5% (e.g., Psaltis et al. 2015) and it would be difficult to be identified using the current EHT. In addition, the dynamics of accretion flows and launching mechanism of jets are still enigmatic. It would be, therefore, important to explore these problem using the additional feature, e.g., polarization (Broderick & Loeb 2006;Shcherbakov et al. 2012;Dexter 2016;Mościbrodzka et al. 2017;Gold et al. 2017;Tsunetoe et al. 2020Tsunetoe et al. , 2021Event Horizon Telescope Collaboration et al. 2021) and/or multi-wavelength spectra (Ohsuga et al. 2005;Mościbrodzka et al. 2009Mościbrodzka et al. , 2016Chael et al. 2016Chael et al. , 2018, where the latter approach is the motivation of the code development in this paper. The multi-wevelength radiative transfer calculations, whose energy band ranging from radio to very-highenergy (VHE) gamma-ray, potentially constrain the various parameters with respect to BH spacetimes and theoretical models of the accretion flows and jets via the simultaneous comparison between theoretical and observed radiative flux in the various energy bands. In order to calculate the spectral energy distributions (SEDs) of accretion flows and jets around BHs above the X-ray band, it is necessary to incorporate the effects of the Compton and inverse-Compton scattering (e.g., Rybicki & Lightman 1979). For the calculations of the Comptonized SEDs, Monte-Carlo radiative transfer (e.g., Pozdnyakov et al. 1977Pozdnyakov et al. , 1983Gorecki & Wilczewski 1984;Canfield et al. 1987;Dolence et al. 2009) are suitable for accurate computations of the Compton/inverse-Compton scattering without introducing approximations (e.g., the Kompaneets approximation, which is applicable when the radiative field are isotropic in the fluid-rest frame (Kompaneets 1957)), although the Monte-Calro calculations require high computational costs.
Despite the most of the previous works on GRRT calculations focused on the thermal electrons, the observations of BHs strongly demonstrates the evidences of the acceleration of electrons. For example, the VHE gamma-ray observation of M87 by MAGIC Telescope indicates the generation of nonthermal electrons at the distance ∼100r g from the BH in M87 (MAGIC Collaboration et al. 2020). Recent simultaneous multiwavelength observation campaign for EHT 2017 found that the gamma-ray emission is not generated at the same place with the ring-emission in the vicinity of the BH if the isotropic distribution of the nonthermal electrons is assumed (EHT MWL Science Working Group et al. 2021). On the other hand, IC 310 shows the subday scale, rapidly varying VHE gamma-ray emission, which is interpreted as the emission orginated from the nonthermal electrons accelerated from the gap electric field in the vinicity of the BH (Aleksić et al. 2014). In the lower frequency bands, the synchrotron emission via nonthermal electrons can reproduce the radio emission at ν < 100 GHz (Yuan et al. 2003;Davelaar et al. 2018) and it will be important to consider the images and spectra at 230 GHz, infrared , and X-ray bands in Sgr A* (e.g., Ball et al. 2016;Chael et al. 2017;Mao et al. 2017;Chatterjee et al. 2020;Scepi et al. 2021).
Calculations of the Comptonized SEDs are important for the general studies on BHs, whose mass ranging from stellar-mass BHs to SMBHs, with the mass accretion rate ranging from sub-Eddington to super-Eddington accretion-rate. For example, in the super-Eddington accretion flows, the Compton/inverse-Compton scattering in jets and winds drastically modifies the SEDs (Kawashima et al. 2012;Kitaki et al. 2017;Narayan et al. 2017;Curd & Narayan 2019;Kawanaka & Mineshige 2021), which will be important to explain the SEDs of ultraluminous X-ray sources (ULXs) and tidal disruption events (TDEs). The origin of the hard X-ray photons from the BH accretion flows with marginally Eddington accretion rate, which is an open question and is proposed to be Comptonized in hot coronae formed above the standard disks or inside the truncated accretion disks, can be also approached (Kawanaka et al. 2008;Davis et al. 2009;. In this paper, we describe our newly developed general relativistic, multi-wavelength radiative transfer code RAIKOU (来 光), by which the multi-wavelength SEDs from radio to VHE gamma-ray bands and images of the BH shadow, accretion flows and outflows can be calculated. The electron distribution functions of the thermal electrons with the Maxwell Jüttner distribu-tion and the nonthermal electrons with single/broken power-law distributions are incorporated. The VHE gamma-ray emission above TeV energy bands can be computed because of the implement of the inverse-Compton scattering processes via nonthermal electrons. The code is applicable to calculate the images and spectra of accreting BHs by post-processing general relativisitic (radiation-)magnetohydrodynamic simulations, hereafter GR(R)MHD, as well as by using semi-analytic models. In §2-5, the formalism and numerical method are described. In §6, numerical tests are shown. In §7, we demonstrate an application of the code to calculations of multi-wavelength images and spectra of a radiatively inefficient accretion flow onto Sgr A* based on GRMHD simulations. We devote §8 as the summary and discussion.

OVERVIEW OF METHOD
We describe the formalism of equations and algorithms implemented to RAIKOU. For the various purposes of the radiative transfer calculations, two types of calculation algorithm are employed. One is the observerto-emitter algorithm, by which the radiative transfer equations are integrated from the observer-screen to the emitting plasma. This method enables us to efficiently compute the accurate images and spectra in the limit of plasmas being optically thin against the inverse-Compton scattering. The other is the emitterto-observer algorithm, by which the superphotons are transported from the emitting plasma to the observer with a Monte-Carlo method. By using the emitter-toobserver algorithm, one can compute broadband spectra from radio to VHE gamma-rays which is significantly affected by Compton/inverse-Compoton scattering.
The same formalism and algorithm of the ray-tracing part, i.e., the integration module of the geodesic equations of photons, are used for both of the oberver-toemitter and the emitter-to-observer algorithm. The details of the radiative transfer part including various radiative processes, i.e., the emission, abosorption, and scattering, are presented in the subsequent subsections.

GEODESICS
We integrate the Hamilton's canonical equations of motion, which describe time evolution of the position and momentum of photons r, θ, ϕ, p r and p θ in Boyer-Lindquist (BL) coordinates, to solve the geodesic equations for photons . We do not solve the time evolution of conserved variable of photons, i.e., time and azimuthal components of photons p t and p ϕ during the calculation of the geodesics, while these are updated when the photons are scattered via the Compton process.
The metric of the BL coordinate is written as follows: where, in the geometrical unit G = c = 1, The Hamiltonian in the BL coordinate is described as The Hamilton equations dx i /dt = ∂H/∂p i and dp i /dt = −∂H/∂x i are then written as follows: As is mentioned above, dp ϕ /dt is not solved during the calculation of geodesics as p ϕ is conserved variable in the axisymmetric spacetime. The coefficients in the above Hamilton equations are described as follows: ∂α ∂r = 1 2α ∂ ∂r ∂ ∂θ ∂ ∂r ∂ ∂θ ∂ ∂r ∂ ∂θ 1 ρ 2 = 2 ρ 4 a 2 sin θ cos θ.
The above equations are integrated using an 8th-order embedded Runge-Kutta method with the adaptive stepsize control (e.g., Press et al. 2002).

RADIATIVE TRANSFER (I): OBSERVER-TO-EMITTER ALGORITHMM
The covariant radiative transfer equations, which are expressed by a suitable formula for the time-reversed calculations, are solved when we use the observer-toemitter algorithm.
The covariant radiative transfer equations without scatterings are described as following: where τ ν is the optical depth for photons with their frequency ν measured in the observer frame. The invariant specific intensity I and source function S is wrtten as where J and A are the invariant emissivity and absorption coefficient, respectively. The covariant radiative transfer equation (25) leads to the equation which is convenient to be integrated from the observer to emitter (Younsi et al. 2012;Pu et al. 2016): where τ ν(a) is the optical depth for absorption calculated from the observer to emitter. In practice, the infinitesimal changes of the invariant quantity τ ν(a) calculated in our code such as where dl is the distance of photon propagation measured in the zero angular momentum observer (ZAMO; Bardeen et al. 1972) frame. In the ZAMO frame, the infinitesimal dxμ can be transformed from that in the BL coordinate dx µ as where Then, dl is obtained by dl = (dxr) 2 + (dxθ) 2 + (dxφ) 2 . Since the radiative transfer equation is not scale free, we should calculate dl in the physical unit, i.e., cgs unit in the code, while the geodesic equations can be calculated with being normalized by the gravitational radius. We note that ν f and ν z are the photon frequency observed in fluid-rest and ZAMO frame, respectively. The covariant 4-momentum of photon in the ZAMO frame is calculated as where p µ is the covariant 4-momentum in the BL coordinates and The covariant source function Eq. (27) is calculated by using total emissivity and absorption coefficient described in the following subsection.

Emission and Absorption
Various emission mechanism can be simultaneously incorporated by summing up each emissivity in the fluidrest frame: where the detailed formulation of the emissivity of each radiative processes via thermal and nonthermal electrons are described in Appendix B.1 The absorption coefficient can be obtained by summing up those for each radiative process: where the detailed description of the absorption coefficients for the radiative processes via thermal and nonthermal electrons are shown in Appendix C. We integrate the covariant radiative transfer equation from the observer to the emitter described as Eq (30), by substituting the total emissivity and absorption coefficient to the covariant form of radiative quantities desribed in Eq (27), (28), and (29).

Observer Screen
The observer screen is located at the distance far from the black hole, e.g., 10 4 r g in default. Each ray incoming to the screen is assumed to be perpendicular to the surface of the screen. We inversely trace the photon which arises from the center of the pixel of the screen. We solve the radiative transfer with various photon frequency at observer screen simultaneously, during the ray-tracing calculations.

RADIATIVE TRANSFER (II): EMITTER-TO-OBSERVER ALGORITHM
The emitter-to-observer algorithm with a Monte Carlo method is used to transport the photons taking into account the effects of Compton/inverse-Compton scattering. For the Monte Carlo calculations, pseudo random numbers are generated by using Mersenne-Twister method.
The four reference frames are used in this method: (i) observer frame (BL frame) is used for the ray tracing, (ii) ZAMO frame is basically used for analysis, (iii) fluid-rest frame is used for the calculation of emission and absorption coefficient, and (iv) electron-rest frame is used for computation of Compton scattering. For the transformation among (ii), (iii) and (iv), the Lorentz transformation is used because these are the tetrad frames.

Emission
To transport the photons with the Monte Carlo method, computational particles representing the photons so called superphotons, are generated in the computation box. The weight of superphotons w, which is a dimensionless quantity, represents the ratio of number density of real photons to superphotons in the phasespace generated in the unit time (e.g., Dolence et al. 2009): where N real and N super are the number of real photons and superphotons, respectively. Therefore, w can be described as whereN super is the number of superphotons per photonfrequency bin per computational cell. It should be noted that √ −g, ∆t, ∆ 3 x dν are defined in the observer frame, while dΩ (f) is defined in the fluid-rest frame 1 .

Absorption
The absorption processes reduce the weight of the superphoton as where w new is the updated value of superphoton weight w after one tiemstep, and ∆τ ν(a) is the optical depth for absorption per timestep of ray tracing. The total absorption coefficient described by (37) and formulation of the optical depth (31) are used to evaluate ∆τ ν(a) .

Scattering
The overview of the procedure for calculating the Compton scattering is as follows: (i) we calculate the optical depth for the Compton scattering dτ ν(s) between the current and previous position of the superphoton, (ii) we evaluate the probability of Compton scattering from dτ ν(s) , and examine whether the scattering takes place at the current position by using a Monte Carlo method, (iii) the Compton scattering with the Monte Carlo method is calculated if the scattering occurs, and (iv) we go back to (i) after integrating the geodesic equation for one timestep. Our procedure for the solving Compton scattering with a Monte Carlo method is similar to that of Dolence et al. (2009), except that we incorporate the effect of the nonthermal electrons in addition to the thermal electrons.

Optical Depth for Compton Scattering
The optical depth for the Compton scattering is computed taking into account the fluid motion and kinetic motion of thermal and nonthemal electrons: where σ KN,th and σ KN,nth are the Klein-Nishina cross section modified by taking into account the effect of kinetic motion of thermal and nonthermal electrons in the fluid-rest frame, respectively. We note that the former is the same as the hot cross section in Dolence et al. (2009). These modified cross section is formulated as: where dn e,th /dγ e is the Maxwell-Jüttner distribution function described in eq. (A1), dn e,nth /dγ e is the distribution function of nonthermal electron assumed to be isotropic in the momentum space, and σ KN is the Klein-Nishina cross section defined in the electron-rest frame, where e is the photon energy, which is normalized by the rest-mass energy of electrons, in the electron-rest frame.
We tabulate the modified cross section for thermal and nonthermal electrons before the GRRT calculations. It should be noted that the nonthermal electrons are divided to subgroups by their Lorentz factor in the fluid rest frame, which is explained in the next subsection.

Scattering Probability
It should be determined whether the Compton scattering occurs at each timestep of integration of geodesic equation. In order to accelerate the calculation, it is useful to introduce the the bias parameter b ): Due to the introduction of the bias parameter b, each superphoton is split into two superphotons with dividing the weight w, when it is determined to be scattered. The weight of the scattered photons w s and the the other one are described as follows: so that the number of photons is conserved via the scattering processes with the split of superphotons 2 . When we include the nonthermal electrons, we further extend this method. At first, we set the bias parameter b for nonthermal electrons as follows: where b th /b nth,i and ∆τ ν(s)th /∆τ ν(s)nth,i are the bias parameter and the optical depth for the Compton scattering for thermal/nonthermal electrons, respectively. The subscript "i" identifies a subgroup of nonthermal electrons, which is introduced to accelerate the calculation. In astrophysical plasmas, the nonthermal electrons with power-law distribution γ −p e usually has the power-law index with p > 0, i.e., the number density of nonthermal electorns with high Lorentz factor is very small. It requires, therefore, a huge computational time to sample the nonthermal electrons with high Lorentz factor with Monte-Carlo method. In order to avoid this problem, we divide the energy distribution of nonthermal electrons into subgroups with their Lorentz factor.
We summarize our practical method to determine whether the scattering occur at each timestep and which of the thermal and nonthermal electrons are sampled if the scattering takes place. We tabulate the modified cross section for each subgroup of the nonthermal electrons: as well as that for thermal electrons (which are not divided to subgroups). Then the optical depth for the Compton scattering is calculated as The scattering probability is calculated as We generate a pseudo-random number ξ 1 to determine whether the superphotons are scattered at each timestep. If the following condition is satisfied, the superphoton is scattered: When this condition is satisfied, we examine which of the thermal and nonthermal electrons scatters the superphoton. We assume that the thermal electrons scatter the superphotons if the following relation is satisfied: where, ξ 2 is a newly generated random number. The nonthermal electrons in the subgroup i = i (< i max ) is selected when the following condition is satisfied: In the other case, the nonthermal electrons in the subgroup i = i max is selected.

Solving Compton Scattering
When the condition for the scattering is satisfied, we solve the Compton scattering by using a Monte Carlo method. First of all, (i) a scattering electron in the fluidrest frame is sampled, and then, (ii) the four-momentum of the scattered photon is determined using another pseudo-random number.
Before the Compton scattering, the scattering electrons are sampled as follows: 1. The Lorentz factor of thermal or nonthermal electrons in the fluid-rest frame is determined. For thermal electrons, the Lorentz factor of the electrons is sampled to reproduce the Maxwell-Jütter distribution (eq. A1). We employ a procedure based on a rejection method described in Canfield et al. (1987), which is summarized in Appendix A.2.1. On the other hand, for nonthermal electrons with a single/broken power law distribution, the Lorentz factor of the scattering electrons can be sampled much easier by the inverse function method descibed in Appendix A.2.2 and A.2.3.
2. The direction of the velocity of the scattering electron, i.e., the angle between the momenta of scattering electron and the incident photon, is determined. When we sample the scattering electrons, the relativistic motion of electron should be taken into consideration during the MC sampling. This is because (i) the scattering cross section in the fluid-rest frame is (1 − µ e β e )σ KN and (ii) σ KN is significantly reduced when the Doppler-boosted photon energy in the electron-rest frame is significantly larger than the rest-mass energy of electrons (Gorecki & Wilczewski 1984). In order to take into account these effects, we follow the algorithm in Canfield et al. (1987) as follows. The probability distribution function is anisotropic by the factor (1 − µ e β e ). The cosine of the angle between the incident photon and the scattering electron µ e is then sampled by the inverse function method using a pseudo random number ξ 3 : There still remains the degree of freedom of the other angle of four-momentum of the electron, and this angle is sampled by a uniform random number. At this step, the four-momentum of the electron is tentatively determined.
3. Finally, we examine the validity of the sampled electrons with their four-momentum using a rejection method. We accept the sampled electron if it satisfies where σ KN shown in equation (44) is evaluated by the photon energy in the electron-rest frame e , which is calculated by the Lorentz transformation of the four-momentum of the incident photon from the fluid-rest frame to the electron-rest frame of the sampled electron.
Next, we sample the four-momentum of the scattered photons. As described above, the four-momentum vector of the incident photons is Lorentz transformed from the fluid-rest frame to the electron-rest frame. Using a pseudo-random number ξ 5 , a tentative value of the scattered photon energy can be represented as where e is related to e and the cosine of the angle between the electron and the scattered photon µ e = 1 + 1/ e + 1/ e , so that e,min = (1 + 2 e )/ e and e,max = e . Once e is determined, we obtain µ e . The remaining angle, whose probability distribution function uniformly distributes [0, 2π), is determined by generating another pseudo-random number. Finally, in order to examine whether the tentative value of e is presumable, we accept the value which satisfies the following relation: If the tentative e is rejected, the procedure for sampling the scattered photon is repeated until it is accepted through the rejection method. Once the sampled scattered photon is accepted, the four-momentum of the scattered photon in the electron-rest frame is converted to that in the fluid-rest frame and in the ZAMO frame via the Lorentz transformation. The covariant four-momentum of photon in the BL frame is also calculated via the coordinate transformation in order to solve the geodesic equations at the next time step.

Calculation of SED
The SED is calculated after the superphoton escape to the outer boundary of the computational domain. We divide the bins of the photon-frequency in the observer frame ν and the viewing-angle θ ob in such a way that ∆ ln ν and ∆ cos θ ob are constant. Here viewing angle θ ob is defined as the angle between the rotation axis of the Kerr BH and the line-of-sight of observers.
For an observer at jth viewing-angle bin, the value of νL ν at ith frequency bin can be calculated as follows: (63) where w n (i, j) is the weight of n-th photons leaving the computational domain from the outer boundary with belonging to the i-th frequency bin and the j-th viewingangle bin. We set ∆Ω = 2π∆ cos θ ob and ∆t = 1 in default.

Trajectory of photon
We test the geodesics solver by comparing the numerical results and the analytic solution on the equatorial plane of a Kerr BH for the case with p ϕ = −ap t (Chandrasekhar 1983): where r + and r − are the outer and inner horizon of the Kerr BHs r ± = M ± √ M 2 − a 2 . Figure 1 shows the result with setting a * = 0.998, where a * ≡ a/M . The BH rotates in the counterclockwise direction. The photons are initiated at r = 10r g , θ = π/2, and ϕ given by eq. (64). The four-momentum of photon is p t = −E, p θ = 0, p ϕ = a * E, and p r is derived from the the constraint p µ p µ = 0. As the photon gets close to the BH, its trajectory is strongly affected by the BH gravity and the frame-dragging effect. The result computed by our code (red points) well reproduces the analytic solution (black line).

Black Hole Shadow
As another tests of the geodesic equation solver, the shape of the BH shadows are calculated with the viewing angle 30 • and 90 • for a * = 0, 0.5, and 0.998, and shown in Figure 2. The photon trajectory is solved from the observer screen with the time-reversing way. The black color displays the region where the photon is captured by the black hole, while the light-gray color represents the area where the photons escape to infinity. The dotted, orange lines represent the shape of black-hole shadow given by an analytic solution, whose radius and angle in the observer screen is described as follows (Bardeen 1973;Johannsen 2013;Wong 2021): general relativistic, multi-wavelength radiative transfer . where x and y are the cartesian coordinate on the observer screen in which the BH spin vector projected onto the screen is on the y-axis, and θ ob is the viewing angle of the observer. The variables ξ and η are parameterized by the radius of the corresponding photon orbit r ph as follows: We note that r ph continuously exists between their minimum and maximum vales r ph,min and r ph,max , which are the roots of the equation where the region formed by these unstable circular photon orbits is so called photon shell (Johnson et al. 2020; Event Horizon Telescope Collaboration 2019e; Kawashima et al. 2021). Our numerical calculations successfully reproduce the analytical shape of the BH shadow for the various inclination angle and BH spin. We note that, as the BH spin increases, the location of the BH shadow shifts in the direction perpendicular to the BH spin vector projected onto the observer screen, because of the frame-dragging effect (for the application to a possible estimate of the BH spin, see Kawashima et al. 2019;Chael et al. 2021).
More tests of computations of geodesic equations and black hole shadows using RAIKOU are also carried out by GRRT code comparison projects of EHT (Gold et al. 2020).

Compton scattering via thermal electron
In this section, we test the Compton scattering of blackbody photons via thermal electrons in a uniform plasma sphere. For the sake of the simplicity, the spacetime is assumed to be flat (i.e., the Minkowski spacetime) and the absorption processes are ignored here. The radiative temperature of blackbody photons T rad is set to be k B T rad /m e c 2 = 10 −8 , where k B , m e and c are the Boltzmann constant, electron rest mass, and speed of light, respectively. The blackbody photons are injected from the inner sphere at r = r in = 1 cm. The outer boundary of the computational domain locates at r = r out = 10 5 cm, so that the finite radius of r in negligibly affect the results.
We examine the Compton scattering via thermal electrons as in Dolence et al. (2009). The temperature and number density of electrons are set to be Θ e = k B T e /m e c 2 = 4 and n e,th = τ (s) /σ T (r out − r in ), respectively. Here, τ (s) is the optical depth for Thomson scattering and we set τ (s) = 10 −4 , 10 −1 , and 3. Figure 3 displays the resulting SEDs, which reproduces the results shown in Dolence et al. (2009). In Fig 3, the SEDs are also decomposed into those of photons without being scattered, with being scattered once, twice, and three or more times. The optically thin cases (i.e., τ (s) = 10 −4 and 10 −1 ) more clearly show the feature of individual Comptonization processes. Spectral bumps at ν ∼ 10 15 , 10 17 , and ∼ 10 19 Hz are formed via the inverse-Compton scattering, which is consistent with the fact that the change of photon frequency per scattering ∆ν is ∆ν/ν ∼ 16Θ 2 e ∼ 10 2 (e.g., Rybicki & Lightman 1979). One can also find that the L ν at each spectral bump reduces by ∼ 10 −4 and ∼ 10 −1 for the case with τ (s) = 10 −4 and 10 −1 , respectively, since the L ν at the frequency of scattered photon reduced by the scattaring probability 1−exp(−τ (s) ) ∼ τ (s) in the optically thin limit. In the optically thick case τ (s) = 3, the number of scattering is ∼ 10 and the Compton y-parameter is significantly greater than unity y = 16Θ 2 e max(τ (s) , τ 2 (s) ) ∼ 2 × 10 3 1, which results in the formation of the Wien-like peak around 3k B T e /h ∼ 10 21 Hz. The resulting SEDs therefore, well reproduce the spectral features of Comptonization both in the opticallty thin and thick plasmas.

Compton scattering via nonthermal electron
We examine the Comptonization of blackbody photons via nonthermal electrons with single and broken power-law distribution functions in the uniform spherical plasma. The radial optical depth against the Thomson scattering is 10 −4 . In this subsection, we show the SEDs of νL ν rather than L ν to clearly present the high energy photons Compton-upscattered by electrons with a high Lorentz factor. Figure 4 displays the SEDs of photons Comptonized via nonthermal electrons with a single power-law distribution n e,nth ∝ γ −p e for γ e,min ≤ γ e ≤ γ e,max , where γ min = 2 and γ max = 2 × 10 4 . The SED peak of photons, which experienced one scattering event, appears at ∼ 10 21 Hz. This peak is formed by inverse-Compton scattering with electrons whose Lorentz factor γ e γ e,max , i.e., the photons with the peak frequency of injected blackbody photons ν bb 4 × 10 12 Hz are upscattered to the photon frequency γ 2 e,max ν bb 10 21 Hz. The additional peak at ν ∼ 10 23 Hz is formed by two scattering events, where the peak frequency appears at γ e,max m e c 2 /h ∼ 10 24 Hz because of the reduce of the scattering cross section (i.e., Klein-Nishina effect). Figure 5 shows the Comptonized SEDs with the same model setup except for the higher maximum Lorentz factor of nonthermal electrons γ e,max = 2 × 10 6 . In the gamma-ray energy bands, only one remarkable spectral peak appears at ν ∼ 10 25 Hz, although two peaks appear in the gamma-ray energy band for the model with the lower maximum Lorentz factor (see Fig. 4). This is because the spectral peak formed by the one-scattering component γ 2 e,max ν bb ∼ 10 25 Hz and the possible maximum peak of Comptonized SEDs γ e,max m e c 2 /h ∼ 10 26 Hz are similar. For the component with scattering twice or more, the Compton scattering with high Lorentz fac- Figure 3. Comptonized SEDs of blackbody photons, which are scattered by the hot electrons Θe = 4 in the plasma sphere with optical depth for Thomson scattering τ (s) = 10 −4 (top), 10 −1 (middle), and 3(bottom). The resulting SEDs are displayed with the red solid lines. Black lines represent the decomposed SEDs, whose photons are escaped towards the observer without being scattered (dashed), with being scattered once (dotted), twice (dash-dot-dash), and more than three times (dash-dot-dot-dash).
tor of electrons are highly suppressed by the Klein-Nishina effect and the spectral peak appears at ν γ e,max m e c 2 /h ∼ 10 26 Hz. . Comptonized SEDs of blackbody photons, which are scattered by nonthermal electrons with power-law index p = 2.5 in the plasma sphere. The minimum and maximum Lorentz factor of nonthermal electrons are γe,min = 2 and γe,max = 2 × 10 4 , respectively. The optical depth for Thomson scattering is τ (s) = 10 −4 . The resulting SEDs are displayed with the red solid lines. Black lines represent the decomposed SEDs, whose photons are escaped towards the observer without being scattered (dashed), with being scattered once (dotted), twice (dash-dot-dash), and more than three times (dash-dot-dot-dash). In figure 6, we show a result of Comptonized SEDs calculated for broken power-law model of nonthermal electrons. The power-law distribution extends from γ e,min = 2 × 10 2 to γ e,max = 2 × 10 4 and the powerlaw index is broken at γ br = 2 × 10 2 . The power-law indices are p 1 = 2.5 and p 2 = 3.5 in γ e,min ≤ γ e ≤ γ e,br and γ e,br < γ e ≤ γ e,max , respectively. As a consequence of the break of power-law electron destribution Figure 6. The same as Fig. 4, except for broken power-law distribution of nonthermal electrons. The power-law indices are p1 = 2.5 and p2 = 3.5 between 2 ≤ γe ≤ 2 × 10 2 and 200 < γe ≤ 2 × 10 4 , respectively.

Synchrotron Self-Compton
In this subsection, the synchrotron self-Compton processes via nonthermal electrons are tested. We set a plasma sphere with uniform spatial distribution of number density of nonthermal electrons n e,nth = 10 3 cm −3 and magnetic-field strength B = 10 −2 G, where the orientation of the magnetic field is isotropic. The radius of the plasma sphere is set to be 10 15 cm. The bulk speed of the plasma sphere is assumed to be zero. We examine the models with single-and broken-power law distribution function of nonthermal electrons.
First, we present the results for the model with a single power-law distribution of electrons, where the distribution function is shown in eq. A2. The power-law indices is p = 2.5 for the Lorentz factor of nonthermal electrons γ e,min ≤ γ e ≤ γ e,max , where γ e,min = 10 2 and γ e,max = 10 7 . Figure 7 shows the resulting SED calculated by RAIKOU (red line). As a reference, the SED calculated by a numerical code without based on MC methods (AGN SED tool, Massaro et al. 2006;Tramacere et al. 2009Tramacere et al. , 2011, in which one-zone models are assumed, is also plotted as black filled circles. It is found that our result well agrees with the reference. The spectral features in Fig. 7 is briefly mentioned here. The SED is steep in ν 10 9 Hz because of the synchrotron self-absorbed process. In 10 9 ν 10 17.5 Hz, the SED formed by the optically-thin, synchrotron Figure 7. SEDs of synchrotron Self-Compton via nonthermal electrons with a single power-law distribution. The power-law indices are p = 2.5 for the electrons with the Lorentz factor γe,min ≤ γe ≤ γe,max. Here, we set γe,min = 10 2 and γe,max = 10 7 . The red line and the black circles present the results obtained by RAIKOU and a public code AGN SED tool, respectively.  Fig. 7, but for a broken power-law distribution. The power-law indices are p1 = 1.5 and p2 = 3.5 for the electrons with the Lorentz factor γe,min ≤ γe ≤ γ e,br and γ e,br < γe ≤ γe,max, respectively. Here, we set γe,min = 10 2 , γ e,br = 10 4 , and γe,max = 10 7 . emission via the nonthermal electrons with p = 2.5, i.e., νL ν ∝ (3 − p)/2 = 0.25 appears. As is well known, the SEDs due to the synchrotron emission via electrons with a power-law distribution is the same as those due to the inverse Compton scattering of photons via nonthermal electrons with the same distribution function (e.g., Rybicki & Lightman 1979). In ν 10 19 Hz, the spectral bump resulted by the inverse Compton scattering of synchrotron photons via nonthermal electrons (i.e., synchrotron self-Compton) appears. Here, it should be noted that the SED formed by the synchrotron self-Compton processes show a bump rather than the simple power-law shape, because the seed photons (i.e., synchrotron photons) have the broadband SEDs which consequently result in the deviation from the simple powerlaw SED of the Comptonized photons.
Next, we present the results for the model with a broken power-law distribution of electrons, where the distribution function is shown in eqs. (A4) and (A5). The power-law indices are p 1 = 1.5 and p 2 = 3.5 for γ e,min ≤ γ e ≤ γ e,br and γ e,br ≤ γ e ≤ γ e,max , respectively, where γ e,min = 10 2 , γ e,br = 10 4 , and γ e,max = 10 7 . Figure 8 shows the resulting SED. As is the calculation with single power-law model, our computed SED well reproduce that calculated by AGN SED tool. Here, we briefly note the spectral features of calculated SED in Figure  8 . As is the same as the single power-law model, the SED is steep in 10 9 Hz because of the synchrotron self-absorbed process. In 10 9 ν 10 12.5 Hz and 10 12.5 ν 10 17 Hz, the optically-thin, synchrotron emission via the nonthermal electrons with p 1 and p 2 appears, respectively. The spectral indices of νL ν are (3 − p 1 )/2 = 0.75 and (3 − p 2 )/2 = −0.25, respectively. In ν 10 17 Hz, the spectral bump resulted by the synchrotron self-Compton via the nonthermal electrons appears.

APPLICATION TO ACCRETING KERR BLACK HOLE BASED ON GRMHD SIMULATIONS
We show examples of full code tests with applying RAIKOU to the calculation of images and SEDs of the accretion flows, wind, and jets of GRMHD simulation data. We demonstrate the results of a radiatively insufficient accretion flow around a SMBH where the accretion flows are simulated with by using a GRRMHD code UWABAMI (Takahashi et al. 2016). The key parameters for the GRRT calculations are described bellow and summarized in Table 1. The GRMHD simulation data are briefly summarized in appendix D.
We consider the synchrotron emission/abosrption and Compton/inverse-Compton scattering processes via thermal and nonthermal electrons. The dimensionless spin parameter a * is set to be 0.9375. We need to give the scale of the system, because GRRT calculations are not scale free. Here, we set the parameters suitable for one of the main targets of EHT, Sgr A*: the BH mass M BH = 4.1 × 10 6 M , the mass accretion ratė M ∼ 10 −8 M y −1 , and Distance of the observer screen for image calculation D = 8.1 kpc. We set the black hole mass, accretion rate, and distance suitable for one of the main targets of EHT, Sgr A*: The field-of view (FoV) of the screen is assumed to be (10 3 µas, 10 3 µas), which is divided by (2.5 × 10 3 , 2.5 × 10 3 ) pixels.
Because of the uncertainty of the electron temperature in the GRMHD simulations, we use a formula describing T p /T e as a function of the plasma β, which is the same as Mościbrodzka et al. (2016); Event Horizon Telescope Collaboration (2019e): where R high = 5 and R low = 1. The distribution function of nonthermal electrons is assumed to be a single power-law with a power-law index p = 3.5. The minimum and maximum Lorentz factor of the nonthermal electrons are γ e,min = 30 and γ e,max = 10 6 , respectively. The number density of nonthermal electrons are set in such a way that their energy density u e,nth is proportional to that of the magnetic energy density u B at each computational cell: where a efficiency parameter η nth is set to be 0.03 and to be uniform in space in this work. Figure 9 demonstrates images of accretion flow and jet at 22 (left), 86 (center), and 230 GHz (right), and the broadband SEDs from radio to VHE gamma-ray. The top/middle panels show the images with/without taking into account the effects of nonthermal electrons. The size of bright region becomes large as a consequence of the inclusion of nonthermal electrons, which is consistent with the previous interpretation of the Sgr A* using simple toy models (e.g., Özel et al. 2000;Cho et al. 2021). This is because the emissivity via the synchrotron processes is still high in the outer region, while the plasma temperature is low and the thermal synchrotron emission is very faint there. Although the effect of nonthermal electrons is less significant at higher frequency in the radio band, the image taking into account the nonthermal electorns still shows a slightly extended emission image including compared with that without nonthermal electrons at 230 GHz (i.e., at the photon frequency observing black hole shadows by EHT) 3 . Note-A set of parameters for GRMHD simulations is summarized in Appendix D.
The broadband SEDs shown in the bottom panel in Figure 9 present significant differences between the models with and without nonthermal electrons. As is described above, the luminosity in the lower radio band 100 GHz is increased via the nonthermal synchrotron emission. Importantly, the significant amount of X-ray photons are generated by the synchrtron emission via nonthermal electrons in addition to the inverse Compton scattering via the thermal electrons. The SED with thermal and nonthermal electrons (red solid line) at ∼ 10 17−20 Hz is the superposition of the nonthermal synchrotron (red dashed line) plus the SSC via thermal electrons (black solid line). In the gamma-rays band 10 20 , the significant amount of the photons up to ∼ TeV are generated by the inverse-Compton scattering due to the nonthermal electrons (red solid line), while a poor amount of gamma-ray photons with ν > 10 22 are emitted in the model with thermal electrons only (black solid line). The power-law shaped SED with the spectral index ∼ (3 − p)/2 = −0.25 expected for synchrotron emmission and inverse-Compton scattering via nonthermal electron appears in the photon-frequency range 10 16 ν 10 22 Hz and 10 23 ν 10 25 Hz, respectively, while the SED shows slightly complicated structure because of the contribution of the thermal electron processes.
We check the consistency between the two algorithms, i.e., observer-to-emitter and emitter-to-observer algorithms. Figure 10 shows the SEDs calculated by using the observer-to-emitter and emitter-to-observer algorithms, where the Compton effects are not included for the comparison because the observer-to-emitter algorithm cannot solve the Compton scattering. The SEDs with the former algorithms are averaged in the azimuthal direction to compare with the latter ones, since our broadband SEDs are computed by summing up over the azimuthal angle ϕ and averaged as in grmonty ). One may find that the SEDs computed by the different methods are mutually-consistent, except a small difference in ν ∼ 10 15−16 Hz, i.e., the SED are slightly deviated due to the relativistic Doppler effects caused by the relativistic bulk motion in the jet between 10r g ≤ r < 30r g . It should be also noted that the deviation is still (several) × 10% at highest in the all frequency range of the SEDs, which will be enough to compare the theoretical SEDs with observation data. The radio emission, which is crucial for the comparison with the EHT data, is nearly identical between the two different calculation algorithms.
The results of the self-convergence test for our broadband SEDs calculations are shown in Figure 11. We used the SEDs with N 4.9 × 10 11 shown in Figure 9 for the references, where N is the total number of superphotons. As is expected, the relative errors accumulated over the photon-frequency bins decrease roughly as N −1/2 for both of the models with and without the nonthermal electrons.

SUMMARY AND DISCUSSION
We have developed a multi-wavelength, general relativistic radiative transfer code RAIKOU (来光) to compute the multi-wavelength images and SEDs of accreting Kerr BHs in the broadband energy bands from radio to VHE gamma-ray, based on a ray-tracing method. The two different algorithms are implemented to RAIKOU; an observer-to-emitter algorithm for efficient calculations of images and an emitter-to-observer algorithm based on a Monte-Carlo method for computations of broadband SEDs by incorporating the inverse-Compton scattering of soft photons via themal and nonthermal electrons.
The radiative processes of cyclo-synchrotron emission/absorption via thermal electrons, synchrotron emission/absorption via thermal and nonthermal electrons, bremsstrahlung emission/absorption via thermal electrons, and inverser-Compton/Compton scattering via thermal and nonthermal electrons are incorporated. The implementation of these important radiative processes to the code enabled us to calculate the detailed images and broadband SEDs with the broad energy band ranging from radio to VHE gamma-ray. The nonthermal synchrotron emission will enlarge the bright region in the resulting images, and X-ray and gammaray emissions are significantly increased due to the combination of the nonthermal synchrotron emission and the inverse-Compton scattering via the nonthermal electrons, respectively. These will be useful to constrain the physical parameters in the main target of EHT, i.e., Sgr A* and M87, and other accreting black holes. Al-   though we focused on the application to the calculations for radiatively inefficient accretion flow in this paper, our code can be applicable to the studies on the super-Eddington accretion flow around black holes and neutron stars, which will be realized in ultraluminous Xray sources and tidal disruption events which we leave as future works.
In the current version of our code, the e ± creation and annihilation via the γ-γ collision, which will be im-portant to study the formation of the e ± jet ) and the SEDs in VHE gamma-ray bands absorbed by the intrinsic photons and the extragalactic background light (EBL). While the latter effect, i.e., the absorption of VHE gamma-ray via the pair creation due to the interaction with EBL, is negligible to study the Galactic Center and nearby low luminosity AGN, which are the main targets of EHT (e.g., Sgr A* and M87), it is important to be incorporated when we extend our study to the SEDs of distant AGNs (e.g., MAGIC Collaboration et al. 2008). It should be also mentioned that we focused on the leptonic radiative mechanism. The incorporation of the hadronic processes, e.g., the p-p and p-γ collision and successive production of gamma-ray photons via the decay of pions, as well as the production of high energy neutrinos via the decay of the charged pions. The direct synchrotron emission via accelerated protons will also contribute to the VHE gamma-ray emissions. In our demonstration, the fast-light approximation (i.e., the snapshot data is fixed during the radiative transfer calculation) is used. For more precise calculations and the study of the time-variable images and SEDs, it will be important to use the snapshot depending on the physical time of during the GRRT computation. The inclusion of these effects will be addressed in the near future.
We thank M. Kino, I. Cho, K. Kawaguchi, and K. Asano for fruitful discussion. The numerical simulations were carried out on the XC50 at the Center for Computational Astrophysics, National Astronomical Observatory of Japan. This work also used computational resources of the supercomputer Fugaku provided by RIKEN through the HPCI System Research Project (Project ID: hp120286). This work was supported by JSPS KAKENHI grant Nos. JP18K13594, 19H01908, 19H01906 (T.K.), 18K03710, 21H04488 (K.O.), 20K11851, 20H01941, 20H00156 (H.R.T). This work was also supported in part by MEXT SPIRE, MEXT as "Priority Issue on post-K computer" (Elucidation of the Fundamental Laws and Evolution of the Universe) and as "Program for Promoting Researches on the Supercomputer Fugaku"(Toward a unified view of the universe: from large scale structures to planets), and JICFuS. The distribution function of relativistic thermal electrons, i.e., the Maxwell-Jüttner distribution, is described as: where n e,th is the number density of thermal electrons, γ e is the Lorentz factor of the thermal electron motions, K 2 is the modified Bessel function of the second kind, Θ e = k B T e /m e c 2 is the dimensionless electron temperature, T e is the electron temperature, k B is the Boltzmann constant, m e is the electron mass, and c is the speed of light.

A.1.2. Single Power Law Distribution
The single power-law distribution dn e,nth /dγ e ∝ γ −p e between γ e,min ≤ γ e ≤ γ e,max is written as Using the same procedure as Canfield et al. (1987), we sample the electrons with the Maxwell-Jütter Distribution as following: 1. We generate a pseudo random number ξ 1 to select i(= 3, 4, 5, 6) of the function π i in the following equation: π 6 (Θ e ) = 1 S 3 (Θ e ) where S 3 (Θ e ) is described as follows.
Obviously, 6 i=3 π i (Θ e ) = 1 and i can be determined by using the pseudo random number 0 ≤ ξ 1 ≤ 1. 2. Using the chosen i, we generate the random number ξ 2 reproducing g i , which is the χ 2 distribution of order i where ξ 2 distributes −∞ < ξ 2 < ∞.
3. The generated ξ 2 is acceptable if it satisfies where ξ 3 (0 ≤ ξ 3 ≤ 1) is a newly generated random number for the rejection method. The Lorentz factor of electron is given by If the above condition is not satisfied, then we go back to the step 1. The procedure repeats until the condition is satisfied.

A.2.2. Single Power-Law Distribution
On the other hand, for nonthermal electrons with a single/broken power law distribution, the Lorentz factor of the scattering electrons can be sampled much easier by the inverse function method. For a single power law distribution, the Lorentz factor can be represented by generating a pseudo random number ξ: γ e = γ min 1 + ξ γ e,max γ e,min γ e = γ min γ e,max γ e,min ξ (for p = 1).

A.2.3. Broken Power-Law Distribution
For a broken power law distribution, the Lorentz factor can be written as follows for the case with p 1 = 1 and p 2 = 1: where ξ br ≡ A γ 1−p1 e,br − γ 1−p1 e,min For the case with p 1 = 1 and p 2 = 1, γ e = γ e,br exp γ p1−p2 e,br where For plasma with mildly relativistic temperature, it is important to consider synchrotron emission taking into account the cyclotron emission. It is convenient to use the fitting formula for the angle-averaded, cyclo-synchrotron emissivity (Mahadevan et al. 1996): M (x cs ) ≡ 4.0505 x 1/6 cs 1 + 0.4 x 1/4 cs + 0.5316 where K 2 is the modified Bessel function of the second kind, n e,th is the thermal electron density measured in the fluidrest frame, B is the magnetic field strength measured in the fluid-rest frame, Θ e (= k B T e /m e c 2 ) is the dimensionless electron temperature, T e is the electron temperature, k B is the Boltzmann constant, m e is the rest mass of electron, and e is the elementary charge.

B.1.3. Bremsstrahlung via Thermal Electron
For relativistic temperature plasma, both of the electron-proton and electron-electron bremsstrahlung are important. Here, we describe the bremsstrahlung emissivity via thermal electrons (see, e.g., Svensson 1982;Stepney & Guilbert 1983;Narayan & Yi 1995;Manmoto et al. 1997). The total bremsstrahlung emissivity j (f) ν f (brms) can be written as where f Gaunt is the Gaunt factor given by and the bremsstrahlung cooling rate per unit volume q − brms can be written as q − brms = q − brms,ep + q − brms,ee .

B.2. MC sampling
In general, there are mainly two ways to sample the photon emissions with MC methods: (i) superphotons with an uniform weight are generated in such a way that their number is proportional to the emissivity of plasma among computational mesh, or (ii) the weight of superphotons is proportional to the emissivity with generating a uniform number of superphotons among computational meshes.
We employ the method (ii) because we consider various emission processes, i.e., thermal/nonthermal synchrotorn and bremsstrahlung emission. We can simply sum up the weight of superphotons when they are generated. The direction of the generated superphotons is sampled to be isotropic using uniform pseudo-random numbers. It should be emphasized that the dependence of the synchrotron emissivity on the direction are automatically incorporated by this method, because the weight of superphotons is set to be proportional to the emissivity.
Our procedure is summarized as follows: 1. Superphotons are generated in such a way that their direction becomes isotropic in the fluid-rest frame, using uniform pseudo-random numbers.
Initially, we set an isentropic hydroequilibrium torus rotating around the Kerr BH (Fishbone & Moncrief 1976), which is emmeded in a hot, static, uniform, and non-magnetized ambient gas. The position of the inner edge and the pressure maximum of the torus are set at r = 20r g and r = 33r g on the equatorial plane, respectively. The specific heat ratio is assumed to be γ heat = 13/9. After the growth of the magneto-rotational instability (Balbus & Hawley 1991), the accretion flow is formed via the angular momentum transport. We use the snapshot at t = 1.2 × 10 4 r g /c, in which the accretion flow in the quasi-steady state. The structure of the accretion flow at t = 1.2 × 10 4 r g /c is displayed in Fig.12.
The accretion flow is roughly categorized into two states by the magnitude of their magnetization: SANE (Standard And Normal Evolution) and MAD (Magnetically Arrested Disk), which are weakly and strongly magnetized states, respectively (e.g., Narayan et al. 2012;SÄdowski et al. 2013;Tchekhovskoy 2015; Event Horizon Telescope Collaboration 2019e, and references therein). The SANE and MAD are defined by the dimensionless magnetic flux threading the event horizon, which is described as where the magnetic flux Φ BH and the mass accretion rate onto the BHṀ BH is Here, dA θϕ = √ −gdθdϕ is an area element in the θ-ϕ plane and g is the determinant of the Kerr-Schild metric. In SANE and MAD states, Φ BH ∼ 1, and ∼ 15, respectively. The normalized magnetic flux is Φ BH ≈ 5 in our simulation, so that the accretion flow is categorized to the state between SANE and MAD, which is sometimes referred as semi-MAD 4 although the magnetic flux is slightly high.