Tidal Forces in Majumdar-Papapetrou Spacetimes

Tidal disruption events occur when astrophysical objects are destroyed by black holes due to strong tidal force effects. Tidal forces have been studied in a variety of black hole spacetimes, including Reissner-Nordstr\"om and Kerr spacetimes. Despite the vast literature on the subject, tidal forces around black holes in static equilibrium have never been investigated before. The aim of this work is to fill in this gap and explore tidal forces in the Majumdar-Papapetrou spacetime describing two extremely charged binary black holes in equilibrium. We focus on tidal forces associated with radial and circular geodesics of massive neutral particles moving on the plane equidistant to the black holes. In particular, we study the behavior of the tidal forces as a function of the distance from the black holes and as a function of the energy of the geodesics. We also investigate the numerical solutions of the geodesic deviation equation for different initial conditions.


I. INTRODUCTION
Recent observational results [1][2][3] have increased the interest of the physics community and the general public in black holes (and also in the wide range of phenomena associated with them).The possibility of detecting and studying effects like quasinormal ringing, shadows, and superradiance through radio telescopes and gravitational wave detectors has become a reality in the last decade.The prospect of new groundbreaking results of gravitational physics related to black holes and other compact objects is very promising.
Since the first detection of gravitational waves, more than 90 events consisting of the inspiralling and merger of binary systems of compact objects have been observed by the LIGO-Virgo collaboration in its first three observing runs [4].Theoretical and numerical studies of such binaries are crucial for a greater understanding of current and future observations.Unfortunately, there are no known analytical solutions of Einstein's equations that describe astrophysical binary systems.Nevertheless, exact and simple solutions corresponding to the metrics of non-coalescing pairs of black holes exist.One of such solutions is the Majumdar-Papapetrou (MP) metric [5,6] of extremely charged black holes in static equilibrium.The MP spacetime has been successfully used in the past as a toy model to investigate how the presence of a second black hole influences quasinormal ringing, shadows, Penrose energy extraction, and other typical black hole effects [7][8][9][10][11][12][13][14].
Building on recent articles [15][16][17][18][19][20][21][22][23] that have analyzed tidal forces in black hole spacetimes (including Reissner-Nordström and Kerr), the objective of this work is to investigate tidal forces acting on neutral massive particles that move along geodesics in the spacetime of two static black holes described by the MP metric.Tidal forces arise, for instance, when one considers a body falling towards the Earth.Due to the gravitational interaction, the body suffers stretching in the direction of motion and compression in transverse directions.In particular, tidal forces in the Schwarzschild spacetime can be used to describe the orbital motion of satellites around Earth and estimate deviations in their trajectories [24].
This work is organized as follows.In Section II, after introducing the MP spacetime for two black holes in equilibrium, we obtain the equations of motion for massive particles moving along geodesics that lie on the plane that is equidistant from the black holes.In Section III, we determine a freely falling frame, which is adapted to the geodesics discussed in Section II.In Section IV, we review the geodesic deviation equation of General Relativity and write it down explicitly for the MP metric.Using this equation, we analyze the properties of the tidal forces acting on radial and circular geodesics.In Section V, we solve the geodesic deviation equation numerically and analyze the behavior of the corresponding solutions as a function of the radial distance.Finally, in Section VI, we review the main results of our work and discuss future research directions.

II. THE MAJUMDAR-PAPAPETROU SPACETIME
The MP spacetime describes a collection of maximally charged black holes that are in static equilibrium due to the balance of electromagnetic and gravitational interactions [5,25,26].The MP metric, in Weyl cylindrical coordinates (t, ρ, ϕ, z), is where U (ρ, z) denotes the associated electromagnetic potential.In this work, we focus on a binary system of equalmass black holes and assume that the black holes are located along the z-axis, at z = ±b (the quantity 2b, hence, measures the separation between the black holes).With these assumptions, the electromagnetic potential takes the explicit form: where M denotes the mass of each black hole (each black hole has electric charge Q = M ).In particular, when b = 0, the metric describes a single extremal Reissner-Nordström black hole of mass 2M .In our numerical analyses and figures, we have chosen M = 1, which corresponds to normalizing the variables and parameters with respect to the mass of the black holes.
We consider timelike geodesics that correspond to the trajectories of massive particles subjected exclusively to the gravitational interaction.Geodesics, being curves that maximize the proper time between two events in a spacetime, can be determined from the Euler-Lagrange equations for the following Lagrangian L: where Taking advantage of the fact that t and ϕ are cyclic coordinates of the Lagrangian, one can determine two invariants of motion [11,27]: the energy E and the angular momentum L, defined, respectively, as These constants of motion can also be obtained from the fact that ∂ t and ∂ ϕ are Killing vector fields of the spacetime associated, respectively, with the stationarity and the axisymmetry of the MP metric (2.1).We restrict our analysis to timelike geodesics confined to the plane z = 0. Note that these planar geodesics only exist when the masses of the black holes are the same (as we have assumed in this work).The restriction z = 0, together with (2.5)-(2.7),implies the following equation of motion: where V eff denotes the associated effective potential.This is a first-order ordinary differential equation for ρ, which can be solved after one chooses an initial condition.From the solution ρ(τ ) of (2.8), one can determine t(τ ) using Equation (2.6) and ϕ(τ ) using Equation (2.7) to find the complete trajectory (t(τ ), ρ(τ ), ϕ(τ ), 0).To investigate tidal forces in the z = 0 plane, we shall assume that the geodesics are either radial (in which case, φ = 0 and, consequently L = 0) or circular (in which case, ρ = 0 and L ̸ = 0).Note that, according to Equation (2.8), radial geodesic motion on the z plane, corresponding to ϕ(τ ) = ϕ 0 = constant, is possible only if the absolute value of the energy is greater than a minimum value E min given by In particular, if E min ≤ |E| < 1, the trajectory is bounded and has a turning point at (2.10) Additionally, according to Equation (2.8), circular geodesic motion on the z plane, corresponding to ρ(τ ) = ρ 0 = constant, is only possible when E, L, and ρ 0 satisfy the following requirements: and Given L, one can first use (2.11) to determine the value of ρ 0 corresponding to a circular geodesic.The associated energy E can then be directly determined from (2.12).From Equation (2.7), we find that the angle changes linearly with respect to the proper time for circular geodesics: where ϕ 0 is the angle at τ = 0.In Figure 1, we show the values of E and L as a function of the radius of the circular geodesic and the separation between the black holes.We have observed that the values of E and L grow without bound when b and ρ 0 approach the white dashed line in Figure 1.Inside the black region, Equations (2.11) and (2.12) do not yield real solutions, meaning that circular geodesics are not allowed for the corresponding parameters b and ρ 0 .

III. FREELY FALLING FRAMES IN THE MAJUMDAR-PAPAPETROU SPACETIME
To investigate tidal forces along a geodesic, it is convenient to define a freely falling frame, i.e., a tetrad basis, which is adapted to the corresponding trajectory (meaning that the elements of the tetrad are parallel propagated along the geodesic).A tetrad basis {e â} = {e 0, e 1, e 2, e 3} is a non-coordinate basis that follows a prescribed normalization [28][29][30].We denote each component (according to the coordinate basis ∂ µ ) of the tetrad element e â as e µ â.The components of the corresponding dual vectors are e µâ .We assume that the tetrad basis is normalized according to e µ âe µ b = η âb , where η âb = diag (−1, 1, 1, 1).We set the first element of the tetrad basis as the tangent vector of the timelike geodesic we are interested in, i.e., e µ 0 = v µ .The remaining elements of the tetrad basis are determined using the normalization prescribed above and the requirement that the tetrad elements are parallel transported along the trajectory, i.e., where D/Dτ = v µ ∇ µ = e µ 0∇ µ is the directional derivative along the geodesic.We derive below a general expression for a freely falling frame associated with geodesics confined to the plane z = 0 of the MP spacetime.Equations (2.6)-(2.8)determine the tangent vector of the geodesic and, therefore, also determine the first vector of the tetrad basis: By definition, λ µ 0 is parallel propagated along the geodesic.It is straightforward to find another vector, which is parallel propagated along the geodesic and which satisfies the requirements of the tetrad basis.In fact, the symmetries of the problem, together with the fact that the geodesic is confined to the z = 0 plane, lead to: The difficult part is to find two extra vectors that are parallel propagated along the geodesic.In order to do so, we follow the idea used in Ref. [31] to investigate tidal effects in the Kerr metric.We start by determining two vectors, λ µ 1 and λ µ 2, which, although not parallel propagated along the geodesic, are such that λ µ âλ µ b = η âb : To simplify the notation, we have omitted the dependence (ρ, 0) from U (ρ, 0) in the expression above.Note that one can rotate the vectors λ µ 1 and λ µ 2 by an arbitrary angle while maintaining the orthonormalization λ µ âλ µ b = η âb .Using this property, we can determine under which conditions the rotated vectors become parallel transported along the geodesic.More precisely, if Ψ = Ψ(τ ) denotes the rotation angle, the rotated vectors are given by By requiring that D(e µ 1)/Dτ = D(e µ 2)/Dτ = 0, we find that the angle Ψ must satisfy the equation Hence, if we choose Ψ according to (3.6) and define e µ 0 = λ µ 0, e µ 3 = λ µ 3, the tetrad basis {e â} will be parallel propagated along the geodesic, whose tangent vector is (3.2).
Note that, in general, the angle Ψ will change along the geodesic.Since E cannot be zero, Ψ will remain constant along the geodesic only when L = 0.In fact, if we let L = 0 and Ψ = 0, Equations (3.2)-(3.5)reduce to: which, therefore, constitutes a freely falling frame for radial geodesics on the z plane.
The angle Ψ also takes a simple form for circular geodesics.Considering the fact that ρ = ρ 0 is constant along such geodesics, one finds that dΨ/dτ is also constant.Taking into account Equations (2.11)-(2.12),one finds the following freely falling frame for circular geodesics: where The parameter Ψ 0 , which is independent of ϕ, represents the rotation angle between the tetrads {e â} and {λ â} when ϕ = 0.In particular, if we choose Ψ 0 = 0, then the tetrads will coincide when ϕ = 0.More generally, we can always set Ψ 0 so that the tetrads coincide at a particular angle ϕ (or, equivalently, at a particular proper time τ ).

IV. TIDAL FORCES IN THE MAJUMDAR-PAPAPETROU SPACETIME
Tidal forces in General Relativity are associated with the geodesic deviation equation: which determines the evolution of the geodesic deviation vector ξ µ that connects infinitesimally close geodesics in terms of the Riemann tensor R µ νρδ and the tangent vector v µ (τ ) of the geodesic.In simple terms, a non-vanishing Riemann tensor indicates the existence of acceleration between neighboring geodesics and, consequently, the presence of tidal forces.
To derive the geodesic deviation equation, one chooses a one-parameter family of geodesics around a reference geodesic x µ (τ ) [29,32].The family of geodesics defines a smooth two-dimensional surface x µ (s, τ ) on the spacetime in such a way that, for each s = s 0 , the curve x µ (s 0 , τ ) is a geodesic parametrized by the affine parameter τ with corresponding tangent vector v µ (s 0 , τ ) = (∂ x µ /∂τ )(s 0 , τ ).The parameter s, defined in a neighborhood of s = 0, parametrizes deviations from the reference geodesic x µ (0, τ ) = x µ (τ ).The deviation vector of a given geodesic in the family is defined by ξ µ (s 0 , τ ) = (∂ x µ /∂s)(s 0 , τ ).In particular, we denote the deviation vector ξ µ (0, τ ) of the reference geodesic as ξ µ (τ ).Since (s, τ ) can be used as a coordinate system on the two-dimensional surface, the Lie derivative of the vector field ξ µ along v µ must vanish, meaning that the vector fields ξ µ and ṽµ commute [29,32]: By taking the derivative of the equation above along the direction of v µ , one obtains the geodesic deviation equation after some algebraic manipulation.Note that there is a critical difference between Equations (4.1) and (4.2).While Equation (4.1) is a second-order differential equation that is well defined independently of the choice of the one-parameter family of geodesics around x µ (τ ), Equation (4.2) is a first-order differential equation that only makes sense if the family of geodesics is given a priori.This difference can be understood in the following sense: Equation (4.1) has a larger space of solutions in comparison to Equation (4.2) since the arbitrariness of choosing the family of geodesics around x µ (τ ) has not been used in Equation (4.1).In fact, the extra degrees of freedom in the initial conditions of Equation (4.1) (in comparison to the initial conditions of Equation (4.2)) correspond to this freedom of choosing the one-parameter family of geodesics.In the language of Riemannian and pseudo-Riemannian manifolds, the geodesic equation is typically referred to as the Jacobi equation.A key result in the mathematical literature is that every solution of the geodesic Equation (4.1) corresponds to a one-parameter family of geodesics defined around x µ (τ ) (see, e.g., Proposition 10.4 of Ref. [33]).
When analyzing geodesic deviations, it is convenient to apply the projection operator h µ ν = δ µ ν − v µ v ν on the deviation vector ξ µ to obtain the so-called orthogonal connection vector η µ : 3) The components of the orthogonal connection vector in the tetrad basis, η â, can be related to the components in the coordinate basis through An important characteristic of the orthogonal connection vector is the fact that the component along the direction of the associated geodesic vanishes, i.e., η 0 = 0.With the help of the projection operator, the geodesic deviation equation in terms of η â becomes [29] where K âb denotes the tidal tensor associated with the geodesic.Note that we can replace D/Dτ in the equation above by d/dτ = e µ 0∂ µ .To determine the explicit form of the tidal forces in the MP spacetime, we need the nonzero components of the Riemann tensor, which are (apart from the trivial symmetries R µνϵδ = R ϵδµν = −R νµϵδ ): (4.6)

A. Radial Geodesics
For radial geodesics in the z = 0 plane, we employ the freely falling frame given by (3.7) and find that the associated tidal tensor is diagonal.The non-trivial components of the geodesic deviation equation are given by η1 We remark that the potential U and its derivatives in the expression above must be evaluated at z = 0.The tidal forces in each direction correspond to the ratios between ηâ and η â.The tidal forces along the radial (ρ), angular (ϕ), and vertical (z) directions are shown, respectively, in Figures 2-4.For the sake of comparison, we have included in the left panels of each figure the corresponding tidal force when the MP metric describes an extremal Reissner-Nordström black hole (i.e., when b = 0).Figure 2 tells us that, for sufficiently large values of the radial coordinate ρ, the radial tidal force is positive, meaning that it is a stretching force.As the origin ρ = 0 is approached, the tidal force component attains a maximum value and then decreases, becoming negative (compression force) for sufficiently small values of ρ, attaining a minimum value at ρ = 0.According to (4.7), we see that the radial tidal force is independent of the energy E. The dependence on the separation b, however, is non-trivial.We note that increasing the separation b between the black holes typically decreases the magnitude of the tidal forces at z = 0.This is expected, since the black holes move away from the geodesics at z = 0 when b is increased.
Unlike the radial tidal force, the angular force is always compressive for radial geodesics on the z = 0 plane.This compressive force increases in magnitude with decreasing ρ.The magnitude of the compressive force also increases when the separation b between the black holes decreases and when the energy E associated with the geodesics increases.Finally, the vertical component of the tidal force is compressive for large values of the radial coordinate ρ and stretching for small values of ρ.The stretching force is most intense at the origin and, then, decreases monotonically with increasing ρ, changes sign, and reaches a local minimum, where the compressive force is most intense.As in the case of the angular component, the tidal force along the z direction increases with the energy E and decreases with the separation b.We note that the explicit form of Equation (4.7), after the substitution of the potential U given in Equation (2.2), does not provide any useful analytical insights for the angular and vertical components.The expression for the radial component, on the other hand, is sufficiently simple, allowing one to determine an analytical expression of the coordinate where the tidal force vanishes.Explicitly, one has meaning that the radial tidal force vanishes when

B. Circular Geodesics
For circular geodesics in the z = 0 plane, we employ the freely falling frame given by (3.8) and find that the corresponding tidal tensor mixes the radial and angular components.The non-trivial components of the geodesic deviation Equation (4.5) are The coefficients κ 1 , κ 2 , κ 3 , α, and β of the matrix K , which are constant along the circular geodesic, are given by We remark that the potential U and its derivatives in (4.11) must be evaluated at (ρ, z) = (ρ 0 , 0).According to (4.10), the component of the tidal force in the z direction is always decoupled from the other components, while the components in the other directions are coupled and oscillate as Ψ = Ψ(ϕ) changes along the geodesic.The tidal forces decouple when Ψ = Ψ(ϕ) is an integer multiple of π.
As mentioned previously, we have the freedom to align the tetrads {e â} and {λ â} at a specific point of the circular geodesic.In other words, we can choose Ψ 0 so that Ψ = Ψ(ϕ) vanishes at a specific angle ϕ of the circular geodesic (or, equivalently, at a specific proper time τ ).If we do that, then the tidal tensor becomes diagonal at that specific angle.For instance, if we set ϕ 0 = 0 in (2.13) and choose Ψ 0 = 0 in (3.9), then the tidal forces in the radial and angular directions decouple at time τ = 0.In fact, whenever Ψ is an integer multiple of 2π, the tidal forces decouple according to: Hence, at the points where Ψ is an integer multiple of 2π, the tidal forces along the directions e 1, e 2, and e 3 are given, respectively, by κ 1 , κ 2 , and κ 3 .We plot such tidal forces as a function of the radius ρ 0 of the circular geodesic in Figures 5-7.For the sake of comparison, we have included in the left panels of each figure the corresponding tidal force when the MP metric describes an extremal Reissner-Nordström black hole (i.e., when b = 0).These figures tell us how the tidal forces change when we move from one circular geodesic to another.In contrast, Figures 2-4 exhibit how the tidal forces change along a given radial geodesic.Note that tidal forces associated with radial geodesics do not change as we move from one radial geodesic to another while keeping the distance ρ constant (since Equation (4.7) does not depend on the angle ϕ).
We see in Figure 5 that the tidal force in the radial direction is stretching for circular geodesics of sufficiently large radius, while for small radii, it is compressive.On the other hand, according to Figure 7, the behavior of the vertical tidal force is the opposite: stretching for small radii and compressive for large radii.From the explicit form of the potential U , given in Equation (2.2), we find that the vertical tidal force vanishes when the radius of the geodesic is The angular force shown in Figure 6, in contrast, is always compressive.We have also observed that, as the white dashed line in the right panels of Figures 5-7 is approached, the radial and vertical components of the tidal force diverge while the angular component remains finite.We close this section by noting that the eigenvalues of the matrix K are exactly κ 1 , κ 2 , and κ 3 .The corresponding eigenvectors are, respectively, Consequently, if we rotate back the vectors e µ 1 and e µ 2 by the angle Ψ, we obtain a new basis where the matrix K is diagonal.In other words, if we rotate the tetrad {e â} back to {λ â} at every point along the geodesic (consequently, undoing the rotation described in (3.5)), the tensor K âb becomes diagonal.However, since {λ â} is not parallel propagated along the geodesic (i.e., it is not a freely falling frame), the left-hand side of Equation (4.10) would not be the second derivative of the rotated vectors.

V. NUMERICAL SOLUTIONS OF THE GEODESIC DEVIATION EQUATION
We can also investigate tidal effects by determining the components of the deviation vector η â.We proceed as in [15][16][17][18][19][20][21][22][23] by directly solving the geodesic deviation Equation (4.1) with appropriate initial conditions.Being a second-order differential equation, the initial-value problem for Equation (4.1) requires an initial condition for both the deviation vector and its first-order derivative.We impose boundary conditions at the starting time τ = τ 0 , corresponding to the spacetime event x 0 = x(τ 0 ) = (t 0 , ρ 0 , ϕ 0 , 0).In our calculations, we have set ρ 0 = 10 and ϕ 0 = 0 (we also set Ψ 0 = 0 in Equation (3.9) and leave t 0 unspecified since the results do not depend on it).We consider two different initial conditions (denoted by IC-I and IC-II) to solve the geodesic deviation equation: IC-II : η â(τ 0 ) = 0, ηâ (τ 0 ) = 1. (5.2) Initial condition IC-I can be understood as representing the deviation between two geodesics that are initially parallel (in the sense that the rate of change of the their deviation vanishes at τ = τ 0 ).Initial condition IC-II, on the other hand, corresponds to two initially diverging geodesics that intercept at the initial time τ 0 (at the point x µ (τ 0 )).

A. Radial Geodesics
When analyzing radial geodesics, it is convenient to replace derivatives with respect to the proper time τ by derivatives with respect to the radial coordinate ρ.To accomplish this, we use the relation: which follows from Equation (2.8), to transform Equation (4.7) into: The equations for the radial and the angular components of the geodesic deviation can be integrated by quadratures: ) where K 1 , K 2 , K 3 , and K 4 are constants of integration that can be associated with the initial conditions for the differential equations.The component of the geodesic deviation along the z direction, on the other hand, cannot be cast in terms of simple integrals.
In order to solve the differential Equations (5.4)-(5.6),we impose the initial conditions (5.1) and (5.2).Note that, in terms of the radial coordinate, the initial condition for the derivative becomes: Regarding the vertical z direction, we see that, while the qualitative behavior of η 3 is sensitive to the specific values of b and E for IC-I, it is monotonically decreasing with ρ for IC-II.Figure 8 also shows that, when the energy of the geodesic increases, the radial deviation typically decreases as well.The same behavior with respect to the energy is observed for the angular η 2(ρ) and vertical η 3(ρ) components of the orthogonal connection, shown, respectively, in Figures 9 and 10, when IC-II is considered.On the other hand, for IC-I, we see that an increase of the energy typically increases the deviation in the angular and vertical directions.

B. Circular Geodesics
When considering circular geodesics, it is convenient to replace the derivatives with respect to the proper time τ by the derivatives with respect to the angular coordinate ϕ.To accomplish this, we use the relation: which follows from Equation (2.7).As a consequence, the derivatives in the left-hand side of Equation (4.10) are replaced by the derivatives with respect to the angle ϕ, while the matrix K is multiplied by ρ 4 0 U 4 (ρ 0 , 0)/L 2 .Analogously, in terms of ϕ, the initial condition for the derivatives in (5.1) and (5.2) becomes: In Figures 11 and 12, we show the numerical solutions of Equation (4.10) for both initial conditions IC-I (left panels) and IC-II (right panels). 1 In each panel, we indicate the separation b between the black holes and the radius of the corresponding circular geodesics (the associated values of energy E and angular momentum L are given in Figure 1).We have observed that the geodesic deviation in the directions e 1 and e 2 will diverge with the angle ϕ unless b is sufficiently large.Conversely, the deviation in the vertical direction e 3 will diverge unless b is sufficiently small.Such a behavior can be understood analytically for the vertical component η 3 (since η 1 and η 2 are coupled, a similar analysis for the radial and angular components is not possible).After plugging in the explicit form of U into Equation (4.10) and using Equation (2.11), we find that the differential equation for η 3 becomes (5.12)  where K 5 and K 6 are constants associated with the initial conditions.It is straightforward to see that the solutions will be exponentially growing if b is sufficiently large (b > ρ 0 / √ 2) and oscillatory otherwise (b < ρ 0 / √ 2).In fact, considering that the vertical tidal force, shown in Figure 7, vanishes exactly when b = ρ 0 / √ 2, we conclude that the deviation vector in the vertical direction will be bounded if the tidal force in the vertical direction is compressive.

VI. DISCUSSION
In this work, we investigated tidal forces acting on geodesics of the MP spacetime that describes two extremely charged black holes in equilibrium.We focused on radial geodesics (vanishing angular momentum L = 0) and circular geodesics (constant ρ = ρ 0 ) that live on the plane equidistant to two black holes of equal masses.In Section IV, we analyzed the radial, angular, and vertical components of the tidal forces, while in Section V, we analyzed the numerical solutions of the geodesic deviation equation for two classes of initial conditions.
In particular, we have seen that the tidal forces in the radial direction, for both radial and circular geodesics, can be either compressive or stretching: near the origin, they are compressive, and far away, they are stretching.A similar behavior has been observed for Reissner-Nordström black holes in Ref. [15].However, unlike Reissner-Nordström black holes, the tidal forces in the angular direction are always compressive for the geodesics of the MP that were analyzed.The tidal forces in the vertical direction, as in the case of the radial component, exhibit a sign change with respect to the distance.Nevertheless, the qualitative behavior is quite the opposite of their radial counterparts: for small radii, the force is stretching, and for large radii, it is compressive.Finally, we solved numerically the differential equations that determine the geodesic deviation vector as a function of the radial distance.We employed two types of initial conditions and analyzed the behavior of the corresponding solutions.

FIG. 1 .
FIG. 1. Contour plots of the energy E (left) and of the angular momentum L (right) associated with circular geodesics in the z = 0 plane, as a function of the distance b between the black holes and the radius ρ0 of the geodesic.The black region corresponds to the parameters for which circular geodesics are not allowed.

FIG. 2 .
FIG. 2. Tidal force along the radial direction for radial timelike geodesics confined in the z = 0 plane.(Left) Plots for several black hole separations b. (Right) Contour plot as a function of b and ρ (the dashed contour corresponds to the curve given by Equation (4.9) and indicates the transition from compression to stretching).

FIG. 4 .
FIG. 4. Tidal force along the vertical z direction for radial timelike geodesics confined in the z = 0 plane.(Left) Plots for several black hole separations b and energies E. (Right) Contour plot as a function of b and ρ for E = 1 (the red dashed contour indicates the transition from compression to stretching).

FIG. 5 .
FIG. 5. Tidal force along the radial direction for circular timelike geodesics confined in the z = 0 plane when Ψ is an integer multiple of 2π.(Left) Plots for several black hole separations b as a function of the radius ρ0 of the geodesic.(Right) Contour plot as a function of b and ρ0 (the dashed black contour indicates the transition from compression to stretching).

FIG. 6 .
FIG. 6. Tidal force along the angular ϕ direction for circular timelike geodesics confined in the z = 0 plane when Ψ is an integer multiple of 2π.(Left) Plots for several black hole separations b as a function of the radius ρ0 of the geodesic.(Right) Contour plot as a function of b and ρ0.

. 9 )
Additionally, due to Equation (2.8), the energy of the associated geodesics must satisfy the constraint |E| ≥ 1/U (ρ 0 , 0).The numerical solution of Equations (5.4)-(5.6),corresponding to the radial η 1(ρ), angular η 2(ρ), and vertical η 3(ρ) components of the orthogonal connection are shown in Figures8-10, respectively.In each plot, we indicate the separation b between the black holes and the energy E of the corresponding geodesic.The left panel in each figure corresponds to the initial condition IC-I, while the right panel corresponds to IC-II.We see that the radial deviation is largest at ρ = 0 and decreases as one moves away from the origin for both IC-I and IC-II.On the other hand, the angular deviation for IC-I increases with ρ, while the qualitative behavior for IC-II depends on the choices of b and E.

FIG. 8 .FIG. 9 .
FIG. 8. Component η 1 of the connection vector associated with radial geodesics in the z = 0 plane as a function of ρ for IC-I (left) and for IC-II (right).Each curve corresponds to a radial geodesic specified by the parameters b and E indicated in the panels.

FIG. 10 .
FIG. 10.Component η 3 of the connection vector associated with radial geodesics in the z = 0 plane as a function of ρ for IC-I (left) and for IC-II (right).Each curve corresponds to a radial geodesic specified by the parameters b and E indicated in the panels.

FIG. 12 .
FIG. 12. Components η 3 of the connection vector associated with circular geodesics in the z = 0 plane as a function of ϕ for IC-I (left) and for IC-II (right).Each curve corresponds to a circular geodesic specified by the parameters b and ρ0 indicated in the panels.