© ESO 2022

1. Introduction

Accretion disks around compact objects – neutron stars (NSs) and black holes (BHs) – are among the most efficient energy conversion engines in the Universe. As the matter spirals in towards the compact object, it releases the excess of gravitational energy via radiation. Understanding of the physical mechanisms leading to energy liberation and geometrical properties of the emitting medium is an important goal of the modern high-energy astrophysics. Recent investigations have been focused on the innermost parts of the accretion disk in the regime of strong gravity. X-ray spectroscopy (Reynolds 2014; Bambi et al. 2021) and timing techniques (Revnivtsev et al. 1999; Gilfanov et al. 2003; Uttley et al. 2014; Axelsson & Veledina 2021) have been exploited to trace the geometry of the inner parts of the accretion flow and to identify contributions of various components that are emitting in the X-ray energy range.

Polarimetry is known to be a fine measure of the geometry and radiative processes operating in the accretion flow. Numerous efforts are being aimed at predicting and finding distinct signatures of the accretion disk in a strong gravity regime from the polarimetric information (e.g., Dovčiak et al. 2008; Li et al. 2009; Ingram et al. 2015). Polarimetry in the infrared and millimeter bands was recently shown to be a powerful tool for studying the structure of accretion flows in the vicinity of supermassive BHs in the Milky Way (GRAVITY Collaboration 2018; Bower et al. 2018) and in M 87 (Event Horizon Telescope Collaboration 2021a,b). For studies of NSs and stellar-mass BHs, radiating mostly in the X-ray range, the full capacity of polarimetry will be used with the upcoming launch of dedicated satellites, such as Imaging X-ray Polarimetry Explorer (IXPE, Weisskopf et al. 2021) and enhanced X-ray Timing and Polarimetry mission (eXTP, Zhang et al. 2019). For the BH X-ray binaries, in the absence of a solid surface of the compact object, accretion is the only source of the observed X-ray emission, and the obtained signatures can be directly connected to the innermost geometry of the accretion disk. For the NS X-ray binaries, such as accreting millisecond pulsars, the disk can serve as a source of constant energy-dependent polarimetric background, which has to be subtracted to obtain polarimetric profiles of pulsations (Viironen & Poutanen 2004).

Rees (1975) discussed the polarimetric signatures of the accretion disks with the dominant role of electron scattering (Chandrasekhar 1960; Sobolev 1963) in the Newtonian approximation. Observing such polarization signatures has been considered crucial for confirming the existence of an accretion disk in the first place (Lightman & Shapiro 1975). A number of deviations from this simple model were discussed, such as the role of true absorption, for which Monte-Carlo estimates have been presented by Lightman & Shapiro (1975) and the analytic results are discussed in Loskutov & Sobolev (1979) and Loskutov & Sobolev (1981), along with the alignment of the inner parts of the disk with the BH spin (Bardeen & Petterson 1975).

Yet another important addition to the early works on the disk polarization signatures is inclusion of the effects of special and general relativity (GR and SR; Connors & Stark 1977; Stark & Connors 1977; Pineault & Roeder 1977a,b). Relativistic aberration and light deflection lead to a rotation of the polarization angle (PA) and alter the viewing angle of different segments of the disk, the latter effect leading to a different polarization degree (PD) of a given segment; while the PD is Lorentz invariant, the difference with respect to the non-relativistic case appears because of the different angle between the disk normal and the observer direction. Frame dragging effects, relevant to spinning BHs described by Kerr (1963) metric, lead to the additional rotation of the polarization plane along the photon trajectory. Calculations of all these effects involve the parallel transport of the polarization vector along null geodesics, which requires calculations of the Walker & Penrose (1970) constant of motion. For every null geodesics leaving the disk, the polarization at infinity can be evaluated by solving linear equations for the components of the polarization vector (e.g., Connors & Stark 1977; Connors et al. 1980; Dovčiak et al. 2008; Li et al. 2009; Ingram et al. 2015).

In this implicit formulation, the calculation of Stokes parameters is related to the computationally expensive numerical integration of the equation of geodesic, known as the ray-tracing technique. With the launch of X-ray polarimetric satellites, the data fitting procedures will require fast routines for calculating the polarization from accretion disks. A high time gain can be achieved, for instance, by tabulating the observed polarization characteristics of the disk rings for different BH and orbit parameters, so that the minimization routines would only need to proceed via interpolation between the precomputed models. On the other hand, the tabulated models would have to be recalculated with any change to the local model, such as an alteration of the angular dependence of emission or radial energy dissipation profile. Similar problems are faced within the analysis and theoretical modeling of polarimetric images of accretion flows around supermassive BHs. To accelerate the calculations, Narayan et al. (2021) recently developed a fast code for evaluating polarimetric images in Schwarzschild metric, which is based on evaluation of Walker-Penrose constants and an approximation of the light bending formula from Beloborodov (2002). Understanding the role of different relativistic effects and their separation in the total image, however, remains elusive, as explicit analytical expressions for the PA and PD transformation – as a result of the joint action of GR and SR – have not been presented thus far.

Using such explicit formulae to directly relate the local and the observed PA and PD can also substantially speed up the minimization routine. Connors et al. (1980) give an explicit analytical expression for the rotation angle of the polarization plane due to solely SR effects (their Eq. (18)). Pineault (1977) pointed out that the total rotation of PA in Schwarzschild metric and after accounting for relativistic motions in the disk is not a simple sum of rotations caused by the relativistic motion in the flat space (${\chi }_{\mathrm{flat}}^{\mathrm{SR}}$) and the light bending (χ^GR).

In this work, we derive explicit analytical expressions for the rotation of the PA accounting for the relativistic motion of matter in the accretion disk and light bending in Schwarzschild metric. Corresponding formulae were previously derived for polarization properties of rapidly rotating NSs (Viironen & Poutanen 2004; Poutanen 2020a; Loktev et al. 2020). We use the laws of geometrical optics and exploit the fact that the light trajectories are flat in Schwarzschild metric, that is, the orientation of the polarization vector is fixed with respect to the angular momentum of the ray trajectory (see Pineault 1977). The light bending is computed using the recent analytical approximation derived in Poutanen (2020b), which allows us to achieve high accuracy, with the deviations from the exact solution being smaller than the existing or expected statistical errors of observations. This removes the need for using computationally expensive ray-tracing algorithms when computing the polarization of the accretion disks around NSs or low-spin BHs.

2. Polarized radiation from the accretion disk

2.1. Observed flux

We consider the emission of a disk surface element in an axially symmetric, geometrically (infinitely) thin accretion disk in Schwarzschild metric. We compute the observed Stokes parameters following the approach described in Poutanen (2020b) and repeat here the basic notations and formulae for completeness.

We choose a Cartesian coordinate system with the z-axis coinciding with the normal to the disk and x-axis lying along with the projection of the line of sight on the disk (see Fig. 1). In this system, the normal to the disk $\widehat{\mathit{n}}$, the unit vector in the observer direction $\widehat{\mathit{o}}$, and the radius-vector of the surface element $\widehat{\mathit{r}}$ have the following coordinates:

$$\begin{array}{cc}& \widehat{\mathit{n}}=(0,0,1),\hfill \\ \hfill & \widehat{\mathit{o}}=(sini,0,cosi),\hfill \\ \hfill & \widehat{\mathit{r}}=(cos\phi ,sin\phi ,0),\hfill \end{array}$$(1)

Fig. 1. Geometry of a flat accretion disk ring. The emitting surface element is described by the radius-vector r and velocity v. Vector $\widehat{\mathit{o}}$ pointing in the observer’s direction makes angle i with the disk normal.

where φ is the azimuth of the radius-vector measured from the x-axis and i is the disk inclination to the line of sight. The radius-vector of the surface element makes angle ψ to the line of sight:

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

The photon trajectories are planar in Schwarzschild metric, hence, the direction of the photon momentum close to the disk surface can be described as a linear combination of the observer vector and the radius-vector of the emission point:

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

where

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

The relation between the angles α and ψ is described by the light-bending integral or can be approximated using a simple analytical formula (Pechenick et al. 1983; Beloborodov 2002; Poutanen & Beloborodov 2006; Salmi et al. 2018; Poutanen 2020b).

We consider purely Keplerian motion of disk matter, with the velocity unit vector being parallel to the azimuthal vector:

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

The dimensionless velocity, β = v/c, relative to a static observer at the circumferential radius, r = R/R_S, measured in units of the Schwarzschild radius R_S = 2GM/c² of the central object of mass M, is (see e.g., Luminet 1979):

$$\begin{array}{c}\hfill \beta =\sqrt{\frac{u}{2(1-u)}},\end{array}$$(6)

where u = 1/r is the compactness. The corresponding Lorentz factor is

$$\begin{array}{c}\hfill \gamma =\frac{1}{\sqrt{1-{\beta }^2}}=\sqrt{\frac{1-u}{1-3u/2}}·\end{array}$$(7)

The photon momentum makes angle ξ with the velocity vector:

$$\begin{array}{c}\hfill cos\xi =\widehat{\mathit{v}}·{\widehat{\mathit{k}}}_0=\frac{sin\alpha }{sin\psi }\widehat{\mathit{v}}·\widehat{\mathit{o}}=-\frac{sin\alpha }{sin\psi }sinisin\phi ,\end{array}$$(8)

and angle ζ with the disk normal,

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

The Doppler factor is:

$$\begin{array}{c}\hfill \delta =\frac{1}{\gamma (1-\beta \widehat{\mathit{v}}·{\widehat{\mathit{k}}}_0)}=\frac{1}{\gamma (1-\beta cos\xi )}·\end{array}$$(10)

The unit vector of the photon momentum in the frame comoving with the surface element (fluid frame hereafter) is computed using 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}$$(11)

From this, we get the angle between the photon momentum and the local normal in the fluid frame

$$\begin{array}{c}\hfill cos{\zeta }^{\prime }=\delta cos\zeta .\end{array}$$(12)

The projection of the photon momentum on the disk plane in that frame has azimuth ϕ′ (as measured from the radial direction) and can be determined through:

$$\begin{array}{cc}& sin{\varphi }^{\prime }=\frac{{\widehat{\mathit{k}}}_0^{\prime }·\widehat{\mathit{\phi }}}{sin{\zeta }^{\prime }}=\delta \gamma \frac{cos\xi -\beta }{sin{\zeta }^{\prime }},\hfill \end{array}$$(13)

$$\begin{array}{cc}& cos{\varphi }^{\prime }=\frac{{\widehat{\mathit{k}}}_0^{\prime }·\widehat{\mathit{r}}}{sin{\zeta }^{\prime }}=\delta \frac{cos\alpha }{sin{\zeta }^{\prime }}·\hfill \end{array}$$(14)

Specific flux observed from the surface element at photon energy E is:

$$\begin{array}{c}\hfill \mathrm{d}{F}_{\mathrm{E}}=I_{\mathrm{E}}\mathrm{d}\mathrm{\Omega },\end{array}$$(15)

where dΩ is the solid angle occupied by the surface element on the observer’s sky, I_E is the specific intensity of radiation at infinity, which is related to the intensity in the fluid frame:

$$\begin{array}{c}\hfill I_{\mathrm{E}}={\left(\frac{E}{E^{\prime }}\right)}^3{I}_{\mathrm{E}\prime }^{\prime }({\zeta }^{\prime }).\end{array}$$(16)

Here, we assume that I′ does not depend on the azimuthal angle ϕ′ and for brevity omit its dependence on radius r. The total redshift factor (Luminet 1979; Chen et al. 1989),

$$\begin{array}{c}\hfill g\equiv \frac{E}{E^{\prime }}=\delta \sqrt{1-u}=\frac{\sqrt{1-3u/2}}{1+\beta sinisin\phi sin\alpha /sin\psi },\end{array}$$(17)

combines the effects of the gravitational redshift and the Doppler effect. The solid angle occupied by the surface element of area $\text{d}S=R\text{d}R\text{d}\phi /\sqrt{1-u}$ is given by:

$$\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{1-u}}\mathcal{L}cos\zeta .\end{array}$$(18)

Here, dS cos ζ is the projection of the element area on the plane of the sky, D is the distance to the source and

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

is the lensing factor (Beloborodov 2002; Poutanen 2020b). We obtain the final expression for the observed spectral flux in the form:

$$\begin{array}{c}\hfill \mathrm{d}{F}_{\mathrm{E}}(r,\phi )=g^3{I}_{\mathrm{E}\prime }^{\prime }({\zeta }^{\prime })\frac{\mathrm{d}Scos\zeta }{D^2}\mathcal{L}.\end{array}$$(20)

The total observed flux from the disk can then be obtained by integrating Eq. (20) over the radius and azimuthal angle:

$$\begin{array}{c}\hfill F_{\mathrm{E}}=\frac{R_{\mathrm{S}}^2}{D^2}{\displaystyle \underset{r_{\mathrm{in}}}{\overset{r_{\mathrm{out}}}{\int }}\frac{r\mathrm{d}r}{\sqrt{1-u}}{\int }_0^{2\mathrm{\pi }}\mathrm{d}\phi g^3\mathcal{L}cos\zeta I_{\mathrm{E}\prime }^{\prime }({\zeta }^{\prime }).}\end{array}$$(21)

For a chosen inclination i and for every pair (r, φ), we compute ψ using Eq. (2), which is then used to compute α and ℒ – either exactly from elliptical integrals (as described, e.g., in Salmi et al. 2018), or using analytical formulae (see Poutanen 2020b and Sect. 2.4). Then, ξ and ζ are obtained from Eqs. (8) and (9), respectively. Using the Keplerian velocity and the Lorentz factor given by Eqs. (6) and (7), we then get the Doppler factor δ from Eq. (10). Further, from Eqs. (12) and (17), we get the photon energy E′ and the zenith angle ζ′ in the fluid frame, which are needed to obtain $I_{{\text{E}}^{\prime }}^{\prime }({\zeta }^{\prime })$ and, if needed, the azimuth ϕ′ from Eqs. (13) and (14).

2.2. Computing polarized flux

We consider only the linear polarization, which is described by the three-component Stokes vector¹. The direction of the polarization vector changes along the photon trajectory because of the gravitational light bending and aberration. The key property of the Schwarzschild metric that allows us to derive simple analytical formulae for the PA rotation is that the angle between the polarization vector and the trajectory plane is conserved.

In order to describe polarized radiation in the fluid frame, it is convenient to introduce the polarization basis formed by the vector of the local normal and the photon momentum ${\widehat{\mathit{k}}}_0^{\prime }$:

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

The radiation field emitted at angle ζ′ from the disk surface with coordinates (r, φ) is then described by the Stokes vector

$$\begin{array}{c}\hfill {\mathit{I}}_{{\mathrm{E}}^{\prime }}^{\prime }({\zeta }^{\prime })=I_{{\mathrm{E}}^{\prime }}^{\prime }({\zeta }^{\prime })\left[\begin{array}{c}1\\ p({\zeta }^{\prime })cos2{\chi }_0\\ p({\zeta }^{\prime })sin2{\chi }_0\end{array}\right],\end{array}$$(23)

where p is the linear PD of radiation, which is invariant, namely, it does not change along photon trajectory. The PA χ₀ is defined as the angle between the polarization vector and the basis vector ${\widehat{\mathit{e}}}_1^{\prime 0}$ measured in the counterclockwise direction. The transformation of vector ${\mathit{I}}_{{\text{E}}^{\prime }}^{\prime }({\zeta }^{\prime })$ to the observed Stokes vector involves rotation by the corresponding Mueller matrix:

$$\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}$$(24)

where χ^tot is the rotation angle of the polarization plane. Combining this matrix with Eq. (20), we get the expression for the observed Stokes vector (in terms of fluxes):

$$\begin{array}{c}\hfill \mathrm{d}{\mathit{F}}_{\mathrm{E}}(r,\phi )=g^3\mathbf{M}(r,\phi ){\mathit{I}}_{{\mathrm{E}}^{\prime }}^{\prime }({\zeta }^{\prime })\frac{\mathrm{d}Scos\zeta }{D^2}\mathcal{L}.\end{array}$$(25)

Integrating over azimuth and radius, we get the observed Stokes vector from the entire disk surface:

$$\begin{array}{c}\hfill {\mathit{F}}_{\mathrm{E}}\equiv \left[\begin{array}{c}{F}_I\\ F_Q\\ F_U\end{array}\right]=\frac{R_{\mathrm{S}}^2}{D^2}{\displaystyle \underset{r_{\mathrm{in}}}{\overset{r_{\mathrm{out}}}{\int }}\frac{r\mathrm{d}r}{\sqrt{1-u}}{\int }_0^{2\mathrm{\pi }}\mathrm{d}\phi g^3\mathcal{L}cos\zeta \mathbf{M}(r,\phi ){\mathit{I}}_{{\mathrm{E}}^{\prime }}^{\prime }({\zeta }^{\prime }).}\end{array}$$(26)

If the angular distribution of radiation escaping from the disk surface (as measured in the fluid frame) does not depend on the azimuth ϕ′, but only on the zenith angle ζ′, the polarization vector is parallel to one of the basis vectors (22) and the Stokes U-parameter vanishes, namely, χ₀ = 0 or π/2. Equation (25) is transformed to:

$$\begin{array}{c}\hfill \mathrm{d}{\mathit{F}}_{\mathrm{E}}(r,\phi )=g^3{I}_{{\mathrm{E}}^{\prime }}^{\prime }(r,{\zeta }^{\prime })\left[\begin{array}{c}1\\ \pm p({\zeta }^{\prime })cos2{\chi }^{\mathrm{tot}}\\ \pm p({\zeta }^{\prime })sin2{\chi }^{\mathrm{tot}}\end{array}\right]\frac{\mathrm{d}Scos\zeta }{D^2}\mathcal{L}.\end{array}$$(27)

Integrating over r and φ, we easily get the total observed Stokes vector. We note that the transformation of the Stokes vectors does not depend on the nature of radiation escaping from the disk or on the assumption about the direction of the polarization vector in the fluid frame; it is fully determined by the rotation of the polarization plane, χ^tot. In order to obtain an analytical solution relating the emitted and observed Stokes vectors, we need to find an explicit formula for χ^tot.

2.3. Polarization angle

Here, we consider the physical reasons for the rotation of the polarization plane. The polarization vector is determined by the predominant direction of the electric field oscillations. The angle the polarization vector makes with respect to the projection of the disk normal, $\widehat{\mathit{n}}$, on the plane perpendicular to the photon trajectory is called the polarization angle (PA). As the photon propagates to the observer, the polarization vector remains perpendicular to the direction of motion. The light bending imposes changes to the photon propagation direction. In the Schwarzschild metric, each trajectory lies in a plane, and the normal to this plane remains the same along the whole trajectory. Hence, the polarization vector rotates around the normal to the trajectory plane. In a general case, the normal to the trajectory plane does not coincide with the disk normal, $\widehat{\mathit{n}}$, hence, the angle between the polarization vector and $\widehat{\mathit{n}}$ changes. Accordingly, this causes changes in the PA.

The Lorentz transformation between the lab frame and the fluid frame also constitutes a rotation of the photon momentum vector, and the polarization vector rotates correspondingly. As long as the normal $\widehat{\mathit{n}}$ is not the axis of the rotation, the PA acquires additional rotation.

The PA seen by the observer is the sum of the PA of emitted photon and its total rotation along the trajectory:

$$\begin{array}{c}\hfill \chi ={\chi }_0+{\chi }^{\mathrm{tot}}.\end{array}$$(28)

The PA rotation can be computed directly from the rotation of the polarization plane and the projection of the polarization vector onto the sky, as viewed by the observer. The rotation by the general relativistic (GR) light bending effects is (see derivation in Appendix A.2):

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

where

$$\begin{array}{c}\hfill \stackrel{\sim }{a}=\frac{1-cos\alpha cos\psi }{cos\alpha -cos\psi }·\end{array}$$(30)

This expression accounts only for GR (light bending) effects, while in reality the matter in the disk moves around the compact object at relativistic velocities, hence we need to account also for the SR effects, namely, for the relativistic aberration.

If the surface element moves at a velocity $\mathit{\beta }=\beta \widehat{\mathit{v}}$, the photon vector, ${\widehat{\mathit{k}}}_0$, as measured in the observers (static) frame, is related to the vector ${\widehat{\mathit{k}}}_0^{\prime }$, as measured in the fluid frame, via the Lorentz transformation. This transformation leads to the additional rotation of the PA, which can be written as (see Appendix A.3 for derivation)

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

We note that this expression takes into account the fact that the photon reaching the observer experiences light bending and, hence, the outgoing photon in the static reference frame close to the disk surface propagates along with the vector ${\widehat{\mathit{k}}}_0$ – and not $\widehat{\mathit{o}}$.

For the case of a flat space, when α = ψ, the expression above is reduced to:

$$\begin{array}{c}\hfill tan{\chi }_{\mathrm{flat}}^{\mathrm{SR}}=-\beta \frac{cosicos\phi }{sini+\beta sin\phi },\end{array}$$(32)

which is equivalent to Eq. (18) of Connors et al. (1980) for the case of the disk polarization perpendicular to meridional plane made by the photon momentum and the local normal, namely, χ₀ = π/2.

The total rotation of the PA is the sum of two rotations (Eqs. (29) and (31)):

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

This expression is different from the simple addition of rotations due to light bending, χ^GR, and rotation in flat space, ${\chi }_{\mathrm{flat}}^{\mathrm{SR}}$ (Eq. (32)), as originally noted in Pineault (1977). The reason for this is that the Lorentz transformation has to be applied to the photon momentum after accounting for the bending effect.

2.4. Light bending

The formulae presented above allow us to directly transform the Stokes parameters of radiation escaping from the surface element to the observed ones. To obtain the complete analytical transformation, we need to have an expression for α(ψ, r). The exact integral equation for the inverse function, ψ(α, r) (Pechenick et al. 1983; Beloborodov 2002; Poutanen 2020b), can be used to compute the original expression at high accuracy. For this, the function ψ(α, r) has to be tabulated at a dense grid of arguments and then the function α(ψ, r) can be obtained by interpolation (Salmi et al. 2018).

Alternatively, an approximate analytical expression (Beloborodov 2002) can be used:

$$\begin{array}{c}\hfill cos\alpha \approx 1-(1-u)y,\end{array}$$(34)

where y = 1 − cos ψ. This approximation gives ℒ = 1 for the lensing factor in Eq. (19). This relation was recently used in Narayan et al. (2021) to compute polarized images of an accretion disk in Schwarzschild metric observed nearly face-on. The approximate relation (34) is accurate only when ψ ≲ 90° and therefore cannot be used for high inclinations.

Higher accuracy can be achieved by using the recently proposed relation (Poutanen 2020b):

$$\begin{array}{c}\hfill cos\alpha \approx 1-y(1-u)\{1+\frac{u^2{y}^2}{112}-\frac{\mathit{euy}}{100}[ln(1-\frac{y}{2})+\frac{y}{2}]\},\end{array}$$(35)

where e is the base of natural logarithm. This relation can also be used for high inclinations. In particular, this relation gives a relative accuracy in α better than 0.06% at ψ < 120° and any radius exceeding 1.5 R_S. At these radii, the error does not exceed 0.2% for ψ < 162° (which is valid for any azimuthal angle φ and inclinations of i < 72°). Using Eq. (35), we get a similarly accurate analytical representation for the lensing factor:

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

Equations (35) and (36) can be used instead of the exact relations for fast and accurate calculations of the observed Stokes parameters, as well as the polarized images.

2.5. Images

The formalism developed above allows us to obtain images of the accretion disk on the sky in polarized light. For any point on the disk with polar coordinates (r, φ), we first compute the angle between the observer direction and the radius vector of the element ψ using Eq. (2). Then, using either exact or approximate relations for light bending, we get angle α, which is then used to get the impact parameter in units of R_S (Pechenick et al. 1983; Beloborodov 2002):

$$\begin{array}{c}\hfill b=\frac{r}{\sqrt{1-u}}sin\alpha .\end{array}$$(37)

Owing to the fact that photon trajectories are planar is Schwarzschild metric, we get the position angle Φ (measured counterclockwise from the projection of the disk axis on the sky) of the point where photon hits the plane of the sky:

$$\begin{array}{cc}& sin\mathrm{\Phi }=-\frac{sin\phi }{sin\psi },\hfill \end{array}$$(38)

$$\begin{array}{cc}& cos\mathrm{\Phi }=-\frac{cosicos\phi }{sin\psi }·\hfill \end{array}$$(39)

The error on b is fully determined by the accuracy of the relation α(ψ), while Φ is exact.

3. Applications

In this section, we show examples of calculations of the observed polarization signatures affected by the GR and SR effects. We start with the simple decomposition of the action of these effects on the PA rotation (χ^GR and χ^SR). We then proceed to a comparison between PA rotations obtained using the exact and approximate α(ψ) relations. This relation is also used when obtaining the angle of the outgoing photon to the disk normal, which ultimately affects the observed flux and PD. Finally, we show the calculations of the polarization signatures of the geometrically thin, optically thick (Shakura & Sunyaev 1973; Novikov & Thorne 1973) accretion disk.

3.1. Stokes vector for a narrow ring

Relativistic effects on polarization are expected to be of the highest importance at the smallest disk radii. First, we consider radiation produced at the radius of r = 3, corresponding to the innermost stable circular orbit for a Schwarzschild BH. In Fig. 2, we show the angles χ^SR, χ^GR, and χ^tot as a function of azimuth in the disk for different system inclinations. We define PA in the range [− 90° ;90°]. In general, both the GR and SR rotations decrease with increasing inclination. This leads to higher depolarizing effects at smaller inclinations, which is a known result (Dovčiak et al. 2008).

Fig. 2. Rotation angles of the polarization plane, χGR (red dot-dashed lines), χSR (blue dashed lines), ${\chi }_{\mathrm{flat}}^{\mathrm{SR}}$ (green dotted lines), and χtot (black solid lines) for a ring at r = 3 at different viewing angles (a) i = 30°, 60° (b), and i = 80° (c).

The curve χ^GR (red dot-dashed in Fig. 2) is (anti-)symmetric with respect to the azimuth φ = 0° and 180°. For r = 3, the rotation is highest at $\phi =180°\mp 17\stackrel{°}{.}6$. The position of the extrema can be obtained analytically if we apply the Beloborodov (2002) approximation for the light-bending angle to the expression (29) for the rotation angle (see Appendix A.2 for derivation):

$$\begin{array}{c}\hfill cos{\phi }_{\mathrm{ext}}^{\mathrm{GR}}=-\frac{r{sin}^{2}i+{cos}^{2}i}{rsini}·\end{array}$$(40)

This approximation gives ${\phi }_{\mathrm{ext}}^{\mathrm{GR}}=180°\mp 16°$.

The SR effects are computed for the Keplerian rotation given by Eq. (6). The curve χ^SR (blue dashed curve in Fig. 2) is not symmetric, which adds to the asymmetry of the total PA rotation χ^tot. At sufficiently large inclinations, χ^SR has a maximum very close to φ = 180°. Interestingly, the shape of the similar curve in flat space, ${\chi }_{\mathrm{flat}}^{\mathrm{SR}}$, is very different (compare dashed blue with dotted green curves). In flat space, the maximum rotation angle is reached at (see Eq. (32)):

$$\begin{array}{c}\hfill sin{\phi }_{max}^{\mathrm{SR},\mathrm{flat}}=-\frac{\beta }{sini}·\end{array}$$(41)

Obviously, the extrema do not exist if β > sin i and then ${\chi }_{\mathrm{flat}}^{\mathrm{SR}}$ is monotonic. A general formula for the rotation angle χ^SR does not allow us to obtain a simple analytical expression for the position of the extrema. We note that the numerator of Eq. (31) is a smooth function of φ, therefore, the rotation angle reaches maximum close to the azimuth where the denominator reaches minimum:

$$\begin{array}{c}\hfill sin{\phi }_{max}^{\mathrm{SR}}\approx -\frac{{sin}^2\zeta }{\beta sinisin\alpha /sin\psi }·\end{array}$$(42)

By noticing that the corresponding azimuth is close to 180° and that α and ψ vary slowly with φ, and putting ψ ≈ 90° +i, we get:

$$\begin{array}{c}\hfill sin{\phi }_{max}^{\mathrm{SR}}\approx -\frac{{cos}^2\alpha cosi}{\beta sinisin\alpha },\end{array}$$(43)

where we can substitute cos α = u − (1 − u) sin i using Beloborodov (2002) approximation for the light-bending angle. For r = 3 and i = 60° (80°), the exact calculations give the maximum at ${\phi }_{max}^{\mathrm{SR}}=189\stackrel{°}{.}7$ ($182\stackrel{°}{.}7$), while Eq. (43) gives $184\stackrel{°}{.}1$ ($182\stackrel{°}{.}2$).

The total rotation, χ^tot, is shown with the black solid lines in Fig. 2. Notably, for inclinations of i ≲ 30°, the PA rotation is monotonic with azimuth and crosses 90°. This means that the polarization direction significantly changes, making two full cycles as the azimuth varies from 0 to 360°. Although in Fig. 2a, we see that only χ^SR makes two cycles, the total effect is not entirely caused by the relativistic rotation. Namely, for small inclinations, all the light rays coming to the observer were originally emitted along the direction of radius-vectors, that is, ${\widehat{\mathit{k}}}_0·\widehat{\mathit{r}}>0$, for any azimuth. It means that the observer sees the ring from its outer parts as if every disk element was located between the observer and the BH in the flat space. This is taken into account in χ^SR, when we apply Lorentz transformation to these rays. The direction of matter motion, with respect to the observer, varies between these disk segments, leading to variations of the angles between vectors ${\widehat{\mathit{k}}}_0^{\prime }$ and ${\widehat{\mathit{k}}}_0$. Hence, the SR effects only enhance the PA rotation caused by the light bending. If GR effects alone do not cause the full rotation of PA (such as in Fig. 2a), the additional rotation coming from relativistic motion results in the crossing of the 90° boundary. For the same reason, we see two full rotations of PA for the case when the bright spot orbits the BH. For inclination 180° −i, both the GR and SR rotation angles change sign.

In Fig. 3, we show the total PA rotation χ^tot for different radii and inclinations, as well as the error on the PA computed using the analytical expression for light bending (35), as compared to the exact, numerical solution. We see that the relativistic effects decrease with radius, as expected. The largest rotation of the PA is observed at smaller inclinations. The difference between the PA computed using analytical formulae and numerically is always smaller than $\sim 0\stackrel{°}{.}6$.

Fig. 3. Upper panels: PA χtot as a function of azimuth for disk inclinations i = 30° (a), 60° (b) and 80° (c) and different r = 3 (black solid line), 5 (red dashed), 15 (green dot-dashed) and 50 (blue dotted). Lower panels: difference in the PA rotation computed using analytical formula (35) for light bending and via the exact calculations.

The topology of χ^tot is shown in Fig. 4 as a map on the accretion disk plane in the form of contours of constant values. It is similar to that shown in Fig. 3 of Dovčiak et al. (2008), but defined here on the range of angles [− 90° ,90°] (rather than the previously used [− 180° ,180°], which included two full cycles of PA). We see that for i = 30°, there exists a critical point (at r ≈ 4.06, φ ≈ 230°) where photons are emitted along the disk normal in the comoving frame (i.e., ζ′ = 0) resulting in zero PD and not defined PA. At large inclinations, the strongest variations in PA are seen around φ = 180°. Calculations using the analytical formula for light bending produce nearly identical pictures (see the dotted lines barely separable from the solid lines in Fig. 4).

Fig. 4. Contours of constant χtot at the accretion disk plane (r, φ). The observer is situated at inclinations i = 30° (panel a), 60° (panel b), and 80° (panel c). The coordinates are in units of RS. The polar angle φ is measured from the projection of the direction to the observer on the disk plane. The disk rotates in the counterclockwise direction. The innermost stable circular orbit at r = 3 is shown with a black circle. The dotted lines (wherever visible) show corresponding contours computed using the approximate formula for light bending (35).

The general topology of χ^tot is easier to understand by looking separately at the topology of χ^GR and χ^SR, which are presented in Figs. 5 and 6, respectively. The anti-symmetry of χ^GR in φ and very small rotation angles around φ ≈ 0 are clearly seen. Also at φ ≈ 180° the GR rotation is small at large radii, while closer to r = 3, the contours of constant values become highly packed, corresponding to very fast changes of PA with azimuth. At small inclinations the GR rotation of PA is large, while, for instance, at i = 80°, significant rotation happens only close to φ = 180°. For the χ^SR, we see the existence of the critical point at small inclinations outside r = 3. The peak in SR rotation is reached close to φ = 180° and the effect of SR rotation decreases with the increasing inclination, similarly to the GR rotation angle.

Fig. 5. Same as Fig. 4, but for χGR.

Fig. 6. Same as Fig. 4, but for χSR.

Interestingly, in the specific case of zero inclination i = 0, the expressions for the GR and the SR rotation angles can be considerably simplified:

$$\begin{array}{cc}& {\chi }^{\mathrm{GR}}(r,i=0)=\phi ,\hfill \end{array}$$(44)

$$\begin{array}{cc}& tan{\chi }^{\mathrm{SR}}(r,i=0)=-\beta tan\alpha ,\hfill \end{array}$$(45)

resulting in the total rotation:

$$\begin{array}{c}\hfill {\chi }^{\mathrm{tot}}(r,i=0)=\phi -\mathrm{atan}(\beta tan\alpha ).\end{array}$$(46)

Using Beloborodov (2002) approximation and assuming Keplerian rotation, we further get:

$$\begin{array}{c}\hfill {\chi }^{\mathrm{tot}}(r,i=0)\approx \phi -\mathrm{atan}\left(\sqrt{\frac{1+u}{2u}}\right)=\phi -\mathrm{atan}\left(\sqrt{\frac{1+r}{2}}\right).\end{array}$$(47)

For r → ∞, this transforms to χ^SR(r, i = 0) = φ − π/2. This is different by π from the expression for the rotation angle in flat space given by Eq. (32) in the limit i → 0:

$$\begin{array}{c}\hfill {\chi }_{\mathrm{flat}}^{\mathrm{SR}}(r,i=0)=\frac{\mathrm{\pi }}{2}+\phi ,\end{array}$$(48)

that is, both limits give the same position angle of the polarization (pseudo-)vector. On the other end, at the innermost stable orbit, r = 3, according to Eq. (47) the PA rotates by $\approx \phi -54\stackrel{°}{.}74$ (exact calculations give $\approx \phi -54\stackrel{°}{.}87$).

Light bending and aberration do not only rotate the polarization plane but also affect the flux and PD, as they alter the angle at which we see the surface element. In Fig. 7 we show, for the photon reaching the observer, the zenith angle ζ′ between the local normal and photon propagation direction in the fluid frame, as well as the difference between the angles computed numerically and using approximate analytical formula (35). This angle differs from the observer inclination i, because of the light bending and aberration. The photon rays originating in the parts of the disk behind the BH (around φ = 180°) experience the most pronounced bending, because the light trajectory lies above the BH at closer distances than the emission radius. This results in pronounced dips in ζ′ (seen in Fig. 7c), where the disk is seen more face-on (i.e. ζ′< i). The SR effects make the surface seen more face-on at places, where matter moves towards the observer, that is, φ > 180° and otherwise more edge-on for φ < 180°, where the matter of the disk moves away from the observer. This effect is noticeable for moderate and low inclinations, thus, the change of ζ′ in Fig. 7a is mainly due to SR effects. The error in ζ′, when using the approximate bending formula, grows with inclination, but the largest difference is still below $\sim 0\stackrel{°}{.}4$ (for i = 80°, see Fig. 7c).

Fig. 7. Upper panels: zenith angle ζ′ the photon reaching the observer makes to the local normal in the fluid frame as a function of the azimuth for disk inclinations i = 30° (panel a), 60° (b), and 80° (c) and different emission radii r = 3 (black solid line), 5 (red dashed), 15 (green dot-dashed), and 50 (blue dotted). Lower panels: difference between ζ′ computed using the analytical formula (35) for light bending and the one using exact formula.

There is a general decrease with regard to the importance of GR and SR effects with increasing distance from the BH, hence, the difference between the PAs computed using exact and approximate lensing formulae decreases with radius (Fig. 3a). However, counter-intuitively, as we see in Figs. 7b and c that the error in ζ′ grows with the radius of the ring around azimuthal angle φ = 180°. For high disk inclination, the impact parameter of photons propagating towards the observer – as well as the distance of closest approach to the central compact object – are much smaller than the radius of the corresponding ring. This leads to a larger, although still within $\sim 0\stackrel{°}{.}5$, inaccuracy in ζ′ for the approximation formula at high inclinations.

In Fig. 8, we show the quantity g³ℒ cos ζ, which is proportional to the observed flux from a surface element (see Eq. (20)). Here, we clearly see the effect of Doppler boosting for part of the ring at φ > 180° and deboosting for φ < 180°. For a high-inclination observer, gravitational lensing strongly amplifies the flux from the part of the disk behind the BH at φ ∼ 180°. The azimuthal dependence of the flux is not symmetric because of the increasing role of the Doppler boosting towards φ ∼ 270°. The relative error on the flux arising from the approximate bending and lensing formulae is largest for high inclinations, at azimuths of φ ∼ 180°, when the angle ψ is large. The largest error at r = 3 reaches 26% for i = 80°, while for i = 30° and 60°, it does not exceed 0.6%.

Fig. 8. Upper panels: combination g3ℒ cos ζ which is proportional to the observed flux from a disk surface element as a function of the azimuth for inclinations i = 30° (a), 60° (b), and 80° (c) and different emission radii r = 3 (black solid line), 5 (red dashed), 15 (green dot-dashed), and 50 (blue dotted). Lower panels: relative error (in %) on that combination computed using analytical formulae (35) and (36) for the light bending angle and lensing factor as compared to the exact calculations.

3.2. Polarization of accretion disk in the Schwarzschild metric

We compute the polarization signatures of the optically thick, geometrically (infinitely) thin accretion disk using the analytical formulae derived above. The intensity of the disk in the fluid frame depends on the radius r, energy E′, and the angle ζ′ between the photon vector and the disk normal, $cos{\zeta }^{\prime }={\widehat{\mathit{k}}}_0^{\prime }·\widehat{\mathit{n}}$. As an illustration, we consider the simple case of the standard accretion disk in Newtonian gravity (Shakura & Sunyaev 1973) and pure electron scattering atmosphere for polarization properties. Substituting relevant dependencies in Eq. (23), we obtain the Stokes vector in the fluid frame:

$$\begin{array}{c}\hfill {\mathit{I}}_{{\mathrm{E}}^{\prime }}^{\prime }(r,{\zeta }^{\prime })=\frac{1}{f_{\mathrm{c}}^4}{B}_{{\mathrm{E}}^{\prime }}\left(T_{\mathrm{c}}(r)\right)a_{\mathrm{es}}({\zeta }^{\prime })\left[\begin{array}{c}1\\ -p_{\mathrm{es}}({\zeta }^{\prime })\\ 0\end{array}\right],\end{array}$$(49)

where

$$\begin{array}{c}\hfill a_{\mathrm{es}}({\zeta }^{\prime })\approx 0.421+0.868cos{\zeta }^{\prime }\end{array}$$(50)

approximates the angular dependence of the outgoing intensity (Suleimanov et al. 2020),

$$\begin{array}{c}\hfill p_{\mathrm{es}}({\zeta }^{\prime })\approx 0.1171\frac{1-cos{\zeta }^{\prime }}{1+3.5cos{\zeta }^{\prime }}\end{array}$$(51)

approximates the angular dependence of the PD, and the PA of χ₀ = π/2 is used corresponding to the direction of the polarization vector perpendicular to the meridional plane (Chandrasekhar & Breen 1947; Chandrasekhar 1960; Sobolev 1949, 1963).

The spectral shape is described by the Planck function B_E′ of the color temperature, T_c = f_cT_eff, with the color correction factor assumed to be f_c = 1.7 (Shimura & Takahara 1995). For the radial dependence of the effective temperature, we use a simple expression for a standard Newtonian disk (Shakura & Sunyaev 1973):

$$\begin{array}{c}\hfill T_{\mathrm{eff}}^4(r)=\frac{3GM\dot{M}}{8\pi {\sigma }_{\mathrm{SB}}{R}^3}(1-\sqrt{\frac{3R_{\mathrm{S}}}{R}})=T_{\ast }^4{t}^4(u),\end{array}$$(52)

where

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

$$\begin{array}{cc}& t(u)={\left[u^3(1-\sqrt{3u})\right]}^{1/4}.\hfill \end{array}$$(54)

Figure 9 shows images of the accretion disk as viewed at three different inclinations. The colors reflect the bolometric intensity, which is given by the product g⁴t⁴(u)a_es(ζ′). The black and white contours represent the sky images of lines of equal radii r and equal azimuths φ, computed using exact bending relation and its approximation (35), respectively. The difference is only visible for i ≥ 60° at the accretion disk side behind the black hole, φ ≈ 180°, corresponding to ψ ≈ 150° −170° and photon trajectories lying close to the black hole. We also plot the polarization pseudo-vectors. The green sticks correspond to the exact calculations, while the blue one to the approximate light bending formula. The PA and PD are nearly identical for these two cases, as they are defined by χ^tot and ζ′, which are very well approximated by the analytical formulae (see Figs. 4 and 7). Small inaccuracy of the light bending formula leads to minor differences in the positions of the sticks on the plane of the sky, while the difference in their angles (PA) and lengths (PD) cannot be seen by eye.

Fig. 9. Images of an accretion disk as viewed at three inclinations i = 30° (a), 60° (b), and 80° (c). The coordinates on the sky are given in units of RS. Contours show the images of rings of equal radii r (from 3 to 15, with a step of 3) and equal azimuths φ (every 30°). Black solid contours correspond to the exact calculations of α(ψ), while the white dashed contours (almost fully overlapping with the black ones) are for the approximation (35). The colors reflect the logarithm of g4t4(u)aes(ζ′), which is proportional to the bolometric intensity. The blue sticks give polarization computed using approximate light bending expression (35), with their length being proportional to the PD and their position angle is given by the PA. Exact calculations shown by green sticks are nearly indistinguishable from the approximation.

Using Eqs. (26) and (27), we integrate over radius and azimuth to get the observed Stokes vector. The Stokes vector of the apparent luminosity is computed by multiplying the Stokes vector F_E by 4πD². The luminosity can be represented in a dimensionless form by scaling it to ${\sigma }_{\text{SB}}{T}_{*}^4{R}_{\text{S}}{}^2$, multiplying by photon energy and using the dimensionless photon energy x = E/kT_* as an argument:

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

where x = gx′. With such scaling of photon energy and luminosity, spectra of accretion disk become independent of the BH mass and accretion rate.

The positively defined observed PD is then

$$\begin{array}{c}\hfill p(x)=\frac{\sqrt{l_{x,Q}^2+l_{x,U}^2}}{l_{x,I}}·\end{array}$$(56)

The observed PA is computed as an argument of the complex quantity formed by the Stokes parameters:

$$\begin{array}{c}\hfill \chi (x)=\frac{1}{2}arg(l_{x,Q}+\mathrm{i}{l}_{x,U}).\end{array}$$(57)

In Fig. 10, we show the resulting spectral energy distributions of x l_x, PD (p) and PA (χ) for the accretion disk extending from r_in = 3 to r_out = 3000 in the Schwarzschild metric, as seen by a distant observer at different inclinations: i = 30°, 60°, and 80°. Figure 10a shows the dimensionless luminosity as a function of photon energy at three inclinations. It is much higher at low inclinations due to the strong beaming of radiation along with the normal in the electron-scattering dominated atmosphere. At higher inclinations, the Doppler effects shift the peak of emission to higher energies. The lower subpanel shows the relative error on luminosity when computations are done using approximate formulae for the light-bending angle (35) and lensing factor (36). In spite of the fact that the error on the lensing factor reaches nearly 30% at r = 3 for i = 80° (see Fig. 8), the integral flux has at most 1.5% relative error, owing to the fact that the error in ℒ at different radii has a different sign.

Fig. 10. Polarized spectrum of a standard accretion disk (3 < r < 3000) viewed at different inclinations: i = 30° (red solid line), 60° (blue dashed), 80° (black dotted) as a function of dimensionless photon energy x = E/kT*. Panel a: normalized luminosity x lx. Panel b: polarization degree. Panel c: polarization angle. The lower subpanels show the errors in the corresponding quantities when computations are done with approximate expressions Eqs. (35) and (36) for the bending angle and lensing factor. We note that the error on x lx is relative, while the errors on PA and PD are absolute.

The PD (see Fig. 10b) is much higher for high inclinations. At lower energies, the emission is dominated by large radii where relativistic effects are not important and PD follows the angular dependence given by p(i)≈p_es(ζ′). In agreement with previous findings (Stark & Connors 1977; Connors et al. 1980; Dovčiak et al. 2008), we observe the increase of the value of the rotation angle and the increasing role of depolarization effects at higher energies (see Figs. 10b and c), where the emission from the innermost radii is most important. Because the observed flux is dominated (due to the Doppler effect) by the part of the disk at φ ≈ 270°, where the rotation angle is large, the integrated rotation angle is high at small inclinations. The range of rotation angles is wide in this case (see Figs. 2–4), hence the depolarization effect is strongest at higher energies. As described, for example, in Dovčiak et al. (2008), the recovery of PD at energies higher than kT_* is caused by the fact that only a small area of the disk contributes to this range of energies and, therefore, the depolarizing effects are smaller. Calculations using approximate formulae for light bending and lensing factor give very accurate results in these cases, too. For example, the absolute error on PD is smaller than 0.03% for all inclinations and the error on the PA barely exceeds $0\stackrel{°}{.}1$.

Our calculations of one image shown in Fig. 9 were performed in 0.01 s when using the approximate analytical formula, while the calculations of one polarized spectrum from Fig. 10 take about 0.05 s. Similar calculations using interpolation in table α(ψ, r) obtained from the exact implicit formula for ψ(α, r) took about 1 s. On the other hand, direct calculations using a ray-tracing code ARCMANCER (Pihajoki et al. 2018) took about 10³ s to compute an image 400 × 400 pixels.

3.3. Extensions of the model

The fully analytical formalism developed above can naturally be applied to a number of problems on the production of polarized radiation near relativistic objects. First, the assumption that radiation is produced in a plane-parallel electron-scattering atmosphere can be relaxed. In this specific model, the Stokes vector defined in a polarization basis connected to the local normal contains only two non-zero parameters: I and Q. Properties of the local emission can be very different in other setups. For example, the accretion flows around BHs at low accretion rates are optically thin and polarization at low photon energies (radio to sub-millimeter for super-massive BHs and optical–infrared for X-ray binaries) may be related to synchrotron radiation of relativistic electrons (Poutanen & Veledina 2014; Yuan & Narayan 2014). The magnetic field direction in this situation will not likely be aligned with the local normal. In this case, the Stokes vector will contain three (or four, if we also consider circular polarization) components. In the X-ray domain, the PD and PA may also be very different from the optically thick case, as the photons are produced by Comptonization in a geometrically thick, optically thin hot flow. However, these deviations will change only the Stokes vector of the radiation escaping the disk, but will not affect the computational scheme of the observed polarized flux Stokes vector, which still will follow Eq. (26), with the expressions for the PA rotation angle remaining the same.

Another extension of the model is related to the velocity field of the accreting matter. In this paper, we assume that the gas velocity has only the azimuthal component, while, for example, an optically thin hot flow at low accretion rates, as seen by the EHT in M 87, may have a significant radial component (Event Horizon Telescope Collaboration 2021a,b; Narayan et al. 2021). We note that rotation of the polarization plane due to light bending χ^GR does not depend at all on the velocity field. To calculate the rotation χ^SR caused by special relativity effects in this case, we can use a more general expression (A.27) in the vector form.

A further extension of the model concerns the properties of the polarized flux if the emission is confined to a small region (hot spots) corotating with the flow (e.g., Pineault 1977; Connors et al. 1980). In this case, Eq. (25) for the observed Stokes vector as a function of emission azimuthal angle φ should be modified to account for the light travel time delays. It will depend on the photon arrival time (corresponding to the phase of the orbit, φ_obs):

$$\begin{array}{c}\hfill \mathrm{d}{\mathit{F}}_{\mathrm{E}}(r,{\phi }_{\mathrm{obs}})=g^3\mathbf{M}(r,\phi ){\mathit{I}}_{{\mathrm{E}}^{\prime }}^{\prime }({\zeta }^{\prime })\frac{\mathrm{d}{S}^{\prime }cos{\zeta }^{\prime }}{D^2}\mathcal{L},\end{array}$$(58)

where dS′ is the emission region area in the fluid frame, which is now multiplied by cos ζ′, instead of cos ζ. The relation between the emission azimuth φ and the arrival time can be obtained from the time delay integral (see e.g., Pechenick et al. 1983; Poutanen & Beloborodov 2006; Salmi et al. 2018).

4. Summary

In this work, we derive explicit analytical expressions describing the rotation of the polarization plane in Schwarzschild metric. We show that the total rotation is a sum of two effects: the rotation caused by the light bending (Eq. (29)), which is a pure GR effect, and the relativistic aberration. The latter effect acts on the PA in a way that is different from the relativistic aberration in flat space, because it has to be applied to the photons whose path will subsequently be altered by the bending (Eq. (31)).

We show test cases for the observed PA from a disk ring, as a function of azimuth for various emission radii and observer inclinations. In combination with the recently derived analytical formula for the light bending, it is possible to compute both the zenith angle of emission in the fluid frame (that will affect the observed flux and PD) and PA with an accuracy better than 1° for any inclination and emission radius. This opens up a possibility to produce high-precision polarized images of the accretion disks in the Schwarzschild metric in an unprecedentedly small computing time.

Integration over the disk surface gives the full observed Stokes vector as a function of energy. Utilizing the analytical formulae allows reducing the computing time compared to the ray-tracing calculations by a factor of a hundred giving accuracy better than 1% on the flux, 0.03% on the integrated PD, and about $0\stackrel{°}{.}1$ for the PA. Our analytical technique can be used for a detailed comparison of theoretical models with the polarimetric data on accretion disks around NSs and BHs.

Acknowledgments

This research (specifically in Sect. 3.2) was supported by the Russian Science Foundation grant 20-12-00364. We also acknowledge support from the Jenny and Antti Wihuri foundation (VL) and the Academy of Finland grants 309308 (AV), 322779 and 333112 (JP).

References

Appendix A: Rotation of the polarization plane

A.1. Rotation of PA for two polarization bases

For a photon moving along a unit vector $\widehat{\mathit{k}}$, we can define a polarization basis of two unit vectors perpendicular to each other as well as to $\widehat{\mathit{k}}$ as:

$$\begin{array}{c}\hfill \mathcal{B}(\widehat{\mathit{k}},\mathit{l})=\{{\widehat{\mathit{e}}}_1^{k,l},{\widehat{\mathit{e}}}_2^{k,l}\}=\{\frac{\mathit{l}-(\widehat{\mathit{k}}·\mathit{l})\widehat{\mathit{k}}}{|\widehat{\mathit{k}}\times \mathit{l}|},\frac{\widehat{\mathit{k}}\times \mathit{l}}{|\widehat{\mathit{k}}\times \mathit{l}|}\},\end{array}$$(A.1)

where l is an arbitrary vector not collinear with $\widehat{\mathit{k}}$. The angle the polarization plane of such a photon makes with the vector ${\widehat{\mathit{e}}}_1^{k,l}$ measured counterclockwise as viewed by the observer defines the polarization angle (PA) of the photon in this basis.

Two bases, $\mathcal{B}(\widehat{\mathit{k}},{\mathit{l}}_1)$ and $\mathcal{B}(\widehat{\mathit{k}},{\mathit{l}}_2)$, built around the same photon vector $\widehat{\mathit{k}}$, differ only by rotation around $\widehat{\mathit{k}}$ by an angle $\chi (\widehat{\mathit{k}},{\mathit{l}}_2\to {\mathit{l}}_1)$ (see Fig. A.1) that may be defined by trigonometric functions as:

$$\begin{array}{cc}\hfill cos\chi (\widehat{\mathit{k}},{\mathit{l}}_2\to {\mathit{l}}_1)& ={\widehat{\mathit{e}}}_2^{k,l_1}·{\widehat{\mathit{e}}}_2^{k,l_2}={\widehat{\mathit{e}}}_1^{k,l_1}·{\widehat{\mathit{e}}}_1^{k,l_2},\hfill \end{array}$$(A.2)

$$\begin{array}{cc}\hfill sin\chi (\widehat{\mathit{k}},{\mathit{l}}_2\to {\mathit{l}}_1)& ={\widehat{\mathit{e}}}_2^{k,l_1}·{\widehat{\mathit{e}}}_1^{k,l_2}=-{\widehat{\mathit{e}}}_1^{k,l_1}·{\widehat{\mathit{e}}}_2^{k,l_2}.\hfill \end{array}$$(A.3)

Fig. A.1. Polarization bases, $\mathcal{B}(\widehat{\mathit{k}},{\mathit{l}}_1)$ and $\mathcal{B}(\widehat{\mathit{k}},{\mathit{l}}_2)$, polarization vector, p, and corresponding PAs.

In this case, the PA (i.e., position angle of the electric vector of a photon moving along $\widehat{\mathit{k}}$) in basis 1 can be computed from the PA in basis 2 as χ₁ = χ₂ + χ. The angle χ can be expressed through its tangent

$$\begin{array}{cc}\hfill tan\chi (\widehat{\mathit{k}},{\mathit{l}}_2\to {\mathit{l}}_1)& =\frac{{\widehat{\mathit{e}}}_2^{k,l_1}·{\widehat{\mathit{e}}}_1^{k,l_2}}{{\widehat{\mathit{e}}}_2^{k,l_1}·{\widehat{\mathit{e}}}_2^{k,l_2}}=\frac{(\widehat{\mathit{k}}\times {\mathit{l}}_1)·({\mathit{l}}_2-(\widehat{\mathit{k}}·{\mathit{l}}_2)\widehat{\mathit{k}})}{(\widehat{\mathit{k}}\times {\mathit{l}}_1)·(\widehat{\mathit{k}}\times {\mathit{l}}_2)}\hfill \\ \hfill & =\frac{\widehat{\mathit{k}}·({\mathit{l}}_1\times {\mathit{l}}_2)}{{\mathit{l}}_1·{\mathit{l}}_2-(\widehat{\mathit{k}}·{\mathit{l}}_1)(\widehat{\mathit{k}}·{\mathit{l}}_2)}·\hfill \end{array}$$(A.4)

A.2. Rotation of PA due to light bending

We go on to quantify the effect of the gravitational light bending only. We assume that the Stokes vector of emitted radiation is defined at the disk surface in the lab (non-rotating) frame in the basis related to the disk normal $\widehat{\mathit{n}}$, that is, $\mathcal{B}({\widehat{\mathit{k}}}_0,\widehat{\mathit{n}})$. On the other hand, the observer defines PA in the basis formed by $\widehat{\mathit{n}}$ and the photon momentum at infinity $\widehat{\mathit{o}}$, that is, $\mathcal{B}(\widehat{\mathit{o}},\widehat{\mathit{n}})$. The way to find the rotation of the polarization plane between those two bases is to note that photon trajectories are flat in the Schwarzschild metric and the polarization vector is parallel transported along the trajectory. Since the photon trajectory plane is defined by the vectors $\widehat{\mathit{r}}$ and $\widehat{\mathit{o}}$ (as well as ${\widehat{\mathit{k}}}_0$), the PA measured in bases $\mathcal{B}(\widehat{\mathit{o}},\widehat{\mathit{r}})$ and $\mathcal{B}({\widehat{\mathit{k}}}_0,\widehat{\mathit{r}})$ is the same. Thus, the GR effect on PA consists of two terms only. The first one is the rotation of the polarization vector due to transformation from the basis $\mathcal{B}({\widehat{\mathit{k}}}_0,\widehat{\mathit{n}})$ to $\mathcal{B}({\widehat{\mathit{k}}}_0,\widehat{\mathit{r}})$, that is, $\chi ({\widehat{\mathit{k}}}_0,\widehat{\mathit{n}}\to \widehat{\mathit{r}})$. The second one is the rotation angle $\chi (\widehat{\mathit{o}},\widehat{\mathit{r}}\to \widehat{\mathit{n}})$ from basis $\mathcal{B}(\widehat{\mathit{o}},\widehat{\mathit{r}})$ to $\mathcal{B}(\widehat{\mathit{o}},\widehat{\mathit{n}})$. The total effect is then

$$\begin{array}{c}\hfill {\chi }^{\mathrm{GR}}=\chi ({\widehat{\mathit{k}}}_0,\widehat{\mathit{n}}\to \widehat{\mathit{r}})+\chi (\widehat{\mathit{o}},\widehat{\mathit{r}}\to \widehat{\mathit{n}}),\end{array}$$(A.5)

where

$$\begin{array}{cc}\hfill tan\chi ({\widehat{\mathit{k}}}_0,\widehat{\mathit{n}}\to \widehat{\mathit{r}})& =\frac{{\widehat{\mathit{k}}}_0·(\widehat{\mathit{n}}\times \widehat{\mathit{r}})}{({\widehat{\mathit{k}}}_0·\widehat{\mathit{r}})({\widehat{\mathit{k}}}_0·\widehat{\mathit{n}})}=-\frac{tanisin\phi }{cos\alpha },\hfill \end{array}$$(A.6)

$$\begin{array}{cc}\hfill tan\chi (\widehat{\mathit{o}},\widehat{\mathit{r}}\to \widehat{\mathit{n}})& =\frac{\widehat{\mathit{o}}·(\widehat{\mathit{r}}\times \widehat{\mathit{n}})}{(\widehat{\mathit{o}}·\widehat{\mathit{n}})(\widehat{\mathit{o}}·\widehat{\mathit{r}})}=\frac{tanisin\phi }{cos\psi },\hfill \end{array}$$(A.7)

where we used the facts that $\widehat{\mathit{n}}\times \widehat{\mathit{r}}=\widehat{\mathit{\phi }}$, $\widehat{\mathit{r}}·\widehat{\mathit{\phi }}=0$ and $\widehat{\mathit{r}}·\widehat{\mathit{n}}=0$, with the rest of scalar products obtained directly from Eqs. (1)–(5):

$$\begin{array}{cc}\hfill \widehat{\mathit{o}}·\widehat{\mathit{r}}& =cos\psi ,\hfill \end{array}$$(A.8)

$$\begin{array}{cc}\hfill {\widehat{\mathit{k}}}_0·\widehat{\mathit{r}}& =cos\alpha ,\hfill \end{array}$$(A.9)

$$\begin{array}{cc}\hfill \widehat{\mathit{o}}·\widehat{\mathit{n}}& =cosi,\hfill \end{array}$$(A.10)

$$\begin{array}{cc}\hfill {\widehat{\mathit{k}}}_0·\widehat{\mathit{n}}& =\frac{sin\alpha }{sin\psi }cosi=\frac{sin\alpha }{sin\psi }\widehat{\mathit{o}}·\widehat{\mathit{n}},\hfill \end{array}$$(A.11)

$$\begin{array}{cc}\hfill \widehat{\mathit{o}}·\widehat{\mathit{\phi }}& =-sinisin\phi ,\hfill \end{array}$$(A.12)

$$\begin{array}{cc}\hfill {\widehat{\mathit{k}}}_0·\widehat{\mathit{\phi }}& =-\frac{sin\alpha }{sin\psi }sinisin\phi =\frac{sin\alpha }{sin\psi }\widehat{\mathit{o}}·\widehat{\mathit{\phi }}.\hfill \end{array}$$(A.13)

We thus get the tangent of the sum of the two rotation angles

$$\begin{array}{cc}\hfill tan{\chi }^{\mathrm{GR}}& =tan(\chi ({\widehat{\mathit{k}}}_0,\widehat{\mathit{n}}\to \widehat{\mathit{r}})+\chi (\widehat{\mathit{o}},\widehat{\mathit{r}}\to \widehat{\mathit{n}}))\hfill \\ \hfill & =-\frac{(cos\alpha -cos\psi )(\widehat{\mathit{o}}·\widehat{\mathit{\phi }})(\widehat{\mathit{o}}·\widehat{\mathit{n}})}{cos\alpha cos\psi {(\widehat{\mathit{o}}·\widehat{\mathit{n}})}^2+{(\widehat{\mathit{o}}·\widehat{\mathit{\phi }})}^2}\hfill \\ \hfill & =\frac{-(\widehat{\mathit{o}}·\widehat{\mathit{\phi }})(\widehat{\mathit{o}}·\widehat{\mathit{n}})}{\stackrel{\sim }{a}[1-{(\widehat{\mathit{o}}·\widehat{\mathit{n}})}^2]+(\widehat{\mathit{o}}·\widehat{\mathit{r}})}\hfill \\ \hfill & =\frac{cosisin\phi }{\stackrel{\sim }{a}sini+cos\phi },\hfill \end{array}$$(A.14)

where

$$\begin{array}{c}\hfill \stackrel{\sim }{a}=\frac{1-cos\alpha cos\psi }{cos\alpha -cos\psi }\end{array}$$(A.15)

and we used the relation between directional cosines

$$\begin{array}{c}\hfill {(\widehat{\mathit{o}}·\widehat{\mathit{n}})}^2+{(\widehat{\mathit{o}}·\widehat{\mathit{r}})}^2+{(\widehat{\mathit{o}}·\widehat{\mathit{\phi }})}^2=1.\end{array}$$(A.16)

We see that χ^GR = 0 at φ = 0 where it changes the sign. Using Beloborodov (2002) approximation for the light bending angle, cos α = u + (1 − u) cos ψ, we get

$$\begin{array}{c}\hfill \stackrel{\sim }{a}=r+(r-1)cos\psi ,\end{array}$$(A.17)

which allows us to obtain a simple approximate expression for the rotation angle due to the light bending only:

$$\begin{array}{c}\hfill tan{\chi }_{\mathrm{B}02}^{\mathrm{GR}}\approx \frac{cosisin\phi }{rsini+(r{sin}^{2}i+{cos}^{2}i)cos\phi }·\end{array}$$(A.18)

The extrema are reached at

$$\begin{array}{c}\hfill cos{\phi }_{\mathrm{ext}}=-\frac{r{sin}^{2}i+{cos}^{2}i}{rsini}·\end{array}$$(A.19)

This approximation works whenever the absolute value of the right-hand-side of Eq. (A.19) is smaller than unity; for instance, for r = 3, the inclination should be i > 30°. For r = 3 and i = 60°, the extrema appear at $cos{\phi }_{\mathrm{ext},\mathrm{B}02}=-5/(3\sqrt{3})$, i.e. at φ_ext, B02 ≈ 180° ∓16°, corresponding to $tan{\chi }_{\mathrm{ext},\mathrm{B}02}^{\mathrm{GR}}=1/\sqrt{2}$ (or ${\chi }_{\mathrm{ext},\mathrm{B}02}^{\mathrm{GR}}=\pm 35\stackrel{°}{.}3$), while the exact calculations give ${\phi }_{\mathrm{ext}}\approx 180°\mp 17\stackrel{°}{.}6$ and ${\chi }_{\mathrm{ext}}^{\mathrm{GR}}\approx \pm 31\stackrel{°}{.}9$.

A.3. Rotation of PA due to aberration

Next we consider the effect of the special relativity (SR) on the PA. In the local frame, we assume the element of the disk symmetric with respect to the local normal vector to the disk surface, $\widehat{\mathit{n}}$. Therefore, the properties of the outgoing radiation in the local frame are only defined by the angle it makes with the vertical direction. An atmosphere model under such an assumption would give the polarization vector in $\mathcal{B}({\widehat{\mathit{k}}}_0^{\prime },\widehat{\mathit{n}})$ to be collinear to one of the basis vectors. The rotation of PA due to aberration can be decomposed into three rotations. First, from the basis $\mathcal{B}({\widehat{\mathit{k}}}_0^{\prime },\widehat{\mathit{n}})$ to $\mathcal{B}({\widehat{\mathit{k}}}_0^{\prime },\widehat{\mathit{v}})$. Second, from $\mathcal{B}({\widehat{\mathit{k}}}_0^{\prime },\widehat{\mathit{v}})$ to $\mathcal{B}({\widehat{\mathit{k}}}_0,\widehat{\mathit{v}})$. Third, from $\mathcal{B}({\widehat{\mathit{k}}}_0,\widehat{\mathit{v}})$ to $\mathcal{B}({\widehat{\mathit{k}}}_0,\widehat{\mathit{n}})$. We note that the second rotation is actually zero, because the three vectors ${\widehat{\mathit{k}}}_0^{\prime },{\widehat{\mathit{k}}}_0$, and $\widehat{\mathit{v}}$ lie in the same plane, with the photon momenta unit vectors connected via the Lorenz transformation (11). Thus, the total rotation of the PA due to the SR may be written as:

$$\begin{array}{c}\hfill {\chi }^{\mathrm{SR}}=\chi ({\widehat{\mathit{k}}}_0^{\prime },\widehat{\mathit{n}}\to \widehat{\mathit{v}})+\chi ({\widehat{\mathit{k}}}_0,\widehat{\mathit{v}}\to \widehat{\mathit{n}}).\end{array}$$(A.20)

Noting that $\widehat{\mathit{v}}·\widehat{\mathit{n}}=0$, each angle can be computed using Eq. (A.4):

$$\begin{array}{cc}& tan\chi ({\widehat{\mathit{k}}}_0^{\prime },\widehat{\mathit{n}}\to \widehat{\mathit{v}})=\frac{{\widehat{\mathit{k}}}_0^{\prime }·(\widehat{\mathit{n}}\times \widehat{\mathit{v}})}{({\widehat{\mathit{k}}}_0^{\prime }·\widehat{\mathit{n}})({\widehat{\mathit{k}}}_0^{\prime }·\widehat{\mathit{v}})},\hfill \end{array}$$(A.21)

$$\begin{array}{cc}& tan\chi ({\widehat{\mathit{k}}}_0,\widehat{\mathit{v}}\to \widehat{\mathit{n}})=\frac{{\widehat{\mathit{k}}}_0·(\widehat{\mathit{v}}\times \widehat{\mathit{n}})}{({\widehat{\mathit{k}}}_0·\widehat{\mathit{v}})({\widehat{\mathit{k}}}_0·\widehat{\mathit{n}})}·\hfill \end{array}$$(A.22)

Using the relation $\widehat{\mathit{v}}\times \widehat{\mathit{n}}=\widehat{\mathit{r}}$ valid in our case, and expression (11) for the Lorentz transformation, we get

$$\begin{array}{cc}& \frac{{\widehat{\mathit{k}}}_0^{\prime }·\widehat{\mathit{n}}}{{\widehat{\mathit{k}}}_0·\widehat{\mathit{n}}}=\frac{{\widehat{\mathit{k}}}_0^{\prime }·\widehat{\mathit{r}}}{{\widehat{\mathit{k}}}_0·\widehat{\mathit{r}}}=\delta =\frac{1}{\gamma (1-\beta {\widehat{\mathit{k}}}_0·\widehat{\mathit{v}})},\hfill \end{array}$$(A.23)

$$\begin{array}{cc}& {\widehat{\mathit{k}}}_0^{\prime }·\widehat{\mathit{v}}=\frac{{\widehat{\mathit{k}}}_0·\widehat{\mathit{v}}-\beta }{1-\beta {\widehat{\mathit{k}}}_0·\widehat{\mathit{v}}}·\hfill \end{array}$$(A.24)

The tangent of the sum of two angles is then

$$\begin{array}{cc}\hfill tan{\chi }^{\mathrm{SR}}& =\frac{({\widehat{\mathit{k}}}_0·\widehat{\mathit{r}})({\widehat{\mathit{k}}}_0^{\prime }·\widehat{\mathit{n}})({\widehat{\mathit{k}}}_0^{\prime }·\widehat{\mathit{v}})-({\widehat{\mathit{k}}}_0^{\prime }·\widehat{\mathit{r}})({\widehat{\mathit{k}}}_0·\widehat{\mathit{v}})({\widehat{\mathit{k}}}_0·\widehat{\mathit{n}})}{({\widehat{\mathit{k}}}_0^{\prime }·\widehat{\mathit{n}})({\widehat{\mathit{k}}}_0^{\prime }·\widehat{\mathit{v}})({\widehat{\mathit{k}}}_0·\widehat{\mathit{n}})({\widehat{\mathit{k}}}_0·\widehat{\mathit{v}})+({\widehat{\mathit{k}}}_0^{\prime }·\widehat{\mathit{r}})({\widehat{\mathit{k}}}_0·\widehat{\mathit{r}})}\hfill \\ \hfill & =\frac{({\widehat{\mathit{k}}}_0·\widehat{\mathit{r}})({\widehat{\mathit{k}}}_0·\widehat{\mathit{n}})[({\widehat{\mathit{k}}}_0·\widehat{\mathit{v}}-\beta )-(1-\beta {\widehat{\mathit{k}}}_0·\widehat{\mathit{v}})({\widehat{\mathit{k}}}_0·\widehat{\mathit{v}})]}{({\widehat{\mathit{k}}}_0·\widehat{\mathit{v}}-\beta )({\widehat{\mathit{k}}}_0·\widehat{\mathit{v}}){({\widehat{\mathit{k}}}_0·\widehat{\mathit{n}})}^2+(1-\beta {\widehat{\mathit{k}}}_0·\widehat{\mathit{v}}){({\widehat{\mathit{k}}}_0·\widehat{\mathit{r}})}^2}\hfill \\ \hfill & =\frac{\beta ({\widehat{\mathit{k}}}_0·\widehat{\mathit{r}})({\widehat{\mathit{k}}}_0·\widehat{\mathit{n}})({({\widehat{\mathit{k}}}_0·\widehat{\mathit{v}})}^2-1)}{{({\widehat{\mathit{k}}}_0·\widehat{\mathit{v}})}^2{({\widehat{\mathit{k}}}_0·\widehat{\mathit{n}})}^2+{({\widehat{\mathit{k}}}_0·\widehat{\mathit{r}})}^2-\beta ({\widehat{\mathit{k}}}_0·\widehat{\mathit{v}})[{({\widehat{\mathit{k}}}_0·\widehat{\mathit{r}})}^2+{({\widehat{\mathit{k}}}_0·\widehat{\mathit{n}})}^2]}\hfill \\ \hfill & =-\frac{\beta ({\widehat{\mathit{k}}}_0·\widehat{\mathit{r}})({\widehat{\mathit{k}}}_0·\widehat{\mathit{n}})}{1-{({\widehat{\mathit{k}}}_0·\widehat{\mathit{n}})}^2-\beta ({\widehat{\mathit{k}}}_0·\widehat{\mathit{v}})}=-\beta \frac{cos\alpha cos\zeta }{{sin}^2\zeta -\beta cos\xi },\hfill \end{array}$$(A.25)