Issue 
A&A
Volume 685, May 2024



Article Number  A84  
Number of page(s)  10  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/202347821  
Published online  08 May 2024 
ARTPOL: Analytical raytracing method for spectropolarimetric properties of accretion disks around Kerr black holes
^{1}
Tuorla Observatory, Department of Physics and Astronomy, University of Turku, 20014 Turku, Finland
email: vladislav.loktev@utu.fi
^{2}
Nordita, KTH Royal Institute of Technology and Stockholm University, Hannes Alfvéns väg 12, 10691 Stockholm, Sweden
^{3}
Physics Department and Columbia Astrophysics Laboratory, Columbia University, 538 West 120th Street, New York, NY 10027, USA
^{4}
Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA
^{5}
Institut für Astronomie und Astrophysik, Universität Tübingen, Sand 1, 72076 Tübingen, Germany
Received:
28
August
2023
Accepted:
21
February
2024
Spectropolarimetric signatures of accretion disks in Xray binaries and active galactic nuclei contain information on the masses and spins of their central black holes, as well as the geometry of matter in proximity to the compact objects. This information can be extracted by means of Xray polarimetry. In this work, we present a fast analytical raytracing technique for polarized light (ARTPOL) that helps us to obtain the spinning black hole parameters from the observed properties. This technique can replace the otherwise timeconsuming numerical raytracing calculations for any optically thick or geometrically thin accretion flow. For the purposes of illustration, we considered a standard optically thick, geometrically thin accretion disk in the equatorial plane of the Kerr black hole. We show that ARTPOL proves accurate for dimensionless spin parameter a ≤ 0.94 with a speed that is over four orders of magnitude faster than direct raytracing calculations. This approach opens up broader prospects for direct fittings of the spectropolarimetric data from the Imaging Xray Polarimetry Explorer.
Key words: accretion, accretion disks / gravitational lensing: strong / polarization / methods: analytical / stars: black holes / Xrays: binaries
© The Authors 2024
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1. Introduction
The spin of a black hole (BH) is a fundamental parameter that controls the behavior of the inflowing matter in the accretion disk, along with the properties of the outflowing material, namely, the relativistic jets. The magnitude of spin determines the curvature of spacetime close to the BH and the energy dissipation profile, while also affecting the spectral energy distribution of the observed emission. The spin values for both Galactic BHs in Xray binaries and their supermassive counterparts have been probed using the distinct imprints in spectral and timing properties via the continuumfitting method and the iron line reflectionreverberation method (Miller et al. 2009; McClintock et al. 2014; Uttley et al. 2014; Reynolds 2021). The methods are based on obtaining the radius of the innermost stable circular orbit (ISCO), which, in turn, is related to the BH spin. For Xray binaries, an additional constraint on the parameters comes from the relativistic precession model (Motta et al. 2014, 2022), which links the mass, spin, and radius with the characteristic frequencies found in the Xray light curves. The statistical distribution of spins probed by these methods, as well as spin values obtained by different methods for the same source, do not match (Draghis et al. 2023). This calls for an alternative method to verify spin determination measurements.
The polarization of radiation escaping from the BH vicinity can be used as a fine tool to determine the curvature of spacetime. In this context, it has been anticipated for some time that Xray polarimetric signatures of accretion disks can carry important information about the BH spin. The launch of Imaging Xray Polarimetry Explorer (IXPE; Weisskopf et al. 2022) has opened up these exciting research possibilities. At the same time, the observed polarimetric signatures of BH Xray binaries have proved many previously proposed models to fail in describing the data (Krawczynski et al. 2022; Rodriguez Cavero et al. 2023; Ratheesh et al. 2024). Thus, our basic assumptions on the geometry and radiative mechanisms producing local spectra have required updates, however, it is often impossible to make such alterations in the datafitting models. It is therefore necessary to adopt a fast tool that relates the local spectra to the observed spectropolarimetric signatures.
The accretion disk polarization is produced by multiple scatterings in the upper layers of its atmosphere. The first predictions of the disk polarization (Rees 1975) have been made using the results of calculations in the case of pure electronscattering, semiinfinite planeparallel atmospheres (Chandrasekhar & Breen 1947; Sobolev 1949; Chandrasekhar 1960; Sobolev 1963) and were limited to Newtonian approximation. The polarization degree (PD) may be altered due to the presence of absorption effects in the atmosphere, which were considered using MonteCarlo (Lightman & Shapiro 1975) and analytical (Loskutov & Sobolev 1979, 1981) means. Further, the effects of general and special relativity introduce important modifications to both PD and polarization angle (PA; Connors & Stark 1977; Stark & Connors 1977; Pineault & Roeder 1977a,b). Aberration and light deflection lead to a rotation of the PA and alter the viewing angle of different parts of the accretion disk, affecting the PD. Additional rotation of PA along the photon trajectory is expected for the case of spinning BHs, described by the Kerr (1963) metric, thanks to the framedragging effects. Therefore, the observed spectral dependence of polarimetric signatures can act as an independent probe of the BH spin (Connors et al. 1980; Dovčiak et al. 2008).
Precise computations of the effects of general and special relativity on the polarization properties can be done using the parallel transport of the polarization vector along null geodesics, which often involves computations of the constant of motion from Walker & Penrose (1970; Connors et al. 1980; Dovčiak et al. 2008; Ingram et al. 2015). The geodesics, in turn, are computed using the raytracing techniques (e.g., Dexter 2016; Zhang et al. 2019; Chan et al. 2013; Pihajoki et al. 2018). It is also essential to keep track of the convergence of the computed flux, PD and PA down to the characteristic scale of the observed errors on these quantities – this often means that several simulations with increased resolution must be performed for each parameter set. An additional source of computational errors comes from the assumptions of a small outer radius of the disk, which is enforced by the high computational costs. Typically, the value ≲100 R_{S} (where R_{S} is the Schwarzschild radius) is considered, which allows for the (highly polarized) secondary images of the disk to be visible in the region surrounding the outer radius of the simulated disk; in reality, a much higher extend of the disk, ≳10^{5}R_{S}, completely covers those from the line of sight.
The raytracing technique is too computationally expensive for a direct data fitting. Instead, precomputed geodesics have been used to accelerate calculations (e.g., Li et al. 2005; Krawczynski 2012; Krawczynski & Beheshtipour 2022). Alternatively, analytical and semianalytical approaches can be applied to solve geodesic equations (Dexter & Agol 2009; Yang & Wang 2014; CárdenasAvendaño et al. 2023). Their applicability is, however, limited; for instance, the semianalytical expressions for geodesics can only be used for the equatorial plane of the BH. Finally, the calculation of geodesics may be omitted in the Schwarzschild metric and the photon trajectory can be treated using an approximation to the lightbending relation (Beloborodov 2002; Poutanen 2020a). The latter approach gave reliable results for the lowenergy synchrotron emission observed from the supermassive BH M87*, with an almost faceon disk (Narayan et al. 2021).
The physical understanding of the modifications to polarization signatures caused by the curvature of spacetime and fast motions of matter is difficult to achieve if we use the implicit WalkerPenrose constant. For this purpose, we can use the explicit analytical expression for the rotation of the polarization plane along the photon path, which we call the analytical raytracing technique for polarized light (ARTPOL hereafter). This approach has been used to extract polarimetric properties of spinning spherical (Poutanen 2020b) and oblate (Loktev et al. 2020) neutron stars, as well as accretion disks around Schwarzschild BHs (Loktev et al. 2022). This method was first proposed in Pineault (1977) and applied, in the context of accretion disks, to the PA rotation caused by general and special relativity individually, while the expression for their combined effects was first derived in Loktev et al. (2022).
In this work, we apply ARTPOL to the accretion disks around spinning BH. In Sect. 2, we describe the formalism that can be used in spectropolarimetric modeling and imaging of accretion disks around Kerr BHs. In Sect. 3, we compare the results of the ARTPOL technique to those obtained using explicit raytracing calculations. We show that the PD and PA computed using this method remain accurate to the level of current observational (IXPE) uncertainties for BH spin parameters up to a = 0.94. We summarize our findings and discuss a broad range of applications of the technique in Sect. 4.
2. Methods
In this section, we describe the procedure for computing polarization spectra from a standard optically thick, geometrically thin accretion disk in the equatorial plane of a Kerr BH. This example is selected for our analysis as it represents one of the most commonly assumed scenarios for Kerr BH accretion and has been chosen for its simplicity. The method for acquiring polarization parameters described here is generalizable and adaptable to different flow geometries using the vector formalism described in Loktev et al. (2022).
2.1. Local model of the disk emission
We considered a BH with mass M and dimensionless spin a = Jc/GM^{2}, where J is the angular momentum. For the sake of simplicity, we considered a standard equatorial geometrically thin accretion disk (Novikov & Thorne 1973). The BH is situated at the origin with the spin directed along the zaxis and is orthogonal to the disk plane. We assume that the state of the fluid only depends on the BoyerLindquist dimensionless radial coordinate, r = R/R_{S}, which expresses the distance from the central object, R, in units of the BH Schwarzschild radius R_{S} = 2GM/c^{2}. At a given radius, r, the matter moves, relative to the locally nonrotating observer, with Keplerian velocity, which we express in the units of the speed of light c as (Kato et al. 2008):
where
The disk is considered to be optically thick. The energy flux from a surface element is described by the effective temperature, which has a radial profile that takes the form:
where
combines the BH mass and the accretion rate, Ṁ, while f(r, a) is the factor accounting for the relativistic and boundary condition corrections (Page & Thorne 1974; Eq. (3.191) in Kato et al. 2008). In the Newtonian case, assuming the inner disk edge at r = 3, the correction factor is:
The disk extends from the ISCO up to an outer edge at r_{out} = 3000, which was chosen to reproduce the correct spectra down to the photon energy E ∼ 0.01kT_{*}. The radius of the ISCO depends on the spin as:
where the minus sign corresponds to a corotation of the disk and the BH and the plus sign is for the case of the retrograde rotation and:
The local spectrum of the disk is assumed to be a diluted blackbody with the color correction f_{col} = 1.7 (Shimura & Takahara 1995) with the flux being:
where E′ is photon energy measured in the frame comoving with the matter of the disk (where all quantities are denoted with primes) and B_{E′} is the Planck function. The angular distribution of the specific intensity and the polarization is assumed to correspond to the case of electronscatteringdominated semiinfinite atmosphere (Chandrasekhar & Breen 1947; Chandrasekhar 1960; Sobolev 1949, 1963). The threecomponent Stokes vector fully describes a linearly polarized radiation field as a function of photon energy E′ and zenith angle ζ′ (measured in the comoving frame):
where the angular distribution can be approximated as (Pihajoki et al. 2018)
and the PD (Viironen & Poutanen 2004)
The Stokes parameter U is zero in the comoving frame because of the azimuthal symmetry. Furthermore, in our case, the escaping radiation is polarized perpendicular to the meridional plane formed by the normal to the local surface and the photon momentum, resulting in the local PA of χ_{0} = π/2 (Stokes Q < 0). We omitted the fourth component of the Stokes vector, describing the circular polarization, from our treatment here. However, it can be easily added, because the circular polarization degree is conserved along a light trajectory in a vacuum.
Once the model for the disk structure (i.e., radial dependence of the velocity profile) and the energy and angular dependence of the local Stokes vector are established, we need to consider how the polarized radiation is modified toward the observer in the curved spacetime. One way to account for that is to compute a set of geodesics from the vicinity of the BH to the observer (we focus on this in more detail in Sect. 2.3). The other possibility is to use a faster approach that exploits analytical approximations, which we describe in Sect. 2.2.
2.2. Analytical raytracing in the Schwarzschild metric
We previously tested the analytical method of calculating the rotation of the PA for the case of the planar light trajectories (Loktev et al. 2022). This assumption is fulfilled in the Schwarzschild metric but does not hold for the Kerr metric. Nevertheless, we can use this approach as the firstorder approximation for the case of a spinning BH.
The considered geometry is shown in Fig. 1. The disk is assumed to be flat and located in the equatorial plane of the spinning BH. An element of the accretion disk surface, located at the tip of the radius vector, r, is described in BoyerLindquist coordinates by an azimuth angle, φ, and the length, r. We utilized a threedimensional Cartesian coordinate system to describe the polarization frame’s rotation. The normal to the disk, , is aligned with the z axis and the direction to the observer, , is in the xz plane making an inclination, i, to the normal:
Fig. 1. Geometry of a flat accretion disk ring. An element of the flat equatorial disk, defined by the unit radiusvector, at azimuthal angle, φ, moves along the direction, . The coordinate system is based on the normal of the disk, and the observer vector, . A photon is emitted from the disk element along . The important angles used in this work are also shown. 
In these coordinates, the unit vector of the disk element,
makes an angle ψ with the observer direction:
Close to the disk surface, photon escapes along a unit vector,
where α is the angle between the radius vector and the observer direction,
and is related to ψ through the light bending formula (see below). The photon momentum makes an angle ζ with the normal
Here we assume that the fluid velocity is aligned with the azimuthal unit vector
The matter moves with Keplerian speed given by Eq. (1).
The unit vector of photon momentum in the laboratory frame and the one in the comoving frame are related by the Lorentz transformation:
where is the Lorentz factor and
is the Doppler factor and
We then also get
The geodesics are not computed explicitly; that is, we only need the relation between the angles α and ψ, for which we can use the following analytical approximate formula (Poutanen 2020a):
where y = 1 − cosψ and u = 1/r. The formula is not applicable for r < 1 (within the Schwarzschild radius); hence, for cases with r_{ISCO} < 1, namely, a ≳ 0.943, ARTPOL cannot be used.
Under the Schwarzschild metric assumption, the geodesics are flat; therefore, the parallel transport of the Lorentz frame along the geodesics is unnecessary. The rotation of the polarization basis is a sum of several simple rotations due to the gravitational light bending (the general relativity effect) and Lorentz aberration (the special relativity effect):
Analytical expressions for those have been derived in Loktev et al. (2022):
where and
Expressions for χ^{GR} and χ^{SR} in vector form can be found in Loktev et al. (2022).
The PA is measured in the main polarization basis, formed by the disk normal and the observer vector:
The total PA in this basis is given by the sum of the relativistic rotations and the intrinsic PA χ_{0},
The PA χ_{0} in the comoving frame is described in the polarization basis that is formed by the local normal vector, , and the photon momentum vector in that frame, :
In our case of an electronscattering atmosphere, the intrinsic PA is χ_{0} = π/2.
The area of an element of the disk surface is expressed as
It occupies the solid angle dΩ on the sky and can be computed in Schwarzschild approximation as
where D is the distance to the source and the lensing factor ℒ is defined for Schwarzschild spacetime as (Beloborodov 2002)
It can be computed analytically (Poutanen 2020a) following Eq. (28):
Next, the intensity in the comoving frame is related to the observed one as follows:
where g is the redshift factor, which can be expressed in the case of a disk in the equatorial plane of the BH as (see Eq. (C.13) in Li et al. 2005):
where
In the case of a = 0, the redshift factor is reduced to .
The total disk flux is obtained by integration over the surface:
where M is the rotation matrix defined at every surface element dr dφ of the disk:
The PD and PA are then defined from the Stokes parameters of the total flux. The PD is obtained as
while the observed PA can be computed from either
or
Following Loktev et al. (2022), we can use dimensionless energy, x = E/kT_{*} (and x′=x/g), and scale the luminosity to to get the dimensionless luminosity, l_{x}, in the following form:
where t(r, a) = r^{3/4}f^{−1/4}(r, a) is related to the disk temperature (Eq. (5)). In this notation, the spectral shape and normalization of xl_{x} are independent of the BH mass and accretion rate.
Image of the disk can be reproduced using Cartesian coordinates X and Y (expressed in units of R_{S}) on the plane of the sky:
where Φ is the position angle of the point where the photon hits the plane of the sky measured counterclockwise from the projection of the disk axis on the sky and
is the approximation for the impact parameter, which is exact in the Schwarzschild case (Pechenick et al. 1983; Beloborodov 2002).
2.3. Numerical raytracing in Kerr metric
We went on to compare our results to numerical raytracing calculations to determine the accuracy of ARTPOL. We constructed the image of the BH accretion disk in the Kerr metric using the ARCMANCER code (Pihajoki et al. 2018), designed to integrate the exact geodesic equation numerically. The resulting trajectories were parametrized by the point at which they intersect the (observer) image plane.
To obtain an image of the accretion disk, we first defined the region of interest on the image plane to include the relevant region of the disk. The image plane was positioned at a distance r_{img} = 2500, where the effects of general relativity on the photon trajectory are negligible. We then introduced a grid of nodes and calculate the observed radiation intensity from the accretion disk at each node. The rectangular grid on the image plane is orthogonal to the vector pointing toward the observer. At each point on the grid, we define a local Lorentz frame with the basis vectors along BoyerLindquist coordinate grid. The Lorentz frames constitute the initial conditions for the geodesic curves with the tangent vectors normal to the image plane. Essentially, a set of parallel rays are propagated from the points of the grid towards the BH, backward in time, until they either intersect the disk surface or reach a distance of > r_{img} (as measured from the BH). When all the geodesics were calculated, we computed the local intensity and polarization for each point that intersects the disk surface, which we assume lies in the equatorial plane of the Kerr BH. We considered the velocity and energy dissipation profile following Novikov & Thorne (1973).
The Stokes parameters are paralleltransported in the “polarization frame” along the geodesic in Kerr spacetime. The polarization frame 𝒫 = {v^{a}, h^{a}} is a pair of orthogonal spacelike vectors that are orthogonal to the geodesic itself and the fourvelocity of the observer. Initially, at the image plane, the polarization frame consists of the corresponding two vectors of the local Lorentz frame. Then, the polarization frame is paralleltransferred along the geodesic to the disk surface to obtain the polarization parameters. Then, the frame is projected to the rest frame of the fluid in the disk (comoving frame). The propagation of the geodesics, transporting of the Lorentz frames, and the projection of the screen to the comoving frame are automatically performed by ARCMANCER.
All the necessary values, namely, the redshift g, the rotation of the polarization basis, χ^{tot}, and the emission angle, ζ′, can be defined in the comoving frame, where the geodesic hits the disk surface. The gravitational redshift factor g is computed using the curve tangents of the initial and the endpoints of the light trajectory. For more details on this numerical procedure, we refer to Sect. 4 of Pihajoki et al. (2018). We also corrected the redshift by a factor of to mitigate the redshift between the observer at r_{img} and the one at infinity. The total redshift agrees with the one given by the analytical expression given by Eq. (40).
To compute the overall flux from the disk in the case of numerical raytracing, we summed the intensities over the whole image grid as
where dx dy is the size of the pixel at r_{img}, measured in Schwarzschild radii, and M(x, y) is the rotation (Mueller) matrix along the trajectory from the point (x, y) on the image plane to the point on the disk surface with a coordinate r < r_{out}.
The secondary and higher order images of the disk and its bottom side are visible through the gap between the r_{ISCO} and the BH event horizon. From there, a segment of the disk surface can be seen multiple times, depending on how many revolutions its photons make around the BH before escaping toward the observer. We do not account for these trajectories in calculations with ARTPOL because we find the contribution of the secondary images to the observed flux are below 0.5% (see Sect. 3.3 below) for any inclination and spin values (see also Zhou et al. 2020). Lastly, we note that it is also possible to compute polarization rotation using Walker–Penrose theorem (e.g., see the appendix in Li et al. 2009); however, it does not significantly decrease the computing time as long as the geodesics are still computed numerically.
3. Applications
In this section, we show the applications of the ARTPOL technique to spectra and polarization signatures of matter near a Kerr BH. We first verify the accuracy of the ARTPOL for the PA rotation χ and the emission zenith angle ζ′ (measured in the comoving frame) for the case of the narrow disk ring located at different distances from the BH. We proceed to the comparison of polarized images of the accretion disk. Finally, we show the spectropolarimetric energy distributions of the accretion disk obtained with analytical and numerical approaches. In all cases, we assume that the matter (ring or disk) rotates counterclockwise. Positive spin values correspond to the prograde rotation (BH spin vector aligned with the orbital momentum of the accretion disk).
3.1. Narrow ring
The difference between Stokes parameters computed via ARTPOL and numerical ray tracing is expected to be high at the r_{ISCO}, as the framedragging effects – omitted in ARTPOL – are most important here. In Fig. 2, we show the total rotation angle χ^{tot} and the emission zenith angle ζ′ at the ISCO computed by two methods for different spin values and observer inclinations. For ARTPOL, both PA rotation and ζ′ are systematically shifted, over the azimuth, with respect to the numerical values. The difference is highest at φ ∼ 180° and it is generally larger for higher inclinations. The effect is caused by the close BH approach of the photon trajectories starting from regions of the disk located behind the BH. Naturally, the difference is larger for higher spin values. For example, in the case of a = 0.8 (r_{ISCO} ≈ 1.45) and i = 30° the maximum Δχ^{tot} is just 10° and Δζ′ is only 9°. We also note here that the largest error on χ^{tot} corresponds to a small emission angle cosζ′∼1 at which the PD is minimal (see Eq. (14)); therefore, the polarized flux is affected very slightly.
Fig. 2. Influence of relativistic effects on the PA and the emission angle. The rotation of the PA (a sum of the GR and SR effects) χ^{tot} (upper panels) and the emission angle ζ′ (lower panels) computed at r = r_{ISCO} (given in the upper left corner of each panel) for spins a = 0.2 (left panels), 0.5 (middle panels), and 0.8 (right panels) at three inclinations i = 30° (red), 60° (green), and 80° (blue lines). The results using ARTPOL code are shown with the solid lines, while the results from the raytracing in Kerr metric are presented with the dashed lines. The difference between the two values is shown under each panel. 
While the accuracy of ARTPOL can be low for certain parts of the narrow ring of matter at the ISCO, for the case of the zero torque boundary condition (Novikov & Thorne 1973), the total flux emerging from this ring is zero. The error is then weighted with the dissipation profile, which achieves maximum at distances further than ISCO. In Fig. 3, we show χ^{tot} and ζ′ for a = 0.8 at different radii r = 2, 3, and 5, computed with ARTPOL and numerical raytracing. The shapes of the lines are similar, but the deviations between the two methods are much smaller. For example, the maximum Δχ^{tot} and Δζ′ for r = 2 and i = 30° are 5° and 4°, respectively, and those are the peak values only at φ ≈ 180°. We note that for r > 5, ARTPOL gives χ^{tot} and ζ′ that are nearly indistinguishable from the numerical results, even for BHs with extreme spins.
Fig. 3. Same as Fig. 2 for the spin of a = 0.8 and three ring radii of r = 5 (left panels), 3 (middle panels), and 2 (right panels). 
3.2. Disk imaging
ARTPOL can be used to construct images of the disk (see more details in Loktev et al. 2022). In Fig. 4, we show the observed intensity and polarization of an accretion disk for BH spin parameters a = 0.2, 0.5, and 0.8, and observer inclinations of i = 30°, 60°, and 80°. The observed intensity for the innermost parts of the standard disk with the electron scattering dominating atmosphere is colormapped. The sticks depict the polarization vector computed with the raytracing and analytical techniques in the Schwarzschild metric. The difference in polarization parameters can barely be seen only in the innermost regions for high spin cases. Otherwise, the effect of BH rotation is so tiny that the sticks visually coincide. The sticks are displaced predominantly in the azimuthal direction, mainly caused by the framedragging effect, while the radial displacement due to the approximate impact parameter estimation is negligible. The most distorted region of the disk appears to be directly behind the black hole, where the surface of the disk appears to be dragged with the BH spin. This effect is, again, especially pronounced for higher observer inclinations.
Fig. 4. Polarization parameters of radiation from the disk at inclinations i = 30° (upper row), 60° (middle row), and 80° (lower row), and different spin parameters a = 0.2 (left column), 0.5 (middle column), and 0.8 (right column). The black lines outline an even polar grid on the disk with rays spaced by 30° in azimuth and contours spaced by 5R_{S} in radius (the corresponding values are denoted in each ring). The sticks represent the polarization parameters derived for the center of each grid segment. The polarized disk emission were calculated using the exact numerical interaction of the geodesics in the Kerr metric (with the ARCMANCER code; green lines) and using the approximate analytical formulae in the Schwarzschild metric (with the ARTPOL code; cyan lines). The coordinates of the cyan sticks were computed using Eq. (50). The length of the sticks is proportional to the observed PD from the region and the orientation shows the observed PA. The 5% polarization stick is shown in the bottom right panel to scale. The colors reflect the logarithm of the bolometric intensity relative to the maximal value across all panels. We note that differences between the exact numerical solution and the approximate ARTPOL method are visible only at the inner radii for systems with the highest spin values. 
3.3. Disk spectra and polarization
Using the ARTPOL technique, we perform fast computations of the spectral properties, which are essential to extract the spin information in both reflection spectroscopy and continuumfitting methods and spectral dependence of the polarization signatures. Comparisons of the computed spectral energy distribution, PD, and PA for the case of the optically thick, razorthin accretion disk to the exact raytracing calculations are shown in Fig. 5. The local model for the disk is described in Sect. 2.1.
Fig. 5. Spectra of the luminosity xl_{x} (upper panels), PD (middle panels) and PA (bottom panels) for different inclinations i = 30° (left), 60° (middle), and 80° (right column) and BH spin a = −1 (pink lines), 0 (red), 0.5 (orange), 0.8 (green), and 0.94 (blue). The solid lines correspond to the results obtained with the approximate analytical formulae (ARTPOL), while the dashed lines correspond to the exact numerical integration of the geodesics (ARCMANCER). The accretion disk extends from r_{ISCO} to r_{out} = 3000 for each case. Narrow panels show the difference between the results given by the two methods. The black dotted lines denote the results for a Newtonian disk from Loktev et al. (2022). 
In the upper panels of Fig. 5, we show the multicolor blackbody accretion disk spectra viewed at inclinations i = 30° (left column), 60° (middle column), and 80° (right column). We consider different BH spins from a retrograde one a = −1 to the prograde a = 0.94. The solid lines correspond to the calculations using ARTPOL, and the dashed lines correspond to numerical raytracing calculations. The differences between the methods are shown below each panel. The deviations increase with the spin value, energy, and inclination, mainly due to the increased importance of the framedragging effects. These are more important for the innermost radii, contributing to higher energies and the photon trajectories passing close to the BH, which appear mostly at higher inclinations. For higher spin values, the ISCO radii are smaller, and the framedragging effect is more pronounced. At the same time, in the case of the retrograde BH rotation (BH rotates clockwise in our case), a = −1, smaller deviations are caused by the more distant location of the ISCO (r_{ISCO} = 4.5). We find that for all the considered spin values, the deviations in the flux values are smaller than 20% at energies where the flux is relatively high (i.e., where it drops by less than a factor of 10 from the peak value).
At this point, we may consider the question of whether the difference between the fluxes computed with the two methods has its origin in the presence of secondary images. Using ARCMANCER, we computed the accretion disk spectra for two extreme BH spins and two inclinations. The results shown in Fig. 6 demonstrate that the secondary images contribute at most 0.5% of the flux for the high spin BH and large inclination. At lower inclinations or smaller spin values, their contributions are even smaller.
Fig. 6. Contribution of secondary images to the observed accretion disk flux for two inclinations i = 30° (left panel) and 80° (right panel) and three BH spins a = −1 (pink line), 0 (red), and 0.95 (blue) as computed with the numerical raytracing code ARCMANCER. The contribution of the secondary images for i = 30° and a = 0.95 is below 10^{−5}. The solid lines correspond to the full luminosity and the dashed lines to the contributions of the secondary images. 
The middle and lower panels of Fig. 5 show the observed PD and PA. The trend of increasing deviations towards higher spins, energies, and inclinations is also evident in these quantities. For all considered cases, however, the deviations remain small as compared to the typical observational errors achieved in IXPE observations. The maximum deviations in the predicted PD are about 0.4% for a = 0.94 and are much smaller than the expected accuracy of measurements with IXPE for smaller spins. The error in PA is below 1° for a ≤ 0.5, and even for a = 0.94, it is at most 5° at energies where the flux is relatively high. This indicates that fitting the Stokes spectra performed with our analytical technique will work fairly well for all inclinations and all considered spins.
We note the difference between the values obtained for the Newtonian disk (Shakura 1973) in the Schwarzschild metric in Loktev et al. (2022) and the zerospin case considered here. The difference comes solely from the relativistic temperature profile in Eq. (5), which is different from the Newtonian temperature profile from Eq. (7) assumed in Loktev et al. (2022). Notably, the Newtonian disk polarization profiles more closely resemble relativistic profiles with the spin parameter of a ∼ 0.2 than a = 0.
The reduction in calculation time using ARTPOL compared to ARCMANCER is consistent with that reported in Loktev et al. (2022), as we are comparing the same two methods. Specifically, our calculations of one image using the approximate analytical formulae are performed in ∼0.01 s at a MacBook Pro laptop using single processor, while the calculations of one polarized spectrum take about 0.05 s. This has to be compared to calculating an image 400 × 400 pixels using a raytracing code ARCMANCER, which takes about 1000 s. Our method offers flexibility, allowing for its application to arbitrary disk geometries, different energy dissipation profiles, and various local spectra and polarization.
4. Summary
We present the analytical technique ARTPOL for computing the images and spectropolarimetric characteristics of relativistic accretion disks around Kerr BHs. While the velocity and energy dissipation profiles used in the calculations correspond to the spacetime around spinning BHs, the photon paths have been considered in the Schwarzschild metric. We showed that ARTPOL technique is highly efficient, allowing for reduction of the computing time by a factor of 10^{4}, and remains accurate, generally within 10% in flux, 0.2% in PD, and 1° in PA, for spin values a ≤ 0.5. This enables ARTPOL technique to be used for a broad range of BH spin parameters. The deviations from the exact results obtained by numerical raytracing techniques systematically increase towards the highest spins and inclinations, reducing the accuracy down to 20% in flux, 0.5% in PD, and 7° in PA for a = 0.8–0.94 and i = 80°. The systematic discrepancy arises from the framedragging effect, which is omitted in our analytical raytracing calculations.
Applications of the analytical technique include fast computations of static spectropolarimetric signatures of accretion disks of Xray binaries in the soft, intermediate, and very high states (Ratheesh et al. 2024; Rodriguez Cavero et al. 2023; Podgorný et al. 2023), while accounting for the relativistic effects on spectra and polarimetric signatures of Comptonization observed in the hard state of Xray binaries and Seyfert galaxies (Krawczynski et al. 2022; Marinucci et al. 2022; Gianolli et al. 2023; Ingram et al. 2023; Tagliacozzo et al. 2023). The longanticipated detection of the quasiperiodic oscillations from BH Xray binaries in the polarized light would open new prospects for ARTPOL in terms of fast calculations of the phaseresolved characteristics. The method is also well suited for the warped disks whose spectropolarimetric signatures have not been studied in detail so far.
Acknowledgments
We acknowledge support from the Jenny and Antti Wihuri foundation (VL), the Academy of Finland grants 347003, 355672 (AV), 333112, 349906 (JP), and 349144 (VL, AV, JP), and the German Academic Exchange Service (DAAD) project 57525212 (VFS). VFS thanks Deutsche Forschungsgemeinschaft (DFG) grant WE 1312/591 for financial support.
References
 Beloborodov, A. M. 2002, ApJ, 566, L85 [Google Scholar]
 CárdenasAvendaño, A., Lupsasca, A., & Zhu, H. 2023, Phys. Rev. D, 107, 043030 [Google Scholar]
 Chan, C.K., Psaltis, D., & Özel, F. 2013, ApJ, 777, 13 [Google Scholar]
 Chandrasekhar, S. 1960, Radiative Transfer (New York: Dover) [Google Scholar]
 Chandrasekhar, S., & Breen, F. H. 1947, ApJ, 105, 435 [NASA ADS] [CrossRef] [Google Scholar]
 Connors, P. A., & Stark, R. F. 1977, Nature, 269, 128 [Google Scholar]
 Connors, P. A., Piran, T., & Stark, R. F. 1980, ApJ, 235, 224 [Google Scholar]
 Dexter, J. 2016, MNRAS, 462, 115 [Google Scholar]
 Dexter, J., & Agol, E. 2009, ApJ, 696, 1616 [Google Scholar]
 Dovčiak, M., Muleri, F., Goosmann, R. W., Karas, V., & Matt, G. 2008, MNRAS, 391, 32 [Google Scholar]
 Draghis, P. A., Miller, J. M., Zoghbi, A., et al. 2023, ApJ, 946, 19 [Google Scholar]
 Gianolli, V. E., Kim, D. E., Bianchi, S., et al. 2023, MNRAS, 523, 4468 [Google Scholar]
 Ingram, A., Maccarone, T. J., Poutanen, J., & Krawczynski, H. 2015, ApJ, 807, 53 [Google Scholar]
 Ingram, A., Ewing, M., Marinucci, A., et al. 2023, MNRAS, 525, 5437 [Google Scholar]
 Kato, S., Fukue, J., & Mineshige, S. 2008, BlackHole Accretion Disks– Towards a New Paradigm (Kyoto: Kyoto Univ. Press) [Google Scholar]
 Kerr, R. P. 1963, Phys. Rev. Lett., 11, 237 [Google Scholar]
 Krawczynski, H. 2012, ApJ, 754, 133 [Google Scholar]
 Krawczynski, H., & Beheshtipour, B. 2022, ApJ, 934, 4 [Google Scholar]
 Krawczynski, H., Muleri, F., Dovčiak, M., et al. 2022, Science, 378, 650 [NASA ADS] [CrossRef] [Google Scholar]
 Li, L.X., Zimmerman, E. R., Narayan, R., & McClintock, J. E. 2005, ApJS, 157, 335 [Google Scholar]
 Li, L.X., Narayan, R., & McClintock, J. E. 2009, ApJ, 691, 847 [CrossRef] [Google Scholar]
 Lightman, A. P., & Shapiro, S. L. 1975, ApJ, 198, L73 [Google Scholar]
 Loktev, V., Salmi, T., Nättilä, J., & Poutanen, J. 2020, A&A, 643, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Loktev, V., Veledina, A., & Poutanen, J. 2022, A&A, 660, A25 [Google Scholar]
 Loskutov, V. M., & Sobolev, V. V. 1979, Astrofizika, 15, 241 [NASA ADS] [Google Scholar]
 Loskutov, V. M., & Sobolev, V. V. 1981, Astrofizika, 17, 97 [NASA ADS] [Google Scholar]
 Marinucci, A., Muleri, F., Dovciak, M., et al. 2022, MNRAS, 516, 5907 [Google Scholar]
 McClintock, J. E., Narayan, R., & Steiner, J. F. 2014, Space Sci. Rev., 183, 295 [Google Scholar]
 Miller, J. M., Reynolds, C. S., Fabian, A. C., Miniutti, G., & Gallo, L. C. 2009, ApJ, 697, 900 [Google Scholar]
 Motta, S. E., Belloni, T. M., Stella, L., MuñozDarias, T., & Fender, R. 2014, MNRAS, 437, 2554 [Google Scholar]
 Motta, S. E., Belloni, T., Stella, L., et al. 2022, MNRAS, 517, 1469 [Google Scholar]
 Narayan, R., Palumbo, D. C. M., Johnson, M. D., et al. 2021, ApJ, 912, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Novikov, I. D., & Thorne, K. S. 1973, Black Holes (Les Astres Occlus) (New York: Gordon& Breach), 343 [Google Scholar]
 Page, D. N., & Thorne, K. S. 1974, ApJ, 191, 499 [NASA ADS] [CrossRef] [Google Scholar]
 Pechenick, K. R., Ftaclas, C., & Cohen, J. M. 1983, ApJ, 274, 846 [Google Scholar]
 Pihajoki, P., Mannerkoski, M., Nättilä, J., & Johansson, P. H. 2018, ApJ, 863, 8 [Google Scholar]
 Pineault, S. 1977, MNRAS, 179, 691 [Google Scholar]
 Pineault, S., & Roeder, R. C. 1977a, ApJ, 212, 541 [Google Scholar]
 Pineault, S., & Roeder, R. C. 1977b, ApJ, 213, 548 [Google Scholar]
 Podgorný, J., Marra, L., Muleri, F., et al. 2023, MNRAS, 526, 5964 [Google Scholar]
 Poutanen, J. 2020a, A&A, 640, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Poutanen, J. 2020b, A&A, 641, A166 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ratheesh, A., Dovčiak, M., Krawczynski, H., et al. 2024, ApJ, 964, 77 [Google Scholar]
 Rees, M. J. 1975, MNRAS, 171, 457 [Google Scholar]
 Reynolds, C. S. 2021, ARA&A, 59, 117 [NASA ADS] [CrossRef] [Google Scholar]
 Rodriguez Cavero, N., Marra, L., Krawczynski, H., et al. 2023, ApJ, 958, L8 [Google Scholar]
 Shakura, N. I. 1973, Sov. Astron., 16, 756 [Google Scholar]
 Shimura, T., & Takahara, F. 1995, ApJ, 445, 780 [Google Scholar]
 Sobolev, V. V. 1949, Uch. Zap. Leningrad Univ, 16 [Google Scholar]
 Sobolev, V. V. 1963, A treatise on Radiative Transfer (Princeton: Van Nostrand) [Google Scholar]
 Stark, R. F., & Connors, P. A. 1977, Nature, 266, 429 [NASA ADS] [CrossRef] [Google Scholar]
 Tagliacozzo, D., Marinucci, A., Ursini, F., et al. 2023, MNRAS, 525, 4735 [Google Scholar]
 Uttley, P., Cackett, E. M., Fabian, A. C., Kara, E., & Wilkins, D. R. 2014, A&ARv, 22, 72 [NASA ADS] [CrossRef] [Google Scholar]
 Viironen, K., & Poutanen, J. 2004, A&A, 426, 985 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Walker, M., & Penrose, R. 1970, Commun. Math. Phys., 18, 265 [NASA ADS] [CrossRef] [Google Scholar]
 Weisskopf, M. C., Soffitta, P., Baldini, L., et al. 2022, JATIS, 8, 026002 [NASA ADS] [Google Scholar]
 Yang, X.L., & Wang, J.C. 2014, A&A, 561, A127 [Google Scholar]
 Zhang, W., Dovčiak, M., & Bursa, M. 2019, ApJ, 875, 148 [NASA ADS] [CrossRef] [Google Scholar]
 Zhou, M., Ayzenberg, D., Bambi, C., & Nampalliwar, S. 2020, Phys. Rev. D, 101, 043010 [Google Scholar]
All Figures
Fig. 1. Geometry of a flat accretion disk ring. An element of the flat equatorial disk, defined by the unit radiusvector, at azimuthal angle, φ, moves along the direction, . The coordinate system is based on the normal of the disk, and the observer vector, . A photon is emitted from the disk element along . The important angles used in this work are also shown. 

In the text 
Fig. 2. Influence of relativistic effects on the PA and the emission angle. The rotation of the PA (a sum of the GR and SR effects) χ^{tot} (upper panels) and the emission angle ζ′ (lower panels) computed at r = r_{ISCO} (given in the upper left corner of each panel) for spins a = 0.2 (left panels), 0.5 (middle panels), and 0.8 (right panels) at three inclinations i = 30° (red), 60° (green), and 80° (blue lines). The results using ARTPOL code are shown with the solid lines, while the results from the raytracing in Kerr metric are presented with the dashed lines. The difference between the two values is shown under each panel. 

In the text 
Fig. 3. Same as Fig. 2 for the spin of a = 0.8 and three ring radii of r = 5 (left panels), 3 (middle panels), and 2 (right panels). 

In the text 
Fig. 4. Polarization parameters of radiation from the disk at inclinations i = 30° (upper row), 60° (middle row), and 80° (lower row), and different spin parameters a = 0.2 (left column), 0.5 (middle column), and 0.8 (right column). The black lines outline an even polar grid on the disk with rays spaced by 30° in azimuth and contours spaced by 5R_{S} in radius (the corresponding values are denoted in each ring). The sticks represent the polarization parameters derived for the center of each grid segment. The polarized disk emission were calculated using the exact numerical interaction of the geodesics in the Kerr metric (with the ARCMANCER code; green lines) and using the approximate analytical formulae in the Schwarzschild metric (with the ARTPOL code; cyan lines). The coordinates of the cyan sticks were computed using Eq. (50). The length of the sticks is proportional to the observed PD from the region and the orientation shows the observed PA. The 5% polarization stick is shown in the bottom right panel to scale. The colors reflect the logarithm of the bolometric intensity relative to the maximal value across all panels. We note that differences between the exact numerical solution and the approximate ARTPOL method are visible only at the inner radii for systems with the highest spin values. 

In the text 
Fig. 5. Spectra of the luminosity xl_{x} (upper panels), PD (middle panels) and PA (bottom panels) for different inclinations i = 30° (left), 60° (middle), and 80° (right column) and BH spin a = −1 (pink lines), 0 (red), 0.5 (orange), 0.8 (green), and 0.94 (blue). The solid lines correspond to the results obtained with the approximate analytical formulae (ARTPOL), while the dashed lines correspond to the exact numerical integration of the geodesics (ARCMANCER). The accretion disk extends from r_{ISCO} to r_{out} = 3000 for each case. Narrow panels show the difference between the results given by the two methods. The black dotted lines denote the results for a Newtonian disk from Loktev et al. (2022). 

In the text 
Fig. 6. Contribution of secondary images to the observed accretion disk flux for two inclinations i = 30° (left panel) and 80° (right panel) and three BH spins a = −1 (pink line), 0 (red), and 0.95 (blue) as computed with the numerical raytracing code ARCMANCER. The contribution of the secondary images for i = 30° and a = 0.95 is below 10^{−5}. The solid lines correspond to the full luminosity and the dashed lines to the contributions of the secondary images. 

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.