Simulation of geodesic trajectory of charged BTZ black holes in massive gravity

In order to classify and understand the spacetime structure, investigation of the geodesic motion of massive and massless particles is a key tool. So the geodesic equation is a central equation of gravitating systems and the subject of geodesics in the black hole dictionary attracted much attention. In this paper, we give a full description of geodesic motions in three-dimensional spacetime. We investigate the geodesics near charged BTZ black holes and then generalize our prescriptions to the case of massive gravity. We show that electric charge is a critical parameter for categorizing the geodesic motions of both lightlike and timelike particles. In addition, we classify the type of geodesics based on the particle properties and geometry of spacetime.

In order to classify and understand the spacetime structure, investigation of the geodesic motion of massive and massless particles is a key tool. So the geodesic equation is a central equation of gravitating systems and the subject of geodesics in the black hole dictionary attracted much attention. In this paper, we give a full description of geodesic motions in three-dimensional spacetime. We investigate the geodesics near charged BTZ black holes and then generalize our prescriptions to the case of massive gravity. We show that electric charge is a critical parameter for categorizing the geodesic motions of both lightlike and timelike particles. In addition, we classify the type of geodesics based on the particle properties and geometry of spacetime.
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 [126], 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 [56]. Later on, many people studied various aspects of BTZ black hole as a typical laboratory of high energy physics [57][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 fourdimensional Schwarzschild-de Sitter [100], Kerr [101,127], Kerr-Newman [102], and Kerr-de Sitter [103] spacetime. Moreover, the higher dimensional solution of Schwarzschild, Schwarzschild-(anti-) de Sitter, Reissner-Nordström, Reissner Nordström-(anti-) de Sitter [104] and Myers-Perry spacetimes [105] 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 [106][107][108][109][110], charged dilatonic black holes [111], the singly spinning and (charged) doubly spinning black ring [112,113], (rotating) black string [114], Schwarzschild, and Kerr pierced by black string spacetimes [115,116].
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 Sec. II, and define acceptable regions of motion. After that, we investigate possible regions of motion and orbits for charged and uncharged black holes in Sec. III. Then, we completely classify geodesic motions and solve the geodesic equation around the (charged) BTZ black holes in Sec. IV. The neutral black holes in massive gravity are treated in Sec. V and the charged ones in Sec. VI. Finally, we end the paper with some concluding remarks.

II. THREE DIMENSIONAL LINE ELEMENT
In (2 + 1)−dimensional spacetimes, the line element takes the following form where ψ(r) is an arbitrary function of coordinate r which should be obtained based on the gravitational field equation.
In this paper, we will consider four different cases for the metric function ψ(r). The first one describes static BTZ black holes [56] where Λ denotes the cosmological constant and m 0 > 0 is an integration constant proportional to the total mass of the black hole. Note that here Λ < 0 is necessary to have a black hole solution. For a discussion of the singularity at r = 0 see [57]. Geodesic motions of massive/massless particles around the rotating uncharged BTZ black holes have been investigated in Refs. [117,118]. In the mentioned papers, it is 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] 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, Λ 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 (Λ ∝ −l −2 ), however, in [125], 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 [106]. Since the charge term of the metric function in Ref. [106] 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] where W(x) denotes the principal branch of the Lambert W function and W(−1, x) the branch with W(−1, x) ≤ −1.
The only physically acceptable solution of Eqs. (4) and (5) is The metric function of BTZ black holes in massive gravity is [90] ψ 3 (r) = −Λr 2 − m 0 + m 2 cc 1 r, 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 ′ to combine all massive parameters, as follow Note that in the massive case Λ < 0 is not necessary for black hole solutions and we may as well consider Λ > 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 ′ 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.
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 The Lagrangian L of a test particle is given by where ǫ takes the values 1 and 0 for massive and massless particles, respectively, and λ is the affine parameter for massless particles and the proper time for massive particles. From this, we then find the following geodesic equations dr dϕ In addition, considering Eq. (12), an effective potential can be introduced as 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) ≥ 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 eff . Equivalently, from (14) the acceptable region of motion is E 2 ≥ V 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].

III. 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 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 section II, the function P (r) has a double root r = 0 for all parameters. Therefore, we can introduce a new function P (r) by Instead of searching for roots of P (r) we may solve P (r) = 0 and d 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 → ∞. In the limit of large r we find For negative Λ, massive particles may therefore never reach infinity.

A. 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, 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 ∞, then approaches a periapsis r = r min and goes back to ∞.
• Terminating bound orbits: r starts in (0, r m ] for 0 < r m < ∞ and it falls into the singularity at r = 0.
• Terminating escape orbits: r comes from ∞ and falls into the singularity at r = 0.
The polynomial 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) → +∞ for r → ∞ always corresponds to an even number of roots whereas P (r) → −∞ corresponds to an odd number. Therefore, only the following five different orbit configurations are possible: • In Region (0), P (r) has no positive real roots and P (r) > 0 for r ≥ 0. The possible orbit types are terminating escape orbits.
• In Region (1), P (r) has one positive real root r 1 with P (r) ≥ 0 for 0 ≤ r ≤ r 1 and possible orbit types are terminating bound orbits.
• In Region (2), P (r) has two positive real roots r 1 ≤ r 2 with P (r) ≥ 0 for 0 ≤ r ≤ r 1 and r 2 ≤ r and possible orbit types are flyby and terminating bound orbits.
• In Region (3), P (r) has three positive real zeros r 1 ≤ r 2 ≤ r 3 with P (r) ≥ 0 for 0 ≤ r ≤ r 1 and r 2 ≤ r ≤ r 3 and possible orbit types are bound and terminating bound orbits.
• In Region (4), P (r) has four positive real zeros r 1 ≤ r 2 ≤ r 3 ≤ r 4 with P (r) ≥ 0 for 0 ≤ r ≤ r 1 and r 2 ≤ r ≤ r 3 and r 4 ≤ r and possible orbit types are bound, terminating bound and flyby orbits .

B. 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 P (r) approaches −∞ for r → 0. Therefore, based on Bolzano's theorem, since P (r) is continuous in (0, +∞), if Λ > 0 there is always a real and positive root.
Since the logarithmic charge term changes the asymptotic behaviour of P (r) for r → 0 from a positive finite value to negative infinity, this suggests the presence of a small additional root r 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 called two-world escape orbit and many-world bound orbit, respectively (see [104,114,119,120] for more detail).
Numerical calculations show that the following kind of orbits can be identified for the charged cases • Flyby orbits: r starts from ∞, then approaches a periapsis r = r min for r min > r + , and goes back to ∞.
• Bound orbits: r changes between two boundary values r 1 ≤ r ≤ 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.
For the charged case we can find six regions with different number of roots and, therefore, different possible orbits, with the following properties: • In Region (0), P (r) has no real positive root with P (r) < 0 and the motion is impossible .
• In Region (1), P (r) has one real positive root 0 < r 1 < r − with P (r) ≥ 0 for r 1 ≤ r and possible orbit types are two-world escape orbits.
• In Region (2), P (r) has two real positive roots 0 < r 1 < r − < r + < r 2 with P (r) ≥ 0 for r 1 ≤ r ≤ r 2 and possible orbit types are many-world bound orbits.
• In Region (3), P (r) has three real positive roots 0 < r 1 < r − < r + < r 2 ≤ r 3 with P (r) ≥ 0 for r 1 ≤ r ≤ r 2 and r 3 ≤ r and possible orbit types are many-world 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 Λ < 0, m ′ > 0.

IV. 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 [117,118], 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.
A. Uncharged BTZ black holes

General classification of motion
Substituting Eq. (2) into Eq. (13), we obtain Let us first consider the case of massive test-particles with ǫ = 1. Upon inspection of the effective potential (18) and keeping in mind that Λ < 0, we see that it diverges to infinity for r → ∞ and to minus infinity for r → 0. Its derivative 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 / √ −Λ. The same result can of course be inferred from (19). By Descartes' rule of signs, P (r) posses at most one positive real zero r 0 , and P (r) → −∞ for r → ∞. Massive test particles are therefore bound to the region 0 ≤ r ≤ r 0 . Now let us turn to massless particles, ǫ = 0. Here, it is convenient to rescale the affine parameter λ such that the equation of motion simplifies to where b = L/E is the impact parameter. Then, the effective potentialV eff approaches −Λb 2 from below for r → ∞ and diverges to minus infinity for r → 0. As the derivative ofV eff , is positive for r > 0, this implies that a photon may reach infinity if −Λ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.

Analytic solution of geodesic equations
We use the substitution r 2 = 1 u and slightly rewrite Eq. (19) to obtain where c 1 = E 2 L 2 m0 + Λ m0 + ǫ L 2 , c 2 = ǫΛ m0L 2 . This integral has the solution ln Equivalently, we find for u as function of ϕ, where ξ is related to u 0 as Therefore, r(ϕ) is 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 Sec. III.

B. Charged BTZ black holes
For the line element (1), the linearly charged solution of BTZ black holes has been obtained as [95] so from Eq. (13), P (r) and V eff are In contrast to the uncharged case, here we can find circular orbits, which are given by dr dλ = 0 and d 2 r dλ 2 = 0. These two conditions are equivalent to P (r) = 0 and dP dr = 0. We can solve these two equations for the squared energy and angular momentum as for massive particles (ǫ = 1). Let us discuss these two equations. From Eq. (31), it is evident that E 2 > 0 is valid only for a negative denominator and thus In addition, the numerator of Eq. (32) has to be negative, and therefore, keeping in mind that Λ < 0 Considering the constraints (33) and (34), simultaneously, we find that circular orbits exist if the following relation is satisfied, It is notable that such inequality implies that there exist no horizons 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 = 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 Sec. III.  Let us now turn to massless particles. In this case (ǫ = 0) circular orbits exist for where b = L/E is the impact parameter. Due to the fact that Γ is real, in order to have physically acceptable motion for the massless particles, the denominator of Eq. (36) must be negative. Therefore, the following constraint must hold, keeping again in mind that Λ < 0, This results in the same inequality (35) as for the case of massive particles, that is incompatible with the existence of an event horizon.
The many-world bound orbit of null geodesics are shown in Fig. 3.

A. 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 Let us first discuss massive particles (ǫ = 1). From the form of P (r), according to the Descartes rule of signs, it is clear that for Λ < 0, there may be one or three positive real zeros, whereas for Λ > 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 ′ . This points to the existence of circular orbits, which are given by dr dλ = 0 and d 2 r dλ 2 = 0. These two conditions are equivalent to P (r) = 0 and dP dr = 0. We may solve these two equations for E 2 and L 2 , which gives for massive particles (ǫ = 1) It is clear from the expression for E 2 that circular orbits can only exist if 2m 0 − m ′ r < 0, which is for r > 0 only possible if m ′ > 0, which then gives us r > 2m 0 /m ′ . The equation (42) for L then implies that m ′ − 2Λr > 0 has to hold, which is automatically fulfilled for Λ < 0. For Λ > 0 we find r < m ′ /(2Λ).
If the radius of the circular orbit corresponds to a maximum of P (r), i.e. if the second derivative of P (r) is negative, it is stable. The second derivative d 2 P dr 2 together with Eqs. (41) and (42) reads This expression can be solved for the radius of the circular orbit, For Λ < 0, we find one positive zero (for the negative sign before the root in (44)) 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 ′ for small Λ and 8m 0 /3m ′ for large negative Λ. We plotted the ISCO for fixed m 0 in Fig. 4.
Interestingly, for Λ > 0 stable orbits only exist if r c from (44) is real, that is, if and only if (m ′ ) 2 > 16Λm 0 . Note that (m ′ ) 2 < 4Λm 0 can be rewritten as m ′ /(2Λ) < 2m 0 /m ′ , which implies that no circular orbits exist according to our discussion of (41) and (42) above. Summarized, for m ′ > 0 and Λ > 0 circular orbits exist for r ∈ [2m 0 /m ′ , m ′ /(2Λ)], and may be stable if m ′ /(2Λ) > 8m 0 /m ′ . Then the circular orbits are stable in between the two radii given by (44) and we have an innermost and an outermost stable circular orbit. We show the important radii for the case m ′ > 0, Λ > 0 in Fig. 5. Now let us turn to massless particles (ǫ = 0), where such as before, 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 = −Λ. If 1/b 2 < −Λ the polynomial P (r) diverges to −∞ for r → ∞ 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 ≤ r ≤ r 0 . On the other hand, if 1/b 2 > −Λ, then P (r) is positive for r → ∞, 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 dP dr = 0 for the impact parameter b and the radius of the circular orbit r c , we find From Eq. (45) for b we infer that 4Λm 0 − m ′2 < 0, or equivalently Λ < m ′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 we find with the values above, which implies that the circular photon orbit is unstable. Interestingly, the radius of the photon orbit only depends on the ratio m 0 /m ′ , but not on Λ or the individual parameters m 0 and m ′ . For Λ 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 ′ . However, we correctly recover that the circular orbit vanishes (its radius shifts to infinity) for m ′ → 0. The results of Eqs. (41), (42) and (45)

B. Analytic solution of geodesic equations
Here, we discuss the analytical solution of Eq. (13). The function P (r) in given in Eq. (40) is a polynomial of degree 6, which can be written in the following form where the coefficients a i are 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 and the coefficients b j can be calculated as where α j = b j | ǫ=0 . The simple roots of Eq. (52) lead to an elliptic type differential equation. Another substitution u = 1 α3 (4y − α2 3 ) transforms P 3 (u) into the following form in which the coefficients g 2 and g 3 are It is known that the solution of Eq. (53) is the Weierstrass function ℘(ϕ) in the following form where in which r 0 and ϕ 0 are the initial values of the differential equation. Therefore, the general solution of Eq. (52) is 2. Timelike geodesics (ǫ = 1) By substituting (ǫ = 1) into Eq. (50), it can be transformed to where β j = b j+2 | ǫ=1 . Equation (58) is of hyperelliptic type and its analytical solution is given by the Kleinian sigma function [101,105,121] in which σ i (z) denotes the derivative of the Kleinian sigma function with respect to the i th component of z where 2ω, 2ω ′ is the period matrix, 2η, 2η ′ is the period matrix of the second kind, C is a constant that can be given explicitly but does not enter (59) and τ = ω −1 ω ′ . The theta function is defined as follow in which g, h 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 ϕ σ in (59) is defined as where f is given by the condition σ (ϕ σ ) = 0. The constant ϕ in reads and only depends on the initial values ϕ 0 and u 0 . The general solution of the radial coordinate r is given by Equations (57) and (64) 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-11. Types of orbits mentioned here have been introduced in the total classification section III.

A. General classification of motion
In the charged black holes in massive gravity, the metric function ψ(r) is obtained as [90] ψ(r) = −Λr 2 − m 0 + m ′ r − 2q 2 ln r r 0 .
From Eq. (13), we then find P (r) as Again, we can look for circular orbits by solving P (r) = 0 and dP (r) dr = 0 for the squared energy E 2 and angular momentum L 2 . For massive particles (ǫ = 1) we find As expected, for q = 0 Eqs. (67) and (68) reduce to Eqs. (41) and (42), respectively, and for m ′ = 0 they reduce to Eqs. (31) and (32). From equation (67) we infer the inequality For m ′ > 0 the function R diverges to infinity at r = 0 and r = ∞, and has only a single extrema, a minimum at r = 4q 2 /m ′ . It may have two zeros Γ m,0 < Γ m,−1 given by If m ′ ≤ 0 the denominator R is monotonically decreasing on (0, ∞) and introduces therefore an upper bound on r given by r < Γ 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. ψ 4 (r H ) = 0. Then we find from this Therefore, R is positive at the inner horizon and negative at the event horizon. For Λ > 0 also a cosmological horizon exists, where R is positive again. We conclude that r H,inner < Γ m,0 < r H,event < Γ m,−1 < r H,cosmo , if the respective horizons and zeros of R exist.
We can directly infer that for Λ < 0 and m ′ ≤ 0 circular orbits outside of a black hole cannot exist, as for charged BTZ black holes discussed in section IV B and uncharged BTZ black holes in massive gravity discussed in section V.
In addition, the nominator in (68) needs to be negative. For Λ < 0 this implies a lower bound on r given by Note that this bound is the minimum of ψ 4 and, therefore, is smaller than the event horizon. For Λ < 0 and m ′ > 0 circular orbits around a black hole exist at r > Γ m,−1 , as in the case of uncharged BTZ black holes in massive gravity discussed in section V.
From the condition that the nominator in (68) has to be negative, in the case Λ > 0 circular orbits are not allowed for m ′ ≤ 0, and have only a limited range for m ′ > 0 given by Note that the two bounds correspond to the minimum and the maximum of the metric function ψ 4 . It therefore only remains to show that Γ m,−1 is smaller than the upper bound in (74): as r * := Γ m,−1 is a zero of R we find ln It is clear that ψ 4 is positive at r * = Γ m,−1 which implies that r * is smaller than the upper bound in (74). We therefore find circular orbits for Λ > 0 and m ′ > 0 in the range Let us now turn to massless particles (ǫ = 0). We find for the circular orbits where r c is a solution of R = 0, see Eq. (69). Therefore, again for m ′ < 0 there is no circular photon orbit outside of a black hole. If m ′ > 0 a circular photon orbit exists at r = Γ m,−1 . In the limit m ′ = 0 this reduces to the results of section IV B, and for q = 0, the zero of R is r = 2m 0 /m ′ as in section V. The circular photon orbit is always unstable: if we insert b 2 and r into the second derivative of P , we find stability for which is however incompatible with a real Γ m,−1 .
The results of Eqs. (67), (68) and (78) for both massive and massless particles are given in Figs. 12-13. (For more details we refer the reader to Sec. III). 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 Λ < 0 and Λ > 0. In the case Λ < 0, it is impossible to reach radial infinity and bound orbits outside the horizons do not exist. On the other hand, for Λ > 0 all parameter combinations allow for orbits reaching infinity and bound orbits outside the horizons are possible. 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-15 (see also Sec. III for more details). We restricted to positive Λ 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 manyworld 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 figure 14(b) and (e), where a massive test particle approaches the black hole rather straight and is reflected back without revolving around the black hole.

VII. SUMMARY AND CONCLUSIONS
In this paper, we analyzed the geodesic motion in the three-dimensional (charged) BTZ black hole spacetime 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 black hole spacetimes with charge, where an analytical solution could not be found.
While the uncharged BTZ black hole only has a rather poor variety of orbital motion, 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 ′ introduced in Eq. (8) further enriches the structure of geodesic motion. We found that m ′ > 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 Λ > 0) the outermost stable circular orbit. The results are plotted in figures 4 and 5. In particular, this implies that stable bound orbital motion outside the event horizon is possible, and we showed an example in figure 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.