Soddy Sphere Packing

DEFINITION

(b₁+b₂+b₃+b₄+b₅)² = 3(b₁²+b₂²+b₃²+b₄²+b₅²)
b₅′ = b₁+b₂+b₃+b₄ − b₅, and the same on (b·x, b·y, b·z)
section of a sphere at distance d: radius √(r² − d²)

NOTES

Five spheres can all touch each other at once, and when they do their curvatures are locked together: the square of the sum of all five is three times the sum of their squares. Frederick Soddy set that down in Nature in 1936, as the sphere stanza of a poem that first restates the four-circle theorem of Descartes; Thorold Gosset answered in the same journal the following January with the n-dimensional form, where the three becomes n. Read as a quadratic in one of the five curvatures the relation has two roots, so four mutually tangent spheres admit two fifths, and the two roots sum to the four. That is the whole construction: the other sphere is b₁+b₂+b₃+b₄ minus b₅, with no square root taken and no sign to choose. In the plane the same reflection carries a coefficient of two; in space it carries a coefficient of one, which is why a packing that begins with whole-number curvatures keeps them forever. The centres come free with the curvatures. Write a sphere as the four numbers b, b·x, b·y and b·z — curvature first, then curvature times centre — and the identical reflection applies to all four at once; dividing the last three by the first gives the centre back. The packing here starts from the curvatures −1, 2, 2, 3, 3: a unit sphere holding two of radius one half and two of radius one third. Every curvature the recurrence makes from it is a whole number, and the six spheres of curvature 3 turn out to form a closed ring at sixty degree spacing — the hexlet, from the paper Soddy wrote the year after the poem. No curvature in this packing is ever one more than a multiple of three, and the reason is short: the five root curvatures sum to nine, the reflection moves one of them by that sum less three times itself, and the sum stays a multiple of three at every step — so every curvature keeps the remainder of the one it replaced, and the remainders of the five it started with were only ever nought and two. What the plate draws is not the packing but a plane cut through it, drifting. Every sphere the plane enters contributes one circle, of radius the square root of r squared minus d squared, and those circles are almost never tangent to one another: two spheres touch at a single point, and a plane that misses that point meets both of them in circles that miss. So the section is a visibly different object from the gasket one dimension down — a scatter of separate rings rather than a chain of tangencies — and it breathes as the plane travels. Three dimensions is where the construction stops. Graham, Lagarias, Mallows, Wilks and Yan proved the group is discrete for two and three dimensions and not discrete above them, and the spheres it generates in four dimensions and higher overlap: space is the last dimension in which this is a packing at all.

PROVENANCE

Origin

F. Soddy, "The Kiss Precise", Nature 137 (20 June 1936), 1021 — the five-sphere relation is its closing stanza. The n-dimensional generalisation is a note of the same title in Nature 139 (9 January 1937), 62, doi:10.1038/139062a0. It is attributed to Thorold Gosset throughout the literature, and the atlas says so, but the record it is checked against does not: the Crossref entry for that note carries a title, a volume, a page and a date and NO author field at all. The attribution here rests on the secondary literature, not on the identifier.

The one it is not

The relation for four mutually tangent circles is the Descartes circle theorem, which the same poem restates a stanza earlier, and it is what PL. 36 draws. This plate is the sphere case: five spheres rather than four circles, three where the circle case has two, and a reflection of coefficient one where the circle case has coefficient two. The integer packings and the hexlet are named in a separate paper, Nature 139 (1937), 77–79, doi:10.1038/139077a0, which Graham, Lagarias, Mallows, Wilks and Yan call perhaps the first observation of integrality in these packings.

Curvature-center coordinates

R. L. Graham, J. C. Lagarias, C. L. Mallows, A. R. Wilks and C. H. Yan, "Apollonian circle packings: number theory", J. Number Theory 100 (2003), 1–45, doi:10.1016/S0022-314X(03)00015-5, and the geometry series in Discrete and Computational Geometry: part I, DCG 34 (2005), 547–585, doi:10.1007/s00454-005-1196-9; parts II and III, DCG 35 (2006), 1–36 and 37–72, doi:10.1007/s00454-005-1195-x and doi:10.1007/s00454-005-1197-8. A dating caveat, because both forms circulate: Crossref records parts II and III as published online 17 October 2005 and printed in volume 35, issue 1, dated January 2006. The atlas cites the printed volume, so part I is 2005 and parts II and III are 2006. The four numbers per sphere are Definition 3.2 of part III and the reflection is its equation 4.1.

The last packing

Part III, Theorem 4.1: the Apollonian group is a discrete subgroup for n = 2 and n = 3 and is not discrete for any n ≥ 4. The consequence stated in the introduction of the same paper is that in four dimensions and above the orbits are no longer sphere packings, because the spheres overlap.

Standing

Public domain — a theorem of 1936 and its nineteenth-century antecedents. No patent ever applied.

Constants

depth is the reflection depth and the only cost governor, locked from regenerate on the apollonian precedent: each level multiplies the sphere count by about 3.1 and only ever resolves smaller circles into the same section. It starts at 6 because it was measured and not guessed — at depth 4 one plane position in twenty draws fewer than eight circles, and the shallow end of a slider that goes empty is the same defect as a slider that does nothing. tilt sets the angle of the cutting plane to the axis, travel how far along its own normal it drifts, rate how fast. All three are live across their whole declared ranges, measured the way regenerate jitters them, which is all of them at once rather than one at a time.

Source

doi:10.1038/1371021a0

HOUDINI · VEX

The same published mathematics as a Detail Wrangle body. Paste it into a Wrangle with Run Over set to Detail; every constant is the published value plus a tweak channel, so Create Spare Parameters gives a slider that starts where the paper does.

// FORMA — PL. 152 · SODDY SPHERE PACKING — Frederick Soddy, 1936 · Thorold Gosset, 1937
//   (b₁+b₂+b₃+b₄+b₅)² = 3(b₁²+b₂²+b₃²+b₄²+b₅²)
//   b₅′ = b₁+b₂+b₃+b₄ − b₅, and the same on (b·x, b·y, b·z)
//   section of a sphere at distance d: radius √(r² − d²)
// Paste into a Detail Wrangle (Run Over: Detail), no inputs needed.
// Written from the published mathematics, not adapted from any code.
// Constants arrive at their published values. Press the node's Create
// Spare Parameters button and every tweak becomes a slider — starting
// at 0, the published figure, and moving in the constant's own units.
// https://forma-gen.com/#plate=soddy

float p_depth  = 7 + chf('depth_tweak');        // reflection depth · live 6 .. 8
float p_tilt   = 34 + chf('tilt_tweak');        // plane tilt · live 0 .. 90
float p_travel = 0.62 + chf('travel_tweak');    // plane travel · live 0.15 .. 0.85

// The plate's own colour: FORMA's LATTICES accent as a cosine ramp,
// brightest near t = 0 and t = 1, near-black around t = 0.5.
vector forma_ramp(float t){
  return set(
    0.46 + 0.5 * cos(6.28318530718 * (t + 0)),
    0.4185 + 0.4549 * cos(6.28318530718 * (t + 0.05)),
    0.1389 + 0.151 * cos(6.28318530718 * (t + 0.1)));
}

// The Soddy-Gosset relation (b1+b2+b3+b4+b5)^2 = 3(b1^2+...+b5^2) read as a
// quadratic in one curvature: its two roots sum to the other four, so the
// OTHER sphere tangent to a given four is b1+b2+b3+b4 - b5 with no root taken
// and no sign to choose. In curvature-center coordinates (b, b*x, b*y, b*z)
// -- Graham, Lagarias, Mallows, Wilks and Yan, part III, Definition 3.2 -- the
// identical reflection applies to all four numbers at once (their equation
// 4.1), so the centres come out of the same line the curvatures do and are
// recovered by dividing by b. The plane case has coefficient 2 over a sum of
// three; space has coefficient 1 over a sum of four, which is why an integral
// packing in space stays integral.
//
// The root is the quintuple of curvatures -1, 2, 2, 3, 3: a unit sphere at the
// origin, two of radius 1/2 on the z axis, two of radius 1/3 on the waist.
// 9^2 = 3 * 27 exactly. VEX has no recursion, so the reflection tree is an
// explicit queue of quintuples, the same idiom apollonian.vex uses one
// dimension down.
//
// What is emitted is the plate's own picture: ONE PLANE SECTION of the
// packing, as closed polylines, in the section plane's own two axes with y
// negated the way every 2-D canvas port here negates it. The spheres
// themselves are right there in the arrays should a reader want them --
// centre set(bxs[i], bys[i], bzs[i]) / float(cbs[i]), radius 1/cbs[i] -- but
// the plate is the section, and the port is the plate.
//
// waived: rate — it paces the drift of the cutting plane, and a cook has no
// clock. The plane is parked at p_travel * sin(2*pi*PHASE), which is where the
// page puts it every time its own drift cycle comes back round to PHASE --
// once per drift period, whatever the rate is, so the cook shows a frame the
// page really renders rather than an independently chosen slice. forma_phase
// below is the computed PHASE_OF value for the id "soddy", by the same FNV
// hash the page uses.

float TAU          = 6.28318530718;
float forma_phase  = 0.596096069086343;
int   forma_seg    = 64;
float forma_minr   = 0.006;    // smallest sphere kept, as the plate keeps it
float forma_mincut = 0.002;    // and the plate's sub-pixel floor on a section
                               // circle, in the same units: 0.55 px against a
                               // 270 px bounding radius at drawer size

int depth = int(rint(p_depth));

// the packing, as parallel arrays: curvature and curvature-times-centre.
// The curvature is an int on purpose -- it comes out an exact whole number at
// every step, because the recurrence is a sum and difference of integers and
// never divides, and that is what makes it a safe bucket key below.
int   cbs[];
float bxs[], bys[], bzs[];
int   lvls[];
float rt3 = sqrt(3.0);
push(cbs, -1); push(bxs, 0.0); push(bys,  0.0); push(bzs,  0.0); push(lvls, 0);
push(cbs,  2); push(bxs, 0.0); push(bys,  0.0); push(bzs,  1.0); push(lvls, 0);
push(cbs,  2); push(bxs, 0.0); push(bys,  0.0); push(bzs, -1.0); push(lvls, 0);
push(cbs,  3); push(bxs, rt3); push(bys,  1.0); push(bzs,  0.0); push(lvls, 0);
push(cbs,  3); push(bxs, rt3); push(bys, -1.0); push(bzs,  0.0); push(lvls, 0);

// Spheres bucketed by curvature, as head-of-chain plus next-link, so the
// duplicate test costs a short walk instead of a scan of the whole list.
// Every curvature is a positive integer no larger than 1/forma_minr once the
// floor below has passed, so the table is sized from the floor itself and
// cannot be overrun.
int slots = int(ceil(1.0 / forma_minr)) + 2;
int bhead[], bnext[];
for (int i = 0; i < slots; i++){
    push(bhead, -1);
}
for (int i = 0; i < 5; i++){
    int k = cbs[i] > 0 ? cbs[i] : 0;   // the bounding sphere parks in slot 0
    push(bnext, bhead[k]);
    bhead[k] = i;
}

// The quintuple queue: five sphere indices and the level they sit at, walked
// with a read cursor rather than popped. That is not a style choice. A stack
// walks the reflection tree depth first, and depth first plus a duplicate test
// plus a depth bound builds a SMALLER packing: a sphere first reached down a
// deep branch is recorded at that deep level, and when the short branch reaches
// it later the duplicate test throws the short path away, so the sphere never
// gets the levels of expansion it was owed. Measured against the plate, which
// walks breadth first: 2,501 spheres instead of 6,340 at the published depth.
int q0[] = {0};
int q1[] = {1};
int q2[] = {2};
int q3[] = {3};
int q4[] = {4};
int ql[] = {0};
int head = 0;

while (head < len(q0)){
    int i0 = q0[head], i1 = q1[head], i2 = q2[head], i3 = q3[head], i4 = q4[head];
    int lvl = ql[head];
    head += 1;
    if (lvl >= depth) continue;
    int quint[] = array(i0, i1, i2, i3, i4);
    int sb = 0;
    float sx = 0.0, sy = 0.0, sz = 0.0;
    for (int j = 0; j < 5; j++){
        int m = quint[j];
        sb += cbs[m];  sx += bxs[m];  sy += bys[m];  sz += bzs[m];
    }
    for (int j = 0; j < 5; j++){
        int m = quint[j];
        int nb = sb - 2 * cbs[m];
        // curvature 0 is a plane and negative curvature is an enclosing
        // sphere: neither is a child inside a bounded packing, and the
        // bounding sphere is already in the list. These are the degeneracies
        // of the quadratic, refused here by name.
        if (nb <= 0) continue;
        if (1.0 / float(nb) < forma_minr) continue;
        float nx = sx - 2.0 * bxs[m];
        float ny = sy - 2.0 * bys[m];
        float nz = sz - 2.0 * bzs[m];
        // Duplicate test. Two distinct spheres of curvature b have centres at
        // least 2/b apart, so in these coordinates they differ by at least 2 --
        // far outside the tolerance below, and far outside anything the
        // arithmetic could drift: VEX floats are 32-bit, which at the largest
        // magnitude these coordinates reach is a few parts in 10^5, four
        // orders under the 0.5 the test allows.
        int dup = 0;
        for (int k = bhead[nb]; k >= 0 && dup == 0; k = bnext[k]){
            if (abs(bxs[k] - nx) + abs(bys[k] - ny) + abs(bzs[k] - nz) < 0.5) dup = 1;
        }
        if (dup == 1) continue;
        int ni = len(cbs);
        push(cbs, nb);  push(bxs, nx);  push(bys, ny);  push(bzs, nz);
        push(lvls, lvl + 1);
        push(bnext, bhead[nb]);
        bhead[nb] = ni;
        push(q0, j == 0 ? ni : i0);
        push(q1, j == 1 ? ni : i1);
        push(q2, j == 2 ? ni : i2);
        push(q3, j == 3 ? ni : i3);
        push(q4, j == 4 ? ni : i4);
        push(ql, lvl + 1);
    }
}

// The cutting plane: unit normal at p_tilt from the z axis on the bearing this
// plate takes from its own PHASE, offset along that normal. uax and vax are
// the other two directions of the same spherical frame, so they span the
// plane and the foot of the perpendicular from a centre lands at
// (dot(c, uax), dot(c, vax)) with no projection left to compute.
float tilt = p_tilt * M_PI / 180.0;
float psi  = forma_phase * TAU;
float ct = cos(tilt), st = sin(tilt), cs = cos(psi), sn = sin(psi);
vector nrm = set(st * cs, st * sn, ct);
vector uax = set(ct * cs, ct * sn, -st);
vector vax = set(-sn, cs, 0.0);
float off = p_travel * sin(TAU * forma_phase);

for (int i = 0; i < len(cbs); i++){
    float b = float(cbs[i]);
    float r = 1.0 / b;
    vector c = set(bxs[i], bys[i], bzs[i]) / b;
    // The bounding sphere needs no special case: its curvature is -1, so r is
    // negative, r*r is the same number every other sphere gives, and the
    // section radius comes out right.
    float d = dot(c, nrm) - off;
    float rr = r * r - d * d;
    if (rr <= 0.0) continue;                 // the plane misses this sphere
    float rad = sqrt(rr);
    if (rad < forma_mincut) continue;
    // sparse strokes on a dark ground, so the whole band map stays inside the
    // bright lobe of the ramp -- 0.82 to 1.06, never crossing the trough at
    // 0.5. The bounding sphere is the one negative curvature in the packing
    // and takes the accent itself, which is ramp(0).
    vector col = i == 0 ? forma_ramp(0.0)
                        : forma_ramp(0.82 + 0.24 * float(lvls[i]) / float(depth));
    float alpha = i == 0 ? 0.9 : 0.85;
    float px = dot(c, uax), py = dot(c, vax);
    int prim = addprim(0, "polyline");
    for (int k = 0; k <= forma_seg; k++){
        float th = TAU * float(k) / float(forma_seg);
        // canvas y runs down; negated so the section sits as the plate shows it
        int pt = addpoint(0, set(px + rad * cos(th), -(py + rad * sin(th)), 0.0));
        setpointattrib(0, "Cd", pt, col);
        setpointattrib(0, "Alpha", pt, alpha);
        addvertex(0, prim, pt);
    }
}

AFTER EFFECTS · DECLINED

Circles that mostly never touch: a plane section of a sphere packing is a scatter of separate rings, and no continuous path visits them.