1 Introduction

Black holes are a miracle of gravity. Although their existence used to be controversial in the last century, they are now believed to play important roles in the life cycle of a massive star, in the center of a galaxy and in the theory of quantum gravity. During the recent decade, more evidence for black holes has been obtained [1,2,3], including the amazing shadow and ring of M87* imaged by the Event Horizon Telescope (EHT) at a wavelength of 1.3 millimeters [3]. In the center of M87, there is a supermassive spinning black hole embedded in a geometrically thick, optically thin accretion disk. The shadow and ring is formed by strong gravitational lensing of synchrotron emission from the hot plasma in the disk. In this scenario, the plasma plays the part of a light source and the central black hole works as a very strong lens.

In the future, hopefully images of higher resolution and of other black holes will be generated, then it will become necessary to study the formation of the shadow and ring in more details. For example, the success of the EHT [3] relied on the condition that the accretion disk of M87* is optically thin at millimeter wavelengths. In this way, the plasma in the disk is treated as a pure emission source, and the small opacity has negligible influence on the fuzzy image. However, for other black holes, in longer baselines [4] or at smaller wavelengths [5], the plasma cannot always be transparent, then it will become necessary to consider the scattering and absorption of photons by the plasma.

Besides opacity, the plasma has two other effects on the shadow and ring. First, as a refractive medium the plasma can change the deflection angle of light rays in the gravitational field of a black hole [6,7,8,9]. Second, the energy and momentum of the plasma particles can modify the gravitational field outside the black hole. For a Schwarzschild black hole surrounded by a spherically symmetric plasma, the former refractive effect has been investigated by Ref. [10]. In this paper, we aim to study the latter gravitational effect in a similar toy model. In our model, we will assume that the plasma is static and its density profile falls as \(r^{-3/2}\) outside the event horizon. This model is far from realistic, but it enables us to assess the gravitational effect of plasma particles following the method introduced in Ref. [10] and make a comparison with the refractive effect of the plasma medium.

The rest of the paper is organized as follows. In Sect. 2, we briefly review some useful formulae in Ref. [10] and establish our convention of notations. In Sect. 3, we extend the toy model of Ref. [10] by incorporating the plasma’s gravity into the metric. The density profile is cut off inside a sphere of certain critical radius. By specifying the cutoff radius, we consider two concrete models:

  1. (A)

    The plasma density profile is cut off at the event horizon;

  2. (B)

    The plasma density profile is cut off at the innermost stable circular orbit.

Applying the general formulae to model A, influences of the plasma on black hole shadows are investigated in Sect. 4. In Sect. 4.1 we work out the corrections of refractive and gravitational effects to radii of the photon sphere and the shadow. Inserting observational data of some active galactic nuclei (AGNs) [11, 12], we find these corrections are negligibly small, though the gravitational correction can overtake the refractive correction for AGNs more massive than about \(10^9M_{\odot }\). In Sect. 4.2, we illustrate the influences of refractive and gravitational effects on the trajectories of light rays and the observed intensity image. This is done by assigning exaggerated values to model parameters. For model B, a parallel investigation is done in Sect. 5. In Sect. 6 we summarize the our results and discuss the implications.

2 Previous results and convention of notations

We are interested in the shadow of a Schwarzschild black hole surrounded by a spherically symmetric plasma. In this situation, the shadow seen by a distant observer has the shape of a circular disk. To determine the radius of the shadow, one should study the equations of motion of light rays outside the black hole. This has been investigated in Ref. [10] by including the refractive effect of the plasma medium but neglecting the gravitational field of the plasma particles. In this paper, we will take the gravitational effect of the plasma particles into account. In the current section, as a preparation, we will collect some useful formulae from Ref. [10] (see also Ref. [13]) which hold for spherical spacetimes generally.

In general, a static spherical spacetime is described by a metric of the form

$$\begin{aligned} ds^2= & {} -A(r)c^{2}dt^{2}+B(r)dr^{2}\nonumber \\&+D(r)\left( d\theta ^{2}+\sin ^{2}\theta d\phi ^{2}\right) \end{aligned}$$
(1)

in which c represents the light velocity. In accordance with the symmetry, the electron number density N(r) is a function of radius only. If we denote the charge and mass of the electron as e and \(m_e\) respectively, then the plasma frequency is given by

$$\begin{aligned} \omega _p(r)^2=\frac{4\pi e^2}{m_e}N(r) \end{aligned}$$
(2)

in Gaussian units. It is useful to introduce the function

$$\begin{aligned} h(r)^2=\frac{D(r)}{A(r)}\left[ 1-A(r)\frac{\omega _p(r)^2}{\omega _0^2}\right] \end{aligned}$$
(3)

with a constant of motion \(\omega _0\). On the equatorial plane, the orbit equation of light rays can be written as

$$\begin{aligned} \frac{dr}{d\phi }=\pm \sqrt{\frac{D(r)}{B(r)}\left[ \frac{h(r)^2}{b^2}-1\right] } \end{aligned}$$
(4)

with the parameter b being a constant of motion. In terms of \(h(r)^2\), the radius of the outermost photon sphere \(r_\mathrm {ph}\) is the largest root of the equation

$$\begin{aligned} \frac{d}{dr}h(r)^2=0. \end{aligned}$$
(5)

For an observer at a distance of r, the opening angle of the shadow is determined by

$$\begin{aligned} \sin ^2\alpha _\mathrm {sh}=\frac{h(r_\mathrm {ph})^2}{h(r)^2}, \end{aligned}$$
(6)

and the radius-squared of the shadow is accordingly

$$\begin{aligned} R_\mathrm {sh}^2=r^2\sin ^2\alpha _\mathrm {sh}=\frac{r^2h(r_\mathrm {ph})^2}{h(r)^2}. \end{aligned}$$
(7)

More generally, it can be verified that in the image seen at a distance of r, photons with the same value of b form a circle of radius-squared

$$\begin{aligned} R^2=r^2\sin ^2\alpha =\frac{b^2r^2}{h(r)^2}. \end{aligned}$$
(8)

If we restrict to asymptotically flat spacetimes, \(\omega _0\) can be interpreted as the photon frequency at infinity, b can be named as an impact parameter, while \(R_\mathrm {sh}\) will be identical to the impact parameter of the photon sphere.

Equations (5) and (7) are key formulae for us to compute radii of the photon sphere and the shadow, while Eqs. (4) and (8) are crucial for tracing the trajectory of light rays. We will assume \(h(r)\ge 0\) and \(0\le \alpha \le \pi /2\) for practical reasons. To trace the trajectory of light rays, one should integrate the orbit equation (4), usually using a numerical method. For the convenience of numerical integration, we will introduce a new variable \(u=r_g^{1/2}/r^{1/2}\).

Throughout this paper, we work in Gaussian units. The mass of the black hole will be denoted as M, and its mass accretion rate will be denoted as \({\dot{M}}_A\). For brevity, we introduce the gravitational radius \(r_g=2GM/c^2\) with G being the Newtonian constant and two dimensionless parameters

$$\begin{aligned} \beta&=\frac{e^2{\dot{M}}_A}{m_em_p\omega _0^2cr_g^2}\approx 1.0\times 10^{-7}\nonumber \\&\quad \times \left( \frac{\lambda _0}{0.1\,\mathrm {cm}}\right) ^2\left( \frac{10^9M_{\odot }}{M}\right) ^2\left( \frac{{\dot{M}}_A}{0.1M_{\odot }\mathrm {yr}^{-1}}\right) , \end{aligned}$$
(9)
$$\begin{aligned} \gamma&=\frac{4G{\dot{M}}_A}{3c^3}\approx 2.1\times 10^{-14}\times \left( \frac{{\dot{M}}_A}{0.1M_{\odot }\mathrm {yr}^{-1}}\right) . \end{aligned}$$
(10)

The wavelength of the photon is related to its frequency by \(\lambda _0=2\pi c/\omega _0\).

3 Toy models

It is well known that the metric of the Schwarzschild black hole

$$\begin{aligned} ds^2= & {} -\left( 1-\frac{r_g}{r}\right) c^{2}dt^{2}\nonumber \\&+\left( 1-\frac{r_g}{r}\right) ^{-1}dr^{2}+r^2\left( d\theta ^{2}+\sin ^{2}\theta d\phi ^{2}\right) \end{aligned}$$
(11)

is a static spherically symmetric solution of the vacuum Einstein equations. In Ref. [10], the refractive effect of a cold plasma on the shadow of Schwarzschild black holes has been evaluated in a concrete model [14], where the plasma has a mass density

$$\begin{aligned} \rho (r)=\left\{ \begin{array}{ll} 0, &{}\quad \mathrm {for}~r\le r_c,\\ \frac{{\dot{M}}_A}{4\pi cr_g^{1/2}r^{3/2}}, &{}\quad \mathrm {for}~r>r_c. \end{array}\right. \end{aligned}$$
(12)

Here we have introduced a critical or cutoff radius \(r_c\) for clarity. For a fully ionized hydrogen plasma, the plasma mass density is related to the electron number density by \(\rho (r)=m_pN(r)\).

The gravitational field of the plasma has been neglected in Ref. [10]. If we take the plasma’s gravity into consideration, the Schwarzschild metric will be modified into a time-dependent metric. The reasons are as follows. First, owing to the accretion, the mass or equivalently radius of the black hole must increase with time. Second, even if we ignore the accretion, the plasma cannot keep static against the gravitational attraction of the central black hole. Consequently, the modified metric ought to be time-dependent both inside and outside the event horizon. This will prevent us from applying the formulae shown in Sect. 2. In order to circumvent this obstacle, it is customary to neglect the time dependence and assume [15, 16]

$$\begin{aligned} A(r)=\frac{1}{B(r)},~~~~D(r)=r^2 \end{aligned}$$
(13)

in Eq. (1). Making use of one of the Einstein equations, one can verify that

$$\begin{aligned} \frac{1}{B(r)}&=1-\frac{2Gm(r)}{c^2r}, \end{aligned}$$
(14)
$$\begin{aligned} m(r)&=\left\{ \begin{array}{ll} M, &{}\quad \mathrm {for}~r\le r_c,\\ M+4\pi \int _{r_c}^{r}\rho (r')r'^2dr', &{}\quad \mathrm {for}~r>r_c. \end{array}\right. \end{aligned}$$
(15)

Working out the integral with the mass density (12), we find

$$\begin{aligned} m(r)=\left\{ \begin{array}{ll} M, &{}\quad \mathrm {for}~r\le r_c,\\ M+\frac{2{\dot{M}}_A}{3cr_g^{1/2}}\left( r^{3/2}-r_c^{3/2}\right) , &{}\quad \mathrm {for}~r>r_c. \end{array}\right. \end{aligned}$$
(16)

Outside the sphere of critical radius, \(r>r_c\), the above metric coincides with the case of \(w=-0.5\) in Ref. [17]. However, one should be cautious that the physical interpretations are completely different [18, 19], and here we are studying a different region of parameter space to estimate the gravitational effect of plasma particles.

In following sections, we will focus on two special models. All of the above equations are assumed to hold for both models. Their difference is merely the value of \(r_c\), which is specified to \(r_g\) in model A and to \(3r_g\) in model B. Let us study them separately in Sects. 4 and 5.

4 Model A: \(r_c=r_g\)

In Ref. [10], it was implicitly assumed that the plasma density does not vanish on the photon sphere, which implies \(r_c<3r_g/2\). In model A, we set \(r_c=r_g\) to meet this condition concretely. Then Eq. (16) becomes

$$\begin{aligned} m(r)=\left\{ \begin{array}{ll} M, &{}\quad \mathrm {for}~r\le r_g,\\ M+\frac{2{\dot{M}}_A}{3cr_g^{1/2}}\left( r^{3/2}-r_g^{3/2}\right) , &{}\quad \mathrm {for}~r>r_g. \end{array}\right. \end{aligned}$$
(17)

Corresponding to (17), the spacetime has two horizons. One is the event horizon of black hole located at \(r=r_g\), the other is the cosmological horizon at

$$\begin{aligned} r=\frac{r_g(1-\gamma )}{4\gamma ^2}\left( \sqrt{1-\gamma }+\sqrt{1+3\gamma }\right) ^2 \end{aligned}$$
(18)

with a small parameter \(\gamma \) defined in Eq. (10). Both the photon sphere and the static observer live in the patch of spacetime between the two horizons.

4.1 Influences on the radius

With the model specified, we are now ready to compute the plasma’s refractive and gravitational effects on the shadow accurately. Before entering on computations, let us make a rough comparison of the two effects. The refractive effect is induced intrinsically by the electromagnetic interaction. It is textbook knowledge that in a hydrogen atom, the ratio of the gravitational force to the electromagnetic force is \(Gm_em_p/e^2=4.4\times 10^{-40}\). Therefore, one would naively expect that the gravitational effect is suppressed by this factor in comparison with the refractive effect.

Let us take a closer look at model A. For this model, we can express the function \(h(r)^2\) in terms of \(\beta \) and \(\gamma \) explicitly,

$$\begin{aligned} h(r)^2= & {} r^2\left[ 1-\frac{r_g}{r}-\frac{\gamma }{r_g^{1/2}r}\left( r^{3/2}-r_g^{3/2}\right) \right] ^{-1}\nonumber \\&-r^2\beta \left( \frac{r_g}{r}\right) ^{3/2}. \end{aligned}$$
(19)

According to Eqs. (5) and (7), radii of the photon sphere and the shadow are determined simply by this function. It is clear that the refractive effect is controlled by \(\beta \), while the gravitational effect is dictated by \(\gamma \). For instance, one can switch off the gravitational effect by setting \(\gamma =0\), then the above expression will reduce to Eq. (50) in Ref. [10]. From this point of view, one would expect that the gravitational effect is suppressed by a factor \(\gamma /\beta \sim Gm_em_pr_g^2/(e^2\lambda _0^2)\) in comparison with the refractive effect. Intriguingly, here is an enhancement factor \(r_g^2/\lambda _0^2\) in competition with \(Gm_em_p/e^2\).

As a further step, one can substitute Eq. (19) into Eq. (5) and search for the largest root, but the resulted equation cannot be solved exactly. To get some analytical results, we solve it perturbatively for small \(\beta \) and \(\gamma \). In this way, we find the second-order solution

$$\begin{aligned} r_\mathrm {ph}\approx r_{00}+r_{10}-r_{01}+r_{20}-r_{11}-r_{02}, \end{aligned}$$
(20)

where

$$\begin{aligned} r_{00}&=\frac{3}{2}r_g,~~~~r_{10}=\frac{\sqrt{6}}{108}\beta r_g,~~~~r_{01}=\left( \frac{3}{2}-\frac{9\sqrt{6}}{16}\right) \gamma r_g, \end{aligned}$$
(21)
$$\begin{aligned} r_{20}&=\frac{7}{5832}\beta ^2r_g,~~~~r_{11}=\left( \frac{1}{16}-\frac{\sqrt{6}}{216}\right) \beta \gamma r_g,~~~~r_{02}\nonumber \\&=\left( \frac{27\sqrt{6}}{32}-\frac{243}{128}\right) \gamma ^2r_g. \end{aligned}$$
(22)

With the help of Eq. (7), we can proceed to compute the radius of shadow observed at a distance \(r_\mathrm {o}\). To the first order, the result is

$$\begin{aligned} R_\mathrm {sh}\approx R_{00}-R_{10}-R_{01} \end{aligned}$$
(23)

in which

$$\begin{aligned} R_{00}&=\frac{3\sqrt{3}}{2}r_g\left( 1-\frac{r_g}{r_\mathrm {o}}\right) ^{1/2}, \end{aligned}$$
(24)
$$\begin{aligned} R_{10}&=\left[ \frac{\sqrt{2}}{6}-\frac{3\sqrt{3}}{4}\left( \frac{r_g}{r_\mathrm {o}}\right) ^{3/2}\left( 1-\frac{r_g}{r_\mathrm {o}}\right) \right] \nonumber \\&\quad \times \beta r_g\left( 1-\frac{r_g}{r_\mathrm {o}}\right) ^{1/2}, \end{aligned}$$
(25)
$$\begin{aligned} R_{01}&=\left[ \frac{3\sqrt{3}}{4}\left( \frac{r_\mathrm {o}}{r_g}\right) ^{1/2}+\frac{3}{8}\left( 4\sqrt{3}-9\sqrt{2}\right) \right. \nonumber \\&\quad \left. +\frac{3\sqrt{3}r_g^{1/2}}{4\left( r_\mathrm {o}^{1/2}+r_g^{1/2}\right) }\right] \gamma r_g\left( 1-\frac{r_g}{r_\mathrm {o}}\right) ^{1/2}. \end{aligned}$$
(26)

It is straightforward to derive higher order terms. They are very lengthy but negligibly small, and thus not shown here. Remarkably, the gravitational correction to the radius of the black hole shadow is enhanced further by \((r_\mathrm {o}/r_g)^{1/2}\) in Eq. (26). In summary, the ratios of the gravitational corrections to the refractive corrections are

$$\begin{aligned} \frac{r_{01}}{r_{10}}&\approx 36\pi ^2\left( 4\sqrt{6}-9\right) \times \frac{Gm_em_pr_g^2}{e^2\lambda _0^2}\nonumber \\&\approx 1.1\times 10^{-6}\times \left( \frac{0.1\,\mathrm {cm}}{\lambda _0}\right) ^2\left( \frac{M}{10^9M_{\odot }}\right) ^2, \end{aligned}$$
(27)
$$\begin{aligned} \frac{R_{01}}{R_{10}}&\approx 12\sqrt{6}\pi ^2\times \frac{Gm_em_pr_g^{3/2}r_\mathrm {o}^{1/2}}{e^2\lambda _0^2}\nonumber \\&\approx 0.4\times \left( \frac{0.1\,\mathrm {cm}}{\lambda _0}\right) ^2\left( \frac{M}{10^9M_{\odot }}\right) ^{3/2}\left( \frac{r_\mathrm {o}}{10\mathrm {Mpc}}\right) ^{1/2}. \end{aligned}$$
(28)
Table 1 Masses [11], distances and mass accretion rates [12] for a sample of AGNs
Fig. 1
figure 1

Influences of the plasma on the radius of the black hole shadow for model A in Sect. 4. Left panel: the radius of the photon sphere and corrections to it according to Eq. (21) and Table 1 with \(\lambda _0=0.13\,\mathrm {cm}\). Right panel: the radius of the shadow and corrections to it according to Eqs. (24), (25), (26) and Table 1. In both panels, black asterisks represent the uncorrected radii, red diamonds denote the refractive corrections, and blue circles mark the gravitational corrections

To convert the above analytical expressions to numbers, we need the values of \(\lambda _0\), \(r_\mathrm {o}\), M and \({\dot{M}}_A\). Table 1 is a sample of AGNs selected from Refs. [11, 12]. The data of masses are extracted from Ref. [11], while the data of distances and mass accretion rates are obtained from Ref. [12]. Making use of the data and ignoring the spins of black holes, we evaluate Eq. (21) for each AGN and draw the results in the left panel of Fig. 1 in logarithmic coordinates, and Eqs. (24), (25), (26) in the right panel. In both panels, the uncorrected radii are drawn as asterisks, the refractive corrections are depicted by diamonds, and the gravitational corrections are denoted by circles. For all of them, we have set \(\lambda _0=0.13\,\mathrm {cm}\) [3]. This figure can be regarded as a semi-realistic estimation of refractive and gravitational effects of plasma on the radii of the photon sphere and the shadow. We can see clearly that these effects are negligible despite of the diversity of AGNs in the sample. The right panel is in good agreement with Eq. (28), which shows that the two effects on the radius of shadow are comparable for black holes of mass \(M\sim 10^9M_{\odot }\). For Sgr A*, the gravitational corrections are smallest because of its slowest accretion rate.

4.2 Influences on the intensity

From the conclusion of the previous subsection, one can infer that the refractive and gravitational effects on the intensity image should be undetectable for realistic AGNs. This is due to the smallness of \(\beta \) and \(\gamma \) in Eqs. (9), (10). In the current subsection, to make the effects noticeable, we will assign unrealistically large values to these parameters.

For the toy model A, the orbit equation (4) of light rays can be reformed as

$$\begin{aligned} \frac{r_g}{D(r)}\frac{dr}{d\phi }=\pm \sqrt{\frac{r_g^2}{b^2}\left[ 1-A\frac{\omega _p(r)^2}{\omega _0^2}\right] -\frac{r_g^2A}{D(r)}}. \end{aligned}$$
(29)

In terms of a new variable \(u=r_g^{1/2}/r^{1/2}\), we can write down \(\omega _p(r)^2/\omega _0^2=\beta u^3\), \(D(r)=r_g^2u^{-4}\) and

$$\begin{aligned} A=u^{-1}\left[ u-(1-\gamma )u^3-\gamma \right] \end{aligned}$$
(30)

outside the event horizon. For clockwise light rays \(du/d\phi >0\), the orbit equation can be reexpressed in terms of u as

$$\begin{aligned} 2u\frac{du}{d\phi }=\sqrt{\frac{r_g^2}{b^2}\left( 1-\beta u^3A\right) -u^4A}. \end{aligned}$$
(31)

When performing the numerical integration, exaggerated values are assigned to the model parameters as labeled in Fig. 2, \(r_\mathrm {o}=10r_g\), \(\beta =0\) or 0.5, \(\gamma =0\) or 0.05.

Let us consider an observer located outside the photon sphere. According to the ratio \(b/h(r_\mathrm {ph})\), the light rays arriving at this observer can be classified into three types [20] which are depicted by different line-types in Fig. 2: (i) Light rays with \(b<h(r_\mathrm {ph})\) can travel through the photon sphere. As illustrated by green solid curves in Fig. 2, all orbits of these rays start near the event horizon of the black hole. (ii) A critical light ray has \(b=h(r_\mathrm {ph})\). It propagates in an unstable circular orbit of radius \(r_\mathrm {ph}\) and has a chance to escape to the observer under radial perturbations. The circle is a great circle of the photon sphere. In every left panel of Fig. 2, we plot such a critical light ray as a red dotted curve. (iii) Light rays with \(b>h(r_\mathrm {ph})\) are always outside the photon sphere. As depicted by blue dashed curves in Fig. 2, each orbit of such rays is symmetric with respect to a diametrical line through its pericenter. The radial coordinate \(r_\mathrm {min}\) of the pericenter is dictated by \(h(r_\mathrm {min})=b\).

Fig. 2
figure 2

Influences of the plasma on the intensity of the black hole shadow for model A in Sect. 4. Left panels: trajectories of the light rays observed at the point \((x,y)=(10r_g,0)\). Green solid curves, red dotted curves and blue dashed curves depict trajectories with \(b<h(r_\mathrm {ph})\), \(b=h(r_\mathrm {ph})\) and \(b>h(r_\mathrm {ph})\) respectively. The interval of adjacent orbits is \(\Delta b/r_g=0.1\). From top left to bottom left, the critical trajectories \(b=h(r_\mathrm {ph})\) correspond to \(r_\mathrm {ph}/r_g=1.5,1.511647,1.493443\) and \(h(r_\mathrm {ph})/r_g=2.598076,2.477186,2.714083\) in turn. Right panels: the observed intensity image for \(r_\mathrm {o}=10r_g\), \(\beta =0\), \(\gamma =0\) (top right) and its differences with images for \(\beta =0.5\), \(\gamma =0\) (middle right) and \(\beta =0\), \(\gamma =0.05\) (bottom right). All Axes are rescaled by \(r_g\)

It is reasonable to expect that the radiant energy density is proportional to the plasma density. This implies a specific emissivity in the rest-frame of the emitter

$$\begin{aligned} j(\nu ,r)\propto \frac{1}{r^{3/2}}. \end{aligned}$$
(32)

For simplicity we assume the emission is uniform in frequency from a static source. Then \(r^{3/2}j(\nu ,r)\) is independent of \(\nu \) and r, and thus the specific intensity at the point (xy) of the observed image is [17, 21, 22]

$$\begin{aligned} I_{\beta ,\gamma }=\frac{1}{2}{\mathcal {I}}r_g^{1/2}\int _{\mathrm {ray}}\frac{A(r)^{3/2}}{r^{3/2}}\sqrt{B(r)+r^2\left( \frac{d\phi }{dr}\right) ^2}dr. \end{aligned}$$
(33)

Here the normalization \({\mathcal {I}}\) is unimportant in our simulations, and the subscripts are used to reminding us the dependence of intensity on \(\beta \) and \(\gamma \). Numerically it is more convenient to perform this integration in terms of u defined above,

$$\begin{aligned} I_{\beta ,\gamma }={\mathcal {I}}\int _{\mathrm {ray}}A^{3/2}\sqrt{\frac{1}{A}+\frac{u^2}{4}\left( \frac{d\phi }{du}\right) ^2}du. \end{aligned}$$
(34)

The integration is performed using the backward ray shooting method [20, 21]. For type (i) light rays, we integrate Eq. (34) from the observer to the event horizon of the black hole. For type (iii) light rays, the integration is performed from the observer to the pericenter and then to the cosmological horizon (18). In practice, telescopes do not collect light rays with \(\pi /2<\alpha \le \pi \), thus we assume \(0\le \alpha \le \pi /2\) in our simulations.

By virtue of the spherical symmetry of the toy model, the image is rotationally invariant, and thus the intensity is dependent simply on \(R=\left( x^2+y^2\right) ^{1/2}\). By definition Eq. (8), it is easy to show

$$\begin{aligned} R^2=\frac{b^2A}{1-\beta u^3A} \end{aligned}$$
(35)

with \(u=r_g^{1/2}/r_\mathrm {o}^{1/2}\). By continuously changing the value of b, we have computed the dependence of intensity on R and simulated the images with \(r_\mathrm {o}=10r_g\), \(\beta =0\) or 0.5, \(\gamma =0\) or 0.05. The slight differences between these images are not easy to notice, especially near the ring. In the top right panel of Fig. 2, we present the image for \(I_{0,0}\). In the middle right and bottom right panels, we subtract \(I_{0,0}\) from \(I_{0.5,0}\) and \(I_{0,0.05}\) to illustrate the slight differences. In comparison with the case of \(\beta =0,\gamma =0\), the intensity outside the ring decreases more significantly in the case of \(\beta =0,\gamma =0.05\) than the case of \(\beta =0.5,\gamma =0\), but the changes in radius are about the same.

5 Model B: \(r_c=3r_g\)

In Ref. [14], dropping corrections to the metric from the external plasma, it was shown that no physically acceptable solution exists if \(r_c<3r_g\). In other words, the critical radius ought to be slightly larger than the radius of the innermost stable circular orbit. In model B, we consider the limiting case and set \(r_c=3r_g\). Then Eq. (16) becomes

$$\begin{aligned} m(r)=\left\{ \begin{array}{ll} M, &{}\quad \mathrm {for}~r\le 3r_g,\\ M+\frac{2{\dot{M}}_A}{3cr_g^{1/2}}\left( r^{3/2}-3\sqrt{3}r_g^{3/2}\right) , &{}\quad \mathrm {for}~r>3r_g. \end{array}\right. \nonumber \\ \end{aligned}$$
(36)

Corresponding to (36), the event horizon of black hole is located at \(r=r_g\), and the cosmological horizon is located at

$$\begin{aligned} r= & {} \frac{r_g}{9\gamma ^2}\bigg [1+\left( C+\sqrt{C^2-1}\right) ^{1/3}\nonumber \\&+\left( C-\sqrt{C^2-1}\right) ^{1/3}\bigg ]^2 \end{aligned}$$
(37)

with

$$\begin{aligned} C=\frac{1}{2}\left( 2-27\gamma ^2+81\sqrt{3}\gamma ^3\right) . \end{aligned}$$
(38)

Both the photon sphere and the static observer live in the patch of spacetime between the two horizons.

5.1 Influences on the radius

In model A, the density profile of the plasma is a continuous function outside the event horizon. In model B, however, the density profile has a step at the innermost stable circular orbit. Thereby the function \(h(r)^2\) is discontinuous, taking the form

$$\begin{aligned} h(r)^2=\left\{ \begin{array}{ll} r^2\left( 1-\frac{r_g}{r}\right) ^{-1}, &{}\quad \mathrm {for}~r\le 3r_g,\\ r^2\left[ 1-\frac{r_g}{r}-\frac{\gamma }{r_g^{1/2}r}\left( r^{3/2}-3\sqrt{3}r_g^{3/2}\right) \right] ^{-1}-r^2\beta \left( \frac{r_g}{r}\right) ^{3/2}, &{} \quad \mathrm {for}~r>3r_g. \end{array}\right. \end{aligned}$$
(39)

Applying Eq. (5) to this function, we find the radius of the photon sphere can be derived from the function in the region \(r\le 3r_g\) as long as \(\beta \) and \(\gamma \) are small. The result is

$$\begin{aligned} r_\mathrm {ph}=\frac{3}{2}r_g, \end{aligned}$$
(40)

the same as the Schwarzschild black hole without a plasma. Indeed, Eq. (5) indicates that \(r_\mathrm {ph}\) is determined locally by \(h(r)^2\) near the photon sphere as long as the equation has no root outside the sphere. Reversing this logic, one can infer that the plasma density does not vanish on the photon sphere in Ref. [10], because \(r_\mathrm {ph}^1\ne 0\) in Eq. (53) of Ref. [10].

Although the value of \(r_\mathrm {ph}\) is trivial in model B, the value of \(R_\mathrm {sh}\) is nontrivial. According to Eq. (7), the radius of the shadow depends on both the radius of the photon sphere and the radial coordinate of the observer. As an analytical result, the radius of shadow observed at a distance \(r_\mathrm {o}\) is

$$\begin{aligned} R_\mathrm {sh}\approx R_{00}+R_{10}-R_{01} \end{aligned}$$
(41)

to the first order of \(\beta \) and \(\gamma \), in which

$$\begin{aligned}&R_{00}=\frac{3\sqrt{3}}{2}r_g\left( 1-\frac{r_g}{r_\mathrm {o}}\right) ^{1/2}, \end{aligned}$$
(42)
$$\begin{aligned}&R_{10}=\frac{3\sqrt{3}}{4}\left( \frac{r_g}{r_\mathrm {o}}\right) ^{3/2}\left( 1-\frac{r_g}{r_\mathrm {o}}\right) ^{3/2}\beta r_g, \end{aligned}$$
(43)
$$\begin{aligned}&R_{01}=\left[ \frac{3\sqrt{3}}{4}\left( \frac{r_\mathrm {o}}{r_g}\right) ^{1/2}+\frac{3\sqrt{3}r_g^{1/2}}{4\left( r_\mathrm {o}^{1/2}+r_g^{1/2}\right) }\right] \nonumber \\&\qquad \quad \gamma r_g\left( 1-\frac{r_g}{r_\mathrm {o}}\right) ^{1/2}. \end{aligned}$$
(44)

In this model, the ratio of the gravitational correction to the refractive correction is

$$\begin{aligned} \frac{R_{01}}{R_{10}}\approx & {} \frac{16\pi ^2}{3}\times \frac{Gm_em_pr_\mathrm {o}^2}{e^2\lambda _0^2}\nonumber \\\approx & {} 2\times 10^{15}\times \left( \frac{0.1\,\mathrm {cm}}{\lambda _0}\right) ^2\left( \frac{r_\mathrm {o}}{10\mathrm {Mpc}}\right) ^2. \end{aligned}$$
(45)

Making use of the data in Sect. 4.1, we evaluate Eqs. (42), (43), (44) for each AGN and draw the results in Fig. 3.

Fig. 3
figure 3

Influences of the plasma on the radius of the black hole shadow for model B in Sect. 5. Black asterisks represent the uncorrected radius according to Eq. (42), red diamonds denote the refractive corrections according to Eq. (43), and blue circles mark the gravitational corrections according to Eq. (44). Data in Table 1 have been used here

5.2 Influences on the intensity

In terms of \(u=r_g^{1/2}/r^{1/2}\), the orbit equation for clockwise light rays can be expressed as

$$\begin{aligned} 2u\frac{du}{d\phi }=\left\{ \begin{array}{ll} \sqrt{\frac{r_g^2}{b^2}\left( 1-\beta u^3A\right) -u^4A}, &{}\quad \mathrm {for}~u<1/\sqrt{3},\\ \sqrt{\frac{r_g^2}{b^2}-u^4\left( 1-u^2\right) }, &{}\quad \mathrm {for}~1/\sqrt{3}\le u\le 1 \end{array}\right. \nonumber \\ \end{aligned}$$
(46)

in model B, where

$$\begin{aligned} A=u^{-1}\left[ u-(1-3\sqrt{3}\gamma )u^3-\gamma \right] . \end{aligned}$$
(47)

Considering an observer located at \(r_\mathrm {o}=10r_g\), we illustrate the trajectories of light rays in left panels of Fig. 4 for \(\beta =0\) or 0.5, \(\gamma =0\) or 0.05. The innermost stable circular orbit is depicted by a big black circle. Noticeably, when \(\beta \ne 0\), the behaviors of light rays in the neighborhoods of the innermost stable circular orbit are distinct. This can be understood from

$$\begin{aligned}&\lim _{u\rightarrow 1/\sqrt{3}}\sqrt{\frac{r_g^2}{b^2}\left( 1-\beta u^3A\right) -u^4A}\nonumber \\&\quad =\sqrt{\frac{r_g^2}{b^2}\left( 1-\frac{2\beta }{9\sqrt{3}}\right) -\frac{2}{27}}, \end{aligned}$$
(48)

which suggests that the bending of light is weakened by the refractive effect near the innermost stable circular orbit.

In model B, the plasma is the same as in model A when \(r>3r_g\) but vanishes when \(r\le 3r_g\). Therefore, the specific emissivity has the form

$$\begin{aligned} j(\nu ,r)\propto \frac{1}{r^{3/2}}H(r-3r_g). \end{aligned}$$
(49)

Here the Heaviside step function is defined as

$$\begin{aligned} H(x)=\left\{ \begin{array}{ll} 0, &{}\quad \mathrm {for}~x\le 0,\\ 1, &{}\quad \mathrm {for}~x>0. \end{array}\right. \end{aligned}$$
(50)

When computing the specific intensity of the observed image in this model, we should insert the Heaviside unit step function into Eq. (34),

$$\begin{aligned} I_{\beta ,\gamma }={\mathcal {I}}\int _{\mathrm {ray}}A^{3/2}\sqrt{\frac{1}{A}+\frac{u^2}{4}\left( \frac{d\phi }{du}\right) ^2}H\left( \frac{1}{\sqrt{3}}-u\right) du\nonumber \\ \end{aligned}$$
(51)

and take A as Eq. (47). The images for model B are presented in right panels of Fig. 4. In this model, the image of photon sphere makes negligible contributions to the total brightness. Similar results have been reported previously in Refs. [17, 23, 24].

Fig. 4
figure 4

Influences of the plasma on the intensity of the black hole shadow for model B in Sect. 5. Left panels: trajectories of the light rays observed at the point \((x,y)=(10r_g,0)\). Green solid curves, red dotted curves and blue dashed curves depict trajectories with \(b<h(r_\mathrm {ph})\), \(b=h(r_\mathrm {ph})\) and \(b>h(r_\mathrm {ph})\) respectively. The interval of adjacent orbits is \(\Delta b/r_g=0.1\). Right panels: the observed intensity image for \(r_\mathrm {o}=10r_g\), \(\beta =0\), \(\gamma =0\) (top right) and its differences with images for \(\beta =0.5\), \(\gamma =0\) (middle right) and \(\beta =0\), \(\gamma =0.05\) (bottom right). All Axes are rescaled by \(r_g\)

6 Discussion

In this paper, we have investigated the gravitational effect of plasma particles on the radius and the intensity of the black hole shadow. We did this in a toy model of a Schwarzschild black hole surrounded by a static plasma of density \(\rho (r)\propto r^{-3/2}\) outside the event horizon. Inserting observational data for a sample of AGNs into our analytical results, we find the gravitational correction \(R_{01}\) to the shadow radius is in the range \(10^{-8}\sim 10^{-6}\). In contrast, we find the correction \(R_{10}\) from the refractive effect of plasma medium decreases prominently as the black hole mass increases, covering the range \(10^{-4}\sim 10^{-10}\). Alternatively, turning to the scenario that the space is empty between event horizon and innermost stable circular orbit, we find the shadow is not determined by the photon sphere, but by the inner boundary of the accretion disk.

Strictly speaking, our analytical results are reliable only for the static toy model [25], so the numbers above are very rough estimations for realistic AGNs which usually have nonnegligible spins [26, 27]. However, as a concrete example, they demonstrate that the gravitational effect of plasma particles can exceed the refractive effect of the plasma medium under some circumstances. The gravitational effect is more complicated to compute than the refractive effect. If the refractive effect is proven to be detectable in a more realistic model, then the gravitational effect should be scrutinized, especially for AGNs of high masses.

We did not consider synchrotron self-absorption [5], which is desirable when studying the background emission of the plasma. It would also be interesting to take the black hole spin [28,29,30,31,32,33,34,35] and the QED effect [36] into account.