1 Introduction

The images of M87* and Sgr A* black holes photographed by the Event Horizon Telescope (EHT) [1, 2] are direct evidences for supporting the existence of black holes in the universe. The astronomical observations including the images reveal strong magnetic fields near the event horizons of M87* and Sgr A* [3, 4]. The generation of such external magnetic fields is closely related to the jets and energy extraction and magnetic accretion disks.

When a magnetic field in the vicinity of a black hole with mass m has its strength much smaller than the order \(B\sim B_{max}\sim 10^{19}m_{\odot }/m\) [5], it does not change the spacetime geometry. It has a negligible impact on the motion of neutral test particles and photons around the black hole. However, it can drastically affect the motion of charged test particles subjected to the large Lorentz force. A possible dramatic effect is that the external magnetic field causes the dynamics of charged particles to be nonintegrable. In fact, the chaotic motions of charged particles around one of the Schwarzschild, Reissner–Nordström (RN), Kerr and Kerr–Newman black holes immersed in uniform external magnetic fields were shown in Refs. [6,7,8,9,10,11,12,13,14,15,16,17,18,19,20]. However, the motions of neutral particles in these spacetimes should be integrable and regular if the external magnetic fields are absent. The chaotic or regular orbits are also allowed to the timelike bound orbits, which neither fall into the black holes nor escape to infinity. The authors of [21] showed that no negative precession of timelike bound orbits but positive precession of timelike bound orbits is allowed for the S-stars in Schwarzschild or Kerr spacetimes. The existence of these bound orbits is owing to the effective potentials of massive particles on the equatorial plane having closed pockets or potential wells with local minimum values. Nevertheless, the equatorial photon effective potentials have no closed pockets but potential barriers with local maximum values outside the horizons in the four spacetimes. Thus, bound photon orbits involving stable circular photon orbits are not allowed in general. An exceptional case is the presence of stable equatorial circular photon orbits within the inner horizon in Kerr–Newman spacetimes [22]. Instead, the motions of photons around the black holes have three possible regions for the photons falling to the center, scattering to infinity and unstably circling in the center [23]. The three types of photon orbits are not necessarily restricted to the equatorial plane. The dark region for the light rays falling into the event horizon of a black hole is a shadow of the black hole, seen by an observer inside a celestial sphere. All the unstable circular or spherical photon orbits correspond to a region in the celestial sphere as the black hole shadow boundary. In this sense, the black hole shadows in these spacetimes can be calculated through an analytical method [23].

If an external magnetic field is close to the upper limit of magnetic field with the order \(B_{max}\), such a strong magnetic field would significantly distort the spacetime geometry near a black hole [5]. It makes the black hole strongly magnetized and is included in the spacetime geometry. Ernst [24] gave magnetized Schwarzschild, RN, Kerr, and Kerr–Newman black hole spacetimes to stationary and axisymmetric solutions of the Einstein–Maxwell equations, where the Schwarzschild, RN, Kerr, and Kerr–Newman black holes are nonlinearly coupled with the Melvin’s magnetic universe. A full magnetized Kerr–Newman solution with an electric charge and a magnetic charge was derived in [25]. This magnetized Kerr–Newman spacetime is not asymptotic to the static Melvin metric unless the total electric charge of the black hole satisfies the condition \(Q=jB(1+ j^2B^4/4)\) or a purely electric charge of the original seed black hole is \(q=-amB\). Here, \(j=ma\) is the angular momentum of the original seed solution and a is the spin parameter of the black hole. In addition, m is the mass of the black hole and B denotes the magnetic field strength. The external magnetic fields seeded in the Schwarzschild, RN, Kerr, and Kerr–Newman black holes make these magnetized spacetimes be similar to AdS spacetimes that are not asymptotically flat. They act as gravitational effects in the magnetized spacetimes and strongly affect the radii of the innermost stable circular orbits for neutral test particles around the Schwarzschild and RN black holes [26, 27]. They lead to the nonintegrability of these spacetimes, and therefore they can even induce the occurrence of chaos of neutral particles under some circumstances, as was reported in Refs. [28,29,30,31]. Of course they can also strongly influence the dynamics of charged particles around these magnetized black holes. These chaotic orbits of both neutral and charged particles should be bound orbits.

On the other hand, the motions of photons in the vicinity of these magnetized spacetimes can contain imprints of a chaotic pattern. Using the backward ray-tracing method, Junior et al. [32] observed the shadows of Schwarzschild–Melvin black hole possessing self-similar fractal structures originating from the chaotic lensing. One scattered orbit in the chaotic region on the celestial sphere was also taken as an example to show the chaotic pattern. In addition, bound photon orbits can exist around the stable circular photon orbits on the equatorial plane, and can also exist outside the equatorial plane in the Schwarzschild–Melvin black hole spacetime. Wang et al. [33] claimed that there are equatorial stable circular photon orbits as well as equatorial unstable ones in the Kerr–Melvin spacetime. They further pointed out the existence of some self-similar fractal structures in the shadow of Kerr–Melvin black hole. Hou et al. [34] reported that chaos occurs in the Kerr–Melvin image with thin accretion disk. The chaotic features in the shadows and lensing of magnetized black holes arise from the chaotic motion of photons. Note that these magnetized spacetimes are not asymptotically flat, therefore, their shadows are not computed by choosing a reference frame with zero angular momentum at spatial infinity. An observer at finite distance should be considered. Recently, a new method was used to determine physical parameters of Kerr black holes by shadow observation at finite distance [35]. The method is applicable to various type of black holes such as the magnetized black holes [32,33,34] and superentropic AdS black holes [36]. Besides those in general relativity, the shadows of black holes in modified gravity theories (see e.g. [37]) are widely focused on. Some of the modified gravity theories can cope with the two drawbacks of the theory of general relativity on the presence of curvature singularities and the failure of explaining the existence of the dark sector. An example of the modified gravity theories is that the scalar–tensor–vector gravity theory has successfully explained the solar system observations, cosmological observations, galaxy rotation curves and the dynamics of galaxy clusters [37]. Another example is that renormalization group improved Schwarzschild black holes in quantum gravity eliminate the existence of curvature singularities by the gravitational constant modified as a position-dependent quality. Their parameters were constrained in the allowed ranges based on the shadow results of Sgr A* from EHT [38]. Some other modified gravity theory black hole spacetimes like Kerr black holes with scalar hair [39, 40] and a Hartle–Thorne spacetime including a quadrupole deformation from the Kerr spacetime [41] have chaotic lensing or self-similar structures in black hole shadows. RN diholes allow for the existence of bound regular and chaotic photon orbits outside the equatorial plane [22]. The presence of such bound photon orbits involving the stable circular photon orbits is due to the shape of the equatorial photon effective potential having closed pockets. Although the equatorial photon effective potential is simple, it plays an important role in the study of black hole shadows because it can provide some information on the unstable or stable circular photon orbits and classification of photon orbits.

A notable point is that the chaotic photon properties were shown through chaotic light scattering and the self-similar fractal structures of black hole shadows in Refs. [32,33,34, 39, 40, 42, 43]. No bound chaotic photon orbits but several equatorial stable photon orbits were given. A problem is that the common chaos detection methods like Lyapunov exponents [44, 45] and fast Lyapunov indicators (FLIs) [46, 47] are not very suitable for identifying the chaotic behavior of a scattered orbit because the scattered orbit may exponentially escape to infinity in the nonchaotic case. Another problem is that the fractal does not necessarily mean chaos [48]. In this case, it is vital to use more chaotic indicators to distinguish photon chaotic dynamics from photon regular dynamics. In fact, only a few authors focused on the dynamical structures of bound regular and chaotic photon orbits. Dolan and Shipley exhibited some bound regular and chaotic photon orbits on Poincaré sections in Majumdar–Papapetrou diholes [22]. Employing three tools consisting of Poincaré sections, power spectrum of time series and rotation number, Kostaros and Pappas studied the properties of bound photon orbits in the Hartle–Thorne spacetime [41].

Another notable point is that traditional numerical integration schemes such as Runge–Kutta integrators were used to construct the backward ray-tracing method in Refs. [32,33,34]. They have an artificial excitation or damping and cause a nonphysical increase in energy errors with time. As a result, the qualitative result of longtime evolution becomes quite unreliable. Symplectic integration schemes [49, 50] are particularly suitable for studying the qualitative long-term evolution of Hamiltonian dynamics because they preserve the symplectic structure of phase space of Hamiltonian dynamics and show no secular drift in errors of the integrals of motion in the Hamiltonian system. Explicit symplectic integrators are superior to implicit ones at same orders in computational efficiency. A main path for the construction of explicit symplectic integrators is based on splitting and composing and requires that the Hamiltonian should have the separation of variables or be split into two or more explicitly integrable terms. However, the Hamiltonian for a curved spacetime is not separable or has no two explicitly integrable splitting terms. In this case, the explicit symplectic integrators seem to be difficulty available for curved spacetimes. In fact, the Hamiltonians for some curved spacetimes such as the Schwarzschild black hole metric allow for the splitting of more than two explicitly integrable terms, and therefore allow for the application of explicit symplectic integrators [17,18,19]. Unfortunately, such a direct split does not exist for some other curved spacetimes, e.g. the Kerr black hole metric. Nevertheless, this problem can be addressed via appropriate time transformations to the Hamiltonians for these spacetimes and explicit symplectic integrators are still applicable [51, 52].

Noticing the above points, we shall investigate the motion of photons in the vicinity of the magnetized Kerr–Newman black hole, where the total electric charge of the black hole is \(Q=jB(1+ j^2B^4/4)\) or the purely electric charge of the original seed black hole is \(q=-amB\) [25]. We carry out our work along the equatorial photon effective potential due to its importance. The explicit symplectic integrators [17,18,19, 51, 52] combined with the techniques of Poincaré maps, Lyapunov exponents, FLIs and 0–1 test correlation method [53] were used to explore the dynamics of neutral and charged test particles near the magnetized Kerr–Newman black hole in [31]. Now, they are also adopted to give insight into the dynamics of photons around the magnetized Kerr–Newman black hole in the present paper. In this way, several types of photon orbits on and outside the equatorial plane will be found. In particular, many bound regular and chaotic photon orbits in the spacetime are given. These photon orbit regions are taken into account in the phase space rather than in the regions on the celestial sphere. In spite of this, the regions on the celestial sphere, which correspond to the photon orbits in the phase space, should be given in this similar way. That is, one of the explicit symplectic integrators combined with the fast Lyapunov indicators or the 0–1 test correlation method will be used as a new ray-tracing method to find the related photon regions on the observer’s plane in the future research of black hole shadows.

2 Photon effective potentials and circular orbits in a magnetized Kerr–Newman spacetime

A magnetized Kerr–Newman spacetime is introduced simply. Then photon effective potential and circular orbits are discussed mainly.

2.1 Magnetized Kerr–Newman spacetime

Consider that a rotating black hole has mass m, rotating angular momentum a and purely electric charge q. This rotating black hole is immersed in an external asymptotic uniform magnetic field with the strength B. In the Boyer–Lindquist coordinates \(x^{\alpha }=(t, r, \theta , \phi )\), Gibbons et al. [25] gave such a rotating black hole spacetime in the following metric

$$\begin{aligned} ds^2= & g_{\alpha \beta }dx^\alpha dx^\beta \nonumber \\= & g_{tt}dt^2+g_{rr}dr^2+g_{\theta \theta }d\theta ^2 \nonumber \\ & +g_{\phi \phi }d\phi ^2+2g_{t\phi } dt d\phi . \end{aligned}$$
(1)

The nonzero covariant metric components are

$$\begin{aligned} g_{tt}= & \frac{\omega ^2\Sigma \sin ^2\theta }{H\mathcal {R} ^2}-fH, ~~~~ g_{t\phi } = -\frac{\omega \Sigma \sin ^2\theta }{H\mathcal {R}^2}, \nonumber \\ g_{rr}= & \frac{H\mathcal {R}^2}{\Delta }, ~~~~~~~~~~~~~~~~~~ g_{\theta \theta } = H\mathcal {R}^2, \nonumber \\ g_{\phi \phi }= & \frac{ \Sigma \sin ^2\theta }{H\mathcal {R}^2}. \end{aligned}$$
(2)

The above notations are defined by

$$\begin{aligned} \mathcal {R}^2= & r^2+a^2\cos ^2\theta , \end{aligned}$$
(3)
$$\begin{aligned} \Delta= & (r^2+a^2)-2mr+q^2, \end{aligned}$$
(4)
$$\begin{aligned} f= & \frac{\mathcal {R}^2 \Delta }{\Sigma }, \end{aligned}$$
(5)
$$\begin{aligned} \Sigma= & (r^2+a^2)^2-a^2 \Delta \sin ^2\theta ,\end{aligned}$$
(6)
$$\begin{aligned} H= & 1+(H_{1}B+H_{2}B^2 +H_{3}B^3+H_{4}B^4)/\mathcal {R}^2. \end{aligned}$$
(7)

The expressions of functions \(H_{1}\), \(H_{2}\), \(H_{3}\) and \(H_{4}\) are

$$\begin{aligned} H_{1}= & 2 a q r \sin ^2 \theta , \end{aligned}$$
(8)
$$\begin{aligned} H_{2}= & \frac{1}{2}\left[ \left( r^2+a^2\right) ^2-a^2 \Delta \sin ^2 \theta \right] \sin ^2 \theta \nonumber \\ & +\frac{3}{2} q^2\left( a^2+r^2 \cos ^2 \theta \right) , \end{aligned}$$
(9)
$$\begin{aligned} H_{3}= & \frac{a q}{2 r}\left( r^2+a^2\right) ^2\left( 1+\cos ^2 \theta \right) \nonumber \\ & +\frac{aq^3}{2 r}\left[ \left( 2 r^2+a^2\right) \cos ^2 \theta +a^2\right] -\frac{aq\Delta }{2 r} \nonumber \\ & \times \left[ r^2\left( 3-\cos ^2 \theta \right) \cos ^2 \theta +a^2\left( 1+\cos ^2 \theta \right) \right] , \end{aligned}$$
(10)
$$\begin{aligned} H_{4}= & \frac{\mathcal {R}^2}{16}(r^2+a^2)^2 \sin ^4 \theta +\frac{mr}{4}a^2(r^2+a^2)\sin ^6\theta \nonumber \\ & +\frac{1}{4} m a^2 q^2 r(\cos ^2 \theta -5) \sin ^2 \theta \cos ^2 \theta \nonumber \\ & +\frac{m^2}{4} a^2\left[ r^2(\cos ^2 \theta -3)^2 \cos ^2 \theta +a^2(1+\cos ^2\theta )^2\right] \nonumber \\ & +\frac{q^2}{8} (r^2+a^2)(r^2+a^2+a^2 \sin ^2 \theta ) \sin ^2 \theta \cos ^2 \theta \nonumber \\ & +\frac{q^4}{16} \left[ r^2 \cos ^2 \theta +a^2\left( 1+\sin ^2\theta \right) ^2\right] \cos ^2 \theta . \end{aligned}$$
(11)

Function \(\omega \) reads as

$$\begin{aligned} \omega= & \frac{1}{\Sigma }[\left( 2 m r-q^2\right) a+\omega _{1} B] \nonumber \\ & +\frac{ B^2}{\Sigma }[\omega _{2}+\omega _{3} B+\omega _{4} B^2], \end{aligned}$$
(12)

where \(\omega _1\), \(\omega _2\), \(\omega _3\) and \(\omega _4\) are of the following functions

$$\begin{aligned} \omega _{1}= & -2 q r\left( r^2+a^2\right) , \end{aligned}$$
(13)
$$\begin{aligned} \omega _{2}= & -\frac{3}{2} a q^2\left( r^2+a^2+\Delta \cos ^2 \theta \right) , \end{aligned}$$
(14)
$$\begin{aligned} \omega _{3}= & 4 q m^2 a^2 r-2a^2 m q^3 \nonumber \\ & +\frac{1}{2} q r\left( r^2+a^2\right) \left[ r^2-a^2+\left( r^2+3 a^2\right) \cos ^2 \theta \right] \nonumber \\ & +\frac{1}{2} q^3 r\left[ \left( r^2+3 a^2\right) \cos ^2 \theta -2 a^2\right] \nonumber \\ & +q m\left[ r^4-a^4+r^2\left( r^2+3 a^2\right) \sin ^2 \theta \right] , \end{aligned}$$
(15)
$$\begin{aligned} \omega _{4}= & \frac{1}{2} a^3 m^3 r\left( 3+\cos ^4 \theta \right) -\frac{1}{16} a q^6 \cos ^4 \theta \nonumber \\ & -\frac{1}{8} a q^4 r^2\left( 2+\sin ^2 \theta \right) \cos ^2 \theta \nonumber \\ & -\frac{1}{8} a q^4a^2\left( 1+\cos ^4 \theta \right) \nonumber \\ & +\frac{1}{16} a q^2 r^2 \left( r^2+a^2\right) \left( 1-6 \cos ^2 \theta +3 \cos ^4 \theta \right) \nonumber \\ & -\frac{1}{16} a^3 q^2 \left( r^2+a^2\right) \left( 1+\cos ^4 \theta \right) \nonumber \\ & -\frac{1}{4} a^3 m^2 q^2\left( 3+\cos ^4 \theta \right) +\frac{1}{8} amq^4 r \cos ^4 \theta \nonumber \\ & +\frac{1}{4} a m^2\left[ r^4\left( 3-6 \cos ^2 \theta +\cos ^4 \theta \right) \right] \nonumber \\ & +\frac{1}{4} a m^2\left[ 2 a^2 r^2\left( 3 \sin ^2 \theta -2 \cos ^4 \theta \right) \right] \nonumber \\ & +\frac{1}{4} a m^2\left[ -a^4\left( 1+\cos ^4 \theta \right) \right] \nonumber \\ & +\frac{1}{4} a m q^2 r\left[ 2 r^2\left( 3-\cos ^2 \theta \right) \cos ^2 \theta \right] \nonumber \\ & +\frac{1}{4} a m q^2 r\left[ -a^2\left( 1-3 \cos ^2 \theta -2 \cos ^4 \theta \right) \right] \nonumber \\ & +\frac{1}{8} a m r^3 \left( r^2+a^2\right) \left( 3+6 \cos ^2 \theta -\cos ^4 \theta \right) \nonumber \\ & -\frac{1}{8} a^3 m r \left( r^2+a^2\right) (1-6 \cos ^2 \theta -3 \cos ^4 \theta ). \end{aligned}$$
(16)

Equations (2)–(16) are obtained from Eqs. (B.2), (B.4)–(B.6), (B.8) and (B.9) of Ref. [25], where magnetic charge \(p=0\) is considered. Here, the geometric unit systems \(C=G=1\) are given to the speed of light C and the gravitational constant G. When \(a^2+q^2<m^2\), the spacetime (1) has two event horizons

$$\begin{aligned} r^{\pm }_{eh}=m\pm \sqrt{m^2-(a^2+q^2)}. \end{aligned}$$
(17)

If \(a^2+q^2=m^2\), then the spacetime (1) has only one event horizon \(r_{eh}=m\).

For \(m=a=q=0\), the metric (1) is the Melvin magnetic universe as a static solution of the Einstein–Maxwell equations [54]. If \(a=q=0\), then the metric (1) reduces to the Schwarzschild–Melvin metric obtained from an exact solution of the Einstein–Maxwell equations [24]. However, the metric (1) for \(q=0\) corresponds to the Kerr–Melvin metric that does not exactly satisfies the Einstein–Maxwell equations in general. Only when the black hole takes the total electric charge Q as the Wald charge \(Q_W=2jB\) [55], i.e. \(Q=Q_w\), is such a Kerr–Melvin metric an exact solution of the Einstein–Maxwell equations [56]. Here, \(j=am\) represents the angular momentum of the original seed black hole. Similarly, the metric (1) contains an ergoregion extending to infinity and is not asymptotic to the Melvin solution for generic values of the mass m, charge q and rotation parameter a of the original seed Kerr–Newman solution, as claimed by the authors of Ref. [25]. These authors further pointed out that the metric can be asymptotic to the static Melvin metric when the total electric charge of the black hole satisfies \(Q=jB(1+j^2B^4/4)\), or the original seed Kerr–Newman black hole carries purely electric charge \(q=-amB\). In fact, the total electric charge consists of three parts: \(Q=q+Q_{bc}+Q_W\), where \(Q_{bc}=-q^3B^2/4\) is a charge resulting from the gravitational back-reaction of the external magnetic field on the geometry. If a rotating black hole has vanishing purely electric charge, then Q is the Wald charge generated by the rotating body as a rotating conductor immersed in the magnetic field [57]. When the horizon radius of a charged rotating black hole is much smaller than the Melvin radius 2/B, the back-reaction of the electromagnetic field is neglected, and then the Wald charge can be recovered [25].

The magnetized Kerr–Newman black hole solution (1) resembles the Melvin universe that is not asymptotically flat. However, the asymptotically flat case exists in the absence of a magnetic field. The magnetic field in the metric (1) leads to the highly curved geometry near infinity. Such a curved geometry near infinity also occurs in the case of AdS black holes.

For simplicity, dimensionless operations are given to the related variables and parameters via scale transformations: \(t\rightarrow mt\), \(r\rightarrow mr\), \(a\rightarrow ma\), \(q\rightarrow mq\), \(Q\rightarrow mQ\) and \(B\rightarrow B/m\). In this way, the black mass m in the above expressions is an eliminated factor; equivalently, it is also taken as one geometric unit \(m=1\). Because the magnetized Kerr–Newman black hole solution (1) is required to be asymptotic to the static Melvin metric, we let the black hole carry purely electric charge \(q=-aB\) throughout the present paper.

A photon moves around the magnetized Kerr–Newman black hole according to the Hamiltonian system

$$\begin{aligned} K = \frac{1}{2}g^{\alpha \beta }p_\alpha p_\beta , \end{aligned}$$
(18)

where the nonzero contravariant components of the metric are

$$\begin{aligned} g^{tt}= & \frac{g_{\phi \phi }}{g_{tt}g_{\phi \phi }-g^2_{t\phi }},\\ g^{\phi \phi }= & \frac{g_{tt}}{g_{tt}g_{\phi \phi }-g^2_{t\phi }}, \\ g^{t\phi }= & \frac{-g_{t\phi }}{g_{tt}g_{\phi \phi }-g^2_{t\phi }}, \\ g^{rr}= & \frac{1}{g_{rr}}, ~~~~ g^{\theta \theta } = \frac{1}{g_{\theta \theta }}. \end{aligned}$$

The covariant 4-momentum \(p_{\alpha }\) is governed by one set of the canonical equations of the Hamiltonian system: \(\dot{x}^{\alpha }=\partial K/\partial p_{\alpha }\). Because the Hamiltonian system does not explicitly depend on the coordinates t and \(\phi \), another set of the Hamiltonian canonical equations \(\dot{p}_{\alpha }=-\partial K/\partial x^{\alpha }\) imply that the covariant 4-momentum has two constant components \(p_t=-E\) and \(p_\phi =L\). The two constants of motion satisfy the relations

$$\begin{aligned} E= & -g_{tt}\dot{t} -g_{t\phi }\dot{\phi }, \end{aligned}$$
(19)
$$\begin{aligned} L= & g_{t\phi }\dot{t}+ g_{\phi \phi }\dot{\phi }. \end{aligned}$$
(20)

Note that \(\dot{t}\) and \(\dot{\phi }\) are derivatives of t and \(\phi \) with respect to the affine parameter \(\tau \) rather than the proper time. L and \(\tau \) are given dimensionless quantities after the scale transformations \(L\rightarrow mL\) and \(\tau \rightarrow m\tau \). E and L are the formal energy and angular momentum of the light ray. However, E in general is not the energy measured by a static observer at spatial infinity because the spacetime (1) is not asymptotically flat. To show this fact, we introduce the tetrad basis of an observer at \((r,\theta )\) (see e.g. Ref. [33]): \(e^{\nu }_{\hat{t}}=(\zeta ,0,0,\sigma )\), \(e^{\nu }_{\hat{r}}=(0,D^r,0,0)\), \(e^{\nu }_{\hat{\theta }}=(0,0,D^{\theta },0)\) and \(e^{\nu }_{\hat{\phi }}=(0,0,0,D^{\phi })\). Here, \(\zeta \), \(\sigma \), \(D^r\), \(D^{\theta }\) and \(D^{\phi }\) have the following expressions:

$$\begin{aligned} D^{r}= & \frac{1}{\sqrt{g_{rr}}}, ~~~~D^{\theta } = \frac{1}{\sqrt{g_{\theta \theta }}}, ~~~~ D^{\phi } = \frac{1}{\sqrt{g_{\phi \phi }}}, \\ \zeta= & \sqrt{\frac{g_{\phi \phi }}{g^2_{t\phi }-g_{tt}g_{\phi \phi }}}, ~~~~~~ \sigma =-\frac{g_{t\phi }}{g_{\phi \phi }}\zeta =\omega \zeta . \end{aligned}$$

The four-momentum is locally measured in this frame by \(p^{\hat{\alpha }}=e^{\nu }_{\hat{\alpha }}p_{\nu }\). Thus, the measured energy and angular momentum of the photon are

$$\begin{aligned} E_{obs}= & p^{\hat{t}}=-p_{\hat{t}}=-e^{\nu }_{\hat{t}}p_{\nu }=\zeta E-\sigma L, \end{aligned}$$
(21)
$$\begin{aligned} L_{obs}= & p^{\hat{\phi }}=p_{\hat{\phi }}=e^{\nu }_{\hat{\phi }}p_{\nu }=\frac{L}{\sqrt{g_{\phi \phi }}}. \end{aligned}$$
(22)

Equation (22) is rewritten as

$$\begin{aligned} L_{obs}= & L\sqrt{\frac{\mathcal {R}^2+H_1B+H_2B^2+H_3B^3+H_4B^4}{\Sigma \sin ^2\theta }}, \end{aligned}$$
(23)
$$\begin{aligned} H_4= & \frac{\sin ^4\theta }{16}\left( r^2+a^2\cos ^2\theta \right) \left( r^2+a^2\right) ^2+\cdots . \end{aligned}$$
(24)

For the case of \(B\ne 0\) and \(\theta \ne 0,\pi \), \(L_{obs}\rightarrow \infty \) as \(r\rightarrow \infty \). Namely, the observed angular momentum is nonzero at spatial infinity. Obviously, \(E_{obs}\ne E\) near infinity. In fact, the authors of [32] showed that for the case of \(a=q=0\), \(B\ne 0\) and \(\theta \ne 0,\pi \), \(E_{obs}\rightarrow 0\) as \(r\rightarrow \infty \). They also showed that if \(\theta =0,\pi \), then \(E_{obs}=E\) at spatial infinity. This means that E is the photon’s energy measured by an observer along the symmetric axis at spatial infinity in the Schwarzschild–Melvin spacetime. Therefore, the two constants in Eqs. (19) and (20) are not generally the observed energy and angular momentum in the magnetized Kerr–Newman spacetime (1). On the contrary, \(E_{obs}= E\) and \(L_{obs}\rightarrow 0\) in the case of \(B=0\) as \(r\rightarrow \infty \). Due to the photon with zero rest mass, there is a third constant of motion

$$\begin{aligned} K=0. \end{aligned}$$
(25)

This means that the Hamiltonian is always identical to zero. As Yang and Wu [31] claimed, no fourth constant of motion exists for the particle motion around the black hole because the Hamilton-Jacobi equation is not separable to the variables. The result is also suitable for the photon motion.

2.2 Photon effective potentials and circular orbits

The Hamiltonian K is rewritten as

$$\begin{aligned} K=\frac{1}{2}\left( g^{rr}p_{r}^2+g^{\theta \theta }p_{\theta }^2\right) +\frac{1}{2} {E^2}V_{eff}, \end{aligned}$$
(26)

where the effective potential \(V_{eff}\) on the equatorial plane \(\theta =\pi /2\) is expressed as

$$\begin{aligned} V_{eff}=g^{\phi \phi }\left( \eta ^2+2\eta \frac{g_{t\phi }}{g_{tt}}+\frac{g_{\phi \phi }}{g_{tt}}\right) , \end{aligned}$$
(27)

where \(\eta =L/E\) is a constant parameter. The effective potential can be found in Ref. [33].

If the effective potential satisfies the relations

$$\begin{aligned} \theta= & \frac{\pi }{2}, \end{aligned}$$
(28)
$$\begin{aligned} V_{eff}= & 0, \end{aligned}$$
(29)

we have

$$\begin{aligned} \eta _{1}= & -\frac{g_{t\phi }}{g_{tt}}-\frac{\sqrt{{g_{t\phi }^2}-g_{tt}g_{\phi \phi }}}{g_{tt}}, \end{aligned}$$
(30)
$$\begin{aligned} \eta _{2}= & -\frac{g_{t\phi }}{g_{tt}}+\frac{\sqrt{{g_{t\phi }^2}-g_{tt}g_{\phi \phi }}}{g_{tt}}. \end{aligned}$$
(31)

If the condition

$$\begin{aligned} \frac{ \partial V_{eff} }{ \partial r}=0 \end{aligned}$$
(32)

is also satisfied, then r is the radius of a circular photon orbit or light ring on the equatorial plane. The light ring with its radius \(r_1\) given by Eqs. (30) and (32) is prograde, while the ring with its radius \(r_2\) governed by Eqs. (31) and (32) is retrograde. The ring is unstable when

$$\begin{aligned} {\partial ^2 V_{eff}/\partial r^2}<0, \end{aligned}$$
(33)

whereas it is stable if

$$\begin{aligned} {\partial ^2 V_{eff}/\partial r^2}>0. \end{aligned}$$
(34)

In other words, the local maximum value of the effective potential may correspond to an unstable photon circular orbit satisfying the conditions (30), (32) and (33), but the local minimum value of the effective potential may show the presence of a stable photon circular orbit satisfying the conditions (31), (32) and (34). The conditions for the unstable or stable photon circular orbit were given in [33]. In fact, all the unstable light rings are associated with the boundary of black hole shadows.

We take the parameter \(a=0.2\). The purely electric charge of the black hole should be \(q=-aB\), as is above-mentioned. Such a value should always be given to the purely electric charge in the later computations. Figure 1a plots the effective potentials for several values of the parameter \(\eta \) in the nonmagnetized case \(B=0\). Note that these values \(\eta \) are given arbitrarily and do not necessarily satisfy Eq. (29). The effective potentials have potential barriers corresponding to the local maximum values but do not exist closed pockets (or potential wells) corresponding to the local minimum values for \(\eta =-6.0302\), \(\eta =-5.5856\) and \(\eta =4.7832\). Here, the positive and negative signs of \(\eta \) correspond to the prograde and retrograde directions, respectively. However, neither the local maximum value nor the local minimum value can exist in the effective potential for \(\eta =3.5176\). The effective potential runs from the black hole with the magnitude of \(\eta \) increasing. It is worth noting that the local maximum value of the effective potential for \(\eta =-6.0302\) does not represent the existence of an unstable circular orbit because \(V_{eff}\ne 0\). Nevertheless, the local maximum values of the effective potentials for \(\eta =4.7832\) and \(\eta =-5.5856\) correspond to the existence of unstable circular orbits, which have radii \(r=2.7592\) for \(\eta =4.7832\) and \(r=3.2282\) for \(\eta =-5.5856\). When the magnetic field strength \(B=0.13\) is included, the effective potentials are shown in Fig. 1b, c. A typical difference between Fig. 1a–c is that the effective potentials have the local maximum values for \(\eta =4.4551\), \(\eta =-5.1134\), \(\eta =-5.7780\), \(\eta =5.6760\) and \(\eta =-6.0302\) in the magnetized case (Fig. 1b, c), whereas they do not have for \(\eta =-6.0302\), \(\eta =-5.5856\) and \(\eta =4.7832\) in the nonmagnetized case (Fig. 1a). Of course, the local maximum values are present in the two cases. Particular for \(\eta =-6.0302\), the local maximum value of the effective potential has no unstable circular orbit, and the local minimum one does not indicate the presence of stable circular orbit. The results can be seen clearly from \(V_{eff}\ne 0\). Similarly, no stable circular orbits but unstable ones exist for \(\eta =4.4551\) with the radius \(r=2.8866\) and \(\eta =-5.1134\) with the radius \(r=3.5149\). There are no unstable circular orbits but stable ones for \(\eta =-5.7780\) with the radius \(r=7.7135\) and \(\eta =5.6760\) with the radius \(r=7.9986\).

Fig. 1
figure 1

Effective potentials on the equatorial plane \(\theta =\pi /2\) for several values of \(\eta \). The other parameters are \(a=0.2\) and \(q=-aB\). (a) When the magnetic field is \(B=0\), no stable circular orbits but unstable ones exist for the retrograde direction with \(\eta =-5.5856\) and the prograde direction with \(\eta =4.7832\). (b) When the magnetic field is \(B=0.13\), two unstable circular orbits exist for \(\eta =-5.1134\) and \(\eta =4.4551\). (c) When the magnetic field is still \(B=0.13\), two stable circular orbits exist for \(\eta =-5.7780\) and \(\eta =5.6760\)

Three results can be concluded from Fig. 1. (i) The effective potentials outside the horizons can accept the possibility of the local maximum values but cannot allow for the local minimum values in the nonmagnetized Kerr type spacetimes. Because of this point, no stable circular photon orbits but unstable ones exist. No bound photon orbits without falling to the center or escaping to infinity exist, either. Note that stable equatorial circular photon orbits can exist within the inner horizon in Kerr–Newman spacetimes [22]. However, the effective potentials may have the local maximum values and the local minimum ones in the magnetized Kerr type spacetimes. Thus, unstable and stable circular photon orbits may be allowed. The stable circular photon orbits are a class of bound photon orbits. Other bound photon orbits are also possible. (ii) When the spin parameter a is given, two unstable circular photon orbits and two stable ones are present for an appropriately large magnetic field strength B, while only two unstable circular photon orbits are for an extremely small magnetic field strength B or the vanishing magnetic field. (iii) When the parameters B and a are given, the radii of the stable circular photon orbits are larger than those of the unstable circular photon orbits. This is because the unstable circular photon orbits go toward the black hole, while the stable ones run from the black hole. In addition, the radii of the unstable and stable circular photon orbits are larger for the retrograde black holes than those for the prograde black holes.

Considering that the magnetic fields play a dominant role in the existence of the local minimum values of the effective potentials, we pay attention to the dependence of the radii of unstable or stable circular orbits on the magnetic field strength B. For the prograde case, the circular orbits are solved from Eqs. (30) and (32). Choosing an initial value of r in the range \(r^{+}_{eh}<r<5\), we use the Newtonian iteration method to solve Eqs. (30) and (32). In this way, we can obtain the radii of unstable circular orbits for several values of the spin parameter a in Fig. 2a. The radii of stable circular orbits are also shown in Fig. 2b through an initial value of r in the range \(8<r<15\). For a given value of the spin parameter a in the prograde case, the unstable circular orbit has its larger radius as the magnetic field strength B increases. However, the radius of the stable circular orbit decreases. Particularly for \(B=a=0\), the metric (1) is the Schwarzschild black hole spacetime. Here, the radius of the unstable photon ring is 3, but that of the stable circular orbit tends towards infinity, that is, no stable circular photon orbit exists. When \(B=0\) and \(a\ne 0\), the metric (1) is the Kerr spacetime, where the radii of the unstable circular photon orbits exist in the range \(1.2\le r\le 3\), while those of the stable circular orbits still tend to infinity. In a word, no stable photon rings but unstable ones can exist in the Schwarzschild and Kerr spacetimes. The result is what we have well known. When the magnetic field strength B is given a nonzero value, the radius of the unstable circular orbit decreases but that of the stable circular orbit increases with the spin parameter a increasing. On the other hand, the circular orbits for the retrograde case are governed by Eqs. (31) and (32). The dependence of the radii of the circular orbits on the magnetic field strength B for the retrograde case is similar to that for the prograde case in Fig. 2c, d. However, the dependence of the radii of the circular orbits on the spin parameter a for the retrograde case is unlike that for the prograde case. The radius of the retrograde unstable circular orbit increases but that of the retrograde stable circular orbit decreases when B is given and a gets from small to large. See Fig. 2 for more information about the effects of B or a on the circular photon orbits for the prograde and retrograde two cases in the present problem. The results are consistent with those of [33] in the magnetized Kerr spacetime without the purely electric charges.

Fig. 2
figure 2

Radii r of circular orbits at the equatorial plane \(\theta =\pi /2\) depending on the magnetic field parameter B for several values of the rotation parameter a. (a) The unstable circular orbits in the prograde case with \(\eta _{1}\) of Eq. (30). (b) The stable circular orbits in the prograde case. (c) The unstable circular orbits in the retrograde case with \(\eta _{2}\) of Eq. (31). (d) The stable circular orbits in the retrograde case

Fig. 3
figure 3

The radii r of circular photon orbits in the equatorial plane depending on the magnetic field parameter B and the rotation parameter a. (a) The unstable circular orbits in the prograde case. (b) The stable circular orbits in the prograde case. (c) The unstable circular orbits in the retrograde case. (d) The stable circular orbits in the retrograde case

Besides the two-dimensional plots on the effect of the parameter a or B on the radii of equatorial circular photon orbits in Fig. 2, three-dimensional plots in Fig. 3 well describe the dependence of the radii of equatorial circular photon orbits on the parameters a and B. Some examples of the stable or unstable equatorial circular photon orbits in Fig. 3 are listed in Tables 1 and 2. It is clear that the dependence of the radii r on a or B in Fig. 3 is similar to that in Fig. 2.

In general, the circular photon orbits are determined by Eqs. (29), (32) and the following equation

$$\begin{aligned} {\partial V_{eff}/\partial \theta }=\mathbb {V}(r,\theta )\cos \theta =0. \end{aligned}$$
(35)

Obviously, \(\theta =\pi /2\) is a solution of Eq. (35). This is why the above-mentioned circular photon orbits are discussed in the equatorial plane. The equatorial circular photon orbits are relatively useful to compute black hole shadows in many references. In addition, Eqs. (29), (32) and the equation \(\mathbb {V}(r,\theta )=0\) may have a solution \((r,\theta ,\eta )\). In this case, nonequatorial circular photon orbits are possibly present. Such nonequatorial circular photon orbits are too expensive in computations because the equation (32) and the equation \(\mathbb {V}(r,\theta )=0\) must be solved iteratively after Eq. (30) or (31) is substituted into the two equations. The nonequatorial circular photon orbits are seldom considered in the calculation of black hole shadows in the literature. The reason is that the equatorial circular photon orbits correspond to the boundary of black hole shadows, whereas the nonequatorial ones may not.

3 Bound quasiperiodic and chaotic photon orbits

Now, let us consider the motion of photons outside the equatorial plane. The motion is governed by the Hamiltonian (18) with Eq. (25).

3.1 Explicit symplectic integrators

The authors of [31] introduced the time transformation

$$\begin{aligned} d\tau = g(r,\theta )d w, \end{aligned}$$
(36)

where g is a time transformation function

$$\begin{aligned} g(r,\theta )=\frac{H \mathcal {R} ^2}{r^2}. \end{aligned}$$
(37)

They obtained the time transformation Hamiltonian

$$\begin{aligned} \mathcal {J}= g(K+p_0) =\mathcal {J}_1+\mathcal {J}_2+\mathcal {J}_3+\mathcal {J}_4+\mathcal {J}_5, \end{aligned}$$
(38)

where \(p_0=-K=0\) is a new momentum corresponding to a new coordinate \(q_0=\tau \). Obviously, \(\mathcal {J}=0\). Detailed expressions of \(\mathcal {J}_1\), \(\mathcal {J}_2\), \(\mathcal {J}_3\), \(\mathcal {J}_4\) and \(\mathcal {J}_5\) are

$$\begin{aligned} \mathcal {J}_1= & g(r,\theta )(\mathcal {H}_1+p_0), \end{aligned}$$
(39)
$$\begin{aligned} \mathcal {J}_2= & \frac{1}{2}p_r^2, \end{aligned}$$
(40)
$$\begin{aligned} \mathcal {J}_3= & -\frac{1}{r}p_r^2, \end{aligned}$$
(41)
$$\begin{aligned} \mathcal {J}_4= & \frac{p_\theta ^2}{2r^2}, \end{aligned}$$
(42)
$$\begin{aligned} \mathcal {J}_5= & \frac{1}{2r^2}(a^2+q^2)p_r^2, \end{aligned}$$
(43)

where

$$\begin{aligned} \mathcal {H}_1=\frac{E^2}{2}V_{eff}. \end{aligned}$$
(44)

The five sub-Hamiltonian parts have their analytical solutions, which are explicit functions of the new time w.

The solvers of the five sub-Hamiltonians are \(\mathfrak {T} _1\), \(\mathfrak {T} _2\), \(\mathfrak {T} _3\), \(\mathfrak {T} _4\), and \(\mathfrak {T} _5\), respectively. Taking h as a time step, the authors of [31] gave two first-order operators

$$\begin{aligned} \varXi (h)= & \mathfrak {T}_1(h)\times \mathfrak {T}_2(h)\times \mathfrak {T}_3(h)\nonumber \\ & \times \mathfrak {T}_4(h)\times \mathfrak {T}_5(h), \end{aligned}$$
(45)
$$\begin{aligned} \varXi ^{*} (h)= & \mathfrak {T}_5(h)\times \mathfrak {T}_4(h)\times \mathfrak {T}_3(h)\nonumber \\ & \times \mathfrak {T}_2(h)\times \mathfrak {T}_1(h). \end{aligned}$$
(46)

Using the two operators, they established a second-order explicit symplectic method for the Hamiltonian (38)

$$\begin{aligned} S_2(h)=\varXi ^{*}\left( \frac{h}{2}\right) \times \varXi \left( \frac{h}{2}\right) . \end{aligned}$$
(47)

There is a fourth-order method [58]

$$\begin{aligned} S_4(h)=S_2(\gamma h)\times S_2(\delta h)\times S_2(\gamma h), \end{aligned}$$
(48)

where \(\gamma =1/(2-\root 3 \of {2}) \) and \(\delta =1-2\gamma \). An optimized fourth-order partitioned Runge–Kutta explicit symplectic integrator reads [59,60,61,62]

$$\begin{aligned} PRK_64(h)= & \varXi ^{*} (\alpha _{12} h)\times \varXi (\alpha _{11} h)\times \cdots \nonumber \\ & \times \varXi ^{*} (\alpha _2 h) \times \varXi (\alpha _1 h), \end{aligned}$$
(49)

where the time coefficients are given by

$$\begin{aligned} \alpha _1 = \alpha _{12}= & 0.079203696431196, \\ \alpha _2 = \alpha _{11}= & 0.130311410182166, \\ \alpha _3 = \alpha _{10}= & 0.222861495867608, \\ \alpha _4 = \alpha _9= & -0.366713269047426, \\ \alpha _5 = \alpha _8= & 0.324648188689706, \\ \alpha _6 = \alpha _7= & 0.109688477876750. \end{aligned}$$
Table 1 The radii of unstable (U) and stable (S) circular photon orbits for several values of the magnetic field strength B and a in the prograde case of Fig. 3a and the prograde case of Fig. 3b
Table 2 Same as those in Table 1 but those in Fig. 3c, d

\(S_2\), \(S_4\) and \(PRK_64\) are the explicit symplectic integrators for the time transformation Hamiltonian (38). Fixed time steps are used for the new time w, but variable time steps seem to be considered for the original time \(\tau \). In fact, these methods are not adaptive, and adaptive explicit symplectic integrators for the curved spacetimes are developed recently in Ref. [63].

Let the new time step be \(h=1\). The parameters are given as follows: the magnetic field strength \(B=0.10\), black hole rotating parameter \(a=0.01\) and photon energy \(E=0.995\). The photon parameter is \(\eta _{1} =7.0550\) for Orbit 1, and \(\eta _{1} =6.599\) for Orbit 2. The initial conditions are \(\theta =\pi /2\) and \(p_{r}=0\). The initial separations are \(r=9\) for Orbit 1, and \(r=15\) for Orbit 2. The initial values \(p_{\theta }>0\) are solved from Eq. (25). The three methods give no secular drifts to the Hamiltonian errors \(|\mathcal {J}|\) for Orbit 1 in Fig. 4a and Orbit 2 in Fig. 4b. The accuracy of \(S_4\) is superior to that of \(S_2\), but inferior to that of \(PRK_64\). In view of better accuracy and computational efficiency, the \(S_{4}\) method is used in our later calculations.

Fig. 4
figure 4

Hamiltonian errors \(|\mathcal {J}|\) in Eq. (38) for the three sympletic methods \(S_{2}\), \(S_{4}\) and \(PRK_{6}4\). The parameters are \(B=0.1\), \(a=0.01\) and \(E=0.995\). The initial conditions are \(p_r=0\) and \(\theta =\pi /2\). (a) \(\eta _{1}=7.0550\) and the initial separation \(r=9\) for Orbit 1. (b) \(\eta _{1}=6.599\) and the initial separation \(r=15\) for Orbit 2

3.2 Chaos indicators

Several methods are used to investigate the regular or chaotic dynamical features of the two tested orbits in Fig. 4.

3.2.1 Poincaré section

The phase space structures of the two orbits in Fig. 4 can be examined through the Poincaré section map in the two-dimensional plane r-\(p_r\).

There are many intersections of the photons’ trajectories with the surface of section \(\theta =\pi /2\) and \(p_\theta >0\). Suppose that \((r,p_r,\theta )\) correspond to a set of solutions \((X_1,Y_1,\Theta _1)\) at some integration step before a photon’ trajectory crosses the section surface \(\theta =\pi /2\), where \(\Theta _1<\pi /2\). \((r,p_r,\theta )\) have another set of solutions \((X_2,Y_2,\Theta _2)\) in the next integration step after this trajectory crosses the section surface, where \(\Theta _2>\pi /2\). Using the linear interpolation method, one can get a line with the two points \((X_1,\Theta _1)\) and \((X_2,\Theta _2)\). This line is described by a point \((X,\Theta )\) satisfying the function

$$\begin{aligned} X=\frac{X_2-X_1}{\Theta _2-\Theta _1}(\Theta -\Theta _1)+X_1. \end{aligned}$$
(50)

Letting \(\Theta =\pi /2\), one has X corresponding to r as a coordinate of an intersection point on the section surface. Similarly, the two points \((Y_1,\Theta _1)\) and \((Y_2,\Theta _2)\) determine the function of point \((Y,\Theta )\):

$$\begin{aligned} Y=\frac{Y_2-Y_1}{\Theta _2-\Theta _1}(\Theta -\Theta _1)+Y_1. \end{aligned}$$
(51)

For \(\Theta =\pi /2\), Y is known and corresponds to \(p_r\) as another coordinate of the intersection point on the section surface. In this way, the intersection point \((r,p_r)\) on the section surface is obtained and the Poincaré section map can be plotted.

Fig. 5
figure 5

Dynamical features of Orbits 1 and 2 in Fig. 4. (a), (c) Poincaré section maps on the plane \(\theta =\pi /2\). (b) Evolution of \(\theta \) with time w for Orbit 1 in panel (a)

The intersections \((r,p_r)\) form a closed curve, which shows the regularity of Orbit 1 in Fig. 5a. As is above-mentioned, the intersections show that the photon orbits do frequently cross the equatorial plane. This point is also shown clearly in Fig. 5b, which relates to the evolution of \(\theta \) running from 0 to \(\pi \). However, the intersections \((r,p_r)\) form an area with many randomly points distributed and indicate the chaoticity of Orbit 2 in Fig. 5c.

3.2.2 Lyapunov exponents and fast Lyapunov indicators

The dynamical features of the two orbits are confirmed by the maximal Lyapunov exponents in Fig. 6a. With the aid of the two-nearby-orbits method in a curved spacetime, the Lyapunov exponents are calculated in [45] by

$$\begin{aligned} \lambda =\lim _{w\rightarrow \infty } \frac{1}{w}\ln \frac{d(w)}{d0}, \end{aligned}$$
(52)

where d0 and d(w) are the proper distances of two nearby orbits at the starting time and time w, respectively. If the maximal Lyapunov exponent \(\lambda \) is zero, then the considered bounded orbit is regular. When the maximal Lyapunov exponent \(\lambda \) is positive, the considered bounded orbit is chaotic.

Fig. 6
figure 6

Fig. 5 (continued). (a) The largest Lyapunov exponents \(\lambda \). (b) The fast Lyapunov indicators (FLIs). (c), (d) The values of Kc and the visuals of point \((P_c, Q_c)\). These methods provide the consistent results on regular and chaotic dynamics of the two orbits

The dynamical features of the two orbits are also shown through the fast Lyapunov indicators (FLIs) in Fig. 6b. Based on a modified version of the FLI in [46], the FLI with two nearby orbits was defined in an invariant form [47]:

$$\begin{aligned} FLI=\log _{10}\frac{d(w)}{d0}. \end{aligned}$$
(53)

Orbit 1 is regular because its FLI grows slowly with time \(\log _{10}w\), whereas Orbit 2 is chaotic because its FLI grows exponentially.

3.2.3 The 0–1 test correlation method

In addition to the above three methods finding chaos, the 0–1 test correlation method of Gottwald and Melbourne [53] can well characterize the chaotic or regular nature of the two orbits. It is described as follows.

Let \(\psi (j)\) with \(j=1,2,3,\ldots ,n_1\) be a time series. Two translation variables are given by

$$\begin{aligned} P_c(n_1)= & \sum _{j= 1}^{n_1} \psi (j)\cos (cj), \end{aligned}$$
(54)
$$\begin{aligned} Q_c(n_1)= & \sum _{j= 1}^{n_1} \psi (j)\sin (cj), \end{aligned}$$
(55)

where \(c\in (0,\pi )\) is a parameter and \(n_1=1,2,\ldots ,N_1\). The mean square displacement between the translation variables \(P_c\) and \(Q_c\) is computed by

$$\begin{aligned} M_c(n_2)= & \lim _{N_2 \rightarrow \infty } \frac{1}{N_2-n_2}\sum _{j = 1}^{N_2-n_2} [P_c(j+n_2)-P_c(j)]^2 \nonumber \\ & +[Q_c(j+n_2)-Q_c(j)]^2, \end{aligned}$$
(56)

where \(n_2\ll N_2\).

A modified mean square displacement is considered as

$$\begin{aligned} D_c(n_2) = M_c(n_2) -V_{osc}(c,n_2), \end{aligned}$$
(57)

where \(V_{osc}(c,n_2)\) is an oscillatory term

$$\begin{aligned} V_{osc}(c,n_2)= & \Lambda ^2 \frac{1-\cos (cn_2)}{1-\cos (c)}, \end{aligned}$$
(58)
$$\begin{aligned} \Lambda= & \lim _{N_2 \rightarrow \infty } \frac{1}{N_2}\sum _{j = 1}^{N_2}\psi (j). \end{aligned}$$
(59)

\(D_c(n_2)\) and \(M_c(n_2)\) have the same asymptotic growth rate with time \(n_2\). However, \(D_c(n_2)\) has better convergence properties than \(M_c(n_2)\) because the oscillatory term is eliminated in \(D_c(n_2)\), while is not in \(M_c(n_2)\).

Fig. 7
figure 7

Photon motion orbits with the magnetic field B and the initial separation r. Black: the photons falling to the center, Purple: the photons scattering to infinity, Yellow: unstable circular photon orbits, Green: bound stable circular photon orbits, Blue: bound quasiperiodic orbits and Red: bound chaotic orbits. The parameter is \(a=0.80\). (a) \(\eta _{1}\) in the prograde situation. (b) \(\eta _{2}\) in the retrograde situation

The vectors \(\xi =(1,2,\ldots ,\tilde{n})\) and \(\Delta _c =(D_c(1),D_c(2), \ldots ,D_c(\tilde{n}))\) correspond to mean values

$$\begin{aligned} \bar{\xi }= & \frac{1}{\tilde{n}}\sum ^{\tilde{n}}_{j=1} j=\frac{1}{2}(\tilde{n}+1), \\ \bar{\Delta }_c= & \frac{1}{\tilde{n}}\sum ^{\tilde{n}}_{j=1}D_c(j). \end{aligned}$$

Here, \(n_2\le \tilde{n}\le N_2/10\). Based on the two vectors, the variances are

$$\begin{aligned} var(\xi )= & \frac{1}{\tilde{n}}\sum ^{\tilde{n}}_{j=1} (j-\bar{\xi })^2, \end{aligned}$$
(60)
$$\begin{aligned} var(\Delta _c)= & \frac{1}{\tilde{n}}\sum ^{\tilde{n}}_{j=1} (D_c(j)-\bar{\Delta }_c)^2. \end{aligned}$$
(61)

For \(\xi \) and \(\Delta _c\), their covariance is

$$\begin{aligned} cov(\xi ,\Delta _c)= \frac{1}{\tilde{n}}\sum ^{\tilde{n}}_{j=1} (j-\bar{\xi })(D_c(j)-\bar{\Delta }_c). \end{aligned}$$
(62)

The correlation coefficient of the two vectors is defined as

$$\begin{aligned} K_c = corr(\xi ,\Delta _c) =\frac{cov(\xi ,\Delta _c)}{\sqrt{var(\xi )~var(\Delta _c)}}. \end{aligned}$$
(63)

The mean square displacement \(D_c(n_2)\) is bounded as time \(n_2\) increases for a nonchaotic time series, whereas grows linearly when time number \(n_2\) increases for a chaotic time series. Namely, \(K_c = 0\) shows regular dynamics, but \(K_c = 1\) indicates chaotic dynamics. \(K_c\) is the 0–1 chaos test from the correlation method of Gottwald and Melbourne [53].

Orbit 1 with the coefficient \(K_c=0.00035\approx 0\) is regular, but Orbit 2 with the coefficient \(K_c=0.99915\approx 1\) is chaotic. The visual plot on the evolution behavior of \(P_c\) and \(Q_c\) that is a bounded region in Fig. 6c underlies regular dynamics. However, the visual plot that is an asymptotic unbounded Brownian motion in Fig. 6d characterizes chaotic dynamics.

Table 3 The radii of unstable (U) and stable (S) circular photon orbits for several values of the magnetic field strength B in the prograde case of Fig. 7a
Table 4 Same as those in Table 3 but those in the retrograde case of Fig. 7b

Besides the afore-mentioned four techniques, many other tools can demonstrate the origins of chaos in the system. For instance, homoclinic or heteroclinic intersections of the stable and unstable manifolds in the Poincaré sections or the horseshoe dynamics show the onset of chaos [64, 65].

Fig. 8
figure 8

(a)–(c) The evolution of the photon motion radius r with the new time w. The parameters are \(E=0.998\), \(\eta _{1}\) and \(a=0.80\). (a) Three escaping orbits. (b) Three falling orbits. (c) Three stable circular photon orbits. (d) The visuals of \(P_{c}\) and \(Q_{c}\) with the value of \(K_{c}\) for the three orbits in panel (c)

Fig. 9
figure 9

(a)–(e) Poincaré maps. (a) \(\eta _{1}=4.0881\), \(E=0.998\); (b) \(\eta _{1}=4.2499\), \(E=0.998\); (c) \(\eta _{1}=4.2527\), \(E=0.998\); (d) \(\eta _{1}=4.1321\), \(E=0.998\); (f) \(\eta _{1}=3.8004\), \(E=0.998\). The other parameters are \(B=0.175\) and \(a=0.80\). (f) Lyapunov exponents \(\lambda \). It is clear that bound regular and chaotic photon orbits can exist in the vicinity of the magnetized Kerr–Newman black hole

Fig. 10
figure 10

(a)–(c) Poincaré maps. (a) \(\eta _{1}=4.2813\), \(E=0.998\); (b) \(\eta _{1}=4.5703\), \(E=0.998\); (c) \(\eta _{1}=4.4239\), \(E=0.998\). The other parameters are \(B=0.16\) and \(a=0.80\). (df) The visuals of \(P_{c}\) and \(Q_{c}\) with the values of \(K_{c}\), which correspond to panels (a)–(c). Clearly, bound regular and chaotic photon orbits exist in the vicinity of the magnetized black hole

Fig. 11
figure 11

(a)–(e) The visuals of \(P_{c}\) and \(Q_{c}\) with the values of \(K_c\). (a) \(\eta _{1}=3.7835\), \(E=0.998\); (b) \(\eta _{1}=3.9701\), \(E=0.998\); (c) \(\eta _{1}=3.9699\), \(E=0.998\); (d) \(\eta _{1}=3.7425\), \(E=0.998\); (f) \(\eta _{1}=3.4738\), \(E=0.998\). The other parameters are \(B=0.19\) and \(a=0.80\). (f) FLIs. There are bound regular and chaotic photon orbits in the vicinity of the magnetized black hole

In a word, the above demonstrations sufficiently show the quasiperiodicity of Orbit 1 and the chaoticity of Orbit 2. The two photon orbits stay at the allowed regions without falling to the center or escaping to infinity. In this sense, Orbit 1 is a bound quasiperiodic orbit, and Orbit 2 is a bound chaotic orbit. The two orbits stay on the equatorial plane at the starting time, but run away from the equatorial plane in many other times. However, such bound photon orbits do not generally exist in the nonmagnetized black hole spacetime (1) with \(B=0\). These different results are caused by distinct shapes of the effective potentials in the two cases. It is possible to exist the local maxima but impossible to have the local minima for the effective potentials outside the horizon in the spacetime with \(B=0\), as shown in Fig. 1a. Therefore, there are three possible photon orbits around the rotating charged black hole, which correspond to the photons falling to the center, the photons escaping to infinity and the photons unstably circling in the center. On the other hand, the effective potentials in the magnetized spacetime (1) with \(B\ne 0\) may have the local maxima and the local minima in Fig. 1b, c. This leads to the photon motions having the additionally allowed bound orbits without falling to the center or escaping to infinity as well as the three cases with falling to the center, scattering to infinity or unstably circling in the center. The bound photon orbits include the stable circular orbits, bound quasiperiodic orbits and bound chaotic orbits. The five types of photon orbits (i.e. falling orbits, escaping orbits, unstable circular orbits, stable circular orbits, bound quasiperiodic orbits) are not necessarily restricted to the equatorial plane. However, the bound chaotic photon orbits must be allowed outside the equatorial plane.

3.3 More bound photon orbits

To clearly know the presence of more bound photon orbits, we trace the FLI depending on the initial radius r and the magnetic field strength B. Here, the parameter \(a=0.8\) is fixed. The initial values \(\theta =\pi /2\) and \(p_r=0\) are still used. The parameter \(\eta _{1}\) is chosen in the prograde situation and \(\eta _{2}\) is adopted in the retrograde situation. The initial radius r ranges from 0 to 20. The magnetic field strength B is varied in an interval (0, 0.45] for the prograde case, and another interval (0, 0.16] for the retrograde case. A photon with the radial distance \(r>10,000\) is considered in the scattering region. If the radial distance r approaches \(r^{+}_{eh}\), the photon is believed to fall to the center. Except the two types of photon orbits, each of the FLIs for other photon orbits is obtained after the integration time \(10^6\). The threshold of the FLIs between regular dynamics and chaotic dynamics is 5 for the orbits in the bound regions without falling to the center or escaping to infinity. The FLIs larger than this threshold mean chaotic dynamics, while those smaller than or equal to this threshold show regular dynamics.

Figure 7a describes six types of photon orbits in the plane of the initial radius r and the magnetic field strength B for the prograde case. A large region colored black is the falling orbits. The purple region corresponds to the scattering orbits. The unstable circular orbits are existent in the yellow region, whereas the stable ones are given in the green region. The blue region shows bound regular dynamics, but the red one represents bound chaotic dynamics. In short, the 6 types of photon motion orbits are shown clearly. Some of the unstable and stable circular orbits are listed in Table 3. The similar results can be also seen from Fig. 7b on the photon motion orbits for the retrograde case. There are some of the unstable and stable circular orbits in Table 4.

Figure 8a is used to check the three orbits with the magnetic field \(B=0.01\) and the initial separations \(r=4\), 10, 16 in the purple area of Fig. 7a. The three orbits escape quickly from the black hole after short integration times. Figure 8b shows that the three orbits with the values of (Br) corresponding to (0.25, 15), (0.32, 16) and (0.35, 17) in the black area of Fig. 7a go towards the black hole. Figure 8c confirms that the radii remain invariant for the values of (Br) corresponding to (0.16, 6.62243), (0.175, 6.600962) and (0.19, 5.49422) in the green area. The three stable circular orbits are also shown by the 0–1 test method with \(K_c=0.00039 \approx 0\) and the visual plots of \(P_c\) and \(Q_c\) with closed curves in Fig. 8d.

The Poincaré maps and the Lyapunov exponents in Fig. 9 show that the initial separations \(r=4.5\), 7.5, 9.0 correspond to the onset of chaos for \(B=0.175\) in the blue and red regions of Fig. 7a. However, the initial separations \(r=5.5\), 6.5 allow for the regularity. For \(B=0.175\) in Fig. 10, chaotic dynamics occurs for the initial separations \(r=4.5\), 8.5, but regular dynamics appears for the initial separation \(r=7.5\). When \(B=0.19\) is chosen in Fig. 11, the visual plots \((P_c, Q_c)\), the correlation coefficients \(K_c\) and the FLI show the chaotic dynamics for \(r=4.0\), 7.5, 8.5, but the regular dynamics for \(r=5.2\), 5.8.

In brief, Figs. 5, 6, 8, 9, 10 and 11 sufficiently show the existence of the bound photon orbits, including the stable circular orbits, regular quasi-periodic orbits and chaotic orbits. The bound regular orbits and chaotic orbits lie outside the equatorial plane.

4 Summary

When the external magnetic field surrounding the Kerr–Newman black hole is extremely strong, it can make the black hole be magnetized and change the spacetime geometry. When the condition on the total electric charge or the purely electric charge of the black hole is satisfied, the magnetized Kerr–Newman black hole metric is asymptotic to the static Melvin metric but is not asymptotically flat. The external magnetic field greatly affects the motion of neutral or charged particles in the vicinity of the black hole, as was reported in [31]. It also exerts a typical influence on the motion of photons around the black hole. This is because the strong magnetic field in the magnetized spacetime gives a gravitational attraction to the photons and leads to a lack of the fourth motion constant.

The photon effective potentials outside the horizons on the equatorial plane in the magnetized Kerr–Newman spacetime with the condition on the total electric charge or the purely electric charge of the black hole the are unlike those in the nonmagnetized one. They can allow for the existence of potential barriers corresponding to the local maximum values but exclude that of closed pockets (or potential wells) corresponding to the local minimum values in the nonmagnetized case. In this case, no stable circular photon orbits but unstable ones exist. There are no bound photon orbits without falling to the center or escaping to infinity, either. However, the photon effective potentials are possible to have both the local maximum values and the local minimum ones in the magnetized case. Thus, two unstable circular photon orbits and two stable ones can exist for the given parameters including an appropriately large magnetic field strength. Besides the stable circular photon orbits, other bound photon orbits such as quasiperiodic orbits and chaotic orbits are also possible. In sum, there are six possible types of photon orbits: the photons falling to the center, scattering to infinity and unstably circling in the center, and the bound photon orbits including the stable circular orbits, bound quasiperiodic orbits and bound chaotic orbits. The former five types of photon orbits are allowed on and outside the equatorial plane, but the sixth type orbits (i.e. the bound chaotic photon orbits) are only allowed outside the equatorial plane. These photon orbits are confirmed numerically. In particular, the techniques of Poincaré maps, Lyapunov exponents, fast Lyapunov indicators and 0–1 test correlation method show the existence of bound quasiperiodic orbits and bound chaotic orbits for an appropriately large magnetic field strength.

It is worth pointing out that the six types of photon orbits for an appropriately given magnetic field strength (e.g. \(B=0.2\) in Fig. 5a) can correspond to those on the observer’s screen. In other words, the celestial coordinates in the observer’s position locate at the six regions. In fact, a new backward ray-tracing method is used to find the six regions on the observer’s plane. This method includes the time transformation explicit symplectic integrator and the fast Lyapunov indicators or the 0–1 test correlation method. However, the technique on homoclinic or heteroclinic intersections of the stable and unstable manifolds in the Poincaré sections is not conveniently included in the new backward ray-tracing method. Another notable point is that the existing backward ray-tracing methods (see e.g. [32]) adopt traditional numerical integrators such as Runge–Kutta methods to find two regions on the observer’s plane, where one region corresponds to the photons falling to the center and the other region corresponds to the photons not falling to the center. The traditional integrators would provide unreliable results if the integration time is long enough. Clearly, the new ray-tracing method will provide more dynamical information to the photons on the observer’s plane. All points corresponding to the photons falling to the center form a black hole shadow on the observer’s plane. The unstable photon rings are at the boundary of the black hole shadow. In this way, the chaotic behavior may be observed through the black hole image. The new idea for finding the six regions on the observer’s plane is a preliminary to the study of black hole shadows in our future works.