Gravitational Potential and Newton's Law

The Big Picture

Drop an apple. Why does it fall straight down? Why does the moon stay in orbit instead of flying off into space? Why does light bend around the sun? Three questions, one answer at the Newtonian level: a single scalar field called the gravitational potential $\Phi(\mathbf{r})$ encodes the entire effect of every mass in the universe. Just like the electric potential $V$ of the previous section, it satisfies a Poisson equation:

Here $\rho(\mathbf{r})$ is the mass density at the point $\mathbf{r}$, $G \approx 6.674 \times 10^{-11}\,\mathrm{N\,m^2/kg^2}$ is Newton's gravitational constant, and $\nabla^2$ is once again the Laplacian. The equation says the same thing in gravity-speak that it said in electrostatic-speak: mass density tells the potential how to curve; the gradient of the potential is the acceleration any test particle feels.

Why this is exhilarating. Newton wrote down $F = -GMm/r^2$ for two point masses. That equation, paired with the principle of superposition, is enough to derive everything in the rest of this section — the potential, Gauss's law for gravity, Poisson's equation, the shell theorem, planetary orbits, escape velocity, even (with general-relativistic corrections) the perihelion of Mercury. One inverse-square law unfolds into the entire structure of celestial mechanics.

Compare side-by-side with the electrostatic version you just learned. The equations are almost identical, but the sign and the constant tell a story:

The key difference: mass has no negatives. Every source pulls. The wells of the gravitational potential never invert — and that is why everything eventually clumps under gravity, while charge distributions can shield each other out. Cosmology, planet formation, structure formation: all consequences of that single sign.

Newton's Inverse-Square Law

In 1687 Newton published the law that started it all: any two point masses $M$ and $m$ separated by $r$ attract each other along the line between them with magnitude

The minus sign and the unit vector $\hat{\mathbf{r}}$ (pointing from $M$ toward $m$) make the force attractive: it always pulls $m$ back toward $M$. Three things to chew on before we generalize.

The inverse-square shape

Why the exponent 2? Geometrically, a point source "spreads its influence" over the surface of a sphere of radius $r$, and that surface has area $4\pi r^2$. Whatever the total outward-pointing influence is (call it "flux"), it must thin out as $1/(4\pi r^2)$ to keep the integral conserved. The inverse-square law is, at its heart, a statement about the dimensionality of space.

Analogy — paint on a balloon. Take a fixed amount of paint and spray it uniformly outward from a tiny nozzle. At distance $r$ from the nozzle, the paint covers a sphere of area $4\pi r^2$. So the surface density of paint (mass per unit area) is total mass divided by $4\pi r^2$. That is the inverse-square law. Gravity isn't paint, but the geometric reasoning is the same.

Superposition: gravity is linear

If many masses act on a single test particle, the total force is the vector sum of the individual forces:

Linearity is profound: it means we can decompose any extended mass distribution into infinitesimal point masses, compute each pull separately, and add them up. Every continuous distribution becomes an integral. Every numerical galaxy simulation in the last 60 years is built on this.

Mass attracts mass — there is no "anti-gravity"

Unlike electric charge, mass has only one sign. There is no "negative mass" in classical Newtonian gravity. Two consequences:

• You cannot shield gravity. Pile up a wall of lead between you and the Earth and the Earth still pulls just as hard on you.
• Gravitational potential energy is always negative (taking $\Phi(\infty) = 0$). A bound system has lower energy than its separated parts — this is why stars and planets form spontaneously out of diffuse gas.

Why a Potential? From Vectors to a Scalar

Newton's law gives you the force. But carrying around a three-component vector for every particle gets tedious — and worse, it hides a beautiful structure. The fix: introduce a scalar potential $\Phi$ such that the gravitational acceleration is its (negative) gradient,

The acceleration $\mathbf{g}$ (also called the gravitational field strength; units m/s²) is just the force per unit mass: $\mathbf{g} = \mathbf{F}/m$. For a point mass $M$ at the origin the potential turns out to be

which you can verify by differentiating: $-\nabla(-GM/r) = -GM\,\hat{\mathbf{r}}/r^2 = \mathbf{g}$ — Newton's law recovered. So $\Phi(r) = -GM/r$ carries exactly the same information as the force law but as a single number per point.

Three reasons the potential is the right object

1. Scalars superpose more cleanly than vectors. For $N$ masses, summing $N$ potentials is one addition per point; summing $N$ force vectors is three additions. More importantly, the potential of a continuous distribution is a single triple integral, while the field is three.
2. Energy is potential, not force. The work done against gravity to move a unit mass from $A$ to $B$ is exactly $\Phi(B) - \Phi(A)$ — a path-independent difference. Try writing that statement using only the force: you need a line integral and a curl-free assumption.
3. Boundary-value problems are tractable in the potential. The equation $\nabla^2 \Phi = 4\pi G\rho$ is a single linear PDE with well-developed analytic and numerical machinery — the same machinery we used for the electric potential.

Mental model — the bowl. Picture a bowl whose depth at a point is the local potential. A marble placed anywhere rolls in the direction of steepest descent — that is $-\nabla\Phi$. The bottom of the bowl is the deepest well; the rim is far from any mass where $\Phi \to 0$. Every mass digs a well in the bowl, and wells from many masses superpose to form the cosmic landscape that planets, stars, and galaxies live in.

Interactive: The Gravitational Well

Drag the sliders. The left panel shows the inverse-square force per unit mass $F(r) = -GM/r^2$; the right panel shows the corresponding potential well $\Phi(r) = -GM/r$. Notice how the force diverges at $r \to 0$ much faster than the potential does — that is the singularity all of classical gravitation hides under the rug, and is one of the reasons general relativity was eventually needed.

Things to try. Set $M$ to Earth's value (5.97) and slide $r$ to 6.371 — you should see $g \approx 9.81\,\text{m/s}^2$ and $\Phi \approx -6.25\times 10^7\,\text{J/kg}$. That potential value is the famous "escape energy per kilogram": to leave Earth's gravitational well from the surface you must supply at least $|\Phi|$ of kinetic energy per kg, giving an escape speed $v_{\text{esc}} = \sqrt{2|\Phi|} \approx 11.2\,\text{km/s}$.

Gauss's Law for Gravity

Just as the electric field has a Gauss law, so does gravity. The reasoning is identical because both forces are inverse-square. Pick any closed surface $S$ in space; let $M_{\text{enc}}$ be the total mass enclosed inside. Then

The minus sign says: the gravity-flux through any closed surface points inward, proportional to the enclosed mass. Compare with the electric version $\oint \mathbf{E}\cdot d\mathbf{A} = Q_{\text{enc}}/\varepsilon_0$ — same form, with $-4\pi G$ replacing $1/\varepsilon_0$ and the minus sign reflecting the attractive nature of gravity.

Why this matters: symmetry shortcuts

Gauss's law turns hard volume integrals into easy surface integrals when symmetry permits. Three workhorse examples:

The spherical case alone is responsible for: orbits, escape velocity, planetary masses inferred from satellites, and Newton's elegant proof that a uniform spherical shell exerts no net force on anything inside it (the shell theorem). That theorem is why Earth can be treated as a point mass when computing the moon's orbit, even though Earth is obviously not a point.

From Gauss to Poisson

The Poisson form of gravity is one substitution away. We take Gauss's law in integral form, convert it to a divergence equation, then substitute $\mathbf{g} = -\nabla\Phi$.

1. Apply the divergence theorem to Gauss's integral form. The flux integral on the left becomes a volume integral of the divergence: $\oint_S \mathbf{g}\cdot d\mathbf{A} = \int_V \nabla\cdot\mathbf{g}\,dV$. On the right, $M_{\text{enc}} = \int_V \rho\,dV$.
2. Equate the integrands. Since the volume is arbitrary, the integrands must agree pointwise: $\nabla\cdot\mathbf{g} = -4\pi G\,\rho.$
3. Use $\mathbf{g} = -\nabla\Phi$. Then $\nabla\cdot\mathbf{g} = -\nabla^2 \Phi$. Substituting into step 2 and flipping the sign:

That is Poisson's equation for gravity. Outside any mass, where $\rho = 0$, it reduces to Laplace's equation $\nabla^2 \Phi = 0$, and the potential is determined entirely by boundary data — just like in electrostatics outside the charged region.

One picture, two physics. The electric and gravitational Poisson equations are mathematically the same problem. Any algorithm that solves one — finite differences, multigrid, fast multipole, neural-network PDE solver — solves the other with one sign flip. That is why graduate courses in plasma physics, geophysics, astrophysics, and electrical engineering all spend a week or two on the same linear PDE.

Worked Example: Earth as a Uniform Sphere

Treat Earth as a ball of radius $R$ with uniform density $\rho_0 = M/(\tfrac{4}{3}\pi R^3)$. By spherical symmetry $\Phi$ depends only on $r$. The Laplacian in spherical coordinates for a radial-only function is

Two regions, two ODEs. Inside the ball the source is constant $\rho_0$; outside the source is zero. We solve each region, then glue them at $r=R$ so that $\Phi$ and $d\Phi/dr$ are continuous (no infinite-density surface layer). Add the natural boundary condition $\Phi(\infty) = 0$ and the answer is uniquely determined.

The piecewise result is the gravitational mirror of the charged-sphere formula from the previous section:

• Outside the ball, gravity is identical to a point mass at the centre — Newton's shell theorem in one line. Earth's 6,000-km radius is invisible to the moon: only the total mass matters.
• Inside, the field grows linearly: $g(r) = -d\Phi/dr = GM r/R^3$ for $r \le R$. So if you fell down a frictionless shaft through the centre of Earth you would oscillate in simple harmonic motion with period $T = 2\pi\sqrt{R/g_{\text{surf}}}$ ≈ 84 minutes.
• The well is deepest at the centre: $\Phi(0) = -3GM/(2R)$, which is exactly $1.5\times$ the surface value $\Phi(R) = -GM/R$.

Interactive: Mass Distributions to Potential Fields

Click anywhere to plant a mass. The 64×64 grid behind the canvas solves $\nabla^2 \Phi = 4\pi G\rho$ by Jacobi relaxation — every cell updates to the average of its neighbours minus a source bump. The colormap goes from white (potential ≈ 0) at the edges to indigo at the bottom of the wells, and the arrows show the gravitational acceleration $\mathbf{g} = -\nabla\Phi$, which always points toward mass.

What to play with. Two masses far apart: each digs its own well and the saddle between them is a Lagrange-like point where $\mathbf{g}$ vanishes. A tight cluster: wells merge into a single deeper well — exactly how galaxies look in gravitational potentials. Place a ring of masses: at the centre, the radial arrows from each mass cancel out by symmetry (the shell theorem in flat-land form). Crank up the iterations to converge tightly; drop them low to watch Jacobi propagate information one cell per sweep.

Computing It: Direct N-Body in Python

Before reaching for a grid-based PDE solver, you should know the most direct way to compute the gravitational potential: sum it. For a set of point masses, Newton's law gives $\Phi(\mathbf{r}) = -G\sum_i m_i/|\mathbf{r}-\mathbf{r}_i|$ as a literal arithmetic sum. The implementation below evaluates this at a single field point. The trick: every step is direct addition — no PDE, no iteration. It is the rawest possible expression of the law.

Compare: G ≈ 6.67e-11. Coulomb's k ≈ 8.99e+9. Ratio ≈ 1e20 in favour of electromagnetism.

phi = 0.0  →  after particle 1: phi = -3.99e7  →  after particle 2: phi = -7.97e7

point = (0,0,0), pos = (-5e6, 0, 0)  →  dx=+5e6, dy=0, dz=0

dx=5e6, dy=0, dz=0  →  r = sqrt(25e12) = 5.0e6 m

term = -6.674e-11 * 5.972e24 / 5.0e6 = -3.987e7 J/kg per particle.

dx=5e6,dy=0,dz=0  →  r² = 25e12,  r³ = 25e12 · 5e6 = 1.25e20

At midpoint, particle 1 at -5e6: dx=+5e6, scale = -6.67e-11·5.97e24/1.25e20 = -3.19e-6.
Contribution to gx = scale·dx = -3.19e-6 · 5e6 = -15.95 m/s² ... but particle 2 contributes +15.95, exact cancellation.

particles[0] = (5.972e24 kg, (-5e6, 0, 0))
particles[1] = (5.972e24 kg, (+5e6, 0, 0))

phi_mid = -7.9716e7 J/kg.

g_mid ≈ (0, 0, 0) m/s².

1import math
2
3G = 6.674e-11   # N m^2 / kg^2
4
5def potential_at(point, particles):
6    """Gravitational potential Phi at a single point in space.
7
8    Phi(r) = -G * sum_i  m_i / |r - r_i|
9
10    Args:
11        point:     (x, y, z) location where we want Phi
12        particles: list of (mass_kg, (x, y, z)) tuples
13    Returns:
14        Phi in J/kg (negative number)
15    """
16    phi = 0.0
17    for mass, pos in particles:
18        dx = point[0] - pos[0]
19        dy = point[1] - pos[1]
20        dz = point[2] - pos[2]
21        r = math.sqrt(dx*dx + dy*dy + dz*dz)
22        if r == 0.0:
23            continue                       # skip self-interaction
24        phi += -G * mass / r
25    return phi
26
27def field_at(point, particles):
28    """Gravitational acceleration g = -grad Phi at a point.
29
30    g(r) = -G * sum_i  m_i (r - r_i) / |r - r_i|^3
31    """
32    gx, gy, gz = 0.0, 0.0, 0.0
33    for mass, pos in particles:
34        dx = point[0] - pos[0]
35        dy = point[1] - pos[1]
36        dz = point[2] - pos[2]
37        r2 = dx*dx + dy*dy + dz*dz
38        if r2 == 0.0:
39            continue
40        r3 = r2 * math.sqrt(r2)
41        scale = -G * mass / r3             # negative -> points TOWARD mass
42        gx += scale * dx
43        gy += scale * dy
44        gz += scale * dz
45    return (gx, gy, gz)
46
47# --- Demo: two Earth-mass particles 1e7 m apart ---
48ME = 5.972e24
49particles = [
50    (ME, (-5.0e6, 0.0, 0.0)),
51    (ME, (+5.0e6, 0.0, 0.0)),
52]
53
54midpoint = (0.0, 0.0, 0.0)
55phi_mid  = potential_at(midpoint, particles)
56g_mid    = field_at(midpoint, particles)
57
58print(f"Phi at midpoint  = {phi_mid:.4e} J/kg")
59print(f"g   at midpoint  = ({g_mid[0]:.3e}, {g_mid[1]:.3e}, {g_mid[2]:.3e}) m/s^2")

Why direct summation matters. The Poisson PDE and the direct sum agree because the kernel $-1/r$ is literally the Green's function of the 3D Laplacian (up to a constant) — that is the content of Green's functions you saw in Section 3. So for free-space problems with no boundary conductors, direct sum and grid-based Poisson solve are two equivalent computations of the same $\Phi$. Use direct sum when you have many point particles and few field points; use grid Poisson when you have a continuous $\rho(\mathbf{r})$ and want $\Phi$ everywhere.

Scaling Up: PyTorch on the GPU

Direct summation is $O(N K)$ — one term per particle per field point. With a million particles and a million field points the Python version above would take days. PyTorch lets us trade nested for-loops for tensor broadcasting and run the whole computation as one parallel kernel on a GPU. Below we build a toy "galaxy" of 200 stars and evaluate $\Phi(x, 0, 0)$ along a line slicing through the galactic plane.

K=3 points, N=2 particles  →  disp shape = (3, 2, 3).  disp[0,0] = grid_points[0] - positions[0].

Without eps: at r=0, phi = -inf. With eps=1e3: phi capped at -G·m/1e3, finite.

disp[0,0] = (5e15, 0, 0)  →  r²=2.5e31, r = sqrt(2.5e31 + 1e6) ≈ 5.00e15 m (eps invisible here).

After reduction: phi shape = (K,) = (200,).  Multiplied by -G the values become negative.

On a Colab GPU: device='cuda'. On a laptop with no GPU: device='cpu'. Code path identical.

masses = [2.0e30, 2.0e30, ..., 2.0e30]  shape (200,).  Total = 200 M_sun.

Uniform u in [0,1]:  fraction below u  vs fraction inside r = R·sqrt(u):
  u=0.25  →  r=R/2  →  area = πR²/4 = ¼ of total ✓

positions[0] might be (3.2e15, -1.8e15, 6.4e14)  → 3.7e15 m from the centre.

1import torch
2
3G = 6.674e-11
4
5def potential_grid_3d(masses, positions, grid_points):
6    """Gravitational potential on a batch of evaluation points.
7
8    Phi(r_k) = -G * sum_i  m_i / ||r_k - r_i||
9
10    Args:
11        masses:      (N,)        masses of N particles
12        positions:   (N, 3)      particle positions
13        grid_points: (K, 3)      points where Phi is wanted
14    Returns:
15        phi: (K,) potential values (J/kg, all <= 0)
16    """
17    # Pairwise displacements:  (K, 1, 3) - (1, N, 3) -> (K, N, 3)
18    disp = grid_points.unsqueeze(1) - positions.unsqueeze(0)
19    # Distances with a small softening eps to avoid r=0:
20    eps  = 1.0e3
21    r    = torch.sqrt((disp * disp).sum(dim=-1) + eps * eps)   # (K, N)
22    # Per-pair contribution -G m / r, then sum over particles:
23    phi  = -G * (masses.unsqueeze(0) / r).sum(dim=-1)          # (K,)
24    return phi
25
26# --- Build a small galaxy: 200 stars in a flat disk ---
27torch.manual_seed(0)
28device = "cuda" if torch.cuda.is_available() else "cpu"
29
30N = 200
31masses    = torch.full((N,), 2.0e30, device=device)             # 1 solar mass each
32radii     = 1.0e16 * torch.rand(N, device=device).sqrt()         # 0..1e16 m, uniform in area
33angles    = 2 * torch.pi * torch.rand(N, device=device)
34positions = torch.stack([
35    radii * torch.cos(angles),
36    radii * torch.sin(angles),
37    0.10 * radii * torch.randn(N, device=device),                # thin disk
38], dim=-1)                                                       # (200, 3)
39
40# --- Evaluate Phi along the x-axis through the disk ---
41K = 200
42xs = torch.linspace(-2e16, 2e16, K, device=device)
43grid_points = torch.stack([xs, torch.zeros_like(xs), torch.zeros_like(xs)], dim=-1)
44
45phi = potential_grid_3d(masses, positions, grid_points)
46
47print(f"min Phi (well bottom) = {phi.min().item():.3e} J/kg")
48print(f"Phi at x=0            = {phi[K//2].item():.3e} J/kg")
49print(f"Phi at x=+2e16        = {phi[-1].item():.3e} J/kg")

With a softening length $\epsilon = 10^3\,\text{m}$ the potential is bounded everywhere — no infinities even when a probe lands on a star. The minimum value of $\Phi$ along the slice gives the depth of the cluster's gravitational well, and the rate at which $\Phi \to 0$ at large $|x|$ tells you the total mass (because at large distances the cluster looks like a single point of mass $M_{\text{tot}}$, and $\Phi \approx -GM_{\text{tot}}/|x|$).

From toy to real. Real cosmological simulations (Millennium, IllustrisTNG, FIRE) follow billions of particles through cosmic time. Direct summation at $N=10^9$ would be impossible — those codes use tree codes (Barnes-Hut, $O(N\log N)$) or particle-mesh hybrids that solve Poisson on a grid for the long-range force and sum directly only for nearby pairs. Every one of those algorithms is a smart approximation of the inner loop you just wrote.

From Tides to Cosmology to ML

The gravitational potential is doing more work in modern science and engineering than you might guess. A few highlights to anchor the idea.

Tides — the second derivative of $\Phi$

Standing on Earth, you don't feel the moon's pull as a uniform shift — you feel its gradient. The tidal force on a body of size $L$ at distance $d$ from a mass $M$ is $F_{\text{tidal}} \sim G M L/d^3$ — the second derivative of $\Phi$. The reason ocean tides exist at all is that $\Phi$ isn't flat across Earth: the near side is pulled more strongly than the far side. The same $1/d^3$ falloff is why a planet too close to a black hole gets shredded (spaghettification).

Orbits as motion in a 1D effective potential

For a body in orbit around a point mass, conservation of angular momentum $\ell = m r^2 \dot\theta$ lets you eliminate the angular variable and write the energy as if the body lived in a 1D potential

The two terms compete: the gravitational well pulls inward, the centrifugal barrier pushes outward. The minimum of $\Phi_{\text{eff}}$ sits at the circular-orbit radius; energies above the minimum give elliptic orbits, exactly at zero give parabolic escape, above zero give hyperbolic flybys. Every Kepler problem in the universe lives on this one curve.

Cosmology and dark matter

Galaxy rotation curves measure $v(r)$ for stars orbiting a galactic centre. From Newton's law $v^2/r = GM_{\text{enc}}(r)/r^2$, so $M_{\text{enc}}(r) = v^2 r/G$. Observed curves stay flat far out, implying $M_{\text{enc}}(r) \propto r$ well beyond the luminous matter. The Poisson equation for the inferred mass distribution doesn't add up unless there is invisible matter — that is the observational origin of dark matter. Everything here uses nothing more advanced than $\nabla^2 \Phi = 4\pi G\rho$.

Machine learning connections

• Score-based generative models & Poisson flow networks. The 2022/2023 PFGM family (Xu, Jiang, Karras, et al.) explicitly trains a neural network to learn the gradient of a 3D Poisson potential whose source is the training data distribution. Sampling = following $-\nabla \Phi$ back toward the wells (the data). Same equation as gravity, with images as "masses."
• Graph neural networks via the Laplacian. Replace the continuous Laplacian by the graph Laplacian $L = D - A$. Smoothing signals over a graph by $(I - \alpha L) x$ is exactly one Jacobi sweep of $\nabla^2 \Phi = 0$ on the discrete graph. Many GNN message-passing schemes can be read as "solve Laplace on the data graph."
• Physics-informed neural networks (PINNs). Same idea as in the electrostatic section: train a net $\Phi_\theta(x,y,z)$ by minimizing $\|\nabla^2 \Phi_\theta - 4\pi G\rho\|^2$ over sampled points plus boundary loss. The Laplacian inside the loss is computed by autograd. PINNs for Newton-gravity reproduce planetary potentials with parameters trained from a handful of orbit samples.
• N-body emulators. Modern cosmology uses convolutional neural networks trained on small high-resolution simulations to predict the density-to-potential map for huge volumes, replacing FFT-based Poisson solvers. The training target is exactly the inverse Laplacian.

Summary

The takeaway in one sentence. Newton's inverse-square law, repackaged as $\nabla^2 \Phi = 4\pi G\rho$, lets the whole celestial sphere — planet, star, galaxy, dark matter halo — be solved by the same averaging-plus-source idea that solved electrostatics, with one difference: mass has no negatives, so every well goes only one way, and the universe (eventually) wants to clump.