In this paper, we study the geodesic motion in spherically symmetric electro-vacuum Euclidean solutions of the Einstein equation. There are two kinds of such solutions: the Euclidean Reissner–Norström (ERN) metrics, and the Bertotti–Robinson-like (BR) metrics, the latter having constant Kretschmann scalar.

First, we derive the motion equations for the ERN spacetime and we generalize the results of Battista–Esposito, showing that all orbits in as ERN spacetime are unbounded if and only if it has an event horizon. We also obtain the Weierstrass form of the polar radial motion, providing an efficient tool for numerical computations.

We then study the angular deflection of orbits in the Euclidean Schwarzschild spacetime which, in contrast to the Lorentzian background, can be either positive or negative. We observe the presence of a null and a maximal deflection rings for particles with velocity at infinity $v>1$ and we give approximate values for their size when $v\gtrsim1$.

For BR spacetimes, we obtain analytic solutions for the radial motion in proper length, involving (hyperbolic) trigonometric functions and we deduce that orbits either exponentially go to the singularity or are periodic.

Finally, we apply the previous results and use algorithms related to Weierstrass' elliptic functions to produce a Python code to plot orbits of the spacetimes ERN and BR, and draw “shadows” of the first ones, as it was already done before for classical black holes.

§ INTRODUCTION AND MOTIVATION

Inspired by <cit.>, a gravitational instanton can be defined as a Euclidean (i.e. a Riemannian) metric on a 4-dimensional manifold which satisfies the Einstein field equation. Such metrics were first introduced in quantum gravity by Hartle and Hawking <cit.> in order to make some path integral converge, hence defining the so-called Hartle–Hawking propagator. See <cit.> for more details.

The simplest example of a gravitational instanton is the Euclidean Schwarzschild metric <cit.>, given in Schwarzchild coordinates $(\tau,r,\theta,\phi)$ on $\mathcal{M}:=\R\times]2M;+\infty[\times\Sph^2\simeq\R^2\times\Sph^2$ by

\begin{equation}\label{ES}\tag{ES}
{\rm d}s^2=\left(1-\frac{2M}{r}\right){\rm d}\tau^2+\left(1-\frac{2M}{r}\right)^{-1}{\rm d}r^2+r^2({\rm d}\theta^2+\sin^2\theta{\rm d}\phi^2),
\end{equation}

($M\ge0$ being the mass of the central body) and the link with the Lorentzian Schwarzschild metric is given by setting the Euclidean time $\tau=it$, with $t$ being the (Lorentzian) coordinate time. As for any metric, it is natural to study the geodesic motion associated to it. Regarding the metric (<ref>) above, this question was addressed thoroughly in <cit.>.

As suggested in <cit.>, one can also look for the Reissner–Nordström analogue of the Euclidean Schwarzschild solution; a metric that was used in <cit.> and <cit.>, for instance. Assuming the central body has an electric charge $Q\in\R$, the Euclidean Reissner–Nordström (ERN) metric is given by

\begin{equation}\label{ERN}\tag{ERN}
{\rm d}s^2=\left(1-\frac{2M}{r}+\frac{Q^2}{r^2}\right){\rm d}\tau^2+\left(1-\frac{2M}{r}+\frac{Q^2}{r^2}\right)^{-1}{\rm d}r^2+r^2({\rm d}\theta^2+\sin^2\theta{\rm d}\phi^2).
\end{equation}

In this paper, we generalize the approach of <cit.> to this metric and study the motion of a test particle in the ERN spacetime. Specifically, we prove that the method of <cit.> still applies to the ERN metric, hence obtaining analytic solutions for (non-purely radial) geodesics in terms of Weierstrass' elliptic functions. As we will see, one of the remarkable results obtained in <cit.> extends to the ERN spacetime with horizon (i.e. such that $Q^2<M^2$), namely the fact that the energy of an exterior geodesic is confined in the open interval $]-1,1[$. In particular, no elliptic-like geodesics exist in the sub-extremal case. This fails in the super-charged case $Q^2>M^2$, where arbitrary high energy is allowed and attained by a circular geodesic. As mentioned in <cit.>, these facts show that these Euclidean metrics present substantial differences in their dynamics, when compared to their usual Lorentzian avatars.

Moreover, we will see how the polar motion equation simplifies in the Schwarzschild case $Q=0$ and we retrieve the results from <cit.> using only the elementary geometric properties of the real elliptic curve describing the phase portrait in (affinely transformed) Binet variable.

Another remarkable dynamical distinction between Lorentzian and Euclidean Schwarzschild spacetimes relies in the gravitational deflection of orbits. As is well-known, given an orbit coming from and to infinity, the deflection angle $\delta\phi$ (in the motion plane) between its two asymptotic directions is always positive in Schwarzschild geometry, for photons as well as for massive particles. However, in Euclidean Schwarzschild geometry, the deflection angle can be positive or negative, depending on the velocity at infinity $v$ and the perihelion $r_{\rm min}$ of the orbit. More precisely, we observe that, at fixed $v<1$, we have $\delta\phi<0$ for all values of $r_{\rm min}$, while for $v>1$, the deflection $\delta\phi$ vanishes (resp. is maximal positive) at some perihelion $r_{\rm min}=\rho_0$ (resp. $r_{\rm min}=\rho_{\rm max}$). In (<ref>) and (<ref>), we give approximate values for $\rho_0$ and $\rho_{\rm max}$ when $v\gtrsim1$. The existence of $\rho_{\rm max}$ gives rise to a visible maximal deflection ring in the shadow of such a spacetime. Moreover, the fact the $\delta\phi<0$ at small perihelia indicates that particles are repelled by the central mass and thus the horizon becomes invisible to the observer, see Figures <ref> and <ref>.

The only other possible type of spherically symmetric instanton is given by a Bertotti–Robinson-like metric, whose line element has the form

\begin{equation}\label{BRr}\tag{BR}
{\rm d}s^2=Q^2\left(\frac{1-2mr+q^2r^2}{r^2}{\rm d}\tau^2+\frac{{\rm d}r^2}{r^2(1-2mr+q^2r^2)}+{\rm d}\theta^2+\sin^2\theta{\rm d}\phi^2\right),
\end{equation}

where $Q\ne0$ and $m,q\in\R$ are some constants. This is the Euclidean analogue of the general Bertotti–Robinson electro-vacuum Lorentzian solution. As in the Lorentzian case, this metric essentially differs from (<ref>) is this sense that its Kretschmann invariant is constant. To the knowledge of the author, the dynamics of this solution doesn't appear in the literature. This is treated in <ref>, where we provide a full analytic solution of the geodesic equation, in terms of (hyperbolic) trigonometric functions; see the Theorem <ref> and the Corollary <ref>.

Finally, the efficient algorithms available to approximate Weierstrass' elliptic functions <cit.> are used to produce a Python code[available at <https://github.com/arthur-garnier/euclidean_orbits_and_shadows.git>], designed to draw orbits in the two types of instantons discussed here, as well as to obtain the “shadow” of an ERN space by ray-tracing, as it was already done for black holes (see <cit.>). Since there are no null geodesics in Euclidean geometry, photons are replaced by particles with a velocity at infinity that should be provided by the user. Conformally to what was mentioned above concerning the deflection angle, we observe the presence of a maximal deflection ring when $v>1$ and we notice that the horizon is not visible. In particular, the optical difference between the cases $Q^2\lesssim M^2$ and $Q^2\gtrsim M^2$ is not as obvious as in the Lorentzian background. We still observe that the size of the maximal deflection ring diminishes as the charge increases, see Figure <ref>.

The layout of the paper is as follows: first, we prove that the metrics (<ref>) and (<ref>) are the only spherically symmetric solutions of the Einstein–Maxwell field equation with complex vector potential $A_\mu=-iQr^{-1}{\rm d}\tau$. Then, we derive the motion equations and the motion constants for the metric (<ref>), and we prove that the energy $E$ of a geodesic satisfies $E^2<1$ when $Q^2<M^2$, as mentioned above. We then obtain the Weierstrass equation from of the polar radial motion equation and we investigate the particular case where $Q=0$. Next, we study the gravitational deflection in Euclidean Schwarzschild spacetime and provide the aforementioned approximations for the deflection angle, as well as for the null and maximal deflection rings. Concerning the Bertotti–Robinson family (<ref>), we derive the motion equations and obtain analytic solutions with (hyperbolic) trigonometric functions. Finally, we give some explanations on how the Python code is constructed and we finish with some figures illustrating our results and programs.

§ THE TWO TYPES OF SPHERICALLY SYMMETRIC ELECTRO-VACUUM INSTANTONS

We systematically use Stoney units where $G=c=4\pi\epsilon_0=1$. Let $(M,Q)\in\R_+\times\R$ and consider the numbers $r_+,r_-$ defined by

\[r_{\pm}:=\left\{\begin{array}{ll}M\pm\sqrt{M^2-Q^2} & \text{if }Q^2<M^2, \\ 0 & \text{otherwise.}\end{array}\right.\]

Let also $\mathcal{M}:=\R\times ]r_+,+\infty[\times\Sph^2\simeq\R^2\times\Sph^2$, with coordinates $x^\mu=(\tau,r,\theta,\phi)$, the pair $(\theta,\phi)$ describing spherical coordinates on $\Sph^2$. Since we work on the Euclidean section, we restrict our study to the open subset $r>r_+$, just as in <cit.>. This is a reasonable restriction, as we are interested by the exterior region of the spacetime.

We denote by ${\rm d}\Omega^2:={\rm d}\theta^2+\sin^2\theta{\rm d}\phi^2$ the usual round metric on $\Sph^2$ and we have the following result:

Let ${\rm d}s^2=g_{\mu\nu}{\rm d}x^\mu{\rm d}x^\nu$ be a spherically symmetric solution of the Einstein–Maxwell equation with complex vector potential

\[A_\mu=-iQr^{-1}{\rm d}\tau,\]

defined for $r\gg0$.

If the Kretschmann invariant $K=R^{\alpha\beta\mu\nu}R_{\alpha\beta\mu\nu}$ associated to $g_{\mu\nu}$ is independent of $r$, then there are constants $m,q\in\R$ such that the metric takes the Bertotti–Robinson form

\[{\rm d}s^2=Q^2\left[\frac{1-2mr+q^2r^2}{r^2}{\rm d}\tau^2+\frac{{\rm d}r^2}{r^2(1-2mr+q^2r^2)}+{\rm d}\Omega^2\right]\]

in which case the Kretschmann invariant is $K=8Q^{-4}$.

Otherwise, there are coordinate transformations of the form $\tilde{\tau}=C\tau$ and $R=r/(Ar+B)$ (with $A\in\R$ and $B,C\in\R^*$) as well as constants $\tilde{M},\tilde{Q}\in\R$ such that $g_{\mu\nu}$ takes the Reissner–Nordström form

\[{\rm d}s^2=\left(1-\frac{2\tilde{M}}{R}+\frac{\tilde{Q}^2}{R^2}\right){\rm d}{\tilde{\tau}}^2+\left(1-\frac{2\tilde{M}}{R}+\frac{\tilde{Q}^2}{R^2}\right)^{-1}{\rm d}R^2+R^2{\rm d}\Omega^2,\]

whose Kretschmann invariant is $K=8R^{-8}(6\tilde{M}^2R^2-12\tilde{M}\tilde{Q}^2R+7\tilde{Q}^4)$.

Let ${\rm d}s^2=g_{\mu\nu}{\rm d}x^\mu{\rm d}x^\nu$ be a spherically symmetric solution of the field equation. As proved in <cit.>, a Lorentzian spherically symmetric solution of the electro-vacuum Einstein equations is necessarily static and the proof can be adapted verbatim to the Euclidean case. Therefore, the metric is diagonal: ${\rm d}s^2=g_{\mu\mu}({\rm d}x^\mu)^2$ and moreover, the functions $g_{\tau\tau}$ and $g_{rr}$ only depend on $r$. Furthermore, since the restriction to a hypersurface with constant coordinates $(\tau,r)$ must be an $SO(3)$-invariant metric on $\Sph^2$, there is a positive smooth function $\rho : ]r_+,+\infty[\longto\R^*_+$ such that $g_{\theta\theta}{\rm d}\theta^2+g_{\phi\phi}{\rm d}\phi^2=\rho(r)r^2{\rm d}\Omega^2$. Thus the metric may be reduced to the form

\[{\rm d}s^2=\rho(r)\left[u(r){\rm d}\tau^2+v(r){\rm d}r^2+r^2{\rm d}\Omega^2\right].\]

The Christoffel symbols ${\Gamma^\alpha}_{\mu\nu}=\tfrac12g^{\alpha\beta}(g_{\beta\mu,\nu}+g_{\beta\nu,\mu}-g_{\mu\nu,\beta})$ of this metric read

\[\left\{\begin{array}{ll}
{\Gamma^\tau}_{\mu\nu}=\frac{(u\rho)'}{2u\rho}\left(\begin{smallmatrix}0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0\end{smallmatrix}\right), & {\Gamma^r}_{\mu\nu}=\frac{1}{2v\rho}\left(\begin{smallmatrix}-(u\rho)' & 0 & 0 & 0 \\ 0 & (v\rho)' & 0 & 0 \\ 0 & 0 & -r(r\rho'+2\rho) & 0 \\ 0 & 0 & 0 & -r\sin^2\theta(r\rho'+2\rho)\end{smallmatrix}\right), \\[2em]
{\Gamma^\theta}_{\mu\nu}=\left(\begin{smallmatrix}0 & 0 & 0 & 0 \\ 0 & 0 & \frac{r\rho'+2\rho}{2r\rho} & 0 \\ 0 & \frac{r\rho'+2\rho}{2r\rho} & 0 & 0 \\ 0 & 0 & 0 & -\sin\theta\cos\theta\end{smallmatrix}\right), & {\Gamma^\phi}_{\mu\nu}=\left(\begin{smallmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{r\rho'+2\rho}{2r\rho} \\ 0 & 0 & 0 & \cotan\theta \\ 0 & \frac{r\rho'+2\rho}{2r\rho} & \cotan\theta & 0\end{smallmatrix}\right),\end{array}\right.\]

where we have dropped the dependence in the variable $r$ for simplicity.

On the other hand, the electro-magnetic field tensor associated to the potential $A_\mu$ is given by

\[F_{\mu\nu}=A_{\mu,\nu}-A_{\nu,\mu}=\frac{-iQ}{r^2}\left(\begin{smallmatrix}0 & 1 & 0 & 0 \\ -1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0\end{smallmatrix}\right)\]

so that the Maxwell equations

\[0={F^{\mu\nu}}_{;\mu}={F^{\mu\nu}}_{,\mu}+{\Gamma^\mu}_{\mu\lambda}F^{\lambda\nu}+{\Gamma^\nu}_{\mu\lambda}F^{\mu\lambda}=\frac{1}{\sqrt{\det g}}\left(\sqrt{\det g} F^{\mu\nu}\right)_{,\mu}\]

are satisfied for all $\nu\ne\tau$ and we have ${F^{\mu\tau}}_{;\mu}={F^{r\tau}}_{,r}+{\Gamma^\mu}_{\mu r}F^{r\tau}=\frac{iQ(uv)'}{(ruv\rho)^2}$, so that the Maxwell equations hold if and only if there is a constant $k\in\R^*$ such that $v=k/u$. Therefore, the stress-energy tensor $T_{\mu\nu}$ reads

\begin{align*}
T_{\mu\nu}&=\frac{1}{\mu_0}\left(g^{\alpha\beta}F_{\alpha\mu}F_{\beta\nu}-\frac14g_{\mu\nu}F_{\alpha\beta}F^{\alpha\beta}\right)=\frac{1}{4\pi}\left(\frac{Q^2}{2kr^4\rho^2}g_{\mu\nu}-\frac{1}{g_{\alpha\alpha}}F_{\mu\alpha}F_{\alpha\nu}\right) \\
&=\frac{Q^2}{8\pi r^2\rho}{\rm diag}\left(\frac{-u}{kr^2},\frac{-1}{ur^2},\frac{1}{k},\frac{\sin^2\theta}{k}\right).
\end{align*}

Next, the Ricci tensor $R_{\mu\nu}={\Gamma^\alpha}_{\mu\nu,\alpha}-{\Gamma^\alpha}_{\mu\alpha,\nu}+{\Gamma^\alpha}_{\beta\alpha}{\Gamma^\beta}_{\mu\nu}-{\Gamma^\alpha}_{\nu\beta}{\Gamma^\beta}_{\mu\alpha}$ and Einstein tensor $G_{\mu\nu}=R_{\mu\nu}-\tfrac12Rg_{\mu\nu}$ are diagonal and have

\[\left\{\begin{array}{l}
2kr\rho R_{\tau\tau}=-u(ru\rho''+r\rho u''+2\rho'(ru'+u)+2\rho u'), \\[1em]
2ru\rho^2R_{rr}=-3ru\rho\rho''-r\rho^2u''+3ru(\rho')^2-2\rho\rho'(ru'+u)-2\rho^2u', \\[1em]
2k\rho R_{\theta\theta}=-r^2u\rho''-\rho'(r^ 2u'+4ru)+2\rho'(k-u-ru'),\end{array}\right.\]

as well as

\[\left\{\begin{array}{l}
kr^2\rho^2 G_{\tau\tau}=u(r^2u\rho\rho''-3r^2u(\rho')^2/4+r\rho\rho'(ru'+4u)/2+\rho^2(k-u-ru')),\\[1em]
k\rho^2 G_{\theta\theta}=r(ru\rho\rho''+r\rho^2u''/2-3ru(\rho')^2/4+\rho\rho'(ru'+u)+\rho^2u'),\end{array}\right.\]

and we have $R_{\phi\phi}=\sin^2\theta R_{\theta\theta}$ and $G_{\phi\phi}=\sin^2\theta G_{\theta\theta}$. Thus, if the Einstein equation $G_{\mu\nu}=8\pi T_{\mu\nu}$ holds, then we have

\[0=\frac{2\rho^2}{u^2}\left(u^2(G_{rr}-8\pi T_{rr})-k(G_{\tau\tau}-8\pi T_{\tau\tau})\right)=3(\rho')^2-2\rho\rho'',\]

then the function $\rho$ satisfies $(\rho/\rho')'=-1/2$, so that there are constants $\alpha,\beta\in\R$ such that $\rho(r)=(\alpha r+\beta)^{-2}$. The metric now reads

\[{\rm d}s^2=(\alpha r+\beta)^{-2}(u(r){\rm d}\tau^2+ku(r)^{-1}{\rm d}r^2+r^2{\rm d}\Omega^2).\]

We then distinguish two cases:

If $\beta\ne0$, then the field equations reduce to the single equation $G_{rr}=8\pi T_{rr}$, which may be re-written as

\[\beta r^3(\alpha r+\beta)u'(r)+\beta r^2u(r)(\beta-2\alpha r)+(\alpha r+\beta)^2((\alpha^2Q^2-k)r^2+2Q^2\alpha\beta r+Q^2\beta^2)=0,\]

an equation which is equivalent to

\[\left(\frac{ru(r)}{(\alpha r+\beta)^3}\right)'=\frac{(k-\alpha^2Q^2)r^2-2Q^2\alpha\beta r-Q^2\beta^2}{\beta r^2(\alpha r+\beta)^2}\]

and we integrate the right-hand side by decomposing it into simple rational fractions. For $\alpha\ne0$, we obtain the solution

\[u(r)=\frac{(\alpha r+\beta)^3}{r}\left({\gamma}+\frac{1}{\beta}\left(\frac{Q^2}{r}-\frac{k}{\alpha(\alpha r+\beta)}\right)\right),\]

for some additional constant $\gamma\in\R$. If $\alpha=0$, then we get the equation

\[r^3u'(r)+r^2(u(r)-k)+Q^2\beta^2=0,\]

whose solution is

\[u(r)=k+\tilde{\alpha}/r+Q^2\beta^2/r^2\]

for some $\tilde{\alpha}\in\R$. Observe now that since $\beta\ne0$, the map $r\mapsto r(\alpha r+\beta)^{-1}$ is a diffeomorphism onto its image for any $\alpha$ and if we let

\[\tilde{\tau}:=\frac{\tau\sqrt{k}}{\beta},~R:=\frac{r}{\alpha r+\beta},~\tilde{Q}:=\frac{Q}{\sqrt{k}},~\tilde{M}:=\left\{\begin{array}{ll} \frac{1}{2}\left(\frac1\alpha+\frac{\alpha Q^2-\beta^2\gamma}{k}\right) & \text{if $\alpha\ne0$},\\[.75em] -\frac{\tilde{\alpha}}{2k\beta} & \text{otherwise,}\end{array}\right.\]

then we obtain the Reissner–Norström form of the statement.

Assume now that $\beta=0$ (forcing $\alpha\ne0$), in which case the equation $G_{rr}=8\pi T_{rr}$ reduces to $k=\alpha^2Q^2$. Moreover, the remaining non-trivial component of the field equation is $G_{\theta\theta}=8\pi T_{\theta\theta}$, which yields

\[r^2u''(r)-2ru'(r)+2(u(r)-\alpha^2 Q^2)=0,\]

whose solution reads $u(r)=\check{\beta}r^2+\check{\gamma}r+\alpha^2 Q^2$ for constants $\tilde{\beta},\tilde{\gamma}\in\R$. Letting

\[\check{\tau}:=\alpha Q\tau,~q:=\frac{\sqrt{\check{\beta}}}{\alpha Q},~m:=-\frac{\check{\gamma}}{2\alpha^2Q^2}\]

we obtain the form

\[{\rm d}s^2=\frac{1}{\alpha^2}\left[\frac{1-2mr+q^2r^2}{r^2}{\rm d}\tilde{\tau}^2+\frac{{\rm d}r^2}{r^2(1-2mr+q^2r^2)}+{\rm d}\Omega^2\right].\]

But in the new coordinates $(\tilde{\tau},r,\theta,\phi)$, the vector potential becomes $A_\mu=-i\alpha^{-1}r^{-1}{\rm d}\tilde{\tau}$, so that we may assume that $\alpha=1/Q$, thus obtaining the stated form of the metric.

To conclude, it remains to compute the Kretschmann scalar $K=R^{\alpha\beta\mu\nu}R_{\alpha\beta\mu\nu}$ of each metric of the statement. To do this, we use the symmetry and give the non-zero components $R_{\alpha,\beta,\mu,\nu}=g_{\alpha\lambda}({\Gamma^\lambda}_{\beta\nu,\mu}-{\Gamma^\lambda}_{\beta\mu,\nu}+{\Gamma^\lambda}_{\sigma\mu}{\Gamma^\sigma}_{\beta\nu}-{\Gamma^\lambda}_{\sigma\nu}{\Gamma^\sigma}_{\beta\mu})$ of the Riemann tensor for indices $\alpha<\beta$ and $\mu<\nu$. First, for the Bertotti–Robinson metric we have

\[R_{\tau r\tau r}=-\frac{Q^2}{r^4},~R_{\theta\phi\theta\phi}=Q^2\sin^2\theta~\Longrightarrow~K=\frac{8}{Q^4}\]

is indeed independent of $r$. For the Reissner–Nordström form, we have

\[\left\{\begin{array}{ll}
R_{\tau\theta\tau\theta}=\frac{R_{\tau\phi\tau\phi}}{\sin^2\theta}=\frac{(\tilde{Q}^2-\tilde{M}R)(R^2-2\tilde{M}R+\tilde{Q}^2)}{R^4}, & R_{\tau r\tau r}=\frac{2\tilde{M}R-3\tilde{Q}^2}{R^4},\\[1em]
R_{r\theta r\theta}=\frac{R_{r\phi r_\phi}}{\sin^2\theta}=\frac{\tilde{Q}^2-\tilde{M}R}{R^2-2\tilde{M}R+\tilde{Q}^2}, & R_{\theta\phi\theta\phi}=(2\tilde{M}R-\tilde{Q}^2)\sin^2\theta.\end{array}\right.\]

and thus

\[K=\frac{8(6\tilde{M}^2R^2-12\tilde{M}\tilde{Q}^2R+7\tilde{Q}^4)}{R^8},\]

which depends on $R$ and thus the Kretschmann scalar expressed in the original coordinates depends on $r$ as well, as claimed.

We make the following observations:

* The Ricci scalar of any of the above solutions vanishes.

* The proof shows in particular that for a vector potential $A_\tau=\mathcal{Q}r^{-1}$, a spherically symmetric solution of the field equation is Euclidean (resp. Lorentzian) if and only if $\mathcal{Q}$ is purely imaginary (resp. is real).

* Recalling the notation from the proof, we observe that in the Reissner–Nordström case, the new potential is

\[A'_\mu=-i\tilde{Q}R^{-1}{\rm d}\tilde{\tau}=-iQ(r^{-1}+\alpha\beta^{-1}){\rm d}\tau=A_\mu-i\nabla_\mu f,\]

where $f:=Q\alpha\beta^{-1}\tau$. Therefore, the coordinate transformation $(\tau,r)\mapsto(\tilde{\tau},R)$ induces a gauge transformation $A_\mu\mapsto A_\mu-i\nabla_\mu f$.

* The second metric of the statement with $m=q=0$ gives the Euclidean version of the original Bertotti–Robinson line element derived in <cit.>

\[{\rm d}s^2=\frac{Q^2}{r^2}\left[{\rm d}\tau^2+{\rm d}r^2+r^2{\rm d}\Omega^2\right].\]

Observe moreover that in Binet variable $u=1/r$, the general Bertotti–Robinson metric has an even simpler form

\[{\rm d}s^2=Q^2\left[(u^2-2mu+q^2){\rm d}\tau^2+\frac{{\rm d}u^2}{u^2-2mu+q^2}+{\rm d}\Omega^2\right].\]

The Euclidean Reissner–Nordström metric is the only spherically symmetric, asymptotically (Euclidean) flat metric satisfying the electro-vacuum Einstein–Maxwell equation associated to the complex vector potential

\[A_\mu:=-iQr^{-1}{\rm d}\tau.\]

More precisely, if ${\rm d}s^2=g_{\mu\nu}{\rm d}x^\mu{\rm d}x^\nu$ is a such a metric, defined for $r\gg0$, then there exists a constant $\tilde{M}\in\R$ such that

\[{\rm d}s^2=\left(1-\frac{2\tilde{M}}{r}+\frac{Q^2}{r^2}\right){\rm d}\tau^2+\left(1-\frac{2\tilde{M}}{r}+\frac{Q^2}{r^2}\right)^{-1}{\rm d}r^2+r^2{\rm d}\theta^2+r^2\sin^2\theta{\rm d}\phi^2.\]

In particular, the metric (<ref>) is the only spherically symmetric solution of the Einstein–Maxwell equation associated to $A_\mu$ defined on $\mathcal{M}$, which reduces to the Euclidean Schwarzschild metric (<ref>) when $Q\to0$.

Such a solution ${\rm d}s^2=g_{\mu\nu}{\rm d}x^\mu{\rm d}x^\nu$ is diagonal and since it is assumed to be asymptotically flat, the Kretschmann scalar should vanish as $r\to+\infty$, thus only the Reissner–Nordström form from the previous theorem is allowed, with $\tilde{\tau}=C\tau$ and $R=r/(Ar+B)$. But the asymptotic conditions

\[\lim_{r\to+\infty}(Ar+B)^{-2}=\lim_{r\to +\infty}\frac{g_{\theta\theta}}{r^2}=1=\lim_{r\to+\infty}g_{\tau\tau}=\lim_{r\to+\infty}C^2(1-2\tilde{M}/R+\tilde{Q}^2/R^2)\]

impose $R=r$ and $\tilde{\tau}=\tau$. Finally, the electromagnetic tensor $F_{\mu\nu}$ has $iF_{\tau r}=\tilde{Q}R^{-2}=Qr^{-2}$, so that $\tilde{Q}=Q$, hence the result.

§ GEODESIC MOTION IN EUCLIDEAN REISSNER–NORSTRÖM INSTANTONS

§.§ Motion equations and energy of orbits

Recall the notation from the beginning of the previous section and consider a non-constant geodesic $\gamma=(\tau,r,\theta,\phi)$ in $\mathcal{M}=\R\times ]r_+,+\infty[\times\Sph^2$ for the metric (<ref>), with affine parameter $\lambda$. We will analyse the geodesic equation in the same fashion as in <cit.>.

By spherical symmetry, we may assume that $\theta\equiv\pi/2$ and letting

\[\Delta(r):=1-\frac{2M}{r}+\frac{Q^2}{r^2},\]

the relativistic Lagrangian $\mathcal{L}=\tfrac12 g_{\mu\nu}\dot{\gamma}^\mu\dot{\gamma}^\nu$ reads

\[2\mathcal{L}=\Delta(r)\dot{\tau}^2+\Delta(r)^{-1}\dot{r}^2+r^2\dot{\phi}^2.\]

Thus, the temporal and angular Euler–Lagrange equations provide constants $C,J\in\R$ such that

\begin{equation}\label{tandphi}
\dot{\tau}=\frac{C}{\Delta(r)},~\dot{\phi}=\frac{J}{r^2}
\end{equation}

and thus the scalar

\[\mathcal{H}:=2\mathcal{L}=\Delta(r)^{-1}(C^2+\dot{r}^2)+\frac{J^2}{r^2}>0\]

is conserved along $\gamma$ and the proper length $s$ satisfies ${\rm d}s^2=\mathcal{H}{\rm d}\lambda^2$ so that we get

\begin{equation}\label{mot_ERN}
\left(\frac{{\rm d}\phi}{{\rm d}s}\right)^2=\frac{\dot{\phi}^2}{\mathcal{H}}=\frac{L^2}{r^4},~\left(\frac{{\rm d}r}{{\rm d}s}\right)^2=\frac{\dot{r}^2}{\mathcal{H}}=\Delta(r)\left(1-\frac{L^2}{r^2}\right)-E^2.
\end{equation}

where $E$ (resp. $L$) is the energy per unit mass (resp. angular momentum per unit mass) of $\gamma$, defined by analogy with the Lorentzian framework as the proper temporal (resp. azimuthal) conjugate momentum

\[E:=p_{\tau(s)}=\tfrac{1}{\sqrt{|\mathcal{H}|}}p_{\tau(\lambda)}=\frac{g_{\tau\mu}\dot{\gamma}^\mu}{\sqrt{\mathcal{H}}}=\frac{\Delta(r)\dot{\tau}}{\sqrt{\mathcal{H}}}=\frac{C}{\sqrt{\mathcal{H}}},~\left(\text{resp. } L:=\tfrac{1}{\sqrt{|\mathcal{H}|}}p_\phi=\frac{J}{\sqrt{\mathcal{H}}}\right)\]

For the rest of this section, we assume that $\gamma$ is non-purely radial and we denote the differentiation with respect to the Euclidean time $\tau$ with a dot.

We want to apply the Weierstrass analysis of this equation and since it has degree 4, one first needs to choose a real root of the quartic right-hand side. Such a real root is guaranteed to exist as soon as $E^2<1$, a property that we shall prove to always hold, provided that the metric (<ref>) presents an event horizon (that is, when $Q^2<M^2$). To do this, we need expressions for the energy and angular momentum, as functions of the initial conditions $\gamma(0)=:(\tau_0,r_0,\pi/2,\phi_0)$ and $\dot{\gamma}(0)=:(1,\dot{r}_0,0,\dot{\phi}_0)$. By symmetry, we may assume that $\tau_0=\phi_0=0$. We have the following result:

The energy and angular momentum of a non-constant equatorial geodesic $\gamma$ parametrized by coordinate time, with $\gamma(0)=(r_0,0)$ and $\dot{\gamma}(0)=(\dot{r}_0,\dot{\phi}_0)$ are given by

\[E=\frac{\Delta(r_0)}{\sqrt{\Delta(r_0)+\Delta(r_0)^{-1}\dot{r}_0^2+r_0^2\dot{\phi}_0^2}},~L=\frac{r_0^2\dot{\phi}_0}{\sqrt{\Delta(r_0)+\Delta(r_0)^{-1}\dot{r}_0^2+r_0^2\dot{\phi}_0^2}}.\]

Let $\alpha:=\left.\frac{{\rm d}\tau}{{\rm d}\lambda}\right|_{\lambda=0}$. The constant $C$ reads $C=\alpha\Delta(r_0)$ and we also have $J=\alpha r_0^2\dot{\phi}_0$. Then,

\[\mathcal{H}=\Delta(r_0)^{-1}\left(C^2+\left(\frac{{\rm d}r}{{\rm d}\lambda}\right)^2\right)+\frac{J^2}{r_0^2}=\alpha^2\left(\Delta(r_0)+\Delta(r_0)^{-1}\dot{r}_0^2+r_0^2\dot{\phi}_0^2\right),\]

hence the result.

If $Q^2\le M^2$, then any (exterior) non-constant geodesic $\gamma$ for the metric (<ref>) has $E^2<1$. Otherwise, there are circular orbits with arbitrary energy.

Assume first that $Q^2\le M^2$. For simplicity, we may rescale the radius by $1/M$ and assume that $M=1$. From the Lemma <ref>, we have

\[E^2\ge1~\Longleftrightarrow~\Delta(r_0)^2\ge\Delta(r_0)+\Delta(r_0)^{-1}\dot{r}_0^2+r_0^2\dot{\phi}_0^2~\Longleftrightarrow~P(r_0)\le0,\]

where $P=P_{Q^2,\dot{r}_0^2,\dot{\phi}_0^2}\in\R_8[x]$ is the following polynomial of degree at most 8:

\begin{align*}
\end{align*}

We will inductively prove that $P$ is positive on the interval $I=]r_+,+\infty[$ (with $r_+=1+\sqrt{1-Q^2}$), hence contradicting the above inequality and establishing the result. Observe first that $P^{(8)}/8!=\dot{\phi}_0^2\ge0$ so $P^{(7)}$ is non-decreasing on $I$ and since

\[\tfrac{1}{7!}P^{(7)}(r_+)=2\dot{\phi}_0^2(3+4\sqrt{1-Q^2})\ge0\]

we find that $P^{(7)}\ge0$ so that $P^{(6)}$ is non-decreasing on $I$. We repeat the process, with

\[\tfrac{1}{6!}P^{(6)}(r_+)=3\dot{\phi}_0^2\left(5+9(1-Q^2)+14\sqrt{1-Q^2}\right)+\dot{r}_0^2,\]

which is positive since $\gamma$ is non-constant. Thus, we get $P^{(6)}>0$ on $I$ and $P^{(5)}$ is increasing (not only non-decreasing) on $I$. Next, we have

\begin{align*}
\tfrac{1}{5!}P^{(5)}(r_+)=&2\sqrt{1-Q^2}\left[45\dot{\phi}_0^2+25\dot{\phi}_0^2(1-Q^2)+3\dot{r}_0^2\right] \\
&+20\dot{\phi}_0^2+120\dot{\phi}_0^2(1-Q^2)+2(1+3\dot{r}_0^2)>0, \\[1em]
\tfrac{1}{4!}P^{(4)}(r_+)=&5(1+\sqrt{1-Q^2})^2\left[3(Q^2\dot{\phi}_0^2+\dot{r}_0^2)+14\dot{\phi}_0^2(1+\sqrt{1-Q^2})\sqrt{1-Q^2}\right] \\
&+10\sqrt{1-Q^2}+2-Q^2>0, \\[1em]
\tfrac{1}{3!}P^{(3)}(r_+)=&\sqrt{1-Q^2}\left[12\dot{\phi}_0^2(3Q^4+2)+60\dot{r}_0^2+4+4(1-Q^2)(78\dot{\phi}_0^2+5\dot{r}_0^2+1)\right] \\
\end{align*}

\begin{align*}
\tfrac12 P''(r_+)=&\sqrt{1-Q^2}\left[60\dot{r}_0^2+4(1-Q^2)(2+15\dot{r}_0^2)+2\dot{\phi}_0^2(33Q^4-136Q^2+112)\right]+15\dot{r}_0^2Q^4 \\
&+\dot{\phi}_0^2(-13Q^6+174Q^4-384Q^2+224)+4(1-Q^2)(30\dot{r}_0^2+(1-Q^2)+1)\ge0, \\[1em]
\end{align*}

and thus, $P'$ is increasing on $I$ and has $P'(r_+)\ge0$, so $P$ is increasing on $I$. We conclude by observing that $P(r_+)=\dot{r}_0^2r_+^6\ge0$.

To prove the second statement, assume that $Q^2>M^2$, choose $e>0$ and consider a circular (equatorial) geodesic $\gamma=(\tau,r,\pi/2,\phi)$ with energy $E=e$. Introducing the potential

\[V(r):=\sqrt{\Delta(r)\left(\frac{L^2}{r^2}-1\right)},\]

the condition that $\gamma$ is circular reads

\[V(r)^2+e^2=0=\frac{\partial V}{\partial r}.\]

The first equation imposes

\[L=\pm r\sqrt{1-\frac{e^2}{\Delta(r)}},\]

while the second yields

\begin{align*}
0=\frac{\partial V}{\partial r}&=\frac{1}{2V}\left(\left(\frac{L^2}{r^2}-1\right)\frac{\partial\Delta}{\partial r}-\frac{2L^2\Delta}{r^3}\right) \\
\end{align*}

This quartic has a positive root, since it evaluates to $-Q^4<0$ at $r=0$ and is positive when $r\gg0$.

When $Q=0$, we recover the result from <cit.>.

The metric (<ref>) has an event horizon if and only if no orbit is bounded.

§.§ Reduction of the polar radial equation to Weierstrass' form

Let $\gamma=(\tau,r,\phi)$ be a non-purely radial equatorial geodesic. Then $J\ne0$ so that the map $s\mapsto\phi(s)$ is a diffeomorphism onto its image and from (<ref>) we find the polar radial equation

\begin{equation}\label{polmot_ERN}
\left(\frac{{\rm d}r}{{\rm d}\phi}\right)^2=\frac{1-E^2}{L^2}r^4-\frac{2M}{L^2}r^3+\left(\frac{Q^2}{L^2}-1\right)r^2+2Mr-Q^2=:F(r)
\end{equation}

In this section, the dot denotes differentiation with respect to the polar variable $\phi$. We use the same trick as in <cit.> to reduce the degree of the above equation and then re-write it in Weierstrass form. Let $\overline{r}\in\C$ be a root of the quartic $F$ (which, in view of <ref>, is guaranteed to be real positive when $Q^2\le M^2$) and consider the shifted Binet variable

\[u:=\frac{1}{r-\overline{r}}.\]

Then, the equation (<ref>) becomes

\begin{align*}
\dot{u}^2=&\frac{\dot{r}^2}{(r-\overline{r})^4}=u^4F(\overline{r}+1/u) \\
\end{align*}

It is now straightforward to re-write this in Weierstrass form. We summarize the discussion in the following result:

Let $\gamma=(\tau,r,\phi)$ be a non-purely radial equatorial geodesic in the metric (<ref>), with energy $E$ and angular momentum $L$. The map $\lambda\mapsto\phi(\lambda)$ being a diffeomorphism onto its image, we re-parametrize $\gamma$ using the angular variable $\phi$ and we abusively still denote by $r$ the re-parametrized coordinate $\phi\mapsto r(\phi)$. Let $\overline{r}\in\C$ be a root of the quartic

\[\frac{1-E^2}{L^2}r^4-\frac{2M}{L^2}r^3+\left(\frac{Q^2}{L^2}-1\right)r^2+2Mr-Q^2=0,\]

which we can choose to be real positive if $Q^2\le M^2$, and if we let

\[\left\{\begin{array}{l}
\delta=\frac{1-E^2}{L^2},\\[.5em]
\gamma=2\left(2\delta\overline{r}-\frac{M}{L^2}\right),\\[.5em]
\beta=6\overline{r}\left(\delta\overline{r}-\frac{M}{L^2}\right)+\frac{Q^2}{L^2}-1,\\[.5em]
\alpha=2\left(2\delta\overline{r}^3-\frac{3M}{L^2}\overline{r}^2+\left(\frac{Q^2}{L^2}-1\right)\overline{r}+M\right),\end{array}\right.~~\text{as well as}~~\left\{\begin{array}{l} g_2:=\frac14\left(\frac{\beta^2}{3}-\alpha\gamma\right),\\[.5em]
\wp:=\frac{\alpha}{4(r-\overline{r})}+\frac{\beta}{12},\end{array}\right.\]

then the function $\wp$ satisfies the Weierstrass equation

\[\dot{\wp}^2=4\wp^3-g_2\wp-g_3.\]

In other words, the polar radial motion is given by

\[r(\phi)=\overline{r}+\frac{\alpha}{4\wp(\phi)-\beta/3},\]

where $\wp=\wp_{g_2,g_3}$ is the Weierstrass function associated to the pair $(g_2,g_3)\in\C^2$.

More practically, given an initial radius $r_0:=r(\phi_0)$, we have to find some $z_0\in\C$ such that $\wp(z_0)=\frac{\alpha}{4(r_0-\overline{r})}+\frac{\beta}{12}$, a task that can be achieved using Carlson's integrals <cit.>

\begin{equation}\label{carlsonRF}\tag{$R_F$}
R_F(x,y,z):=\frac12\int_0^\infty\frac{{\rm d}\zeta}{\sqrt{(\zeta+x)(\zeta+y)(\zeta+z)}}.
\end{equation}

These integrals can be approximated efficiently by the Carlson algorithm <cit.> and we can approach $\wp$ using the Coquereaux–Grossmann–Lautrup algorithm <cit.>. This is the method we use to approximate ERN orbits and produce a Python code. We summarize the method in the following result:

Fix $(r_0,L,E)\in]r_+,+\infty[\times\R^*\times\R^*_+$ and let $\gamma$ be the unique maximal non-purely radial equatorial ERN geodesic with energy $E$, angular momentum $L$ and such that $r(0)=r_0$ in polar parametrization $r=r(\phi)$. Recalling the constants $\overline{r},\alpha,\beta,\gamma,\delta,g_2,g_3$ from the Proposition <ref>, if $g_2^3-27g_3^2\ne0$, then the function $r$ is given by

\[r(\phi)=\overline{r}+\frac{\alpha}{4\wp_{g_2,g_3}(z_0+\phi)-\beta/3},\text{ with }z_0:=R_F(\wp_0-z_1,\wp_0-z_2,\wp_0-z_3)\in\C,\]

where $z_{1,2,3}\in\C$ are the roots of the Weierstrass cubic $4z^3-g_2z-g_3$ and $\wp_0:=\frac{\alpha}{4(r_0-\overline{r})}+\frac{\beta}{12}$.

§.§ The special case of the Euclidean Schwarzschild geometry

When $Q=0$, we have $F(0)=0$ so that we may take $\overline{r}=0$ so that the polar equation in Binet variable simplifies to

\[\dot{u}^2=u^4F(1/u)=2Mu^3-u^2-\frac{2M}{L^2}u+\frac{1-E^2}{L^2}\]

which is <cit.>. Therefore, the Weierstrass form

\[\dot{\wp}^2=4\wp^3-g_2\wp-g_3\]

is obtained by letting $\wp:=M/(2r)-1/12$, as well as

\begin{equation}\label{ges}\tag{$g_{12}^{\rm ES}$}
g_2:=\frac{1}{12}+\frac{M^2}{L^2}\text{ and }g_3:=\frac{1}{216}-\frac{M^2}{12L^2}(2-3E^2).
\end{equation}

These expression can be compared to the Lorentzian Schwarzschild case, where for a test-particle of mass $\mu=2\mathcal{L}$ (twice the Schwarzschild Lagrangian), the constants $g_2$ and $g_3$ read (cf <cit.>)

\begin{equation}\label{gls}\tag{$g_{12}^{\rm LS}$}
g_2=\frac{1}{12}+\frac{\mu M^2}{L^2}\text{ and }g_3=\frac{1}{216}-\frac{M^2}{12L^2}(2\mu+3E^2).
\end{equation}

Therefore, one may view a Euclidean Schwarzschild geodesic as a space-like Lorentzian Schwarzschild geodesic with complex energy $E_{{\rm Euc}}=iE_{{\rm Lor}}$. Notice that this last equality is expected since the energy is defined as the temporal momentum and because the Euclidean time $\tau$ and the Lorentzian time $t$ are related by the relation $\tau=it$.

The Weierstrass formulation also permits to derive a shorter proof of the fact that $E<1$ for every geodesic $\gamma$ in Euclidean Schwarzschild geometry, with initial radius $r_0>r_+=2M$. As above, we may simplify the notation by rescaling the radius and assuming that $M=1$. If $L=0$, then the second equation (<ref>) reduces to $({{\rm d}r}/{{\rm d}s})^2=1-E^2-{2}/{r}$ so $1-E^2\ge 2/r>0$, as claimed. Now if $L\ne0$, then we may use Weierstrass' form and the discriminant of the equation $4z^3-g_2z-g_3$ reads

\[\Delta:=16(g_2^3-27g_3^2)=L^{-6}\left[(1-E^2)L^4-(27E^4-36E^2+8)L^2+16\right]=\frac{1-E^2}{L^2}+\bigo\left(\frac{1}{L^4}\right),\]

so that if, for the sake of contradiction, we assume $E^2>1$, then $\Delta<0$ for $L\gg0$, something which can be achieved by rescaling the initial azimuthal angular velocity. Therefore, we assume that $\Delta<0$ and look for an absurdity.

Consider the Weierstrass cubic $\mathfrak{q}(x):=4x^3-g_2x-g_3$; then the phase portrait in Weierstrass variable $x=1/(2r)-1/12$ describes a portion of the (real) elliptic curve

\[\mathcal{E}_{\mathfrak{q}}:=\{(x,y)\in\R^2~|~y^2=\mathfrak{q}(x)\}.\]

Since $\lim_{x\to+\infty}\mathfrak{q}(x)=+\infty$, we have $\mathcal{E}_\mathfrak{q}\cap\{x=x_0\}\ne\emptyset$ for $x_0\gg0$ but because $\mathfrak{q}(1/6)=-(E/2L)^2<0$, the non-compact connected component of $\mathcal{E}_\mathfrak{q}$ lies in the open half-plane $\{x>1/6\}$ (which corresponds in radial variable to $\{r<2\}$), so that the considered phase portrait cannot describe a portion of this component. However, to say that $\Delta<0$ amounts to say that $\mathcal{E}_\mathfrak{q}$ is connected, a contradiction.

At this point, we know that $\Delta\ge0$. If $\Delta>0$, then the elliptic curve $\mathcal{E}_\mathfrak{q}$ has an additional compact connected component and since $\mathfrak{q}(-1/12)=(1-E^2)/(4L^2)>0$, this component intersects the subset $\{x=-1/12\}$. Therefore, the corresponding orbit is indeed unbounded. In the case where $\Delta=0$, the curve is singular but connected so that the orbit imposes $0\le\mathfrak{q}(1/6)=-(E/2L)^2\le0$ so $r=2$, a new contradiction.

We summarize the discussion in the following result:

Any Euclidean Schwarzschild orbit has $E^2<1$, and is unbounded.

Moreover, the phase portrait of a non-radial equatorial Euclidean Schwarzschild orbit in Weierstrass variable $\wp=M/(2r)-1/12$ (with polar argument) describes a portion of the unique compact connected component of the associated (real) elliptic curve $\mathcal{E}$, included in $\mathcal{E}\cap\{-1<12\wp<2\}$. In particular, the discriminant of $\mathcal{E}$ is positive.

§ GRAVITATIONAL DEFLECTION OF EUCLIDEAN SCHWARZSCHILD ORBITS

Consider an equatorial orbit $\gamma=(\tau,r,\phi)$, with energy $-1<E<1$ and angular momentum $L\in\R^*$. Since $\gamma$ is unbounded, we may consider its velocity at infinity, defined by $v^2:=\lim_{r\to\infty}({\rm d}r/{\rm d}\tau)^2$. Using equations (<ref>) and (<ref>), we find the expression

\[v^2=\lim_{r\to\infty}\left(\frac{{\rm d}r}{{\rm d}\tau}\right)^2=\lim_{r\to\infty}\left(\frac{\dot{r}}{\dot{\tau}}\right)^2=\lim_{r\to\infty}\frac{\Delta(r)^2}{C^2}\left(\Delta(r)\left(\mathcal{H}-\frac{J^2}{r^2}\right)-C^2\right)=\frac{1}{E^2}-1\]

The equation (<ref>) can be recast in Binet variable $u=1/r$ and yields

\begin{equation}\label{polES}
\left(\frac{{\rm d}u}{{\rm d}\phi}\right)^2+u^2=2Mu^3-\frac{2M}{L^2}u+\frac{1}{b^2}=2Mu\left(u^2-\frac{1}{L^2}\right)+\frac{1}{b^2},
\end{equation}

where the constant $b:=\sqrt{\frac{L^2}{1-E^2}}$ is the impact parameter of $\gamma$, satisfying $L=\pm bvE$. Since

\[\lim_{r\to\infty}\frac{{\rm d}\phi}{{\rm d}\tau}=\lim_{r\to\infty}\frac{J\Delta(r)}{Cr^2}=0,\]

the orbit admits asymptotic lines, and we are first interested in the deflection angle at infinity $\delta\phi$ between these two asymptotic directions, as a function of the orbit's perihelion $r=r_{\rm min}$.

(z) at (0,0);
(zp) at (0,4);
(zm) at (0,-4);

[color=red,domain=-1.55*pi/4:1.55*pi/4,scale=1,samples=200] plot (deg():2/sqrt(3^2*cos(r)^2-1));

(pu) at (7*1/3,7*sqrt(1-1/3^2));
(pup) at (7.5*1/3,7.5*sqrt(1-1/3^2));
(pd) at (7.5*1/3,-7.5*sqrt(1-1/3^2));
(mpu) at (-2*1/3,-2*sqrt(1-1/3^2));
(mpd) at (-2*1/3,2*sqrt(1-1/3^2));
(a) at (2/sqrt(3^2-1),0);
(ap) at (2/sqrt(3^2-1),5);
(am) at (2/sqrt(3^2-1),-5);
(f1) at (2*3/sqrt(3^2-1),0);

(Pb) at (2*3/sqrt(3^2-1)+1.15*(5.5*1/3),0+1.15*(5.5*sqrt(1-1/3^2)));
(Pbp) at (2*3/sqrt(3^2-1)+1.25*(5.5*1/3),0+1.25*(5.5*sqrt(1-1/3^2)));

[fill=black] (f1) circle (6pt);

[dotted] (mpu)–(pup) (mpd)–(pd);
(z)–(a) (ap)–(am);
[dotted] (f1)–(Pbp);
[ultra thick] (a)–(f1);
[<->,dashed,very thick] (Pb)–(pu);

($(Pb)!0.5!(pu)$) node[right]$b$;
($(a)!1.1!(f1)$) node[below right]$0$;
($(f1)!0.5!(a)$) node[right]$r_{\rm min}$;

($(z)!-0.75!(f1)$) coordinate (A) (z) coordinate (B) (pd) coordinate (C) pic ["$\phi$",draw,->,angle radius=0.8cm]angle;
(pu) coordinate (A) (z) coordinate (B) (mpd) coordinate (C) pic ["$\delta\phi$",draw,->,angle radius=1.5cm]angle;

(de) at (0.71,2);
[fill=black] (de) circle (1.5pt);
[fill=black] (z) circle (1.5pt);
(pu) coordinate (A)  (de) coordinate (B)  (ap) coordinate (C) pic ["$\delta\phi/2$",draw,->,angle radius=2.5cm]angle;

Schematics of an orbit with angle deflection $\delta\phi=2\phi-\pi$.

§.§ Analytic expression of the deflection angle using Carlson's integrals

As illustrated in the Figure <ref>, the deflection (at infinity) $\delta\phi$ is given by

\[\delta\phi=2|\phi(r=\infty)-\phi(r=r_{\rm min})|-\pi\]

and using (<ref>) again leads to the expression

\[\delta\phi=2\int_{r_{\rm min}}^\infty{\rm d}\phi-\pi=2\int_{r_{\rm min}}^\infty\frac{{\rm d}\phi}{{\rm d}r}{\rm d}r-\pi=2\int_{r_{\rm min}}^\infty\frac{1}{\sqrt{\frac{1}{b^2}-\frac{1}{r^2}\left(1-\frac{2M}{r}+\frac{2Mr}{L^2}\right)}}\frac{{\rm d}r}{r^2}-\pi.\]

This expression can be simplified using the Weierstrass variable $\wp=Mu/2-1/12$. Indeed, using the constants $g_2,g_3$ given by equation (<ref>), we have

\[\delta\phi=2\left|\int_0^{1/r_{\rm min}}\frac{{\rm d}u}{\sqrt{\frac{1}{b^2}-u^2+2Mu\left(u^2-\frac{1}{L^2}\right)}}\right|-\pi=2\int_{-\tfrac{1}{12}+\tfrac{M}{2r_{\rm min}}}^ {-\tfrac{1}{12}}\frac{{\rm d}p}{\sqrt{4p^3-g_2p-g_3}}-\pi.\]

This expression can be rewritten in terms of elliptic integrals, as in the pioneer work <cit.>. However, it is both easier and numerically more adequate to express it with the integrals (<ref>), which we numerically approximate using Carlson's algorithm <cit.>. Observe first that because $r=r_{\rm min}$ is a turning point of $\gamma$, we have

\[0=\left.\frac{{\rm d}\wp^2}{{\rm d}\phi^2}\right|_{r=r_{\rm min}}=4\wp(r=r_{\rm min})^3-g_2\wp(r=r_{\rm min})-g_3.\]

In other words, the point $\wp_{\rm max}:=M/(2r_{\rm min})-1/12$ is a root of the Weierstrass cubic and this leads to the factorization

\[4p^3-g_2p-g_3=(p-\wp_{\rm max})\left(4p^2+\left(\frac{2M}{r_{\rm min}}-\frac13\right)p+\left(\frac{1}{36}-g_2-\frac{M}{3r_{\rm min}}+\frac{M^2}{r_{\rm min}^2}\right)\right),\]

which allows to find the other two roots $\wp_{\pm}\in\C$ of the cubic. We then have

\begin{align*}
\int^{-\tfrac{1}{12}}_{\wp_{\rm max}}\frac{{\rm d}p}{\sqrt{4p^3-g_2p-g_3}}&=\int^0_{\tfrac{M}{2r_{\rm min}}}\frac{{\rm d}\zeta}{\sqrt{4\left(\zeta-\frac{1}{12}\right)^3-g_2\left(\zeta-\frac{1}{12}\right)-g_3}} \\
&=\frac12\int^0_{\tfrac{M}{2r_{\rm min}}}\frac{{\rm d}\zeta}{\sqrt{\left(\zeta-\wp_{\rm max}-\tfrac{1}{12}\right)\left(\zeta-\wp_--\tfrac{1}{12}\right)\left(\zeta-\wp_+-\tfrac{1}{12}\right)}} \\
&=\frac12\left[\int_{\tfrac{M}{2r_{\rm min}}}^\infty-\int_0^\infty\right]\frac{{\rm d}\zeta}{\sqrt{\left(\zeta-\wp_{\rm max}-\tfrac{1}{12}\right)\left(\zeta-\wp_--\tfrac{1}{12}\right)\left(\zeta-\wp_+-\tfrac{1}{12}\right)}} \\
&=R_F\left(e_{\rm max}+\tfrac{M}{2r_{\rm min}},e_-+\tfrac{M}{2r_{\rm min}},e_++\tfrac{M}{2r_{\rm min}}\right)-R_F(e_{\rm max},e_-,e_+),
\end{align*}

where $e_{{\rm max},\pm}:=-\wp_{{\rm max},\pm}-1/12$. Thus, we obtain the following expression for the deflection:

\begin{equation}\label{deltaphi_carlson}
\delta\phi=2\left(R_F\left(e_{\rm max}+\tfrac{M}{2r_{\rm min}},e_-+\tfrac{M}{2r_{\rm min}},e_++\tfrac{M}{2r_{\rm min}}\right)-R_F(e_{\rm max},e_-,e_+)\right)-\pi
\end{equation}

This is an efficient formula for numerical calculations (see below), but it is natural to ask for an approximation of $\delta\phi$ when $r_{\rm min}\to\infty$. This is the goal of the next subsection. Observe moreover that the two terms of (<ref>) involving the integrals (<ref>) can be complex, but their difference is real and non-negative.

§.§ Approximation of the deflection angle with perturbed solution

The aim of this section is to obtain an expansion of $\delta\phi$ in powers of the perihelion $r_{\rm min}$, up to order 3. This choice of order will become transparent later, when we study the null and maximal deflection rings.

To do so, we could use for instance the well-known perturbed solution method <cit.>. However, we will avoid complicated calculations with a simple observation. First, using (<ref>), the turning point condition at $r=r_{\rm min}$ gives the following relation

\[0=\left(1-\frac{2M}{r}\right)\left(1-\frac{L^2}{r^2}\right)-E^2,\]

or, in terms of $b$ and $v$,

\[r_{\rm min}^2\left(1-\frac{2M}{r_{\rm min}}\left(1+\frac{1}{v^2}\right)\right)=b^2\left(1-\frac{2M}{r_{\rm min}}\right),\]

that is,

\begin{equation}\label{b_rm}
b=r_{\rm min}\sqrt{\frac{1-\frac{2M}{r_{\rm min}}\left(1+\frac{1}{v^2}\right)}{1-\frac{2M}{r_{\rm min}}}}=r_{\rm min}\sqrt{1-\frac{2M}{v^2({r_{\rm min}}-2M)}}.
\end{equation}

Recall also that $L^2=b^2v^2(1+v^2)^{-1}$ and that the equation (<ref>) in terms of $b$ and $L$ reads

\[\left(\frac{{\rm d}r}{{\rm d}\phi}\right)^2=r^4\left(\frac{1}{b^2}-\frac{1}{r^2}\left(1-\frac{2M}{r}+\frac{2Mr}{L^2}\right)\right).\]

On the other hand, in Lorentzian Schwarzschild geometry, the same equation

\[\left(\frac{{\rm d}r}{{\rm d}\phi}\right)^2=r^4\left(\frac{1}{b_{\rm Sch}^2}-\frac{1}{r^2}\left(1-\frac{2M}{r}-\frac{2Mr}{L_{\rm Sch}^2}\right)\right)\]

differs from the above one by just a sign in the term in $L^{-2}$, and we have in this case

\[b_{\rm Sch}^2=r_{\rm min}\sqrt{1+\frac{2M}{v^2({r_{\rm min}}-2M)}},\]

as well as $L_{\rm Sch}^2=b^2v^2(1-v^2)^{-1}$. This means that the expression of the deflection angle in Euclidean Schwarzschild geometry with squared velocity $v^2$ is the same as the expression of the Lorentzian deflection angle with same perihelion and “squared velocity” $-v^2$. This is mathematically well-defined since the data only depend on the squared quantities $(b^2,v^2)$, and not on the pair $(b,v)$ itself, so no complex number is involved. Moreover, this interpretation is physically consistent with the relation $\tau=it$ between the Euclidean time $\tau$ and the Lorentzian time $t$ and the definition of the Euclidean velocity $v^2=\lim_{r\to\infty}({\rm d}r/{\rm d}\tau)^2$, while in the Lorentzian case, $v^2=\lim_{r\to\infty}({\rm d}r/{\rm d}t)^2$.

After <cit.> or <cit.> for instance, up to order 3, in Lorentzian Schwarzschild geometry we have

\begin{equation}\label{def_sch}
\delta\phi=\frac{2M}{b}\left(1+\frac{1}{v^2}\right)+\frac{3\pi M^2}{4b^2}\left(1+\frac{4}{v^2}\right)+\frac{2M^3}{3b^3}\left(5+\frac{45}{v^2}+\frac{15}{v^4}-\frac{1}{v^6}\right)+\bigo\left(\frac{M^4}{b^4}\right).
\end{equation}

Therefore, changing the sign of the terms in $v^{-2}$ yields the following estimate, in Euclidean Schwarzschild geometry,

\begin{equation}\label{estimdeltaphi}\tag{$\delta_3^b$}
\delta\phi=\frac{2M}{b}\left(1-\frac{1}{v^2}\right)+\frac{3\pi M^2}{4b^2}\left(1-\frac{4}{v^2}\right)+\frac{2M^3}{3b^3}\left(5-\frac{45}{v^2}+\frac{15}{v^4}+\frac{1}{v^6}\right)+\bigo\left(\frac{M^4}{b^4}\right).
\end{equation}

Using now the relation (<ref>), we arrive at the expression

\begin{align*}\label{estimdeltaphir}\tag{$\delta_3^r$}
\delta\phi=&\frac{2M}{r_{\rm min}}\left(1-\frac{1}{v^2}\right)+\frac{M^2}{r_{\rm min}^2}\left(\frac{3\pi}{4}\left(1-\frac{4}{v^2}\right)+\frac{2}{v^2}\left(1-\frac{1}{v^2}\right)\right) \\
&+\frac{M^3}{r_{\rm min}^3}\left(\frac{3\pi}{2v^2}\left(1-\frac{4}{v^2}\right)+\frac{10}{3}-\frac{26}{v^2}+\frac{9}{v^4}-\frac{7}{3v^6}\right)+\bigo\left(\frac{M^4}{r_{\rm min}^4}\right).
\end{align*}

Concerning the metric (<ref>) with non-zero charge $Q\ne0$, the perturbed solution method gives, up to order 2 and after elementary calculations we omit,

\[\delta\phi=\frac{2M}{b}\left(1-\frac{1}{v^2}\right)+\frac{3\pi M^2}{4b^2}\left(1-\frac{4}{v^2}\right)-\frac{\pi Q^2}{4b^2}\left(1-\frac{2}{v^2}\right)+\bigo(b^{-3}),\]

in agreement with <cit.>. Expressing the impact parameter $b$ in terms of the perihelion $r_{\rm min}$ as

\[b=r_{\rm min}\sqrt{1-\frac{2Mr_{\rm min}-Q^2}{v^2(r_{\rm min}^2-2Mr_{\rm min}+Q^2)}},\]

yields the expansion

\begin{equation}\label{def_with_charge}
\delta\phi=\frac{2M}{r_{\rm min}}\left(1-\frac{1}{v^2}\right)+\frac{M^2}{r_{\rm min}^2}\left(\frac{3\pi}{4}\left(1-\frac{4}{v^2}\right)+\frac{2}{v^2}\left(1-\frac{1}{v^2}\right)\right)-\frac{\pi Q^2}{4r_{\rm min}^2}\left(1-\frac{2}{v^2}\right)+\bigo(r_{\rm min}^{-3}).
\end{equation}

As a sanity check, observe that replacing $v^2$ by $-v^2$ in the previous expression and letting $v\to1$ leads to the light deflection formula of <cit.>

\[\delta\phi=\frac{4M}{r_{\rm min}}+\frac{M^2}{r_{\rm min}^2}\left(\frac{15\pi}{4}-4\right)-\frac{3\pi Q^2}{4r_{\rm min}^2}+\bigo(r_{\rm min}^{-3}).\]

§.§ Null and maximal deflection rings

We start by observing that, at lowest order, the deflection angle for a usual (Lorentzian) massive Schwarzschild orbit with velocity at infinity $v$ is given by

\[\delta\phi_{\rm Lor}\approx\frac{2M}{r_{\rm min}}\left(1+\frac{1}{v^2}\right)>0,\]

while in the Euclidean background, the estimation (<ref>) reads

\[\delta\phi_{\rm Euc}\approx \frac{2M}{r_{\rm min}}\left(1-\frac{1}{v^2}\right).\]

This suggests that for the Euclidean Schwarzschild solution, the deflection $\delta\phi$ may vanish for some values of $r_{\rm min}$. This can also be noticed from the motion equation itself. Indeed, recall from <cit.> that a polar curve, parametrized in Binet variable $u$ is concave with respect to its pole (hence, has a positive deflection) if and only if $u+\frac{{\rm d}^2u}{{\rm d}\phi^2}>0$. Now, this quantity in Schwarzschild geometry is given by

\[\ddot{u}+u=M(3u^2-\epsilon L^{-2}),\]

where $\epsilon$ is the signature of the metric. Hence, this is always positive in the Lorentzian background $\epsilon=-1$, while it may change of sign in the Euclidean world $\epsilon=1$.

To numerically appreciate this phenomenon, in the Figure <ref> we depict pencils of orbits with different velocities. Their deflection at infinity is computed using the formula (<ref>) and Carlson's algorithm. Observe that the deflection stays negative when $v<1$. The more extreme case where $v=5$ is displayed in the Figure <ref>, along with the graph representing the deflection angle as a function of the perihelion. Observe the presence of a null deflection at $r_{\rm min}\approx 2.16$ and a maximal deflection at $r_{\rm min}\approx2.5$. The Figure <ref> depicts the deflection as a function of the perihelion, for several velocities at infinity.

Orbits with different velocities and their perihelion (in units of $M$) and deflection at infinity.

Orbits with $2.1<r_{\rm min}<5$.

The deflection $\delta\phi$ as a function of $r_{\rm min}$.

Orbits with $v=5$ and corresponding deflection graph.

Representations of $\delta\phi=\delta\phi(r_{\rm min})$ for several values of $v$.

Observe that in the Lorentzian (resp. Euclidean) Schwarzschild spacetime, the impact parameter $b$ of a massive particle (resp. of any particle) satisfies $bv=L/E$. Moreover, a photon in the usual Schwarzschild metric has $b=L/E$ and this motivates the following terminology: an orbit in a Euclidean Schwarzschild spacetime with velocity at infinity $v$ will be called of sub-photon type (resp. of sup-photon type) if $v\lesssim1$ (resp. $v\gtrsim1$).

We have numerically observed that sub-photon orbits have a negative deflection (i.e. are repelled by the central mass) and that the absolute value of $\delta\phi$ behaves like in the usual Schwarzschild spacetime. This can also be noticed using the lowest order approximations from (<ref>) and (<ref>) recalled at the beginning of the subsection.

For orbits with fixed velocity $v>1$, one could use numerical methods on the expression (<ref>) to find an approximation of the critical value $\rho_0$ (resp. $\rho_{\rm max}$) such that the deflection angle of the orbit with perihelion $r_{\rm min}=\rho_0$ (resp. $r_{\rm min}=\rho_{\rm max}$) vanishes (resp. is maximal). However, we can also use the estimate (<ref>). Indeed, as it is a third order approximation of the deflection, to solve $\delta\phi=0$ at this order amounts to solve a quadratic equation. Similarly, we can differentiate (<ref>) with respect to $r_{\rm min}$ to obtain an approximation of $\rho_{\rm max}$ using the quadratic formula again. This is precisely the reason why we chose this order at the first place. Though the resulting expressions for $\rho_0$ and $\rho_{\rm max}$ are not very enlightening, it is worth noticing that they both diverge when $v\to1$. This suggests to expand these expressions in powers of $v-1$, yielding estimations of $\rho_0$ and $\rho_{\rm max}$ for orbits of sup-photon type. We find

\begin{equation}\label{Rm0}
\frac{\rho_0}{M}=\frac{9\pi}{16}(v-1)^{-1}+1-\frac{21\pi}{32}+\frac{64}{9\pi}+\bigo(v-1)\approx\frac{1.767}{v-1} +1.202+\bigo(v-1)
\end{equation}

\begin{equation}\label{Rm1}
\frac{\rho_{\rm max}}{M}=\frac{9\pi}{8}(v-1)^{-1}+1-\frac{21\pi}{16}+\frac{32}{3\pi}+\bigo(v-1)\approx\frac{3.534}{v-1}+0.272+\bigo(v-1).
\end{equation}

Observe that, at lowest order, we have the remarkable relation $\rho_{\rm max}=2\rho_0$. Moreover, at this order, the maximal value of the deflection is given by

\begin{equation}\label{deltamax}\tag{$\delta_{\rm max}$}
\delta\phi_{\rm max}\approx\frac{16(v-1)}{9\pi}\left(1-\frac{1}{v^2}\right).
\end{equation}

However, this estimate badly fails when $v\gg1$. We illustrate our approximations in the Figure <ref>.

Accuracy of the approximations (<ref>) and (<ref>) for several velocities.

§.§ Deflection of orbits passing through a given point

Consider now an orbit starting from a given point $(\phi,r)=(0,r_0)$, for some fixed $r_0>2M$, whose velocity vector at this point makes an angle $0<\alpha\le\tfrac{\pi}{2}$ with respect to the radial direction. Thinking of the point $(\phi,r)=(0,r_0)$ as an observer receiving particles coming from infinity, we are interested in the total deflection angle $\delta_\alpha\phi$ between the tangent line to the orbit at $(0,r_0)$ and the asymptotic direction. The situation may be visualized using the Figure <ref>.

(z) at (0,0);
(zp) at (0,4);
(zm) at (0,-4);

[color=red,domain=-1.17*pi/4:1.17*pi/4,scale=1,samples=200] plot (deg():2/sqrt(1.7^2*cos(r)^2-1));

(pu) at (8*1/1.7,8*sqrt(1-1/1.7^2));
(pup) at (8.5*1/1.7,8.5*sqrt(1-1/1.7^2));
(pd) at (8*1/1.7,-8*sqrt(1-1/1.7^2));
(mpu) at (-3.5*1/1.7,-3.5*sqrt(1-1/1.7^2));
(mpd) at (-3.5*1/1.7,3.5*sqrt(1-1/1.7^2));
(a) at (2/sqrt(1.7^2-1),0);
(ap) at (2/sqrt(1.7^2-1),4.5);
(am) at (2/sqrt(1.7^2-1),-4.5);
(f1) at (2*1.7/sqrt(1.7^2-1),0);

(rph) at (2*cosh(1.4)/sqrt(1.7^2-1),2*sinh(1.4));
(rphp) at (2*cosh(1.4)/sqrt(1.7^2-1)-1.05*2*sinh(1.4)/sqrt(1.7^2-1),2*sinh(1.4)-1.05*2*cosh(1.4));
(rphm) at (2*cosh(1.4)/sqrt(1.7^2-1)+0.5*2*sinh(1.4)/sqrt(1.7^2-1),2*sinh(1.4)+0.5*2*cosh(1.4));

(R) at (2*1.7/sqrt(1.7^2-1)+1.5*(2*cosh(1.4)/sqrt(1.7^2-1)-2*1.7/sqrt(1.7^2-1)),0+1.5*(2*sinh(1.4)-0));
(L) at (2*1.7/sqrt(1.7^2-1)-1.5*(2*cosh(1.4)/sqrt(1.7^2-1)-2*1.7/sqrt(1.7^2-1)),0-1.5*(2*sinh(1.4)-0));
(Pb) at (2*1.7/sqrt(1.7^2-1)+1.2*(5.5*1/1.7),0+1.2*(5.5*sqrt(1-1/1.7^2)));
(Pbp) at (2*1.7/sqrt(1.7^2-1)+1.3*(5.5*1/1.7),0+1.3*(5.5*sqrt(1-1/1.7^2)));

[fill=black] (rph) circle (3pt);
[fill=black] (f1) circle (6pt);

[dotted] (mpu)–(pup) (mpd)–(pd);
[very thick] (rph)–(rphp) (rph)–(rphm);

(rph)–(R) (f1)–(L);
[dotted] (f1)–(Pbp);
[ultra thick] (a)–(f1) (f1)–(rph);
[<->,dashed,very thick] (Pb)–(pu);

($(Pb)!0.5!(pu)$) node[right]$b$;
($(f1)!0.5!(rph)$) node[below]$r_0$;
($(a)!1.1!(f1)$) node[below right]$0$;
($(f1)!0.5!(a)$) node[right]$r_{\rm min}$;

(rphm) coordinate (A)  (rph) coordinate (B)  (R) coordinate (C) pic ["$\alpha$",draw,->,angle radius=1.5cm]angle;
(L) coordinate (A)  (f1) coordinate (B)  ($(a)!2.5!(f1)$) coordinate (C) pic ["$\Phi$",draw,->,angle radius=1cm]angle;
($(z)!-0.75!(f1)$) coordinate (A) (z) coordinate (B) (pd) coordinate (C) pic ["$\phi$",draw,->,angle radius=0.8cm]angle;
(mpu) coordinate (A) (z) coordinate (B) (pd) coordinate (C) pic ["$\delta\phi=2\phi-\pi$",draw,->,angle radius=3.2cm]angle;

(De) at (0.36,-0.485);
(de) at (1.46,2);
[fill=black] (De) circle (1.5pt);

[fill=black] (z) circle (1.5pt);

(rph) coordinate (A) (De) coordinate (B) (mpd) coordinate (C) pic ["$\delta_\alpha\phi$",draw,->,angle radius=1.5cm,very thick]angle;

Schematics of an orbit passing through the point $(\phi,r)=(0,r_0)$, with total deflection $\delta_\alpha\phi=\Phi+\alpha+\tfrac12(\delta\phi-\pi)=\alpha-\pi+\phi+\Phi$.

We first proceed as in the first subsection, to obtain a numerically efficient analytic formula for $\delta_\alpha\phi$, using Carlson's integrals. Recall the Weierstrass variable $\wp=M/(2r)-1/12$, satisfying the Weierstrass equation with constants $g_{2,3}=g_{2,3}(v,r_{\rm min})$ given by (<ref>). Recall also that $\wp_{\rm max}=M/(2r_{\rm min })-1/12$ is the maximal value of $\wp$ along the geodesic. We have

\[\phi=\int_{\infty}^{r_{\rm min}}\frac{{\rm d}r}{r^2\sqrt{\frac{1}{b^2}-\frac{1}{r^2}\left(1-\frac{2M}{r}+\frac{2Mr}{L^2}\right)}}=\int_{\wp_{\rm max}}^{-\frac{1}{12}}\frac{{\rm d}p}{\sqrt{4p^3-g_2p-g_3}},\]

\[\Phi=\int_{r_0}^{r_{\rm min}}\frac{{\rm d}r}{r^2\sqrt{\frac{1}{b^2}-\frac{1}{r^2}\left(1-\frac{2M}{r}+\frac{2Mr}{L^2}\right)}}=\int_{\wp_{\rm max}}^{\wp_0}\frac{{\rm d}p}{\sqrt{4p^3-g_2p-g_3}},\]

where $\wp_0:=M/(2r_0)-1/12$, so that

\[\delta_\alpha\phi=\alpha-\pi+\phi+\Phi=\alpha-\pi+\left[\int_{\wp_{\max}}^{-\frac{1}{12}}+\int_{\wp_{\rm max}}^{\wp_0}\right]\frac{{\rm d}p}{\sqrt{4p^3-g_2p-g_3}}.\]

We now introduce the two remaining roots $\wp_{\pm}\ne\wp_{\rm max}$ of the cubic $4p^3-g_2p-g_3$, as well as the constants $e_{{\rm max},\pm}:=-\wp_{{\rm max},\pm}-1/12$ as above and arrive at the expression

\begin{align*}\label{Deltaphi_carlson}\tag{$11$}
\delta_\alpha\phi=&R_F\left(e_{\rm max}+\tfrac{M}{2r_{\rm min}},e_-+\tfrac{M}{2r_{\rm min}},e_++\tfrac{M}{2r_{\rm min}}\right)+R_F\left(e_{\rm max}+\tfrac{M}{2r_0},e_-+\tfrac{M}{2r_0},e_++\tfrac{M}{2r_0}\right) \\
&+\alpha-\pi-2R_F(e_{\rm max},e_-,e_+).
\end{align*}

It may be observed that, at fixed $\alpha\in]0;\pi/2]$, when $r_0\gg0$, we have $\Phi\approx\tfrac{\pi}{2}-\alpha$, so $r_{\rm min}\approx r_0\cos(\Phi)\approx r_0\sin(\alpha)$ and $\delta_\alpha\phi\approx\tfrac12\delta\phi(r_{\rm min}=r_0\sin\alpha)$. On the other hand, at fixed $r_0$, when $\alpha=\tfrac{\pi}{2}$, we have $r_{\rm min}=r_0$ and $\Phi=0$. We thus obtain the boundary conditions

\[\delta_\alpha\phi\stackrel[r_0\to\infty]{}\approx \frac12\delta\phi(r_{\rm min}=r_0\sin\alpha)\text{ and }\delta_{{\pi}/{2}}\phi=\frac12\delta\phi(r_{\rm min}=r_0).\]

In particular, at lowest order,

\[\delta_\alpha\phi\stackrel[r_0\to\infty]{}\approx\frac{M}{r_0\sin\alpha}\left(1-\frac{1}{v^2}\right).\]

We illustrate the situation in the Figure <ref>.

Orbits with $0.05<\alpha<0.5$.

Total deflection $\delta_\alpha\phi$ and deflection at infinity $\delta\phi$, as functions of $\alpha$.

Orbits with $r_0=10M$ and $v=1.5$, and the corresponding graph.

For fixed values of $r_0$ and $v$, we now aim to find the critical angle $\alpha_0$ (resp. $\alpha_{\rm max}$) at which the total deflection $\delta_\alpha\phi$ vanishes (resp. is maximal). First recall the following relation from <cit.> (which is just the Pythagorian theorem in curved space)

\[\tan^2\alpha=\left.\frac{p_\phi\dot{\gamma}^\phi}{p_r\dot{\gamma}^r}\right|_{r=r_0}=\left.r^2\left(1-\frac{2M}{r}\right)\left(\frac{{\rm d}\phi}{{\rm d}r}\right)^2\right|_{r=r_0}=\frac{1-\frac{2M}{r_0}}{\frac{r_0^2}{b^2}-1+\frac{2M}{r_0}-\frac{2Mr_0}{L^2}}\]

\begin{equation}\label{tanalpha}
\alpha=\arctan\sqrt{\frac{1-\frac{2M}{r_0}}{\frac{r_0^2}{b^2}-1+\frac{2M}{r_0}-\frac{2Mr_0}{L^2}}}.
\end{equation}

The value of the impact parameter $b$ and the angular momentum $L$ are obtained using the relations (<ref>) and $L=bv(1+v^2)^{-1/2}$. Therefore, we obtain the value of $\alpha$ as a function of the triple $(v,r_0,r_{\rm min})$. With the notation of the previous subsection, we are thus interested in the values of $\alpha$ at $(v,r_0,\rho_0)$ and $(v,r_0,\rho_{\rm max})$.

For sup-photon type orbits with high perihelion, we may inject the approximations (<ref>) and (<ref>) into (<ref>), expand in powers of $r_0$ up to order 2, and then expand in powers of $v-1$ up to order 0 and obtain the following estimates

\begin{equation}\label{alpha0}
\alpha_0=\frac{9\pi M}{16r_0}\left(1+\frac{M}{r_0}\right)(v-1)^{-1}+\frac{M}{r_0}\left(\frac{64}{9\pi}-\frac{21\pi}{32}\right)-\frac{M^2}{r_0^2}\left(\frac{57\pi}{32}-\frac{64}{9\pi}\right)+\bigo(v-1,r_0^{-3})
\end{equation}

\begin{equation}\label{alphamax}
\alpha_{\rm max}=\frac{9\pi M}{8r_0}\left(1+\frac{M}{r_0}\right)(v-1)^{-1}-\frac{M}{r_0}\left(\frac{21\pi}{16}-\frac{32}{3\pi}\right)-\frac{M^2}{r_0^2}\left(\frac{57\pi}{16}-\frac{32}{3\pi}\right)+\bigo(v-1,r_0^{-3}).
\end{equation}

We illustrate our approximations in the Figure <ref>.

Accuracy of the approximations (<ref>) and (<ref>) for several velocities and $r_0=30M$.

Finally, concerning the observable sizes $R_0$ and $R_{\rm max}$ of the null and maximal deflection rings, we have $R_0=r_0\tan\alpha_0$ and $R_{\rm max}=r_0\tan\alpha_{\rm max}$. For an orbit of sup-photon type with $r_0\gg0$, we may use the relation (<ref>) and the estimates (<ref>) and (<ref>) to obtain

\begin{equation}\label{R0}
R_0=\frac{9\pi M}{16}\left(1+\frac{M}{r_0}\right)(v-1)^{-1}+{M}\left(\frac{64}{9\pi}-\frac{21\pi}{32}\right)-\frac{M^2}{r_0}\left(\frac{57\pi}{32}-\frac{64}{9\pi}\right)+\bigo(v-1,r_0^{-2})
\end{equation}

\begin{equation}\label{Rmax}
R_{\rm max}=\frac{9\pi M}{8}\left(1+\frac{M}{r_0}\right)(v-1)^{-1}-{M}\left(\frac{21\pi}{16}-\frac{32}{3\pi}\right)-\frac{M^2}{r_0}\left(\frac{57\pi}{16}-\frac{32}{3\pi}\right)+\bigo(v-1,r_0^{-2}).
\end{equation}

Observe that this is consistent with the approximations (<ref>) and (<ref>) since, at fixed velocity, we have $\lim_{r_0\to\infty}\alpha_0=0$ and thus $R_0=r_0\tan\alpha_0\sim r_0\alpha_0$ and similarly, $R_{\rm max}\sim r_0\alpha_{\rm max}$.

We may proceed similarly for the charged case where $Q\ne0$, but the calculations are a bit tougher. At the lowest order, we have

\[\rho_{\rm max}\stackrel[v\to1]{}\sim\frac{\pi M}{8}\left(9-\frac{Q^2}{M^2}\right)(v-1)^{-1},\]

as well as

\[\alpha_{\rm max}\stackrel[\substack{r_0\to\infty \\ v\to1}]{}\sim\frac{\pi M}{8r_0}\left(9-\frac{Q^2}{M^2}\right)\left(1+\frac{M}{r_0}\right)(v-1)^{-1}\]

\[R_{\rm max}\stackrel[\substack{r_0\to\infty \\ v\to1}]{}\sim\frac{\pi M}{8}\left(9-\frac{Q^2}{M^2}\right)\left(1+\frac{M}{r_0}\right)(v-1)^{-1}.\]

Furthermore, just as in the Schwarzschild case, at lowest order, we have $\rho_{\rm max}=2\rho_0$ and similarly for $\alpha_0$ and $R_0$. Note that for these approximations to make sense, the velocity at infinity $v=\sqrt{E^{-2}-1}$ should be real and in view of the Proposition <ref>, this is ensured only when $Q^2<M^2$.

§ GEODESIC MOTION IN BERTOTTI–ROBINSON SPACETIMES

In this section we investigate the dynamics in Euclidean and Lorentzian Bertotti–Robinson spacetimes and provide in particular analytic solutions for the geodesic equation. For $m,q\in\R$, consider the line element in Binet radial variable

\begin{equation}\label{BR}\tag{BRb}
{\rm d}s^2=Q^2\left[\epsilon \Delta(u){\rm d}\tau^2+\Delta(u)^{-1}{\rm d}u^2+{\rm d}\Omega^2\right],
\end{equation}

where $\Delta(u):=u^2-2mu+q^2$. This is an electro-vacuum solution of Einstein's equation for $A_\mu=\sqrt{-\epsilon}Qu{\rm d}\tau$.

Let $\gamma=(\tau,u,\theta,\phi)$ be a geodesic with respect to the metric (<ref>), parametrized by an affine parameter $\lambda$. By spherical symmetry, we may assume that $\theta\equiv\pi/2$, in which case the Lagrangian reads

\[\mathcal{L}=\frac12 g_{\mu\mu}(\dot{\gamma}^\mu)^2=\frac{Q^2}{2}\left[\epsilon \Delta(u)\dot{\tau}^2+\Delta(u)^{-1}\dot{u}^2+\dot{\phi}^2\right].\]

Thus, the temporal and angular Euler–Lagrange equations immediately yield constants $J,C\in\R$ such that $\dot{\tau}=C\Delta(u)^{-1}$ and $\dot{\phi}=J$ and so the quantity

\[\mathcal{H}:=2\mathcal{L}=Q^2\left[\Delta(u)^{-1}(\epsilon C^2+\dot{u}^2)+J^2\right]\]

is conserved along $\gamma$. Thus, the proper length $s$ satisfies ${\rm d}s^2=\mathcal{H}{\rm d}\lambda^2$ and we get

\begin{equation}\label{mot_BR}
\left(\frac{{\rm d}u}{{\rm d}s}\right)^2=\frac{\dot{u}^2}{\mathcal{H}}=\Delta(u)\left(\frac{1}{Q^2}-L^2\right)-\epsilon E^2,
\end{equation}

where we denote $L:=J/\sqrt{\mathcal{H}}$ and $E:=C/\sqrt{\mathcal{H}}$. From now on, differentiation with respect to $s$ will be denoted by a dot. Since the right-hand side of (<ref>) is quadratic in $u$, we may solve it explicitly.

Before doing this, we need a technical lemma on ODEs:

For any $(\alpha,\beta)\in\R^*\times\R$, their is a unique smooth maximal non-constant solution of the incomplete initial value problem

\[\left\{\begin{array}{rcl}
\dot{y}^2 & = & \alpha(y^2-\beta^2),\\ y(0) & = & \beta.\end{array}\right.\]

Moreover this solution is globally defined and given, for $s\in\R$, by

\[y(s)=\beta\cosh(s\sqrt{\alpha}).\]

The following argument was communicated to the author by Olivier Goubet. First, up to exchanging $y\leftrightarrow -y$, we may assume that $\beta\ge0$. Let $y$ be a maximal non-constant solution of the problem, defined on a (maximal, open) interval $I$. Consider the open subset

\[U:=\{s\in I~|~y(s)>\beta\}\subset\R^*\]

and suppose that $]a;b[\subset U$ is a connected component of $U$, of finite length. By continuity of $y$, we have $\beta\le \lim_{s\to a^+}y(s)=y(a)\le\beta$, so that $y(a)=\beta=y(b)$ and by Rolle's theorem, we may choose $a<t<b$ such that $0=y'(t)=\alpha(y(t)^2-\beta^2)$. This leads to $y(t)\le |y(t)|=\beta$: a contradiction.

This proves that $U$ has no connected component of finite length, leaving four possibilities:

* If $U=\emptyset$, then $y$ is constant.

* If $U=\R_+^*$, then $y=\beta$ on $\R_-$ and solving the well-defined ODE $\dot{y}=\sqrt{\alpha(y^2-\beta^2)}$ on $U$ leads to

\[y(s)=\left\{\begin{array}{cc}\beta\cosh(s\sqrt{\alpha}) & \text{if $s>0$}, \\ 1 & \text{otherwise}.\end{array}\right.\]

* The case $U=\R_-^*$ is the symmetric of the previous one.

* if $U=\R^*$, then we solve $\dot{y}=\pm\sqrt{\alpha(y^2-\beta^2)}$ on $\R^*_\pm$ and obtain the stated solution, the only one being non-constant and smooth.

We may now prove the main result:

Let $\gamma=(\tau,u,\phi)$ be an exterior equatorial geodesic in the spacetime (<ref>), with proper angular momentum $L=\dot{\phi}(0)$. The Binet component $u$ of $\gamma$, as a function of the proper length $s$, is given by

\[u(s)=m+(u(0)-m)\cos(s\sqrt{\ell})+\dot{u}(0)\frac{\sin(s\sqrt{\ell})}{\sqrt{\ell}},\]

where $\ell:=L^2-1/Q^2$.
Moreover, the only circular orbit has $u\equiv m$.

Let $u_0:=u(0)$ and $\dot{u}_0:=\dot{u}(0)$ and observe first that, up to dilating the affine parameter, we may assume that $\left.\left(\frac{{\rm d}\tau}{{\rm d}\lambda}\right)\right|_{\lambda=0}=1$. Then, we have

\[\mathcal{H}=Q^2\left[\frac{\epsilon\left(\Delta(u_0)\left.\left(\frac{{\rm d}\tau}{{\rm d}\lambda}\right)\right|_{\lambda=0}\right)^2+\left.\left(\tfrac{{\rm d}u}{{\rm d}\lambda}\right)\right|_{\lambda=0}^2}{\Delta(u_0)}+J^2\right]=Q^2\left[\epsilon \Delta(u_0)+\mathcal{H}\frac{\dot{u}_0^2+\Delta(u_0)L^2}{\Delta(u_0)}\right],\]

\[\mathcal{H}=\frac{\epsilon Q^2\Delta(u_0)^2}{\Delta(u_0)(1-Q^2L^2)-Q^2\dot{u}_0^2}~\Longrightarrow~E^2=\frac{C^2}{{\mathcal{H}}}=\frac{\Delta(u_0)(1-Q^2L^2)-Q^2\dot{u}_0^2}{\epsilon Q^2}=-\frac{\ell \Delta(u_0)+\dot{u}_0^2}{\epsilon}\]

and thus the equation (<ref>) becomes

\begin{equation}\label{new_mot_BR}
\dot{u}^2=\ell(\Delta(u_0)-\Delta(u))+\dot{u}_0^2=\ell(u_0-u)(u_0+u-2m)+\dot{u}_0^2.
\end{equation}

If $\dot{u}_0=0$, then the Lemma <ref> applied to $y=u-m$, $\beta=u_0-m$ and $\alpha=-\ell$ leads to the stated expression for $u$. Otherwise, on a neighbourhood of $0$, we have $\dot{u}\ne0$ and differentiating (<ref>) with respect to $s$ yields

\[\ddot{u}=-\frac{\ell \Delta'(u)}{2}=\ell(m-u),\]

a linear ODE of order 2 whose solution reads

\[u(s)=m+(u_0-m)\cos(s\sqrt{\ell})+\dot{u}_0\frac{\sin(s\sqrt{\ell})}{\sqrt{\ell}},\]

and this function indeed satisfies (<ref>) and is globally defined.

Now, if $\gamma$ is circular, then $\dot{u}_0=0$ and by the Lemma <ref> again, we find that $u(s)=m+(u_0-m)\cos(s\sqrt{\ell})$ is constant, so that $\ell=0$ or $u_0=m$. This can also be seen by analysing the potential $V=\sqrt{\ell\Delta}$. But suppose now that $u\ne m$, then $\ell=0$ and since $\dot{u}(0)=0$, we get

\[\ell=0\Longleftrightarrow L^2=\frac{1}{Q^2}\Longleftrightarrow \frac{Q^2\dot{\phi}_0^2}{\mathcal{H}}=1\Longleftrightarrow\Delta(u_0)=0,\]

contradicting the fact that $\gamma$ is exterior.

This proposition implies a couple of results, which we summarize in the following corollary:

* The calculations of the previous proof allow to re-write the motion equation (<ref>) as

\[\left(\frac{{\rm d}u}{{\rm d}s}\right)^2=(1/Q^2-\dot{\phi}_0^2)(u-u_0)(u+u_0-2m)+\dot{u}_0^2.\]

In particular, the proper spatial dynamics does not depend on the parameter $q$, nor on the signature of the metric (<ref>). Moreover, the translated variable $v:=u-m$ satisfies the equation

\[\dot{v}^2+\ell v^2=\ell(u_0-m)^2+\dot{v}_0^2.\]

Therefore, in the non-degenerate case $\ell\ne0$, the phase portrait in Binet variable is either an ellipse if $\ell>0$, and a hyperbola otherwise.

* A non-circular geodesic $\gamma=(\tau,u,\theta,\phi)$ with angular momentum $L=\dot{\phi}(0)$ is bounded if and only if $Q^2L^2>{1}$, in which case it has periodic radial component with proper period

\[\omega_\gamma=\frac{2\pi Q}{\sqrt{Q^2L^2-1}}=\frac{2\pi}{L}+\bigo\left(\frac{1}{Q^2L^3}\right).\]

When $Q^2L^2=1$, the geodesic is affine in proper length and when $Q^2L^2<1$, we have

\[u(s)\stackrel[s\to\infty]{}=\bigo\left(e^{s\sqrt{1-Q^2L^2}}\right),\]

so the Binet variable $u$ blows-up exponentially in proper length.

* There are exterior circular orbit only when $q^2>m^2$, in which case the only such orbit has $r=1/m$. In particular, there are no circular orbits in the original Bertotti–Robinson space <cit.>.

* Any (exterior) Bertotti–Robinson spacetime is geodesically complete in Binet variable.

From this we can deduce the analytic expression for the geodesic motion in terms of affine parameter. We arrive at the following result:

Let $\gamma=(\tau,r,\phi)$ be an exterior equatorial geodesic in the spacetime (<ref>), parametrized by the affine parameter $\lambda$, chosen so that $\dot{\tau}_0=1$. Then the constants $\mathcal{H}$ and $L$ are given by

\[\mathcal{H}=Q^2\left(\frac{\epsilon(1-2mr_0+q^2r_0^2)}{r_0^2}+\frac{\dot{r}_0^2}{r_0^2(1-2mr_0+q^2r_0^2)}+\dot{\phi}_0^2\right)\text{ and }L=\frac{\dot{\phi}_0}{\sqrt{\mathcal{H}}},\]

and letting $\ell:=L^2-Q^{-2}$, we have

\[r(\lambda)=\left(m+\left(\frac{1}{r_0}-m\right)\cos(\lambda\sqrt{\ell})-\frac{\dot{r}_0}{r_0^2}\frac{\sin(\lambda\sqrt{\ell})}{\sqrt{\ell\mathcal{H}}}\right)^{-1},\text{ as well as }\phi(\lambda)=\phi_0+\lambda\dot{\phi}_0.\]

Moreover, the geodesic $\gamma$ has energy $E=\frac{1-2mr_0+q^2r_0^2}{r_0^2\sqrt{\mathcal{H}}}$.

§ IMPLEMENTATION OF ORBITS AND RAY-TRACING

We have developed a package[<https://github.com/arthur-garnier/euclidean_orbits_and_shadows.git>] under Python, for drawing orbits in the spacetimes (<ref>) and (<ref>), as well as for ray-tracing the first ones, using a standard backward ray-tracing method (see for instance <cit.>). As mentioned in <ref>, the Carlson and Coquereaux–Grossmann–Lautrup algorithms <cit.> are used to compute geodesics in the Reissner–Nordström instanton.

Let us briefly recall the backward ray-tracing procedure we use for shadowing a spacetime. First, we consider an artificial celestial hemisphere on which we project our original image, seeing it as a portion of its tangent plane parallel to the screen (and on the other side of the black hole). As a projection, we simply choose the standard and widely used equirectangular projection, which has the advantage of taking the celestial hemisphere to a square, which we may rescale to fit our image.

Next, for each pixel of the screen, we consider the geodesic for the spacetime (null in the Lorentzian and normalized in the Euclidean case) starting at this point and with velocity directed by the line from the point observer. We then solve the geodesic equations (backwards) and we see if the ray ends in (came from) the black hole or touches the sphere somewhere. If so, the RGB value of the pixel on the screen is given by the value of the landing pixel on the sphere and we carry this process on until every pixel has been worked out.

To simplify calculations, we make heavy use of the spherical symmetry of the spacetimes considered here: given an initial datum, use a linear rotation to bring the initial velocity (and hence the full orbit) in the plane $\{\theta=\pi/2\}$. Then, we give values to the various constants involved in the expression of the radial geodesic and, instead of computing the full orbit, we simply solve the equation $r=r_{\Sph}$ where $r_{\Sph}$ is the radius of the celestial sphere, using the Weierstrass function. This can be done rather easily, precisely and quickly: we compute some values until we cross the sphere and the first such point is used as an initial value for the Newton method. We finally rotate the result back and find our landing pixel. Thus, no full orbit calculation is required. For more details and illustrations, see <cit.>. We should mention however that since, in contrast to the Lorentzian framework, photons (i.e. null geodesics) do not properly exist in the Euclidean world, we have to trace orbits with some prescribed velocity at infinity $v$: an additional input to the program. For each pixel (i.e each particle), we provide its initial position and the direction of its velocity vector, whose norm is adjust to match the velocity $v$.

We finish our discussion by providing some figures illustrating the programs and our results.

First, concerning the orbits drawings, it should be mentioned that, while the functions of the package are designed to draw orbits in 3D, we chose to plot planar orbits here, to make the figures more readable. In each case, the mass is set to unity and we vary the charge in the different plots, distinguishing between the sub-extremal and sup-extremal cases for the charge. In each figure, the legend gives the values of the energy and angular momentum of each displayed orbit.

The Figure <ref> depicts some orbits in the spacetime (<ref>). In particular, we illustrate the Proposition <ref>: when there is a horizon, the squared energy is smaller than unity and there are circular orbits with arbitrary energy when we have a naked singularity. Notice also the presence of bounded orbits in the sup-extremal case.

In the Figures <ref> and <ref>, we show some Bertotti–Robinson orbits, with $\epsilon=\pm1$, $Q=1/2$ and $m=1$. The Figure <ref> displays orbits in the original Bertotti–Robinson spaces where $m=q=0$.

Finally, some particular unstable orbits may be observed in the sup-extremal Euclidean cases. In the Figure <ref> are displayed some “flower orbits” in a horizon-less ERN spacetime, while in the Figure <ref>, we plot some “star orbits” in a horizon-less BR space. These can be seen as analogues of “leaf orbits” in the usual Lorentzian solutions, see <cit.>.

Regarding the shadows, just as in <cit.>, we use the color grid shown in Figure <ref> as our base picture for the shadows. In the Figure <ref>, we depict the shadows of an ERN spacetime with various charges. The render time of each figure is approximately 350 seconds, on an 8-core 3.00 GHz CPU with 16 Go of RAM[We have tested the efficiency of the shadowing program by running it on (randomly generated) images with $10\times10$ to $500\times500$ pixels and we made an exponential regression on the render time. We found

\[{\rm time}[s]\approx e^{-10.6}\times{\rm pix}^{2.5},\]

with a regression coefficient $r>0.99$.].

To illustrate our results on deflection, we also wrote a code similar to the shadowing program, but which rather displays the total angular deflection $\delta_\alpha\phi$ (here, $\alpha$ is the angle between the initial velocity vector and the radial direction passing through the common converging point of the “light rays”. The results are depicted in the Figures <ref> and <ref>. Observe that, in accordance with the Section <ref>, the deflection is always positive in the usual Schwarzschild spacetime, with a visible event horizon, while in the Euclidean world, the horizon disappears and the presence of a null and maximal deflection rings is manifest.

Lastly, in the Figure <ref>, we give the shadows of ERN spacetimes on a celestial background. The original picture is from the NASA and the render time for each shadow is about 2800 seconds, the resolution of the picture being of $1080$ pixels. While the null deflection ring is impossible to pinpoint on such a picture, the maximal deflection ring is still quite noticeable. Observe moreover that, again in contrast to the usual Schwarzschild metric, the horizon-full case $Q^2<M^2$ and the naked singularity case $Q^2>M^2$ produce similar pictures, the most noticeable difference being the size $R_{\rm max}$ of the maximal deflection ring, which decreases as the charge increases. This is in accordance with the lowest order approximation of $R_{\rm max}$ provided in the Remark <ref>.

Orbits with $Q=0.998$

Orbits with $Q=1.098$

Circular orbits with $Q=1.098$

Some Euclidean Reissner–Nordström orbits with $M=1$.

Orbits in Lorentzian Bertotti-Robinson spacetimes with $Q=1/2$ and $m=1$.

Orbits in Euclidean Bertotti-Robinson spacetimes with $Q=1/2$ and $m=1$.

Orbits in Bertotti–Robinson spaces with $Q=1$ and $m=q=0$.

Flower orbits in ERN space with $M=1$ and $Q=1.098$.

Star orbits in BR space with $\epsilon=-1$, $Q=1/2$, $m=1$ and $q=1.098$.

The original $600\times600$ color grid.

The original $600\times600$ pixels color grid.

Picture obtained with the code ($M=Q=0$, $v=1$).

The color grid used in the program.

Shadows of ERN spacetimes with $M=1$ and $v=1.5$.

Deflection maps in Euclidean Schwarzschild spacetime ($M=1$) with several velocities at infinity.

$v=1$ (photon)

Deflection maps in Lorentzian Schwarzschild spacetime ($M=1$) with several velocities at infinity.

The base picture (Arches cluster).

Shadows of ERN spacetimes with $M=1$ and $v=1.5$ on a celestial background (original image: <https://images.nasa.gov/details/GSFC_20171208_Archive_e000717>).

§ CONCLUSIONS

§.§ Summary

In this work, we study the geodesic motion in the two different kinds of spherically symmetric electro-vacuum Euclidean solutions of Einstein's equation with complex vector potential, namely the Reissner–Nordström and Bertotti–Robinson-like instantons. More precisely, after proving that these are indeed the only two possible such solutions (and that these two are incompatible), we derive the motion equations and give the main properties of the test-particle orbits, using their motion constants such as the energy and angular momentum.

We start with the Euclidean Reissner–Nordström solution, for which we prove that if the spacetime has a horizon, then all exterior (non-constant) geodesics have bounded energy and, otherwise, there are circular orbits with arbitrary energy. In particular, the spacetime features a horizon if and only if there are no bounded orbits. This generalizes the results from <cit.> on the Euclidean Schwarzschild solution and shows that the differences in the dynamics of the Lorentzian and Euclidean solutions is deeper than the apparent mere sign change in the line element. Then, we show that the polar radial motion equation may be put in Weierstrass form, something that we can take advantage of to numerically solve the motion equation. Moreover, this approach shows that, in the special case of the Euclidean Schwarzschild metric, the (polar) phase portrait in Binet variable[specifically, in an affine transform of the Binet variable] describes a real elliptic curve. This allows us to derive a new proof of the results from <cit.>, using the elementary geometry of the curve. Furthermore, we observe that this elliptic curve is always disconnected, while in the usual Lorentzian Schwarzschild solution, it can be either connected or not.

We then study the gravitational bending of geodesics in Euclidean Schwarzschild spacetimes. More precisely, given a geodesic with closest approach radius $r=r_{\rm min}$ and velocity at infinity $v^2=\lim_{r\to\infty}({\rm d}r/{\rm d}\tau)^2$, we first provide an analytic formula for the deflection angle $\delta\phi$ that occurs between the two asymptotic directions of the geodesic, in terms of Carlson's elliptic integrals, which are used to numerically compute the deflection. Then, we give approximations for $\delta\phi$ when $r_{\rm min}\to\infty$ using the previous results of <cit.>. We observe that when $v>1$, there are particular values $\rho_0$ and $\rho_{\rm max}$ for which $\delta\phi=0$ for $r_{\rm min}=\rho_0$ and $\delta\phi$ is maximal for $r_{\rm min}=\rho_{\rm max}$. Estimates for $\rho_0$ and $\rho_{\rm max}$ are provided when $v\gtrsim1$. We then do the same for the deflection $\delta_\alpha\phi$ of geodesics passing through a point at fixed radius $r_0$ and with varying angle $0<\alpha<\pi/2$, modelling orbits coming from infinity to the eye of an observer. We give approximations, when $r_0\gg2M$ and $v\gtrsim1$, for the critical angles $\alpha_0$ and $\alpha_{\rm max}$ corresponding to the null and maximal deflection rings, as well as for the observable size of these rings. We use numerical evaluations of the analytic formulas to check the accuracy of our estimates.

Next, we do the same for the Bertotti–Robinson solution for which, after a technical lemma on ordinary differential equations, we give a general analytic solution, in terms of (hyperbolic) trigonometric functions. In particular, orbits are either periodic or unbounded and there is a unique circular orbit, which is exterior exactly when the metric has no horizon.

Finally, we provide some details on a Python code[available at <https://github.com/arthur-garnier/euclidean_orbits_and_shadows.git>], designed to plot orbits in the aforementioned spacetimes, as well as to draw shadows of the Euclidean Reissner–Nordström family, by the usual backward ray-tracing method. The numerical computation takes advantage of the Weierstrass form of the polar equation, coupled with the well-known Carlson <cit.> and Coquereaux–Grossmann–Lautrup <cit.> algorithms, to produce an efficient and rather fast code. This method is the same as the one used by the author in <cit.> to ray-trace Reissner–Nordström–(anti)de Sitter black holes. We illustrate the code and our results by providing some figures.

§.§ Discussion and perspectives

In the Figures <ref> and <ref>, we depict some particularly shaped unstable periodic orbits, that can be seen as analogues of the so-called “leaf orbits” that occur in the Lorentzian framework <cit.>. However, we have no theoretical interpretation for these orbits and this could be a matter of interest for future works.

In <ref>, we have observed that a Euclidean Schwarzschild orbit can mathematically be interpreted as a Lorentzian space-like Schwarzchild orbit with purely imaginary energy. To decide whether this is merely a mathematical curiosity or is a manifestation of a deeper physical phenomenon is beyond the knowledge of the author and could be an interesting subject to investigate as well.

In <ref>, we give approximate values for the null and maximal deflection rings of a sup-photon type particle (i.e. with velocity at infinity $v\gtrsim1$). Moreover, we have numerically observed that such particular rings occur exactly when $v>1$. It would be suitable to provide a rigorous proof of this fact.

Besides the theoretical results, one aim of this work is to provide an open, transparent, user-friendly and customizable Python code to plot orbits and draw shadows of spherically symmetric (asymptotically flat) electro-vacuum instantons. This has a negative side though, which is that the code is under-optimal and still slow in comparison to other widely-used ray-tracing codes, such as ,  or the more recent  <cit.>. Therefore, it would be interesting to improve it by using a GPU parallelization, for instance.

As in <cit.>, the shadowing program is designed to work with a common plane background image and in order to avoid distortions, we took the compromise of projecting the image on a celestial hemisphere. However, to be able to catch any possible particle that is ray-traced, we had to take the whole celestial sphere into account and we arbitrarily chose to project a mirrored version of the original image onto it. Thus, the code could be completed by producing a panoramic version of it.

§.§ Acknowledgments

I am much grateful to Emmanuele Battista for bringing the subject to my attention and for many fascinating discussions. I also warmly thank him for his careful reading of the present work and his useful comments.

Missing 'biblatex' package
The bibliography requires the 'biblatex' package.

journaltitleClass. Quantum Grav.

journaltitleClass. Quantum Grav.
relatedstringCorrected in
titleGravitational deflection of massive particles in classical and semiclassical general relativity

journaltitleEur. Phys. J. C
titleGeodesic motion in Euclidean Schwarzschild geometry

titleGravitational bending of light near compact objects

journaltitlePhys. Rev.
titleUniform electromagnetic field in the theory of general relativity

titleDetermining the dimensionality of spacetime by gravitational lensing

journaltitleNumerical Algorithms
titleNumerical computation of real or complex elliptic integrals

The Americal Astronomical Society

title: A massively parallel GPU-based code for ray tracing in relativistic spacetimes

Oxford University Press

journaltitleIMA Journal of Numerical Analysis
titleIterative method for calculation of the Weierstrass elliptic function

journaltitleProc. R. Soc. Lond. A
titleThe gravity field of a particle

The American Astronomical Society

titleA fast new public code for computing photon orbits in a Kerr spacetime

journaltitleClass. Quantum Grav.
titleQuantum vacuum energy near a black hole: the Maxwell field

seriesLecture notes in Physics Monographs
titleQuantum gravity, quantum cosmology and Lorentzian geometries

IOP Publishing

journaltitleClass. Quantum Grav.
titleMotion equations in a Kerr–Newman–de Sitter spacetime: some methods of integration and application to black holes shadowing in Scilab

journaltitlePubl. Astron. Dep. Eotvos Univ.
titleWeak gravitational lensing in brane-worlds

Cambridge University Press

booktitleGeneral relativity - An Einstein centenary survey
titleQuantum field theory in curved spacetime

American Physical Society

journaltitlePhys. Rev. D
titleAction integrals and partition functions in quantum gravity

IOP Publishing

journaltitleClass. Quantum Grav.
titleThe application of Weierstrass elliptic functions to Schwarzschild null geodesics

journaltitleJapanese Journal of Astronomy and Geophysics
titleTheory of the relativistic trajectories in a gravitational field of Schwarzschild

American Physical Society

journaltitlePhys. Rev. D
titlePath-integral derivation of black-hole radiance

booktitleRecent Developments in Gravitation: Cargèse 1978
titleEuclidean Quantum Gravity

journaltitlePhys. Lett. A
titleGravitational instantons

journaltitleEur. Phys. J. C
titleGravitational deflection of massive particles in Schwarzschild–de Sitter spacetime

journaltitleThe Quarterly Journal of Mathematics
titleOn the spherically symmetric field in relativity

American Physical Society

journaltitlePhysical Review D
titleA periodic table of black hole orbits

journaltitleCommun. Theor. Phys.
titleGravitational deflection of massive particles by a Schwarzschild black hole in radiation gauge

American Book Company

seriesThe Cornell Mathematical Series
titleElements of the differential calculus

IOP Publishing

journaltitleClass. Quantum Grav.
titleBlack holes and gravitational instantons

American Physical Society

journaltitlePhys. Rev. D
titleNegative modes and the thermodynamics of Reissner–Nordström black holes

journaltitleClass. Quantum Grav.
titleGravitational lensing of massive particles in Reissner–Nordström black hole spacetime

The American Astronomical Society

title: A public GPU-based code for general-relativistic radiative transfer in Kerr spacetime

journaltitleBull. Acad. Polon. Sci
titleA solution of the Maxwell-Einstein equations

editionSecond edition
seriesGraduate texts in Physics
titleGeneral Relativity

journaltitleEur. Phys. J. C
title: a new code for ray tracing around compact objects

IOP Publishing

journaltitleClass. Quantum Grav.
title: a new general relativistic ray-tracing code