Opposition surge refers to the significant increase in brightness observed on many objects as their phase angle approaches zero (i.e. at opposition), when the observer is located directly between the object and the sun. Opposition surge can be observed on essentially all solid bodies in the Solar system, but is most prominent on rocky airless bodies and Saturn's rings. The effect is caused by a combination of shadow-hiding (particles on an object's surface blocking the view of their own shadows at zero phase), retroreflection (objects preferentially reflecting light back to the source due to their material properties and geometry) and coherent backscattering (constructive interference of reflected light waves back towards the light source). Shadow-hiding seems to be the main contributor in most cases.

Scroll down for real-life images of opposition surge.

Moon

Opposition surge on the lunar surface, centered on the shadow of Buzz Aldrin
Photo credit: NASA

In this video of an Earthrise from Japan's Kaguya lunar orbiter, opposition surge is visible during the first 18 seconds as a bright moving spot on the left side of the frame
Video credit: JAXA/NHK

Mars

Opposition surge is visible near the center of this low-orbit photo mosaic from the Mars Global Surveyor spacecraft
Photo credit: NASA/JPL/MSSS

Opposition surge is visible just above the center of this colorized image taken by the Perseverance rover during its descent to the martian surface
Photo credit: NASA/JPL-Caltech/Simeon Schmauß/AndreaLuck

Saturn

Saturn's rings display an especially strong opposition surge in this photo from the Cassini spacecraft
Photo credit: NASA/JPL/Space Science Institute

Asteroids and Comets

Opposition surge is clearly visible on this photo of comet 67P from the Rosetta spacecraft
Photo credit: ESA/Rosetta/MPS for OSIRIS Team MPS/UPD/LAM/IAA/SSO/INTA/UPM/DASP/IDA

More Examples

Some great examples of opposition surge on Earth: [1] [2] [3] [4]
A few of the many great images processed by Thomas Appéré showing opposition surge on Mars: [1] [2] [3]
Opposition surge is visible in the descent videos captured by NASA's Curiosity and Perseverance Mars rovers

For a quick introduction to the SpaceEngine climate model and its main features, see our blog post here. In this article we provide a closer look at the physics and calculations behind planetary temperatures.

The temperatures of the SpaceEngine climate model are based on energy transport calculations and account for planetary albedo, presence of an atmosphere, atmosphere properties (including wind speeds, radiation and advection, and greenhouse effects), internal planetary heating, day sides, night sides, planet obliquity (for seasons and varying daylight hours, polar days and polar nights), eccentric orbits, tidal locking, and incident light of all stars in the system. Additionally, the altitude dependence had a major overhaul from its previous implementation. Now, we use real temperature-pressure profile data for different types of atmospheres, allowing for exciting behavior like temperature inversions (like we see at the tropopause, stratopause, and mesopause for Earth).

All planets are born hot, but without a sustained energy source like the fusion taking place inside of stars, they cool quickly. So unless planets are very young, almost all of their energy comes from stellar irradiation. Therefore, the main focus of SpaceEngine’s climate model is on stellar irradiation. We also account for internal planetary heating and greenhouse effects when an atmosphere is present.

The stellar part of the model is implemented in three main parts: 1) a starting temperature, based on the distance to and properties of the host star(s); 2) a longitudinal dependence based on the type of planet; and 3) a latitudinal dependence that also accounts for obliquity (axial tilt).

Where do we start? Calculating a Planetary Temperature

We start with the Planetary Equilibrium Temperature (the theoretical temperature of a planet if it were a blackbody heated only by its parent star):

\begin{equation}
T_{\mathrm{equilibrium}} = \left ( F_*\frac{(1-A)}{\sigma_{\mathrm{SB}} f} \right )^{1/4}= \left ( \frac{L_*}{4\pi d^2}\frac{(1-A)}{\sigma_{\mathrm{SB}} f} \right )^{1/4} ,
\end{equation}

where \(F_*\) and \(L_*\) are the flux and luminosity of the host star, respectively, \(d\) is the distance between the star and planet, \( A \) is the bond albedo of the planet, \( \sigma_{SB} \) is the Stefan-Boltzmann constant, and \( f \) is a redistribution factor. Usually, \( f=4 \) for planets where all incident energy is uniformly distributed, and \( f=2 \) for planets where energy is distributed over the starlit hemisphere only, like a tidally locked planet (see, for example, Selsis et al. 2007 A&A 476).

This value is computed for each star in the system at the current point in the orbit (\(d\)) changes as time progresses for elliptical orbits) and is used to calculate a fiducial temperature (\(T_{\rm fiducial}\)) that a surface temperature map is built out from.

Day/Night Cycles and Longitude

We split planets into three categories: terrestrial planets, gas giant planets, and tidally locked planets (that can be terrestrial or gas giant). Each category has a different base temperature model in SpaceEngine. Our first step in all three categories was to model how the temperature changes around the planet with longitude, i.e.: how should the temperature change from day to night?

Starting with terrestrial planets with atmospheres, we used a thermal transport model to calculate temperatures across planet surfaces to account for day and night. We considered the temperature of a gas parcel as it travels over the planet surface and the planet rotates away from the host star and then back towards it. This gas parcel absorbs incident flux from the host star on the day side and radiates heat away during both day and night. We made the simple assumption that this gas parcel radiates as a blackbody. To determine the temperature at any time during the day or night, we set up and solved a thermodynamic energy equation (e.g.: Showman et al. 2009, Crossfield course notes 2019):

\begin{equation}
\frac{\mathrm{d}T}{\mathrm{d}t} = \frac{1}{c_p \rho H} \Delta F = \frac{1}{c_p \rho H} \left( \sigma_{\mathrm{SB}} {T_{\rm fiducial}}^4 T_{\rm latitude} ~{\rm max}(\cos \theta, 0) - \sigma_{\mathrm{SB}} T^4 \right) ,
\end{equation}

where \( T \) is the temperature we are solving for, \( t \) is time, \( c_p \) is the heat capacity (under constant pressure, since this gas parcel moves along a constant surface on the planet), \( \rho \) is the atmospheric density, \( H \) is the thickness of the atmosphere, and \( \Delta F \) is the net flux on the gas parcel. \( T_{\rm fiducial} \) (in Kelvin) is a function of starting planetary equilibrium temperature, \( T_{\rm latitude} \) is a unitless latitude scaling (described below), \( \theta \) is the longitude measured relative to the subsolar point (where the host star is directly overhead), and "max" is a function that gives the maximum of the two arguments (to account for the fact that when the parcel is on the night side of the planet is absorbs zero flux, not negative flux).

We rewrote this equation and substituted a longitude dependence for the time dependence so solving the equation resulted in a longitudinal temperature map for the planet, where \( \theta \) is a proxy for time of day. We used a radiative time scale (\(\tau_{\rm rad}\)) and an advective time scale (\(\tau_{\rm adv}\)) to convert time to longitude:

\begin{equation}
\tau_{\rm rad} = \frac{c_p \rho H}{\sigma_{\mathrm{SB}} {T_0}^3},
\end{equation}

and

\begin{equation}
\tau_{\rm adv} = \frac{2 \pi R_p}{v_{\rm wind}} ,
\end{equation}

where \(T_0\) is the fiducial temperature combined with the latitude scaling, \(R_p\) is the planet’s radius, and \(v_{\rm wind}\) is a global wind speed.

The ratio of these two timescales tells us how heat is carried around the planet. When \(\tau_{\rm rad} / \tau_{\rm adv} \gg 1\), winds carry the heat around the planet before it is radiated away, the planet stays warmer overall. If \(\tau_{\rm rad} / \tau_{\rm adv} \ll 1\), radiation is more efficient than advection, so energy is radiated away before it is carried around the planet; the planet stays cooler overall, and we get more “averaged” temperatures over the surface (the amplitude between the daytime maximum and nighttime minimum temperatures is smaller). This also produces a phase offset of the maximum temperature in our model. This means the hottest point on a planet’s surface is not necessarily at noon, when the host star is directly overhead (note that on Earth, winds are not the reason for this offset — on Earth it is warmer shortly after noon than directly at noon due to “thermal inertia.” We have not modeled that effect here, but our phase-offset due to winds makes a nice approximation for this in SpaceEngine).

The resulting longitudinal temperature map is composed of two solutions of the thermodynamic energy equation: a daytime function where the gas parcel absorbs incident flux and radiates away heat, and a night time function where the gas parcel only radiates away heat. The longitudinal temperature map of an arbitrary planet with varying wind speeds is shown in the following figure:

A longitudinal temperature map for a random terrestrial planet with an atmosphere. The temperature along the equator of a planet with 0 obliquity (axial tilt) for different global wind speeds is shown. There are three curves for a planet with identical atmospheric and orbital parameters in SpaceEngine and they only differ in global wind speed. \(\theta = 0^\circ, 360^\circ\) is the subsolar point, or noon. \(\theta = 90^\circ\) and \(\theta = 270^\circ\) are the relative longitudes of sunset and sunrise, respectively (for a planet with 0 obliquity), and \(\theta = 180^\circ\) is the relative longitude of midnight, when the host star is on the opposite side of the planet. Higher wind speeds means heat is carried around the planet, before it radiates away. This results in more average temperatures over the surface, and an offset of the maximum temperature from the position where the host-star is directly overhead.

For terrestrial planets without atmospheres, heat is not carried around the planet by moving air. For such planets, we removed the atmospheric component of the model and the resulting longitudinal temperature maps are similar to the “terrestrial with atmospheres” case, but with the hottest point being always at the subsolar point, and without the “averaging” effect over the surface of the planet. This results in a more extreme temperature variation from day to night.

Gas giant planets have massive atmospheres with significant convection, resulting in longitudinal temperature functions that are smoother, with a less sharp boundary at sunrise and sunset points. We adopted a sinusoidal function for the temperatures of gas giants, where the wind speed corresponds to a phase offset between the subsolar point and the hottest point on the surface of the planet (where, as per convention, the “surface” is the level in the atmosphere where the pressure is equal to 1 atm).

A longitudinal temperature map for the northern hemisphere of a gas giant planet on its northern summer solstice (when the host star’s subsolar point is at its highest latitude —see the next section for details on accounting for latitude!). As with the terrestrial temperature map, \(\theta = 0^\circ, 360^\circ\) is the subsolar point, or noon. \(\theta = 90^\circ\) and \(\theta = 270^\circ\) are the relative longitudes of sunset and sunrise, respectively, at the equator, and \(\theta = 180^\circ\) is the relative longitude of midnight, when the host star is on the opposite side of the planet. Following convention, the “surface” of a gas giant planet is the level in the atmosphere where the pressure is equal to 1 atm.

A longitudinal temperature map for a tidally locked planet. As with the other temperature maps, \(\theta = 0^\circ, 360^\circ\) is the subsolar point, or noon and \(\theta = 180^\circ \) is the relative longitude of midnight, where the host star is on the opposite side of the planet. For tidally locked planets, \(\theta = 90^\circ\) and \(\theta = 270^\circ\) are the relative longitudes of sunset and sunrise, respectively, for all latitudes.

Since tidally locked planets, as the name suggests, always have the same side facing the host star, our longitudinal temperature map for tidally locked planets is based on the Lambert cosine law. This law states that the relative intensity of light on a surface is proportional to the cosine of the angle of incidence. This results in a hot spot at the subsolar point, with temperatures dropping as angular distance from the subsolar point increases, with a fairly constant temperature across the (permanent) night side of the planet. Our model also incorporates a phase offset calculated from the planet’s wind speed. For tidally locked planets, the temperature, on the night side of the planet, is typically 1/3 the day time high (e.g.: Lunine & Lorenz 2002).

For a tidally locked planet on an eccentric orbit, the exact position of the subsolar point will oscillate over the planet’s surface with time. The subsolar point in SpaceEngine is not fixed to a particular latitude/longitude position on the planet’s surface, and this oscillation is included in our model.

Obliquity, Seasons, and Latitude

Our model includes a latitudinal dependence that fixes the maximum temperature at the subsolar latitude (for a planet with zero obliquity/no axial tilt, this is the equator) and a minimum temperature at the maximum relative latitude from the subsolar point (for a planet with zero obliquity/no axial tilt, this would be the poles).

We started with Lambert’s cosine law, which states that the relative intensity of light on a surface is proportional to the cosine of the angle of incidence. The angle in our case is the latitude (\(\phi\)) on the planet relative to the host star’s subsolar latitude. So, if a planet has zero obliquity, at the equator (\(\phi = 0^\circ\)) where the surface normal is parallel to the incident light, this is \(\cos 0^\circ = 1\), and at the poles (\( \phi = \pm 90^\circ \)) we have \(\cos \pm 90^\circ = 0\).

We multiplied our latitudinal dependence with our longitudinal dependence to scale the temperature appropriately, but wanted to modify this slightly. We didn’t want the temperature to go to 0 at the poles, and we wanted to be able to account for obliquities different from zero.

The first modification we made was to re-scale the cosine so we could accommodate obliquities that are not zero, and for latitude we used the absolute value of the latitude relative to the subsolar latitude. This means that on the solstices, when the subsolar latitude is equal to the obliquity (for Earth on the June solstice this would be \(\phi = 23.44^\circ \)), the winter pole is at a relative latitude of \(90^\circ + obliquity\) and the summer pole is at \(90^\circ - obliquity\). So our latitude dependence needed to go to \(90^\circ + obliquity\), not just \(90^\circ\) . Our next modification was to prevent the minimum value from reaching zero. We defined a unitless minimum temperature fraction for the winter pole (\(T_{\rm pole}\)) to scale our latitude dependence for temperature. This minimum temperature is based on the relationship between surface pressure, rotation period, and equator-to-pole temperature difference (e.g.: Kaspi & Showman 2015 ApJ 804), and is calculated for each planet/moon.

We applied this scaling to our longitudinal temperature map to account for the latitude in our temperature calculations.

Left panel: the cosine of \(\phi\) goes to 0 at \(\phi = 90^\circ\). Right panel: Our latitude dependence as a function of latitude relative to the subsolar latitude for a planet with obliquity = \(23.44^\circ\). This function goes to a minimum value of \(T_{pole}\) at \(\phi = 90^\circ + 23.44^\circ\).

The next aspect related to latitude we accounted for is the varying daylight hours throughout the year for planets with non-zero obliquity/axial tilt. For example, on Earth’s June solstice, we know that the northern hemisphere experiences their longest day of the year, and the southern hemisphere experiences their shortest day of the year. At this time, the sun never sets inside of the northern polar circle (polar day), and the sun never rises inside of the southern polar circle (polar night).

The amount of time a point on the planet surface experiences incident sunlight affects the temperature, and is needed to set the boundaries between our day and night longitudinal temperature maps. To include this in our climate model, we needed to compute the fraction of the day a surface point is in sunlight at a given latitude. This is a straightforward calculation that can be done a few different ways, typically with the true anomaly (an orbital parameter). It can also be computed using the subsolar latitude:

\begin{equation}
fraction~of~day~receiving~sunlight = \frac{\arccos \left ( -\tan \phi \tan \phi_{SS} \right )}{\pi} ,
\end{equation}

where \( \phi \) is the latitude of interest and \( \phi_{SS} \) is the subsolar latitude (e.g.: Rauscher 2017 ApJ 846).

The number of hours a given latitude experiences incident sunlight on Earth. Ignoring atmospheric effects and the angular size of the sun, at the equator there are equal hours of daylight and night on every day of the year. Other latitudes experience different hours of daylight throughout the year. On the vernal and autumnal equinoxes, all points on the Earth’s surface experience equal hours of daylight and night (ignoring the angular size of the Sun).

A Complete Surface Temperature Map

By combining the fiducial temperature, the longitudinal and latitudinal dependencies, and the varying length of day, we obtained a complete surface temperature map for a planet. Our model accounts for the properties of the host star, the planet’s orbit, obliquity/axial tilt (resulting in seasons!), rotation, and atmospheric properties (if an atmosphere is present). At this point, we also accounted for any internal heating from the planet itself, and greenhouse effects in the atmosphere.

We added a \( 20^\circ \) obliquity to the planet from the terrestrial longitudinal temperature map figure and applying the latitude dependence for a few sample latitudes, we obtain a more complete picture of the surface temperatures on a planet. The time of local sunrises and sunsets are marked for each latitude curve. For this planet, the atmospheric properties result in a small offset between noon and the hottest temperature of the day.

The surface temperatures for a planet with a \( 20^\circ \) obliquity on the northern summer solstice in SpaceEngine, as determined from our model, plotted on a sphere.

Tatooine and Beyond: Systems with Multiple Stars

The SpaceEngine climate model also accounts for systems with multiple stars! Whether your favorite planet is orbiting a star orbiting a star (an S-type orbit), orbiting two stars (a P-type orbit), or some other configuration with even more stars, it is accounted for in the model. Temperatures are computed as time progresses and as stars and planets move along their orbits. Binary planets and moons (and moons around moons!) are properly accounted for as well.

We treated each star in the system independently, calculating a longitudinal temperature map of the planet for each star. We then combined the effects of all stars in the system. Energies are additive, and energy is proportional to temperature to the power of 4, so our final temperature at a point on the surface will be:

\begin{equation}
T_{\mathrm{final,~current~position}} = [~( T_{\mathrm{current~position,~star~1}})^4 + ( T_{\mathrm{current~position,~star~2}})^4 + ( T_{\mathrm{current~position,~star~3}})^4 + \dots~]^{1/4}.
\end{equation}

Any internal heating from the planet itself, and greenhouse effects in the atmosphere are only added to the model after all stars are accounted for.

Accounting for Altitude and Layered Atmospheres

The final part of calculating temperature was to account for the altitude. Once the surface temperature calculation included the incident light from all stars in the system, the longitude and latitude, the internal planetary heating, and any greenhouse effects due to an atmosphere, we calculated the temperature as a function of altitude.

The main challenge in accounting for altitude is that temperature depends on both altitude and pressure, but pressure also depends on altitude and temperature! In addition, vertical temperature profiles have varying slopes and inversions with altitude, because pressure and atmospheric chemical composition change with altitude.

To avoid the pressure-temperature mutual-dependence, in the SpaceEngine climate model we assumed pressure decreases exponentially with altitude (e.g.: Wikipedia - Scale Height).

Note: Following convention, zero altitude on gas giant planets is set to pressure=1 atm in SpaceEngine.

Then, to allow for temperature inversions and to account for different layers of atmospheres, we used real and simulated pressure-temperature profiles to determine the temperature as a function of pressure. The majority of the pressure-temperature profiles we use in SpaceEngine were computed using the NASA Planetary Spectrum Generator (Villanueva et al. 2018). Other profiles and simulated profiles were sourced from assorted publications (Thorngren et al. 2019, Piette & Madhusudhan et al. 2020, Zhang 2020, Ohno & Fortney 2023).

The suite of pressure-temperature profiles included in SpaceEngine. Surfaces of terrestrial solar system planets are marked with dots. When these profiles are selected for planets with higher surface pressures, or differing surface temperatures, the profiles are scaled appropriately in SpaceEngine.

For each planet, a suitable profile is selected, scaled for the surface pressure and temperature (calculated as described above), and then interpolated at the user’s current altitude/pressure position to determine the local temperature.

And that’s the final step in the temperature component of the new SpaceEngine Climate Model! We hope you enjoyed this peek inside the model.

For a look at how these quantities are displayed in-game, be sure to check out our related blog post here.

Note: None of the figures shown in this article are displayed "in-program": this is simply a look “under the hood” of the SpaceEngine climate model.

When thinking about renders of things like warp drives and black holes we usually just expect to see a simple approximation or an artist rendition, assuming the math required to pull off something accurate would require someone with at least a PhD in Mathematical Physics. Which I won't tell that its completely untrue, but in this blog post I'll try to explain a way to do quite accurate visualizations within a 100 or so lines of code, for basically any kind of space-time for which you can write its metric as code. Although, the detailed mathematical derivation of this approach might be somewhat math heavy.

The main ingredient of any GR render is figuring out how the rays of light move around. Knowing how light moves we can trace rays from the camera into the scene and see where the light came from. So, to render the simplest scene without objects we simply trace a ray for each pixel and assign the color of the pixel to the color of the skybox in the direction in which the ray ends up pointing to.

• What are geodesics?
• Mathematical description of shortest path
• Lagrangian description of a geodesic
• Hamiltonian description of a geodesic
• Writing a Hamiltonian geodesic tracer in GLSL
• Conclusions
• References

What are geodesics?

How exactly do we trace rays in curves space? Any object inside a curved space follows something called a geodesic.

A geodesic is just a fancy word for, in some sense, a path of shortest length between 2 points inside a space.

I should note that there could be multiple of such paths, which are locally minimal, in the sense that you can't nudge the path to make it shorter, while globally there might be a shorter path. Also, in Minkowski space-time the definition is a bit more complicated, because of the imaginary distances (when \( ds^2 < 0 \) ).

In our case, instead of paths between 2 points, we are only interested in finding how a ray moves given an initial point and direction, but the definition above will still prove useful when deriving the equations describing a geodesic, which we will use here.

Mathematical description of shortest path

I'll try to quickly go through the derivation. If you wish to skip the math part, jump to Writing a Hamiltonian geodesic tracer in GLSL.

Mathematically speaking we have some coordinate system, a path, and a way to compute distances between 2 points.

A coordinate system being a set of several numbers labeling each point in the space. A path is a function that takes in the path parameter and outputs a coordinate. In General Relativity the path parameter is usually proper time (like a clock moving with the object, labeling each point), but it can be anything really. The way to compute distances is called a metric (and it's the main source of scary math here).

In physics, or more generally differential geometry, a metric is defined as an integral ("sum") of something called the metric tensor. A metric tensor is a bilinear form \( g(a, b) \), it essentially maps pairs of vectors to real numbers, and is a generalization of dot product for curved spaces. So using a metric tensor we can find the length of a vector in space, and also distances \( ds \) between infinitely close points in space.

\begin{equation}
ds^2 = g(dx, dx)
\end{equation}

In our case, where we describe vectors as a set of numbers, a metric is simply a matrix product of some matrix \( g_{\mu \nu} \) times the vectors. For our infinitesimal distance \( ds \) we get this expression:

\begin{equation}
ds^2 = \sum_{\mu \nu}^N g_{\mu \nu} dx^\mu dx^\nu
\end{equation}

Usually the sum is just implicitly assumed by Einstein notation [1].

\begin{equation}
ds^2 = g_{\mu \nu} dx^\mu dx^\nu
\end{equation}

Here we can actually see that for some simple choice of \( g_{\mu \nu} \) we can get the distances by Pythagoras' theorem. Specifically for the case when the metric tensor matrix is a unit matrix.

\begin{equation}
ds^2 = dx_1^2 + dx_2^2 + dx_3^2
\end{equation}

For a flat space-time like space, we get something similar but with the exception that the time coordinate component is with a negative sign.

\begin{equation}
ds^2 = - dx_0^2 + dx_1^2 + dx_2^2 + dx_3^2
\label{flat}
\end{equation}

Here I used the \( (- + + +) \) signature, but signs can actually be flipped without changing the geodesics, and in some cases, like for particle physics, it makes more sense to use the opposite \( (+ - - -) \) signature.

Going back to the main question of computing distances, to compute the length between 2 points, A and B, along some path \( x^i(t) \) we simply need to sum the infinitesimal distances together using an integral:

\begin{equation}
l = \int_A^B \sqrt{g_{\mu \nu} dx^\mu dx^\nu} = \int_A^B \sqrt{g_{\mu \nu} dx^\mu dx^\nu} \frac{dt}{dt} = \int_A^B \sqrt{g_{\mu \nu} \frac{dx^\mu}{dt} \frac{dx^\nu}{dt}} dt
\end{equation}

Where \( \frac{dx^i}{dt} \) is simply how fast the coordinate x changes with respect to the path parameter ("clock"), in some sense can be interpreted as the velocity.

Now our main question is how do we minimize the path length? Here is where we introduce a thing called calculus of variations, which is roughly speaking a way to find how a functional (distance) changes by infinitesimally small variations of its input function (path). Such derivatives have similar properties to normal function derivatives. And in fact, similarly to calculus, to find the extremum of a function (min, max or stationary point), we simply need to equate the variation to 0.

Lagrangian description of a geodesic

There is an entire branch of physics related to variational principles, which states that any kind of physical system has a value it likes to minimize (or more generally make unchanging under small variations of path, i.e. "stationary"). That value is called action, and the function under the integral is called the Lagrangian function of the system. The branch of physics studying Lagrangian functions of systems is called Lagrangian mechanics.

In our case the Lagrangian can be written like this:

\begin{equation}
L = \sqrt{g_{\mu \nu} \frac{dx^\mu}{dt} \frac{dx^\nu}{dt}}
\end{equation}

Turns out we don't need the square root for the minimum of the functional to be a geodesic, and we can simply remove it from our geodesic Lagrangian:

\begin{equation}
L = g_{\mu \nu} \frac{dx^\mu}{dt} \frac{dx^\nu}{dt}
\end{equation}

The proof of this you can find here [2]. The only difference such simplification makes is that the parametrization of the path might be different, but the path itself will be the same. Also notably, the equations for this specific case turn out to be the same, meaning the parametrization is also the same.

Additionally, we want this with a 1/2 factor, to simplify the equations down the line.

\begin{equation}
L = \frac{1}{2} g_{\mu \nu} \frac{dx^\mu}{dt} \frac{dx^\nu}{dt}
\end{equation}

So, our goal right now is to minimize this functional:

\begin{equation}
S = \int_A^B \frac{1}{2} g_{\mu \nu} \frac{dx^\mu}{dt} \frac{dx^\nu}{dt} dt
\end{equation}

In general the minimum of a functional like this can be found by applying the Euler-Lagrange equations [3]:

\begin{equation}
\frac{\partial L}{\partial x^i} - \frac{d}{dt} \frac{\partial L}{\partial \frac{dx^i}{dt}} = 0
\label{el}
\end{equation}

You can find a more detailed derivation here [4]

Let's derive the Euler-Lagrange equations for our geodesic Lagrangian (keep in mind that there is an equation for each coordinate \( x^i \) ):

\begin{equation}
\frac{\partial L}{\partial \frac{dx^i}{dt}} = \frac{1}{2} \frac{\partial }{\partial \frac{dx^i}{dt}} g_{\mu \nu} \frac{dx^\mu}{dt} \frac{dx^\nu}{dt} = \frac{1}{2} g_{i \nu} \frac{dx^\nu}{dt} + \frac{1}{2} g_{\mu i} \frac{dx^\mu}{dt} = g_{i \nu} \frac{dx^\nu}{dt} \label{el0}
\end{equation}

Then we take the derivative with respect to the path parameter:

\begin{equation}
\frac{d}{dt} \left( g_{i \nu} \frac{dx^\nu}{dt} \right) = \frac{d g_{i \nu} }{dt} \frac{dx^\nu}{dt} + g_{i \nu} \frac{d^2x^\nu}{dt^2} = \frac{d g_{i \nu} }{dx^\mu} \frac{dx^\mu}{dt} \frac{dx^\nu}{dt} + g_{i \nu} \frac{d^2x^\nu}{dt^2} \label{el1}
\end{equation}

And lastly:

\begin{equation}
\frac{\partial L}{\partial x^i} = \frac{1}{2} \frac{d g_{\mu \nu} }{dx^i} \frac{dx^\mu}{dt} \frac{dx^\nu}{dt} \label{el2}
\end{equation}

Substituting \eqref{el1} and \eqref{el2} into Euler-Lagrange equations \eqref{el} leads us to the equation of a geodesic:

\begin{equation}
\frac{1}{2} \frac{d g_{\mu \nu} }{dx^i} \frac{dx^\mu}{dt} \frac{dx^\nu}{dt} - \frac{d g_{i \nu} }{dx^\mu} \frac{dx^\mu}{dt} \frac{dx^\nu}{dt} - g_{i \nu} \frac{d^2x^\nu}{dt^2} = 0
\end{equation}

Multiplying by the metric tensor inverse \( - g^{i \nu} \) we get:

\begin{equation}
\frac{d^2x^i}{dt^2} + g^{i \nu} \left( \frac{d g_{i \nu} }{dx^\mu} - \frac{1}{2} \frac{d g_{\mu \nu} }{dx^i} \right) \frac{dx^\mu}{dt} \frac{dx^\nu}{dt} = 0
\end{equation}

And that's our final system of equations for a geodesic, we could of course also substitute the Christoffel symbols here, but for our application there is no difference. Of course, we could just use these equations for tracing geodesic rays and call it a day, but unfortunately this would require computing a whole lot of derivatives (in 4d space-time it's 64 of them to be specific), either manually, or by using numerical differentiation. Thankfully there is a way to avoid this, and in fact simplify the entire algorithm! (At a slight performance cost)

Hamiltonian description of a geodesic

So here comes the star of the show - Hamiltonian mechanics. Hamiltonian equations of motion have a really nice form which allows to easily write a computer program that integrates them by using Euler integration.

\begin{equation}
\frac{dp^i}{dt} = - \frac{\partial H}{\partial x^i}
\end{equation}
\begin{equation}
\frac{dx^i}{dt} = \frac{\partial H}{\partial p^i}
\end{equation}

Where \( p \) is the so called generalized momentum, it's the derivative of the Lagrangian with respect to the coordinate path parameter derivative:

\begin{equation}
p_i = \frac{\partial L}{\partial \frac{dx^i}{dt} }
\label{momentumdef}
\end{equation}

And to get the Hamiltonian itself you need to apply the Legendre Transform [6] on the Lagrangian:

\begin{equation}
H = \sum_{i}^N p^i \frac{dx^i}{dt} - L
\label{legandre}
\end{equation}

A different derivation of Hamilton's equations of motion can be found here [5].

And for our case the momentum would be the following, which we already computed when writing down the Euler-Lagrange equations \eqref{el0}, and given the definition of generalized momentum \eqref{momentumdef}:

\begin{equation}
p_i = g_{i j} \frac{dx^j}{dt}
\label{momentum}
\end{equation}

To get the "time" derivatives you simply need to multiply both sides by the metric tensor inverse:

\begin{equation}
\frac{dx^i}{dt} = g^{i j} p_j
\label{dxdt}
\end{equation}

And the Hamiltonian itself:

\begin{equation}
H = \sum_{i}^N \frac{dx^i}{dt} p i - L = g{i j} \frac{dx^i}{dt} \frac{dx^j}{dt} - \frac{1}{2} g{i j} \frac{dx^i}{dt} \frac{dx^j}{dt} = \frac{1}{2} g{i j} \frac{dx^i}{dt} \frac{dx^j}{dt} = L
\end{equation}

Turns out that for this simple choice of a geodesic Lagrangian, the Hamiltonian is equal to the Lagrangian!

Also, we want to know the Hamiltonian as a function of the generalized momentum by substituting \eqref{dxdt} into the Hamiltonian equation:

\begin{equation}
H = \frac{1}{2} g_{i j} \frac{dx^i}{dt} \frac{dx^j}{dt} = \frac{1}{2} g^{i j} p_i p_j
\label{hamiltonian}
\end{equation}

While the equations of motion will simply be:

\begin{equation}
\frac{dp_i}{dt} = - \frac{\partial H}{\partial x^i}
\label{eqmotion1}
\end{equation}

\begin{equation}
\frac{dx^i}{dt} = g^{i j} p_j
\label{eqmotion2}
\end{equation}

This is all we need to write a numerical geodesic integrator!

Writing a Hamiltonian geodesic tracer in GLSL

You might have noticed that in the final Hamilton's equations of motion I didn't write out \( \frac{\partial H}{\partial x^i} \), this is actually important! We want to keep the derivative of the Hamiltonian as is, because then instead of computing the 64 derivatives of the metric tensor, we only need 4 to find the Hamiltonian gradient. This is the main simplification of the geodesic tracing algorithm.

Here we will use the GLSL shading language, since it has variables and functions which map quite well to the mathematical operations we will perform here. On top of that we can easily then make a real time GR visualization shader.

First of all, we need a function that evaluates the metric tensor at a 4d point in space and time. Let's use the Alcubierre warp drive [7] metric as an example, since it is quite simple.

mat4 Metric(vec4 x)
{
  //Alcubierre metric
  const float R = 1.0;
  const float sigma = 35.0; 
  const float v = 1.1;

  float x = v*x.x;
  float r = sqrt(sqr(x.y - x) + x.z*x.z + x.w*x.w);
  float f = 0.5*(tanh(sigma*(r + R)) - tanh(sigma*(r - R)))/tanh(sigma*R);
  float gtt = v*v*f*f - 1.0;
  float gxt = -v*f;

  return mat4(gtt, gxt,  0,  0,
              gxt,   1,  0,  0,
                0,   0,  1,  0,
                0,   0,  0,  1);
}

In our case x is a 4D vector representing position. The first component x.x or x[0] being time. As an output we get a 4 by 4 matrix represented by mat4 in GLSL.

Then we need to write down the Hamiltonian \eqref{hamiltonian}. The Hamiltonian is a function that takes 2 things, the position in space-time, and the 4d momentum vector, and outputs a scalar.

float Hamiltonian(vec4 x, vec4 p)
{
  mat4 g_inv = inverse(Metric(x));
  return 0.5*dot(g_inv*p,p);
}

As a bonus here is the Lagrangian

float Lagrangian(vec4 x, vec4 dxdt)
{
  return 0.5*dot(Metric(x)*dxdt,dxdt);
}

Surprisingly enough that's it, GLSL already has a matrix inverse function inverse(), on top of it the Hamiltonian is just the dot product (in GLSL sense) of g_inv*p and p, which are the contravariant and covariant momentum vectors respectively. The contravariant momentum actually just being the time derivative of the coordinate dxdt, i.e. dot(dxdt,p).

After this we need to compute the 4D gradient of the Hamiltonian. We can do this by using a forward numerical difference in all 4 spacial directions, using some small value eps:

vec4 HamiltonianGradient(vec4 x, vec4 p)
{
  const float eps = 0.001;
  return (vec4(Hamiltonian(x + vec4(eps,0,0,0), p),
               Hamiltonian(x + vec4(0,eps,0,0), p),
               Hamiltonian(x + vec4(0,0,eps,0), p),
               Hamiltonian(x + vec4(0,0,0,eps), p)) - Hamiltonian(x,p))/eps;
}

Now that we have the Hamiltonian gradient, we can finally write down the equation of motion \eqref{eqmotion1} \eqref{eqmotion2} integration code

vec4 IntegrationStep(inout vec4 x, inout vec4 p)
{
  const float TimeStep = 0.1;
  p = p - TimeStep * HamiltonianGradient(x, p);
  x = x + TimeStep * inverse(Metric(x)) * p;
}

You might ask, "wait, that's it?", and indeed that is all you need to integrate the geodesic. Of course, it is quite slow since we do a whopping 6 matrix inverse evaluations, which can be optimized down to 1, by replacing most Hamiltonians with Lagrangians which don’t have inverses, since they are equal. Even better is to have the metric inverse already computed analytically, but it’s not possible for every metric, especially for an implicitly defined one.

There is of course the last problem, while initializing the space-time position is easy, how do we initialize the value of the momentum vector p when starting to trace?

Before tracing the geodesic, you can use the equation \eqref{momentum}

p = Metric(x) * dxdt;

But what is dxdt? It's nothing more than the 4D direction the ray moves inside space-time. There are 3 categories the directions can fall into:

• Time-like, when \( A < 0 \)
• Null, when \( A = 0 \)
• Space-like, when \( A > 0 \)

Where \(A\) is
\begin{equation}
A = g_{\mu \nu} \frac{dx^\mu}{dt} \frac{dx^\nu}{dt}
\end{equation}

or in GLSL, A = dot(Metric(x) * dxdt, dxdt)

When simulating how light travels we just want null directions, which lead to null geodesic solutions. On the other hand, if you want to model an object moving slower than light you need a time-like geodesic. And space-like geodesics for tachyonic stuff, which doesn't happen in real life, so we ignore it.

So, assuming a flat space metric \eqref{flat} and some 3D direction in space our p for a light ray would be

vec4 GetNullMomentum(vec3 dir)
{
  return Metric(x) * vec4(1.0, normalize(dir));
}

And the inverse of this operation

vec3 GetDirection(vec4 p)
{
  vec4 dxdt = inverse(Metric(x)) * p;
  return normalize(dxdt.yzw);
}

So, in the end the final simple algorithm will look like this:

void TraceGeodesic(inout vec3 pos, inout vec3 dir, inout float time)
{
  vec4 x = vec4(time, pos);
  vec4 p = GetNullMomentum(dir);

  const int steps = 256;
  for(int i = 0; i < steps; i++)
  {
    IntegrationStep(x, p);
    //you can add a stop condition here when x is below the event horizon for example
  }

  pos = x.yzw;
  time = x.x;
  dir = GetDirection(p);
}

Essentially this is just a 4D ray marching algorithm where the direction of the ray changes every step. In this specific case the size of the step also changes, which can be avoided by normalizing the momentum p = normalize(p). This only changes the step length, and doesn't change the geodesic path, i.e., it works just like a dynamic reparameterization of the path. The time step of the integration can also be varied depending on the metric used. For example, in the case of black holes I change the time step proportionally to the distance to the event horizon, so that the accuracy of the geodesic is roughly proportional to the curvature of space. This is an important optimization to get accurate results, while keeping the computational cost relatively small.

Example shadertoy code:

mat4 diag(vec4 a)
{
    return mat4(a.x,0,0,0,
                0,a.y,0,0,
                0,0,a.z,0,
                0,0,0,a.w);
}

mat4 Metric(vec4 x)
{
    //Kerr-Newman metric in Kerr-Schild coordinates 
    const float a = 0.8;
    const float m = 1.0;
    const float Q = 0.0;
    vec3 p = x.yzw;
    float rho = dot(p,p) - a*a;
    float r2 = 0.5*(rho + sqrt(rho*rho + 4.0*a*a*p.z*p.z));
    float r = sqrt(r2);
    vec4 k = vec4(1, (r*p.x + a*p.y)/(r2 + a*a), (r*p.y - a*p.x)/(r2 + a*a), p.z/r);
    float f = r2*(2.0*m*r - Q*Q)/(r2*r2 + a*a*p.z*p.z);
    return f*mat4(k.x*k, k.y*k, k.z*k, k.w*k)+diag(vec4(-1,1,1,1));
}

float Hamiltonian(vec4 x, vec4 p)
{
    mat4 g_inv = inverse(Metric(x));
    return 0.5*dot(g_inv*p,p);
}

/*
float Lagrangian(vec4 x, vec4 dxdt)
{
    return 0.5*dot(Metric(x)*dxdt,dxdt);
}
*/

vec4 HamiltonianGradient(vec4 x, vec4 p)
{
    const float eps = 0.001;
    return (vec4(Hamiltonian(x + vec4(eps,0,0,0), p),
                 Hamiltonian(x + vec4(0,eps,0,0), p),
                 Hamiltonian(x + vec4(0,0,eps,0), p),
                 Hamiltonian(x + vec4(0,0,0,eps), p)) - Hamiltonian(x,p))/eps;
}

void IntegrationStep(inout vec4 x, inout vec4 p)
{
    const float TimeStep = 0.15;
    p = p - TimeStep * HamiltonianGradient(x, p);
    x = x + TimeStep * inverse(Metric(x)) * p;
}

vec4 GetNullMomentum(vec4 x, vec3 dir)
{
    return Metric(x) * vec4(1.0, normalize(dir));
}

vec3 GetDirection(vec4 x, vec4 p)
{
    vec4 dxdt = inverse(Metric(x)) * p;
    return normalize(dxdt.yzw);
}

void TraceGeodesic(inout vec3 pos, inout vec3 dir, inout float time)
{
    vec4 x = vec4(time, pos);
    vec4 p = GetNullMomentum(x, dir);

    const int steps = 256;
    for(int i = 0; i < steps; i++)
    {
        IntegrationStep(x, p);
        //you can add a stop condition here when x is below the event horizon for example
    }

    pos = x.yzw;
    time = x.x;
    dir = GetDirection(x, p);
}

void mainImage(out vec4 fragColor, in vec2 fragCoord) {
    vec2 uv = 2.0 * (fragCoord - 0.5 * iResolution.xy) / max(iResolution.x, iResolution.y);

    vec3 RayPos = vec3( 0.0,  0.0,  32.0);
    vec3 RayDir = vec3(uv.x, uv.y, -1.0);
    float  Time = 0.0;

    RayDir = normalize(RayDir);

    TraceGeodesic(RayPos, RayDir, Time);

    fragColor = vec4(texture(iChannel0, RayDir).rgb, 1.0);
}

You can check out this Shadertoy implementation to see some of the optimizations, like variable timestep, replacing Hamiltonians with Lagrangians, using a symmetric matrix inversion function (a bit faster), reusing some of the computed values (restart if the Shadertoy is black):

mat4 diag(vec4 a)
{
    return mat4(a.x,0,0,0,
                0,a.y,0,0,
                0,0,a.z,0,
                0,0,0,a.w);
}

mat4 Metric(vec4 x)
{
  //Kerr-Newman metric in Kerr-Schild coordinates 
  const float a = 0.8;
  const float m = 1.0;
  const float Q = 0.0;
  vec3 p = x.yzw;
  float rho = dot(p,p) - a*a;
  float r2 = 0.5*(rho + sqrt(rho*rho + 4.0*a*a*p.z*p.z));
  float r = sqrt(r2);
  vec4 k = vec4(1, (r*p.x + a*p.y)/(r2 + a*a), (r*p.y - a*p.x)/(r2 + a*a), p.z/r);
  float f = r2*(2.0*m*r - Q*Q)/(r2*r2 + a*a*p.z*p.z);
  return f*mat4(k.x*k, k.y*k, k.z*k, k.w*k)+diag(vec4(-1,1,1,1));    
}

This is, of course, not the limit for optimization. The main other optimization is computing analytical inverses of the metric. For a large class of metrics you can use the Sherman–Morrison formula [9]

Some useful things to keep in mind

Since we used numerical finite differences, the results can actually depend quite a lot on the relative values of the float numbers. For example, the accuracy of the numerical derivatives is a lot lower far from the coordinate system center, so you'd need to vary the size of eps to avoid excessive numerical noise. Also, metrics quite often have numerical singularities which you should avoid, unless you want to get NaN results.

I usually avoid metrics in spherical coordinates due to their polar axis singularity, which has strong visual effects which are extremely hard to avoid even with a tiny varying timestep, although such metrics are usually mathematically simpler and allow for larger timesteps without breaking the look of the Black hole. For spherically symmetric metrics, like non-spinning black holes and wormholes there is a trick to avoid the polar singularity altogether! The thing about spherical symmetry is that the geodesic is always moving inside a 2d plane, which can be mapped to the equatorial plane of the coordinate system, basically reducing the 3d + time problem to 2d + time. (Scott Manley has a video explaining how he rendered wormholes, in there he used this trick + precomputing a lookup table to simplify the computation by a lot)

I've also used the dimensionality reduction trick in my shadertoy wormhole:

And finally you could always derive the equations analytically, while this is the most annoying method it is usually the fastest performance-wise. A compromise solution would be derive the equations automatically, this approach is used by geodesic_raytracing made by James Berrow (you should follow him on Twitter, he has a lot of cool stuff on this topic).

Figuring out if the ray has fallen inside the event horizon is actually not trivial, and there is no universal method, and while you could just set the color to 0 if the ray is below the event horizon surface, this is incorrect when viewing things from inside the black hole. Tracing the rays should also be done backwards in time, since we trace the rays from the camera, not to the camera, this has a noticeble effect on the resulting render, if not done this also results in completely dark renders inside of black holes, even though light does exist under the event horizon, and can reach from the outside.

Redshift in General Relativity can be computed simply from the ratio of the dot products of the velocity of the object at emission/absorption times the momentum of the photon, where the momentum of the photon is parallel to the photons direction of movement. The momentum needs to be parallel transported along the geodesic from emission to absorption, and the good news is that we can just use the geodesic 4-velocity instead, since it also is technically parallel transported along the geodesic and is pointed in the direction of movement. (Just in case, this also means we need to avoid renormalizing the generalized momentum when integrating the equations)

Also notably, the dot product in General Relativity is simply defined from the metric tensor \( g(u,v) \)

\begin{equation}
u \cdot v = g_{\mu \nu} u^\mu v^\nu
\end{equation}

Finally, if you only want to compute geodesics for Schwarzschild black holes, you can simply use the equation for a particle with mass 1 in a certain classical force field, details are in
this blog post.

Conclusions

Kerr black holes

Alcubierre warp fields and wormholes

Volumetric accretion disks around a Kerr black hole

Fast volumetric ray tracing with geodesics is quite difficult, and we needed to separate the ray marching loop into 2 loops, main loop being the geodesic steps, and the second loop being the volumetric substeps. Since we also use blue noise, it was necessary to keep the steps uniform along the geodesic, otherwise there would be clear artifacts in the volume, which required a few tricks with having a variable number of substeps per geodesic step.

Combining this with SDF's is somewhat easier, you need to vary the geodesic step to be the min() between the current step size and the SDF. Using this I've also tried to make a really simple path tracer in Unity with a Kerr black hole, naturally it was quite slow.

Note that rendering of moving objects is waaay harder, and requires either to have a space-time SDF, or some insane acceleration structure for triangles. And on top of that the entire history of the scene's past needs to be kept in memory, the only simple cases are when the moving objects can be represented as analytical functions you can sample in space and time, like the volumetric accretion disk in SpaceEngine.

References

• [1] Einstein notation
• [2] Equivalence of squared Lagrangian to Lagrangian
• [3] Euler-Lagrange equations
• [4] Euler-Lagrange equations derivation
• [5] Hamiltonian equations derivation
• [6] Legendre Transform
• [7] Alcubierre metric
• [8] Kerr–Newman metric in Kerr-Schild coordinates
• [9] Sherman–Morrison formula
• [10] Space-time path tracer

This is a continuation of this discussion.

Watsisname: I did try directly computing orbits with and without perturbation just to try to get some intuition for how the perturbation effect works.  For starters I assumed co-planar, circular orbits, and also that only Neptune perturbs Uranus and not vice versa (which probably matters very little).

Well, sorry for reviving this Neptune discovery thing again, but I felt I needed more comprehension on how Uranus's residuals (produced by the gavitational perturbations of Neptune) actually behaved. Since I'm not as hardworking (nor as clever) as Watsisname I didn't create my own N-body integrator (you are mad Watsisname, awesome mad). So I tried Rebound, an N-body integrator library used in Python and used it to program a bunch of scenarios and graphs with that. I wanted also to learn how to use Rebound so I could also experiment myself in my modest "coding-lab" on whatever topic comes next, and this was a good excuse to learn.

What I did? Just reproduced Watsisname's residuals graph but now using the actual orbital parameters (Watsisname assumed two circular obits without inclination for Uranus nad Neptune and I've used JPL's current best orbital parameters for the integration). As he understood there is basically no large difference in the result. Rebound handles many different kind of N-body integrators but I choose the IAS15 which is fast and has adaptative time-steps (so that it get's more detailed when larger interactions are possibly occuring). I wanted to know how exactly does this curve was generated (Why that periodicity between oppositions? for example) so either way I want to explain what exactly is going on here physically, because it's quite interesting.

I did this video with the entire explanation but I will write some part of it here for the sake of clarity (you have to stop the video to many times and YouTube's lower progress bar might get in the middle of the text when pausing).

What are we seeing?

Let me explain the parts of the simulation and the color-coding here. So we have an orbital diagram of the Solar System as viewed from above (with both axis measured in astronomical units) to the left and 4 graphs to the right (with the current year represented in the x-axis for all of them). Let's enumerate the graphs on the right from top to bottom; In the 1st graph we have two curves, the cyan for Uranus's residuals in its longitude and the darker blue one for Neptunes residuals. What is longitude? is just the position angle of a planet as viewed from the Sun. What are the resiudals in longitude? Well, they are the angular difference between Uranus's longitude if there was no Neptune (the theoretical position angle predicted for Uranus by scientists before the discovery of the perturbing planet: Neptune) and Uranus's longitude if Neptune exists (the actual empirical observation of the longitude of Uranus at the time). So the residuals are just a measure of how far theoretical predictions are from the actual observations of the time (as Neptune created that difference by its influence). The residuals for Neptune are just plotted in dark blue for comparison (no historical situation needed those but I wanted to know how Uranus's perturbs Neptune also). Even if the orbital diagram in the video does not show the other planets of the Solar System (just Uranus and Neptune) they where used to compute the longitude of Uranus and Neptune (in both the non-Neptune and Neptune scenarios needed for calculating the residual at any given time) since historically Le Verrier computed the residuals of the perturbing unknown body after taking into account the perturbations produced by the other 6 planets. This first graph has both arcsecond and arcminute tick marks in the y-axis. The 2nd graph shows the intensity of the gravitational force between Uranus and Neptune (in billions of Newtons or Giga-Newtons) as they get close and go away (in dark grey), but the important thing is the tangential component of this force (tangential to Uranus's orbit) whos intensity is represented by the magenta curve. Why we care for this particular component of the force vector? Because that force is the one generating an acceleration in the direction of motion and therefore is the one contributing to change the orbital speed of Uranus and in consequence is the one we should pay atention if we want to understand the behaviour of the residuals. The perpendicular component is basically pointing in the Sun-Uranus line (not entirely because the orbit is not a perfect circle) which means that the acceleration generated by this force does not contribute to change the speed of the object but only the direction at which the velocity vector is pointing. The 3th graph shows in green the excentricity of Uranus's orbit (yep it turns out it is very important for the explanation). The 4th graph shows the difference between the orbital speed of Uranus if there was no Neptune and the orbital speed considering the influencing presence of Neptune. The speed difference of the 4th graph is tiny (can be measured in meters per second) but it can accumulate over time to separate both versions of Uranus (the perturbed and unperturbed ones). All the graphs show a dotted vertical line for each opposition year (when Neptune and Uranus have their closest approach).

The orbital diagram is very easy to understand; we have both the Sun as a yellow star icon and the orbits of Uranus and Neptune. But we have two dots for Uranus that are separating from each other; those are the perturbed/real/observed Uranus (magenta dot) and the unperturbed/Neptuneless/predicted version of Uranus (cyan dot). The angle between them as seen from the Sun is the value of Uranus's residual (that's why I drew lines coming from the Sun to both versions of Uranus). The only tricky thing here is that I have exaggerated this angle by a factor of 100 in the orbital diagram (as you can see the residuals fluctuate in the range of the 10ths and 100ths of arcseconds, not sufficiently large for our viewing purpouses in the orbital diagram). Finally I added some arrows; a dark grey arrow (it is nearly invisible) which show the gravitational force vector, and a magenta arrow showing the projection of this vector to the tangent of the orbit at Uranus (the tangential component of the gravitational force). A black filled circle in Uranus's orbit is marking the perihelion of the perturbed Uranus (the opposite side of the orbit is where the aphelion has to be then).

What is going on?

The simulation starts the 1st of January of 1777 and evolves through 550 years of history to the year 2327. You can see how the residuals are as big as 200 arcseconds (more than 3 arcminutes) by the time Le Verrier and Galle discover Neptune. We can see that these residuals turn as big as to be naked eye visible after a few more years and if they tried to wait for some centuries their prediction would have been extremely accurate and based on clear evidence as the residuals get higher and higher. But no wait was needed, just in a few decades since Uranus's discovery the residuals were large enought to be detected by the telescopes of the time above any shadow of doubt (above what was reasonable by accounting uncertainties and instrumental errors).

Why did the residuals got so large so quickly?
We talked about this previously and two easy missconceptions can be shown;

#1 We would expect to see the larger perturbing force when Uranus and Neptune are closer (at opposition).

This turns out to be false. Why? because at opposition nearly all the gravitational force is directed perpendicular to the direction of motion of Uranus and thus even if the gravitational tug of Neptune is the strongest it has no tangential component and the accelerations does not contribute to change the orbital speed of Uranus (only contributes to change the direction of the velocity vector slightly).

#2 Then we should expect the perturbation to be maximized when the entire gravitational force is directed tangentially (as in the next image).

False again. Yes, in this configuration Neptune would be pulling Uranus in the exact direction of motion but this happens in a configuration where Uranus is quite far away from Neptune so the pull, even if 100% tangential, is not as strong.

The truth lies in a place between both situations, where there is plenty of force directed in Uranus's direction of motion since both planets are relatively close and in a sufficiently low incidence as to better influence each other. You can see this in the force graph of the simulation; the peak in tangential force does not occur when the net force is tangent to the orbit of Uranus nor when we are at opposition (at opposition we can see how the tangential force turns to zero but the gravitational force peaks).

Le Verrier only used observations of Uranus made between 1781 and 1821 (zoomed version of the graphs for that period below). But he had some luck since in those 40 years the peak in Neptunes perturbing tangential force happened (in 1812 as can be seen in the simulation). If not he would have had to wait for around 172 years to see the next. Also opposition happenned in 1820. This is the reason humanity had such strong of a signal in that time and the residuals where evidently anomalous.

Why the residuals oscillate between each Neptune's opposition?
Well let's imagine a completely different problem. Let's suppose there is only one planet, Uranus. And let's create two different scenarios; in one Uranus has a low orbital excentricity and in the other Uranus's orbit is slightly more excentric. The other orbital parameters are mantained exactly the same. So, how the difference between the longitude of the low-excentricity Uranus and the longitude of the excentric Uranus is going to look like? Let's start both at aphelion. Here the excentric Uranus has to be farther from the Sun (because its aphelion is elevated with respect to the low-excentricity Uranus), this means that the low-excentricity Uranus will be faster at this point (aphelion) and has to traverse smaller distances. This makes the low-excentricity Uranus advance the excentric Uranus for a little (the residual angle get's larger). But there is another effect here. Both are falling towards perihelion (since both are at aphelion), so both Uranus are gaining speed, but the excentric one is gaining more of that extra speed as it falls from a higher elevation (a place with a higher gravitational potential). So, at some point (at perihelion) this effect would be so large that the excentric Uranus will have shortened the distance to the other Uranus and gained much more speed. Both coincide at perihelion but the excentric Uranus is traveling faster than the other and because of inertia it now surpasses the other (the residual angle start to go negative here). Then the reverse holds. As both Uranus climb the gravitational dwell of the Sun they lose speed in the way to aphelion, but the excentric Uranus gets quicker to higher potentials (since is going faster than the other) so it loses speed faster than the low-excentricity Uranus (the residual angle which is negative here starts to close until it reaches zero again). We arrive to aphelion again and the low-excentricity Uranus has climbed a lot and slowed down below the speed of the other Uranus since a while. The other Uranus is now profiting his larger speed to overcome the excentric Uranus. And the cycle repeats. This creates a sinusoidal behaviour for the residuals we have calculated here (the angular separation between the planets increases, then decreases, turns negative and increases in the negative values, finally it returns to zero after one cycle is completed). This oscillation of the residual has the same period as Uranus's orbital period. But why there are two oscillations for each Neptune opposition? Well, because Uranus has nearly a perfect 2:1 ratio in orbital periods with Neptune (a resonance). Just by chance each time Neptune completes a cycle Uranus has done 2 of them.

What does Uranus excentricity has to do with the perturbation of Neptune?
Nothing. Nothing at all. But Neptune is the responsible for the change of the eccentricity of Uranus's orbit and thus the cause of that half-period oscillation of the residuals between oppositions. Neptunes role here, therefore, is to excite that oscillation only (increasing it in amplitude). Neptune is just making the non-perturbed Uranus having to compete with a new Uranus with an orbit with higher excentricity (just as the problem we talked before). Each Neptune's opposition the perturbed Uranus gets a new more excentric orbit and thus the oscillation increases in the residuals (thanks to Neptune's kicks). Le Verrier would not have seen the oscillation itself (since that was accounted as part of the predicted longitude using the known excentricity of the time), but he noticed the arrival of Neptune as a "excentricity excitator", and only after that, with a new more excentric Uranus, the oscillations in the residuals are revealed.

But how is Neptune really exciting Uranus's orbital excentricity?
Now that we know that the half-period oscillations of the residuals are due to a different excentricity between the non-perturbed and perturbed Uranus we can adress how Neptune is making one more excentric that the other. Here comes my third missconception;
#3 The pull of Neptune just before opposition is of the same magnitude but opposite direction as the pull of Neptune just after the opposition, thus the perturbing effect of the later should cancel the perturbing effect of the former. Perturbations are specular symetric with respect to opposition and thus cancel out.
Here we are claiming there is some symetry in the dynamics of the system for which an effect counteracts the other and no global perturbation is observed. It is also suggested by the odd aspect of the tangential force curve (the curve seems to be two times mirrowed with respect to opposition) But this is not what happens obviously since the residuals increase and perturbation is noticiable long after opposition. So there has to be some symetry break here. And in fact there is one;

Look in the video at the close encounter locations. Each opposition is happening when Uranus is on the way to aphelion (the half of its orbit where it is climbing the potential dwell of the Sun). This means that before opposition Uranus is traveling faster than after opposition (and not the other way around) because it is losing orbital speed as it goes to the aphelion. Neptune (which is always slower than Uranus) has the chance to influence Uranus slightly more after opposition than before opposition (here we have the symetry break!) because it has more time to influence Uranus in this part than in the other. While the before-opposition-influence tries to speed up Uranus the after-opposition-interaction tries to slow down it; since the second has more time to actuate we conclude that Neptune is able to get Uranus's speed decreased more than the previous increase and therefore the result long after Neptune's opposition is a slower Uranus.

But what exactly a slower Uranus means in this part of the orbit? Well, since it is slower than it should be before aphelion this means that it has a bit less kinetic energy to climb in the potential dwell of the Sun; which translates to an decrease in the elevation of the aphelion of its orbit; which is just a complicated way of saying that the excentricity has decreased!

But wait a moment, the excentricity is actually increasing (not decreasing). So there has to be an other effect here that counteracts the one just mentioned. Could the fact that Uranus's perihelion and Neptune's perihelion are located in opposite sides of the Solar System has something to do with this? Indeed this is a strange coincidence we havent mentioned yet; Neptune's "argument of perihelion" is at 276.3º while Uranus's argument of perihelion is 97.0º, so there are 179.3º of separation between them (suspiciously close to 180º which would mean they are opposite). Let's see what implications does it has; since the perihelions of both planets are located in opposite directions we have that Neptune's aphelion coincides with Uranus's perihelion (there is minimal interaction then in this region compared to any other configuration), but in exchange we have that Neptune's perihelion is oriented where Uranus's aphelion is located (which implies closer interactions than any other rotation of the ellipse shape of the orbits). So, if opposition is happening for Uranus when it is in the track between perihelion and aphelion opposition is happening for Neptune when it is in the track between its aphelion and perihelion (and viceversa). So, not only Uranus is quicker before opposition than after opposition because of the way the ellipse of its orbit is oriented but also Neptune is slower before opposition than after opposition because of the way his ellipse is oriented. Which means that both the fact that oppositions happens for Uranus in the perihelion-aphelion path and the fact that opposition also happens in Neptune's aphelion-perihelion path, are both responsable of Neptune having more time for influencing (slowing down) Uranus. Both effects are breaking the symetry, but both would tend to circularize the orbit of Uranus and not to increase its excentricity (as we are seeing in the graphs).

So, we are still left with the problem of how Neptune is increasing Uranus's orbital excentricity exactly? Since none of the previous effects plays a role other than circularizing Uranus's orbit, there has to be an effect that counteracts both and allows for plenty of room for a slight net increase in excentricity.

Finally the answer is; the perpendicular component of the force! Yep, we have disregarded this component because we thought that the important was the tangential force (applied in the direction of motion). We assumed that, because the perpendicular component changes the velocity vector in direction but not in magnitude, and we were focusing our attention in the increase or decrease of orbital speed without realizing that the change in appearence of the orbit can also be acomplished without touching the orbital speed but only changing direction (something that the perpendicular component of the force does). Let me explain in more detail; since oppositions always occur (at least in the time window represented in the video) when Uranus is climbing the potential dwell of the Sun, a "kick to the outside" (perpendicular to its motion and outside its orbit), like the one performed by Neptune (which is external to Uranus) will modify the velocity vector of Uranus, but only in direction. Now, with the same speed as if Neptune didn't existed, but pointing to a slighly outward direction, Uranus is able to have more of its motion pointing away from the Sun which means that if would climb farther its gravitational dwell, which means that the aphelion is been elevated, which means that the orbit is gaining excentricity! Finally! Think about it like this if you want: imagine the infinite line defined by the velocity vector of Uranus when he is in the perihelion-aphelion track and without Neptune present (and try to picture the distance of the Sun to this line). Now imagine the same line for the same situation except Neptune is present and at opposition; the force excerted by Neptune is directly related to the acceleration felt by Uranus. As we know acceleration is just the change in velocity but since the change is perpendicular to velocity in this case the velocity vector is just going to rotate a certain small angle. The infinite line created by this slightly rotated vector is now closer to the Sun. If you are able to take a moment to picture this it would be very important for the explanation. Since now the line of the extended velocity vector gets closer to the Sun, and since this line is tangent to the trajectory at that point we can conclude that the new trajectory gets closer to the Sun at perihelion (the perihelion has lowered and thus we have higher excentricity). I tried to make a scheme with photoshop to ilustrate the idea. As you can see opposition in this perihelion-aphelion track makes Uranus's perihelion go down and aphelion go up (excentricity increased).

As a brief summary for this section:

We are dealing with 3 effects, that have different consequences and happen because of different reasons:

1.  Since opposition is always happening in the perihelion-aphelion track of Uranus's orbit, Neptune has more time to slow down Uranus after opposition than before opposition where he's trying to increase Uranus's speed (using the tangential force). The net result is that Uranus gets slowed down by Neptune and because this happens in this part of the orbit. Uranus's orbit has a decrease in excentricity.
2. Since opposition is always happening in the aphelion-perihelion track of Neptune's orbit, Neptune has more time to slow down Uranus after opposition than before opposition where he's trying to increase Uranus's speed (using the tangential force). The net result is that Uranus gets slowed down by Neptune and because this happens in this part of the orbit. Uranus's orbit has a decrease in excentricity.
3. Since at opposition Neptune is pulling Uranus always to the right (according to a viewer placed in Uranus and facing in the direction of motion) by means of the perpendicular component of the force, Uranus is changing the direction of its velocity vector (without changing its magnitude) away from the Sun. The result is that, with the same speed pointing against a higher slope of the potential dwell of the Sun, Uranus is able to climb further to a new higher aphelion, which means that Uranus's orbital excentricity increases.

So the first 2 effects try to circularize Uranus's orbit, but the third is the one exciting the excentricity, and it does so as to balance and overcome the other two effects. Why? probably because the intensity of the perpendicular force during opposition is much higher than the intensity of the tangential force in any conceivable moment. Since the third effect is the only one caused by this perpendicular force he es able to overcome the others.

I want also to show you that in fact the first two effects break the symetry of the tangential force curve by mirror-reversing it (with a redder color);

Now you can see that it is not perfectly symetrical and in fact you can see that the tangential force in the after-opposition is mantained for much longer than before-opposition (and it is even stronger).

A last interesting note: to discard the fact that Uranus's perihelion is opposite to Neptune's is in any way playing a role here I tried to change the parameters of the simulation by establishing the same argument of perihelion for both Uranus and Neptune. In this way Uranus should be always at a safe distance from Neptune and the second effect would have no reason to exist. And indeed when I did this the only change was that the amplitude of the excentricity increased (but the increase in excentricity happened the same way in the same period and no decrease was shown as expected). This is because by doing this rotation of Neptune's orbit so that both planets are closer to the Sun in the same direction we have inhibited the second effect (which only damped the excentricity) so now the third effect has only to fight the first effect and thus the excentricity was more excited (the third effect was working more easily). Why this is interesting? Why I mentioned this first and second effect If I knew the one to blame was the third? Well, because I think that having both planets perihelia pointing in opposite directions is a controler for the out-of-control excentricity increase of Uranus's orbit and also because I suspect that this item has some important role in another future topic, Laplace resonance (since Io and Europa also have their perihelia pointing in opposite directions).

How is it possible that the Neptune-Uranus system is then stable if the interaction keeps increasing Uranus's orbital excentricity?
Aha! this is where it got me really interested. This has to do with the fact that Uranus and Neptune are close to the 2:1 resonance as I explained. This resonance means that opposition will always happen in the same part of the orbit and Neptune's influence will accumulate over time until it strips Uranus away from that region of the Solar System. I mean, that is what looks like to be happening in the video where the excentricity is just getting higher and higher. But remember this! in reality they are not in 2:1 resonance but just close to it (surprisingly this has a long-term shielding effect). The actual ratio between the periods is 1.96:1. That means that, for each orbit of Neptune, Uranus has made 1.96 orbits (nearly two but not), so for the two planets to coincide Uranus has to travel a little more to reach Neptune, meaning that the opposition would ocurr a little ahead of the last time it happened. So the 1.96:1 ratio is making the opposition location shift slowly in the same direction as the planet's motion. The shift is small (because we are too close to the 2:1 resonance where we would be at the exact same spot for all the oppositions) but in the long run you can see it driffting.

Now! remember that the previous arguments where we showed that there was a break in symetry and the influence of Neptune fighted to increase Uranus's orbital excentricity. The entire argument was supported in the idea that oppositions were happening in the way from perihelion to aphelion. But the video is tricky, it doesn't show what happens across milleniums (just centuries); with enought time you would start to notice that the opposition point shifts slowly to the point were it surpasses the aphelion point and then the entire previous argument gets reversed;

Neptune is able to influence Uranus in the falling-to-perihelion-track of its orbit now. The effect is such that if the opposition point reaches that other half of Uranus's orbit, it would damp the orbital excentricity for centuries until it gets the original values. Since I wanted to see this in action and prove it to you I made another animation (this time in gif format):

Here we wait for thousands of years. Each frame corresponds to one opposition (that's why Uranus and Neptune seem to travel togheter). As you can see the near resonance produces the opposition point to travel in the same direction as the planets in their orbits. Within the first centuries (the span of the video) eccentricity is increasing because this opposition point is behind aphelion still, even if it is traveling closer and closer to it. As the opposition point surpasses aphelion around the year 2700 the effect is reversed and exentricity starts to be damped. In the year 3600 it has reached the original excentricity and continues to decrease until, in the year 4700, the opposition point surpasses perihelion and excentricity starts to be excited again. This 4500 year cycle is what in turn mantains the long-term stability of the Solar System (at least for this region), thanks to the fact that Uranus and Neptune are not in a 2:1 resonance but close. Another interesting thing to note in this gif is the wooble of the perihelion (it was changing position itself all this time!).

Lastly I want to mention the fact that for all the graphs and orbital diagrams I took into account the perturbations of all the Solar System planets, except for the excentricity graph. If we include them in it we would have the next picture:

As you can see the rest of the Solar System only adds a bit of "noise" in Neptune's generated excentricity signal. The periods of Jupiter and Saturn can are much lower so they perform higher frequency perturbations in Uranus's excentricity. Those are relatively large but not large enought to mess with general trend dominated by Neptune's influence.

How is it possible that Neptune's influence on Uranus is stronger than Jupiter's or Saturn's?
This might seem a stupid question, but consider the following; Uranus-Neptune distance at opposition is 11 AU but Saturn-Uranus distance is 10 AU. Not only Saturn gets closer to Uranus but it is also much more massive than Neptune, thus more gravitational force should be at play there. Also consider Jupiter; it gets as close as 14 AU from Uranus. Okay Jupiter is 1,4 farther than Neptune at best and since gravitational interactions goes as the inverse square of distance this implies a 1,96 times weaker force than that excerted from Neptune. But wait! Jupiter is 18,7 times more massive than Neptune! So the reality is that Jupiter and Saturn both have stringer pulls on Uranus than Neptune has. This is a call for another missconception.
#4 Since Jupiter and Saturn pull stronger on Uranus than Neptune they should be the dominant perturbing forces in our graphs.
Well, weaker Neptune has two advantages that neither Jupiter nor Saturn have. First is the near resonance; neither Jupiter nor Saturn have such a close orbital ratio with Uranus as Neptune, this is why Neptune is able to increase the excentricity of Uranus continuosly during long time periods and in a consitent manner while Saturn and Jupiter make opposition with Uranus in more randomized places (since the resonance argument is no longer a rational number) increasing and decreasing the eccentricity without any cumulative power in the long-run. The second thing is time; Neptune has plenty of time to influence Uranus at near opposition since it moves slowly, Jupiter and Saturn just pass very quickly as they have the hurry to complete their orbits 13,9 and 5,6 times more rapidly than Neptune has.

Author(s): Propulsion Disk, FastFourierTransform, Watsisname
Original discussion started on the forum: here

Propulsion Disk:

The first theory of Neptune's existence was made by Galileo Galilei, on December 28, 1612, probably one of the most well known astronomers, was once looking at the night sky, searching for new planets. That night, he found something very interesting, it was what seemed to be a star that didn't twinkle, and therefore Galileo thought he had found a planet! He looked at again a year later to figure out that it was starting to go retrograde and then proved that it was indeed a planet, he wrote it in his drawings of the solar system with no name and showed it to the public, astronomers kept looking for Neptune and eventually found out that it was just a star that was in conjunction with Jupiter witch is what made it look like it was moving to Gallileo, and his small telescope, and was never fought of again until 1821.

In 1821, astronomer Alexis Bouvard was observing Uranus, he kept track of it for months and wrote on his tables that it was constantly getting perturbed by what could only be some unknown body, causing him to believe the existence of Neptune.

In 1844, john couch Adams was fascinated by Neptune's mysteries, and wanted to prove it's existence, he gathered a lot of data on the subject by himself, but mainly talked to other famous astronomers to prove Neptune's existence until 1846.

The final time Neptune was searched for was when Urbain Le Verrier started researching it, unaware of John's observations, he made calculations of where Neptune would be by finding Uranus, and watching where it was being pulled, and how much via newton's laws, then he went to make a note to astronomer Johann Gottfried Galle that said to use the refactor at the Berlin observatory and to use his calculations to find Neptune's position, when Galle got the note he looked for Neptune in the night sky and found it! 1 Degree of where Le Verrier fought it would be.

FastFourierTransform:

I would like to point out some statements that I belief are erroneous here, in your historical summary. Just for any future reader to take these into account (don't hate me please).

This discovery (of Galileo's observations) was published in Nature by Charles T. Kowal & Stillman Drake in 1980. The article is interesting on its own if anyone wants to research on the history of "the precovery of Neptune". It seems that Galileo pointed his telescope to Jupiter in a very appropriate moment; our current high-accuracy ephemerids show that Neptune was been occulted by Jupiter in January 1613 (a one in a century event at best). That means that the months before and after that event Neptune was close to Jupiter in the sky and any curious eye trying to learn something about Jupiter would have been able to possibly spot Neptune.

The diagram above shows Jupiter and the relative position of Neptune from 12/28 (December 28th) to 1/30 (January 30th). The second lucky thing is that Jupiter started in January 1613 it's retrograde motion (apparent from Earth). This is why in the diagram (where Jupiter is fixed) Neptune appears to return. The fact Neptune was slowing down and returning close to Jupiter's position is also good since that incremented the chances of Galileo to spot it.

Galileo is known to have been able to spot even magnitude 9 stars. Since Neptune never gets above magnitude 8 that means that Galileo indeed should have been able to see Neptune in some observation that moth. The field in the sky where the first observation was made is depleted of bright stars so this "fixa" must have been Neptune indeed.

Let's look at the actual notes done by Galileo.

Lets zoom in to the upper part of the second page of his notebook:

In the upper part he writes the time when the observation was made. It is December 27th, 1612, 15 hours and 46 minutes after noon. Galileo's way of formating time seems strange to us now. 15 hours and 46 minutes after noon is just 03:46 am of the next day. So the observation was made at 03:46 local time in December 28th of 1612. The first drawing shows Jupiter and 3 moons aligned in the same plane, he also puts "distance" measurements between Jupiter and each moon; the one in the left was at 9 Jupiter radii, and the two on the right where at 9 and 10,3 Jupiter radii respectively. He wrote "fixa" for a star in the upper left part of the image which he would take as a fixed point during this set of Jupiter observations.

The lower drawing shows the situation a few hours after that. The moons have moved a little (one of the four moons has just started to emerge from the glare of Jupiter's disk) and "fixa" still remains more or less fixed in place.

Now compare the first drawing with the actual positions calculated using modern ephemerids:

As you can see "fixa" must be the planet Neptune and not a star. The moons marked by Galileo in the first drawing are (from left to right) Ganymede, Europa and Callisto.

Now let's jump one month to these notes:

Here again we see Jupiter and three moons (from left to right: Ganymede, Europa and Callisto at 5,10 , 8,40 and 20,40 Jupiter radii away from Jupiter according to Galileo). This observation was made the 28th of January 1613, 6 hours after sunset. A star labeled with an "a" can be seen 29 Jupiter radii away in the far left of the image. This new fixed point he used was indeed a star. The star is SAO 119234 (or HD 105374). He used these fixed points as references to understand the motion of Jupiter in the sky. Little he knew that the fixed point used a month before, "fixa", was Neptune and therefore not fixed at all. He wanted to use a secondary fixed point for this new observation so he took another star, "b" which was a little farther away than "a" but basically in the same direction. Since "b" couldn't fit in the space left in that page he splitted the drawing and continued the Jupiter-a-b line in the lower right part of the image. Well, it turns out that "b" was "fixa" again that had displaced in the sky (look at the first image of this post). "b" was Neptune. You can see how well it agrees with our ephemerids of that observation here:

The interesting thing here is that in this observation Galileo wrote that "b", which was also observed in the preceding night, and "a" seemed to be more remote from each other (inter se) than the previous night. This means Galileo was able to correctly detect the motion of Neptune, even if he dismissed the idea and never thought of it again. I want to emphasize "correctly detected" because he was aware that Jupiter and "b" would be moving relative to one another but the fact that "b" separated from "a" (both supposedly fixed celestial marks) was really unexpected.

Le Verrier predicted the mass and orbital parameters of the planet (as well as the current position) so he had a very good estimate as to how massive and far away Neptune was. That said it is true that semi-major axis and mass were educated guesstimates.

By the way it is also worth noting that Adams is currently disregarded as a contributing force for the discovery of Neptune. Nearly all the credit goes to Urbain Le Verrier.

Watsisname:

What I'm getting at is if he simply observed Uranus for a few months and tried to observe how much it got pulled by Neptune during that time, then the effect is too small to be detected with the telescopes available in the 1800s.

Let's suppose Uranus and Neptune were as close as they could possibly get to each other (about 10AU apart). The acceleration of Uranus due to Neptune is

which during a time Δt will shift Uranus' position by

If this shift were oriented in the best possible way relative to Earth (perpendicular to our line of sight), then the change in Uranus' position on the sky during this time is

Let's plug in the numbers. Neptune has a mass of M = 1.02x10²⁶ kg. Let the observing time be 6 months. Then the shift in Uranus' position due to Neptune pulling on it during those 6 months, assuming they were only 10 AU apart, is 370 km. From our vantage point about 20 AU away, this can produce, at most, an angular shift of 0.025 arcsec.

Historically, the largest telescopes available in the 1800s were less than 1.8m in aperture. (The biggest at 1.8 m was the Leviathan of Parsontown, which was not complete until 1845, after the irregularities in Uranus' orbit had already been noticed). The minimum angle that can be resolved with a 1.8 m optical telescope is 0.07 arcsec. (Even this 0.07 arcsec is still unobtainable because atmospheric turbulence limits resolution as well, and rarely is better than 0.5 arcseconds. Today we overcome this with either space-based telescopes or adaptive optics.)

So in the 19th century, no telescope in the world could possibly hope to see Uranus be perturbed by Neptune over the course of a few months. There has to be more to the explanation of how Neptune's influence was observable.

In fact, there is. Try checking the wikipedia article on Neptune's discovery for more information on the perturbations.

FastFourierTransform:

That's actually quite interesting and I didn't realized at all how minuscule the effect should be. Even with a simple model you have show how tiny this is. As I see, in reality it would be even worse since that acceleration could be archived only for short periods of time. After and before the actual opposition, Neptune would be farther away and we know that gravitational force is fairly sensitive to distance.

Also in the opposition of Uranus and Neptune there is no place on Earth's orbit to have a reasonable angle to actually see the displacement (produced when the perturbing influence is at max). Since there's no situation where Earth could be perpendicular to the line between Uranus and Neptune (at their closest approach) the angle subtended by Earth's orbit at Uranus distance has to be taken into account. I've calculated the real displacement to be around 0.052 times the one calculated by Watsisname when we incorporate this to the model. So instead of 370 km of displacement, from the best vantage point in Earth's orbit, we would at most see 19 km!! Even with current technology it would be nearly impossible to notice at all (if we had 19 km resolution we would be mapping Uranus' moons from Earth).

I think the answer here has to do with two things (but only one of them seems to be relevant in the end);
First: a larger observational baseline is needed (not just a few months).
Second: Neptune's tug, even if weaker, should be more constant and oriented tangential to Uranus' orbit when they are still far from opposition.

The first thing is easy. Le Verrier used only observations made between 1781 and 1821. That's 40 years of observations. Let's take 8 years for example as a baseline. Uranus and Neptune are close to their opposition all this time so we can still use the approximation of the gravitational tug staying constant. The good thing of the formulas provided by Watsisname is that time is squared so if we increase a little the observational baseline we have a significant increase in the effect. In this case and even considering my 0,052 factor we get 0,34 arcseconds of separation between the real Uranus and the "predicted" Uranus after 8 years. Barely noticeable still.

So the actual perturbing effect is, maybe (I don't know so please tell me if you have better insights) due to the second option. The actual configuration where the Sun-Uranus-Neptune angle is 90º (in this configurations Neptune's tug is directed tangentially to Uranus' orbit) was archived around the year 1800 and again around 1850 (in the middle there was the opposition). In this configuration Uranus and Neptune were separated by a distance of r = 25,4 AU. But also in this configuration we (the Earth) can be at a good angle to see Uranus' displacement (perpendicular to the Neptune-Uranus line if you want), so here there's no need for my 5,2% reduction of the effect. You can admire the 100% of the displacement in this configuration. The problem is that the tug is weaker since Uranus and Neptune are farther apart than when they are at opposition. But still considering just 6 years it would mean a 0,59 arcseconds discrepancy (surpassing 1 arcsecond after 8 years only), something that starts to be noticeable I guess.

Here we have a graph (taken from Le Verrier: Magnificent and Detestable Astronomer) that shows the "residuals" of Uranus' longitude across 150 years or so.

O stands for the observed longitude of Uranus and C stands for the computed longitude of Uranus after taking into account the perturbations of all the planets. If there was no Neptune this should have been a straight horizontal line at 0 arcseconds, since the computed position of Uranus and the observed one should agree, so the residual (which is the difference between both values) should be zero. But Neptune exists and the residual fluctuates.

Now lets see the actual configuration of the planets through the XIX century (this has been taken from "The Case of the Pilfered Planet" an article published in the Scientific American journal)

So lets compare both images. Around 1800 and 1850 we have the configuration I mentioned (where Neptune's gravitational pull seems to be tangent to Uranus' orbit). Around 1820 we have the opposition. As you can see Uranus' predicted position is lagging behind with respect to Uranus' actual position because Neptune is accelerating Uranus along his track (his orbit). The opposite can be seen in 1850 when the predicted position of Uranus is ahead of the actual position (because Neptune is pulling Uranus in the opposite way so it decelerates). This is also evident in the residuals graph; in 1800 the computed longitude for Uranus was larger than the observed one so the residuals where negative, but close to 1850 the residuals are positive since the observed longitude is larger than the computed one.

Also interesting to note the fact that at opposition the effect is the lowest possible, even if paradoxically the gravitational tug is the greatest, because of the proyection effect of the acceleration I've talked about in my second statement. You can see that the residuals disappear at opposition and Uranus' position does not differ from the predicted one (also, bare in mind that the separation in the second image is greatly exaggerated to make it clearer).

But I still miss a more in depth explanation. I've shown that in 1800 and 1850 the effect is more noticeable than in 1820 when the planets are closest to each other. But still we get around 1 arcsecond displacements after 8 years. Too little in my opinion. Also in the residuals graph they are using one or two orders of magnitude above that (tens or even near a thousand arcseconds discrepancy) for the residuals so probably I am completely lost in terms of what the residuals actually are and how they are actually calculated here. I can only make little sense of the overall trend in the different years. Also I can't understand why the graph of the residuals behaves like that, it seems to increase a lot when the planets are farther apart for some reason and it is not sinusoidal as I intuitively expected. Also the graph seems to indicate that the zero point for the residuals (where the difference between observation and prediction is zero because the prediction was made from this moment) was taken around 1775, right? Why? And why it is so close to the expected time on which Uranus and Neptune should be the farthest apart possible? Coincidence? Those are the kind of things that make me think that I haven't understood the basics of the ideas behind the discovery. I would like to know more. a lot more. It would be crude but beautiful to see the full demonstration just to feel the power of Newtonian vision.

But sadly this is a difficult topic. Jean-Babtiste Biot attempted to explain Le Verrier's methods in six papers he published in 1846 (just before the discovery) and in the third paper he wrote "As I progress in the task I have undertaken, the difficulty of the subject seems to increase". I feel the same. The more you know the more complex it seems to be. I hope there's someone here that has the patience to explain it in detail because it is very exciting subject in my opinion and there should be some place in the Internet that has an explanation like that.

SOME NEWS: I just found a web where there is an interesting explanation of Le Verrier's calculations but it is in French and I'm unable to understand it. If some French fellow forum members want to read it please show the rest of us what is all about

Watsisname

It's amazing what we can do with computers today. I did try directly computing orbits with and without perturbation just to try to get some intuition for how the perturbation effect works. For starters I assumed co-planar, circular orbits, and also that only Neptune perturbs Uranus and not vice versa (which probably matters very little). For the initial condition I put Uranus directly opposite Neptune, and then iterated through 1 hour time steps. I add vertical lines to show when they are closest together (by direct distance, not angle):

I think this region (zoomed in between a conjunction and following opposition) fits with Le Verrier's figure:

There are still significant differences between them (like the amplitude, and width of the minimum vs. maximum), which I guess is mostly due to Uranus' eccentricity. I must factor that in to 'perturbation simulator v0.2'. The shape and overall trend of the curve also depends a lot on the initial position of Uranus relative to Neptune. I have no idea what Le Verrier used in detail. But any rate I think the physics operating here is pretty clear like you said. The variation occurs because while Uranus is approaching Neptune, Neptune causes it to accelerate, and then decelerate after closest approach. It's a very slow effect since the gravitational attraction is weak when they're far apart, but it has a lot of time to accumulate.

I think this is also exactly the explanation behind exoplanet transit timing variations. We see the projection of the planet in its orbit against its star, so these same perturbations will cause transits to appear earlier or later than they otherwise would if there were not additional planets in the system. I'm not sure of the exact details for how we go from the timing observations to predicting where the other planets are, but it must be essentially the same problem, solvable by computer with Monte Carlo methods.

I'm still very curious in how exactly they solved Neptune's orbit without computers. For a computer this is a problem that takes a few seconds. But by hand...? I can't even imagine. Plus they were observing the effect, and had to deduce the cause, by successively adjusting the "guess" for the location of Neptune to get closer and closer to what was observed. This is more in the realm of "mathematical methods of physics", but a very interesting subject.

Author: FastFourierTransform
Original post on the forum: link

I want to try to address this. As Watsisname said; establishing a non-inertial co-rotating frame of reference has many advantages in terms of understanding the actual situation but it also has some disadvantages (non-inertial forces as Coriolis and centrifugal have to be invoked and even then are difficult or impossible to represent in a potential dwell map). Inertial frames of reference are great because of one thing: they are what we are used to picture in our minds frequently and they are the ones where physical laws are coherent without any extra tweaking needed (no fictitious forces needed). The explanation from the point of view of an inertial observer should be equivalent to Watsisname's but it should be given in other terms. I also feel that the actual inertial perspective on the subject should be extremely complex to wrap our heads around without doing many alienating calculations.

I haven't been able to come with a simple discourse for the stability of L4 and L5 nor for Trojan orbits, but I think the fact that Lagrange points exist as "dynamical" equilibrium points is feasible within an inertial framework, so here it goes:
Since we are not co-rotating is useful to remember that Lagrange points are not fixed in space and move themselves. Also very important to remember is Kepler's third law, which here will be expressed by the fact that objects farther from the Sun have lower orbital speeds (the angular speed will be dependent as the -3/2 power of the distance to the Sun). We place ourselves at rest relative to the Sun-Earth barycenter.

So. Lets say we want our probe to orbit the Sun but always been at the same distance from Earth without any additional effort, profiting from the gravitational landscape of the Solar System to accomplish that. If the Earth was a massless object then any point in Earth's orbit would be a suitable place to put our probe since it would travel just as Earth does but leading the march or lagging always the same constant amount. But the Earth has mass and some influence so your probe would be pulled back or ahead by the planet and it would be unable to maintain its heliocentric Earth-like orbit. This is interesting because Watsisname did the reverse, he started with a huge mass and then slowly went to the small Earth mass to show how nearly the entire orbit of Earth was covered with a low-inclination plain potential field. If we went a step further, to a massless Earth, we indeed would have equilibrium points all across Earth's orbit.

L1 point

So the L1 point is the easiest to understand. But first you have to abandon the misconception that the gravitational tug of the Sun equals that of Earth in that point. This is not at all about that (in fact as we saw in previous posts, even the Moon is farther from Earth than the equal-force point). So no, L1 is not the place where the Earth screens the Sun's gravity and you can move freely as if there was no Sun nor Earth until they move and lose that particular vantage point, no! L1 is located closer to the Sun than Earth so an object at that distance in a circular heliocentric orbit should get ahead of the Earth since it moves faster in its orbit, but the Earth here is creating his magic; having the Earth right in the opposite direction of the Sun what your probe feels is that the Sun-Earth system is equivalent to a Sun without Earth but this Sun having a little less gravitational attraction. So the Earth is orbiting the Sun but our probe orbits as if it orbited a lighter Sun, this means that even if it is closer to it our probe can still move as slow as Earth due to the lesser influence of this lighter Sun. So the the probe moves at the same angular velocity around the Sun as the Earth because it feels the Sun less and it is closer (both effects balance out) so the L1 - Sun - Earth configuration is maintained constant. Bare in mind that the vector sum of the gravitational force of Earth and that of the Sun on an object placed here at L1 will have the same direction as to point to the Sun-Earth barycenter (so in fact L1 is orbiting the Sun-Earth barycenter also).

L2 point

L2 obeys a similar interaction. Instead of the tug of Earth screening a small part of the Sun's influence we get the Earth adding to the pull of the Sun (since L2 the Sun and the Earth are both in the same direction). A probe in L2 would, in principle, lag behind Earth because of Kepler's third law since it is farther from the Sun, but at the same time the Sun "feels stronger" (thanks to the extra pull of Earth) so you end up going faster than expected. Both effects balance out in the L2 point. And since you move at the same rate as the Earth here (same period or same angular velocity) you will always maintain the same configuration with respect to the Sun and to Earth and that in turn means that you will maintain this behaviour and the L2 point wouldn't go away from you (you will move with it). Again, as in L1, bare in mind that the vector sum of the gravitational force of Earth and that of the Sun on an object placed here (at L2) will have the same direction as to point to the Sun-Earth barycenter.

L3 point

L3 is similar to L2 but symmetrical to the Sun. At the opposite side of Earth's orbit you notice the pull of the Sun as if Earth nearly didn't existed (you are ~2 AU from it) but Earth still contributes to make you feel the combined pulls of the Earth and the Sun in the same direction as if the Sun had gained some mass and was a little stronger than the Sun we feel from Earth. This means that you can orbit this "more massive Sun" quicker if you are in Earth's orbit but if you get a little farther from the Sun than Earth's orbit you would reach a state of "dynamical equilibrium" again were you are moving at the same speed as Earth does. This is interesting because we often see pictures of L3 been located in Earth's orbital path but in fact it is outside of it, like L2. The difference with L2 is that for you to beat the angular speed increase of the strong gravitational pull of the Sun and Earth you (now in L3) just need to get a little farther but since in L2 the Earth is quite close the influence that has to be overcome is greater and thus you move considerably farther from Earth's orbit than in the case of L3 to archive the same balance. Again, as in L1 and L2, bare in mind that the vector sum of the gravitational force of Earth and that of the Sun on an object placed here (at L3) will have the same direction as to point to the Sun-Earth barycenter.

L4 and L5 points

Since the gravitational field of the Sun-Earth system is mirror-symmetric to the Sun-Earth line then both L4 and L5 points have the exact same underlying explanation. Remember that L1, L2 and L3 all were orbiting the Sun-Earth barycenter because the net force (the result of adding Earth's tug and the Sun's tug) always pointed there? Well this fact was simple to grasp because in the 3 cases we had only vectors in the Sun-Earth line (so the net result would always lie in the same line as the Sun's and Earth's pull). But what if we could have a net force still pointed to the Sun-Earth barycenter but with the Earth's pull and the Sun's not been parallel from our position? Well, that's exactly the L4 and L5 points. Here you will experience the Sun pulling in one direction strongly and the Earth pulling in an other (not parallel to that of the Sun) more weakly but enough to have your net force dragging you to the Sun-Earth barycenter instead of just the Sun. Since this geometrical situation can be archived in many places we have to also include the restriction that you maintain the same orbital speed as Earth (if not you would feel the net force pointing to the Sun-Earth barycenter just for a brief moment and then you would be pulled slowly to a different direction as the distance to the Sun or to the Earth would vary). With both restrictions (having your net force pointing to the barycenter and moving at the same orbital period as Earth as to maintain that net force pointing to the barycenter) you end up with only two places in the entire system, and those are the L4 and L5 positions.

L1 stability

This is difficult to address in an inertial frame but let's try. To understand if L1 is a stable or unstable equilibrium point we are going to imagine our probe placed in 4 different places (all very close to L1). Let's take a slightly higher orbit from L1 (this means that we are in a Sun - system barycenter - L1 - our probe - Earth - L2 configuration). In this situation we have a greater pull from Earth so our net force will be pointing the same way (to the barycenter) but with less intensity (since here the Earth is now screening more of the Sun's influence). Been in a higher orbit alone would mean we are moving slower due to Kepler's 3th law but also we move slower because of the "lighter Sun" we have now. Because of that our probe would lag behind Earth. As this situation continues to evolve the net force would start to have a component pointing outside of the barycenter of the system (as the angle between the pull from the Sun and the pull from Earth gets lower than 180º). This in principle should be nice because it should impart some tangential motion on our probe so that it reaches L1 again, but sadly we are in a to high orbit for this to work. Our probe will gain a little acceleration but would end of still losing the race to Earth. But! this could happen if instead of moving to a higher orbit from L1 we maintain the same orbit as L1 but place our probe behind L1 as it moves. Then the orbit would have a smaller circumference than in the previous case and this new tangential pull would be enough to make our probe reach L1. Once the probe reached L1 it would probably have to much speed to maintain that position so it would come ahead of L1 and the reverse would hold. The probe would be oscillating like in a valley in this case. Lastly let's think about the Sun - system barycenter - our probe - L1 - Earth - L2 configuration (our probe is in a closer orbit than L1). In this situation we are faster than Earth since the pull of Earth is not enough (since we are farther from it) to overcome the effect of a smaller quicker orbit around the Sun. The probe would surpass Earth and it would end pulling us backwards, but not enough strongly. As the tangential acceleration performs its magic the probe would slow down and lower the orbit even farther from L1 making the probe lose its place of near equilibrium.
What does all of this means? It means that any displacement from L1 performed in the Sun-Earth line would end up been amplified and any deviation in the perpendicular direction of that line will be corrected and the equilibrium restored. So here we have for L1 a saddle point as was expected (stable equilibrium in one direction but unstable in other direction).

L2 and L3 stability

Same thing as before but here the situation is more simple. Any small deviation from L2 in the Sun-Earth line can be explained in the same way for higher or lower orbits from L2; we are now considering a "heavier Sun" (due to Earth's contributing to the net pull in the same direction as the Sun). If we are in a lower orbit (with respect to L2) the probe will move faster and will overcome the Earth, but then a tangential force would try to slow down the probe because Earth would be a little behind it. If it slows down the orbit shrinks again so a positive feedback loop is generated where the probe ends been separated more and more from L2. For a higher orbit the the probe would lag behind Earth so a tangential pull would appear that would speed up the probe throwing it to an even higher orbit and escaping L2 more and more. So L2 is an unstable equilibrium point for any direction you choose to displace your probe from it. L3 behaves exactly in the same manner due to the same reasoning but since it is farther away from Earth the instability is less abrupt and we still have a lot of margin to move the spacecraft around L3 without loosing control so quickly (this is a peak with a less inclined slope).

L4 and L5 stability:

Here is where I get lost. The geometrical configuration here is way more complex than in the other cases since it ceases to be a 1-dimensional problem and becomes 2D. I don't know exactly hot to picture this in my mind but any small deviation from L4 has to be diminished since as Watsisname explained this is a stable equilibrium. My guess is that this problem is equivalent to the problem where we only have the Sun and the probe, but this Sun changes mass slightly and harmonically and wobbles from one place to the other in such a way that it pulls the probe as if it was Earth. An animation should make this more visually clear. Let's say, for example, that we are ahead of L4 (still behind Earth). In this situation Earth's pull will increase in magnitude and the direction of the net force will change slightly (from the barycenter towards the Sun). The speed then would increase so the probe would archive a higher orbit and slow down, but this somehow would happen in such a way that L4 would surpass our probe because of that slowdown and we would end up behind L4 quickly, where the reverse would hold. We would be somehow oscillating around L4. How is this the case for displacements from L4 in any direction I don't really know but sure an animated diagram with the forces involved would clarify this a lot (maybe one day I try to do it).

Author: Watsisname
Original post on the forum: link

Introduction

When two masses move in a circular orbit about each other, five points of gravitational equilibrium are generated, which move along with them.  A particle placed at one of these points will have a very small acceleration, as if gravity was in some way cancelled out.  More strangely still, it can be possible to place a particle in a stable orbit around some of these points, despite there being no mass there to attract it!  

These five special locations are known as the Lagrange Points, discovered by Leonhard Euler and Joseph-Louis Lagrange in the 18th century.  Not only are they of special interest in celestial mechanics, they also have important applications for space missions and the placing of satellites.  We even find asteroids caught around some of the Lagrange points: so called "Trojan asteroids".

Why do these points exist?  To investigate, we will again use the concept of a potential, as we saw in the discussion of Mercury's orbital precession.  As a reminder, in physics a "potential" is an imagined landscape of hills and valleys that describes how an object reacts to forces.  We can think of an orbit as being the motion of a particle trapped in a valley of a potential.  This time however, we will generate the potential in a somewhat unusual way: we will adopt a frame of reference which is rotating.  There are pros and cons to this approach.  One one hand, motion in a rotating frame of reference is not very intuitive, and we will need to invoke "non-inertial" forces to understand them.  On the other hand, doing it this way will make life a lot simpler, as then the masses and the Lagrange points themselves will appear stationary.  Motions around the Lagrange points will also be much easier to visualize and interpret.

But, seeing as how much I myself struggled when I was learning about non-inertial forces in rotating frames of reference, I think a quick digression into that subject would be a good idea.

Frames of Reference and Non-Inertial Forces

Most of us are comfortable with physics in "inertial" frames.  An inertial frame is one in which objects subjected to no forces will move in straight lines with constant speed due to inertia.  In other words, an inertial frame is one in which the "law of inertia", or Newton's First law, is valid.  But what if we put ourselves on a turn table, and try rolling a ball back and forth to each another?  We will discover it is quite hard -- the ball will seem to be deflected by some mysterious force.  

Someone watching from the side will still say the ball moves in a straight line relative to the ground.  It is us on the spinning turn-table that accounts for the weird trajectory we see.  However, we can be just as valid in saying that the ball is deflected by something.  We may invoke new forces, unique to our rotating frame.  There is a centrifugal force which tends to fling things away from the axis of rotation, and a Coriolis force which tends to make moving objects turn and follow curling paths.  We all have felt centrifugal force whenever riding in a car going around a turn.  And many are probably familiar with how Coriolis forces on the spinning Earth causes air currents to deflect, and explains the rotation of hurricanes.  These forces are sometimes said to be "fictitious", in the sense that there is no real physical interaction responsible for them.  They are simply consequences of being in a strange frame of reference.  But when we are in that reference frame, they seem totally real.

I have found no more wonderful a demonstration and explanation of these fictitious or non-inertial forces than in this film put together by University of Toronto physicists in 1960.  If you have not seen this, it is an absolute must watch.  Or if you're in a hurry:

-The first 13 minutes build the ideas of different frames of reference.
-13:25 to 17:00 looks at non-inertial (accelerated) frames.  
-17:04 is a beautiful example of motion in a rotating frame.  This segment is extremely relevant to the physics we will be looking at here.

The Generalized Potential for Two Masses

Here we go!  Let's imagine two masses in orbit around each other, and examine the forces and the potential they produce.  For starters, and for symmetry, let's first let the two masses be equal.  We will also "lock them in place", so that (for now) there are no non-inertial forces.  What does the potential look like for this system?

For any point in space there will be a gravitational force.  Actually two of them: one from each of the masses.  The gravitational force from each is given by

where r with a hat over it means "in the direction from the mass to the point in question", since force is a vector, like an arrow, which has magnitude and direction.  The minus sign in front flips it around so that it points to the mass.

In terms of the potential (the landscape whose steepness describes the force), this becomes

This expression lacks any vectors, since each point is assigned only a magnitude (a height in the landscape).  What does this landscape look like?  Below I plot it for the two equal masses, 1 unit apart along the x axis, fixed in place and centered around the origin.  Darker/redder colors indicate more negative values of the potential.  Around each mass is a bottomless hole.  Drop a ball into this landscape, and it will fall into the hole associated with the nearest mass.

As FFT has said earlier, the steepness of the potential (how closely packed the contour lines are) relates to the magnitude of the force at that point.  On a topographic map, close lines means steep terrain.  A ball on steeper terrain will accelerate down the hill faster, and likewise a particle on a steep region of a potential will experience greater force.

To help see how the potential and the force/acceleration are connected, below I'll plot the exact same situation in terms of the acceleration. I'll plot on a log scale, so that details at very small accelerations are visible (whereas on a linear plot they would be overwhelmed by the very large accelerations close to the masses).  Here, darker blues represent greater acceleration.

Unsurprisingly, near the origin the acceleration drops to zero.  There the gravitational force due to the mass at the left exactly cancels the force due to the mass on the right.  We can also see a hint of this from the potential.  This location is a "saddle point", where the landscape rises up if we move along the y axis, but drops down if we move along the x axis.  Trying to balance an object here would be very difficult, with a tendency for it to "roll off".

The Co-Rotating Frame: how do non-inertial forces modify the potential?

Now let's set these two equal masses into a circular orbit around each other, and place ourselves into a rotating frame of reference, such that the masses appear to stay still.  What new forces appear, and how do they affect the potential and the motions?

The first new force this will introduce is centrifugal.  This force acts to fling things away from the center of rotation of the system, about which our frame is rotating with angular velocity Ω (meaning the speed at which the angle is changing).  We can find the angular velocity of the frame by using Kepler's Third Law:

where R is the separation distance between the two masses M₁ and M_2.  With the frame rotating at this angular speed, the centrifugal force is

where r is the distance from the point of interest to the center of mass of the system about which the frame is rotating.  The r-hat in this equation tells us the centrifugal force is directed away from the axis of rotation, as it should be.  The strength of the force also grows in proportion to the distance from the center of the frame, and also with the square of the angular speed that the frame is rotating.  

In terms of the potential, the centrifugal force adds a new term to it, with magnitude

Let's see what this looks like now, again for the two equal masses.  Now the masses are orbiting in a circle, and we are in a frame which rotates with them.

There they are!  Five points where the net force is zero!  One is in the exact center as before, where the gravitational forces cancel and the centrifugal force is zero since it is at the center of rotation.  This is the L1 Lagrange point, in the case of the two masses being equal.  

Two new points are added on either side of the masses.  These are L2 and L3.  Here the centrifugal force is exactly equal and opposite to the sum of the gravitational forces.  In the potential they also look like saddles, and are unstable.

Finally, two new points are added above and below.  These are L4 and L5, where again the centrifugal and combined gravitational forces cancel.

You might notice that the inclusion of centrifugal force makes the potential drop down again at large distances, whereas with gravitation alone, it rose upward.  This is because if you start out far away from the center, and match your speed to rotate along with the system, then you are moving quite fast, and have far too much angular momentum to be in a circular orbit.  Your side-ways speed will carry you off to larger distances.  On the other hand, if you start out too close, then you have too little speed to be in orbit and so you drop closer in.

Let's do something fun!  Earlier I said that L4 and L5 are stable points about which we can place something in orbit.  So let me just drop in a particle very close to the L4 point, and....

Wait, what?  Why did it spiral away?

Well, I had lied a little bit.  L4 and L5 are not always stable.  And you might have noticed that the potential at these points looks like the top of a hill.  Not the bottom of a valley! (Lighter colors instead of darker).  How could there possibly be stable orbits around them?  (Actually I did not understand this until now.  I had suffered the misconception that L4 and L5 are shaped like valleys, until I started working through this in detail!)

What is going on?

L4 and L5 are the tops of hills, not the bottoms of valleys.  So of course an object placed near them will slide away.  But once it does, a new non-inertial or "fictitious" force takes hold: 

The Coriolis Force:  

This is the force which causes the deflection of the puck seen in the demonstration in the above video with the turn table.  It is a velocity dependent force, and it always acts at right angles to the velocity.  This means it will cause a moving object to turn, but not change speed.  The Coriolis force is given by

The cross (X) product between Ω and v means that the direction of the force is given by the right hand rule, and then we flip the direction 180° because of the minus sign in front.  So if the frame of reference is rotating counter-clockwise, then the force on a particle moving upward will be to the right.  That is, it will deflect trajectories clockwise.  This is also why high pressure systems (air flow away from center) spin clockwise in the northern hemisphere, and why low pressure systems (air flow toward the center) spin counter-clockwise.

Because the Coriolis force is velocity dependent, it is impossible to define a potential for it.  If we move faster, the force will be stronger, and we'll take a different path, with a different change in potential energy.  So there is not a unique landscape that we can use to model motions subject to this force.  Instead, we can model the general potential with just gravitation and centrifugal force, and to simulate motions we can include the effect of Coriolis force by direct computation on the computer.  This is how I generated the above plot.

In the case of the masses of M₁ and M₂ being equal, the particle dropped near L4 slides down the potential hill, and then as it speeds up, the Coriolis force whisks it around in a growing spiral.  There are no stable orbits around any of the Lagrange points in this scenario.  

In order for L4 and L5 to be stable, one of the masses must be at least 24.96 times the mass of the other.  Understanding why this is the case (or why this number in particular) is hard, and requires a careful analysis of the equations of motion.  If any readers are interested in the mathematical details, I can refer you to some suggested reading at the end.  But for our purposes, let's just take it as an established fact.  What we can do is look at more examples of the potential where we vary the mass ratio, check what happens visually, and see if we can qualitatively make sense of it.

The generalized potential for a 10:1 mass ratio:

Let's set the mass of one object to be 1/10th the mass of the other.

The acceleration plot is particularly illuminating.  The Lagrange points are starting to look a lot like their commonly-shown configuration, with L4 and L5 leading/lagging the secondary mass by 60°.  This 10:1 ratio is often used in figures that portray the Lagrange points for the Sun and Earth.  Which is an exaggeration of truth, because Earth is of course much less than 1/10th the mass of the Sun!  But it strikes a nice balance, giving an idea of the shape of the potential for the Sun and planet.  Using the true mass ratio of Sun and Earth would be quite unhelpful for visualizing it, as the Sun's influence dominates the landscape.

While we're at the 10:1 ratio, let's drop in another particle near L4:

Once again it spirals away.  10:1 is still too small of a ratio of masses, when we need ~25:1 at least.  But we can see that we're going in the right direction!  The spiral is a bit tighter, and the particle wraps around the L4 point completely before dropping down and hitting the smaller mass.  (Oops.)  

Finally, let's increase the ratio to match the Sun and Jupiter.  That is 1047:1 -- definitely enough for L4 and L5 to be stable.  And since this is a real situation with specific masses and distances, I'll add in units.  The potential is in terms of energy, here shown in gigajoules per kilogram.  Distances are in AU.

It's pretty hard to see the details here (again why most figures online show 10:1 instead of the real thing).  Even the log-scale acceleration plot isn't terribly useful.  But we can at least see that the potential in this case has a crescent ridge along Jupiter's orbit, with Jupiter making a small depression around itself.  The L1, L2, and L3 points still exist and are saddle-shaped, but they are not visible in this plot scale.  Similarly, L4 and L5 are still separate hills along the crescent, but we can't distinguish them from one another.

Let's zoom in to the area around L4, and drop in a particle again.

Here I said "WOW!"  I had no idea that the orbits around Lagrange points look like this.  But they do -- spiraling about in these crazy elongated paths.  I was extremely surprised and excited when the computer output this trajectory.  These orbits are known as "Trojan orbits".  As mentioned before, they have enormous astrophysical significance, and in fact many asteroids are caught in them, especially around Jupiter's L4 and L5 points.  One has also been found in a Trojan orbit of Earth!

Here is an even larger one (starting farther away from the L4 point), with an overview including the Sun and Jupiter but without the mess of the potential.  

As a sanity check, I went to the literature. Check out Figure 3 from this paper "Horseshoe and Trojan orbits associated with Jupiter and Saturn", by Everhart, 1973.  A very close match!

How can we explain this behaviour, qualitatively?

As before, when the particle slides down the potential hill from the Lagrange point, Coriolis force whisks it aside, turning it clockwise.  But soon it is turned so much that it rises back up the potential hill, at a location farther along the orbit -- either closer to or farther from Jupiter.  As it rises up the hill it slows down again, the Coriolis force weakens, and the arc shrinks.  Then the process repeats -- falling down the well, getting deflected to the side by Coriolis, and so on.

But the shape of this potential hill is a like a very broad crescent.  Farther away from the Lagrange point, it gets thinner, and a little less high.  This means the particle has an easier and easier time rising back up to the top of the "ridge".  Finally it reaches a point where it is able to cross over the ridge, and the process repeats on the other side!  It is the Coriolis force that keeps the particle circulating around near the top of this potential hill -- with the hill itself being generated by gravity and centrifugal force.

For larger deviations from the L4/L5 points, the particle will spiral along more and more elongated paths.  Once the departure gets too large, it can even cross over to the other side of the Sun-Jupiter line.  Now we have Horseshoe orbits!   Some asteroids are known to be caught in these orbits as well.

One other thing to add real quick, and then this will be the last for a while.  I wanted to show that it is also possible to orbit around the unstable Lagrange points (L1 for example), but the orbit is extremely delicate and will fall apart if not carefully maintained.

Here it is For Jupiter's L1 point.  The particle starts out a little bit closer to Jupiter than the L1 point, and with an initial velocity downward.  Coriolis force deflects its path clockwise, and with the initial position and velocity tuned "just right", it ends up circulating in an oval.  For the values I chose, it completed about an orbit and a half before falling apart.

To give an idea of the timing involved, this path is plotted over 12.5 years, and the first orbit took about 5.5 years.  Plenty of time for a spacecraft trying to follow this course to make necessary adjustments to stay on it.  This has actually been done with missions to Earth's L1 and L2 points.  Here is the path the LISA Pathfinder used to get to and orbit around the Earth-Sun L1 point.  Again some maintenance is required to follow this orbit.

Concluding Remarks & Suggested Reading

I think I will end it here, as I think this has covered all (and more!) that I had hoped to go through.  I myself have learned an enormous amount while working through this, and came across many surprises and previous misconceptions.  I did not know for example that L4 and L5 are actually repulsive, and that the stability of orbits around them depends on the Coriolis force.  Nor did I really have a sense of what Trojan orbits really look like or why they behave as they do.  This was extremely exciting to work through and discover, and I hope you all find this an interesting and enjoyable exercise as I did.

Cheers!

References & Further reading:

An excellent, more technical treatment of Lagrange points, including deriving them and calculating their stability, can be found here. (written by Neil Cornish of Montana State University).

Another good one can be found here from Tipler's Modern Physics text.

For those interested in a broader scale of classical or celestial mechanics, including deriving equations of motion in non-inertial frames of reference, I can highly recommend the text "Classical Mechanics" by John Taylor.

Background

Often in discussions of Mercury's perihelion advance, the effect is shown greatly exaggerated.  Sometimes even the orbit is shown much more eccentric than it really is.  Before we dive in, let's take a moment to get some perspective on the true shape of the orbit and appreciate the scale of the solar system.

 
All orbits are ellipses, but most planetary orbits are close enough to circular that it can be hard to tell at a glance.  Mercury is the biggest exception with an eccentricity of 0.21.  (Or Pluto at 0.25).  What's interesting, and what this post is all about, is that Mercury's orbit is not stationary.  The ellipse slowly shifts around, it's perihelion point advancing by about 574 arcseconds (or 0.159 degrees) per century.  Much of this (531 arcseconds) can be explained by the perturbations from the other planets.  However, the remaining 43 arcseconds per century are "anomalous", unaccounted for by Newtonian mechanics.  

The puzzle of this additional observed precession was first noted by LeVerrier in 1859.  Following the successful prediction of the existence and location of Neptune from similar perturbations on Uranus' orbit, LeVerrier proposed that Mercury's orbital precession could be caused by another undiscovered planet inside of Mercury's orbit, which he named Vulcan.  Astronomers searched for this missing inner planet during solar eclipses.  Of course, Vulcan was never found... 

The mystery was not solved until Einstein began developing his general theory of relativity.  When he applied his general relativistic equations to the problem in 1915, he found that they exactly predicted the additional 43 arcseconds per century.  He considered this one of his greatest achievements.

Imagine my joy at the recognition of the feasibility of general covariance and at the result that the equations correctly yield the perihelion motion of Mercury. I was beside myself for several days in joyous excitement.

Physics

So, why does it happen? Is there an extra energy involved?  (Spoiler: Yes, in a way!)

For an elliptical orbit, we can think of the orbital motion as being made of two parts: an oscillation in azimuth (angle around the Sun), and an oscillation radially (in and out).  In Newtonian mechanics, the periods of these two oscillations are exactly equal.  That is, each time the planet completes one cycle around the Sun, it also exactly completes one cycle radially in and out.  Therefore in Newtonian mechanics there is no precession of the orbit (besides that which is caused by the influences of the other planets).

In general relativity this is no longer true.  The periods of the two oscillations are different!  To see why, it is easiest to examine the shape of the "effective potential", which you can think of as being just like the potential energy well around a spherical mass, except also taking into account the angular momentum of the orbiting body.  

If you are unfamiliar with the concept of a potential energy well, imagine a rolling landscape of hills and valleys.  If you place a ball at the bottom of a valley, it will just sit there.  It's in a stable equilibrium.  On the other hand, the top of a hill is an unstable equilibrium.  If you displace the ball slightly from the top of the hill and let it go, it will continue to roll down, exchanging gravitational potential energy for kinetic energy as it drops in height.  Finally, if you have some bowl-shaped valley and drop the ball from rest at some initial height, then it will roll down to the bottom, and (ignoring friction) keep going, rolling back up the other side until it reaches its former height, and then roll back again, oscillating back and forth indefinitely.  

Because the gravitational potential energy near Earth's surface is simply proportional to your altitude, if you make a plot of the potential energy of a ball in this landscape as a function of its position, it will be exactly the shape of the landscape.  This form of modelling of the potential is extremely useful.  In general, we can use "potential energy wells" as a way to understand motions of objects subjected to different kinds of attractive and repulsive forces, from balls rolling down hills to the vibrations of atoms bound together in molecules.

An elliptical orbit is an oscillation in an effective potential well, as well.  There is a hill to climb as you move outward due to the gravitational attraction of the star, and there is also a hill to climb as you move inward!  This is because the sideways motion is associated with angular momentum, and angular momentum acts like a repulsive centrifugal force.  By conservation of angular momentum, the sideways velocity increases as your orbit swings inward, increasing the centrifugal repulsion and making it harder to approach the center.  If the angular momentum is large enough, then the inward fall can be halted, and you swing back out again, making an orbit!  

So if we can compute the effective potential well for a planet orbiting the Sun, we can figure out some things about its motion.  The anomalous precession of an orbit can also be understood this way, by comparing the motion from the Newtonian potential with what happens in the general relativistic version of the potential.  Let's try it!

Math

The goal here is to try to understand how things work conceptually, and I will try to emphasize those concepts (especially physical concepts) as we go along.  The hard truth is that much of it will still require going through the math in order to access it, but I'll do my best to turn that math into something visualizable and comprehensible.

In Newtonian mechanics, the gravitational potential energy of a mass m a distance r around a spherical mass M is

This describes a simple curve which plummets downward indefinitely as you move to smaller radii.  As it should.  If you drop something near a massive object and give it no sideways motion (no angular momentum), then it will fall straight into it.

If we instead give the object some angular momentum L, then the potential is modified:

where μ is the reduced mass, which for a planet orbiting the Sun (m << M), reduces to approximately m.

Let's see what this looks like for Mercury:

Again the way to think about this is that the potential describes a landscape that a ball will frictionlessly roll across.  To visualize that, I added the large red dot to represent Mercury's position as a function of time.  The acceleration is determined by the slope of the potential, and I plotted the motion with 2 days per frame.  With 44 frames in all, this traces out one complete period of Mercury's orbit in the radial motion.  I also added a horizontal line to show the total energy of Mercury (its gravitational potential energy + kinetic), which is constant along the orbit.  Where this line is above the effective potential defines the range of distances from the Sun that Mercury's orbit will cover.  The vertical lines represent the extremes (Mercury's perihelion and aphelion distances).

Why is the Newtonian effective potential shaped this way?  
Because of Mercury's amount of angular momentum, the shape of the effective potential that it sees near the Sun is a valley with hills on either side.  The hill at large radii is due to the Sun's gravitational attraction, while the hill at small radii is due to centrifugal repulsion.  Even though the gravitational force grows stronger at closer distances, for an orbit the centrifugal force gets stronger more quickly, due to the conservation of angular momentum which increases your sideways speed.

A useful concept to keep in mind here is that the period of the radial oscillation depends on the curvature of the potential well.  If the well opens up more sharply, then the average acceleration is greater, and the period is shorter. 

And here's one other useful trick.  If the amplitude of the oscillation is not too large (does not reach too far away from the minimum of the well), then we can approximate the well as a parabola around the minimum.  Then for a parabolic well the motion is described very simply as a simple harmonic oscillator.  For Mercury this approximation isn't particularly good (the well is noticeably asymmetric over the region Mercury covers), but it's not terrible either.  

For a simple harmonic oscillator, the frequency is given by

which for the radial motion leads to

The frequency for the angular oscillation on the other hand is given by

Which leads to the exact same expression:

So we see in Newtonian mechanics the periods are exactly equal and there is no anomalous precession.  Now let's see how this gets modified when we move to General Relativity.

The Effective Potential in General Relativity

Mercury's Orbit in General Relativity

Going Further: The Bizarre Orbits in Strongly Curved Space-Time

Don't miss discussion of the original post on the forum: it has a lot of interesting questions and answers.