eq}=2M$, has a global maximum of $(3\sqrt{3}M)^{-1}$ at $r_{\rm eq}=3M$ and then falls back to zero asymptotically.[]{data-label="E-profile"}](E-profile.eps){width="\columnwidth"}

In order to integrate this equation, notice the cubic polynomial inside the square root in the denominator of the rightmost side of equation (\[dtheta/dr\]). Its graph is everywhere concave upward, having two real roots above the horizon, one below and one above $3M$. All our photons start moving perpendicularly from the equatorial plane at $r>3M$ and escape to infinity, which means that the integration is performed just from the outer root $r_{\rm eq}$ (turning point of radial motion) up to a desired radius $r>3M$. Assuming without loss of generality that the photons start moving “upwards" from the equatorial plane (i.e., that $\epsilon^\theta\!=\!-1$ initially), one finds the solution formula obtained by Darwin [@Darwin-59] (see also equation (260) in Chapter 3 of @Chandrasekhar)[^3] $$\begin{aligned}
  -\epsilon^\theta \theta(r)
  &= \frac{\pi}{2}
     -\frac{2\,\sqrt{r_{\rm eq}}}{[(r_{\rm eq}-2M)(r_{\rm eq}+6M)]^{1/4}}
      \left[K(k)-F(\chi,k)\right] \label{theta(r)} \\
  &= \frac{\pi}{2}
     -\frac{2\,\sqrt{r_{\rm eq}}}{[(r_{\rm eq}-2M)(r_{\rm eq}+6M)]^{1/4}}\;
      F(\chi',k) \,,  \label{theta(r)'}\end{aligned}$$ where $F(\chi,k)=\int_0^\chi\frac{{\rm d}\alpha}{\sqrt{1-k^2\sin^2\alpha}}$ is the elliptic integral of the 1st kind, with amplitude $\chi$ and modulus $k$ given by $$\begin{aligned}
  \sin^2\chi
   &= 1-\frac{1}{k^2}\,
        \frac{2M\left(1-\frac{r_{\rm eq}}{r}\right)}{\sqrt{(r_{\rm eq}-2M)(r_{\rm eq}+6M)}}
    \;, \\
  2k^2
   &= 1-\frac{r_{\rm eq}-6M}{\sqrt{(r_{\rm eq}-2M)(r_{\rm eq}+6M)}} \;,  \label{2k2}\end{aligned}$$ and $K(k)=F(\pi/2,k)$ is its 