Skyfield's Github has discussion Jupiter hiccup #815 which then links back to to Non-physical gravitational deflection corrections for Solar System bodies #734.

The script and plot from #815 are shown below. The half amplitude of the 'hiccup' or wiggle in Jupiter's RA is about 0.0001 radians or 1.7 arcminutes, which is quite a sizable gravitational deflection if it were true.

My understanding of space and time is about like this (search for answer by "uhoh") so I can't do any of this myself, but it looks like during this event (between 05:00 and 10:00 on January 1st, 2032) Jupiter passes directly behind the Sun as seen from Earth, and the Skyfield calculation is grappling with a path that passes directly through the Sun. JPL Horizons shows the Sun-Observer-Target (SOT) angle drops to 0.077 degrees during the middle of this pass!

Because of how Skyfield works (simple, easy to use Python interface) and is regularly used to find events through minimizing routines, any observation calculation needs to always return reasonable and continuous position values with respect to time, so no NaNs or sudden jumps.

What I'd like to ask in this question is for an expression for the gravitational deflection (and ideally a light-time) for a ray of light that grazes or even passes through the Sun.

Yes, EM radiation above the Sun's plasma frequency can't go through the Sun,

but since Skyfield will still need something continuous, there still may be an expression that properly handles the mass distribution of the finite-sized Sun rather than lumping it to a single point.

Could one just use the Sun's density profile and integrate a "gravitational tugging" of the trajectory from this mass distribution a bit like models from gravitational lensing by a galaxy do?

Is there a better way, or preexisting, peer-reviewed and vetted method?

Question: Light or (massless) neutrinos graze or pass through the Sun and arrive at Earth - need an expression for Sun's gravitational effect on observation direction/time.

Skyfield calculation RA of Jupiter between 05:00 and 10:00 January 1, 2032 showing a big deflection by the Sun's gravity as it passes directly behind it

from skyfield.api import load
import matplotlib.pyplot as plt

ts = load.timescale()
planets = load('de421.bsp')

earth  = planets['earth']
jupiter = planets['jupiter barycenter']

t0 = ts.ut1(2032, 1, 1, 5, 0, 0)
t1 = ts.ut1(2032, 1, 1, 10, 0, 0)
t = ts.linspace(t0, t1, 301)

ra =  earth.at(t).observe(jupiter).apparent().radec('date')[0]
plt.plot(t.ut1, ra.radians, '-')
  • 1
    $\begingroup$ Yes, I know. Actual GR calculations are above my pay scale. I'm just providing context. $\endgroup$
    – James K
    Commented Dec 4, 2022 at 9:10
  • 1
    $\begingroup$ Try pyAutoLens library $\endgroup$
    – user47732
    Commented Dec 9, 2022 at 15:28
  • 1
    $\begingroup$ Hi @uhoh, Just a quick question, Why is it required for The lights frequency < Plasma resonant frequency, I though The Lights frequency > Plasma resonant frequency to be transparent $\endgroup$
    – user47732
    Commented Dec 10, 2022 at 12:17
  • 2
    $\begingroup$ Wouldn't know what metric to use inside the Sun. However, it will certainly rely on knowing the density profile of the Sun, so there will be no simple analytic expression like for the Schwarzschild metric. $\endgroup$
    – ProfRob
    Commented Jun 2, 2023 at 6:28
  • 3
    $\begingroup$ @PM2Ring well that will do for massless neutrinos, since the pressure term inside the Sun will be quite negligible. Now all you need to do is compute the Shapiro delay and generalise the deflection and lag for neutrinos with mass. $\endgroup$
    – ProfRob
    Commented Jun 2, 2023 at 15:34

2 Answers 2


My solution implements PyAutoLens library which calculates the The Gravitational lensing effect.

My code:

# Hi
#Neccerary Modules
import autolens as al
import autolens.plot as aplt
# Generating grid
grid_2d = al.Grid2D.uniform(shape_native=(50, 50), pixel_scales=0.05)
grid_2d_plotter = aplt.Grid2DPlotter(grid=grid_2d)
# Generating Sersic Profile
sersic_light_profile = al.lp.Sersic(
    centre=(0.0, 0.5),
    ell_comps=(0.1, 0.10),
# Displaying it

image_2d = sersic_light_profile.image_2d_from(grid=grid_2d)
light_profile_plotter = aplt.LightProfilePlotter(light_profile=sersic_light_profile, grid=grid_2d)

#Note: I haven't enter the Suns value

The Image created is the Sersic light profile

Image generated from the above code:

enter image description here

  • 1
    $\begingroup$ Cool! Having a Python package handy to do this calculation is wonderful. I see that the package is well documented, so it can be left as an exercise to plug in the specific values for the Sun's density profile. Thanks! $\endgroup$
    – uhoh
    Commented Dec 10, 2022 at 13:49
  • $\begingroup$ Can you explain what it is your code does, and what your picture represents? $\endgroup$
    – ProfRob
    Commented Dec 21, 2022 at 20:30
  • $\begingroup$ @ProfRob My code shows the Sersic light profile (which is a generalization of De Vaucouleur's law meaning it shows how the Intensity/I of the Galaxy varies with it's Radius/R). The Code is fitted with random values. The image represents the same. The code isn't yet completed it requires the mass profile. Inorder to change the galaxy info You just have to change the ell_comps, sersic_index, centre, effective_Radius and to finally Integrate the mass profile $\endgroup$
    – user47732
    Commented Dec 22, 2022 at 5:31
  • $\begingroup$ How does that address the question? $\endgroup$
    – ProfRob
    Commented Dec 22, 2022 at 7:06
  • $\begingroup$ Well, it simulates the light profile being emitted from the galaxy to be GR (Gravitationally) Lensed $\endgroup$
    – user47732
    Commented Dec 22, 2022 at 11:33

The usual equation for the deflection of light (or a massless particle), used for example for calculations of strong lensing through a cluster is derived in Post-Newtonian theory, which is appropriate for regions in a weak field, i.e., gravity is not too strong. The criteria for a spherical source to be only weakly gravitating is that $\frac{M(r)}{r} << \frac{c^2}{G}$. We know that the event horizon for a solar mass $\frac{2GM_\odot}{c^2}$ = 3 km and the radius of the Sun is $7 \times 10^{5}$ km, so M/R$_\odot$ at the surface of the Sun is $1.5/(7 \times 10^5)$ of $\frac{c^2}{G}$. That means the surface of the sun is well into the weak field limit.

What about the interior of the Sun? The Wikipedia page on Stellar Structure shows the mass distribution for the Sun. The interior values of $\frac{M(r)/M_\odot}{r/R_\odot}$ never rises above about 2 times its value at the surface. We therefore know that the whole interior of the Sun is in the weak limit.

Another issue might be deciding what to integrate to determine the full GR relevant mass. In GR, the total gravitating mass is the sum of the rest mass + internal energy + gravitational potential energy (see The General-Relativistic Theory of Stellar Structure by K. S Thorne ). But, since we are far in the weak limit, the two energy terms, which both depend on M/R (internal energy is set by Virial Theorem), can be ignored. And we again have M(r) depends on rest mass.

The Post-Newtonian calculated gravitational deflection angle is given by: $$ \alpha = \frac{D_{LS}}{D_S}\frac{4GM}{c^2b^2} $$

  • $D_{LS}$ = Distance from lens to source
  • $D_{S}$ = Distance to source
  • b = impact parameter = $D_L(\beta+\alpha) \sim D_L\beta$
  • $D_L$ = Distance to lens
  • $\beta$ = angle to source
  • $M$ is mass in cylinder of radius b You can scale things by knowing that the deflection angle at the limb of the sun is 1.752 arcsec.

There may be multiple solutions, so one may need to search a large area in the lens plane to find all of the solutions where the rays go from the source to viewer.

  • 1
    $\begingroup$ Are you certain that this approximation is valid when $b$ is smaller than the diameter of the Sun? Can you support that with authoritative sources? This looks like it's just the first order approximation valid only for paths which remain outside the object of mass $M$. $\endgroup$
    – uhoh
    Commented Jun 7, 2023 at 21:30
  • 1
    $\begingroup$ As uhoh said, that looks like the standard deflection equation, combined with the small angle approximation, $\theta\approx\sin\theta\approx\tan\theta$. Strictly, it's only valid outside the body, but ProfRob says we can ignore the pressure inside the Sun. Also, that deflection approximation is only valid when $b$ is much larger than $r_s$, as my graphs here show. OTOH, that's probably a safe assumption, even for most paths inside the Sun. $\endgroup$
    – PM 2Ring
    Commented Jun 7, 2023 at 22:25
  • 1
    $\begingroup$ FWIW, I have a Sage / Python script here which calculates the deflection from the impact parameter using an elliptic integral. $\endgroup$
    – PM 2Ring
    Commented Jun 7, 2023 at 22:31
  • 1
    $\begingroup$ @uhoh No, because $M$ is the mass within a cylinder of radius $b$, so the $b^2$ in the numerator & denominator cancel. However, I don't understand why eshaya is using a cylinder instead of a sphere. (Of course, with a sphere, we get $b^3$ in the numerator, once again avoiding divergence). $\endgroup$
    – PM 2Ring
    Commented Jun 8, 2023 at 15:12
  • 1
    $\begingroup$ The mass within a cylinder contributes to the lensing. For a ray moving parallel to the z axis, mass near the poles always contributes. $\endgroup$
    – eshaya
    Commented Jun 8, 2023 at 18:13

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .