Effect of black hole--plasma system on light beams

TL;DR

This paper investigates how a non-magnetized cold plasma around a Kerr black hole affects the propagation of light rays and the observed star distribution near the black hole, using a generalized deviation equation.

Contribution

It derives a generalized deviation equation for light rays in plasma around Kerr black holes and analyzes the plasma's impact on observed star distributions.

Findings

Plasma modifies the angular distribution of stars near the black hole.

The deviation of light rays is affected by plasma density and distribution.

Key features of star distribution are identified and discussed.

Abstract

In the paper we discuss propagation of light around Kerr black hole surrounded by non-magnetized cold plasma with infinite conductivity. For that purpose, we use equations for propagation of light rays obtained within Synge's approach in the approximation of geometrical optics. We derive equation of deviation of a ray propagating close to the reference ray, which is a generalization of the well-known Jacobi equation, and use it to calculate the modification of angular distribution of stars observed close to a black hole surrounded by plasma, compared to the uniform star distribution that would be seen without black hole or plasma. We place the observer in the equatorial plane of the Kerr black hole and try various choices of plasma distributions described by mathematically simple formulae. Key features of star distribution on a local sky near the black hole are identified and the influence of plasma on them is discussed.

Figures

Equations

$$\frac{1}{2}[g^{\mu \nu }+(1-n^2)V_m^{\mu }{V}_m^{\nu }]{\sigma }_{,\mu }{\sigma }_{,\nu }=0$$

$$H=\frac{1}{2}[{\mathbf{k}}^2+(1-n^2)(\mathbf{k}\cdot {\mathbf{V}}_m{)}^2]$$

$$n^2=1-\frac{{\omega }_{pl}^2}{{\omega }^2}$$

$$H=\frac{1}{2}[{\mathbf{k}}^2+{\omega }_{pl}^2(x)]$$

$$S[x]=\frac{1}{2}{\int }_{{\lambda }_1}^{{\lambda }_2}[{\mathbf{k}}^2-{\omega }_{pl}^2(x)]d\lambda$$

$$D\mathbf{k}=-\frac{1}{2}\mathcal{A}$$

$${\mathbf{k}}^2=-{\omega }_{pl}^2$$

$$\displaystyle\begin{array}[]{ll}\bm{U}&=\bm{k}-Q{\bm{V}}_{m}\\ \nabla_{\bm{U}}\bm{k}&=Q\bm{W}(\bm{k},.)+\bm{q},\\ \end{array}$$

$$D^2\mathbf{\xi }=R(\mathbf{k},\mathbf{\xi })\mathbf{k}$$

$$D^2\mathbf{\xi }=R(\mathbf{k},\mathbf{\xi })\mathbf{k}-\frac{1}{2}♯\mathcal{H}(\mathbf{\xi },.)$$

$$\mathbf{k}\cdot D\mathbf{\xi }+\frac{1}{2}\mathbf{\xi }\cdot \mathcal{A}=0$$

$$\tilde{\mathbf{\xi }}=\mathbf{\xi }+\Phi \mathbf{k}\Phi =-{\omega }_{pl}^{-2}\int \mathcal{A}\cdot \mathbf{\xi }d\lambda$$

$$\displaystyle\begin{array}[]{ll}Dx&=\bm{k}\\ D\bm{k}&=-(1/2)\bm{\mathcal{A}}\\ \end{array}\qquad\qquad\qquad\begin{array}[]{ll}D\bm{\xi}&=\bm{\mathcal{V}}\\ D\bm{\mathcal{V}}&=R(\bm{k},\bm{\xi})\bm{k}-(1/2)\sharp\bm{\mathcal{H}}(\bm{\xi},.)\end{array}$$

$$\displaystyle\begin{array}[]{ll}\qquad D\bm{\xi}_{(k)}=\bm{\mathcal{V}}_{(k)}\qquad\qquad\qquad\qquad D\bm{\mathcal{V}}_{(k)}=R(\bm{k},\bm{\xi}_{(k)})\bm{k}-(1/2)\sharp\bm{\mathcal{H}}(\bm{\xi}_{(k)},.)\end{array}$$

$$\displaystyle\begin{array}[]{ll}k_{t}&=\partial L/\partial\dot{t}=-E\\ k_{\varphi}&=\partial L/\partial\dot{\varphi}=\mathcal{L}\end{array}$$

$$\mathcal{L}=g_{\mu \nu }{\dot{x}}^{\mu }{\dot{\xi }}^{\nu }+{\Gamma }_{\mu \nu \rho }{\dot{x}}^{\mu }{\dot{x}}^{\nu }{\xi }^{\rho }-\frac{1}{2}{\xi }^{\mu }{\mathcal{A}}_{\mu }$$

$$\displaystyle\begin{array}[]{llll}P_{t}&=\partial\mathcal{L}/\partial\dot{t}&=(D\xi)_{t}+\Gamma_{\mu t\nu}k^{\mu}\xi^{\nu}&=\bm{\mathcal{V}}\cdot\partial_{t}+\bm{k}\cdot\nabla_{\bm{\xi}}\partial_{t}\\ P_{\varphi}&=\partial\mathcal{L}/\partial\dot{\varphi}&=(D\xi)_{\varphi}+\Gamma_{\mu\varphi\nu}k^{\mu}\xi^{\nu}&=\bm{\mathcal{V}}\cdot\partial_{\varphi}+\bm{k}\cdot\nabla_{\bm{\xi}}\partial_{\varphi}\end{array}$$

$$\mathcal{J}({\omega }_{obs},P,\mathbf{V}):(\bar{\vartheta },\bar{\phi })\mapsto (\vartheta (\lambda =-\infty ),\phi (\lambda =-\infty ))$$

$$dN=\frac{N_s}{4\pi }{n}_s(\bar{\vartheta },\bar{\phi })d\bar{\Omega }$$

$$n_s=\frac{d{\Omega }_{\infty }}{d\bar{\Omega }}$$

$$d{\Omega }_{\infty }=\sin {\vartheta }_{\infty }d{\vartheta }_{\infty }d{\phi }_{\infty }=\sin {\vartheta }_{\infty }|J|d\bar{\vartheta }d\bar{\phi }=\frac{\sin {\vartheta }_{\infty }}{\sin \bar{\vartheta }}|J|d\bar{\Omega }$$

$$n_s=\frac{\sin {\vartheta }_{\infty }}{\sin \bar{\vartheta }}\left|\frac{\partial {\vartheta }_{\infty }}{\partial \bar{\vartheta }}\frac{\partial {\phi }_{\infty }}{\partial \bar{\phi }}-\frac{\partial {\phi }_{\infty }}{\partial \bar{\vartheta }}\frac{\partial {\vartheta }_{\infty }}{\partial \bar{\phi }}\right|$$

$$e_{\hat{t}}=\gamma (B)(E_t-B{E}_{\phi })e_{\hat{r}}=E_r{e}_{\hat{\vartheta }}=E_{\vartheta }{e}_{\hat{\phi }}=\gamma (B)(-B{E}_t+E_{\phi })$$

$$\begin{array}[]{ll}r_{0}&=r_{P}\\ \vartheta_{0}&=\pi/2\\ k^{\hat{r}}_{0}&=\omega_{P}n_{P}c_{\bar{\vartheta}}\\ k^{\hat{\vartheta}}_{0}&=-\omega_{P}n_{P}s_{\bar{\vartheta}}s_{\bar{\varphi}}\\ k^{\hat{t}}_{0}&=\gamma(B)\omega_{P}[1-Bn_{P}s_{\bar{\vartheta}}c_{\bar{\varphi}}]\\ k^{\hat{\varphi}}_{0}&=\gamma(B)\omega_{P}[B-n_{P}s_{\bar{\vartheta}}c_{\bar{\varphi}}]\end{array}\qquad\qquad\qquad\begin{array}[]{ll}\xi^{\hat{t}}_{0}&=\xi^{\hat{r}}_{0}=\xi^{\hat{\vartheta}}_{0}=\xi^{\hat{\varphi}}_{0}=0\\ \mathcal{V}^{\hat{r}}_{0}&=Rb^{\tilde{r}}\\ \mathcal{V}^{\hat{\vartheta}}_{0}&=Rb^{\tilde{\vartheta}}\\ P_{\varphi}&=\bar{\omega}\gamma(B)Rb^{\tilde{\varphi}}\\ \varepsilon&=\alpha k^{\hat{t}}_{0}+\omega\bar{\omega}k^{\hat{\varphi}}_{0}\\ \mathcal{L}&=\bar{\omega}k^{\hat{\varphi}}_{0}\\ \end{array}$$

$$\zeta =\frac{1}{2}\left\{\left[({\xi }_{(2)}^{\hat{\phi }}+{\xi }_{(1)}^{\hat{\vartheta }})+i({\xi }_{(1)}^{\hat{\phi }}-{\xi }_{(2)}^{\hat{\vartheta }})\right]e^{i\sigma }+\left[({\xi }_{(1)}^{\hat{\vartheta }}-{\xi }_{(2)}^{\hat{\phi }})+i({\xi }_{(1)}^{\hat{\phi }}+{\xi }_{(2)}^{\hat{\vartheta }})\right]e^{-i\sigma }\right\}\equiv \frac{1}{2}(\kappa e^{i\sigma }+\mu e^{-i\sigma })$$

$${\delta }_{\pm }=\frac{1}{2}\frac{1}{n_{\infty }{\omega }_{\infty }}|\sqrt{{\kappa }_{\mathcal{V}}{\kappa }_{\mathcal{V}}^{\ast }}\pm \sqrt{{\mu }_{\mathcal{V}}{\mu }_{\mathcal{V}}^{\ast }}|$$

$$d{\Omega }_{\infty }=\frac{\pi }{n_{\infty }^2{\omega }_{\infty }^2}|{\mathcal{V}}_{(1)}^{\hat{\vartheta }}{\mathcal{V}}_{(2)}^{\hat{\phi }}-{\mathcal{V}}_{(2)}^{\hat{\vartheta }}{\mathcal{V}}_{(1)}^{\hat{\phi }}|$$

$$n_s=\frac{n_P^2{\omega }_P^2}{n_{\infty }^2{\omega }_{\infty }^2}|{\tilde{\mathcal{V}}}_{(1)}^{\hat{\vartheta }}{\tilde{\mathcal{V}}}_{(2)}^{\hat{\phi }}-{\tilde{\mathcal{V}}}_{(2)}^{\hat{\vartheta }}{\tilde{\mathcal{V}}}_{(1)}^{\hat{\phi }}|$$

$$\rho =\frac{K}{r^2+r_c^2}$$

$$\rho =\frac{K}{r^2+r_c^2}\exp \left[-\frac{(\vartheta -\pi /2{)}^2}{s^2}\right]$$

Full text

**Effect of black hole–plasma system on light beams

**

Matej Sárený111matej.sareny@gmail.com  and Vladimír Balek222balek@fmph.uniba.sk

*Department of Theoretical Physics, Comenius University, Bratislava, Slovakia

Abstract

In the paper we discuss propagation of light around Kerr black hole surrounded by non-magnetized cold plasma with infinite conductivity. For that purpose, we use equations for propagation of light rays obtained within Synge’s approach in the approximation of geometrical optics. We derive equation of deviation of a ray propagating close to the reference ray, which is a generalization of the well-known Jacobi equation, and use it to calculate the modification of angular distribution of stars observed close to a black hole surrounded by plasma, compared to the uniform star distribution that would be seen without black hole or plasma. We place the observer in the equatorial plane of the Kerr black hole and try various choices of plasma distributions described by mathematically simple formulae. Key features of star distribution on a local sky near the black hole are identified and the influence of plasma on them is discussed.

Keywords: relativistic geometric optics in plasma, ray equation, ray deviation, star distribution function

1 Introduction

The focus of our paper is the behavior of light rays in the vicinity of a rotating uncharged black hole (BH) in the presence of plasma. Since the discovery of Kerr metric [1] in 1963, the behavior of light rays around such BH has been intensively studied by many authors. Geodesics around the Kerr BH were investigated for example in [2, 3] and detailed discussion of spherical orbits can be found in [4]. For an observer at arbitrary position, some of the light rays are arriving “from the BH” (i.e. no rays from distant stars can arrive from that direction), forming the apparent “shadow” of the BH. The shape of the shadow was investigated in [5, 6]. For an observer close to the BH, the look of the night sky is distorted due to the bending of the light rays. Such photographs were simulated in [7, 8] using backward ray-tracing. Also, the appearance of the accretion disc around the BH, as seen by the observer outside the disc, was studied numerically [7, 9, 10, 11].

The inclusion of optical effects of the plasma around the BH, potentially present in the form of accretion disc, or a spherical cloud of dilute material gathered from cosmic environment, can be made by utilizing Synge’s geometrical optics for an isotropic medium with non-unit index of refraction [12]. Lately, this approach is being intensively applied in many papers. Effects of plasma around a non-rotating BH were investigated in [13, 14, 15, 16], including changes in deflection angle, shadow of the BH, magnification factor in gravitational lensing etc., and additional effects arising in Kerr metric were studied in [6, 17, 18, 19]. In particular, the influence of plasma on the apparent shadow of the BH was explored in [6, 19, 20]. As it turns out, among the properties of plasma only one function, the density of electrons, has an effect on light rays, while the 4-velocity of plasma drops out of the equations entirely. Besides, the trajectory of photon depends on its frequency due to the dispersive properties of plasma. The most commonly used “toy models” for plasma distribution are homogeneous plasma [14, 17, 19, 22, 23], singular and nonsingular isothermal sphere [13, 24] and a sphere with power law distribution [13, 16, 17, 19, 20, 22, 23, 24, 25]. Also, it would be interesting to consider disc-like distributions where the disc’s middle plane coincides with equatorial plane of the Kerr BH. Realistic models of accretion discs are not commonly used because of their intricate nature, involving discontinuities of various types and uncertainty in the vertical profile. The models also raise physical questions about the validity of geometrical optics in a disc due to its optical thickness in most cases, as well as its own radiation outshining any light rays from distant stars that would propagate through it. The theoretical work on accretion discs is reviewed in [26]; the most popular models include thin discs [27, 28, 29, 30], advection dominated discs [31, 32] and Polish doughnuts [33].

The body of the paper is divided into four sections. In the second section we rewrite Synge’s equations into Lagrangian form, obtaining equations of motion (EoMs) for the light ray in plasma with such parametrization that ${\dot{x}}^{\mu}$ coincides with the wave 4-vector. We build on this to derive ray deviation equation (RDE), which is a generalization of equation of geodesic deviation known from textbooks on general relativity. We also discuss the conservation laws for the case of Kerr BH surrounded by stationary axially-symmetrically distributed non-gravitating plasma. In the third section we introduce the star distribution function (SDF), which measures the angular dependence of star density as viewed by an observer influenced by light bending in constrast to an observer who sees homogeneous sky. In the fourth section we develop tools for calculating the SDF numerically, following techniques from [7, 34, 35]. Furthermore, we discuss technical details of our calculations, including observer specification, initial conditions, plasma distribution “toy models” and some specifics of our code. In the fifth section we discuss various features of SDF as obtained in our calculations.

Throughout the paper we use the signature $(-+++)$ with Greek spacetime indices and Latin spatial indices. Our other notations follow the textbook [36]: we use the symbols $d$ for exterior derivative, $\mathcal{L}$ for Lie derivative and $\nabla$ for the metric-compatible, torsion-free covariant derivative, and adopt “music notation” $\flat$ and $\sharp$ for abstract lowering and rising of indices. In index notation, these operators are $(df)_{\mu}=f_{,\mu}$, $({\cal L}_{\bm{U}}\bm{V})^{\mu}={V^{\mu}}_{,\nu}U^{\nu}-{U^{\mu}}_{,\nu}V^{\nu}$, $(\nabla_{\bm{U}}\bm{T})^{\rho\ldots\sigma}_{\mu\ldots\nu}=T^{\rho\ldots\sigma}_{\mu\ldots\nu;\lambda}U^{\lambda}\equiv(T^{\rho\ldots\sigma}_{\mu\ldots\nu,\lambda}+\Gamma^{\rho}_{\tau\lambda}T^{\tau\ldots\sigma}_{\mu\ldots\nu}+\ldots-\Gamma^{\tau}_{\mu\lambda}T^{\rho\ldots\sigma}_{\tau\ldots\nu}-\ldots)U^{\lambda}$, $(\flat\bm{V})_{\mu}=g_{\mu\nu}V^{\nu}$ and $(\sharp\alpha)^{\mu}=g^{\mu\nu}\alpha_{\nu}$, where $f$ is scalar, $\bm{U}$, $\bm{V}$ are vectors, $\alpha$ is covector, $\bm{T}$ is tensor of an arbitrary rank and $\Gamma^{\mu}_{\nu\kappa}$ are coefficients of affine connection. (We restrict ourselves to the expressions used in the article, therefore we do not write how $d$ and $\cal L$ act on general tensors, just define their action on scalars and vectors respectively.) We also utilize the geometrized unit convention $G=c=1$, and when considering Kerr geometry we set $M$ (the mass of Kerr black hole) $=1$.

2 Equations of motion for light rays in the presence of plasma

2.1 Equation of motion for a single ray

Geometrical optics for isotropic media is developed in the book by Synge [12] in Hamilton – Jacobi formalism. The basic object is the eikonal $\sigma$, which satisfies the equation

[TABLE]

where $n$ is refraction index of the medium and $\bm{V}_{m}$ is its 4-velocity. In terms of the eikonal, the wave 4-vector at any given point is defined as ${\bm{k}}=\sharp d\sigma$. For an observer carrying a clock, a ruler and a protractor described by an orthonormal tetrad $\bm{e}_{\hat{a}}$, the components $k^{\hat{0}}=-\bm{e}_{\hat{0}}\cdot\bm{k}$ and $k^{\hat{i}}=\bm{e}_{\hat{i}}\cdot\bm{k}$ give us the frequency and the wave vector respectively, hence $\nabla_{{\bm{V}}_{m}}\sigma$ is minus the frequency as observed in the medium rest frame.

The characteristics of the eikonal equation are called light rays and are given by the Hamilton equations following from the Hamiltonian

[TABLE]

Since the Hamiltonian is independent on the parameter of the ray $\lambda$, it must be constant, and from its definition it follows that its only physically acceptable value is $H=0$.

The formula for the index of refraction of cold non-magnetized plasma with infinite conductivity is (see for example [37])

[TABLE]

The plasma frequency $\omega_{pl}$ can be expressed in terms of the number density of free electrons measured in the plasma rest frame $N(x)$ as $\omega^{2}_{pl}=e^{2}N(x)/(\varepsilon_{0}m_{e})$, and as seen from the formula for $n$, it places a lower limit on the frequency $\omega$ of the wave measured in the plasma rest frame. The resulting Hamiltonian for light rays in plasma is

[TABLE]

Thus, the light ray propagating in plasma behaves in the approximation of geometrical optics as a particle with non-zero rest mass that varies from point to point.

The form of the plasma Hamiltonian is familiar from analytical mechanics and could have been obtained from the Lagrangian $L=(1/2)[\bm{k}^{2}-\omega_{pl}^{2}(x)]$, hence the correct action for a light ray is

[TABLE]

where we have chosen the parameter $\lambda$ in such a way that $\bm{k}=d/d\lambda$ (the parameter is not the proper time, the parametrization ensures that $\bm{k}$ is wave 4-vector). Equation of motion for a light ray is obtained from the condition $\delta S=0$ for an arbitrary variation such that the boundary points remain stationary. This yields

[TABLE]

where $\bm{\mathcal{A}}$ is the gradient of $\omega_{pl}^{2}$, $\bm{\mathcal{A}}=\sharp d\omega_{pl}^{2}$, and $D=\nabla_{\bm{k}}$. In index notation, the equation reads $Dk^{\mu}/d\lambda\equiv k^{\nu}{k^{\mu}}_{;\nu}=-(1/2){\omega^{2}_{pl}}{}^{,\mu}$. As we can see, the ray deviates from the geodesic in the direction of maximal decrease of plasma density. The conservation law $H=const$ following from the absence of $\lambda$ in the Lagrangian allows only for $H=0$ due to the eikonal equation. The law describes the normalization of the wave 4-vector,

[TABLE]

In the special case with $\omega_{pl}=const$ we get geodesic motion of a massive particle with constant rest mass. Cyclic coordinates give rise to additional conservation laws.

The analogy between a massive particle and a photon propagating through homogeneous plasma was noted in [13, 38]. In the latter paper, the authors also show that in cold plasma the transversal waves with dispersion relation (6) are accompanied by a longitudinal wave whose frequency with respect to the plasma rest frame is $\omega_{pl}$ irrespective of the wave number.

For an arbitrary index of refraction, the ray equations can be easily obtained by varying the Lagrangian $L=\bm{k}\cdot\bm{U}-H$, $\bm{U}=d/d\lambda$, with respect to $\bm{k}$ and $x$. For the Hamiltonian $H=(1/2)[\bm{k}^{2}+(1-n^{2})\omega^{2}]$, $\omega=-\bm{k}\cdot\bm{V}_{m}$, we find

[TABLE]

where $Q=(1-n^{2})\omega-n(\partial_{\omega}n)_{x}\omega^{2}$, $\bm{q}=n(\sharp dn)_{\omega}\omega^{2}$ and $\bm{W}=\sharp\nabla\flat{\bm{V}}_{m}$ (in index notation, ${W_{\mu}}^{\nu}={(V_{m})_{\mu}}^{;\nu}$). For cold plasma $Q=0$ and $\bm{q}=-\bm{\mathcal{A}}/2$, hence we arrive once again at equation (5).

2.2 Equation of motion for infinitesimally close rays

Geodesic deviation

A powerful tool to describe a light beam (a bundle of nearby light rays) in vacuum is the Jacobi equation a.k.a. equation of geodesic deviation

[TABLE]

The equation describes evolution of the connecting vector $\bm{\xi}$ pointing from the reference geodesic to a geodesic infinitesimally close to it. The vector $\bm{k}$ is the wave 4-vector, however, not timelike as before but null, $\bm{k}^{2}=0$, and $R$ is the curvature operator, $R(\bm{A},\bm{B})=\nabla_{\bm{A}}\nabla_{\bm{B}}-\nabla_{\bm{B}}\nabla_{\bm{A}}-\nabla_{[\bm{A},\bm{B}]}$. The derivation of this equation can be found for example in [40] or [36]. The connecting vector is supposed to be Lie-transported along the reference geodesic, $\mathcal{L}_{\bm{k}}\bm{\xi}\equiv[\bm{k},\bm{\xi}]=0$. The advantage of the Jacobi equation is that it needs to be integrated only along the reference geodesic, but we gain also information about the neighboring geodesics from it.

Ray deviation

Let us consider the action for rays in plasma (4), but for a nearby ray $y^{\mu}=x^{\mu}+\varepsilon\xi^{\mu}$. The generalization of the Jacobi equation in the presence of plasma reads

[TABLE]

where $\bm{\mathcal{H}}$ is the covariant Hessian matrix of $\omega_{pl}^{2}$, $\bm{\mathcal{H}}=\nabla d\omega_{pl}^{2}$; or in components, ${\mathcal{H}}_{\mu\nu}=(\omega^{2}_{pl})_{;\mu\nu}$. The derivation of this equation is discussed in the appendix A. The conservation law due to the absence of $\lambda$ in the Lagrangian governing the evolution of $\xi$ yields

[TABLE]

Actually, the conservation law implies just that the expression on the left hand side is constant, but the differentiation of the constraint (6) restricts the constant to be zero. If we choose an arbitrary observer with 4-velocity $\bm{V}$, we can decompose Jacobi vector into its spatial and time component with respect to that observer, and assuming $\bm{V}\neq C\bm{k}$, the spatial part can further be decomposed into longitudinal and transverse part w.r.t. the directional vector of the reference ray $\bm{k}$ (or alternatively, the spatial part of this vector w.r.t. $\bm{V}$). The transverse part can be, under certain conditions, used to measure the thickness of a narrow light beam.

Light bundles in cold plasma were investigated also in [39]. In the approach developed there, light rays are mapped onto timelike geodesics in a conformally rescaled metric $\tilde{\bm{g}}=\omega_{pl}^{2}\bm{g}$, in which the connecting vector $\bm{\tilde{\xi}}$ obeys standard Jacobi equation (with the 4-velocity normalized to unity in the rescaled metric and the connecting vector Lie-transported along it, i.e. satisfying $[\bm{\tilde{k}},\bm{\tilde{\xi}}]=0$). The authors derive from the equation for $\bm{\tilde{\xi}}$ equations for optical parameters, describing the beam in a basis parallel transported in the conformal metric (Sachs basis). We, on the other hand, use for the description of the beam the connecting vector $\bm{\xi}$ which we determine directly from equation (10). The two vectors differ, because the parameters in the two theories, with the metrics $\bm{g}$ and $\tilde{\bm{g}}$, differ: parameter $\lambda$ in the former theory is related to parameter $\tilde{\lambda}$ in the latter theory by the formula $d\tilde{\lambda}=\omega_{pl}^{2}d\lambda$. As a result, vector $\tilde{\bm{\xi}}$ can be expressed in terms of vector $\bm{\xi}$ as

[TABLE]

Evolution of a narrow elliptical light beam

Consider a light beam with elliptical cross-section. EoMs for a reference ray and an adjacent ray comprise a system of 8 second order ordinary differential equations (ODE), or equivalently 16 first order ODE

[TABLE]

Equation of ellipse can be written as $\vec{r}=\vec{a}\cos\sigma+\vec{b}\sin\sigma$ with general vectors $\vec{a}$ and $\vec{b}$ (they can be made orthogonal by shifting the parameter $\sigma$ by the constant $\rho$ such that $\tan(2\rho)=2\vec{a}\cdot\vec{b}/(a^{2}-b^{2})$). To evolve the whole beam along the reference ray we plug this Ansatz into the equations for ray deviation, i. e. for an arbitrary point at the boundary of the beam we write $\bm{\xi}=\bm{\xi}_{(1)}\cos\sigma+\bm{\xi}_{(2)}\sin\sigma$, $\bm{\mathcal{V}}=\bm{\mathcal{V}}_{(1)}\cos\sigma+\bm{\mathcal{V}}_{(2)}\sin\sigma$ and demand the EoMs to be valid for every $\sigma$. As a result we obtain the equations for two ellipse-generating adjacent rays ($k=1,2$):

[TABLE]

Thus, determining the form of a beam requires integration of 24 equations (if we had at the beginning a beam with elliptical cross section).

2.3 Conserved quantities

Kerr geometry has two cyclic coordinates $t$, $\varphi$, associated with Killing vectors $\partial_{t}$, $\partial_{\varphi}$. As a result, in a stationary, axially symmetric plasma distribution where $N=N(r,\vartheta)$ we obtain two conserved quantities for a light ray described by the action (4),

[TABLE]

These suggestively denoted constants represent (up to a multiplicative constant with the dimension of $\hbar$) energy and angular momentum w.r.t. the axis of symmetry, as measured by the observer at rest at infinity. Additional conservation laws can be derived from the two-ray action (A-1), assuming once again stationary axially symmetric plasma distribution. The cyclic coordinates are explicitly visible when the Lagrangian appearing in the action is written in the coordinate frame:

[TABLE]

The four coordinates that are absent here, hence are cyclic, are $t,\varphi,\xi^{t},\xi^{\varphi}$. The second pair yields the conservation of $k_{t},k_{\varphi}$ again, while the first pair results in additional conserved quantities

[TABLE]

The coordinate-free formula is convenient for recalculation of the conserved quantities into an arbitrary basis. Note that the existence of the second pair of conserved quantities comes from the fact that we have two conserved quantities per ray.

3 Jacobi map from the local sky to the celestial sphere

As an example let us choose the Kerr BH surrounded by plasma with negligible gravitational influence. Relevant formulae for Kerr geometry are found in appendix B. Consider an observer outside of the event horizon, characterized by their position $P=(r_{0},\vartheta_{0},\varphi_{0})$ and 4-velocity $\bm{V}$; moreover, suppose they can observe the radiation only in a narrow band of frequencies centered around $\omega_{obs}$. The set of all directions in which the observer can look is a unit sphere $S^{2}$, which may be called the “local sky” of the observer. The stars located far away from both the observer and the black hole can be thought of as living on a sphere with large radius which we will call the “celestial sphere”. Denote the coordinates on the local sky $(\bar{\vartheta},\bar{\varphi})$ and the coordinates on the celestial sphere $(\vartheta_{\infty},\varphi_{\infty})$ and define Jacobi map $\mathcal{J}$ that assigns a point on the local sky to a point on the celestial sphere by evolving a ray with initial frequency $\omega_{obs}$, which arrives to the observer from the direction given by the angular coordinates $(\bar{\vartheta},\bar{\varphi})$, back in time according to the evolution equation (5):

[TABLE]

The angles $\vartheta,\varphi$ appearing here are Boyer-Lindquist coordinates of the evolved ray. The definition only works if the light ray actually reaches infinity; if it does not for some reason, e.g. because it arrives “from the BH”, the corresponding point of Jacobi map remains undefined.

Suppose the stars are distributed uniformly on the celestial sphere. Consequently, the number of stars seen in the solid angle $d\Omega_{\infty}$ is $dN=(N_{s}/4\pi)d\Omega_{\infty}$, where $N_{s}$ is the total number of stars on the night sky. On the local sky, the distribution of stars will not be uniform due to the bending of the light rays. Thus, it will have the form

[TABLE]

The formula defines the function $n_{s}$ – the star distribution function (SDF) representing relative increase in the angular density of stars in comparison to the uniform case. Rewriting $dN$ in terms of $d\Omega_{\infty}$ we obtain

[TABLE]

In the empty space without plasma, the SDF is identically equal to one. Using the Jacobi map we can relate the two solid angles to each other,

[TABLE]

where $J$ is the Jacobian of the Jacobi map. Therefore, the SDF in the presence of a black hole and plasma can be calculated as

[TABLE]

However, this formula is not very suitable for numerical calculations. First, it contains derivatives which are complicated to calculate numerically, since that would involve many integrations of the ray evolution equation. Second, it contains the sine in the denominator which can be equal to zero. The RDE offers also another way to calculate the SDF, as we will see later.

4 Numerical computation of SDF using LNRF

The following recipe is a modification of the procedure used by Thorne et al. [7], when calculating shots for the movie Interstellar. We have modified the procedure to allow for a non-unit refraction index. Also, we chose to do the calculations in the locally non-rotating reference frame (LNRF), rather than adopting phasors (complex numbers) and adapted reference frame used there. The goal is to find the SDF in Kerr geometry with various plasma distributions. We have written a program integrating previously described equations of ray deviation to calculate the SDF. The program is written in Python 3 and uses the integrator lsode for ODE from the package scipy.integrate.ode [41]. Some SDFs obtained by this method can be found in Appendix C.

4.1 Specifics of the program

We have numerically calculated the SDFs as seen by an observer orbiting the Kerr black hole in the equatorial plane, possibly in the presence of plasma. We have left the following parameters free: the Kerr parameter $a$, the radial coordinate of the observer $r_{0}$, the orbiting velocity of the observer as measured in the LNRF $B$, the observational frequency as measured by the observer $\omega_{P}$, and the plasma distribution $\omega_{pl}^{2}(r,\vartheta)$ (we assumed it is smooth enough, so that the light does not refract by a finite angle after passing an infinitesimal distance along its path). Due to the mirror symmetry w.r.t the equatorial plane, it was sufficient to calculate the star distribution in the lower hemisphere w.r.t. the equatorial plane (the surface $\bar{\varphi}\in[0,\pi]$, see figure 1).

At the observer’s local sky, we have set up a grid of angles from the interval $(\bar{\vartheta},\bar{\varphi})\in[0,\pi]\times[0,\pi]$. Due to the method of calculation used, the accuracy of the values obtained at individual elements of the grid was independent from the size of the grid, which determined only the resolution of the final figure. In our calculation we could set the size of the grid independently in $\bar{\vartheta}$ and $\bar{\varphi}$ directions. To obtain the SDF for each point of the grid, we have numerically integrated the RDE (twice) and then computed the value of the SDF using the formula (29). The numerical integration was the most time-consuming process, therefore the calculation time has increased roughly linearly with the grid size. To decide when (if at all) the ray reached infinity, we have defined the “numerical infinity” at $r_{\infty}=10,000$. We have also set the maximum value of the affine parameter allowed for a ray to $\lambda_{max}=5\times 10^{5}\omega_{P}^{-1}$. Rays that did not reach $r_{\infty}$ before the parameter has assumed this value, or entered the event horizon in the meantime, were assigned blue color in the final picture (meaning undefined $n_{s}$). Note that in vacuum there also exists a more sophisticated criterion determining when the rays escape to infinity, which can be found e.g. in [7].

We have displayed the SDFs in figures consisting of two circles, which represent the “inner” and “outer” hemisphere oriented towards and away from the black hole respectively. The values of SDF were color-coded into a hybrid linear ($n_{S}\in[0,3]$) and segmented ($n_{S}>3$) colormap, since SDF tends to diverge when approaching the BH’s shadow and the resulting large values in a small region around the shadow would mess up the color-coding of the rest of the picture on a purely linear colormap.

4.2 Initial conditions

Initial conditions for the reference ray

Consider an observer at an event $P$ orbiting in the equatorial plane of the Kerr black hole with the velocity with respect to the LNRF $-1<B<1$, and introduce two orthonormal bases along the observer’s worldline, the LNRF basis $(\bm{e}_{\hat{t}},\bm{e}_{\hat{r}},\bm{e}_{\hat{\vartheta}},\bm{e}_{\hat{\varphi}})$ and the observer basis $(\bm{E}_{t}=\bm{V},\bm{E}_{r},\bm{E}_{\vartheta},\bm{E}_{\varphi})$. The components of the 4-velocity of the observer $\bm{V}$ as measured in the LNRF are given by the special relativistic formula $V^{\hat{\mu}}=\gamma(B)(1,0,0,B)$, where $\gamma(B)=(1-B^{2})^{-1/2}$ is the Lorentz factor. Correspondingly, $e_{\hat{t}}^{a}=\gamma(B)(1,0,0,-B)$, which reflects the fact that the LNRF moves w.r.t. the observer $P$ with the velocity $-B$. The complete set of transformation relations is given by the boost formulae (perpendicular vectors do not change and the fourth vector is just a normalized combination of the original vectors orthogonal to the first one):

[TABLE]

In the instantaneous physical 3-space of the observer $P$, standard spherical coordinates $(\bar{\vartheta},\bar{\varphi})$ are introduced using the correspondence $(E_{\varphi},E_{\vartheta},-E_{r})\leftrightarrow(\epsilon_{x},\epsilon_{y},\epsilon_{z})$ (the zenith angle $\bar{\vartheta}$ is measured from the radial direction pointing to the black hole and the azimuthal angle $\bar{\varphi}$ is measured from the vector $E_{\varphi}$ in the clockwise direction w.r.t. $E_{r}$, see figure 1).

In the observer’s frame, the wave 4-vector of the reference ray is initialized as $k^{a}=\omega_{P}(1,-n_{P}\vec{n})$, where $\vec{n}=(-c_{\bar{\vartheta}},s_{\bar{\vartheta}}s_{\bar{\varphi}},s_{\bar{\vartheta}}c_{\bar{\varphi}})$ is the unit vector in the direction in which the observer is looking (the photon is heading in the opposite direction, hence the extra minus sign) and $n_{P}=\sqrt{1-\omega_{pl}^{2}/\omega_{P}^{2}}$ is the index of refraction at the point where the observer is located. After initialization, we have to integrate towards negative values of the affine parameter (into the past).

Initial conditions for the beam

All rays of the beam are focused at the point $P$ at $\lambda=0$ and spread into a circular beam, seen by the observer within the solid angle $d\Omega_{P}$ at $\lambda=-d\lambda$. Focusing of the beam implies initial conditions (ICs) $\bm{\xi}_{(k)}(0)=0$ and spreading of the beam into the solid angle $d\Omega_{P}$ suggests $\bm{\mathcal{V}}_{(k)}(0)=\bm{\dot{\xi}}_{(k)}(0)$ with $\bm{\xi}(-d\lambda)=-(\bm{\mathcal{V}}_{(1)}(0)\cos\sigma+\bm{\mathcal{V}}_{(2)}(0)\sin\sigma)d\lambda$. The spatial part of the 4-vectors $\bm{\mathcal{V}}_{(k)}$ is supposed to be two independent vectors $\vec{V}_{(k)}$ perpendicular to $\vec{n}$ and to each other (in the sense of the induced 3-metric w.r.t. $\bm{V}_{P}$), which have both the magnitude $R$, therefore the radius of the circle as measured from $P$ is333$\epsilon$ comes from our convention that the connecting vector points to an infinitesimally close point, $y^{\mu}=x^{\mu}+\epsilon\xi^{\mu}$ $\epsilon Rd\lambda$. The radius of the circle enters ICs via $\pi\epsilon^{2}R^{2}d\lambda^{2}=dr^{2}d\Omega_{P}=\omega_{P}^{2}n_{P}^{2}d\lambda^{2}d\Omega_{P}$, therefore $R=\epsilon^{-1}\omega_{P}n_{P}\sqrt{d\Omega_{P}/\pi}$. The time component of the vector $\bm{\mathcal{V}}$ is given by the condition $\bm{k}\cdot\bm{\mathcal{V}}=0$ (see (11)), and it vanishes, ${\mathcal{V}}^{t}_{(k)}(0)=0$, for our choice of spatial ICs.

Initialization of differential equations

At the point $P$, 10 dynamic variables entering the evolution equations must be initialized together with 3 required conserved quantities. Position of the central ray is initialized by the choice of the point $P$ somewhere in the equatorial plane (and above the event horizon), $(r_{0}$, $\vartheta_{0})=(r_{P},\pi/2)$. Vector $\bm{\xi}$ is initially zero in the frame of the observer $P$ with our choice of ICs, therefore it is zero in any other frame too. The wave 4-vector, for both the reference ray and the adjacent ray, needs to be transformed from the observer’s frame to LNRF, which can be done by using inverse transformation to the transformation (27). To determine the conserving constants we must further pass from LNRF to the coordinate basis; for that purpose, use (B-10) and the fact that covector components transform like basis vectors. Finally, notice that $P_{\varphi}=\mathcal{V}_{\varphi 0}$ thanks to our choice of ICs. Below we give the full list of ICs expressed in terms of the parameters $\omega_{P},\bar{\varphi}$, $\bar{\vartheta}$, $R=\epsilon^{-1}\omega_{P}n_{P}\sqrt{d\Omega_{P}/\pi}$ and the two ellipse-generating vectors $\vec{b}_{(k)}$, which are unit 3-vectors perpendicular to $\vec{n}(\bar{\vartheta},\bar{\varphi})$ (they can also be chosen to be perpendicular to each other, $\vec{b}_{(1)}\cdot\vec{b}_{(2)}=0$):

[TABLE]

where $\alpha$, $\omega$ and $\bar{\omega}$ are defined in the appendix B. Note that the EoMs (B) are invariant w.r.t. the rescaling $x\mapsto x$, $\bm{k}\mapsto K\bm{k}$, $\omega_{pl}\mapsto K\omega_{pl}$, $\lambda\mapsto K^{-1}\lambda$, $\bm{\xi}\mapsto Q\bm{\xi}$, $\bm{\mathcal{V}}\mapsto KQ\bm{\mathcal{V}}$. This means that we can freely choose the “typical scale” for $\bm{k}$ and $\bm{\xi}$ without the necessity to change the equations in any way. A good choice is $K=\omega_{P}$ and $Q=R/\omega_{P}$, since it removes the superfluous factors from the ICs.

4.3 Evaluation of SDF

Provided the plasma density gradient falls to zero sufficiently fast when $r\rightarrow\infty$, the reference ray 4-velocity at $\lambda=-\infty$ equals the wave 4-vector of radially incoming ray, $k^{\hat{\mu}}=\omega_{\infty}(1,-n_{\infty},0,0)$. The cross-sectional area of the beam is therefore measured in the $Span(\bm{e}_{\hat{\vartheta}},\bm{e}_{\hat{\varphi}})$ plane. Using the Thorne’s approach [7] we write the major and minor angular semiaxis of the ellipse as $\delta_{\pm}=-dY_{\pm}/dr$, where $Y_{\pm}$ is the major and minor linear semiaxis of the ellipse. Since $\bm{\xi}_{(1)}$ is in general not perpendicular to $\bm{\xi}_{(2)}$, it is more convenient (following [7]) to use the description in which the ellipse is encoded in the complex variable $\zeta=\xi^{\hat{\vartheta}}+i\xi^{\hat{\varphi}}$. It holds

[TABLE]

and, consequently, $Y_{\pm}=\frac{1}{2}|\sqrt{\kappa\kappa^{*}}\pm\sqrt{\mu\mu^{*}}|$ (the major and minor semiaxis arises as the constructive and destructive sum of the phasors respectively). Also, if we notice that at infinity $\xi_{(k)}^{\hat{\mu}}\approx{\mathcal{V}}_{(k)}^{\hat{\mu}}\lambda$ and $dr\approx-\omega_{\infty}n_{\infty}d\lambda$, we find

[TABLE]

where the subscript $\mathcal{V}$ means that we have to replace all $\bm{\xi}$’s with $\bm{\mathcal{V}}$’s in the definitions of the corresponding variables. The last step is evaluating $d\Omega_{\infty}=\pi\delta_{+}\delta_{-}$. From the expression for $\delta_{\pm}$ it follows

[TABLE]

and if we insert $d\bar{\Omega}=d\Omega_{P}$ into the definition of SDF and denote $\tilde{\mathcal{V}}=\mathcal{V}/R$, we obtain

[TABLE]

4.4 Plasma distribution

In our calculations we have used various “toy models” of the plasma distribution, all of which can be found in the literature except for the last one which we have proposed ourselves. The advantage of these models is their smoothness, as well as the fact that the derivatives of the plasma density can be expressed analytically. The most trivial case is vacuum with $\rho=0$. Another simple case is the distribution $\rho=const$, when the ray equation coincides with the equation for a massive particle in vacuum. A model often used in the literature [13, 24] is the non-singular isothermal sphere (NSIS), defined as

[TABLE]

where $K$ and $r_{c}$ are positive constants. The same model with $r_{c}=0$ is called the singular isothermal sphere [13, 24]. In Newtonian gravity, such sphere can be obtained as a solution for a static non-rotating self-gravitating gas cloud with the same temperature of the gas everywhere, supported against the collapse solely by the pressure gradient. Finally, as a disc-like distribution we have used a non-singular isothermal sphere modified by a function that flattens the distribution towards the equatorial plane,

[TABLE]

where the parameter $s$ determines the magnitude of flattening.

5 General features of SDFs

We have investigated four plasma distributions, namely no plasma, homogeneous plasma, NSIS (distribution (30)) and flattened NSIS (distribution (31)). The 4-velocity of the plasma was not important since it vanishes from the equations. We have chosen the locally non-rotating family of observers, since they filter out the special relativistic effects as much as possible and exist also near the event horizon (contrary to the static observers who become tachyonic inside ergosphere). In addition to that, we have limited ourselves to the observers located in the equatorial plane. The main properties of the SDFs obtained in this way are:

a) no plasma: it is well known that purely gravitational bending of light is achromatic (frequency independent), hence the SDF is independent of $\omega_{obs}$ in the case with no plasma. Far from the black hole, the SDF is close to one in almost all directions, it changes rapidly only in the vicinity of the black hole shadow, decreasing to zero at Einstein’s ring and diverging at the edge of the shadow. The steep rise of the SDF near the shadow is caused by the multiple image generation by the BH. In fact, the SDF oscillates between zero and large values infinitely many times when approaching the shadow, but in the calculations presented here we do not see this behaviour due to low resolution. As we approach the BH, the behavior of the SDF stays qualitatively the same, but its features become more pronounced. Non-zero spin of the BH introduces asymmetry in SDF w.r.t. the azimuthal direction, with slightly greater values of SDF in the prograde direction. In particular, near the BH the outward-facing maximum of SDF is visibly shifted from $\bar{\vartheta}=\pi$ in the prograde direction if the BH is rotating (see fig. 2). As the observer approaches the event horizon, the BH shadow becomes ever more pronounced, eventually encompassing the whole sky.

b) homogeneous plasma: in this case, the light rays behave as massive test particles in Kerr geometry. Recalling that $k^{\hat{\mu}}=\omega_{obs}(1,-n_{obs}\vec{n})$, we see that $n_{obs}$ plays the role of launching velocity and $\omega_{obs}$ is, modulo factor $1/\omega_{pl}$, the launching Lorentz factor. If at least some rays are to reach infinity, the launching velocity must exceed the escape velocity, $\omega_{obs}>\omega_{min}>\omega_{pl}(r_{obs})$. Thus, the SDF acquires a new property: the night sky is not visible below $\omega_{min}$. In case of rotating BH or a boosted frame, the escape energy exhibits non-trivial dependence on the launching direction. Finally, in the limit of large launching velocity, $\omega_{obs}\gg\omega_{pl}(r_{obs})$, we regain the vacuum picture.

c) non-singular isothermal sphere: the $r$-dependence of NSIS leads to radial acceleration, pushing the light rays away from the BH (see fig. 10 c)). This fights the gravitational pull of the black hole and may result in erasing the black hole shadow from the night sky entirely.

d) flattened NSIS: the disc-like shape of this distribution causes an additional deviation of the acceleration of rays from radial direction: the rays are pushed away from the equatorial plane (see fig. 10 d)). As a consequence, the number of rays the observer receives when looking in equatorial plane directions is hugely enhanced and SDF assumes large values in these directions, fading fast to nearly zero in most of the non-equatorial-plane directions. This model is complicated enough to cause that in the hemisphere facing the BH a quite complex pattern emerges as a result of the struggle between attraction from gravitation and repulsion from plasma (see fig. 5). In both NSIS distributions, the effects of plasma are negligible when $\omega_{obs}\gg\omega_{pl}(r_{obs})$ just like in the constant distribution.

The influence of plasma on the shadow of BH was previously studied as well: analytic approach was cultivated in [19, 21], while numerical calculations were performed in [20]. Analytic approach was restricted to the specific class of plasma distributions, $\omega_{pl}^{2}=\rho^{-2}[f_{r}(r)+f_{\vartheta}(\vartheta)]$ [19]. Numerical calculations, not subject to this constraint, produced the shape of the shadow with purely radial, power law and exponential distributions. Out of these works, our situation is most similar to [19], where the authors placed, like us, the observer at finite distance from the BH ($r=5$ in their case). All papers also included description of the vacuum shadow, either directly or as a high frequency limit of the plasma shadow. Near-extremal BH vacuum shadow exhibits a typical feature – flattened left hand side of the shadow, which was present in all these papers and can be also seen in fig. 2; this feature is called “NHEKline” in some papers, e.g. in [21]. In homogeneous plasma, the increase in the size of the shadow with decreasing frequency was noticed in [19], and is also visible between figures 2 and 3. Paradoxically, distributions that decrease as we are getting further from the BH and vanish as $r$ approaches infinity have the opposite effect on the shadow: its size decreases with decreasing frequency. This feature can be observed in all three papers, as well as in figures 4 and 5. A more detailed discussion of this topic can be found in [19]. Finally, all figures, ours as well as those in the papers [19, 20, 21], suggest that plasma rounds off the BH shadow and suppresses its typical NHEKline signature at frequencies near $\omega_{pl}$.

Besides the plasma distribution, the form of SDF is influenced by other circumstances such as BH rotation, observer’s location and observer’s velocity. Rotation of black hole introduces asymmetry between prograde and retrograde direction at any spatial point outside the axis of rotation, and as a result, in case with no plasma the rays evolving forward (backward) in time are pulled in prograde (retrograde) direction. Thus, in our case the rays become more squashed together in prograde direction and we obtain greater SDF there (this is a naive opinion on the issue which seems to be confirmed by numerical calculation). An observer boosted by a significant fraction of speed of light in a certain direction sees the light rays pulled towards each other in that direction due to the special relativistic aberration of light, hence if the plasma is absent, SDF rises there. Also the BH shadow is pulled in that direction. These effects are modified in the presence of plasma. The situation gets more complicated, since the evolution of rays is now dependent on $\omega_{obs}$ in addition to their direction. Therefore the set of rays used to calculate SDF for unboosted observer does not coincide with the set of rays used to calculate SDF for the boosted observer (this fact can be ignored in vacuum due to achromaticity of the light bending). Thus, there is no direct relation between SDF of two observers moving w.r.t. each other. Let us just note that at large distances from the BH, effects of gravitational pull are severely suppressed except for a narrow region around the BH shadow, while the optical effects of plasma may still be important depending on the plasma distribution. However, in the distributions primarily concentrated around the BH these effects are also suppressed everywhere except nearby the BH shadow. On the other hand, if the observer is very close to the BH, the gravitational pull dominates on a large portion of the night sky, leaving only a small part of the sky with well defined SDF that can be influenced by the presence of plasma. The influence is not easy to predict, the analysis being complicated by the fact that we do not have an analogue of Carter constant for a generic plasma distribution, only for distributions fulfilling the additional condition cited above, which was found in [19] and further discussed or applied in calculations in [6, 17, 18].

6 Conclusion

We have discussed some theoretical aspects of general relativistic geometrical optics in plasma, in particular, we have provided a simple derivation of RDE which was, to our knowledge, not presented in the literature before. We have also discussed conservation laws and normalization constraints, writing them in a covariant way. We have then introduced the concept of SDF and calculated numerically the SDFs in various cases. To do that, we have converted the EoMs for the reference ray, as well as the RDE, into a set of first order ODE. This allowed us to calculate the SDF point by point, in contrast to constructing a network of reference rays and computing the derivatives numerically, which would be certainly less exact, especially for rays that undergo several revolutions around the BH before escaping to infinity. The SDF was calculated for several different plasma distributions: the vacuum (no plasma), the homogeneous plasma, the NSIS and the flattened (disc-like) NSIS. We have discussed also the dependence of SDF on some other parameters, namely the distance of the observer from the BH, the observational frequency and the spin of the BH. In future work there are several directions which appear worth to follow. First, it would be interesting to calculate the damping effect of plasma on the light beam, the most dominant effect being probably the Thompson scattering. Second, we would like to look at the SDF of an observer outside the equatorial plane, which requires only slight modification of the initial conditions. Finally, it would be interesting to calculate the light curves of a distant star seen by a distant observer on the other side of the BH, when the incoming light is influenced simultaneously by the gravitational field of the BH and the refraction index of the surrounding plasma. Obviously, this is only interesting if the plasma itself is not shining significantly, so that it does not outshine the light from the distant star.

Acknowledgements

The work was supported by the grant VEGA 1/0985/16. M. S. is also thankful to the Comenius University for its support via the grants UK/159/2018 and UK/105/2017.

Appendix A Efficient variation procedure for the two-ray action

Consider the action (4) for a reference ray $x^{\mu}(\lambda)$ and a neighboring ray $y^{\mu}(\lambda)$ lying in the close vicinity of it. We can write $y^{\mu}=x^{\mu}+\varepsilon\xi^{\mu}$, where $\varepsilon$ is a positive parameter $\ll 1$ and $\bm{\xi}$ is the connecting vector. Expand the action for the neighboring ray to the first order in $\varepsilon$,

[TABLE]

When rewriting the integral, we started from the equality 444The equality just states that the connecting vector is Lie-transported along the reference ray, $\mathcal{L}_{\bm{k}}\bm{\xi}|_{x(\lambda)}=0$. The same formula is used in the derivation of Jacobi equation.$[\bm{k},\bm{\xi}]=0$ at $x(\lambda)$, which holds because $\bm{\xi}$ closes by definition a 4-sided polygon with $\bm{k}$, and we used its consequence, the equality $\nabla_{\bm{\xi}}\bm{k}=\nabla_{\bm{k}}\bm{\xi}$, implied by vanishing of the torsion tensor. Demanding both $x$ and $y$ to obey the EoM for a ray suggests a variational principle for the second term on the r.h.s., with the boundary values of both $x$ and $\bm{\xi}$ fixed. Thus, we can introduce the action for ray deviation:

[TABLE]

The first term is the action for the equation of geodesic deviation a.k.a. Jacobi equation (similar term is proposed e.g. in [42, 43]), while the second term introduces plasma on the scene. Let us write the $x$-variations as $x+\varepsilon\bm{\mathcal{X}}$ and the $\bm{\xi}$-variation as $\bm{\xi}+\varepsilon\bm{\eta}$, where $\bm{\mathcal{X}},\bm{\eta}$ are vector fields vanishing at the boundaries. Since we need to go only to the first order in $\varepsilon$, we can discuss these variations separately. First consider the $\bm{\xi}$-variation. It holds

[TABLE]

where we introduced the effective equality $\hat{=}$ in which the terms with the total derivative (terms of the form $\nabla_{\bm{k}}(...)$) are dropped. The variation yields yet again the evolution equation for the reference ray,

[TABLE]

Now let us proceed to the variation w.r.t. $x$. We will use the covariant variation rather than the standard one555The difference is that the covariant variation will also cause a change in $\bm{\xi}$, but this change is dependent on the $x$-variation and will drop out eventually due to the ray equation.. Also, it is useful to realize that the commutation relations allow us to swap $\nabla_{\bm{k}}\bm{\xi}$ for $\nabla_{\bm{\xi}}\bm{k}$ and $\nabla_{\bm{k}}\bm{\mathcal{X}}$ for $\nabla_{\bm{\mathcal{X}}}\bm{k}$. The geometric part of the Lagrangian varies as follows:

[TABLE]

where $\bm{A}_{geom}$ is an auxiliary vector defined as $\bm{A}_{geom}=\sharp[(\nabla_{\bm{k}}\bm{k})\cdot(\nabla\bm{\xi})]$. The Riemann curvature operator $R(\bm{U},\bm{V})=\nabla_{\bm{U}}\nabla_{\bm{V}}-\nabla_{\bm{V}}\nabla_{\bm{U}}-\nabla_{[\bm{U},\bm{V}]}$ was introduced in order to rearrange the covariant derivatives, and in the course of calculation the identity $\bm{k}\cdot R(\bm{\mathcal{X}},\bm{k})\bm{\xi}=\bm{\mathcal{X}}\cdot R(\bm{k},\bm{\xi})\bm{k}$, which follows from the symmetry of the Riemann tensor w.r.t. the interchange of the pairs of indices, was used. Let us now determine the $x$-variation of the plasma term:

[TABLE]

where $\bm{A}_{plasma}=\sharp[(1/2)(\sharp d\omega_{pl}^{2})\cdot(\nabla\bm{\xi})]$. Collecting all terms we get:

[TABLE]

where we have isolated the auxiliary term $\bm{A}=\bm{A}_{geom}+\bm{A}_{plasma}$, which appears as a consequence of using the covariant variation and vanishes on the extremal due to the ray equation (A-2). Finally, using the definition of the covariant Hessian matrix we arrive at the ray deviation equation,

[TABLE]

Appendix B Explicit formulae for RDE integration in Kerr geometry in

LNRF

Equations of motion for adjacent ray in LNRF

Here we give the “totally unpacked version” of equations of motion:

[TABLE]

The equations marked with the bullet are algebraic rather than differential and come from the conservation laws. The unknown symbols on the right hand sides encode Kerr geometry and plasma distribution and are defined in the following two paragraphs.

The complete set of equations comprises 10 differential and 4 algebraic equations which are to be solved twice for every reference ray (since the ellipse is generated by two vectors ${\bm{\xi}}_{(k)}$)666Alternatively, we can solve $6+16$ equations simultaneously, should it prove to be faster. In fact, this is how we proceeded.. There are also optional evolution equations for temporal and azimuthal coordinate, $\dot{t}=g^{tA}k_{A}$ and $\dot{\varphi}=g^{\varphi A}k_{A}$, which may be added to the list to obtain the complete set of 16 equations of motion. Note, that Boyer-Lindquist coordinates are ill-defined on north/south pole. This can cause problems, if we try to integrate the system of rays with the reference ray crossing one of the poles (rays with $\mathcal{L}=0$). The equations for reference ray cope with this problem by letting $\vartheta$ escape from the interval $[0,\pi]$, but the equations for neighbouring ray blow up because of term $P_{\varphi}/\bar{\omega}$, so it is generally good idea to avoid integrating these rays when calculating SDF. Rays closely passing the poles may also occasionally suffer from numerical problems, see e.g. blue spots on fig. 5 near values $\bar{\varphi}=\pi/2,3\pi/2$.

Kerr geometry

Kerr metric in the Boyer–Lindquist coordinates, in the notations adopted from Thorne et al. [7], is:

[TABLE]

where $\Delta=r^{2}-2r+a^{2}$, $\rho^{2}=r^{2}+a^{2}\cos^{2}\theta$ and

[TABLE]

In this form an orthonormal covector basis $(\bm{\omega}^{\hat{t}},\bm{\omega}^{\hat{t}},\bm{\omega}^{\hat{\vartheta}},\bm{\omega}^{\hat{\varphi}})$ arises naturally. The corresponding vector basis $(\bm{e}_{\hat{t}},\bm{e}_{\hat{r}},\bm{e}_{\hat{\vartheta}},\bm{e}_{\hat{\varphi}})$, carried by zero angular momentum observers, is called locally non-rotating reference frame (LNRF). The two bases are:

[TABLE]

Connection forms $\bm{\omega}_{\hat{\mu}\hat{\nu}}$ ($(\omega_{\hat{\mu}\hat{\nu}})_{\hat{\lambda}}=\Gamma_{\hat{\mu}\hat{\nu}\hat{\lambda}}$, where hatted indices refer to LNRF) are:

[TABLE]

where $\omega_{\hat{\nu}\hat{\mu}}=-\omega_{\hat{\mu}\hat{\nu}}$ and

[TABLE]

We have used the abbreviations $s=\sin\vartheta$ and $c=\cos\vartheta$. The formulae can be looked up in [44] or obtained directly by solving Cartan’s structure equations (to be found e.g. in [36]).

The LNRF curvature tensor $R_{\hat{\mu}\hat{\nu}\hat{\kappa}\hat{\lambda}}$, as computed in [44], is

[TABLE]

where the pairs of indices $(\hat{\mu}\hat{\nu})$ and $(\hat{\kappa}\hat{\lambda})$ assume values $(\hat{t}\hat{r})$, $(\hat{t}\hat{\vartheta})$, $(\hat{t}\hat{\varphi})$, $(\hat{r}\hat{\vartheta})$, $(\hat{r}\hat{\varphi})$ and $(\hat{\vartheta}\hat{\varphi})$ (in this order), and

[TABLE]

The matrix is symmetric because of the symmetry of $R_{\hat{\mu}\hat{\nu}\hat{\kappa}\hat{\lambda}}$ in the pairs of indices.

Plasma description

In our calculations we have restricted ourselves to stationary, axially symmetric plasma distributions, $\omega_{pl}^{2}=\omega_{pl}^{2}(r,\vartheta)$. Gradient of the function $\omega_{pl}^{2}$ is $d\omega_{pl}^{2}=\mathcal{A}_{\hat{r}}\bm{e^{\hat{r}}}+\mathcal{A}_{\hat{\vartheta}}\bm{e^{\hat{\vartheta}}}$, where $\mathcal{A}_{\hat{r}}=(\sqrt{\Delta}/\rho)(\omega_{pl}^{2})_{,r}$ and $\mathcal{A}_{\hat{\vartheta}}=(1/\rho)(\omega_{pl}^{2})_{,\vartheta}$. The Hessian matrix can be easily calculated using the formula

[TABLE]

The matrix is symmetric and has only 6 nonzero components:

[TABLE]

Appendix C Numerical results

To illustrate the general features of the star distribution function discussed in the text, let us present here some example figures with our results.