Definite integral $\int_{R_0}^{R}\frac{dr}{r^2\sqrt{\frac{R_0-R_S}{R_0^3}-\frac{1}{r^2}\left(1-\frac{R_{s}}{r}\right)}}$

Question

In general relativity, null geodesics (in the unbounded case) can be written under the following form : $$\frac{d\varphi}{dr}=\frac{1}{r^2\sqrt{\frac{R_0-R_S}{R_0^3}-\frac{1}{r^2}\left(1-\frac{R_{s}}{r}\right)}}$$ with:

• $\left(r, \varphi\right)$ the polar coordinates of the photon
• $R_S$ the Schwarschild radius of the central object ($R_S\in\mathbb{R^{+}_{*}}$)
• $R_0$ the distance of closest approach ($R_0 > \frac{3\sqrt{3}}{2}R_S$)

Consequently, to compute the exact trajectory up to a radius $R$ (with $R > R_0$), one can evaluate: $$\varphi\left(R\right)=\int_{R_0}^{R}\frac{dr}{r^2\sqrt{\frac{R_0-R_S}{R_0^3}-\frac{1}{r^2}\left(1-\frac{R_{s}}{r}\right)}}$$

And now comes my question : is there an analytical formula (in terms of special functions for example) corresponding to this integral ? (I would like to compute this integral numerically and an expression in terms of special functions would help a lot).

Answer 1

The integral can be expressed in terms of inverse Weierstrass's elliptic functions and hence in inverse Jacobi elliptic functions. The other answer has mentioned the incomplete elliptic integral of the first kind. That integral is essentially one of the inverse Jacobi elliptic functions in disguise.

I'm going to give the expression in terms of inverse Weierstrass elliptic functions only. Unlike the expression in Jacobi elliptic functions/elliptic integrals, this part is much easier to derive.

Let $\alpha = \frac{R_s}{R_0} = \beta + \frac13$, $u = \frac{R_s}{r} = p + \frac13$, the integral can be rewritten as

$$\varphi(R) = \int_{\frac{R_s}{R}}^\alpha \frac{du}{\sqrt{\alpha^2(1-\alpha) - u^2(1-u)}} = \int_{\frac{R_s}{R}-\frac13}^{\beta}\frac{dp}{\sqrt{p^3 - \frac{p}{3} - ( \beta^3 - \frac{\beta}{3} )}} $$ Let $\wp(z)$ be the Weierstrass elliptic function associated with the ODE

$$\wp'(z)^2 = 4\wp(z)^3 - g_2 \wp(z) - g_3\quad\text{ where }\quad \begin{cases}g_2 &= \frac43\\g_3 &= 4(\beta^3 - \frac{\beta}{3})\end{cases}$$

and substitute $p$ by $\wp(z)$, we have

$$\varphi(R) = 2 \int_{\frac{R_s}{R}-\frac13}^\beta \frac{d\wp (z)}{\sqrt{4\wp(z)^3 - g_2 \wp(z) - g_3}}$$

If one change variable to $z$, the factor $\frac{d\wp}{\sqrt{4\wp^3 - g_2\wp - g_3}}$ in above integrand reduces to either $+dz$ or $-dz$.

Let $\omega$ and $\omega'$ be the half periods of $\wp(z)$ with $\omega \in (0,\infty)$ and $\Im\omega' > 0$. Let $\omega_c = \omega + \omega'$. We will adopt the convention as $\rho$ decreases from $+\infty$ to $-\infty$ along the real axis, $\wp^{-1}(\rho)$ will move along the path

$$0 \to \omega \to \omega_c = \omega + \omega' \to \omega' \to 0$$

For the parameter range $R > R_0 > \frac{3\sqrt{3}}{2} R_s$ we are interested in, $\wp^{-1}(\beta) = \omega_c$ and in general $\wp^{-1}(p)$ will lie on the line segment joining $\omega'$ and $\omega + \omega'$. On this line segment, above factor reduces to $dz$ and we get

$$\varphi(R) = 2 \int_{\wp^{-1}\left(\frac{R_s}{R}-\frac13\right)}^{\wp^{-1}(\beta)} dz = 2 \left[ \omega_c - \wp^{-1}\left(\frac{R_s}{R} - \frac13\right) \right] $$ Computationally, this is not that convenient to use. This is because for the range of $r$ we care, $\wp^{-1}(p)= \wp^{-1}\left(\frac{R_s}{r} - \frac13\right)$ is a complex number! Using following addition formula for $\wp(z)$:

$$\wp(z+\omega_c) + \wp(z) + \wp(\omega_c) = \frac14 \left(\frac{\wp'(z) - \wp'(\omega_c)}{\wp(z) - \wp(\omega_c)}\right)^2$$ We find

$$ \wp(\omega_c-z) = \wp(z+\omega_c) = \frac{\left(p^3-\frac{p}{3}\right) - \left(\beta^3 - \frac{\beta}{3}\right)}{( (p-\beta)^2} - (p + \beta) = \beta + \frac{9\beta^2-1}{3(p-\beta)} $$ and hence $$\varphi(R) = 2\wp^{-1}\left( \beta + \frac{9\beta^2-1}{3(p-\beta)} \right) = 2\wp^{-1}\left(\alpha-\frac13 + \frac{\alpha(3\alpha-2)}{\frac{R_s}{R} - \alpha}\right) \tag{*1}$$ To demonstrate how this work, let us consider the case $(R_s, R_0, R ) = (1,3,4)$ as an example. We have $\alpha = \frac13$ and $\varphi(R)$ becomes

$$\begin{align} \varphi(R) = & \int_3^4 \frac{dr}{r^2\sqrt{\frac{2}{27} - \frac{1}{r^2}\left(1 - \frac{1}{r}\right)}}\\ \approx & 1.00210163826249415497545279509980330387176483490377 \end{align} $$ The numerical value is evaluated using following command on WA (Wolfram Alpha).

$\quad(a)\quad\quad$ N[Integrate[1/(r^2*Sqrt[2/27-1/r^2*(1-1/r)]),{r,3,4}],50]

For this example, $(g_2, g_3) = (4/3,0)$. RHS of $(*1)$ becomes $\wp^{-1}(4)$ and once again, we can evaluate it on WA using another command

$\quad(b)\quad\quad$ N[-2*InverseWeierstrassP[ 4, {4/3,0}],50]

and get $$2\wp^{-1}(4) \approx 1.0021016382624941549754898784512602241681634415061$$

To those with the sharp eyes, you will notice two things.

1. the numerical output of these two commands are very close but not the same. The only thing I can said is WA label the output for $(a)$ as decimal approximation while the output for $(b)$ as result.

2. There is a mysterious minus sign in command $(b)$. It turns out WA uses a different convention for $\wp^{-1}(p)$ than me. For large and possible $p$, InverseWeierstrass[p,{g_2,g_3}] returns a negative instead of a positive number and hence the need of that extra minus sign.

Update

About expressing everything in terms of inverse Jacobi elliptic functions. I finally dig up my reference books and find something useful. Let $e_1, e_2, e_3$ be the 3 roots of the polynomials of $4p^3 - g_2 p - g_3 = 0$. For concreteness, let us choose

$$e_1 = \wp(\omega),\quad e_2 = \wp(\omega_c)\quad\text{ and }\quad e_3 = \wp(\omega')$$

We have

$$\wp'(z)^2 = 4 \wp(z)^3 - g_2 \wp(z) - g_3 = 4(\wp(z) - e_1)(\wp(z) - e_2)(\wp(z) - e_3)$$

When $0 < \alpha < \frac23$, one can show that these 3 roots are real and distinct. Furthermore, they satisfy $e_1 > e_2 = \beta > e_3$. If one define a function $\text{sn}(\cdot)$ by following relation

$$\wp(z) = e_3 + \frac{e_1-e_3}{\text{sn}(z\sqrt{e_1-e_3})^2} \quad\iff\quad \text{sn}(z) = \sqrt{\frac{e_1 - e_3}{\wp\left(\frac{z}{\sqrt{e_1-e_3}}\right) - e_3}}$$

One can show that $\text{sn}(z)$ satisfy the ODE for the Legendre form of Jacobi elliptic sine function:

$$\text{sn}'(z)^2 = (1 - \text{sn}(z)^2)(1 - m\,\text{sn}(z)^2)\quad\text{ where }\quad m = \frac{e_2 - e_3}{e_1-e_3}$$

One can use this relation to convert the expression from inverse Weierstrass elliptic functions to inverse Jacobi elliptic functions.

Answer 2

There is not an analytic pretty closed form to this one, but we can still try and learn a tiny bit more about it. However it is related to Elliptic integrals of the first kind. Note, Jacobi Elliptic functions are inverses for the elliptic integrals.

If you would like to compute it numerically still, let me know, I did so and have codes in mathematica from using NIntegrate. I also computed it using just Integrate in mathematica as well and have this code.

I hope this helps, writing your integral we have $$\varphi\left(R\right)=\int_{R_0}^{R}\frac{dr}{r^2\sqrt{\frac{R_0-R_S}{R_0^3}-\frac{1}{r^2}\left(1-\frac{R_{s}}{r}\right)}}$$ We define $a\equiv (R_0-R_s)/R^3_0$, $b\equiv R_s$ and write $$ \varphi\left(R\right)=\int_{R_0}^{R}\frac{dr}{r^2\sqrt{a-\frac{1}{r^2}\left(1-\frac{b}{r}\right)}}$$

Now we factor our some constants under square root and do algebra to obtain $$ \varphi(r)=a^{-1/2}\int_{R_0}^R \frac{dr}{r^2\sqrt{1-\frac{1}{r^2}\left(a^{-1}-\frac{b}{ar}\right)}}=a^{-1/2}\int_{R_0}^R \frac{dr}{r \sqrt{r^2+\frac{b}{ar}-a^{-1}}}. $$ Now we can briefly look at the indefinite integral of this $$ \varphi(r)=a^{-1/2}\int \frac{dr}{r \sqrt{r^2+\frac{b}{ar}-a^{-1}}} $$ This integral does not have a closed form however can be written in a not so pretty Elliptic Integral of the First Kind (Elliptic-F) Function. So to answer your question, it is expressible in terms of elliptical integrals. It is something like $I\propto r F(\arcsin(...)) $. Mathematica could not compute much else unless you want to use NIntegrate.

Elliptic integrals of first kind in Jacobi's form is given by $$ F(x,k)=\int_0^x \frac{dt}{\sqrt{(1-t^2)(1-k^2t^2)}}, $$ where the usual form you may know is $$ F(\zeta,k)=\int_0^\zeta \frac{d\theta}{\sqrt{1-k^2\sin^2\theta}}, $$ they are related by a change of variables $t=\sin \theta$ and k is related to the eccentricity of the orbit. Sorry I can't find a closed form for your post, but I hope elliptic F functions will do. This is definitely also why they do not put the solution in any book, because it is not so clean.