© The Authors 2024

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 space-time 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 X-ray binaries and their supermassive counterparts have been probed using the distinct imprints in spectral and timing properties via the continuum-fitting method and the iron line reflection-reverberation 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 X-ray 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 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.

The 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, it has been anticipated for some time that X-ray polarimetric signatures of accretion disks can carry important information about the BH spin. The launch of Imaging X-ray Polarimetry Explorer (IXPE; Weisskopf et al. 2022) has opened up these exciting research possibilities. At the same time, the observed polarimetric signatures of BH X-ray 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 data-fitting models. It is therefore necessary to adopt a fast tool that relates the local spectra to the observed spectro-polarimetric 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 electron-scattering, semi-infinite 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 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 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 ≲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⁵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; for instance, the semi-analytical 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 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).

The physical understanding of the modifications to polarization signatures caused by the curvature of space-time and fast motions of matter is difficult to achieve if we use the implicit Walker-Penrose 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 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). 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 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.

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², 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 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². At a given radius, r, the matter moves, relative to the locally non-rotating observer, with Keplerian velocity, which we express in the units of the speed of light c as (Kato et al. 2008):

$$\begin{array}{c}\hfill \beta =\frac{\mathcal{F}}{\mathcal{B}\sqrt{\mathcal{D}}}\sqrt{\frac{1}{2r}},\end{array}$$(1)

where

$$\begin{array}{cc}\hfill \mathcal{B}& =1+\frac{a}{\sqrt{8r^3}},\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill \mathcal{D}& =1-\frac{1}{r}+\frac{a^2}{4r^2},\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill \mathcal{F}& =1-\frac{a}{\sqrt{2r^3}}+\frac{a^2}{4r^3}.\hfill \end{array}$$(4)

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:

$$\begin{array}{c}\hfill T_{\mathrm{eff}}^4(r)=\frac{3GM\dot{M}}{8\pi {\sigma }_{\mathrm{SB}}{R}^3}f(r,a)=T_{\ast }^4\frac{f(r,a)}{r^3},\end{array}$$(5)

where

$$\begin{array}{c}\hfill T_{\ast }^4=\frac{3GM\dot{M}}{8\pi {\sigma }_{\mathrm{SB}}{R}_{\mathrm{S}}^3},\end{array}$$(6)

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:

$$\begin{array}{c}\hfill f(r)=1-\sqrt{\frac{3}{r}}.\end{array}$$(7)

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:

$$\begin{array}{c}\hfill r_{\mathrm{ISCO}}=\frac{1}{2}(3+Z_2\pm \sqrt{(3-Z_1)(3+Z_1+2Z_2)}),\end{array}$$(8)

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:

$$\begin{array}{cc}\hfill Z_1& =1+\sqrt[3]{1-a^2}(\sqrt[3]{1+a}+\sqrt[3]{1-a}),\hfill \end{array}$$(9)

$$\begin{array}{cc}\hfill Z_2& =\sqrt{3a^2+Z_1^2}.\hfill \end{array}$$(10)

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:

$$\begin{array}{c}\hfill F_{E^{\prime }}=\frac{\pi }{f_{\mathrm{col}}^4}{B}_{E^{\prime }}(f_{\mathrm{col}}{T}_{\mathrm{eff}}),\end{array}$$(11)

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-scattering-dominated semi-infinite atmosphere (Chandrasekhar & Breen 1947; Chandrasekhar 1960; Sobolev 1949, 1963). The three-component Stokes vector fully describes a linearly polarized radiation field as a function of photon energy E′ and zenith angle ζ′ (measured in the comoving frame):

$$\begin{array}{c}\hfill {{\mathit{I}}^{\prime }}_{E^{\prime }}({\zeta }^{\prime })=\left[\begin{array}{c}{I^{\prime }}_{E^{\prime }}\\ {Q^{\prime }}_{E^{\prime }}\\ {U^{\prime }}_{E^{\prime }}\end{array}\right]=\frac{1}{\pi }{F}_{E^{\prime }}{a}_{\mathrm{es}}({\zeta }^{\prime })\left[\begin{array}{c}1\\ p_{\mathrm{es}}({\zeta }^{\prime })cos2{\chi }_0\\ p_{\mathrm{es}}({\zeta }^{\prime })sin2{\chi }_0\end{array}\right],\end{array}$$(12)

where the angular distribution can be approximated as (Pihajoki et al. 2018)

$$\begin{array}{c}\hfill a_{\mathrm{es}}({\zeta }^{\prime })=\frac{60}{143}(1+2.3cos{\zeta }^{\prime }-0.3{cos}^2{\zeta }^{\prime })\end{array}$$(13)

and the PD (Viironen & Poutanen 2004)

$$\begin{array}{c}\hfill p_{\mathrm{es}}({\zeta }^{\prime })=0.1171\frac{1-cos{\zeta }^{\prime }}{1+3.582cos{\zeta }^{\prime }}.\end{array}$$(14)

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 χ₀ = π/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 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.

2.2. Analytical ray-tracing 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 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 utilized a three-dimensional Cartesian coordinate system to describe the polarization frame’s rotation. The normal to the disk, $\widehat{\mathit{n}}$, is aligned with the z axis and the direction to the observer, $\widehat{\mathit{o}}$, is in the x-z 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 radius-vector, $\widehat{\mathit{r}\mathbf{,}}$ at azimuthal angle, φ, moves along the direction, $\widehat{\mathit{v}}$. The coordinate system is based on the normal of the disk, $\widehat{\mathit{n}},$ and the observer vector, $\widehat{\mathit{o}}$. A photon is emitted from the disk element along ${\widehat{\mathit{k}}}_0$. The important angles used in this work are also shown.

$$\begin{array}{cc}\hfill \widehat{\mathit{n}}& =(0,0,1),\hfill \end{array}$$(15)

$$\begin{array}{cc}\hfill \widehat{\mathit{o}}& =(sini,0,cosi).\hfill \end{array}$$(16)

In these coordinates, the unit vector of the disk element,

$$\begin{array}{c}\hfill \widehat{\mathit{r}}=(cos\phi ,sin\phi ,0),\end{array}$$(17)

makes an angle ψ with the observer direction:

$$\begin{array}{c}\hfill cos\psi =\widehat{\mathit{r}}·\widehat{\mathit{o}}=sinicos\varphi .\end{array}$$(18)

Close to the disk surface, photon escapes along a unit vector,

$$\begin{array}{c}\hfill {\widehat{\mathit{k}}}_0=[sin\alpha \widehat{\mathit{o}}+sin(\psi -\alpha )\widehat{\mathit{r}}]/sin\psi ,\end{array}$$(19)

where α is the angle between the radius vector and the observer direction,

$$\begin{array}{c}\hfill cos\alpha =\widehat{\mathit{r}}·{\widehat{\mathit{k}}}_0,\end{array}$$(20)

and is related to ψ through the light bending formula (see below). The photon momentum makes an angle ζ with the normal

$$\begin{array}{c}\hfill cos\zeta =\widehat{\mathit{n}}·{\widehat{\mathit{k}}}_0=\frac{sin\alpha }{sin\psi }cosi.\end{array}$$(21)

Here we assume that the fluid velocity is aligned with the azimuthal unit vector

$$\begin{array}{c}\hfill \widehat{\mathit{\phi }}=(-sin\phi ,cos\phi ,0)=\widehat{\mathit{v}}.\end{array}$$(22)

The matter moves with Keplerian speed given by Eq. (1).

The unit vector of photon momentum in the laboratory frame ${\widehat{\mathit{k}}}_0$ and the one in the comoving frame ${\widehat{\mathit{k}}}_0^{\prime }$ are related by the Lorentz transformation:

$$\begin{array}{c}\hfill {\widehat{\mathit{k}}}_0^{\prime }=\delta [{\widehat{\mathit{k}}}_0-\gamma \beta \widehat{\mathit{v}}+(\gamma -1)\widehat{\mathit{v}}(\widehat{\mathit{v}}·{\widehat{\mathit{k}}}_0)],\end{array}$$(23)

where $\gamma =1/\sqrt{1-{\beta }^2}$ is the Lorentz factor and

$$\begin{array}{c}\hfill \delta =\frac{1}{\gamma (1-\beta cos\xi )}=\gamma (1+\beta cos{\xi }^{\prime })\end{array}$$(24)

is the Doppler factor and

$$\begin{array}{cc}\hfill cos\xi & =\widehat{\mathit{v}}·{\widehat{\mathit{k}}}_0=-\frac{sin\alpha }{sin\psi }sinisin\varphi ,\hfill \\ \hfill cos{\xi }^{\prime }& =\widehat{\mathit{v}}·{\widehat{\mathit{k}}}_0^{\prime }=\frac{cos\xi -\beta }{1-\beta cos\xi }.\hfill \end{array}$$(25)

We then also get

$$\begin{array}{c}\hfill cos{\zeta }^{\prime }=\widehat{\mathit{n}}·{\widehat{\mathit{k}}}_0^{\prime }=\delta cos\zeta .\end{array}$$(26)

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):

$$\begin{array}{c}\hfill cos\alpha =1-(1-u)y\{1+\frac{u^2{y}^2}{112}-\frac{\mathrm{e}}{100}uy[ln(1-\frac{y}{2})+\frac{y}{2}]\},\end{array}$$(27)

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):

$$\begin{array}{c}\hfill {\chi }^{\text{tot}}={\chi }^{\text{GR}}+{\chi }^{\text{SR}}.\end{array}$$(28)

Analytical expressions for those have been derived in Loktev et al. (2022):

$$\begin{array}{c}\hfill tan{\chi }^{\text{GR}}=\frac{cosisin\phi }{cos\phi +\stackrel{\sim }{a}sini},\end{array}$$(29)

where $\stackrel{\sim }{a}=(1-cos\alpha cos\psi )/(cos\alpha -cos\psi )$ and

$$\begin{array}{c}\hfill tan{\chi }^{\text{SR}}=-\beta \frac{cos\alpha cos\zeta }{{sin}^2\zeta -\beta cos\xi }.\end{array}$$(30)

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:

$$\begin{array}{c}\hfill {\widehat{\mathit{e}}}_1=\frac{\widehat{\mathit{n}}-cosi\widehat{\mathit{o}}}{sini},{\widehat{\mathit{e}}}_2=\frac{\widehat{\mathit{o}}\times \widehat{\mathit{n}}}{sini}.\end{array}$$(31)

The total PA in this basis is given by the sum of the relativistic rotations and the intrinsic PA χ₀,

$$\begin{array}{c}\hfill \chi ={\chi }_0+{\chi }^{\text{tot}}={\chi }_0+{\chi }^{\text{SR}}+{\chi }^{\text{GR}}.\end{array}$$(32)

The PA χ₀ in the comoving frame is described in the polarization basis that is formed by the local normal vector, $\widehat{\mathit{n}}$, and the photon momentum vector in that frame, ${\widehat{\mathit{k}}}_0^{\prime }$:

$$\begin{array}{c}\hfill {\widehat{\mathit{e}}}_1^{\prime }=\frac{\widehat{\mathit{n}}-cos{\zeta }^{\prime }{\widehat{\mathit{k}}}_0^{\prime }}{sin{\zeta }^{\prime }},{\widehat{\mathit{e}}}_2^{\prime }=\frac{{\widehat{\mathit{k}}}_0^{\prime }\times \widehat{\mathit{n}}}{sin{\zeta }^{\prime }}.\end{array}$$(33)

In our case of an electron-scattering atmosphere, the intrinsic PA is χ₀ = π/2.

The area of an element of the disk surface is expressed as

$$\begin{array}{c}\hfill \mathrm{d}S=\frac{R_{\mathrm{S}}^{2}r\mathrm{d}r\mathrm{d}\phi }{\sqrt{\mathcal{D}}}.\end{array}$$(34)

It occupies the solid angle dΩ on the sky and can be computed in Schwarzschild approximation as

$$\begin{array}{c}\hfill \mathrm{d}\mathrm{\Omega }=\frac{\mathrm{d}Scos\zeta }{D^2}\mathcal{L}=\frac{R_{\mathrm{S}}^2}{D^2}\frac{r\mathrm{d}r\mathrm{d}\phi }{\sqrt{\mathcal{D}}}\mathcal{L}cos\zeta ,\end{array}$$(35)

where D is the distance to the source and the lensing factor ℒ is defined for Schwarzschild space-time as (Beloborodov 2002)

$$\begin{array}{c}\hfill \mathcal{L}=\frac{1}{1-u}\frac{\mathrm{d}cos\alpha }{\mathrm{d}cos\psi }.\end{array}$$(36)

It can be computed analytically (Poutanen 2020a) following Eq. (28):

$$\begin{array}{c}\hfill \mathcal{L}=1+\frac{3u^2{y}^2}{112}-\frac{\mathrm{e}}{100}uy[2ln(1-\frac{y}{2})+y\frac{1-3y/4}{1-y/2}].\end{array}$$(37)

Next, the intensity in the comoving frame $I^{\prime }{}_{E^{\prime }}$ is related to the observed one as follows:

$$\begin{array}{c}\hfill {\mathit{I}}_E=g^3{{\mathit{I}}^{\mathbf{\prime }}}_{E^{\prime }},\end{array}$$(38)

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):

$$\begin{array}{c}\hfill g=E/E^{\prime }=\gamma [\mathcal{X}+\mathcal{Y}\beta +(\mathcal{X}\beta +\mathcal{Y})cos{\xi }^{\prime }],\end{array}$$(39)

where

$$\begin{array}{cc}\hfill \mathcal{X}& =\sqrt{\mathcal{D}/\mathcal{A}},\hfill \end{array}$$(40)

$$\begin{array}{cc}\hfill \mathcal{Y}& =a/\sqrt{4r^4\mathcal{A}},\hfill \end{array}$$(41)

$$\begin{array}{cc}\hfill \mathcal{A}& =1+(r+1)\frac{a^2}{4r^3}.\hfill \end{array}$$(42)

In the case of a = 0, the redshift factor is reduced to $g=\delta \sqrt{1-1/r}$.

The total disk flux is obtained by integration over the surface:

$$\begin{array}{cc}\hfill {\mathit{F}}_E& \equiv \left[\begin{array}{c}{F}_I(E)\\ F_Q(E)\\ F_U(E)\end{array}\right]={\displaystyle \int \mathrm{d}}\mathrm{\Omega }\left[\begin{array}{c}{I}_E\\ Q_E\\ U_E\end{array}\right]\hfill \\ \hfill & =\frac{R_{\mathrm{S}}^2}{D^2}{\displaystyle \underset{r_{\mathrm{ISCO}}}{\overset{r_{\mathrm{out}}}{\int }}}\frac{r\mathrm{d}r}{\sqrt{\mathcal{D}}}{\displaystyle {\int }_0^{2\pi }}\mathrm{d}\phi g^3\mathcal{L}cos\zeta \mathbf{M}(r,\phi ){{\mathit{I}}^{\mathbf{\prime }}}_{E^{\prime }}({\zeta }^{\prime }),\hfill \end{array}$$(43)

where M is the rotation matrix defined at every surface element dr dφ of the disk:

$$\begin{array}{c}\hfill \mathbf{M}(r,\phi )=\left[\begin{array}{ccc}1& 0& 0\\ 0& cos2{\chi }^{\mathrm{tot}}& -sin2{\chi }^{\mathrm{tot}}\\ 0& sin2{\chi }^{\mathrm{tot}}& cos2{\chi }^{\mathrm{tot}}\end{array}\right].\end{array}$$(44)

The PD and PA are then defined from the Stokes parameters of the total flux. The PD is obtained as

$$\begin{array}{c}\hfill p(E)=\frac{\sqrt{F_Q^2(E)+F_U^2(E)}}{F_I(E)},\end{array}$$(45)

while the observed PA can be computed from either

$$\begin{array}{c}\hfill tan2\chi (E)=\frac{F_U(E)}{F_Q(E)}\end{array}$$(46)

or

$$\begin{array}{c}\hfill tan\chi (E)=\frac{p(E)F_I(E)-F_Q(E)}{F_U(E)}.\end{array}$$(47)

Following Loktev et al. (2022), we can use dimensionless energy, x = E/kT_* (and x′=x/g), and scale the luminosity to ${\sigma }_{\text{SB}}{T}_{*}^4{R}_{\text{S}}^2$ to get the dimensionless luminosity, l_x, in the following form:

$$\begin{array}{cc}\hfill x{\mathit{l}}_x& =\frac{60}{{\pi }^4}{\displaystyle \underset{r_{\mathrm{ISCO}}}{\overset{r_{\mathrm{out}}}{\int }}}\frac{r\mathrm{d}r}{\sqrt{\mathcal{D}}}{\displaystyle {\int }_0^{2\pi }}\mathrm{d}\phi \mathcal{L}cos\zeta \frac{g^4}{f_{\mathrm{c}}^4}\hfill \\ \hfill & \times \frac{x^{\prime 4}{a}_{\mathrm{es}}({\zeta }^{\prime })}{{\mathrm{e}}^{{\displaystyle x^{\prime }t(r,a)/f_{\mathrm{c}}}}-1}\left[\begin{array}{c}1\\ p_{\mathrm{es}}({\zeta }^{\prime })cos(2\chi )\\ p_{\mathrm{es}}({\zeta }^{\prime })sin(2\chi )\end{array}\right],\hfill \end{array}$$(48)

where t(r, a) = r^3/4f^−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:

$$\begin{array}{c}\hfill \left[\begin{array}{c}X\\ Y\end{array}\right]=b\left[\begin{array}{c}-sin\mathrm{\Phi }\\ cos\mathrm{\Phi }\end{array}\right]=\frac{r}{\sqrt{\mathcal{D}}}\frac{sin\alpha }{sin\psi }\left[\begin{array}{c}sin\phi \\ -cosicos\phi \end{array}\right],\end{array}$$(49)

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

$$\begin{array}{c}\hfill b=\frac{r}{\sqrt{\mathcal{D}}}sin\alpha \end{array}$$(50)

is the approximation for the impact parameter, which is exact in the Schwarzschild case (Pechenick et al. 1983; Beloborodov 2002).

2.3. Numerical ray-tracing in Kerr metric

We went on to compare our results to numerical ray-tracing 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 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 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 parallel-transported in the “polarization frame” along the geodesic in Kerr space-time. The polarization frame 𝒫 = {v^a, h^a} is a pair of orthogonal space-like vectors that are orthogonal to the geodesic itself and the four-velocity 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 parallel-transferred 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 $g_0=\sqrt{1-1/r_{\mathrm{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 summed the intensities over the whole image grid as

$$\begin{array}{c}\hfill {\mathit{F}}_E=\frac{R_{\mathrm{S}}^2}{D^2}{\displaystyle \int \underset{r<r_{\mathrm{out}}}{\int }}\mathrm{d}x\mathrm{d}y\mathbf{M}(x,y){{\mathit{I}}^{\prime }}_{E^{\prime }}(r,{\zeta }^{\prime }),\end{array}$$(51)

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 spectro-polarimetric 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 frame-dragging 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 = rISCO (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 ray-tracing 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 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 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 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 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.

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 5RS 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 continuum-fitting 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, 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.

Fig. 5. Spectra of the luminosity xlx (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 rISCO to rout = 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 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 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 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).

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

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 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 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 ray-tracing 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 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⁴, 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 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, 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 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 calculations 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.

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/59-1 for financial support.

References

All Figures

Fig. 1. Geometry of a flat accretion disk ring. An element of the flat equatorial disk, defined by the unit radius-vector, $\widehat{\mathit{r}\mathbf{,}}$ at azimuthal angle, φ, moves along the direction, $\widehat{\mathit{v}}$. The coordinate system is based on the normal of the disk, $\widehat{\mathit{n}},$ and the observer vector, $\widehat{\mathit{o}}$. A photon is emitted from the disk element along ${\widehat{\mathit{k}}}_0$. 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 = rISCO (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 ray-tracing 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 5RS 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 xlx (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 rISCO to rout = 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 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.

In the text