Light deflection around a spherical charged black hole to second order. Multivariate Pad\'e approximants

TL;DR

This paper derives higher-order light deflection terms around a charged black hole using perturbation methods and introduces multivariate Padé approximants that better fit numerical data than Taylor series.

Contribution

It develops second-order corrections for light deflection around a charged black hole and introduces multivariate Padé approximants for improved data fitting.

Findings

Second-order light deflection terms derived

Multivariate Padé approximants outperform Taylor series in data fitting

Enhanced understanding of light behavior near charged black holes

Abstract

From the Reissner-Nordstrom metric we obtain the higher-order terms for the deflection of light around a massive-charged black hole using the Lindstedt-Poincar\'e method to solve the equation of motion of a photon around the compact object. The corrections are performed developing the expansion in terms of $\epsilon=\frac{r_{c}}{b}$ and $\frac{\delta}{b}= \frac{Q^{2}}{6 \pi \epsilon_{0}Mc^{2} b}$. Additionally, from the perturbation expansion, we obtain multivariate Pad\'e approximants. We also show how these are a better fit for the numerical data than the original formal Taylor series.

Equations

$${\left(ds\right)}^2=\gamma {\left(cdt\right)}^2-{\gamma }^{-1}{\left(dr\right)}^2-r^2{\left(d\Omega \right)}^2$$

$$(g_{\mu\nu})=\left(\begin{array}[]{cccc}\gamma&0&0&0\\ 0&-\gamma^{-1}&0&0\\ 0&0&-r^{2}&0\\ 0&0&0&-r^{2}\sin^{2}\theta\end{array}\right).$$

$$\hat{t}=t\pm \frac{r_s}{c}ln\left|\frac{r}{r_s}-1\right|.$$

$${\left(ds\right)}^2=c^2\left(1-\frac{r_s}{r}\right){\left(d\hat{t}\right)}^2-\left(1+\frac{r_s}{r}\right){\left(dr\right)}^2\mp 2c\left(\frac{r_s}{r}\right)d\hat{t}dr-r^2{\left(d\Omega \right)}^2.$$

$$\hat{t}=t+\frac{r_s}{c}ln\left|\frac{r}{r_s}-1\right|$$

$$\hat{t}=t-\frac{r_s}{c}ln\left|\frac{r}{r_s}-1\right|$$

$$f\left(M\right)\equiv \frac{{\left(M_X\sin i\right)}^3}{{\left(M_X+M_C\right)}^2}={\left(a\sin i\right)}^3\frac{{\left(\frac{2\pi }{\tau }\right)}^2}{G}$$

$${\left(ds\right)}^2=\Delta {\left(cdt\right)}^2-{\Delta }^{-1}{\left(dr\right)}^2-r^2{\left(d\Omega \right)}^2$$

$$g_{\mu\nu}=\left(\begin{array}[]{cccc}\Delta&0&0&0\\ 0&-\Delta^{-1}&0&0\\ 0&0&-r^{2}&0\\ 0&0&0&-r^{2}\sin^{2}\theta\end{array}\right).$$

$$r_1=\frac{GM}{c^2}-{\left[{\left(\frac{GM}{c^2}\right)}^2-\frac{Q^{2}G}{4\pi {ϵ}_0{c}^4}\right]}^{\frac{1}{2}}$$

$$r_2=\frac{GM}{c^2}+{\left[{\left(\frac{GM}{c^2}\right)}^2-\frac{Q^{2}G}{4\pi {ϵ}_0{c}^4}\right]}^{\frac{1}{2}},$$

$$r_{\infty }=\frac{GM}{c^2}=\frac{Q}{c^2}\sqrt{\frac{G}{4\pi {ϵ}_0}}$$

$$\Delta {\left(cdt\right)}^2-{\Delta }^{-1}{\left(dr\right)}^2-r^2{\left(d\varphi \right)}^2=0.$$

$$\Delta c^2-{\Delta }^{-1}{\left(\frac{dr}{dt}\right)}^2-r^2{\left(\frac{d\varphi }{dt}\right)}^2=0.$$

$${\left(\frac{dr}{dt}\right)}^2=\Delta \left(\Delta c^2-r^2{\left(\frac{d\varphi }{dt}\right)}^2\right)$$

$$M{c}^2=M_0{c}^2+\frac{3}{5}\frac{k{Q}^2}{R}-\frac{3}{5}\frac{G{M}^2}{R},$$

$$M=\frac{1}{2aG}\left({\left(R^2{c}^4+4aG\left(M_{0}R{c}^2+ak{Q}^2\right)\right)}^{\frac{1}{2}}-R{c}^2\right),$$

$$\frac{dM}{dR}=\frac{a\left(G{M}^2-k{Q}^2\right)}{2aGMR+R^2{c}^2}$$

$$(ds{)}^2=c^2(d\tau {)}^2=\Delta c^2(dt{)}^2-{\Delta }^{-1}(dr{)}^2-r^2(d\varphi {)}^2.$$

$$L\left(x^{\mu },\frac{d{x}^{\mu }}{d\sigma }\right)=-g_{\alpha \beta }\left(x^{\mu }\right)\frac{d{x}^{\alpha }}{d\sigma }\frac{d{x}^{\beta }}{d\sigma }-2\frac{q}{m}{A}_{\alpha }\frac{d{x}^{\alpha }}{d\sigma }$$

$$\frac{\partial L}{\partial x^{\mu }}-\frac{d}{d\sigma }\left(\frac{\partial L}{\partial \left(\frac{d{x}^{\mu }}{d\sigma }\right)}\right)=0$$

$$\frac{d{u}_{\mu }}{d\sigma }=\frac{1}{2}\left({\partial }_{\mu }{g}_{\alpha \beta }\right)u^{\alpha }{u}^{\beta }+\frac{q}{m}{F}_{\mu \alpha }{u}^{\alpha }$$

$$F_{01}=-F_{10}=\frac{E^1}{c}=\frac{E\left(r\right)}{c},$$

$$\frac{d}{d\sigma }\left[\Delta c\frac{d\left(ct\right)}{d\sigma }+\frac{qQ}{4\pi {ϵ}_{0}mr}\right]=0$$

$$\frac{d}{d\sigma }\left[m{r}^2\frac{d\varphi }{d\sigma }\right]=0.$$

$$E^{\ast }=\Delta c^2\frac{dt}{d\tau }+\frac{qQ}{4\pi {ϵ}_{0}mr}$$

$$J=r^2\frac{d\varphi }{d\tau }.$$

$$c^2=\Delta c^2{\left(\frac{dt}{d\tau }\right)}^2-{\Delta }^{-1}{\left(\frac{dr}{d\varphi }\right)}^2{\left(\frac{d\varphi }{d\tau }\right)}^2-r^2{\left(\frac{d\varphi }{d\tau }\right)}^2.$$

$${\left(\frac{dr}{d\varphi }\right)}^2=-\Delta \left(\frac{c^2{r}^4}{J^2}+r^2\right)+\frac{r^4}{J^2{c}^2}{\left(E^{\ast }-\frac{qQ}{4\pi {ϵ}_{0}mr}\right)}^2$$

$$\frac{dr}{d\varphi }=-\frac{1}{U^2}\frac{dU}{d\varphi }$$

Full text

light deflection and charged and orbits and Padé

I Introduction

One of the most interesting predictions of the ”General Theory of Relativity” (GTR) is the existence of black holes, regions of space-time from which despising quantum effects, you can not escape. The Physics of a massive spherical black hole with no electric charge and rotational movement is described by the Schwarzschild metric, introduced by Karl Schwarzschild in 1916. For a spherical symmetric space-time with a mass $M$, the invariant interval is Kenyon ; Misner ; Marin2 :

[TABLE]

where $\left(d\Omega\right)^{2}=\left(d\theta\right)^{2}+sin^{2}\theta\left(d\phi\right)^{2}$ , with coordinates $x^{0}=ct$, $x^{1}=r$, $x^{2}=\theta$ and $x^{3}=\phi$. $\gamma=1-\frac{r_{s}}{r}$ ,being $r_{s}=\frac{2GM}{c^{2}}$ the Schwarzschild radius. The arc length $s$ satisfies the relation $ds^{2}=g_{\mu\nu}dx^{\mu}dx^{\nu}$, then, the covariant metric tensor is:

[TABLE]

We have two singularities in this metric. The first when $\gamma=0$ or $r=r_{s}$ is a mathematical singularity that can be removed by a convenient coordinate transformation like the one introduced by Eddington in 1924 or Finkelstein in1958 Kenyon ; Marin2 :

[TABLE]

With this coordinate transformation, the invariant interval can be written as:

[TABLE]

The first transformation in equation (3)

[TABLE]

describes a black hole, while the second one:

[TABLE]

represents what physicists call a ”white hole” emitting material from a singularity in $r=0$ toward space-time.

The other singularity $r=0$ is physical, so it can not be removed. In this singularity, all known physical laws fail, and the curvature of space-time is infinite. If one particle reaches the event horizon ($r=r_{s}$), it will eventually falls to the singularity $r=0$, and it will never escape from the black hole neglecting quantum effects like Hawking radiation.

There is good evidence of the existence of black holes; for example the X-ray source known as Cygnus X-1 Marin2 ; Ramesh . Cygnus X-1 emits X-rays and radio waves in an irregular way that contrasts with the regular emissions that are observed in a binary pulsar such as the PSR 1913 + 16. Cygnus X-1 is the companion of a blue supergiant star (HDE 226868) located at a distance of 2.5 kiloparsecs (8150 light years) from Earth. The mass of the blue supergiant is approximately 20 times that of the Sun. The wavelengths of the spectral lines of the supergiant is too cold to emit X-rays. Material ripped from HDE 226868 forms an accretion disk around Cygnus X-1 that is releasing an immense amount of energy. This energy heats the accretion disk to a temperature at which large amounts of X-rays are emitted. In conclusion, HDE 226868 and Cygnus X-1 form an X-ray binary with an orbital period of 5.6 days.

The mass of Cygnus X-1 can be calculated from the mass function:

[TABLE]

where $M_{X}$ is the mass of the X-ray source, $M_{C}$ is the mass of the blue supergiant, $i$ is the angle between the direction of observation and the normal to the plane of the orbit, $a$ is the length of the semimajor axis and $\tau$ is the orbital period of the X-ray binary. The mass of the source is estimated at approximately 15-16 times that of the Sun, making it a good candidate for a black hole. Other examples of binary systems that can host black holes are LMC X-3 and A0620-00 whose masses are estimated at nine times that of the Sun.

At the core of many galaxies such as the M-87, the active galaxy Centaurus A, the galaxy MCG-6-30-15, etc., there is evidence of the presence of black holes with masses that would range between $10^{6}$ and $3\times 10^{9}$ of times the mass of the Sun GaryHorowitz (supermassive black holes). Observations on the wavelengths corresponding to radio waves and X rays support this statement. These supermassive black holes may have played a very important role in the formation of galaxies in the early universe.

At the core of our galaxy, the Milky Way, there are signs of a swarm of black holes surrounding a supermassive black hole called Sagittarius A * (SgrA *) with a mass of about 4 million of times the mass of the Sun MJReid .

One of the largest black holes observed so far is in the heart of the quasar OJ287 Valtonen1 ; PauliPihajoki ; Valtonen2 at a distance of $3.5\times 10^{9}$ light years from Earth. It has a mass between $17\times 10^{9}$ and $18\times 10^{9}$ times that of the Sun. A smaller black hole with a mass of $100$ million times that of the Sun orbits around that supermassive black hole in an oval orbit with a period of twelve years. The advance of the perihelion due to the movement of translation of the smaller hole around the greater one is of 39 degrees in each orbit Valtonen1 .

The Chandra X-ray telescope Marin2 (name given in honor of the great theoretical astrophysicist Subrahmanyan Chandrasekhar who died in 1995), in orbit around the Earth, has recently obtained evidence of the existence of black holes in intermediate mass, between five hundred and twenty thousand times the mass of the Sun.

In general, depending on factors such as the temperature on the surface and the plasma physics in the solar corona, a star has a small electric charge. So, the contribution of the electromagnetic field in the metric is very small. However, when the star collapses (as nuclear fuel has run out) to become a black hole, one would expect that the entire charge will be ejected out of the star. For this reason, the black hole formed should be electrically neutral. However, there could be some other mechanism (which we do not know) by which some black holes were not electrically neutral. In any case, only observations made in the future will allow us to discard or not this possibility.

In a previous paper MarinRodriguez using the Schwarzschild metric we have obtained the higher-order terms for the deflection of light around a massive object (a black hole) using the Lindstedt-Poincaré method as well as diagonal Padé approximants from the perturbation expansion. In this paper, we obtain from the Reissner-Nordstrom metric, the second-order terms for the deflection of light around a massive-charged black hole using the Lindstedt-Poincaré method to solve the equation of motion of a photon around the compact object. The corrections are performed developing the expansion in terms of $\epsilon=\frac{r_{c}}{b}$ and $\frac{\delta}{b}=\frac{Q^{2}}{6\pi\epsilon_{0}Mc^{2}b}$. Higher order terms can be obtained following a procedure similar to the one we will illustrate in this paper. Also we obtain, in a way similar to the one we developed in the mentioned article MarinRodriguez , multivariate Padé approximants from the perturbation expansion. We make use of Padé approximants on our asymptotic series for the deviation angle to increase its region of validity, and to improve as we shall see, matches the qualitative behavior of the deflection angle.

It is worth mentioning that Padé polynomials were first used in Cosmology with excellent results (see references Christine ; Alejandro ; Salvatore for details).

II The Reissner-Nordström metric

The exact solution to the Einstein field equations for a static space-time with spherical symmetry with charge $Q$ and a mass $M$ in the center of the coordinate system is given by the Reissner-Nordstrom metric. The invariant interval is Misner ; tHooft :

[TABLE]

where $\Delta\equiv 1-\frac{2GM}{rc^{2}}+\frac{Q^{2}G}{4\pi\epsilon_{0}r^{2}c^{4}}$ and $\left(d\Omega\right)^{2}=\left(d\theta\right)^{2}+sin^{2}\theta\left(d\phi\right)^{2}$ , with coordinates $x^{0}=ct$, $x^{1}=r$, $x^{2}=\theta$ and $x^{3}=\phi$.

Then, the covariant metric tensor is:

[TABLE]

The surfaces with infinite red-shift are obtained from $g_{00}=\Delta=0$. If $\frac{GM}{c^{2}}<\frac{Q}{c^{2}}\sqrt{\frac{G}{4\pi\epsilon_{0}}}$ we don’t have such surfaces. In this case the singularity at $r=0$ is naked. However, this is forbidden by the hypothesis of the cosmic censorship Wald . If $\frac{GM}{c^{2}}>\frac{Q}{c^{2}}\sqrt{\frac{G}{4\pi\epsilon_{0}}}$ we will have two surfaces with infinite red-shift:

[TABLE]

[TABLE]

where $\epsilon_{0}=8.854\times 10^{-12}$ $\frac{Coul^{2}}{Nm^{2}}$ is the electric permitivity of the vacuum.

If $\frac{GM}{c^{2}}=\frac{Q}{c^{2}}\sqrt{\frac{G}{4\pi\epsilon_{0}}}$ we will have only one of such surfaces:

[TABLE]

For a light ray the invariant interval is $ds=0$. In the equatorial plane $\theta=\frac{\pi}{2}$, and then $d\theta=0$. Therefore we have:

[TABLE]

or

[TABLE]

From the last equation we can write:

[TABLE]

The event horizons are obtained setting $\frac{dr}{dt}=0$ for any $\theta,\phi$ and so $\Delta=0$. Then, the event horizons coincide with the surfaces with an infinite red-shift, being able to have two, one or no horizons of events.

Now, considering the equivalence between inertial and gravitational mass, the total energy of a spherical object (like a star) with rest mass energy $M_{0}c^{2}$, charge $Q$ (uniformly distributed in the volume of the sphere) and radius $R$ is:

[TABLE]

where $k=\frac{1}{4\pi\epsilon_{0}}$. The second and third terms in (16) represent the Coulomb energy and the gravitational binding energy, respectively. From this equation we can get:

[TABLE]

where $a=\frac{3}{5}$. We note that $M$ is a function of $R$. From equation (17), we can see that $lim_{R\rightarrow\infty}\left(M\left(R\right)\right)=M_{0}$ and $lim_{R\rightarrow 0}\left(M\left(R\right)\right)=\left(\frac{kQ^{2}}{G}\right)^{\frac{1}{2}}$, or $\frac{GM}{c^{2}}\rightarrow\frac{Q}{c^{2}}\sqrt{\frac{G}{4\pi\epsilon_{0}}}$, as expected (see figure 1). The derivative of $M$ with respect to $R$ give us:

[TABLE]

To have collapse, we need that $\frac{dM}{dR}>0$. Then collapse occurs only if: $GM^{2}>kQ^{2}$, or $\frac{GM}{c^{2}}>\frac{Q}{c^{2}}\sqrt{\frac{G}{4\pi\epsilon_{0}}}$, as we had expected.

III Equation of motion

Lets consider a gravitational source (like a black hole) of mass $M$ and charge $Q$ and a massive particle that is moving around the source. The motion of such particle is governed by the Reissner-Nordström metric. Because the orbital angular momentum must be constant, the orbital motion will be performed in a single plane. Then, we can analyze the motion of the particle in the equatorial plane $\theta=\pi/2$:

[TABLE]

The geodesic equation can be written in an alternative form using the Lagrangian

[TABLE]

where $\sigma$ is a parameter of the trajectory of the particle, which is usually taken to be the proper time, $\tau$ for a massive particle. $A_{\alpha}$ is the potential four-vector, $q$ is the charge of the particle and $m$ is its rest-mass. Using the Euler-Lagrange equations:

[TABLE]

we get the geodesic equation for the particle:

[TABLE]

where $u_{\mu}=\frac{dx_{\mu}}{d\sigma}$, and $F_{\mu\alpha}=\partial_{\mu}A_{\alpha}-\partial_{\alpha}A_{\mu}$ is an antisymmetric tensor. In the present case:

[TABLE]

the other $F_{\mu\nu}=0$ (we are considering a static electric field $\vec{E}=E\left(r\right)\vec{e_{r}}$ and we take $\vec{B}=\vec{0}$).

For the coordinates $ct$ ($\mu=0$) and $\phi$ ($\mu=3$) the geodesic equation (22) gives us, respectively :

[TABLE]

and

[TABLE]

Equations (24) and (27) with $\sigma=\tau$ imply that there are two constants of motion.

[TABLE]

and

[TABLE]

The first constant is the energy per unit mass, meanwhile the second is the angular momentum per unit mass.

Returning to equation (19), using $\frac{dr}{d\tau}=\frac{dr}{d\phi}\frac{d\phi}{d\tau}$ , it can be written as:

[TABLE]

From (26) and ( 27) we can express $\frac{dt}{d\tau}$ and $\frac{d\phi}{d\tau}$ in terms of $E^{*}$ and $J$ to obtain:

[TABLE]

This equation can be turned into an equation for $U(\phi)=\frac{1}{r(\phi)}$, noting that

[TABLE]

so we arrive at the following equation for $U(\phi)$:

[TABLE]

By taking the derivative of equation (29) with respect to $\phi$, we get the following differential equation for $U\left(\phi\right)$:

[TABLE]

[TABLE]

This differential equation can be separated into two differential equations for $U(\phi)$. The first one is the equation for a particle traveling inwards or outwards of the source of the gravitational field:

[TABLE]

the other differential equation, applicable for trajectories in which $U(\phi)$ is not constant with respect to $\phi$, is the following:

[TABLE]

This is the equation for the trajectory of a massive particle that travels around the source in the equatorial plane.

For a photon, the geodesic equation is MarinRodriguez

[TABLE]

where $\lambda$ is an afin parameter.

Consider a photon traveling in the equatorial plane ($\theta=\frac{\pi}{2}$) around the source. For the coordinates $ct$ ($\mu=0$) and $\phi$ ($\mu=3$) the geodesic equation (34) give us, respectively, the following constants along the trajectory of the photon around the source:

[TABLE]

[TABLE]

where $E^{**}$ has units of energy per unit mass and $J$ of angular momentum per unit mass (when $\lambda$ has units of time). The invariant interval can be written as:

[TABLE]

Inserting the definitions of $E^{**}$ and $J$ in (37) , with the change of variable $U(\phi)=\frac{1}{r(\phi)}$ , and following a similar procedure to the one used to deduce the orbit equation for a massive particle, we get:

[TABLE]

This differential equation can be separated into two differential equations for $U\left(\phi\right)$. The first one is the equation for a photon that travels directly into or out from the source:

[TABLE]

The other differential equation for trajectories in which $U\left(\phi\right)$ is not constant with respect to $\phi$, is:

[TABLE]

The last one is the equation for the trajectory of a photon that travels around the source in the equatorial plane.

IV Photon Surfaces

Let’s consider the differential equation for a photon traveling around a massive charged object like a charged black hole (equation 39). This equation has exact constant solutions, for the unstable circular orbits of a photon around the source. In fact, for circular orbits we have $\frac{dU}{d\phi}=0$ and $\frac{d^{2}U}{d\phi^{2}}=0$, and then, equation (39) for $U\neq 0$ give us:

[TABLE]

where $r_{1c}$ and $r_{2c}$ are the radii of the so-called photon spheres Claudel , being $r_{2c}$ the outermost photon sphere radius. Equation (40), also can be written as:

[TABLE]

where $r_{c}=\frac{3GM}{c^{2}}$ and $r_{Q}^{2}=\frac{Q^{2}G}{4\pi\epsilon_{0}c^{4}}$.

It can easily be shown that:

[TABLE]

where $r_{c}$ is the radius of the photon sphere for a photon traveling around a compact object (a black hole) in a space-time described by the Schwarzschild metric. In terms of the Schwarzschild radius, we can write $r_{c}=\frac{3}{2}r_{s}$.

In terms of $r_{c}$ the orbit equation (39) can be written as:

[TABLE]

Consider the initial conditions shown in Figure LABEL:fig:defl01. The smallest value of the $r-$coordinate in the trajectory, $r=b$, is taken such that the photon escapes the black hole, $b>r_{2c}$. However because the difference between $r_{2c}$ and $r_{c}$ is very small, we will rewrite Equation (43) in terms of $\epsilon=\frac{r_{c}}{b}<1$, which we will use as a non-dimensional small number for our following perturbative expansions. Note that by multiplying both sides of Equation (43) by $b$, and defining the non-dimensional trajectory parameter

[TABLE]

Equation (43), with the inclusion of the term $\epsilon=\frac{r_{c}}{b}$, then becomes a differential equation in $V(\phi)$:

[TABLE]

where $0<\epsilon=\frac{r_{c}}{b}<1$, $\delta\equiv\frac{Q^{2}}{6\pi\epsilon_{0}Mc^{2}}$ and with initial conditions given by

[TABLE]

Under these conditions, $V(\phi)$ is bounded such that

[TABLE]

V First-order solution for $V(\phi)$

To obtain a solution of Equation (45) we can write $V(\phi)$ as a power series in $\epsilon$:

[TABLE]

Plugging the expansion (48) into Equation (45) results in the following:

[TABLE]

[TABLE]

We can group the powers of $\epsilon$ in Equation (V):

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The initial conditions of $V(\phi)$, applied to the asymptotic expansion in Equation (48), imply (by grouping powers of $\epsilon$):

[TABLE]

[TABLE]

From these differential equations and initial conditions, we can readily obtain $V_{0}$ and $V_{1}$ iteratively:

[TABLE]

[TABLE]

[TABLE]

Thus, we obtain an equation for $V(\phi)$, per Equation (48):

[TABLE]

[TABLE]

According to the coordinate system shown in Figure 1, the photon goes through a total angular deflection of $2\alpha$. This corresponds to setting $V(\phi)=0$ for both $\phi=\pi/2+\alpha$ and $\phi=-\pi/2-\alpha$. From both of these conditions considering that $\alpha$ is very small, to first order in $\epsilon$ we get:

[TABLE]

The total deviation of the photon is then

[TABLE]

Note that $\Omega<\frac{4GM}{bc^{2}}$. Then, the angle of deviation is smaller than the one corresponding to a black hole without electric charge. This result agrees with the calculated by Shchigolev-Bezbatko Shchigolev and Briet-Hobill Briet .

VI Towards a second-order solution for $\Omega(\epsilon)$

We will now see how to obtain a second order solution for $\Omega$.

The term in Equation (V) that goes as $\phi sin(\phi)$ grows without bound. Then, if we naively include these kind of terms (called secular terms Bush ), our solution is no longer bounded. Thus, we have to eliminate any and all secular term that arises to arrive at a well-behaved solution for $V(\phi)$.

One method to do this, due to Lindstedt and Poincaré, is by solving the differential equation in the following strained coordinate Bush :

[TABLE]

Where the $\omega_{k}$ are constants to be determined. In terms of this new strained coordinate $\tilde{\phi}$, Equation (45) becomes

[TABLE]

Assuming an asymptotic expansion on $V(\tilde{\phi})$:

[TABLE]

Plugging the expansion(62) in Equation (61), we obtain:

[TABLE]

We can group the powers of $\epsilon$ in Equation (VI):

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

With some care due to the definitions of the scaled variable and its derivative, we arrive at initial conditions for the $V_{k}(\tilde{\phi})$ from the initial conditions of $V(\phi)$:

[TABLE]

[TABLE]

Solving the differential equation (64) with initial conditions (68), we arrive at the zeroth-order contribution to $V(\tilde{\phi})$:

[TABLE]

Introducing the value of $V_{0}$ given by (70) in equation (65), we get the differential equation:

[TABLE]

The general solution of the equation (71) is:

[TABLE]

[TABLE]

[TABLE]

The conditions (69) for $k=1$ give us:

[TABLE]

and

[TABLE]

and then:

[TABLE]

[TABLE]

[TABLE]

As we can check, a secular term has appeared for $V_{1}(\tilde{\phi})$ (a term that goes as $\tilde{\phi}\sin(\tilde{\phi})$). Then, equation (75) grows without limit. However, we can use the freedom in the definition of $\omega_{1}$ to eliminate this secular term by setting:

[TABLE]

then, the final expression for $V_{1}(\tilde{\phi})$ is:

[TABLE]

Thus, to first order $V_{1}\left(\tilde{\phi};\epsilon\right)$ is:

[TABLE]

We set up $\tilde{\phi}=\frac{\pi}{2}+\tilde{\alpha}$ in equation (78), such that $V\left(\frac{\pi}{2}+\tilde{\alpha}\right)=0$. We obtain the equation:

[TABLE]

To solve this equation let’s assume that $\sin(\tilde{\alpha})$ has the following expansion in $\epsilon$, with a leading order term of order $\epsilon^{1}$:

[TABLE]

where the $\chi_{k}$ are constants to be determined. Inserting this new expansion into Equation (79) leads to the following algebraic equation to first order:

[TABLE]

from which equating to zero the terms with $\epsilon^{1}$ we arrive at:

[TABLE]

Thus, to first order, $sin(\tilde{\alpha})$ is given by:

[TABLE]

To obtain $\tilde{\alpha}$, we employ the Taylor series of $\arcsin(x)$ around $x=0$:

[TABLE]

and obtain

[TABLE]

From the definition of the strained coordinate $\tilde{\phi}$ in (60), it is clear that:

[TABLE]

and to the first order we have:

[TABLE]

and replacing the value of $\tilde{\alpha}$ given by ( 85), we get:

[TABLE]

So, the total deflection angle is to first order:

[TABLE]

Introducing the values of $\epsilon$ and $\delta$, we finally get:

[TABLE]

which matches the value obtained in equation (59).

VII Second-order solution for $\Omega(\epsilon)$

Let´s return to equation (66), but now let´s rewrite it in the form:

[TABLE]

where

[TABLE]

Introducing the expressions of $V_{0}(\tilde{\phi})$ and $V_{1}(\tilde{\phi})$ given in (70) and (77), respectively, we have:

[TABLE]

[TABLE]

The homogeneus equation corresponding to (91 )is:

[TABLE]

with general solution:

[TABLE]

Let´s look for a particular solution of the inhomogeneous equation (91) of the form:

[TABLE]

satisfying:

[TABLE]

where $d_{1}^{{}^{\prime}}=\frac{d(d_{1}(\tilde{\phi}))}{d\tilde{\phi}}$ and $d_{2}^{{}^{\prime}}=\frac{d(d_{2}(\tilde{\phi}))}{d\tilde{\phi}}$. Introducing $V_{2p}(\tilde{\phi})$ in equation (91) we get:

[TABLE]

97 and 98 can be written in matrix form as:

[TABLE]

The solutions to this system of equations are:

[TABLE]

and

[TABLE]

Replacing the value of $V^{*}(\tilde{\phi})$ given by (VII) and performing the integrals, we get:

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Then, the general solution of equation (91) is:

[TABLE]

where $d_{1}(\tilde{\phi})$ and $d_{2}(\tilde{\phi})$ are given by (VII) and (VII), respectively. $c_{1}$ and $c_{2}$ are constants to be determined by the conditions: $V_{2}(\tilde{\phi}=0)=0$ and $\frac{dV_{2}}{d\tilde{\phi}}(\tilde{\phi}=0)=0$. The result is:

[TABLE]

[TABLE]

Introducing the values of $c_{1}$ and $c_{2}$ in (110) we have:

[TABLE]

[TABLE]

[TABLE]

or

[TABLE]

To eliminate the secular term in the last equation (the term containing $\tilde{\phi}\sin(\tilde{\phi})$), we set:

[TABLE]

getting

[TABLE]

[TABLE]

[TABLE]

Then, to second order we can write:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

We set up $\tilde{\phi}=\frac{\pi}{2}+\tilde{\alpha}$ in equation (VII), such that $V\left(\frac{\pi}{2}+\tilde{\alpha}\right)=0$, and use the expansion:

[TABLE]

where the $\chi_{k}$ are constants to be determined. The result is to second order in $\epsilon$:

[TABLE]

implying:

[TABLE]

and

[TABLE]

Thus, $\sin(\tilde{\alpha})$ is given to second order by

[TABLE]

$\tilde{\alpha}$ is obtained from the last equation using the Taylor series of $\arcsin(x)$ around $x=0$:

[TABLE]

to obtain:

[TABLE]

From the definition of the strained coordinate $\tilde{\phi}$ in (60) we have:

[TABLE]

From the last equation, and using the Taylor expansion of $\frac{1}{1+x}=\sum_{n=0}^{\infty}(-1)^{n}x^{n}$ around $x=0$, we get:

[TABLE]

The total deflection angle is then $\Omega=2\alpha$. Replacing the values of $\epsilon$ and $\delta$, we can finally write:

[TABLE]

[TABLE]

When $Q=0$, this result matches the one obtained in Bodenner ; Fischback ; Richter ; Epstein . A third-order solution for $\Omega(\epsilon)$ can be obtained from equation (67) following a procedure similar to the one used in this section. Note that getting a higher-order solution conserves the lower-order terms.

VIII Padé Approximants

The formal Taylor expansion of $\Omega$ around $\epsilon=0$ is:

[TABLE]

We can check easily that as $\epsilon\rightarrow 1$, $\Omega\left(\epsilon\right)\rightarrow\infty$.

As $b\rightarrow r_{c}$, the deflection angle must go to infinity ($\Omega\rightarrow\infty$). This is because the photons starts going into a closed orbit, the photon circular orbit. This implies that $\Omega(\epsilon)$ must diverge at $b=r_{c}$ ( in other words $\Omega(\epsilon)$ has a singularity at $\epsilon=1$). The problem is that a polynomial does not have asymptotes. If we want to recover the asymptote, then we need a rational function MarinRodriguez .

Rational functions of the form:

[TABLE]

are called the Padé approximants of a function $f(x)$, if the equation

[TABLE]

is true for all $k\in[0,M+N]$ Pade1 .

Also, if $f(x)$ can be written as a power series:

[TABLE]

then, the first $M+N$ terms of the Taylor expansion of $f^{[M/N]}$ around $x=0$ must match with the coefficients of equation (131).

To work with the expansion

[TABLE]

it is better to write it in terms of $\epsilon_{1}=\epsilon=\frac{3GM}{bc^{2}}$ and $\epsilon_{2}=\frac{GQ^{2}}{4\pi\varepsilon_{o}b^{2}c^{4}}$. Using this relations, the Taylor expansion of $\Omega$ would be:

[TABLE]

For example, using the second order Taylor polynomial of $\Omega$, one can recover the univariate $[1/1]$ Padé approximant of $\Omega$, for $\epsilon_{1}=\epsilon$, and $\epsilon_{2}=0$ (as in MarinRodriguez ):

[TABLE]

Now, if we want to calculate the first order Padé approximants for equation (133), including the charge term $\epsilon_{2}$, we would need a multivariate expansion.

Working with Padé approximants is very difficult in some cases, because there is not guarantee to obtain accurate results. However, some algorithms have been developed that facilitate the work Pade1 ; Pade2 . This is not the case for multivariate Padé approximants. First, it is difficult to chose the correct numerator and denominator, because not for all cases the system of equations have exact solutions. Also, one can obtain inaccurate results. Second, there is no optimal algorithm that gives us the correct number of equations we need to find the coefficients of the expansion Pade2 ; Pade3 .

Nevertheless, in the present article we use the algorithm showed in Pade4 , and then we analyse the results to find the best fit. Suppose the function $f(x,y)$ can be decomposed in power series as:

[TABLE]

The Padé multivariate problem consist in finding two polynomials of the form $p(x,y)=\sum_{k=0}^{M}A_{k}(x,y)$ and $q(x,y)=\sum_{k=0}^{N}B_{k}(x,y)$, where $A_{k}(x,y)$ and $B_{k}(x,y)$ are polynomials of order $k$; such that $p(x,y)$ and $q(x,y)$ satisfy:

[TABLE]

Or equivalently:

[TABLE]

For simplicity, we define the polynomials $C_{k}(x,y)=\sum_{i+j=k}c_{ij}x^{i}y^{j}$. This polynomials are: $C_{0}=c_{00}$, $C_{1}(x,y)=c_{10}x+c_{01}y$, $C_{2}(x,y)=c_{20}x^{2}+c_{11}xy+c_{02}y^{2}$, etc. Then, form condition (137) we can recover the following system of equations to find the polynomials $A_{k}$ and $B_{k}$:

[TABLE]

[TABLE]

And, the approximant can be written as

[TABLE]

From (133), one can retrieve the coefficients: $C_{0}=0$, $C_{1}(x,y)=\frac{4}{3}\epsilon_{1}-\frac{3\pi}{4}\epsilon_{2}$, $C_{2}(x,y)=\left(\frac{5\pi}{12}-\frac{4}{9}\right)\epsilon_{1}^{2}+\left(\frac{\pi}{2}-\frac{14}{3}\right)\epsilon_{1}\epsilon_{2}+\frac{57\pi}{64}\epsilon_{2}^{2}$. The problem is that setting $C_{0}=0$ retrieves problems like divisions by zero. To avoid that we set $C_{0}=\pi/2$, and then we only need to subtract $\pi/2$ from the final expression.

To see how the algorithm must be used, lets set $M=0$ and $N=1$ to calculate $\Omega^{[0/1]}$:

[TABLE]

Here we are already subtracting the $\pi/2$ term. From equations (138) and (139), one can find the following system of equations:

[TABLE]

[TABLE]

As it can be seen, with these equations we have one degree of freedom, so we will set $B_{o}=1$. Then:

[TABLE]

[TABLE]

Replacing these expressions in the equation of $\Omega^{[0/1]}$, and simplifying we have:

[TABLE]

Using similar steps, we get the following first order Padé approximants of $\Omega(\epsilon_{1},\epsilon_{2})$:

[TABLE]

[TABLE]

In these results, the $\pi/2$ term was already subtracted. Now, we need to test this approximants to analyse what is the best fit for the numerical points.

IX Analysis and numerical tests

Now we need to determine which Padé approximant is the best fit for $\Omega$. In figures (3), (4), (5) are shown the plots of $\Omega$ in function of $\epsilon=\epsilon_{1}$ for the different approximations, including Taylor polynomials and Padé approximants, and for different values of $\epsilon_{2}=\epsilon^{\prime}$.

As it can be seen, no Padé approximant seems to be a good fit for numerical values. This can be solved in the following way: As the numerical data are between the approximation [0/1] and the approximation [0/2], and also the approximation [0/2] does not have a singularity while the [0/1] has the singularity closest to the real one, we can approximate $\Omega$ with the average between [0/1] and [0/2] to find the best fit.

In figures (6), (7), (8) are shown the plots of

[TABLE]

for different values of $\epsilon_{2}=\epsilon^{\prime}$. As it can be seen, this term is a better fit of $\Omega$ for values of $r$ near to the singularity. Nevertheless, if it can be recovered more Taylor expansion terms, it is possible to calculate higher order Padé approximants.

X Conclusions

In this paper we have obtained from the Reissner-Nordström metric the second-order terms for the deflection of light around a massive-charged black hole using the Lindstedt-Poincaré method to solve the equation of motion of a photon around the compact object (equation (VII)). This method allows us to eliminate any and all secular terms that arise (like $\phi\sin(\phi)$) to get a well-behaved solution. The corrections are performed developing the expansion in terms of $\epsilon=\frac{r_{c}}{b}$ and $\frac{\delta}{b}=\frac{Q^{2}}{6\pi\epsilon_{0}Mc^{2}b}$. A third-order solution for $\Omega(\epsilon)$ (and higher order corrections) can be obtained from equation (67) following a procedure similar to the one used in section VII.

Also we have obtained multivariate Padé approximants from the perturbation expansion. By analysing the different Padé expansions, it was determined that the better fit for $\Omega$ is:

[TABLE]

were $\Omega^{[0/1]}$ and $\Omega^{[0/2]}$ are given by equations (142) and (143), respectively. This term is a better fit of $\Omega$ for values of $r$ near to the singularity (see figures (6), (7), (8)). However, if we can recover more Taylor expansion terms, it is possible to calculate higher order Padé approximants.

It is important to mention that Padé polynomials were first used in Cosmology with excellent results Christine ; Alejandro ; Salvatore .

We are convinced that this paper can be very useful for undergraduate and graduate students to learn the use of perturbative techniques as the Lindstedt-Poincaré method or the Padé approximants for solving problems not only within the framework of the General Theory of Relativity, but also in other fields of Physics.