Curvature for Sphere Linkages

Intrinsic and extrinsic curvature of the configuration space, and what each tells us about the linkage

Fourth in a series on spherical configuration spaces of planar linkages, joint work with Aaron Abrams, Dave Bachmann, and Edmund Harriss. The setup is in the first post.

The previous post found a metric on $\mathcal{M}^4_L \cong S^2$ whose geodesics are the linkage’s natural force-free trajectories — what the configuration does when you let go. But linkages are usually not force-free; you grab them, push them, hold one rod still while moving another. What does that feel like?

Here is the picture. The unconstrained position space of three free interior hinges is Euclidean $\mathbb{R}^6$. When you grab the linkage and try to move your hand, you’d like to move freely in $\mathbb{R}^6$, but the rigid rods constrain you to the 2-sphere $\mathcal{M}^4_L$ sitting inside. Your hand gets tugged around because that 2-sphere isn’t flat in $\mathbb{R}^6$ — it bends, and the bending is what you feel.

Two flavors of curvature

Curvature comes in two genuinely different flavors, measuring different things.

Intrinsic curvature is about the surface’s internal geometry — how things behave on the surface, with no reference to the ambient space. The Gaussian curvature $K(p)$ is a single scalar function reading off two related effects:

• Geodesic divergence. Two geodesics start at $p$ with slightly different velocities. Do they converge or diverge? On a flat plane they separate linearly; on a sphere they reconverge; on a hyperbolic surface they fly apart exponentially.
• Parallel transport holonomy. Carry a tangent vector around a closed loop, rotating it as little as possible relative to the surface. When it returns to its starting point, it has rotated by some angle — proportional to $\int K\,dA$ over the enclosed region.

A 2D bug confined to the surface, with no access to the ambient space, can detect $K$ through these effects.

Extrinsic curvature is about how the surface sits inside the ambient space — specifically, how it bends away from its tangent planes as you move across it. This is what you feel when you push the linkage around: stiffness in some directions, floppiness in others, and a constraint force pulling your hand back onto the surface.

The two flavors are independent — different surfaces can have the same intrinsic geometry and entirely different extrinsic geometries. The cleanest way to see this is a small example.

Example: A constrained ball on a plane and a cylinder.. Take a flat plane and a cylinder in $\mathbb{R}^3$. They have the same intrinsic geometry: you can roll the plane up into a cylinder without stretching, and a 2D bug on either surface sees the same flat geometry, with $K = 0$. Geodesics are straight lines on the plane and helices on the cylinder, but in both cases they neither converge nor diverge.

Now drag a ball along each. On the plane the ball coasts in a straight line — no constraint force needed. On the cylinder it wants to fly off tangentially, but the surface won’t let it, and you feel a force pulling it inward as it follows the curve. That force is the centripetal force, and providing it is what extrinsic curvature is about. Move at speed $v$ around the cross-section of a cylinder of radius $r$ and the surface pulls back with force $v^2/r$; move along the axis and there’s no force at all; diagonally, somewhere in between.

Same intrinsic geometry, totally different extrinsic geometry, and dragging a ball detects only the extrinsic.

The Linkage: The linkage is the same kind of object, just bigger: a 2-surface in $\mathbb{R}^6$ instead of $\mathbb{R}^3$. The constraint force on a configuration moving at velocity $\mathbf{v}$ is still the centripetal force you feel as your hand pushes against the rigid rods. What’s new is codimension — the normal space is now 4-dimensional rather than 1-dimensional, so the constraint force has a direction within it as well as a magnitude, and that direction depends on $\mathbf{v}$. The cylinder example couldn’t have hinted at this.

Computational tools

The first and second fundamental forms are tools for organizing the geometric content of the embedding’s derivatives. The first fundamental form packages the first derivatives of $\mathbf{p}(\phi, t) \in \mathbb{R}^6$; the second fundamental form packages the second. The names match the derivative count, and that’s the whole logic. This section sets both up; §2 and §3 read off curvature scalars from them.

The embedding itself we have from the parameterizing post:

function embed(phi, t, L) {
  // Ψ^4_L from the parameterizing post
  const alpha    = Math.acos((L*L - 8) / (2*L));
  const theta    = alpha * Math.sin(phi);
  const p3       = [L - Math.cos(theta), -Math.sin(theta)];
  const d        = Math.hypot(p3[0], p3[1]);
  const alpha_in = Math.acos((d*d - 3) / (2*d));
  const theta_in = alpha_in * Math.sin(t);
  const cti = Math.cos(theta_in), sti = Math.sin(theta_in);

  const rod3x = (p3[0]*cti - p3[1]*sti) / d;
  const rod3y = (p3[0]*sti + p3[1]*cti) / d;
  const p2 = [p3[0] - rod3x, p3[1] - rod3y];

  const p2_abs = Math.hypot(p2[0], p2[1]);
  const sigma  = Math.sign(Math.cos(t)) || 1;
  const nu_in  = sigma * Math.sqrt(Math.max(0, 4 - p2_abs*p2_abs));
  const a      = nu_in / (2 * p2_abs);
  const p1     = [p2[0]/2 - a*p2[1], p2[1]/2 + a*p2[0]];

  return [p1[0], p1[1], p2[0], p2[1], p3[0], p3[1]];
}

The first fundamental form

The first fundamental form is the metric — the inner product on tangent vectors induced by pulling back the ambient inner product through the embedding:

$$I(X, Y) \;=\; \langle d\mathbf{p}(X), d\mathbf{p}(Y)\rangle.$$

In coordinates,

$$I_{ab} \;=\; \langle \partial_a\mathbf{p}, \partial_b\mathbf{p}\rangle, \qquad a, b \in \{\phi, t\}.$$

We computed this in the last post - the metric tensor in our spherical coordinate system - for calculating geodesics. And while that gives exact answers, its also useful to point out that we can compute the first fundamental form numerically via central-difference approximations to the derivative.

This requires four calls to embed, two central differences, three dot products, which may be a comparable amount of work to evaluating our exact expression. But it proves useful when these numerical derivatives are also needed for other calculations, so we record it here.

function firstFundamentalForm(phi, t, L) {
  const eps = 1e-4;
  const pPh = embed(phi + eps, t, L);
  const pPl = embed(phi - eps, t, L);
  const pTh = embed(phi, t + eps, L);
  const pTl = embed(phi, t - eps, L);

  const dPp = pPh.map((x, i) => (x - pPl[i]) / (2*eps));
  const dTp = pTh.map((x, i) => (x - pTl[i]) / (2*eps));

  const dot = (a, b) => a.reduce((s, x, i) => s + x*b[i], 0);
  return {
    E: dot(dPp, dPp),
    F: dot(dPp, dTp),
    G: dot(dTp, dTp),
  };
}

The second fundamental form

The second fundamental form is the same construction one derivative up. The second derivative $\partial_a\partial_b\mathbf{p}$ is a vector in $\mathbb{R}^6$, in general neither tangent nor normal to the surface. Decompose it; the second fundamental form is the normal piece:

$$II_{ab} \;=\; (\partial_a\partial_b\mathbf{p})^\perp \;\in\; N_p\mathcal{M}^4_L.$$

The normal space $N_p\mathcal{M}^4_L \subset \mathbb{R}^6$ is 4-dimensional, so each $II_{ab}$ is a 4-dimensional vector. Symmetry in $a, b$ leaves three independent components $II_{\phi\phi}$, $II_{\phi t}$, $II_{tt}$.

Unlike the first form, we have no already-computed nalytic expression for $II$ — (we could of course compute second derivatives, but the expression would be much longer than what we found for the metric!) Numerical evaluation is straightforward: central differences for $\partial_a\partial_b\mathbf{p}$ on a 9-point stencil, followed by a normal projection.

To extract the normal part of an $\mathbb{R}^6$ vector $\mathbf{v}$, we need its tangential components $(c^\phi, c^t)$ such that $\mathbf{v}^\parallel = c^\phi\partial_\phi\mathbf{p} + c^t\partial_t\mathbf{p}$. These solve

$$\begin{pmatrix} E & F \\ F & G \end{pmatrix}\begin{pmatrix} c^\phi \\ c^t \end{pmatrix} \;=\; \begin{pmatrix} \langle\mathbf{v}, \partial_\phi\mathbf{p}\rangle \\ \langle\mathbf{v}, \partial_t\mathbf{p}\rangle \end{pmatrix}$$

— the matrix is $I$ in the coordinate basis, the right-hand side is $\mathbf{v}$ pulled back through the tangent vectors. Subtracting $\mathbf{v}^\parallel$ from $\mathbf{v}$ gives the normal projection.

function secondFundamentalForm(phi, t, L) {
  const eps = 1e-4, e2 = eps*eps;
  const p   = embed(phi,       t,       L);
  const pPh = embed(phi + eps, t,       L);
  const pPl = embed(phi - eps, t,       L);
  const pTh = embed(phi,       t + eps, L);
  const pTl = embed(phi,       t - eps, L);
  const pPP = embed(phi + eps, t + eps, L);
  const pPN = embed(phi + eps, t - eps, L);
  const pNP = embed(phi - eps, t + eps, L);
  const pNN = embed(phi - eps, t - eps, L);

  const dPp = pPh.map((x, i) => (x - pPl[i]) / (2*eps));
  const dTp = pTh.map((x, i) => (x - pTl[i]) / (2*eps));
  const dPPp = p.map((x, i) => (pPh[i] - 2*x + pPl[i]) / e2);
  const dTTp = p.map((x, i) => (pTh[i] - 2*x + pTl[i]) / e2);
  const dPTp = p.map((_, i) => (pPP[i] - pPN[i] - pNP[i] + pNN[i]) / (4*e2));

  const dot = (a, b) => a.reduce((s, x, i) => s + x*b[i], 0);
  const E = dot(dPp, dPp), F = dot(dPp, dTp), G = dot(dTp, dTp);
  const D = E*G - F*F;

  const projectNormal = v => {
    const bP = dot(v, dPp), bT = dot(v, dTp);
    const cP = (G*bP - F*bT) / D;
    const cT = (E*bT - F*bP) / D;
    return v.map((x, i) => x - cP*dPp[i] - cT*dTp[i]);
  };

  return {
    II_pp: projectNormal(dPPp),
    II_pt: projectNormal(dPTp),
    II_tt: projectNormal(dTTp),
  };
}

The two functions duplicate work — secondFundamentalForm already computes the first-derivative stencils and $E, F, G$ internally for the projection, then throws them away. In applications that want both forms at every point, it’s straightforward to merge the two into a single function that runs the stencil once and returns $E, F, G$ alongside $II_{\phi\phi}, II_{\phi t}, II_{tt}$.

Intrinsic curvature: $K$

Gaussian curvature $K$ is a single scalar function on the surface. One of its main uses is that it controls the behavior of geodesics — and since geodesics here are the linkage’s force-free trajectories, $K$ helps us understand the qualitative dynamics.

Roughly: where $K$ is positive, nearby geodesics converge; where $K$ is negative, they diverge. Precisely, the Jacobi equation

$$\ddot J \;+\; K(\gamma(s))\,J \;=\; 0$$

governs the perpendicular separation $J(s)$ between two infinitesimally-close geodesics, with $K$ evaluated along the reference. Take a configuration, give it two slightly-different initial nudges, integrate both, and the gap is amplified or damped point-by-point by the local sign of $K$. Since $\mathcal{M}^4_L \cong S^2$ topologically, $K$ is positive on average — but it can vary, and regions of small or negative $K$ are where the dynamics is most sensitive to initial conditions.

Computing $K$

Gaussian curvature is famously a function of the first fundamental form alone. The Brioschi formula writes this out explicitly:

$$K \;=\; \frac{1}{(EG - F^2)^2}\,\big(\det M_1 \;-\; \det M_2\big),$$

with $E = h_{\phi\phi}, F = h_{\phi t}, G = h_{tt}$, subscripts denoting partial derivatives, and

$$M_1 \;=\; \begin{pmatrix} -\tfrac{1}{2}E_{tt} + F_{\phi t} - \tfrac{1}{2}G_{\phi\phi} & \tfrac{1}{2}E_\phi & F_\phi - \tfrac{1}{2}E_t \\ F_t - \tfrac{1}{2}G_\phi & E & F \\ \tfrac{1}{2}G_t & F & G \end{pmatrix},$$ $$M_2 \;=\; \begin{pmatrix} 0 & \tfrac{1}{2}E_t & \tfrac{1}{2}G_\phi \\ \tfrac{1}{2}E_t & E & F \\ \tfrac{1}{2}G_\phi & F & G \end{pmatrix}.$$

This is the content of Gauss’s theorema egregium: $K$ depends only on $E, F, G$ and their partials, with no reference to the embedding. Two surfaces with the same intrinsic metric have the same $K$, even if one sits in $\mathbb{R}^6$ and the other lies flat in $\mathbb{R}^2$.

But while it’s true that curvature is a function of the first fundamental form, that’s not usually how it gets computed. In an undergraduate differential geometry course you’re more likely to learn that the Gaussian curvature of a surface in $\mathbb{R}^3$ is the product of the two principal curvatures,

$$K \;=\; \kappa_1 \kappa_2,$$

computed from a priori extrinsic data — the eigenvalues of the shape operator describing how the surface bends in the ambient space.

This works in general. The Gauss equation gives the Gaussian curvature in terms of both fundamental forms:

$$K \;=\; \frac{\langle II_{\phi\phi}, II_{tt}\rangle \;-\; |II_{\phi t}|^2}{EG - F^2}.$$

In codimension 1 this reduces to $\det(II)/\det(I) = \kappa_1\kappa_2$; in codimension 4 the $II_{ab}$ are vectors in the normal space and the scalar products in the numerator become ambient inner products. Either way, theorema egregium guarantees agreement with Brioschi.

function gaussianCurvature(phi, t, L) {
  const { E, F, G } = firstFundamentalForm(phi, t, L);
  const { II_pp, II_pt, II_tt } = secondFundamentalForm(phi, t, L);

  const dot = (a, b) => a.reduce((s, x, i) => s + x*b[i], 0);
  return (dot(II_pp, II_tt) - dot(II_pt, II_pt)) / (E*G - F*F);
}

(This is one of those applications that wants both fundamental forms at every point — the closing remark of §1 applies, and a merged version would skip the duplicated stencil work.)

The first demo colors the sphere $\mathcal{M}^4_L$ by $K$, with the linkage shown alongside; slide $L$ to see how the curvature distribution depends on the rod length.

The second demo makes the Jacobi equation visible: click on the sphere and drag/release to launch ten nearby geodesics. Watch them converge or spread according to the local sign of $K$, with ten superimposed linkages animating in step. Try turning $L$ down to a lower value so there’s some negative curvature to play with!

Extrinsic curvature: $II$

Extrinsic curvature, as §0 set it up, is what your hand feels pushing the linkage around. The second fundamental form $II$ from §1 packages that information — a vector in the 4-dimensional normal space at each configuration, depending quadratically on velocity. To plot it, we extract scalar fields. Two natural ones: total bending $|II|^2$, a single number at each point, and directional stiffness $|II(\mathbf{v}, \mathbf{v})|^2$, a function of the push direction $\mathbf{v}$ at each point.

Both are most easily stated in an orthonormal tangent frame. Gram-Schmidt on the coordinate vectors $\partial_\phi\mathbf{p}, \partial_t\mathbf{p}$ gives

$$\mathbf{e}_1 \;=\; \frac{1}{\sqrt{E}}\,\partial_\phi\mathbf{p}, \qquad \mathbf{e}_2 \;=\; \frac{1}{\sqrt{ED}}\big(E\,\partial_t\mathbf{p} \;-\; F\,\partial_\phi\mathbf{p}\big),$$

and the second fundamental form has components $II_{ab} = II(\mathbf{e}_a, \mathbf{e}_b) \in N_p$ in this frame, derived from $II_{\phi\phi}, II_{\phi t}, II_{tt}$ by change of basis:

$$\begin{aligned} II_{11} \;&=\; \tfrac{1}{E}\,II_{\phi\phi}, \\ II_{12} \;&=\; -\tfrac{F}{E\sqrt{D}}\,II_{\phi\phi} \;+\; \tfrac{1}{\sqrt{D}}\,II_{\phi t}, \\ II_{22} \;&=\; \tfrac{F^2}{ED}\,II_{\phi\phi} \;-\; \tfrac{2F}{D}\,II_{\phi t} \;+\; \tfrac{E}{D}\,II_{tt}, \end{aligned}$$

each a vector in the 4-dimensional normal space.

Total bending $|II|^2$

The squared norm of $II$ as a bilinear form. In any orthonormal tangent frame:

$$|II|^2 \;=\; |II_{11}|^2 \;+\; 2\,|II_{12}|^2 \;+\; |II_{22}|^2,$$

a single scalar at each point summarizing the magnitude of $II$ across all directions. Higher means the linkage is stiffer on average — every push gets pushed back harder.

function totalBending(phi, t, L) {
  const { E, F, G } = firstFundamentalForm(phi, t, L);
  const { II_pp, II_pt, II_tt } = secondFundamentalForm(phi, t, L);
  const D = E*G - F*F, sqrtD = Math.sqrt(D);

  // II in the orthonormal frame (vectors in R^6, lying in N_p)
  const II_11 = II_pp.map(x => x / E);
  const II_12 = II_pp.map((x, i) => -F/(E*sqrtD) * x + II_pt[i]/sqrtD);
  const II_22 = II_pp.map((x, i) =>
    (F*F)/(E*D) * x  -  (2*F/D) * II_pt[i]  +  (E/D) * II_tt[i]
  );

  const dot = (a, b) => a.reduce((s, x, i) => s + x*b[i], 0);
  return dot(II_11, II_11) + 2*dot(II_12, II_12) + dot(II_22, II_22);
}

The figure below colors $\mathcal{M}^4_L$ by $|II|^2$ for several values of $L$.

Directional stiffness $|II(\mathbf{v}, \mathbf{v})|^2$

Total bending averages over push directions; directional stiffness keeps the direction. For a unit tangent

$$\mathbf{v}(\psi) \;=\; \cos\psi\,\mathbf{e}_1 \;+\; \sin\psi\,\mathbf{e}_2,$$

bilinearity of $II$ gives

$$II(\mathbf{v}(\psi), \mathbf{v}(\psi)) \;=\; \cos^2\psi\,II_{11} \;+\; 2\sin\psi\cos\psi\,II_{12} \;+\; \sin^2\psi\,II_{22},$$

an $\mathbb{R}^6$-valued function of $\psi$. Its squared norm $|II(\mathbf{v}(\psi), \mathbf{v}(\psi))|^2$ is the magnitude of the constraint force at unit speed in direction $\mathbf{v}(\psi)$ — a real-valued, degree-4 trigonometric polynomial in $\psi$. Decomposed into Fourier modes $\{1, \cos 2\psi, \sin 2\psi, \cos 4\psi, \sin 4\psi\}$:

• The constant mode is essentially $|II|^2$ (up to a factor): direction-averaged stiffness.
• The $2\psi$ modes give the principal stiffness direction and the anisotropy: max and min stiffness, 90° apart.
• The $4\psi$ modes are higher-order shape that doesn’t fit into a max-and-min summary.

A clean visualization of the directional dependence is to draw, at sample points on the sphere, the ellipse capturing the constant + $2\psi$ part: principal direction is the major axis, anisotropy is the axis ratio. The $4\psi$ part gets dropped, which is fine — the ellipse summary is the standard way to display a symmetric tensor field, and the $4\psi$ remainder is what an ellipse can’t say.

function directionalStiffness(phi, t, psi, L) {
  const { E, F, G } = firstFundamentalForm(phi, t, L);
  const { II_pp, II_pt, II_tt } = secondFundamentalForm(phi, t, L);
  const D = E*G - F*F, sqrtD = Math.sqrt(D);

  const II_11 = II_pp.map(x => x / E);
  const II_12 = II_pp.map((x, i) => -F/(E*sqrtD) * x + II_pt[i]/sqrtD);
  const II_22 = II_pp.map((x, i) =>
    (F*F)/(E*D) * x  -  (2*F/D) * II_pt[i]  +  (E/D) * II_tt[i]
  );

  const c = Math.cos(psi), s = Math.sin(psi);
  const IIvv = II_11.map((_, i) =>
    c*c*II_11[i] + 2*s*c*II_12[i] + s*s*II_22[i]
  );
  return IIvv.reduce((sum, x) => sum + x*x, 0);
}