1 Introduction

Regarding the low energy limit of gravitational interactions, general relativity (GR) is a successful theory describing various phenomena. Unlike the successes of the mentioned theory, some issues such as accelerated expansion of the universe, the existence of dark matter [1], massive gravitons and several other subjects show the necessity of modifying this theory. There are various attempts for modifying GR, such as F(R) gravity [2,3,4,5,6,7], Lovelock gravity [8,9,10,11,12], Horava-Lifshitz gravity [13, 14], brane world cosmology [15,16,17], scalar-tensor theories [18,19,20,21,22,23,24,25,26], massive gravity [27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44], rainbow gravity [45,46,47,48,49,50,51,52], and massive gravity’s rainbow [53,54,55]. The mentioned generalized theories lead to interesting results in various aspects of black objects, especially their geometrical and physical properties. In this regard, we are going to study the effects of massive modification of Einstein gravity on the geometric structure of low dimensional black holes, especially, the geodesic motions of a test particle. The primary motivation of studying the geodesic motions of test particles around the massive objects come from some interesting astrophysical phenomena, such as the precession of the perihelion of Mercury and gravitational lensing. Due to the strong curvature effects, black holes have considerable influence on the geodesic motions.

One of the predictions of GR is the existence of black holes. These mysterious singular solutions can be described by the first static spherically symmetric solution of Einstein’s field equations found by Karl Schwarzschild in 1916. Since the interpretation of the Schwarzschild horizon as a one way membrane by Finkelstein in 1958 [56], considerable efforts have been made to study exact solutions of Einstein gravity with a singularity and their interesting properties. One major achievement in the theoretical analysis of these objects was finding the first three dimensional solution of Einstein field equations by Banados-Teitelboim-Zanelli (BTZ) in 1992 [57]. Later on, many people studied various aspects of BTZ black hole as a typical laboratory of high energy physics [27, 58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76] . BTZ black holes provide good tools for handling some conceptual questions in the context of AdS/CFT correspondence, quantum gravity, string theory models and gauge field theory [77, 78]. In the context of classical gravity, BTZ black holes have been studied in the presence of other gauge fields such as Maxwell and power Maxwell fields [79,80,81], Born–Infeld theory [82,83,84,85,86] and other nonlinear electrodynamics [87,88,89], and also massive gravity [90], gravity’s rainbow [91], massive gravity’s rainbow [92] and dilaton field [81, 93, 94]. Also higher dimensional black hole solutions with BTZ analogy (BTZ like black holes) have been investigated in Refs. [95, 96].

We can investigate the black hole properties in various ways by considering the motion of a test particle in the spacetime with a black hole, described by the geodesic equation. Exact analytical solutions to the geodesic equation provide us with the best basic understanding of geodesic motion, but may not always be possible. We can then employ analytical approximation schemes or numerical solutions. A first seminal paper for studying the geodesic equations and its analytical solutions was by Hagihara in 1931, when he solved the geodesic equations of Schwarzschild black holes in terms of Weierstrass elliptic functions [97]. Afterwards, Darwin solved the geodesic equations by using the Jacobian elliptic functions [98, 99]. Also, the analytical solution of geodesic equations have been studied in four-dimensional Schwarzschild–de Sitter [100], Kerr [101, 102], Kerr–Newman [103], and Kerr–de Sitter [104] spacetime. Moreover, the higher dimensional solution of Schwarzschild, Schwarzschild-(anti-) de Sitter, Reissner–Nordström, Reissner Nordström–(anti-) de Sitter [105] and Myers-Perry spacetimes [106] have been investigated. Studying the geodesic equations have been extended to F(R) gravity and conformal gravity in BTZ and GMGHS (Gibbons-Maeda-Garfinkle-Horowitz-Strominger) black holes [107,108,109,110,111], regular and modified Hayward black holes [112, 113], charged dilatonic black holes [114], the singly spinning and (charged) doubly spinning black ring [115, 116], brane-world black hole [117], (rotating) black string [118], Schwarzschild, and Kerr pierced by black string spacetimes [119, 120]. Moreover, regarding three dimensional spacetime, the stability and existence of circular geodesics in a family of asymptotically AdS black holes in new massive gravity theory [121] and the null geodesics in a static circularly symmetric black hole spacetime [122] have been investigated before. Furthermore, in Ref. [123], the exact solutions of null and timelike geodesics are found and it is shown that the rotating (static) black hole is geodesically complete (incomplete).

In this paper, we consider (charged) BTZ black holes and its generalization to massive gravity. We provide a complete classification of the geodesic motion, and solve the geodesic equation. As far as possible, we follow the method which is introduced in Ref. [100]. For all of the previous works on the solutions of geodesic equations, the metric function was a polynomial function. In this paper, we investigate both polynomial and non-polynomial metric functions. A non-polynomial metric function appears in the charged cases and since there is no known exact analytical solution for this type, we resort to numerical solutions for this case.

The outline of this paper is as follow. First, we present the geodesic equation and effective potential in BTZ and its extension to massive gravity black holes in Sect. 2, and define acceptable regions of motion. After that, we investigate possible regions of motion and orbits for charged and uncharged black holes in Sect. 3. Then, we completely classify geodesic motions and solve the geodesic equation around the (charged) BTZ black holes in Sect. 4. The neutral black holes in massive gravity are treated in Sect. 5 and the charged ones in Sect. 6. Finally, we end the paper with some concluding remarks.

2 Three dimensional line element

In \((2+1)-\)dimensional spacetimes, the line element takes the following form

$$\begin{aligned} ds^{2}=-\psi (r)dt^{2}+\frac{dr^{2}}{\psi (r)}+r^{2}d\varphi ^{2}, \end{aligned}$$
(1)

where \(\psi (r)\) is an arbitrary function of coordinate r which should be obtained based on the field equations. In this paper, we will consider four different cases for the metric function, \(\psi (r)\). The first one describes static BTZ black holes [57]

$$\begin{aligned} \psi _{1}(r)=-\Lambda r^{2}-m_{0}, \end{aligned}$$
(2)

where \(\Lambda \) denotes the cosmological constant and \(m_{0}>0\) is an integration constant proportional to the total mass of the black hole. Note that here \(\Lambda <0\) is necessary to have a black hole solution (for a discussion of the singularity at \(r=0\) see [58]). Geodesic motions of massive/massless particles around the rotating uncharged BTZ black holes were investigated in Refs. [124, 125]. In the mentioned papers, it was shown that for the static case, there is no bound orbit for the positive geometrical mass of the black holes (\(m_{0}\)).

The line element (1) was generalized to the linearly charged solution of BTZ black holes with the metric function [95]

$$\begin{aligned} \psi _{2}(r)=-\Lambda r^{2}-m_{0}-2q^{2}\ln \left( \frac{r}{r_{0}}\right) . \end{aligned}$$
(3)

Both \(m_{0}>0\) and q are integration constants which are, respectively, related to mass and electric charge of the black hole solutions. In order to have an event horizon, here, \(\Lambda \) should be negative. It is also worth mentioning that \(r_{0}\) is an arbitrary constant with length dimension which is necessary to obtain dimensionless argument for the logarithmic function. In general, \(r_{0}\) is different from the length scale related to the cosmological constant (\(\Lambda \propto -l^{-2}\)), however, in [126], it was shown that the equality \(r_{0}=l\) is necessary to avoid an ensemble dependency.

The geodesic equation of charged BTZ black holes have been studied in [107]. Since the charge term of the metric function in Ref. [107] is positive, its related black hole has a Schwarzschild like horizon (spacelike singularity). Here, we consider real valued charged BTZ black hole with two horizon (timelike singularity). The roots of the metric function (3) have been reported in Ref. [95]

$$\begin{aligned} r_{+}= & {} r_{0}\exp \left( \frac{-1}{2}\left[ \mathrm {W}\left( \frac{\Lambda r_{0}^{2}e^{\left( -\frac{m_{0}}{q^{2}}\right) }}{q^{2}}\right) +\frac{m_{0} }{q^{2}}\right] \right) , \end{aligned}$$
(4)
$$\begin{aligned} r_{-}= & {} r_{0}\exp \left( \frac{-1}{2}\left[ \mathrm {W}\left( -1,\frac{ \Lambda r_{0}^{2}e^{\left( -\frac{m_{0}}{q^{2}}\right) }}{q^{2}}\right) + \frac{m_{0}}{q^{2}}\right] \right) , \end{aligned}$$
(5)

where \(\mathrm {W}(x)\) denotes the principal branch of the Lambert W function and \(\mathrm {W}(-1,x)\) the branch with \(\mathrm {W}(-1,x)\le -1\). The only physically acceptable solution of Eqs. (4) and ( 5) is

$$\begin{aligned}&r_{+}>0 \ \ \text {or} \ \ r_{-}>0\ \ \text { so} \ \ m_{0}\ge q^{2}\left( \ln \left( -\frac{\Lambda r_{0}^{2}}{q^{2}}\right) +1\right) \nonumber \\&\quad \text {for} \ \ \Lambda <0. \end{aligned}$$
(6)

The metric function of BTZ black holes in massive gravity is [90]

$$\begin{aligned} \psi _{3}(r)=-\Lambda r^{2}-m_{0}+m^{2}cc_{1}r, \end{aligned}$$
(7)

where \(m_{0}>0\) as before and m, c and \(c_{1}\) are three constants related to the massive gravity (see Ref. [90] for more details). Here, we define a new massive parameter \(m^{\prime }\) to combine all massive parameters, as follow

$$\begin{aligned} m^{\prime }=m^{2}cc_{1}. \end{aligned}$$
(8)

Note that in the massive case \(\Lambda <0\) is not necessary for black hole solutions and we may as well consider \(\Lambda >0\).

Following Ref. [90], we find that by considering specific values for different parameters, the metric function (7) may have two roots or one extreme root (we are not interested in naked singularity). In addition, since the constant c is positive, the sign of \( m^{\prime }\) depends on the positive or negative sign of \(c_{1}\) (see [30] for more details). The uncharged BTZ black holes in massive gravity have a curvature singularity at \(r=0\).

In the charged black holes in massive gravity, the metric function \(\psi (r)\) is obtained as [90]

$$\begin{aligned} \psi _{4}(r)=-\Lambda r^{2}-m_{0}+m^{\prime }r-2q^{2}\ln \left( \frac{r}{ r_{0}}\right) . \end{aligned}$$
(9)

Following Refs. [79,80,81, 90], one finds the mentioned metric functions describe a three dimensional spacetime with a singularity located at \(r=0\). It is also notable that the singularity of the charged solutions is timelike, while it is spacelike for uncharged ones (for more details regarding the horizon and geometry of the mentioned spacetimes, we refer the reader to [79,80,81, 90]).

We will analyze all four cases introduced above in terms of the motion of test particles, both for massless and massive particles. Since we are working in static and spherical symmetric spacetimes, we can immediately obtain two conserved quantities, energy and angular momentum, as

$$\begin{aligned} E= & {} g_{tt}\frac{dt}{d\lambda }= -\psi \left( r\right) \frac{dt}{d\lambda } \,, \nonumber \\ L= & {} g_{\varphi \varphi }\frac{d\varphi }{d\lambda }=r^{2}\frac{d\varphi }{ d\lambda }. \end{aligned}$$
(10)

The Lagrangian \({\mathcal {L}}\) of a test particle is given by

$$\begin{aligned} {\mathcal {L}}= & {} g_{\mu \nu }\frac{dx^{\mu }}{d\lambda }\frac{dx^{\nu }}{d\lambda }=-\epsilon =-\psi \left( r\right) \left( \frac{dt}{d\lambda }\right) ^{2} \nonumber \\&+ \frac{1}{\psi \left( r\right) }\left( \frac{dr}{d\lambda }\right) ^{2}+r^{2}\left( \frac{d\varphi }{d\lambda }\right) ^{2}, \end{aligned}$$
(11)

where \(\epsilon \) takes the values 1 and 0 for massive and massless particles, respectively, and \(\lambda \) is the affine parameter for massless particles and the proper time for massive particles.

From this, we then find the following geodesic equations

$$\begin{aligned} \left( \frac{dr}{d\lambda }\right) ^{2}&=E^{2}-\psi (r)\left( \frac{L^{2}}{ r^{2}}+\epsilon \right) , \end{aligned}$$
(12)
$$\begin{aligned} \left( \frac{dr}{d\varphi }\right) ^{2}&=\frac{r^{4}}{L^{2}}\left\{ E^{2}-\psi (r)\left( \frac{L^{2}}{r^{2}}+\epsilon \right) \right\} :=P(r). \end{aligned}$$
(13)

In addition, considering Eq. (12), an effective potential can be introduced as

$$\begin{aligned} V_{\mathrm {eff}}=\psi \left( r\right) \left( \frac{L^{2}}{r^{2}}+\epsilon \right) . \end{aligned}$$
(14)

Before solving the equations of motion we analyze the structure of possible types of geodesic motion for the test particles. The major point in these analyzes is that Eq. (13) implies \(P(r)\ge 0\) as a necessary condition for the existence of a geodesic. The real roots of P(r) are related to intersection points of \(E^2\) and \(V_{\mathrm {eff}}\). Equivalently, from (14) the acceptable region of motion is \( E^2\ge V_{\mathrm {eff}}\). The number of real roots of P(r) characterize the shape of the orbit [100]. For a given set of parameters, P(r) may have a certain number of roots. By varying E and L, and fixing other parameters the number of real roots and, therefore, the possible types of orbits, will change [101].

3 Classification of geodesic motion

Changes in the possible orbit configurations happen if two real zeros of P(r) merge to a double zero. The corresponding constants of motion can be obtained by solving \(P(r)=0\) and \(\frac{dP(r)}{dr}=0\) for \(E^{2}\) and \(L^{2}\) . Obtaining \(E^{2}\) and \(L^{2}\) for these limit cases is therefore crucial for studying the motion of particles. In all four cases introduced in Sect. 2, the function P(r) has a double root \(r=0\) for all parameters. Therefore, we can introduce a new function \({\widetilde{P}}(r)\) by

$$\begin{aligned} P(r)=r^{2}{\widetilde{P}}(r). \end{aligned}$$
(15)

Instead of searching for roots of P(r), we may solve \({\widetilde{P}}(r)=0\) and \(\frac{d{\widetilde{P}}(r)}{dr}=0\) for \(E^{2}\) and \(L^{2}\), for both massive and massless particles.

The asymptotic behavior of P(r) is very important since it determines whether particles have flyby or bound orbits. A particle reaches infinity if P(r) is positive for \(r\rightarrow \infty \). In the limit of large r we find

$$\begin{aligned} P(r)&\rightarrow \text { sign}\left( \Lambda \right) \infty \quad \text {for }\epsilon =1 \end{aligned}$$
(16)
$$\begin{aligned} P(r)&\rightarrow \text { sign}\left( \frac{E^{2}}{L^{2}}+\Lambda \right) \infty \quad \text {for }\epsilon =0. \end{aligned}$$
(17)

For negative \(\Lambda \), massive particles may therefore never reach infinity.

3.1 Uncharged cases

In the uncharged case, the metric function is a polynomial function. By considering P(r), we can discuss about possible motions. In the uncharged cases that we study in this paper, P(r) is a polynomial in r of order not more than 6. Therefore, \({\widetilde{P}}(r)\) is of order not more than 4.

For all uncharged cases the following kinds of orbits can be identified.

  • Flyby orbits: r starts from \(\infty \), then approaches a periapsis \( r=r_{\min }\) and goes back to \(\infty \).

  • Bound orbits: r changes between two boundary values \(r_{1} \le r \le r_{2}\) with \(0< r_{1} \le r_{2} < \infty \).

  • Terminating bound orbits: r starts in \(\left( 0,r_{m}\right] \) for \(0<r_{m} < \infty \) and it falls into the singularity at \(r=0\).

  • Terminating escape orbits: r comes from \(\infty \) and falls into the singularity at \(r=0\).

The polynomial \({\widetilde{P}}(r)\) can have up to four real zeros, which together with the asymptotic behaviour could give rise to ten different orbit configurations. It will however turn out that \(P(r) \rightarrow +\infty \) for \( r\rightarrow \infty \) always corresponds to an even number of roots whereas \(P(r) \rightarrow -\infty \) corresponds to an odd number. Therefore, only the following five different orbit configurations are possible:

  • In region (0), \({\widetilde{P}}(r)\) has no positive real roots and \({\widetilde{P}}(r)>0\) for \(r\ge 0\). The possible orbit types are terminating escape orbits.

  • In region (1), \({\widetilde{P}}(r)\) has one positive real root \( r_{1}\) with \({\widetilde{P}}(r)\ge 0\) for \(0\le r\le r_{1}\) with possible orbit types being terminating bound orbits.

  • In region (2), \({\widetilde{P}}(r)\) has two positive real roots \(r_{1}\le r_{2}\) with \({\widetilde{P}}(r)\ge 0\) for \(0\le r\le r_{1}\) and \( r_{2}\le r\) with possible orbit types being flyby and terminating bound orbits.

  • In region (3), \({\widetilde{P}}(r)\) has three positive real zeros \(r_{1}\le r_{2}\le r_{3}\) with \({\widetilde{P}}(r)\ge 0\) for \(0\le r\le r_{1}\) and \(r_{2}\le r\le r_{3}\) with possible orbit types being bound and terminating bound orbits.

  • In region (4), \({\widetilde{P}}(r)\) has four positive real zeros \(r_{1}\le r_{2}\le r_{3}\le r_{4}\) with \({\widetilde{P}}(r)\ge 0\) for \(0\le r\le r_{1}\) and \(r_{2}\le r\le r_{3}\) and \(r_{4}\le r\) with possible orbit types being bound, terminating bound and flyby orbits .

Note that the regions (1) and (3) only appear for \(\Lambda <0\), whereas the regions (0), (2), and (4) are only possible if \(\Lambda >0\) in the case of massive particles (\(\epsilon =1\)), and if \((1/b^2)+\Lambda >0\) in the case of massless particles (\(\epsilon =0\)). Here \(b=L/E\) is the impact parameter.

3.2 Charged cases

For the charged case, the presence of the logarithmic function spoils the polynomial form of P(r). Note that although P(r) always has a double zero at \(r=0\), for the charged case, the function \({\widetilde{P}}(r)\) approaches \(-\infty \) for \(r\rightarrow 0\). Therefore, based on Bolzano’s theorem [127], since \({\widetilde{P}}(r)\) is continuous in \( (0,+\infty )\), if \(\Lambda >0\) there is always a real and positive root.

Since the logarithmic charge term changes the asymptotic behaviour of \( {\widetilde{P}}(r)\) for \(r\rightarrow 0\) from a positive finite value to negative infinity, this suggests the presence of a small additional root \(r_{ \mathrm {extra}}\) as compared to the uncharged case. If \(r_{extra}<r_{-}\) (where \(r_{-}\) is inner horizon), then a flyby or bound orbit will be reflected by the singularity and cross the horizons multiple times and entering in a new copy of the universe. These orbits are called two-world escape orbit and many-world bound orbit, respectively (see [105, 118, 128, 129] for more detail).

Numerical calculations show that the following orbits can be identified for the charged cases

  • Flyby orbits: r starts from \(\infty \), then approaches a periapsis \( r=r_{\min }\) for \(r_{\min }>r_{+}\), and goes back to \(\infty \).

  • Bound orbits: r changes between two boundary values \(r_{1}\) \(\le r\) \(\le \) \(r_{2}\), with \(r_{1},r_{2}>r_{+}\) or \(0<r_{1},r_{2}<r_{-}\) where \( r_{+}\) is event horizon and \(r_-\) the inner horizon.

  • Many-world bound orbits: r changes between to boundary values \(r_{1}\) \(\le r\) \(\le \) \(r_{2}\) , with \(0<r_{1}\) \(<r_{-}\) and \(r_{2}>r_{+}\).

  • Two-world escape orbits: with \(r>r_{1}\) and \(0<r_{1}\) \(<r_{-}\).

For the charged case, we can find six regions with different number of roots which leads to following orbits:

  • In region (0), \({\widetilde{P}}(r)\) has no real positive root with \({\widetilde{P}}(r)\) \(<0\) and the motion is impossible .

  • In region (1), \({\widetilde{P}}(r)\) has one real positive root \( 0<r_{1}<r_{-}\) with \({\widetilde{P}}(r)\ge 0\) for \(\,r_{1}\le r\) and possible orbit types are two-world escape orbits.

  • In region (2), \({\widetilde{P}}(r)\) has two real positive roots \(0<r_{1}<r_{-}<r_{+}<r_{2}\) with \({\widetilde{P}}(r)\ge 0\) for \(r_{1}\le r\le r_{2}\) and possible orbit types are many-world bound orbits.

  • In region (3), \({\widetilde{P}}(r)\) has three real positive roots \(0<r_{1}<r_{-}<r_{+}<r_{2}\le r_{3}\) with \({\widetilde{P}}(r)\ge 0\) for \(r_{1}\le r\le r_{2}\) and \(r_{3}\le r\) with possible orbit types being many-world bound and flyby orbits.

  • In region (4), \({\widetilde{P}}(r)\) has four real positive roots \(0<r_{1}<r_{-}<r_{+}<r_{2}\le r_{3}\le r_{4}\) with \({\widetilde{P}} (r)\ge 0\) for \(r_{1}\le r\le r_{2}\) and \(r_{3}\le r\le r_{4}\) with possible orbit types being bound and many-world bound orbits.

  • In region (5), \({\widetilde{P}}(r)\) has five real positive roots \(0<r_{1}<r_{-}<r_{+}<r_{2}\le r_{3}\le r_{4}\le r_{5}\) with \( {\widetilde{P}}(r)\ge 0\) for \(r_{1}\le r\le r_{2}\) and \(r_{3}\le r\le r_{4}\) and \(r_{3}\le r\), with possible orbit types being many-world bound, bound and flyby orbits.

Note that region (2) is the only possible region for the charged BTZ black holes. Region (4) only appears for particle motion in the charged massive BTZ black hole with \(\Lambda <0\), \(m^{\prime }>0\).

4 BTZ black holes

In this section, we study the test particle’s motion in charged BTZ black holes. Although the geodesic motion of a particle around the uncharged black hole has been studied before [124, 125], for the sake of comparison, we will first derive the equations of motion for the uncharged BTZ black holes before proceeding to the charged case.

4.1 Uncharged BTZ black holes

4.1.1 General classification of motion

Substituting Eq. (2) into Eq. (13), we obtain

$$\begin{aligned} V_{\mathrm {eff}}&= -\Lambda \epsilon r^2 - (m_0\epsilon +\Lambda L^2) - \frac{m_0L^2}{r^2}\,, \end{aligned}$$
(18)
$$\begin{aligned} P(r)&=\left( \frac{\epsilon \Lambda }{L^{2}}\right) r^{6}+\left( \frac{ E^{2}}{L^{2}}+\Lambda +\frac{\epsilon m_{0}}{L^{2}}\right) r^{4}+ m_{0} r^{2}. \end{aligned}$$
(19)

Let us first consider the case of massive test-particles with \(\epsilon =1\). Upon inspection of the effective potential (18) and keeping in mind that \(\Lambda <0\), we see that it diverges to infinity for \( r\rightarrow \infty \) and to minus infinity for \(r\rightarrow 0\). Its derivation with respect to r is

$$\begin{aligned} \frac{dV_{\mathrm {eff}}}{dr}=\frac{2m_{0}L^{2}}{r^{3}}-2\Lambda \epsilon r, \end{aligned}$$

which is always positive for \(r>0\). Referring to (12), this implies that all massive particle trajectories have some outer turning point \(r_{0}>0\) and eventually have to cross the black hole horizon at \( r_{+}=m_{0}/\sqrt{-\Lambda }\).

The same result can be inferred from (19). By Descartes’ rule of signs, P(r) posses at most one positive real zero \( r_{0} \), and \(P(r)\rightarrow -\infty \) for \(r\rightarrow \infty \). Massive test particles are therefore bound to the region \(0\le r\le r_{0}\).

Now let us turn to massless particles, \(\epsilon =0\). Here, it is convenient to rescale the affine parameter \(\lambda \) such that the equation of motion simplifies to

$$\begin{aligned} \left( \frac{dr}{d\lambda }\right) ^{2}=1-{\hat{V}}_{\mathrm {eff}}=1-(-\Lambda r^{2}-m_{0})\frac{b^{2}}{r^{2}}, \end{aligned}$$

where \(b=L/E\) is the impact parameter. Then, the effective potential \({\hat{V}}_{\mathrm {eff}}\) approaches \(-\Lambda b^{2}\) from below for \(r\rightarrow \infty \) and diverges to minus infinity for \(r\rightarrow 0\). As for derivation with respect to r of \({\hat{V}}_{\mathrm {eff}}\),

$$\begin{aligned} \frac{d{\hat{V}}_{\mathrm {eff}}}{dr}=\frac{2m_{0}b^{2}}{r^{3}}, \end{aligned}$$

is positive for \(r>0\), this implies that a photon may reach infinity if \( -\Lambda b^{2}<1\), and otherwise it is bounded by an outer turning point \( r_{0}\). From this analysis, it is also clear that circular orbits cannot exist.

4.1.2 Analytic solution of geodesic equations

We use the substitution \(r^{2}=\frac{1}{u}\) and slightly rewrite Eq. (19) to obtain

$$\begin{aligned} \int _{u_{0}}^{u}\frac{du}{\sqrt{u^{2}+c_{1}u+c_{2}}}=2\sqrt{m_{0}}(\varphi -\varphi _{0})\,, \end{aligned}$$
(20)

where \(c_{1}=\frac{E^{2}}{L^{2}m_{0}}+\frac{\Lambda }{m_{0}}+\frac{\epsilon }{L^{2}}\), \(c_{2}=\frac{\epsilon \Lambda }{m_{0}L^{2}}\). This integral has the solution

$$\begin{aligned}&\ln \left( \frac{c_{1}}{2}+u+\sqrt{u^{2}+c_{1}u+c_{2}}\right) \nonumber \\&\qquad -\ln \left( \frac{c_{1}}{2}+u_{0}+\sqrt{u_{0}^{2}+c_{1}u_{0}+c_{2}}\right) \nonumber \\&\quad =2\sqrt{m_{0}} (\varphi -\varphi _{0})\,. \end{aligned}$$
(21)
Fig. 1
figure 1

Null geodesics for the uncharged BTZ black hole with \(\Lambda =-0.1\) and \(m_{0}=2\). The dash dotted line represents the horizon

Equivalently, we find u as a function of \(\varphi \),

$$\begin{aligned} u(\varphi )= & {} \left( \frac{{L}^{4}{\Lambda }^{2}+2{\Lambda L}^{2}\left( {E}^{2}-\epsilon m_{0}\right) +\left( {E}^{2}+\epsilon m_{0}\right) ^{2}}{8{L}^{4}m_{0}^{2}\xi }\right) \nonumber \\&\times e^{2\sqrt{m_{0}}\left( \varphi _{0}-\varphi \right) }+\frac{1}{2}\,e^{-2\sqrt{m_{0}}\left( \varphi _{0}-\varphi \right) }-\frac{c_{1}}{2}, \end{aligned}$$
(22)

where \(\xi \) is related to \(u_{{0}}\) as

$$\begin{aligned} \xi =\frac{c_{1}}{2}+u_{0}+\sqrt{u_{0}^{2}+c_{1}u_{0}+4c_{2}}. \end{aligned}$$
(23)

Therefore, \(r(\varphi )\) is

$$\begin{aligned} r(\varphi )=\frac{1}{\sqrt{u(\varphi )}}, \end{aligned}$$
(24)

which is valid for both timelike and lightlike geodesics.

All the possible types of null geodesics in BTZ black holes are plotted in Fig. 1. Types of orbits mentioned in these plots have been discussed in Sect. 3.

4.2 Charged BTZ black holes

For the line element (1), the linearly charged solution of BTZ black holes has been obtained as [95]

$$\begin{aligned} \psi _2 \left( r\right) =-\Lambda r^{2}-m_{0}-2q^{2}\ln \left( \frac{r}{r_{0}} \right) , \end{aligned}$$
(25)

so from Eq. (13), P(r) and \(V_{\mathrm {eff}}\) are

$$\begin{aligned} P(r)= & {} \left( \frac{\epsilon \Lambda }{L^{2}}\right) r^{6} \nonumber \\&+\left( \frac{E^{2}}{ L^{2}}+\Lambda +\frac{\epsilon }{L^{2}}\left\{ m_{0}+2q^{2}\ln \left( \frac{r }{r_{0}}\right) \right\} \right) r^{4} \nonumber \\&+\left( 2q^{2}\ln \left( \frac{r}{r_{0} }\right) +m_{0}\right) r^{2}. \end{aligned}$$
(26)
$$\begin{aligned} V_{\mathrm {eff}}= & {} \left( -\Lambda r^{2}-m_{0}-2q^{2}\ln \left( \frac{r}{ r_{0}}\right) \right) \left( \frac{L^{2}}{r^{2}}+\epsilon \right) . \end{aligned}$$
(27)

In contrast to the uncharged case, here we can find circular orbits, which are given by \(\frac{dr}{d\lambda }=0\) and \(\frac{d^{2}r}{d\lambda ^{2}}=0\). These two conditions are equivalent to \(P(r)=0\) and \(\frac{dP}{dr}=0\). We can solve these two equations for the squared energy and angular momentum as

$$\begin{aligned} E^{2}&=-\frac{\left( \Lambda r^{2}+2q^{2}\ln \left( \frac{r}{r_{0}}\right) +m_{0}\right) ^{2}}{2q^{2}\ln \left( \frac{r}{r_{0}}\right) -q^{2}+m_{0}}, \end{aligned}$$
(28)
$$\begin{aligned} L^{2}&=\frac{r^{2}\left( \Lambda r^{2}+q^{2}\right) }{2q^{2}\ln \left( \frac{r}{r_{0}}\right) -q^{2}+m_{0}}, \end{aligned}$$
(29)

for massive particles \((\epsilon =1)\). Let us discuss these two equations. From Eq. (28), it is evident that \(E^{2}>0\) is valid only for a negative denominator and thus

$$\begin{aligned}&2q^{2}\ln \left( \frac{r}{r_{0}}\right) -q^{2}+m_{0}<0 \nonumber \\&\quad \Leftrightarrow 0<r<\Gamma =r_{0}\exp \left( \frac{1}{2}\left[ 1-\frac{m_{0}}{q^{2}} \right] \right) . \end{aligned}$$
(30)

In addition, the numerator of Eq. (29) has to be negative, and therefore, keeping in mind that \(\Lambda <0\)

$$\begin{aligned} r^{2}\left( \Lambda r^{2}+q^{2}\right)<0\quad \Leftrightarrow \quad 0< \frac{|q|}{\sqrt{-\Lambda }}<r . \end{aligned}$$
(31)

Considering the constraints (30) and (31), simultaneously, we find that circular orbits exist if the following relation is satisfied,

$$\begin{aligned} \frac{|q|}{\sqrt{-\Lambda }}<\Gamma \ \ \text {so} \ \ m_{0}<q^{2}\left( \ln \left( -\frac{\Lambda r_{0}^{2}}{q^{2}}\right) +1\right) . \end{aligned}$$
(32)

It is notable that such inequality implies the absence of horizon in the domain \(r>0\). This confirms that for the charged BTZ spacetime circular orbits may only exist around naked singularities, but not for the case of a black hole.

Therefore, similar to the previous section, P(r) has the same number of roots in each combination of E, L and constant parameters. For \(r\ne r_{0}\), we can adjust arbitrary combination of E, L and other constant parameters, in which we obtain two real positive roots for P(r). However, for \(r=r_{0}\), the uncharged case is recovered in which P(r) has always one root. All the possible types of timelike orbits have been plotted in the Fig. 2. Properties of the orbits mentioned in these figures have been presented in Sect. 3.

Fig. 2
figure 2

Timelike geodesic in the charged BTZ black hole with \(\Lambda =-0.1\) , \(m_{0}=2\), \(q=1.5\) and \(r_{0}=1.08\). The dash dotted lines represent the horizons

Fig. 3
figure 3

Null geodesic in the charged BTZ black hole with \(\Lambda =-0.1\), \( m_{0}=2\), \(q=1.5\) and \(r_{0}=1.08\). The dash dotted lines represent the horizons

Let us now turn to massless particles. In this case (\(\epsilon =0\)) circular orbits exist for

$$\begin{aligned} b^{2}&=-\frac{\Gamma ^{2}}{\Lambda \Gamma ^{2}+q^{2}}, \end{aligned}$$
(33)
$$\begin{aligned} r&= \Gamma , \end{aligned}$$
(34)

where \(b=L/E\) is the impact parameter. Due to the fact that \(\Gamma \) is real, in order to have physically acceptable motion for the massless particles, the denominator of Eq. (33) must be negative. Therefore, the following constraint must hold, keeping again in mind that \( \Lambda <0\),

$$\begin{aligned} \frac{|q|}{\sqrt{-\Lambda }}<\Gamma . \end{aligned}$$
(35)

This results in the same inequality (32) as for the case of massive particles which indicates absence of an event horizon.

The many-world bound orbit of null geodesics are shown in Fig. 3.

5 Uncharged Massive BTZ black hole

5.1 General classification of motion

The metric function for the BTZ black holes in massive gravity is presented in Eqs. (7) and (8) [90]. By substituting Eqs. (7) and (8) into Eqs. (13) and (14), we have

$$\begin{aligned} V_{\mathrm {eff}}&=(-\Lambda r^{2}-m_{0}+m^{\prime }r)\left( \frac{L^{2}}{ r^{2}}+\epsilon \right) , \end{aligned}$$
(36)
$$\begin{aligned} P(r)&=\left( \frac{\epsilon \Lambda }{L^{2}}\right) r^{6}-\left( \frac{ \epsilon m^{\prime }}{L^{2}}\right) r^{5} \nonumber \\&\quad +\left( \frac{E^{2}}{L^{2}}+\Lambda +\frac{\epsilon m_{0}}{L^{2}}\right) r^{4}- m^{\prime } r^{3}+ m_{0} r^{2}. \end{aligned}$$
(37)

Let us first discuss massive particles (\(\epsilon =1\)). From the form of P(r) , according to the Descartes rule of signs, it is clear that for \( \Lambda <0\), there may be one or three positive real zeros, whereas for \( \Lambda >0\) there are four, two, or no positive real zeros, depending on the values of E and L as well as the sign of \(m^{\prime }\). This points to the existence of circular orbits, which are given by \(\frac{dr}{d\lambda }=0\) and \(\frac{d^{2}r}{d\lambda ^{2}}=0\). These two conditions are equivalent to \(P(r)=0\) and \(\frac{dP}{dr}=0\). We may solve these two equations for \(E^{2}\) and \(L^{2}\), which gives for massive particles \((\epsilon =1)\)

$$\begin{aligned} E^{2}&=-\frac{2\left( \Lambda r^{2}+m_{0}-m^{\prime }r\right) ^{2}}{ 2m_{0}-m^{\prime }r}, \end{aligned}$$
(38)
$$\begin{aligned} L^{2}&=-\frac{\left( m^{\prime }-2\Lambda r\right) r^{3}}{2m_{0}-m^{\prime }r}. \end{aligned}$$
(39)

It is clear from the expression for \(E^{2}\) that circular orbits can only exist if \(2m_{0}-m^{\prime }r<0\), which is for \(r>0\) only possible if \( m^{\prime }>0\), which then gives us \(r>2m_{0}/m^{\prime }\). The equation (39) for L then implies that \(m^{\prime }-2\Lambda r >0\) has to hold, which is automatically fulfilled for \(\Lambda <0\). For \( \Lambda >0 \) we find \(r<m^{\prime }/(2\Lambda )\).

The circular orbit is stable if its radius corresponds to a maximum of P(r), i.e. if the second derivative of P(r) is negative. The second derivative of P(r) together with Eqs.  (38) and (39) reads

$$\begin{aligned} \frac{d^{2}P}{dr^{2}}=2\frac{2[3\Lambda m^{\prime 2}+(8\Lambda m_{0}+(m^{\prime 2})r-3m_{0}m^{\prime }]}{m^{\prime }-2\Lambda r}. \end{aligned}$$

This expression can be solved for the radius of the circular orbit,

$$\begin{aligned} r_c&= \frac{8\Lambda m_0 + (m^{\prime 2 }\pm \sqrt{[(m^{\prime 2}-16\Lambda m_0][(m^{\prime 2}-4\Lambda m_0]}}{6\Lambda m^{\prime }}. \end{aligned}$$
(40)

For \(\Lambda <0\), we find one positive zero (for the negative sign before the root in (40)) with stable orbits for larger radii, which implies that an innermost stable circular orbit (ISCO) exists. The radius of this ISCO approaches \(3m_0/m^{\prime }\) for small \(\Lambda \) and \( 8m_0/3m^{\prime }\) for large negative \(\Lambda \). We plotted the ISCO for fixed \(m_0\) in Fig. 4.

Interestingly, for \(\Lambda >0\) stable orbits only exist if \(r_c\) from (40) is real, that is, if and only if \((m^{\prime 2}>16\Lambda m_0\). Note that \((m^{\prime 2}<4 \Lambda m_0\) can be rewritten as \(m^{\prime }/(2\Lambda )<2m_0/m^{\prime }\), which implies that no circular orbits exist according to our discussion of (38) and (39) above. Summarized, for \(m^{\prime }>0\) and \(\Lambda >0\) circular orbits exist for \(r \in [2m_0/m^{\prime },m^{\prime }/(2\Lambda )]\), and may be stable if \(m^{\prime }/(2\Lambda )>8m_0/m^{\prime }\). Then the circular orbits are stable in between the two radii given by (40) and we have an innermost and an outermost stable circular orbit. We show the important radii for the case \(m^{\prime }>0\), \(\Lambda >0\) in Fig. 5.

Fig. 4
figure 4

Innermost stable circular orbit in the uncharged BTZ black hole in massive gravity with \(m_0=1\). For the left plot we fixed \(m^{\prime }=1\), which implies that circular orbits exist for \(r>2\), and are stable above the solid red line. In the right plot \(\Lambda =-0.01\). Circular orbits exist above the dashed blue line and are stable above the solid red line. The dotted black line indicates the horizon

Fig. 5
figure 5

Regions of stable and unstable circular orbits in the uncharged massive BTZ black hole with \(m_0=1\) and positive \(\Lambda \). The dashed blue lines enclose the region of circular orbits, and the solid red lines indicate the innermost and outermost stable circular orbit. The dotted black lines indicate the horizons. On the left \(m^{\prime }=1\), on the right \( \Lambda =0.04\)

Now let us turn to massless particles \((\epsilon =0)\), where we introduce \(b=L/E\) as the impact parameter. We then see that the leading coefficient of P(r) may change its sign for \(1/b^{2}=-\Lambda \). If \( 1/b^{2}<-\Lambda \) the polynomial P(r) diverges to \(-\infty \) for \( r\rightarrow \infty \) which implies that the photon may not reach infinity. Moreover, P(r) may have at most one positive real zero (\(r_{0}\)), and therefore, it is bound in a region \(0\le r\le r_{0}\). On the other hand, if \(1/b^{2}>-\Lambda \), then P(r) is positive for \(r\rightarrow \infty \), and one finds the photon may reach infinity, and P(r) has two or no positive real zeros. This again points to the existence of a circular orbit for this case. If we solve the two conditions \(P(r)=0\) and \(\frac{dP}{dr}=0\) for the impact parameter b and the radius of the circular orbit \(r_{c}\), we find

$$\begin{aligned} b^{2}= & {} -\frac{4m_{0}}{4\Lambda m_{0}-m^{\prime 2}}, \end{aligned}$$
(41)
$$\begin{aligned} r_{c}= & {} \frac{2m_{0}}{m^{\prime }}. \end{aligned}$$
(42)

From Eq. (41) for b we infer that \(4\Lambda m_{0}-m^{\prime 2}<0\), or equivalently \(\Lambda <m^{\prime 2}/(4m_{0})\) is necessary for the existence of circular orbits. This circular orbit is stable if it corresponds to a maximum of the polynomial P(r). For its second derivative with respect to r, we find

$$\begin{aligned} \frac{d^{2}P}{dr^{2}}=2m_{0}>0, \end{aligned}$$
(43)

which implies that the circular photon orbit is unstable. Interestingly, the radius of the photon orbit only depends on the ratio \(m_{0}/m^{\prime }\), but not on \(\Lambda \) or the individual parameters \(m_{0}\) and \(m^{\prime }\) . For \(\Lambda \) this is to be expected, as it does not enter as an independent parameter in the equation of motion just as in four dimensional Schwarzschild-de Sitter spacetime, but this seems surprising for the independent parameters \(m_{0}\) and \(m^{\prime }\). However, we correctly recover that the circular orbit vanishes (its radius shifts to infinity) for \(m^{\prime }\rightarrow 0\).

The results of Eqs. (38), (39) and ( 41) for both massive and massless particles are given in Figs. 6, 7 and 8. Additional properties of these figures have been presented in Sect. 3.

Fig. 6
figure 6

The regions of different types of geodesic motions for massive particles in the uncharged massive BTZ black hole with \(m_{0}=2\) and \( m^{\prime }=2.54\) for positive \(\Lambda \). The numbers in parentheses indicate the number of real positive roots (vertical dashed line indicates the location of divergency)

Fig. 7
figure 7

The regions of different geodesic motions for massive particles in the uncharged BTZ black hole in massive gravity with \(m_{0}=2\) and \( m^{\prime }=2.54\) for negative \(\Lambda \). The numbers in parentheses indicate the number of real positive roots

5.2 Analytic solution of geodesic equations

Here, we discuss the analytical solution of Eq. (13). The function P(r) in given in Eq. (37) is a polynomial of degree 6, which can be written in the following form

$$\begin{aligned} \left( \frac{dr}{d\varphi }\right) ^{2}=\sum \limits _{i=2}^{6}a_{i}r^{i}=P(r) \end{aligned}$$
(44)

where the coefficients \(a_{i}\) are

$$\begin{aligned} a_{6}= & {} \frac{\epsilon \Lambda }{L^{2}},\quad a_{5}=-\frac{\epsilon m^{\prime }}{L^{2}},\quad a_{4}=\frac{E^{2}}{L^{2}}+\Lambda +\frac{ \epsilon m_{0}}{L^{2}}, \nonumber \\ a_{3}= & {} -m^{\prime },\quad a_{2}=m_{0}. \end{aligned}$$
(45)

By substitution \(r=u^{-1}+r_{M}\) into Eq. (13) , where \(r_{M}\) is a root of P(r), for instance \(r_{M}=0\), we find

$$\begin{aligned}&\left( \frac{du}{d\varphi }\right) ^{2}=P(u)=\sum \limits _{j=-2}^{3}b_{j}u^{j}, \nonumber \\&\quad \left. b_{j}=\frac{1}{ (4-j)!}\frac{d^{(4-j)}P(r)}{dr^{^{(4-j)}}}\right| _{r=r_{M}}, \end{aligned}$$
(46)

and the coefficients \(b_{j}\) can be calculated as

$$\begin{aligned} b_{3}= & {} 6\,{\frac{\Lambda \,\epsilon \,r_{M}^{5}}{{L}^{2}}}-5\,{\frac{ m^{\prime }\,\epsilon \,r_{M}^{4}}{{L}^{2}}} \nonumber \\&+\,4\,{\frac{\left( {L}^{2}\Lambda +{E}^{2}+\epsilon \,m_{0}\right) r_{M}^{3}}{{L}^{2}}}-3\,m^{\prime }\,r_{M}^{2}+2\,m_{0}\,r_{M}, \nonumber \\ b_{2}= & {} 15\,{\frac{\Lambda \,\epsilon \,r_{M}^{4}}{{L}^{2}}}-10\,{\frac{ m^{\prime }\,\epsilon \,r_{M}^{3}}{{L}^{2}}} \nonumber \\&+\,6\,{\frac{\left( {L}^{2}\Lambda +{E}^{2}+\epsilon \,m_{0}\right) r_{M}^{2}}{{L}^{2}}}-3\,m^{\prime }\,r_{M}+m_{0}, \nonumber \\ b_{1}= & {} 20\,{\frac{\Lambda \,\epsilon r_{M}^{3}}{{L}^{2}}}-10\,{\frac{ m^{\prime }\,\epsilon \,r_{M}^{2}}{{L}^{2}}} \nonumber \\&+\,4\,{\frac{\left( {L}^{2}\Lambda +{E}^{2}+\epsilon \,m_{0}\right) r_{M}}{{L}^{2}}}-m^{\prime }, \nonumber \\ b_{0}= & {} 15\,{\frac{\Lambda \,\epsilon \,r_{M}^{2}}{{L}^{2}}}-5\,{\frac{ m^{\prime }\,\epsilon \,r_{M}}{{L}^{2}}}+{\frac{{L}^{2}\Lambda +{E} ^{2}+\epsilon \,m_{0}}{{L}^{2}},} \nonumber \\ b_{-1}= & {} {6\,{\frac{\Lambda \,\epsilon \,r_{M}}{{L}^{2}}}-{\frac{m^{\prime }\,\epsilon }{{L}^{2}}},}\quad b_{-2}={\frac{\Lambda \,\epsilon }{{L}^{2}}}. \end{aligned}$$
(47)

Now, we are in a position to obtain analytic solutions of Eq. (46 ) for massless and massive particles.

Fig. 8
figure 8

The regions of different geodesic motions for massless particles in the uncharged BTZ black hole in massive gravity with \(m_{0}=2\) and \( m^{\prime }=2.54\). The numbers in parentheses indicate the number of real positive roots

Fig. 9
figure 9

Timelike geodesics in the uncharged BTZ black hole in massive gravity with \(\Lambda =0.01\), \(m^{\prime }=2.54\) and \(m_{0}=2\). The dashed dot line represents horizons

Fig. 10
figure 10

Null geodesics in the uncharged BTZ black hole in massive gravity with \(\Lambda =0.01\), \(m^{\prime }=2.54\) and \(m_{0}=2\). The dashed dot line represents horizon

Fig. 11
figure 11

Timelike geodesics in the uncharged BTZ black hole in massive gravity with \(\Lambda =-0.01\), \(m^{\prime }=2.54\) and \(m_{0}=2\). The dashed dot line represents horizons

5.2.1 Null geodesics \((\epsilon =0)\)

By substituting \((\epsilon =0)\) into Eq. (46), it transforms to the following

$$\begin{aligned} \left( \frac{du}{d\varphi }\right) ^{2}=P(u)=\sum \limits _{j=0}^{3}\alpha _{j}u^{j}, \end{aligned}$$
(48)

where \(\alpha _{j}=\left. b_{j}\right| _{\epsilon =0}\).

The roots of Eq. (48) lead to an elliptic type differential equation. Another substitution \(u=\frac{1}{\alpha _{3}}(4y- \frac{\alpha _{2}}{3})\) transforms \(P_{3}(u)\) into the following form

$$\begin{aligned} \left( \frac{dy}{d\varphi }\right) ^{2}=4y^{3}-g_{2}y-g_{3}, \end{aligned}$$
(49)

in which the coefficients \(g_{2}\) and \(g_{3}\) are

$$\begin{aligned} g_{2}= & {} \frac{1}{16}\left( \frac{4}{3}\alpha _{2}^{2}-4\alpha _{1} \alpha _{3}\right) , \nonumber \\ g_{3}= & {} \frac{1}{16}\left( \frac{1}{3}\alpha _{1}\alpha _{2}\alpha _{3}- \frac{2}{27}\alpha _{2}^{3}-\alpha _{0}\alpha _{3}^{2}\right) . \end{aligned}$$
(50)
Fig. 12
figure 12

The regions of different geodesic motion for massless particle in the charged BTZ black hole in massive gravity with \(m_0=1\), \(m^{\prime }=1\) and \(r_0=1\). On the left \(q^2=1\), on the right \(\Lambda =-0.1\). The red dashed line corresponds to \(b^{-2}+\Lambda =0\), where the asymptotic for \(r \rightarrow \infty \) switches sign. The numbers in parentheses indicate the number of real positive roots

Fig. 13
figure 13

The regions of different geodesic motion for massive particle in the charged BTZ black hole in massive gravity with \(m_0=1\), \(m^{\prime }=1\), \(r_0=1\), and \(q^2=1\). On the left \(\Lambda =-0.01\) and on the right \( \Lambda =0.001\). The numbers in parentheses indicate the number of real positive roots

It is known that the solution of Eq. (49) is the Weierstrass function \(\wp (\varphi )\) in the following form

$$\begin{aligned} y=\wp (\varphi -\varphi _{in}), \end{aligned}$$
(51)

where

$$\begin{aligned} \varphi _{in}= & {} \varphi _{0}+\int \nolimits _{y_{0}}^{\infty }\frac{dy}{\sqrt{ 4y^{3}-g_{2}y-g_{3}}}, \nonumber \\ y_{0}= & {} \frac{1}{4}\left( \frac{\alpha _{3} }{r_{0}-r_{M}}-\frac{\alpha _{2}}{3}\right) , \end{aligned}$$
(52)

in which \(r_{0}\) and \(\varphi _{0}\) are the initial values of the differential equation. Therefore, the general solution of Eq. (48) is

$$\begin{aligned} r\left( \varphi \right) =\frac{\alpha _{3}}{4\wp (\varphi -\varphi _{in})- \frac{\alpha _{2}}{3}}+r_{M}, \end{aligned}$$
(53)

which is in agreement with the result of Ref. [107].

Fig. 14
figure 14

Timelike geodesics in the charged BTZ black hole in massive gravity with \(\Lambda =10^{-2}\), \(m_{0}=2\), \(m^{\prime }=2.54\), \(q=1.5\) and \( r_{0}=0.95\). The dashed dot lines represent horizons

Fig. 15
figure 15

Null geodesics in the charged BTZ black hole in massive gravity with \(\Lambda =10^{-5}\), \(m_{0}=2\), \(m^{\prime }=3.02\), \(q=0.6\) and \( r_{0}=0.93\). The dashed dot lines represent horizons

5.2.2 Timelike geodesics \((\epsilon =1)\)

By substituting \((\epsilon =1)\) into Eq. (46), it can be transformed to

$$\begin{aligned}&\left( u\frac{du}{d\varphi }\right) ^{2}=\sum \limits _{j=0}^{5}\beta _{j}u^{j}=R(u), \nonumber \\&\quad \left. \beta _{j}=\frac{1}{(6-j)!}\frac{d^{(6-j)}P(r)}{ dr^{^{(6-j)}}}\right| _{r=r_{M}}, \end{aligned}$$
(54)

where \(\beta _{j}=\left. b_{j+2}\right| _{\epsilon =1}\). Equation (54) is hyperelliptic type and its analytical solution is given by the Kleinian sigma function [101, 106, 130]

$$\begin{aligned} u\left( \varphi \right) =-\frac{\sigma _{1}}{\sigma _{2}}\left( \varphi _{\sigma }\right) , \end{aligned}$$
(55)

in which \(\sigma _{i}\left( z\right) \) denotes the derivative of the Kleinian sigma function with respect to the \(i^{th}\) component of z

$$\begin{aligned} \sigma \left( z\right) =Ce^{-\frac{1}{2}z^{t}\eta \omega ^{-1}z}\vartheta \left[ g,h\right] \left( \left( 2\omega \right) ^{-1}z;\tau \right) , \end{aligned}$$
(56)

where \(\left( 2\omega ,2\omega ^{^{\prime }}\right) \) is the period matrix, \( \left( 2\eta ,2\eta ^{^{\prime }}\right) \) is the period matrix of the second kind, C is a constant that can be given explicitly but does not enter (55) and \(\tau =\omega ^{-1}\omega ^{^{\prime }}\). The theta function is defined as follow

$$\begin{aligned} \vartheta \left[ g,h\right] \left( z;\tau \right) =\sum _{m\in {\mathbb {Z}}^{2}}e^{i\pi \left( m+g\right) ^{t}\left( \tau \left( m+g\right) +2z+2h\right) }, \end{aligned}$$
(57)

in which gh are two dimensional vectors related to the vector of Riemann constants. They are defined as \(g=(1/2,1/2)\), \(h=(0,1/2)\). The argument \( \varphi _{\sigma }\) in (55) is defined as

$$\begin{aligned} \varphi _{\sigma }=\left( \begin{array}{c} f\left( \varphi -\varphi _{in}\right) \\ \varphi -\varphi _{in} \end{array} \right) , \end{aligned}$$
(58)

where f is given by the condition \(\sigma \left( \varphi _{\sigma }\right) =0\). The constant \(\varphi _{in}\) reads

$$\begin{aligned} \varphi _{in}=\varphi _{0}+\int \nolimits _{u_{0}}^{\infty }\frac{udu}{\sqrt{ R(u)}}, \end{aligned}$$
(59)

and only depends on the initial values \(\varphi _{0}\) and \(u_{0}\). The general solution of the radial coordinate r is given by

$$\begin{aligned} r\left( \varphi \right) =-\frac{\sigma _{2}}{\sigma _{1}}\left( \varphi _{\sigma }\right) +r_{M}. \end{aligned}$$
(60)

Equations (53) and (60) completely describe the motion of massive and massless particles in this spacetime. All the possible type of orbits in BTZ black holes of massive gravity have been plotted in Figs. 9, 10 and 11. Types of orbits mentioned here have been introduced in the total classification Sect. 3.

6 Charged black holes in massive gravity

6.1 General classification of motion

In the charged black holes in massive gravity, the metric function \(\psi (r)\) is obtained as [90]

$$\begin{aligned} \psi (r)=-\Lambda r^{2}-m_{0}+m^{\prime }r-2q^{2}\ln \left( \frac{r}{r_{0}} \right) . \end{aligned}$$
(61)

From Eq. (13), we then find P(r) as

$$\begin{aligned} P(r)= & {} \left( \frac{\epsilon \Lambda }{L^{2}}\right) r^{6}-\left( \frac{\epsilon m^{\prime }}{L^{2}}\right) r^{5} \nonumber \\&+\left( \frac{E^{2}}{L^{2}}+\Lambda +\frac{\epsilon }{L^{2}}\left\{ m_{0}+2q^{2}\ln \left( \frac{r}{r_{0}}\right) \right\} \right) r^{4} \nonumber \\&- m^{\prime } r^{3}+\left\{ m_{0} + 2q^{2}\ln \left( \frac{r}{ r_{0}}\right) \right\} r^{2}. \end{aligned}$$
(62)

Again, we can look for circular orbits by solving \(P(r)=0\) and \(\frac{dP(r)}{ dr}=0\) for the squared energy \(E^{2}\) and angular momentum \(L^{2}\). For massive particles \((\epsilon =1)\) we find

$$\begin{aligned} E^{2}&=\frac{2\left( \Lambda r^{2}+2q^{2}\ln \left( \frac{r}{r_{0}}\right) +m_{0}-m^{\prime }r\right) ^{2}}{m^{\prime }r-2m_{0}+2q^{2}-4q^{2}\ln \left( \frac{r}{r_{0}}\right) }, \end{aligned}$$
(63)
$$\begin{aligned} L^{2}&=-\frac{r^{2}\left( 2\Lambda r^{2}-m^{\prime }r+2q^{2}\right) }{ m^{\prime }r-2m_{0}+2q^{2}-4q^{2}\ln \left( \frac{r}{r_{0}}\right) }\,. \end{aligned}$$
(64)

As expected, for \(q=0\) Eqs. (63) and (64) reduce to Eqs. (38) and (39), respectively, and for \(m^{\prime }=0\) they reduce to Eqs. (28) and (29). From equation (63) we infer the inequality

$$\begin{aligned} R := m^{\prime }r-2m_{0}+2q^{2}-4q^{2}\ln \left( \frac{r}{r_{0}}\right) >0. \end{aligned}$$
(65)

For \(m^{\prime }>0\) the function R diverges to infinity at \(r=0\) and \( r=\infty \), and has only a single extrema, a minimum at \(r=4q^2/m^{\prime }\). It may have two zeros \(\Gamma _{\mathrm {m},0} < \Gamma _{\mathrm {m},-1}\) given by

$$\begin{aligned} \Gamma _{\mathrm {m},0}&:= r_0 \exp \left( \frac{1}{2} - \frac{m_0}{2q^2} - \mathrm {W}\left( \frac{-m^{\prime }r_0}{4q^2} e^{\frac{1}{2}-\frac{m_0}{2q^2} } \right) \right) \,, \end{aligned}$$
(66)
$$\begin{aligned} \Gamma _{\mathrm {m},-1}&:= r_0 \exp \left( \frac{1}{2} - \frac{m_0}{2q^2} - \mathrm {W}\left( -1, \frac{-m^{\prime }r_0}{4q^2} e^{\frac{1}{2}-\frac{m_0}{ 2q^2}} \right) \right) \,. \end{aligned}$$
(67)

If \(m^{\prime }\le 0\) the denominator R is monotonically decreasing on \( (0,\infty )\) and introduces therefore an upper bound on r given by \( r<\Gamma _{\mathrm {m},0}\). We can determine the relative position of the zeros of R and the horizons: let \(r_H\) be one of the horizons, i.e. \( \psi _4(r_H)=0\). Then we find from this

$$\begin{aligned} \psi _4(r_H)&= 0 \nonumber \\&\Rightarrow \quad \ln \left( \frac{r_H}{r_0}\right) = \frac{1}{2q^2}(-\Lambda r_H^2 - m_0 + m^{\prime }r_H) \nonumber \\&\Rightarrow \quad R(r_H) = -m^{\prime }r_H + 2q^2 +2\Lambda r_H^2 = - \frac{1}{r_H} \frac{d\psi _4}{dr}\bigg |_{r=r_H}\,. \end{aligned}$$
(68)

Therefore, R is positive at the inner horizon and negative at the event horizon. For \(\Lambda >0\) a cosmological horizon also exists, where R is positive again. We conclude that \(r_{H,\mathrm {inner}}<\Gamma _{\mathrm {m},0}<r_{H,\mathrm {event}}<\Gamma _{\mathrm {m},-1}<r_{H,\mathrm {cosmo}}\), if the respective horizons and zeros of R exist.

We can directly infer that for \(\Lambda <0\) and \(m^{\prime }\le 0\) circular orbits outside of a black hole cannot exist, as for charged BTZ black holes discussed in Sect. 4.2 and uncharged BTZ black holes in massive gravity discussed in Sect. 5.

In addition, the nominator in (64) needs to be negative. For \(\Lambda < 0\) this implies a lower bound on r given by

$$\begin{aligned} r\ge \frac{m^{\prime }- \sqrt{(m^{\prime 2}-16\Lambda q^2}}{4\Lambda }. \end{aligned}$$
(69)

Note that this bound is the minimum of \(\psi _4\) and, therefore, is smaller than the event horizon. For \(\Lambda <0\) and \(m^{\prime }>0\) circular orbits around a black hole exist at \(r>\Gamma _{\mathrm {m},-1}\), as in the case of uncharged BTZ black holes in massive gravity discussed in Sect. 5.

From the condition that the nominator in (64) has to be negative, in the case \(\Lambda >0\) circular orbits are not allowed for \(m^{\prime }\le 0\), and have only a limited range for \( m^{\prime }>0\) given by

$$\begin{aligned} \frac{m^{\prime }- \sqrt{(m^{\prime 2}-16\Lambda q^2}}{4\Lambda } \le r \le \frac{m^{\prime }+ \sqrt{(m^{\prime 2}-16\Lambda q^2}}{4\Lambda }. \end{aligned}$$
(70)

Note that the two bounds correspond to the minimum and the maximum of the metric function \(\psi _{4}\). Therefore, it only remains to show that \(\Gamma _{\mathrm {m},-1}\) is smaller than the upper bound in (70): as \(r_{*}:=\Gamma _{\mathrm {m},-1}\) is a zero of R we find

$$\begin{aligned} \ln \left( \frac{r_{*}}{r_{0}}\right)&=\frac{2q^{2}+m^{\prime }r_{*}-2m_{0}}{4q^{2}} \nonumber \\&\Rightarrow \quad \psi _{4}(r_{*}) =-\Lambda r_{*}^{2}-m_{0}+m^{\prime }r_{*} \nonumber \\&\quad +\frac{1}{2}(2m_{0}-2q^{2}-m^{\prime }r_{*}) \end{aligned}$$
(71)
$$\begin{aligned}&=-\frac{1}{2}(2\Lambda r_{*}^{2}-m^{\prime }r_{*}+2q^{2})=\frac{2}{ r_{*}}\frac{d\psi _{4}}{dr}\bigg |_{r=r_{*}}. \end{aligned}$$
(72)

It is clear that \(\psi _4\) is positive at \(r_*=\Gamma _{\mathrm {m},-1}\) which implies that \(r_*\) is smaller than the upper bound in (70). We therefore find circular orbits for \( \Lambda >0\) and \(m^{\prime }>0\) in the range

$$\begin{aligned} \Gamma _{\mathrm {m},-1} < r \le \frac{m^{\prime }+ \sqrt{(m^{\prime 2}-16\Lambda q^2}}{4\Lambda }. \end{aligned}$$
(73)

Let us now turn to massless particles \((\epsilon =0)\). We find for the circular orbits

$$\begin{aligned} b^{2}=\frac{r_{c}^{2}}{\psi _4(r_c)}, \end{aligned}$$
(74)

where \(r_c\) is a solution of \(R=0\), see Eq. (65). Therefore, again for \(m^{\prime }<0\) there is no circular photon orbit outside of a black hole. If \(m^{\prime }>0\) a circular photon orbit exists at \(r=\Gamma _{\mathrm {m},-1}\). In the limit \(m^{\prime }=0\) this reduces to the results of Sect.4.2, and for \(q=0\), the zero of R is \(r=2m_0/m^{\prime }\) as in Sect. 5. The circular photon orbit is always unstable: if we insert \(b^2\) and r into the second derivative of P, we find stability for

$$\begin{aligned} m^{\prime }> \frac{4q^2}{r_0} e^{\frac{m_0}{2q^2} - \frac{3}{2}}\,, \end{aligned}$$
(75)

which is however incompatible with a real \(\Gamma _{\mathrm {m},-1}\).

The results of Eqs. (63), (64) and (74) for both massive and massless particles are given in Figs. 12 and 13. (For more details we refer the reader to Sect. 3). We observe that a variation of q as well as the value and sign of cosmological constant have considerable effect on the geodesic motion of massless particles. More clearly, we find that increasing the cosmological constant or the charge leads to increasing the possibility of two world escape obits rather than flyby and many world bound orbits. According to these figures, for some values of the charge or the cosmological constant, there is no physical motion for large values of b. For massive particles, we see that the possible types of orbits are rather different for \(\Lambda <0\) and \(\Lambda >0\). In the case \(\Lambda <0\), it is impossible to reach radial infinity and bound orbits outside the horizons do not exist. On the other hand, for \( \Lambda >0\) all parameter combinations allow orbits reaching infinity and bound orbits outside the horizons.

6.2 Numerical solution of geodesic equations

Now, we are in a position to discuss both null and timelike geodesics of charged BTZ solutions in massive gravity. Since the metric function is no longer a polynomial function, to our knowledge, there is no analytical solution for the charged cases. Therefore, we use the Runge–Kutta–Fehlberg numerical method. It is worthwhile to mention that we have checked this method for uncharged cases which we have analytical solutions, and the results were the same with high accuracy.

All the possible type of orbits in charged BTZ black holes in massive gravity with a positive cosmological constant are plotted by using the numerical analysis in Figs. 14 and 15 (see also Sect. 3 for more details). We restricted to positive \(\Lambda \) due to the interesting possibility of bound orbital motion.

As already familiar from the case of the charged BTZ solution, for both timelike and null geodesics, we find many-world bound orbits as well as the two-world escape orbits, that cross both the inner and the event horizon. For a positive cosmological constant, it is well known from four-dimensional black holes like the Schwarzschild-de Sitter solution, that flyby orbits can exist that are deflected from the cosmological barrier. This can be seen in Fig. 14b, e, where a massive test particle approaches the black hole rather straight and is reflected back without revolving around the black hole.

7 Summary and conclusions

In this paper, we analyzed the geodesic motion in (un)charged BTZ black hole and its generalization to massive gravity. We provided a complete classification of the possible types of geodesic motion for the four different spacetimes considered here, thereby investigating the effects of the various parameters of the metric on the geodesics. In particular, to our knowledge for the first time, we investigated the geodesic motion in a spacetime whose equation of motion for test particles has a non-polynomial structure. Moreover, we presented analytical solutions for the equations of motion for the uncharged cases, and used a Runge–Kutta–Fehlberg method for the spacetimes of the charged black holes, where an analytical solution could not be found.

While the uncharged BTZ black hole has rather a poor variety of orbital motions, in particular lacking bound orbital motion outside the horizon, the inclusion of a charge introduces a potential barrier within the inner horizon. This barrier reflects particles and light such that they have to cross the inner horizon for a second time, thereby entering another copy of the universe. This behaviour is known from some four-dimensional spacetimes, and the corresponding orbits are called two-world flyby or many-world bound orbits, respectively.

In three-dimensional massive gravity, the black hole solutions show an even richer structure of geodesic motion. First, we are no longer restricted to a negative cosmological constant, which gives rise to new types of geodesics. Second, the massive gravity parameter \(m^{\prime }\) introduced in Eq.  (8) further enriches the structure of geodesic motion. We found that \(m^{\prime }>0\) is a necessary condition for the existence of circular orbits for both uncharged and charged black holes.

For the uncharged BTZ black hole in massive gravity, in addition to the general classification of geodesic motion, we derived the radius of the innermost stable and (for \(\Lambda >0\)) the outermost stable circular orbit. The results are plotted in Figs. 4 and 5. In particular, it was pointed out that stable bound orbital motion outside the event horizon is possible, and we showed an example in Fig. 9. Moreover, we calculated the radius of the unstable photon sphere, that is an important radius for the calculation of the shadow of the black hole.

Finally, for the most general case discussed in this paper, the charged BTZ black hole in massive gravity, we derived the range of radii where circular orbits may exist. Based on this, we presented a complete classification of geodesic motion. However, due to the complicated structure of the equations, we did not determine the innermost stable circular orbit, and leave this point for future research.