artpol : Analytical ray-tracing method for spectro-polarimetric properties of accretion disks around Kerr black holes

Spectro-polarimetric signatures of accretion disks in X-ray binaries and active galactic nuclei contain information about the masses and spins of their central black holes, as well as the geometry of matter close to the compact objects. This information can be extracted using the means of X-ray polarimetry. In this work, we present a fast analytical ray-tracing technique for polarized light artpol that helps obtain the spinning black hole parameters from the observed properties. This technique can replace the otherwise time-consuming numerical ray-tracing calculations. We show that artpol proves accurate for Kerr black holes with dimensionless spin parameter a ≤ 0 . 94 while being over four orders of magnitude faster than direct ray-tracing calculations. This approach opens broad prospects for directly fitting the spectro-polarimetric data from the Imaging X-ray Polarimetry Explorer .


Introduction
The spin of a black hole (BH) is a fundamental parameter that controls the behavior of the inflowing matter, the accretion disk, and the properties of the outflowing material, the relativistic jets.The magnitude of spin determines the curvature of space-time close to the BH, determines the energy dissipation profile, and affects the spectral energy distribution of the observed emission.The spin values for both Galactic BHs in X-ray 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 reflection/reverberation method (Miller et al. 2009;McClintock et al. 2014;Uttley et al. 2014).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 X-ray binaries, an additional constraint on the parameters comes from the relativistic precession model (Motta et al. 2014(Motta et al. , 2022)), which links the mass, spin, and radius with the characteristic frequencies found in the X-ray 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.
Polarization of radiation escaping from the BH vicinity can be used as a fine tool to determine the curvature of space-time.In this context, X-ray polarimetric signatures of accretion disks have long been anticipated to carry important information about the BH spin, and the launch of Imaging X-ray Polarimetry Explorer (IXPE, Weisskopf et al. 2022) has opened these new exciting possibilities.At the same time, the observed polarimetric signatures of BH X-ray binaries have proved many previously pro-posed models to fail in describing the data (Krawczynski et al. 2022;Rodriguez Cavero et al. 2023;Ratheesh et al. 2023), requiring to alter basic assumptions on the geometry and radiative mechanisms producing local spectra, which are often impossible to alter in the data-fitting models.A fast tool that relates the local spectra with the observed spectro-polarimetric signatures is required.
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 electron-scattering, semiinfinite plane-parallel 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 Monte-Carlo (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 frame-dragging 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 Walker & Penrose (1970) constant of motion (Connors et al. 1980;Dovčiak et al. 2008;Ingram et al. 2015).The geodesics, in turn, are computed using the ray-tracing 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 ≲ 100R 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 ray-tracing technique is too computationally expensive for a direct data fitting.Instead, pre-computed geodesics have been used to accelerate calculations (e.g., Li et al. 2005;Krawczynski 2012;Krawczynski & Beheshtipour 2022).Alternatively, analytical and semi-analytical approaches can be applied to solve geodesic equations (Dexter & Agol 2009;Yang & Wang 2014;Cárdenas-Avendaño et al. 2023).Their applicability is, however, limited, e.g., the semi-analytical expressions for geodesics can only be used for the equatorial plane of the BH.Finally, in the Schwarzschild metric, the calculation of geodesics may be omitted, and the photon trajectory can be treated using an approximation to the light bending relation (Beloborodov 2002;Poutanen 2020a).The latter approach gave reliable results for the low-energy synchrotron emission observed from the supermassive BH M87*, with an almost face-on disk (Narayan et al. 2021).
Physical understanding of the modifications to polarization signatures caused by the curved space-time and fast motions of matter is difficult to achieve if we use the implicit Walker-Penrose constant.For this, one can use the explicit analytical expression for the rotation of the polarization plane along the photon path, which we call the analytical ray-tracing 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).The method was first proposed in Pineault (1977) and applied, in the context of accretion disks, to the PA rotation caused by, separately, general and special relativity, 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 spectro-polarimetric 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 ray-tracing 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.

Local model of the disk emission
We consider a BH with mass M and dimensionless spin a = Jc/GM 2 , where J is the angular momentum.For the sake of simplicity, in this paper, we consider a standard equatorial geometrically thin accretion disk (Novikov & Thorne 1973).The BH is situated at the origin with the spin directed along the z-axis and is orthogonal to the disk plane.We assume that the state of the fluid only depends on the Boyer-Lindquist 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 the radial profile in the form where combines the BH mass and the accretion rate Ṁ, and 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 electron-scatteringdominated semi-infinite atmosphere (Chandrasekhar & Breen 1947;Chandrasekhar 1960;Sobolev 1949Sobolev , 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 (Nättilä & Pihajoki 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 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 space-time.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.

Analytical ray-tracing in 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 first-order 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 Boyer-Lindquist coordinates by an azimuth angle φ and the length r.We utilize a 3-dimensional Cartesian coordinate system to describe the polarization frame's rotation.The normal to the disk n is aligned with the z axis and the direction to the observer ô is in the x-z plane making an inclination i to the normal: In these coordinates, the unit vector of the disk element makes an angle ψ with the 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 k0 and the one in the comoving frame k′ 0 are related by the Lorentz transformation: where γ = 1/ 1 − β 2 is the Lorentz factor and is the Doppler factor and We then also get The geodesics are not computed explicitly -we only need the relation between the angles α and ψ, for which we 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 the cases with r ISCO < 1, i.e., 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 ã = (1 − cos α cos ψ)/(cos α − cos ψ), 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 n and the photon momentum vector in that frame k′ 0 : In our case of electron scattering atmosphere, the intrinsic PA χ 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 L is defined for Schwarzschild space-time as (Beloborodov 2002) It can be computed analytically (Poutanen 2020a) following Eq.( 28).
Next, the intensity in the comoving frame I ′ E ′ is related to the observed one as 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. C13 in Li et al. 2005) where In the case of a = 0, the redshift factor is reduced to g The total disk flux is obtained by integration over the surface: where M is the rotation matrix defined at every surface element of the disk, dr dφ, as following 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 Following Loktev et al. (2022), we can use dimensionless energy x = E/kT * (and x ′ = x/g) and scale the luminosity to σ SB T 4 * R 2 S 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).

Numerical ray-tracing in Kerr metric
Let us now compare our results to numerical ray-tracing calculations to determine the accuracy of artpol.We construct 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 are parametrized by the point at which they intersect the (observer) image plane.
To obtain an image of the accretion disk, we first define the region of interest on the image plane to include the relevant region of the disk.The image plane is positioned at a distance r img = 2500, where the effects of general relativity on the photon trajectory are negligible.We then introduce 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 Boyer-Lindquist 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 > r img (as measured from the BH).When all the geodesics are calculated, we compute 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 consider the velocity and energy dissipation profile following Novikov & Thorne (1973).
The Stokes parameters are parallel transported in the "polarization frame" along the geodesic in Kerr space-time.The polarization frame P = {v a , h a } is a pair of orthogonal space-like 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 correct the redshift by the factor g 0 = 1 − 1/r img 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 ray-tracing, we sum the intensities over the whole image grid as where dx dy is the size of the pixel at r img , measured in Schwarzschild radii, 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.These trajectories are not accounted for in 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 computed numerically anyway.

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/disk) rotates counterclockwise.Positive spin values correspond to the prograde rotation (BH spin vector aligned with the orbital momentum of the accretion disk).

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 frame-dragging effects -omitted in artpol -are most important here.In Fig. 2, we show the total rotation angle χ tot and the 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) and 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 ray-tracing in Kerr metric are presented with the dashed lines.The difference between the two values is shown under each panel.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, in general, is 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 PD is minimal (see Eq. 14); therefore, the polarized flux is affected very little.
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 ray-tracing.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 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.

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 color-mapped.The sticks depict the polarization vector computed with the ray-tracing 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 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 azimuthal direction, mainly caused by the frame-dragging 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.

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 continuum-fitting methods and spectral dependence of the polarization signatures.Comparison of the computed spectral energy distribution, PD, and PA for the case of the optically thick, razor-thin accretion disk to the exact ray-tracing calculations are shown in Fig. 5.The local model for the disk is described in Sect.2.1.
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 ray-tracing 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 frame-dragging effects.These are more im-portant 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 frame-dragging 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).
We can ask here a question 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 contribution is even smaller.
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, 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 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.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.
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 zero-spin 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 resemble relativistic profiles with the spin parameter a ∼ 0.2 more 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, 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 ray-tracing code arcmancer, which takes about 1000 s.Moreover, our analytical method offers flexibility, allowing for its application to arbitrary 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).disk geometries, different energy dissipation profiles, and various local spectra and polarization.

Summary
We presented the analytical technique artpol for computing the images and spectro-polarimetric characteristics of relativistic accretion disks around Kerr BHs.While the velocity and energy dissipation profiles used in the calculations correspond to the space-time 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, in general 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 ray-tracing 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 frame-dragging effect, which is omitted in our analytical ray-tracing calculations.
Applications of the analytical technique include fast computations of static spectro-polarimetric signatures of accretion disks of X-ray binaries in the soft/intermediate/very high states (Ratheesh et al. 2023;Rodriguez Cavero et al. 2023;Rawat et al. 2023a,b;Jana & Chang 2023;Podgorny et al. 2023), and account for the relativistic effects on spectra and polarimetric signatures of Comptonization observed in the hard state of X-ray 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 long-anticipated detection of the quasi-periodic oscillations from BH X-ray binaries in the polarized light would open new prospects for artpol in terms of fast calculation of the phase-resolved characteristics.The method is also well suited for the warped disks whose spectro-polarimetric signatures have not been studied in detail so far.

Fig. 1 .
Fig.1.Geometry of a flat accretion disk ring.An element of the flat equatorial disk, defined by the unit radius-vector r at azimuthal angle φ, moves along the direction û.The coordinate system is based on the normal of the disk n and the observer vector ô.A photon is emitted from the disk element along k0 .The important angles used in this work are also shown.

Fig. 4 .
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 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.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.

Fig. 5 .
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 fromLoktev et al. (2022).

Fig. 6 .
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 ray-tracing 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.