Fundamental photon orbits in the double Schwarzschild space-time

We consider photon orbits in the space-time consisting of two identical Schwarzschild black holes in static equilibrium. In particular, we focus on fundamental photon orbits – bound photon orbits that do not fall into either black hole or go to infinity. It is known that there could exist up to two circular photon orbits lying in the plane of symmetry between the black holes. In this paper, we present more examples of fundamental photon orbits in this space-time. This includes non-circular orbits lying in the plane of symmetry, circular orbits that lie off the plane of symmetry, polar orbits, and non-planar orbits that orbit around one or both of the black holes. These orbits are all unstable. Finally, in an appendix, we present a class of circular time-like orbits that lie off the plane of symmetry.


Introduction
It is well known that photons can orbit around a Schwarzschild black hole in a circle with radius equal to 3M, where M is the mass of the black hole in geometric units. Although this orbit is unstable, it plays an important role in the appearance of the black hole to a distant observer. In general, the black hole will appear as a shadow -a region of the observer's sky devoid of light. The boundary of this shadow is formed by photons which approach the unstable circular photon orbit when traced back asymptotically to the past.
For the rotating Kerr black hole, two circular photon orbits are known to exist in the equatorial plane -one prograde and one retrograde. More interestingly, other bound photon orbits are allowed [1]. These orbits have a constant coordinate radius, but are not necessarily confined to a plane; for this reason, they are known as spherical photon orbits. These orbits play an important role in the imaging of the Kerr black hole, in particular in the reconstruction of its shadow. a e-mail: eteo@nus.edu.sg (corresponding author) In general, a bound photon orbit that does not fall into the black hole or go to infinity is known as a fundamental photon orbit (FPO) [2]. A special but important case of an FPO is a circular photon orbit, also known as a light ring (LR). For the Kerr black hole, it can be shown that the only FPOs are the spherical photon orbits, with the two equatorial LRs as special cases.
For other examples of FPOs, one has to turn to other black hole space-times. It is natural to consider space-times with more than one black hole. In [3], a space-time consisting of two identical extremally charged black holes -known as the double Majumdar-Papapetrou space-time -was considered, and several interesting classes of FPOs were discovered. Amongst them are closed orbits that lie in a plane containing the symmetry axis of the space-time. These orbits can orbit around either one or both of the black holes. One particularly interesting example is an orbit that orbits around both black holes in the shape of a figure 8. The non-planar counterparts of these orbits were also discussed.
The extremal charge of the black holes considered in [3] is necessary to maintain a static equilibrium between them. It might be physically more relevant to consider uncharged black holes, such as Schwarzschild or Kerr black holes. In this paper, we shall consider the double Schwarzschild spacetime. This space-time contains two identical Schwarzschild black holes, held in static equilibrium by a conical singularity stretching between them along the axis of symmetry. This exact solution has been used as a model for the more realistic case of a dynamical binary black hole system in [4,5].
An early study of null geodesics in the double Schwarzschild space-time was made by Coelho and Herdeiro [6]. They showed that there exist two LRs lying in the plane of symmetry between the black holes, when the separation of the black holes is below a certain critical value. Most remarkably, this critical value was found to involve the reciprocal of the golden ratio ϕ ≡ ( √ 5 + 1)/2, also known as the golden ratio conjugate. The two LRs merge into one at this critical value, and cease to exist altogether above this critical value.
In this paper, we will find a number of new classes of FPOs in this space-time, and study their properties. We begin in Sect. 2 with a brief review of the double Schwarzschild spacetime and the relevant geodesic equations. In Sect. 3, we show that when the separation of the black holes is below the abovementioned critical value, there exists a class of non-circular FPOs lying in the plane of symmetry between the black holes. In Sect. 4, we show that when the separation is above the critical value, there exist two LRs lying symmetrically off the plane of symmetry. Thus the LRs found in [6] do not actually cease to exist, but rather they just move off the plane of symmetry.
In Sect. 5, we find a class of closed FPOs lying in a plane containing the symmetry axis of the space-time. We then turn to the existence of non-planar FPOs in Sect. 6, and show that there exists a class of such FPOs which orbit around both black holes, and another class which only orbit around one of the black holes. In Sect. 7, we show how the effective potential method provides an alternative way to visualise the FPOs. The paper concludes in Sect. 8 with some open questions and future directions. Although it is not the main focus of this paper to consider time-like geodesics, we present in the appendix a class of circular time-like orbits that lie off the plane of symmetry between the black holes.

Geodesic equations
It is well known that the general static axisymmetric solution of the vacuum Einstein equations has a metric which can be written in Weyl-Papapetrou coordinates as [7] ds 2 = −e 2U dt 2 + e 2k−2U dρ 2 + dz 2 + ρ 2 e −2U dφ 2 , (2.1) where U = U (ρ, z) and k = k(ρ, z) are functions of ρ and z only. In particular, U satisfies the Laplace equation for three-dimensional flat space in cylindrical polar coordinates, with the metric Once U is specified, k can be derived by solving the set of equations (2.4b) The standard Schwarzschild solution is obtained if U is taken to be the Newtonian potential of an infinitesimally thin rod of mass M and length 2M, placed somewhere along the z-axis [7]. The rod itself corresponds to the horizon of the black hole. To obtain the double Schwarzschild solution, we can take U to be the potential of two such rods, placed along the z-axis in the intervals a 1 ≤ z ≤ a 2 and a 3 ≤ z ≤ a 4 , for some a 1 , . . . , a 4 satisfying a 1 < a 2 < a 3 < a 4 [8]. The two black holes then have masses M 1 = (a 2 − a 1 )/2 and M 2 = (a 4 − a 3 )/2, and are separated by a coordinate distance 2L = a 3 − a 2 . For this solution, the functions U and k are explicitly given by [9] where we have defined Assuming that the azimuthal coordinate φ has the standard periodicity, there is a conical singularity along the z-axis between the two rods, for a 2 < z < a 3 . In this case, it corresponds to an excess angle given by .
The existence of this conical singularity is necessary to counterbalance the gravitational attraction of the two black holes, thus maintaining the system in static equilibrium.
In the rest of this paper, we shall assume that the two black holes have the same mass, i.e., M 1 = M 2 ≡ M. Without loss of generality, we can then choose a 1 = −L − 2M, a 2 = −L, a 3 = L and a 4 = L + 2M.
The geodesic equations can be derived in the standard way using either the Lagrangian or Hamiltonian formalism; we shall use the former here. For the metric (2.1), the geodesic equations are obtained by extremising the Lagrangian (2.8) where the dot denotes derivative with respect to an affine parameter along the geodesic. For the t and φ geodesic equations, we obtain the first-integrals: where E and are constants of motion, interpreted as the energy and angular momentum of the particle respectively. The ρ and z geodesic equations are second-order differential equations, given bÿ The partial derivatives of k can be calculated from those of U using (2.4). For the double Schwarzschild solution, the latter are explicitly For null geodesics, we also have the condition that the Lagrangian (2.8) vanishes. In this case, we can effectively set E = 1 through a rescaling of the affine parameter and ; this is what we will assume from now on. Substituting (2.9) into (2.8), we then have where H (ρ, z) is given by Since the left-hand side of (2.12) is non-negative, we must have | | ≤ H (ρ, z) for physical motion. This result forms the basis of the effective potential method, with H (ρ, z) being the effective potential [3,6]. 1 This method will prove to be very useful for studying certain qualitative properties of FPOs in the space-time. However, we will also need explicit examples of these orbits to understand their properties more fully. They are obtained by numerically integrating the four geodesic equations (2.9) and (2.10).
In this work, the geodesic equations were integrated numerically using the standard fourth-order Runge-Kutta method. The desired orbit was found by treating it as an initial value problem, and using the shooting method. The condition L = 0 was imposed at the initial point, and then used as a consistency check of the subsequent numerical integration. 1 Here, we follow the approach of [3]. The effective potential used in [6] is related to the one used here by V (ρ, z) = H (ρ, z) −2 .
In our simulations, it was found to be always satisfied to an accuracy of at least 10 −10 .
In searching for FPOs, we do not consider photon orbits that fall into either black hole. Thus, they should not approach the two rods along the z-axis. However, they are allowed to approach the remaining parts of the z-axis, which together make up the axis of symmetry of the space-time. In particular, they are allowed to approach arbitrarily close to the conical singularity that lies on the part of the z-axis between the two rods, even though we assume the photons cannot pass through the conical singularity itself.
We remark that the numerical integration might encounter problems near or on the axis of symmetry. This is because R i − ζ i vanishes on the part of the axis where z ≥ a i , leading to potential issues in the computation of e 2U and ∂ ρ U in (2.5a) and (2.11a), respectively. To address this problem, we can replace the exact expressions for e 2U and ∂ ρ U by their corresponding Taylor series expansions in ρ when a 2 < z < a 3 and z > a 4 . This will ensure the stability and accuracy of the numerical integration if the orbit approaches these parts of the symmetry axis.
Once the orbits have been obtained, we will plot them out in the three-dimensional flat space (2.3). Examples of such plots can be found in Figs. 2, 5, 8, 9, 11 and 12 below. In these plots, the x and y axes refer to the usual Cartesian coordinates defined by x = ρ cos φ and y = ρ sin φ. The rods representing the black hole horizons are indicated in red in these plots (except for Fig. 2). However, we have not indicated the conical singularity between them, to reduce clutter in the figures.

Equatorial FPOs
Since the plane of symmetry between the two black holes is a totally geodesic sub-manifold of the space-time [6], it is natural to begin by looking for photon orbits in this plane. It was found in [6] that there exist two LRs lying in this plane when L/M is less than the golden ratio conjugate 2 1/ϕ = ( √ 5 − 1)/2. The two LRs coincide when L/M is equal to this value. In this section, we investigate if there exists any other FPOs in the plane of symmetry. For brevity, we will henceforth refer to this plane as the equatorial plane [3], and orbits lying in it as equatorial FPOs.
The existence of the two LRs was derived using the effective potential method in [6]. In this case, we have to set z =ż = 0 in (2.12) and (2.13). When L < M/ϕ, the general form of H (ρ, 0) is shown in Fig. 1. As can be seen, it has two stationary points. Each stationary point corresponds Note that the general form of H (ρ, 0) in Fig. 1 also allows for the existence of a family of non-circular equatorial FPOs, whose angular momentum ranges between that of the inner and outer LRs. An example is indicated by the horizontal blue line, which describes an orbit whose radius ρ takes the finite range ρ min ≤ ρ ≤ ρ max . The angular momentum of this orbit is then given by the value of ±H (ρ, 0) at either the minimum radius ρ min or the maximum radius ρ max . It is clear from Fig. 1 that ρ max always lies between the two LRs, while ρ min always lies below the inner LR. The smallest value that ρ min can take occurs when ρ max is equal to the radius of the outer LR.
For a fixed value of L in the range 0 < L < M/ϕ, this family of non-circular equatorial FPOs can be conveniently parameterised by the value of ρ max . Two examples of such FPOs are plotted in Fig. 2, for the case L = M/2. As can be seen, these FPOs precess by a certain amount after every complete oscillation in radius. We note that analogues of these orbits have also been found in the double Majumdar-Papapetrou space-time [10], as well as in other black hole space-times [11,12].

Off-equatorial LRs
Since we know that LRs exist in the equatorial plane of the double Schwarzschild space-time, the next question is if LRs can exist off this plane of symmetry. In this section, we show the existence of LRs that lie in a plane of constant z = 0.
To do so, we use the effective potential method again. In this case, we just have to setż = 0 in (2.12). Stationary points of the potential H (ρ, z) in (2.13) then correspond to LRs, with angular momentum equal to the value of ±H (ρ, z) at the point. The conditions for a stationary point, ∂ ρ H (ρ, z) = 0 and ∂ z H (ρ, z) = 0, are equivalent to respectively. The set of all points (ρ, z) satisfying each of these conditions can be found numerically, using the expressions in (2.11). It turns out that the solution of (4.1a) generally consists of two curves in the (ρ, z)-plane. Similarly, the solution of (4.1b) generally consists of two curves in the (ρ, z)-plane, with one of them being the z = 0 line. The qualitative behaviour of these curves will depend on whether L is below, equal to, or above the critical value M/ϕ. We begin by considering the case L < M/ϕ. As an example, the set of points satisfying the conditions in (4.1) is plotted in Fig. 3a for L = M/2. The two orange curves are points at which ∂ ρ H (ρ, z) = 0, while the two blue curves are points at which ∂ z H (ρ, z) = 0. It can be seen that there are only two points -indicated by the black dots -at which one of the orange curves meets one of the blue ones, and they both lie along the z = 0 line. These are precisely the locations of the two equatorial LRs found in [6]. It is clear that there are no off-equatorial LRs in this case.
We next consider the critical case L = M/ϕ in Fig. 3b. In this case, the two orange curves and the two blue curves all meet at a single point. This, of course, describes the case in which the two equatorial LRs have merged into one. Again, it is clear that there are no off-equatorial LRs in this case.
Finally, we consider the case L > M/ϕ. As an example, the set of points satisfying (4.1) is plotted in Fig. 3c for L = M. The qualitative behaviour of the orange curves is now different. While there are still two points at which one of the orange curves meets one of the blue ones, they now both lie symmetrically off the z = 0 line, with the same value of ρ. This shows that there are two off-equatorial LRs with the same radius in this case.
Although we have only illustrated this for one specific value of L > M/ϕ, we have checked that it holds for other values of L in this range. Figure 4 shows how the radius ρ of the LRs and their separation z depends on L. Note that ρ monotonically decreases as L is increased from the critical value, and that it approaches the value √ 3M in the limit L → ∞. By contrast, z monotonically increases as L is increased, and in fact asymptotically approaches the line 2(M + L). These limiting values arise when the two black holes can be regarded as separate Schwarzschild black holes with no interactions between them. The radius √ 3M is that of the circular photon orbit around a Schwarzschild black (b) (a) Fig. 2 Examples of non-circular equatorial FPOs for the case L = M/2, with parameters as indicated. The choice of parameter for a is explained at the end of Sect. 6.2; the parameter for b was chosen to make equal to 10M exactly. The orbits actually start at ρ min in both plots, and terminate after the interval φ = 6π and 12π , respectively hole in Weyl-Papapetrou coordinates. 3 On the other hand, 2(M + L) is the separation between the centres of mass of the two black holes. This is also the separation between the photon orbits around each of the black holes, if the interactions between them are neglected.
We remark that the stability of these orbits can be checked in the usual way by calculating the Hessian . It turns out that all the LRs described in this section -including the equatorial oneshave a negative Hessian determinant, showing that these orbits are unstable. This is consistent with the general result found in [13]. To summarise, we have shown that the two LRs found in [6] for the case L < M/ϕ, continue to exist in the case L > M/ϕ, except that they no longer lie in the equatorial plane. Instead, they lie symmetrically off the equatorial plane, with the same radius. The transition between these two cases occurs at the critical value L = M/ϕ involving the golden ratio, when the two LRs merge into one.

Polar FPOs
The next class of orbits that we shall consider are closed FPOs that lie in a plane containing the z-axis. Since there is a 3 The relation between Weyl-Papapetrou coordinates and the usual spherical polar coordinates is given by ρ = √ r (r − 2M) sin θ and z = (r − M) cos θ [14]. Thus ρ = √ 3M corresponds to r = 3M, when θ = π/2. conical singularity along the z-axis between the black holes, we shall assume that the orbits cannot pass through the inner part of the axis. They can only pass through the outer parts of the axis, where z < a 1 or z > a 4 , as they orbit the black holes. Because of this, we shall refer to such FPOs as polar FPOs. These orbits necessarily have vanishing angular momentum .
Two examples of such FPOs are illustrated in Fig. 5, for the cases L = 2M and 5M. For sufficiently small L, these FPOs resemble ellipses (although we have checked that they are not true ellipses). These orbits become more elongated as L is increased. Two useful characteristics of these orbits are its radius ρ when it passes through the equatorial plane, and the maximum z coordinate reached by the orbit. We denote them by ρ z=0 and z max , respectively. In Fig. 6, we have plotted ρ z=0 and z max − L against L for this class of polar FPOs.
We note that ρ z=0 monotonically increases as L is increased, and that it approaches a constant value in the limit L → ∞. On the other hand, z max monotonically increases as L is increased, and in fact increases linearly with L in the asymptotic limit. The latter can be seen from the fact that z max − L tends to a constant value as L → ∞.
It is, in fact, possible to derive these asymptotic values, using the fact that the two black holes can be regarded as separate Schwarzschild black holes in this limit. The corresponding orbit in this limit is one that approaches a Schwarzschild black hole from infinity, and is deflected by an angle of 180 • before receding back to infinity. The asymptotic value of ρ z=0 is just the impact parameter b of this orbit. The value of b can be calculated to be approximately 5.35696M; see, e.g., The yellow regions denote the points at which ∂ ρ H (ρ, z) < 0 (equivalently, 1 − 2ρ∂ ρ U (ρ, z) < 0); this information will be needed in the appendix when we consider time-like orbits [15,16]. On the other hand, the asymptotic value of z max is the closest point of this orbit to the black hole. This value can be calculated to be approximately 3.52062M, using Eq. (5) of [16]. It is also possible to derive the values of ρ z=0 and z max when L = 0. This limit corresponds to the two black holes merging into a single Schwarzschild black hole with mass 2M. In this case, the polar orbit is just the circular photon orbit around the Schwarzschild black hole with radius r = 6M in spherical polar coordinates. In Weyl-Papapetrou coordinates (c.f. Footnote 3), this corresponds to the values ρ z=0 = 2 √ 3M and z max = 4M.

Non-planar FPOs
We have so far focussed on FPOs that lie in a plane in the three-dimensional space (2.3). In this section, we turn our attention to the existence of non-planar FPOs. We have found a few such classes of FPOs, which we now discuss in turn.

Class A non-planar FPOs
The first class of non-planar FPOs, which we call Class A, contains the polar FPOs described in Sect. 5 as a special case. A common feature of this class of orbits is that they all pass through the equatorial plane z = 0, with the property thaṫ ρ| z=0 = 0. It turns out that a possible choice of parameter for these orbits is the radius ρ of the orbit when it passes through the equatorial plane. As usual, we denote this by ρ z=0 .
Two useful characteristics of these orbits are the inclination angle that the orbit makes with respect to the equatorial plane in the three-dimensional space (2.3), and the maximum z coordinate reached by the orbit. We denote them by θ inc and z max , respectively. In particular, θ inc is given by the general formula: Without loss of generality, we can assumeż| z=0 ≥ 0, so the allowed range of θ inc is from 0 • (equatorial orbits) to 90 • (polar orbits). In Fig. 7, we show how θ inc and z max -as well as -depend on ρ z=0 , for several values of L. It is clear from Fig. 7 that the qualitative behaviour of the orbits is different, depending on whether L is more or less than the critical value M/ϕ. We begin by considering the range L ≤ M/ϕ. In this range, the lower and upper bounds  of ρ z=0 are finite non-zero values, as listed in Table 1. As ρ z=0 is decreased from the upper bound, both θ inc and z max will monotonically decrease. Two example orbits are plotted in Fig. 8a and b for the specific value L = M/2. As ρ z=0 reaches the lower bound, both θ inc and z max will vanish. Thus, we recover a circular orbit lying in the equatorial plane. Indeed, it can be checked that it is just the outer LR discovered in [6]. Next, we consider the range L > M/ϕ. In this range, the lower bound of ρ z=0 is zero while the upper bound is a finite non-zero value, as listed in Table 1. As ρ z=0 is decreased from the upper bound, both θ inc and z max will initially decrease, but may increase again. Two example orbits are plotted in Fig. 8c and d for the specific value L = M. As ρ z=0 → 0, both θ inc and z max will tend to finite non-zero values. Thus, we do not recover an equatorial orbit in this limit. In Fig. 9a, we have illustrated an orbit with a very small value of ρ z=0 . As can be seen, it resembles a figure-8 or lemniscate type of orbit. However, unlike the figure-8 orbit that is known to exist in the double Majumdar-Papapetrou space-time [3], the orbit in Fig. 9a does not lie in a single plane. This is due to the presence of the conical singularity on the z-axis between the black holes. It is known that when a geodesic passes close to a conical singularity, it will undergo a deflection by a certain angle [4].
This angle can be approximately derived as follows. Consider a constant z plane cutting through the conical singularity (so that a 2 < z < a 3 ). The metric of this two-dimensional plane, in the neighbourhood of the conical singularity itself, is given by up to an irrelevant conformal factor. This metric describes a two-dimensional flat space in polar coordinates, with polar angle given by φ = e −k(0,z) φ. Now consider a geodesic which passes through ρ = 0; its φ coordinate will undergo the usual shift φ → φ + π . However, in terms of the coordinate φ, this will correspond to the shift φ → φ + π e k(0,z) . Thus, the angle by which the geodesic will be deflected is given by This formula remains approximately valid for geodesics passing close to the conical singularity; a fact we have verified numerically for our examples. Since we do not allow light to pass through the conical singularity, the exact angle (6.3) will not be reached in practice.

Class B non-planar FPOs
The Class A FPOs that we discussed in the previous subsection have the property that they cross the equatorial plane z = 0 repeatedly. As such, they can be regarded as orbiting around both black holes. We now turn to a second class of orbits, which we call Class B, that are confined to either the region z > 0 or z < 0. Thus, they can be regarded as orbiting around only one of the black holes. Without loss of generality, we can take it to be the upper black hole. Because these orbits do not pass through the equatorial plane, we need to use a different parameter from the Class A FPOs. A possible choice of parameter is the maximum z coordinate reached by the orbit, denoted by z max . Two useful characteristics of these orbits are its radius ρ at z = z max , and the difference between the maximum and minimum z coordinates reached by the orbit. We denote them by ρ z max and z, respectively. In Fig. 10, we show how these two quantities -as well as -depend on z max , for several values of L.
Again, it is clear that the qualitative behaviour of the orbits is different, depending on whether L is more or less than the critical value M/ϕ. We begin by considering the range L > M/ϕ. In this range, the lower and upper bounds of z max are finite non-zero values, as listed in Table 2. We note from Fig. 10b that z = 0 at the lower bound. This in fact corresponds to an off-equatorial LR, of the type described in Sect. 4. As z max is increased, z will monotonically increase and the orbit will start tilting with respect to the z-axis. Two example orbits are plotted in Fig. 11a  will become almost vertical. However, as the lowest point of the orbit nears the conical singularity between the two rods, a deflection of the type described at the end of Sect. 6.1 will occur. In Fig. 9b, we have illustrated an orbit with a value of z max very close to the upper bound. In analogy with the non-planar figure-8 orbit in Fig. 9a, we shall call such an orbit a "non-planar figure-0 orbit". Next, we consider the range L ≤ M/ϕ. In this range, the lower bound of z max is zero while the upper bound is a finite non-zero value, as listed in Table 2. The upper bound again corresponds to a non-planar figure-0 orbit, when the orbit is most vertical. As z max is decreased, the orbit will become less inclined. Two example orbits are plotted in Fig. 11c and d for the specific value L = M/2. When z max = 0, we have z = 0 while ρ z max remains non-vanishing. This corresponds to an equatorial FPO, of the type described in Sect. 3. For L = M/2, we actually obtain the orbit shown in Fig. 2a; indeed, it can be seen that Fig. 2a is like a "flattened" version of Fig. 11d.

Higher-order non-planar FPOs
It turns out that there exist non-planar FPOs that are hybrids of the Class A and Class B FPOs when L > M/ϕ. Like the Class A FPOs, they all pass through the equatorial plane. Two such FPOs are plotted in Fig. 12, for the same specific values of L = M and = 6.5M as the Class A orbit in Fig. 8d; indeed, they can be regarded as higher-order generalisations of the latter orbit. In Fig. 12a, the orbit starts at the equatorial plane z = 0 and initially heads upwards like the Class A orbit in Fig. 8d. However, at the maximum value of z, it then orbits once around the upper black hole like the Class B orbit in Fig. 11b. After returning to the maximum value, it continues the rest of its orbit like the Class A orbit. In Fig. 12b, the first half of the orbit is similar to the previous one. However, at the minimum value of z, it orbits once around the lower black hole, before returning to the equatorial plane.
It is natural to label each of these orbits by a pair of integers (m, n), where m is the number of times it orbits around the upper black hole like a Class B orbit, and n is the number of times it orbits around the lower black hole like a Class B orbit subsequently. The lowest-order case (m, n) = (0, 0) corresponds to the orbit in Fig. 8d. The orbits in Fig. 12a and b have (m, n) = (1, 0) and (1, 1), respectively. One may wonder if there exist generalisations of these orbits to even higher order. Indeed, we have found examples of orbits with (m, n) as high as (2,2). Not surprisingly, the numerical integration has to be carried out to a very high precision to obtain them.
These higher-order (m, n)-orbits can be parameterised by the initial position and direction of the orbit in the equatorial plane. The initial position is, as usual, specified by the radius ρ z=0 . On the other hand, the initial direction is specified by the inclination angle θ inc defined in (6.1), as well as an azimuthal-type angle θ azi defined by  Note that θ azi is non-zero only for (m, n)-orbits with m = n; by symmetry, it vanishes for orbits with m = n. The corresponding parameters for the two higher-order orbits in Fig. 12 are indicated in their respective captions.

Effective potential method for FPOs
The effective potential method is an important and useful tool when analysing geodesics. Indeed, we have so far used it to deduce the existence of the equatorial FPOs in Sect. 3, and the existence and stability properties of the off-equatorial LRs in Sect. 4. In this section, we will expand on this method following [3], and show how it can be applied to understand the general behaviour of the FPOs found in this paper. In particular, it provides an alternative way to visualise the motion of the geodesics in the (ρ, z) plane. This is especially useful for the non-planar FPOs found in Sect. 6. Recall that the effective potential for the null geodesics considered here is given by H (ρ, z) in (2.13), and that the angular momentum of the photon must satisfy | | ≤ H (ρ, z) for physical motion. Equality occurs if and only iḟ ρ =ż = 0, as can be seen from (2.12). Note that, without loss of generality, we can assume ≥ 0; the case < 0 can then be obtained by a reflection φ → −φ. Now, for fixed values of M, L and , one can consider the set of points (ρ, z) satisfying the equation known as level curves of H (ρ, z). They delineate the boundary between the allowed and forbidden regions of the (ρ, z) plane for a null geodesic. Sinceρ =ż = 0 along the level curves, a null geodesic will turn around when it reaches such a curve.
As an example, consider the case L = M/2 and = 10M in Fig. 13a. The shaded regions indicate the parts of the (ρ, z) plane for which < H (ρ, z), where null geodesics are allowed to exist. In this case, there are two allowed regions, one connected to the black holes and the other to infinity.
As the values of the parameters are varied, the level curves will change. The rest of the plots in Fig. 13 illustrate this. Figure 13b shows the situation when is decreased below a certain critical value -approximately equal to 9.78603M in this case. At the critical value itself, the right-most two curves in Fig. 13a will meet up at a single point; this point corresponds to the outer LR in the equatorial plane. Below the critical value, the two curves will intercommute, as shown in Fig. 13b for the case = 8.8M. As a result, the two allowed regions in Fig. 13a will join up to become a single region, connected to both black holes and to infinity.
On the other hand, Fig. 13c shows the situation when is increased above a certain critical value -approximately equal to 10.3154M in this case. At the critical value itself, the leftmost two curves in Fig. 13a will meet up at a single point; this point corresponds to the inner LR in the equatorial plane. Above the critical value, the two curves will intercommute, as shown in Fig. 13c for the case = 10.5M. As a result, there are now three allowed regions: two of which are connected to either black hole, and one to infinity. Figure 13d and e illustrate examples of level curves when the black holes are separated by an intermediate and far distance, respectively. There are still three allowed regions in these cases, although the shapes of the level curves are different. In particular, it is clear from Fig. 13e that the two black holes are beginning to behave like separate Schwarzschild black holes. In such cases, there is a critical value of for which the two left curves will meet up with the right curve. The two points at which they meet up correspond to the two off-equatorial LRs. If is decreased below the critical value, we have the situation shown in Fig. 13f. The three allowed regions in Fig. 13e will join up to become a single region, connected to both black holes and to infinity. We note that in the limit = 0, the three level curves in Fig. 13f will approach the z-axis itself. The allowed region therefore becomes the whole (ρ, z) plane, consistent with the fact that the null geodesics are now allowed to reach the zaxis. This is in fact a general result that applies for any value of L.
The six sets of level curves shown in Fig. 13 can be divided into three topologically distinct types, represented by Figs. 13a, b and c. We call them type (a), (b) and (c), respectively. Note that Fig. 13d and e belong to type (c), while Fig. 13f belongs to type (b). There are no other topologically distinct types. Indeed, it can be checked that the L-parameter space is divided into three parts as shown in Fig. 14, corresponding to types (a), (b) and (c).
As we have already seen, in all three cases, the allowed regions of the (ρ, z) plane are connected to the black holes and/or to infinity. This shows that it is kinematically possible for a null geodesic to reach the black holes or go to infinity. For all the FPOs found in this paper, we have checked that a small perturbation in the appropriate direction will cause them to fall into one of the black holes or move off to infinity. They are thus unstable. 4 Despite being inherently unstable, an FPO will in general trace out a trajectory in the allowed region of the (ρ, z) plane -that neither falls into a black hole nor moves off to infinity. Moreover, it can be shown that these trajectories will touch the level curves at right angles [3], provided one uses the same scale for the ρ and z axes. As examples, we have plotted out three actual FPOs in Fig. 13a and b: (i) the equatorial FPO shown in Fig. 2b; (ii) the Class A FPO shown in Fig. 8a; and (iii) the Class B FPO shown in Fig. 11c.
These examples highlight an immediate advantage of visualising the FPOs in this way: they trace out closed trajectories in the (ρ, z) plane, even though they are in general open trajectories in the three-dimensional (ρ, z, φ) space. This greatly simplifies the search for non-planar FPOs, as the shooting method can be carried out in the (ρ, z) plane instead of the (ρ, z, φ) space. It also provides a way to classify the different types of non-planar FPOs. Note that the two blue curves (ii) and (iii) touch different level curves in Fig. 13b, at least at one end. This suggests an alternative way to distinguish between Class A and Class B FPOs in this topological class by the behaviour of the curves: Class A orbits trace out curves that stretch symmetrically between the two long level curves, whereas Class B orbits trace out curves that stretch between the short level curve and one of the long level curves (without crossing the z = 0 line).
Appropriate generalisations to the higher-order FPOs can also be made. In Fig. 15, we have illustrated the situation corresponding to the two higher-order orbits in Fig. 12. Note that the (1,0)-orbit traces out a curve in the (ρ, z) plane that stretches between the short level curve and the lower long level curve. It also very closely approaches the upper long level curve, although they do not actually touch. On the other hand, the (1,1)-orbit traces out a curve in the (ρ, z) plane that touches the short level curve symmetrically at both ends, while very closely approaching both the upper and lower long level curves.

Conclusion
It has been our main purpose in this paper to investigate the existence of FPOs in the double Schwarzschild space-time.
To this end, we have found a number of new classes of FPOs, extending on the two equatorial LRs discovered in [6]. These FPOs are summarised in Table 3, together with two important properties: the corresponding range of the black hole separation L, and the topological type of the level curves (7.1). It can be seen that the critical separation M/ϕ found in [6] continues to be an important determinant of the existence and properties of the FPOs. The results we have found underscore the rich behaviour of null geodesics in space-times with more than one black hole. We now conclude with a few open questions and directions for future work.
Since the geodesic equations for this space-time cannot be solved analytically, much of our work has been numerical. Although we have tried to be exhaustive in our search for new (a) (b) Fig. 15 The higher-order FPOs in Fig. 12, plotted in the (ρ, z) plane. The insets show that the orbits (blue curves) do not actually touch the level curves (dark grey curves) at the highest point -and by symmetry, the lowest point in b as well Off-equatorial LRs Fig. 11a, b Higher-order FPOs L > M/ϕ (b) Fig. 12a, b classes of FPOs, it is possible that there remains undiscovered classes of FPOs. Also, many of the properties of the FPOs discovered in this paper were studied numerically. It might be worth exploring if it is possible to rederive some of these results analytically.
One assumption made in this work is that the two black holes have equal mass. It is possible to relax this assumption, for example, by varying the mass of one of the black holes while keeping the distance L between them fixed. For small variations, the FPOs discussed in this paper should continue to exist. Indeed, we have verified that the two LRs will continue to exist, although their positions will shift about as the mass of the black hole is increased. The behaviour of the other classes of FPOs remains to be investigated.
Another possible extension of this work is to investigate the nature of time-like orbits in the double Schwarzschild space-time. It is expected that these orbits will exhibit a rich behaviour, similar to their null counterparts. For a start, we have found a class of circular time-like orbits that lie off the plane of symmetry, extending the circular null orbits of Sect. 4. These orbits are described in the appendix. The inves-tigation of other classes of time-like orbits will be left for the future.
It should also be possible to extend this work to the case of more than two black holes, and even to charged black holes. We mention that the existence of equatorial LRs and time-like geodesics in these more general space-times were already considered in the pioneering work of Coelho and Herdeiro [6].
More recently, the authors of [17][18][19] have uncovered a certain topological charge (TC) of a given space-time, equal to the number of stable LRs minus the number of unstable LRs. It is invariant under smooth deformations of the spacetime obeying fixed boundary conditions. So far, most of the work on the TC has been for space-times containing no or one black hole horizon. Thus, it is timely to investigate the TC in space-times with more than one horizon present. We note that TC = −2 for the double Schwarzschild space-time.
It is our hope that the FPOs discovered in this paper would eventually lead to a better understanding of aspects of the double Schwarzschild space-time, such as its shadow. In particular, the existence of the higher-order non-planar FPOs suggests that the shadow might exhibit fractal behaviour [3], and it would be interesting to study this and other chaotic features of the space-time in more detail.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: All data was obtained using the formulae and procedures as described in the article.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 . SCOAP 3 supports the goals of the International Year of Basic Sciences for Sustainable Development.

Appendix: Circular time-like orbits
In this appendix, we shall discuss the existence of circular time-like orbits in the double Schwarzschild space-time. In the case L < M/ϕ, it was found in [6] that such orbits could exist anywhere in the equatorial plane outside the annular region bounded by the two equatorial LRs. When L = M/ϕ, this annular region becomes a circle, but the above result still holds. Here, we shall focus on the case L > M/ϕ. We shall consider time-like orbits in a constant z plane. Substituting (2.9) into (2.8), and using the fact thatż = 0 and that L can always be rescaled to − 1 2 for such orbits, we obtain where the effective potential V (ρ, z) is given by [6] V (ρ, z) = ρ −2 e 4U 2 + e 2U .
Stationary points of this potential then correspond to circular time-like orbits, with energy E equal to the value of √ V (ρ, z) at the point. It can be checked that the conditions for a stationary point, ∂ ρ V (ρ, z) = 0 and ∂ z V (ρ, z) = 0, are equivalent to respectively. Note that the condition (A.3b) is identical to that for the null case in (4.1b). This means that the set of points (ρ, z) satisfying this condition is given by the same two blue curves described in Sect. 4 (c.f. Fig. 3). The condition (A.3a) is however different from the null case; it only reduces to (4.1a) if the right-hand side vanishes, i.e., when → ∞. Since the right-hand side of (A.3a) is in general a positive quantity, points in the (ρ, z) plane satisfying this condition are bounded by the orange curves described in Sect. 4.
In Fig. 3, we have indicated the set of points for which 1 − 2ρ∂ ρ U < 0 by the yellow regions. The parts of the blue curves which lie outside the closure of the yellow regions then indicate the allowed positions of the circular time-like orbits. For the case L ≤ M/ϕ, it agrees with the results found in [6]. For the case L > M/ϕ, it shows that these orbits could lie in the equatorial plane, and also possibly on the curve that stretches between the two off-equatorial LRs.
To find the actual positions of these orbits, one has to solve for the condition (A.3a) for different values of . It turns out that its solution generally consists of three curves in the (ρ, z)-plane. For definiteness, consider the case L = M; other cases will be similar. In Fig. 16, we have illustrated the set of points (ρ, z) satisfying the two conditions (A.3), for four qualitatively different cases. The orange curves are points at which ∂ ρ V (ρ, z) = 0, while the blue curves are points at which ∂ z V (ρ, z) = 0.
When is below a certain critical value -approximately equal to 6.53321M in this case, there will exist a single circular time-like orbit in the equatorial plane. This is illustrated in Fig. 16a for the case = 6.3M; it can be seen that the short orange curve intersects the blue curves at one point, with z = 0. When is equal to this critical value, the orange curve in question will meet the blue curves at the point where the latter two curves intersect. When is increased above this critical value, the short orange curve will intersect the blue curves at three distinct points, one with z = 0 and the other two with z = 0. This is illustrated in Fig. 16b for the case = 6.56M. Thus, we have one circular orbit in the equatorial plane and two circular orbits that are lifted off the equatorial plane. The latter orbits are the time-like analogues of the off-equatorial LRs described in Sect. 4.
It turns out that the time-like case exhibits additional features that are not present in the null case. If we increase to the approximate value 6.56860M, the two long orange curves will meet along the z = 0 line, giving rise to an additional circular orbit in the equatorial plane. Above this value, the two curves will intercommute, as shown in Fig. 16c for the case = 6.8M. There are now a total of three circular orbits in the equatorial plane, in addition to the two off-equatorial ones.
If we continue increasing to the approximate value 7.07266M, the left-most two orange curves in Fig. 16c will meet along the z = 0 line, causing the inner-most two circular orbits in the equatorial plane to merge. Beyond this value, the two curves will intercommute, as shown in Fig. 16d for the  Fig. 3c; the orange curves will always lie outside these regions case = 7.3M. There now remains just one circular orbit in the equatorial plane, in addition to the two off-equatorial ones. The stability of these circular time-like orbits can be checked by calculating the Hessian matrix . The black dots in Fig. 16 indicate orbits that have a negative Hessian determinant, and so are unstable. In particular, this includes all the off-equatorial orbits. On the other hand, the grey dots in Fig. 16 indicate orbits that have a positive-definite Hessian matrix, and so are stable. We conclude that there can exist up to two stable circular time-like orbits in the equatorial plane of this space-time.