Gravitational birefringence of light in Schwarzschild spacetime

TL;DR

This paper investigates how light's polarization affects its trajectory in Schwarzschild spacetime, revealing polarization-dependent offsets and instabilities through numerical and perturbative solutions of modified equations.

Contribution

It derives and solves the Souriau-Saturnini equations in Schwarzschild spacetime to analyze gravitational birefringence effects on light trajectories.

Findings

Polarization-dependent trajectory offsets grow over time.

An instability in the transverse spin component is observed.

Numerical and perturbative methods confirm birefringence effects.

Abstract

We compute the gravitational birefringence of light as it undergoes gravitational lensing. To this end we re-derive the Souriau-Saturnini equations in the Schwarzschild metric and solve them numerically and perturbatively. Our main result is an offset between the trajectories of the photons of opposite polarisations, which grows with time. We also find an intriguing instability of the spin component transverse to the momentum.

Figures

Equations

$${\dot{P}}^{\mu }$$

$${\dot{S}}^{\mu \nu }$$

$${R(S{)}^{\mu }}_{\nu }:={R^{\mu }}_{\nu \alpha \beta }{S}^{\alpha \beta }\mathrm{and}R(S)(S):=R_{\mu \nu \alpha \beta }{S}^{\mu \nu }{S}^{\alpha \beta }.$$

$$\dot{X}$$

$$\dot{P}$$

$$\dot{S}$$

$$\mathrm{g}=-B^2\Vert d\mathbf{x}{\Vert }^2+A^{2}d{t}^2$$

$$A:=\frac{r-a}{r+a},B:={\left(\frac{r+a}{r}\right)}^2,r:=\sqrt{\mathbf{x}\cdot \mathbf{x}},0<a<r.$$

$$\rho =r{\left(1+\frac{GM}{2r}\right)}^2\text{or}r=\frac{1}{2}\left(\rho -GM+\sqrt{\rho (\rho -2GM)}\right).$$

$${{\Gamma }^j}_{ii}=-{{\Gamma }^i}_{ji}=-{{\Gamma }^j}_{jj}=\frac{2a{x}^j}{r^2(r+a)},{{\Gamma }^j}_{44}=\frac{2a{r}^3(r-a)x^j}{(r+a{)}^7},{{\Gamma }^4}_{4j}=\frac{2a{x}^j}{r(r+a)(r-a)},$$

$${R^i}_{jij}$$

$${R^4}_{i4i}$$

$$P=\left(\begin{array}[]{c}\displaystyle\frac{{\mathbf{p}}}{B}\\[10.0pt] \displaystyle\frac{\|{\mathbf{p}}\|}{A}\end{array}\right)$$

$$S=({S^{\mu}}_{\nu})=\left(\begin{array}[]{cc}j({\mathbf{s}})&\displaystyle-\frac{({\mathbf{s}}\times{\mathbf{p}})}{\|{\mathbf{p}}\|}\frac{A}{B}\\[6.0pt] \displaystyle-\frac{({\mathbf{s}}\times{\mathbf{p}})^{T}}{\|{\mathbf{p}}\|}\frac{B}{A}&0\end{array}\right).$$

$$-\frac{1}{2}\mathrm{Tr}(S^2)=s^2,$$

$$s:=\frac{\mathbf{s}\cdot \mathbf{p}}{\Vert \mathbf{p}\Vert },$$

$$\mathrm{Pf}(R(S))$$

$$R(S)(S)$$

$$SR(S)P$$

$$\mathbf{c}$$

$$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad+\left[{\mathbf{x}}\times{\mathbf{p}}\cdot{\mathbf{s}}\right]\left({\mathbf{x}}\cdot\frac{{\mathbf{p}}}{\|{\mathbf{p}}\|}\right){\mathbf{s}}\times\frac{{\mathbf{p}}}{\|{\mathbf{p}}\|}\Bigg{]}\,,$$

$$d$$

$$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad+\|{\mathbf{p}}\|\,(\|{\mathbf{s}}\|^{2}-s^{2})\,r^{2}-2\left[{\mathbf{x}}\times{\mathbf{p}}\cdot{\mathbf{s}}\right]^{2}/\|{\mathbf{p}}\|\Big{]}.$$

$${\left[\mathbf{u}\times \mathbf{v}\cdot \mathbf{w}\right]}^2$$

$$\Psi (Z)=P_{\mu }{Z}^{\mu }+\frac{1}{2}{S}^{\mu \nu }{\nabla }_{\mu }{Z}_{\nu }$$

$$\mathcal{E}=\frac{r-a}{r+a}\Vert \mathbf{p}\Vert +\frac{2ar}{(r+a{)}^4\Vert \mathbf{p}\Vert }\left[\mathbf{x}\times \mathbf{p}\cdot \mathbf{s}\right],$$

$$\mathcal{L}={\left(\frac{r+a}{r}\right)}^2\mathbf{x}\times \mathbf{p}+\frac{r-a}{r+a}\mathbf{s}+\frac{2a}{r^2(r+a)}(\mathbf{s}\cdot \mathbf{x})\mathbf{x},$$

$$\mathbf{s}=\frac{r+a}{r-a}\left[\mathcal{L}-{\left(\frac{r+a}{r}\right)}^2\mathbf{x}\times \mathbf{p}-\frac{2a}{r^2(r+a)}(\mathcal{L}\cdot \mathbf{x})\mathbf{x}\right].$$

$$D:=r^2(\mathbf{s}\cdot \mathbf{p})-3(\mathbf{p}\cdot \mathbf{x})(\mathbf{s}\cdot \mathbf{x}).$$

$$\frac{dt}{d\tau }=\frac{r+a}{r-a}\Vert \mathbf{p}\Vert \left[\frac{sD\Vert \mathbf{p}\Vert }{s^2{r}^2\Vert \mathbf{p}{\Vert }^2-3(\mathbf{s}\cdot \mathbf{x}{)}^2\Vert \mathbf{p}{\Vert }^2+3[\mathbf{x}\times \mathbf{p}\cdot \mathbf{s}{]}^2}\right],$$

$$\frac{d\mathbf{x}}{dt}$$

$$\frac{d\mathbf{p}}{dt}$$

$$\displaystyle\quad\quad\quad\quad-r\,\|{\mathbf{p}}\|^{2}\big{[}D+r\,(r-a)\,({\mathbf{s}}\cdot{\mathbf{p}})\big{]}\,{\mathbf{x}}$$

$$\displaystyle\quad\quad\quad\quad\quad\quad+3\,(r-a)[{\mathbf{x}}\times{\mathbf{p}}\cdot{\mathbf{s}}]\,({\mathbf{p}}\cdot{\mathbf{x}})\,{\mathbf{x}}\times{\mathbf{p}}\Bigg{\}},$$

$$\frac{d\mathbf{s}}{dt}$$

$$\displaystyle+\left(2\|{\mathbf{p}}\|^{2}({\mathbf{x}}\cdot{\mathbf{s}})-({\mathbf{x}}\cdot{\mathbf{p}})({\mathbf{s}}\cdot{\mathbf{p}})\right){\mathbf{x}}\times{\mathbf{p}}\big{]}+2arD\left(({\mathbf{x}}\cdot{\mathbf{s}}){\mathbf{p}}-({\mathbf{x}}\cdot{\mathbf{p}}){\mathbf{s}}\right)$$

$$\displaystyle+2a(r-a)\big{(}-r^{2}({\mathbf{s}}\cdot{\mathbf{p}})^{2}{\mathbf{x}}-3[{\mathbf{x}}\times{\mathbf{p}}\cdot{\mathbf{s}}]^{2}{\mathbf{x}}$$

$$\displaystyle\quad\quad+r^{2}({\mathbf{x}}\cdot{\mathbf{s}})({\mathbf{s}}\cdot{\mathbf{p}}){\mathbf{p}}+3[{\mathbf{x}}\times{\mathbf{p}}\cdot{\mathbf{s}}]({\mathbf{x}}\cdot{\mathbf{s}}){\mathbf{x}}\times{\mathbf{p}}\big{)}\Big{\}}.$$

$$\mathbf{x}\times \mathbf{p}\cdot \mathbf{s}$$

Full text

**Gravitational birefringence of light

in Schwarzschild spacetime111supported by the OCEVU Labex (ANR-11-LABX-0060) funded by the ”Investissements d’Avenir”

 French government program **

Christian Duval, Loïc Marsot222mailto: marsot.loic@gmail.com and Thomas Schücker333mailto: thomas.schucker@gmail.com

Aix Marseille Univ, Université de Toulon, CNRS, CPT, Marseille, France

Abstract

We compute the gravitational birefringence of light as it undergoes gravitational lensing. To this end we re-derive the Souriau-Saturnini equations in the Schwarzschild metric and solve them numerically and perturbatively. Our main result is an offset between the trajectories of the photons of opposite polarisations, which grows with time. We also find an intriguing instability of the spin component transverse to the momentum.

To the memory of Christian Duval

1 Introduction

Birefringence of light is a well kown phenomenon in anisotropic matter like quartz or calcite. On theoretical grounds Federov [1] in 1955 and Imbert [2] in 1972 predicted birefringence in absence of matter but in an electric field with a gradient. This prediction was confirmed experimentally in 2008 by Bliokh [3] et al. and Hosten & Kwiat [4]. They measured an offset between the trajectories of the photons of opposite polarisations, which is of the order of the wavelength of the photons.

Birefringence of light without matter arises in Loop Quantum Gravity [5] and more generally in extensions of Maxwell’s theory or of the $SU(3)\times SU(2)\times U(1)$ model, which violate Lorentz invariance, like the addition of a Chern-Simons term, [6, 7]. Birefringence in these theories is strongly constrained by spectropolarimetry of radio galaxies, quasars, Gamma Ray Bursts and the Cosmic Microwave Background [8, 9, 10, 11, 12].

The starting point of the present work is less ambitious: the Mathisson-Papapetrou-Dixon equations [13, 14, 15] and their geometric derivation in 1974 by Souriau [16]. He and Saturnini [17], using this geometric framework, generalized the null geodesic equation in order to include the spin of photons, thus proposing birefringence of light in a gravitational field with a gradient. In 1976 Saturnini computed this generalized equation in the Schwarzschild metric and obtained first numerical solutions with birefringence.

The usual geodesic equation can be derived exactly from general principles, when neglecting dipole and higher moments of the test particle. See the derivation using conserved quantities by Dixon [15] or, more geometrically by Souriau [16]. This computation can be generalized to include the dipole moment of the test particle. This is how the Mathisson-Papapetrou-Dixon (MPD) equations are derived. They can also be obtained using Souriau’s geometrical approach. With $X=(X^{\mu}(\tau))$ denoting the trajectory of the test particle, $P$ its 4-momentum, and $S$ a skewsymmetric tensor describing its spin state, we have the MPD equations,

[TABLE]

where the dot over the trajectory $X$ denotes the ordinary derivative with respect to its affine parameter, $\dot{X}=dX/d\tau$, while the dot over $P$ and $S$ denotes the covariant derivative with respect to that same parameter.

These equations are very general and can be applied to any kind of test particle, be it a planet with its intrinsique angular momentum or an elementary particle with its spin. However, these equations are not deterministic: there are more unknowns than equations. We need to impose constraints, or “equations of state”, that depend on the test particle we want to describe. In flat spacetime, at least for photons, ${S^{\mu}}_{\nu}P^{\nu}={S^{\mu}}_{\nu}\dot{X}^{\nu}=0$ holds. Souriau shows [18] that the MPD equations together with the previous constraint in flat spacetime for massless particles of spin 1 lead to the Maxwell equations after geometric quantization. Now, in curved spacetime the MPD equations entail that the 4-momentum $P$ may not be parallel to the 4-velocity $\dot{X}$, thus there are different possible equations of state. For extended particles, it is generally accepted to impose the Tulczyjew constraint, ${S^{\mu}}_{\nu}P^{\nu}=0$ [19], to describe uniquely the worldline of the center of mass of a particle [15]. However, for elementary particles, the center of mass gives no criterion. An alternative possibility for the equation of state would then be the Frenkel-Pirani constraint ${S^{\mu}}_{\nu}\dot{X}^{\nu}=0$ [20, 21]. For examples with the Frenkel-Pirani constraint, see [22, 23, 24, 25].

In this paper, we consider the Tulczyjew constraint since this framework can be successfully applied for elementary particles in different situations: in the electromagnetic field, this constraint yields the well known and well used Bargmann-Michel-Telegdi equations [26], with an anomalous velocity, and with correct anomalous magnetic moment [27]; it has also been used [28, 29, 30, 31] to correctly describe the above mentioned birefringence in the Fedorov-Imbert effect.

There are also examples in gravitational fields, with massive particles, see [32, 25]. Let us insist that solving the massless equations in this framework is more involved than for massive equations.

Reference [33] presents the Souriau-Saturnini equations in a generic Robertson-Walker metric and some numerical and some perturbative solutions, showing a striking effect on the photon’s trajectory. Indeed, it travels on a helix, centered around the null geodesic and with a radius of the order of the wavelength.

Here we take up Saturnini’s work in the (outer) Schwarzschild metric. Thanks to present day computing power, we obtain numerical solutions precise enough to lead us to perturbative solutions. These solutions, contrary to the ones in a Robertson-Walker metric, feature an intriguing instability of the spin component transverse to the momentum and an offset between the trajectories of photons with opposite polarisations, which grows linearly with time.

In 2006 Gosselin, Bérard & Mohrbach [34] have published an analysis similar to ours, but starting from the Bargmann-Wigner equations. We compare results at the end of section 5.

Let us note an experimental upper bound on birefringence obtained in 1974 by Very Long Baseline Interferometry of radio sources lensed in the Sun’s gravitational field [35].

The paper is organized as follows. We will start with the Souriau-Saturnini equations in section 2, which are the application of the Mathisson-Papapetrou-Dixon equations with the Tulczyjew constraint to a massless particle of spin 1, in the case of a Schwarzschild spacetime. Section 3 has two aims: to show the similarity between the usual geodesic equations and the equations obtained by taking the photon’s spin into consideration; and to present a different approach to the well known gravitational lensing in Schwarzschild spacetime based on the conserved quantities derived here and on first order differential equations. The following section 4 will be devoted to the numerical study of the complicated equations of motion with spin obtained from section 2. In light of the numerical integrations and the similarity of the spinning equations and the geodesic equations, the last section 5 presents a perturbative solution to the Souriau-Saturnini equations in Schwarzschild spacetime, which will allow us to interpret the differences introduced by the spin of the photon in gravitational lensing.

2 Spinning massless particles

2.1 The Souriau-Saturnini equations

We use the shorthand notations,

[TABLE]

We also suppress indices by using linear maps, i.e. $S=({S^{\mu}}_{\nu})$, and we write $P=(P^{\mu})$. Assuming the consistency condition $R(S)(S)\neq 0$, the equations of motion of photons in space-time [17] read:

[TABLE]

Here $\overline{P}$ denotes the covector associated to the vector $P$ via the metric: $\overline{P}_{\mu}\vcentcolon=g_{\mu\rho}P^{\rho}$. We denote the Pfaffian of a skewsymmetric linear map $F$ as $\mathrm{Pf}(F)$. It is the square root of the determinant of the linear map, noting that the determinant of a skewsymmetric matrix can always be written as a perfect square. An alternative definition of the Pfaffian is $\mathrm{Pf}(F)=-\frac{1}{8}\sqrt{-\det(\mathrm{g}_{\alpha\beta})}\,\varepsilon_{\mu\nu\rho\sigma}F^{\mu\nu}F^{\rho\sigma}$ with $\varepsilon_{\mu\nu\rho\sigma}$ the Levi-Civita symbol such that $\varepsilon_{1234}=1$.

For a derivation of these equations (2.4 - 2.6) in English, see [33].

2.2 Metric

The Schwarzschild metric can be expressed in an isotropic coordinate patch $(X^{\mu})=({\mathbf{x}},t)$ by

[TABLE]

with

[TABLE]

where ${\mathbf{x}}=(x^{1},x^{2},x^{3})$ and $\|\,\cdot\,\|$ is the Euclidean norm. If $(\rho,\theta,\varphi,t)$ are the Schwarzschild coordinates, the isotropic polar ones $(r,\theta,\varphi,t)$ are related by

[TABLE]

Recall that $a=\frac{1}{2}GM$ is the Schwarzschild radius.

The vector product, which abounds in computations involving spin, take a simple form in isotropic coordinates. Therefore we adopt these coordinates. However they have the drawback, that including the cosmological constant, which is straight forward in the Schwarzschild coordinates, becomes difficult. This inclusion will be dealt with in a future publication.

We have the following Christoffel symbols,

[TABLE]

for all $i\not=j=1,2,3$, no summation over repeated indices.

For the Riemann tensor ${R^{\mu}}_{\nu\alpha\beta}=\partial_{\alpha}{\Gamma^{\mu}}_{\beta\nu}-\partial_{\beta}{\Gamma^{\mu}}_{\alpha\nu}+\cdots$ with $i,j$ and $k$ all different, we have

[TABLE]

The Ricci tensor vanishes.

2.3 Momentum and spin

In the above coordinate system, the (future pointing) 4-momentum of the photon is written as

[TABLE]

with ${\mathbf{p}}\in\mathbb{R}^{3}\setminus\{0\}$, the spatial linear momentum, and $\|{\mathbf{p}}\|\vcentcolon=\sqrt{{\mathbf{p}}\cdot{\mathbf{p}}}$ (Euclidean scalar product). We suppose positive energy, $\|{\mathbf{p}}\|>0$. The 4-momentum is light-like, $P^{2}=0$.

The map $S$ is skewsymmetric with respect to the metric: $\mathrm{g}(SV,W)=-\mathrm{g}(V,SW)$ for all vectors $V$ and $W$. Accordingly, the spin tensor is defined by the Tulczyjew constraint SP=0.

For given $P$, the general solution of the Tulczyjew constraint is parametrized by the three components of the spin vector ${\mathbf{s}}\in\mathbb{R}^{3}$ that we suppose non vanishing,

[TABLE]

The vector-product is with respect to the Euclidean metric and we define the linear map $j({\mathbf{s}}):{\mathbf{p}}\mapsto{\mathbf{s}}\times{\mathbf{p}}$. We have,

[TABLE]

with the longitudinal spin, or “scalar” spin,

[TABLE]

which turns out to be a constant of the system [16]. The scalar spin $s$ is not to be confused with the norm $\|{\mathbf{s}}\|$ of the spin vector. The helicity or handedness of the photon is $\mathrm{sign}(s)$.

In the Schwarzschild metric we obtain with the notations (2.3),

[TABLE]

The following vector identity will be useful:

[TABLE]

2.4 Conservation laws

The group of isometries of Schwarzschild spacetime is $\mathrm{O}(3)\times\mathbb{R}$, its generators are the Killing vector fields of the metric (2.7), $Z=\varepsilon^{i}_{\,jk}\,\omega^{j}x^{k}\,\partial/\partial x^{i}+\epsilon\,\partial/\partial t$, where $\bm{\omega}\in\mathbb{R}^{3}$ and $\epsilon\in\mathbb{R}$ stand for infinitesimal rotations and time translations, respectively; the $\varepsilon^{i}_{\,jk}$ are the structure constants of $\mathfrak{so}(3)$. Using the general expression [18]

[TABLE]

of the “moment map”, $\Psi$, associated with a Killing vector field, $Z$, together with the expressions (2.13) and (2.14) for $P$ and $S$, we find in a straightforward fashion $\Psi(Z)=-\bm{{\mathcal{L}}}\cdot\bm{\omega}+{\mathcal{E}}\,\epsilon$ where

[TABLE]

is the conserved energy and

[TABLE]

the conserved angular momentum featuring both an extra spin contribution. The latter equation defines an affine map between spin and angular momentum. We will use its inverse:

[TABLE]

2.5 Specifying the Souriau-Saturnini equations

Now that we have defined the metric together with the objects appearing in the equations of motion and the conserved quantities, we are ready to spell out the Souriau-Saturnini equations (2.4 - 2.6) for the case of Schwarzschild spacetime.

Let us introduce the shorthand,

[TABLE]

To obtain the equations of motion in $3$-space, we trade the affine parameter $\tau$ for the coordinate time $t$ using (2.4),

[TABLE]

which we assume non-vanishing. By abuse of notation we write $\tau(t)$ for the inverse function of $t(\tau)$ and we do not distinguish ${\mathbf{x}}={\mathbf{x}}(t)={\mathbf{x}}(\tau(t))$ and likewise for ${\mathbf{p}}$ and ${\mathbf{s}}$. Then we have, from the Souriau-Saturnini equations (2.4 - 2.6) and (2.28):

[TABLE]

With the equations above, we can verify that the conserved quantities, namely the scalar spin (2.16), the energy (2.24) and the total angular momentum (2.25) are conserved. We have indeed $d{\mathcal{E}}/dt=d\bm{{\mathcal{L}}}/dt=ds/dt=0$.

We can simplify the system by only considering the equations of position and momentum (2.29, 2.30) and by eliminating $[{\mathbf{x}}\times{\mathbf{p}}\cdot{\mathbf{s}}]$ and $({\mathbf{s}}\cdot{\mathbf{x}})$ in favour of the conserved angular momentum $\bm{{\mathcal{L}}}$ using equation (2.26) and by eliminating $({\mathbf{s}}\cdot{\mathbf{p}})$ in favour of the conserved scalar spin $s$ using equation (2.16). We use the following relations

[TABLE]

We are thus left with six equations for six unknown functions of $t$, which will be spelled out later, (5.65, 5.66).

We also have a formula for the norm of ${\mathbf{p}}$ from the conserved quantities (2.24) and (2.25),

[TABLE]

and

[TABLE]

Noticing that this last equation and the three equations for position only depend on ${\mathbf{p}}/\|{\mathbf{p}}\|$ our system effectively reduces to five equations.

The results above can already by found in Saturnini’s thesis [17] of 1976.

2.6 Radial case

The first observation is that in the radial case, i.e. with an initial momentum parallel to the initial position, the equations of motion (2.29 - 2.31) reduce to those of the radial geodesics,

[TABLE]

While the differential equation (2.38) displays the well known redshift effect of light, it is striking that we have the same expression (2.39) for the evolution of the transverse spin. This can be expected when looking at the Souriau–Saturnini equations (2.4 - 2.6) and noticing that the redshift terms in (2.38) and (2.39) come from the covariant derivative. Indeed, when the photon is following the geodesic trajectory, the Souriau-Saturnini equations reduce to the geodesic equations i.e. $\dot{X}=P$ and $\dot{P}=\dot{S}=0$, meaning that both $P$ and $S$ are parallel transported.

We also take the opportunity to note that equation (2.24) tells us that the conserved energy ${\mathcal{E}}$ is modified by the transverse spin in general, but not in the radial case.

3 Null geodesics & spinless gravitational lensing

In this section, we first show that we recover the known spinless, massless case, albeit in a slightly different form than the usual geodesic equations, by putting ${\mathbf{s}}=0$ in our equations, and we rederive the well known deviation angle $\Delta\varphi$.

3.1 Some preliminaries

If we put ${\mathbf{s}}=0$ in (2.24) and (2.25), the Noether quantities are of the form

[TABLE]

From $dX/d\tau=P$, and equation (2.13), we find

[TABLE]

For null geodesics, $P^{2}=0$, we have

[TABLE]

Taking advantage of the conservation of total angular momentum, $\bm{{\mathcal{L}}}$, we compute ${\mathbf{x}}\times\bm{{\mathcal{L}}}$ and end up with

[TABLE]

where the function $\lambda$ satisfies (using (3.41), (3.43) and (3.42))

[TABLE]

with ${\mathcal{L}}=\|\bm{{\mathcal{L}}}\|$. We note that $\lambda^{2}\geq 0$ implies a condition on ${\mathcal{E}},{\mathcal{L}}$ and $r$. By taking the scalar product on both sides of equation (3.43) with ${\mathbf{x}}$, we obtain the simple expression

[TABLE]

We record for further use that

[TABLE]

Let us stress that the latter equations lead precisely to the equations of null geodesics given in terms of the Christoffel symbols (2.10). Here we used, instead, the conservation laws, including a number of computational tricks, to obtain the velocity (3.43). Note that the time-component of the geodesic equation yields (up to a global sign):

[TABLE]

which is clearly non-vanishing. Comparison with the general equation (2.28), which is ill-defined in the limit $s\to 0$, shows a striking similarity with equation (3.47), namely the latter is identical to the former provided we ignore the spin-dependent factor on the RHS.

To make the link with equations (2.29) and (2.30), let us write down the equations of motion of the null geodesic in the form:

[TABLE]

3.2 Lensing in weak fields

We restrict our analysis to geodesics remaining in regions of space where the gravitational field is weak, i.e. where all distances $r(t)$ remain much larger than the Schwarzschild radius $a$,

[TABLE]

and linearize with respect to $\alpha$. We take our initial conditions at $\tau=t=0$:

[TABLE]

To alleviate notations we will write from now on ${\mathbf{x}}=(x_{1},\,x_{2},\,x_{3})$ with lower indices. Following tradition we consider the photon in the $x_{1}$-$x_{2}$ plane with energy $p_{0}>0$ coming in from the left, $x_{0}>0$, with positive impact parameter $b$. We suppose $a\ll b$. Then we have to first order in $\alpha$:

[TABLE]

The equations of motion (3.43) become:

[TABLE]

implying

[TABLE]

Equation (3.45), $\dot{r}=\lambda r$, tells us that the distance of closest approach $r_{p}$ (‘perihelion’) is reached when $\lambda$ vanishes. Therefore

[TABLE]

and $r_{p}\sim b-4\,a+4\,ab/r_{0}$.

Our aim is to compute the scattering angle $\Delta\varphi$ for $x_{0}\rightarrow\infty$. As we have set the cosmological constant to zero, spacetime is flat far away from the mass and there coordinate and physical angles coincide. Denoting by $\varphi_{p}$ the angle of closest approach, we have $\Delta\varphi=\pi-2\varphi_{p}$. We can compute $\varphi_{p}$ by integrating

[TABLE]

between $r_{0}=\infty$ and $r_{p}$. In this interval both $r$ and $\varphi$ decrease and we must choose the positive signs in equation (3.58). Our initial angle is $\varphi_{0}=\pi$ and we obtain with $u\vcentcolon=r/r_{p}$,

[TABLE]

Note the integrable singularity at the perihelion, $u=1$. Finally, in linear approximation, the scattering angle takes its famous value: $\Delta\varphi\sim 4\,GM/r_{p}$.

We thus recover the known geodesic equations in the Schwarzschild metric, and the well known deflecting angle $\Delta\varphi$, from the Souriau-Saturnini formalism and putting ${\mathbf{s}}=0$. The resulting equations of motions (3.48 - 3.49) are first order equations, but are strictly equivalent to the second order geodesic equations. Now, the next step is to consider the spinning case, ${\mathbf{s}}\neq 0$, hence considering the full equations of motions (2.29 - 2.31). This is done in the following sections.

4 Numerical solutions

Since solving the system of equations (2.29, 2.30) is not straightforward, we will use the help of numerical integration to propagate specific initial conditions. These numerical solutions will guide us towards perturbative ones.

The numerical integration meets the usual problem of accuracy errors when computing the difference of two almost identical numbers. It becomes relevant here because the present system of equations involves such computations, especially when conserved quantities are involved, e.g. (2.32). This is why it is better to numerically solve all of the 9 differential equations (2.29 - 2.31), including those of the spin.

Even with such measures, integrating these equations over a long time can be tricky with Mathematica. The step algorithm seems overly cautious and is eager to stop the integration process due to stiffness problems, even though all quantities involved are well defined, finite, and smoothly evolving. We need to select the right precision parameters to keep the step algorithm from stopping the integration. Yet, this does not create instabilities in the trajectory of the simulation and we obtain very precise results.

It is convenient to take the initial conditions not at infinity but at perihelion $r_{0}=r_{p}$ of the trajectory of the photon around the star located at the origin:

[TABLE]

Note that the first component of the initial transverse spin ${\mathbf{s}}^{\perp}_{0}$ vanishes, because at perihelion ${d{\mathbf{x}}}/{dt}|_{0}\cdot{\mathbf{x}}_{0}=0$.

We use SI units here. The photon starts with a wavelength of $\lambda_{0}=600\,$nm and a helicity of $\chi=+1$, the star has a Schwarzschild radius of $a=3\cdot 10^{3}\,$m, and the initial distance from the center of the star to the perihelion is $r_{0}=3\cdot 10^{5}\,$m. The numerical integration runs from [math] to $0.1\,$s. While we have $s=\hbar$ in the initial conditions (4.60), we will put $s_{0}^{\perp}=0$ for the time being, because otherwise the trajectory leaves the neighbourhood of the geodesic when $s_{0}^{\perp}$ is close to $\hbar$. We will come back to the transverse spin in the perturbative analysis in the next section.

Figure 1 shows the trajectory of the spinning photon in the geodesic plane. This trajectory is almost identical to the null geodesic one. Indeed, the difference between the coordinates $x_{1}$ and $x_{2}$ of spinning and spinless photons is of the order of the nanometer at the end of the numerical integration. The main differences are the transverse components $x_{3}$ of the trajectory, and $p_{3}$ of the momentum, pictured in Figures 2 and 3 respectively. While the geodesic trajectory is contained within the plane $(x_{1},x_{2})$, the equations of motion (2.29 - 2.31) imply non-vanishing transverse components $x_{3}$ and $p_{3}$.

The angle $\beta$ of the trajectory going out of the plane is small, but constant. As shown in figure 2, it is about $\beta=-6.3\cdot 10^{-8}\,{}^{\prime\prime}$. The sign of the angle $\beta$ depends directly on the helicity $\chi$. Indeed when changing the helicity from $+1$ to $-1$, the amplitude of the angle stays the same, but the sign switches. We see from numerical integrations that the trajectories of two different helicity photons are symmetric with respect to the null geodesic. The transverse momentum $p_{3}$ also shows the same behavior under helicity changes and its sign is again opposite to that of $x_{3}$.

In the next section, we will confirm and explain these results with a perturbative approach.

5 Perturbative solutions

We wish to compare the behavior of our system (2.29, 2.30) describing the trajectories of photons with their due spin to the behavior of null geodesics.

Now, define two constant small parameters,

[TABLE]

where, for the sun, $\alpha$ is typically of the order of $10^{-6}$ and $\epsilon$ of the order of $10^{-16}$ for photons in the visible spectrum. A small $\epsilon$ corresponds to photons having a wavelength much smaller than its distance to the star, which is a sensible hypothesis. Due to the particularities of this system of equations, namely $D$ (2.27) being of order $\epsilon$, we must consider second order terms in $\epsilon$ to obtain the first order equations. In $\alpha$, linear terms will be sufficient.

Let us redefine the spin by setting

[TABLE]

where $\chi=\pm 1$ is the helicity of the photon and $w$ is finite and dimensionless. We easily obtain the conserved quantities (2.24) and (2.25) from the initial conditions (4.60),

[TABLE]

We define the normalized quantities,

[TABLE]

We can then write the equations (2.29 - 2.30) as

[TABLE]

Let us momentarily forget the physical aspect of this system and set $\epsilon=0$ in (5.63). Then with the initial conditions (4.60) the differential equations (5.65) and (5.66) reduce to those of the null geodesics (3.48) and (3.49). Indeed, from (2.32), we have initially $\widetilde{\mathrm{xps}}|_{0}=0$ and $\widetilde{\mathrm{x{\mathcal{L}}}}|_{0}=0$, reducing the initial system to the geodesic one. If we are on a geodesic trajectory, which is in the plane spanned by ${\mathbf{x}}_{0}$ and ${\mathbf{p}}_{0}$, then $\widetilde{\mathrm{x{\mathcal{L}}}}=0$ and $\widetilde{\mathrm{xps}}=0$ continue to vanish due to geodesic conservation of angular momentum and the photon continues on the geodesic trajectory.

This heuristic argument and our numerical results in the last section motivate the ansatz

[TABLE]

where $x_{1},\,x_{2},\,p_{1},\,p_{2}$ solve the geodesic equations. Define $r_{g}=\sqrt{x_{1}^{2}+x_{2}^{2}}$ and similarly for $p_{g}$. To leading order, we have

[TABLE]

In order to recover the geodesics in the limit $\epsilon\rightarrow 0$, we thus need these two leading terms to be zero implying the initial transverse spin to vanish and some conditions on first order terms in $\epsilon$ that are valid at least to first order in $\alpha$:

[TABLE]

Plugging the ansatz (5.67) into the six scalar equations (5.65) and (5.66) we obtain twelve equations: six in $\epsilon^{0}$ and six in $\epsilon^{1}$. The six equations in $\epsilon^{0}$ are equivalent to the four equations (5.70 - 5.73). The six equations in $\epsilon^{1}$ yield:

[TABLE]

At this point, we may even obtain the terms of order $\alpha\epsilon$ giving us constraints on $z_{1}$ and on the initial transverse spin and we end up with $y_{1}\sim y_{2}\sim q_{1}\sim q_{2}\sim{\mathcal{O}}(\alpha^{2})$ and

[TABLE]

and our perturbative solution reads

[TABLE]

Finally, using (2.26) and ${\mathbf{s}}=\frac{{\mathbf{p}}}{\|{\mathbf{p}}\|}s+{\mathbf{s}}^{\perp}$, we obtain the perturbative solution for the transverse spin,

[TABLE]

The most striking effect of the spin on the trajectory of the photon is that it leaves the geodesic plane, but its projection on this plane coincides up to order $\epsilon\alpha$ with the geodesic. The angle $\beta$ between the trajectory and the geodesic plane is obtained from $\beta\sim d(\epsilon y_{3})/dx_{2}$ at infinity, which is immediate with the help of (5.71),

[TABLE]

with the definition (5.61) for $\epsilon$ and where $\lambda_{0}$ is the wavelength of the photon at perihelion. Notice that this angle depends both on the helicity of the photon $\chi=\pm 1$ and on its wavelength. Photons of the two different helicities follow symmetric trajectories with respect to the geodesic and the dependence on $\lambda_{0}$ produces a rainbow effect. In the case of the sun, with $r_{0}$ its radius, this means that two photons starting at the perihelion with opposite helicity will have an offset given by $2\beta=5.7\cdot 10^{-11}\,{}^{\prime\prime}$. If these two photons then travel to the Earth, the offset between them would be of the order of $41\mathrm{\mu m}$ in perfect conditions. The angle $\beta$ has the curious property of being independent of the mass of the star, at lowest order in $\alpha$. This seems to imply that this angle does not vanish as the mass of the star becomes arbitrarily small. Let us note though, that the limit $\alpha\rightarrow 0$ is ill defined in the equations of motion and therefore in the perturbative solution. Indeed, the first of the Souriau-Saturnini equations (2.4) is independent of $a$ because both $R(S)(S)$ and $SR(S)P$ are proportional to $a$. The introduction of a cosmological constant will regularize this singularity, even at small scales, as will be shown in a forthcoming work.

Also, there is no correction of order $\epsilon\,\alpha$ to the usual deviation angle $\Delta\varphi$ in the plane, computed in subsection 3.2.

Note that the transverse component of the momentum quickly reaches its maximum at a distance of a few $r_{0}$, which is $\epsilon\,{q_{3}}_{max}=2\,\epsilon\,\alpha\,\chi\,p_{0}$. Since the angle $\beta$ comes from a spin-orbit-like effect of the star on the trajectory, we would expect it to only act close to the star. To avoid this problem, we define $\gamma$ to be the angle between the geodesic plane and the momentum carried by the spinning photon. We have:

[TABLE]

This angle does depend on the mass of the star and is even smaller than $\beta$. For the sun we have $2\gamma=4.9\cdot 10^{-16}\,{}^{\prime\prime}$.

Our perturbative results for $y_{3}$ and $q_{3}$ above match our numerical results with a relative error of about $10^{-9}$ and $10^{-4.5}$, respectively. The match is better for $y_{3}$ because it contains terms of order 1 and of order $\alpha$, while $q_{3}$ is of order $\alpha$.

The analysis by Gosselin, Bérard & Mohrbach [34] starting from the Bargmann-Wigner equations shows birefringence with an angle equal to our $\gamma$ but in the geodesic plane.

6 Conclusions

For photons, quantum mechanics teaches us that the longitudinal component $s$ of the spin is $\pm\hbar$. This is in harmony with the conservation of $s$, which follows in general from the Souriau-Saturnini equations. Quantum mechanics also teaches us that the norm of the transverse spin $\|{\mathbf{s^{\bot}}}\|$ is $\hbar$. Two remarks arise from the present work. First, we saw in the radial case that the photon follows the null geodesic trajectory, and that the transverse spin undergoes the same evolution as the momentum: it is parallel transported. However, in our non-radial perturbative solution, equation (5.85), this norm vanishes at perihelion and then grows linearly with time $t$ (to leading order). The linear growth implies that our perturbation theory breaks down for large times. This instability is absent from a generic Robertson-Walker metric where the norm of the transverse spin is proportional to the inverse Hubble parameter [33].

With its continuously varying transverse spin, the instability reminds us of the instability of the classical hydrogen atom and its continuously varying energy. Indeed, the equations we use here are purely classical. The longitudinal spin is a constant of the system. Its definition comes from the co-adjoint representation of the Poincaré group [18, 36]. The transverse spin is not conserved. Its two degrees of freedom come from the introduction of the dipole moment. However, it is not clear from the geometrical derivation of the equations of motion, if these two degrees of freedom are the transverse spin in quantum mechanics. A way to determine their exact meaning would be to derive the Souriau-Saturnini equations (2.4 - 2.6) from quantum mechanics, à la eikonale. This is currently under investigation.

Notice also that the out-of-plane momentum is in the opposite direction with respect to the offset. This means that the star is intrinsicly acting on the photon’s position and momentum, i.e. a spin-orbit effect. Yet, at large time $t$ in the perturbative solution, we see that the trajectory’s offset keeps increasing linearly, while the momentum stays constant and in the opposite direction. We would expect, once we are sufficiently far away from the star, that the star loses grip on the photon. Since spacetime is flat far away, we expect the photon’s momentum to carry the trajectory, which is not what we find. This is in line with the fact that we don’t recover the equations of motion in flat spacetime in the limit $a\rightarrow 0$. We will see in a future work that including the cosmological constant helps us mitigating this problem, namely recovering null geodesics “far away” from the star, even at scales much smaller than those the cosmological constant would suggest. The cosmological constant also puts an upper bound to the transverse spin, whose value still remains in conflict with quantum mechanics.

For us, the most interesting features of birefringence in the Schwarzschild metric are the out-of-plane contributions to trajectory and momentum. First, we have the linearly growing offset – given by an angle $2\beta$, equation (5.86) – between the trajectories of opposite polarisations. Then, the Souriau-Saturnini equations in the Schwarzschild metric become singular far away from the star, a singularity absent in the Kottler metric. Therefore we expect the offset induced by the angle $\gamma$ (5.87) to play a more important role in observations. (The computations with non-vanishing cosmological constant are complicated and will be presented in a later work.) In any case, both angles, $\beta$ and $\gamma$ are wavelength dependent and the offset must feature a rainbow effect.

Despite the mentioned classical instability, we wonder whether this type of gravitational birefringence is accessible to experimental verification. The upper bound of the 1976 Very Long Baseline Interferometry experiment [35] achieves an upper limit for the birefringence angle $\beta$ (or $\gamma$?) of the order of $10^{-3}\,{}^{\prime\prime}$ for $\lambda\sim 10$ cm. With our formulas, we get for this wavelength $\beta\sim 10^{-6}\,{}^{\prime\prime}$ and $\gamma\sim 10^{-11}\,{}^{\prime\prime}$. After 40 years, the needed accuracy for testing our formulas should not be out of reach.

We started this work as a rope team of three. But in September 2018 we lost our friend and guide, Christian Duval. Any error or short coming remaining in this work is ours.

Acknowledgements: This work has been carried out thanks to the support of the OCEVU Labex (ANR-11-LABX-0060) and the A*MIDEX project (ANR-11-IDEX-0001-02) funded by the "Investissements d’Avenir" French government program managed by the ANR.