Part I - Proof of Soddy-Gosset theorem (generalization of Descartes theorem).

For any integer $d \ge 2$, consider the problem of placing $n = d + 2$ hyper-spheres touching each other in $\mathbb{R}^d$. Let $\vec{x}_i \in \mathbb{R}^d$ and $R_i \in \mathbb{R}$ be the center and radius for the $i^{th}$ sphere. The condition for these spheres touching each other can be expressed as:

$$|\vec{x}_i - \vec{x}_j| = | R_i + R_j | \quad\text{ for }\quad 1 \le i < j \le n$$

or equivalently

$$|\vec{x}_i - \vec{x}_j|^2 = (R_i + R_j)^2 - 4R_iR_j\delta_{ij}\quad\text{ for }\quad 1 \le i, j \le n\tag{*1}$$

where $\;\delta_{ij} = \begin{cases}1,&i = j\\0,& i \ne j\end{cases}\;$ is the Kronecker delta.

Since the $n = d+2$ points $\vec{x}_i$ live in $\mathbb{R}^d$, the $d+1$ vectors $\vec{x}_2 - \vec{x}_1, \vec{x}_3 - \vec{x}_1, \ldots, \vec{x}_n - \vec{x}_1$ are linear depedent, this means we can find $n-1$ numbers $\beta_2, \beta_3, \ldots, \beta_n$ not all zero such that $$\sum_{k=2}^n \beta_k (\vec{x}_k - \vec{x}_1 ) = \vec{0}$$ Let $\beta_1 = -(\beta_2 + \ldots + \beta_n)$, we can rewrite this relation in a more symmetric form:

$$ \sum_{k=1}^n \beta_k = 0 \quad\text{ and }\quad \sum_{k=1}^n \beta_k \vec{x}_k = \vec{0} \quad\text{ subject to some }\;\; \beta_k \ne 0 $$

If we fix $j$ in $(*1)$, multiple the $i^{th}$ term by $\beta_i$ and then sum over $i$, we get

$$ \sum_{i=1}^n\beta_i |\vec{x}_i|^2 = \sum_{i=1}^n\beta_i R_i^2 + 2 \left( \sum_{i=1}^n \beta_i R_i \right) R_j - 4R_j^2 \beta_j $$ This leads to $$4R_j^2 \beta_j = 2 A R_j + B \quad\text{ where }\quad\ \begin{cases} A &= \sum\limits_{i=1}^n \beta_i R_i\\ B &= \sum\limits_{i=1}^n\beta_i ( R_i^2 - |\vec{x}_i|^2 ) \end{cases} \tag{*2} $$ Divide $(*2)$ by $R_j$ and sum over $j$, we get

$$4A = 2nA + B\sum_{j=1}^n\frac{1}{R_j}\quad\iff\quad A = -\frac{B}{2d}\sum_{j=1}^n\frac{1}{R_j}\tag{*3}$$

A consequence of this is $B$ cannot vanish. Otherwise $B = 0 \implies A = 0$ and $(*2)$ implies all $\beta_j = 0$ which is clearly isn't the case.

Divide $(*2)$ by $R_j^2$ and sum over $j$, we get

$$0 = 4\sum_{j=1}^n \beta_j = 2A\sum_{j=1}^n\frac{1}{R_j} + B\sum_{j=1}^n\frac{1}{R_j^2}$$

Combine with $(*3)$, the RHS becomes

$$B \left( \sum_{j=1}^n \frac{1}{R_j^2} - \frac{1}{d}\left( \sum_{j=1}^n\frac{1}{R_j}\right)^2\right) = 0 \quad\iff\quad \left( \sum_{j=1}^n\frac{1}{R_j}\right)^2 = d\sum_{j=1}^n \frac{1}{R_j^2}\tag{*4} $$ The RHS of $(*4)$ is sometimes called Soddy-Gosset theorem. When $d = 2$, it reduces to the Descartes four circle theorem, the theorem we wish to prove: $$\left( \frac{1}{R_1} + \frac{1}{R_2} + \frac{1}{R_3} + \frac{1}{R_4} \right)^2 = 2 \left( \frac{1}{R_1^2} + \frac{1}{R_2^2} + \frac{1}{R_3^2} + \frac{1}{R_4^2} \right) $$

Part II - Construction of inner/outer Soddy hyper-spheres

There is an interesting side-product of the proof in Part I. $\beta_k$ is determined up to an overall scaling factor. If we normalize $\beta_k$ such that $B = 4$, $(*2)$ and $(*3)$ together allows us to derive an explicit expression for $\beta_j$

$$\beta_j = \frac{1}{R_j^2} - \left(\frac{1}{d} \sum_{k=1}^n \frac{1}{R_k}\right)\frac{1}{R_j}\tag{*5}$$

We can use this relation to construct the inner and outer Soddy hyper-spheres.

Assume we already have $n-1 = d+1$ hyper-spheres touching among themselves. The inner Soddy hyper-sphere is the sphere outside all these $n-1$ spheres and yet touching all of them.
Let $\vec{x}_k$ and $r_k$ be the center and radius for the $k^{th}$ hyper-sphere for $1 \le k < n$. Let $\vec{x}_{in}$ and $r_{in}$ be the center and radius of the inner Soddy hyper-sphere. If we let $$\vec{x}_n = \vec{x}_{in}\quad\text{ and }\quad R_k = \begin{cases}r_k,& 1 \le k < n\\ r_{in},& k = n\end{cases}$$ discussions in Part I tell us

$$\left( \frac{1}{r_{in}} + \sum_{k=1}^{n-1} \frac{1}{r_k} \right)^2 = d \left( \frac{1}{r_{in}^2} + \sum_{k=1}^{n-1} \frac{1}{r_k^2} \right)\tag{*6a}$$

We can use this to determine $r_{in}$. If the $n-1$ points $\vec{x}_1, \vec{x}_2, \ldots, \vec{x}_{n-1}$ are in general position, i.e. they are vertices of a non-degenerate $d$-simplex, the $d$ vectors $\vec{x}_2 - \vec{x}_1, \ldots, \vec{x}_{n-1} - \vec{x}_1$ will be linearly independent. This implies there exists $d$ coefficients $\gamma_2, \gamma_3, \ldots, \gamma_{n-1}$ such that

$$\vec{x}_{in} - \vec{x}_1 = \gamma_2 (\vec{x}_2 - \vec{x}_1) + \ldots + \gamma_{n-1} ( \vec{x}_{n-1} - \vec{x}_1 )$$

A consequence of this is $\beta_n \ne 0$. This means we can use $(*5)$ and the relation $\sum\limits_{k=1}^n \beta_k \vec{x}_k = \vec{0}$ to compute the center $\vec{x}_{in}$ of the inner Soddy hyper-sphere.

For the outer Soddy hyper-sphere. It is a sphere that contains the original $n-1$ spheres and touching each of them. Let $\vec{x}_{out}$ and $r_{out}$ be the center and radius of the outer Soddy hyper-sphere. The touching condition now takes the form:

$$\begin{array}{ccccl} |\vec{x}_{out} - \vec{x}_j | &=& | r_{out} - r_j |\quad & \text{ for }\quad & 1 \le j < n\\ |\vec{x}_i - \vec{x}_j | &=& | r_i + r_j | \quad & \text{ for }\quad & 1 \le i < j < n \end{array} $$ Once again, if we let $$\vec{x}_n = \vec{x}_{out}\quad\text{ and }\quad R_k = \begin{cases}r_k,& 1 \le k < n\\ -r_{out},& k = n\end{cases}$$ we can repeat discussions in Part I to obtain $$\left( -\frac{1}{r_{out}} + \sum_{k=1}^{n-1} \frac{1}{r_k} \right)^2 = d \left( \frac{1}{r_{out}^2} + \sum_{k=1}^{n-1} \frac{1}{r_k^2} \right)\tag{*6b}$$ We can use this to determine $r_{out}$. Once again, if $\vec{x}_1,\ldots,\vec{x}_{n-1}$ are in general position, we will find $\beta_n \ne 0$. As a result, we can use $(*5)$ to compute $\vec{x}_{out}$, the center of the outer Soddy hyper-spheres, from the remaining centers.

If one compare $(*6a)$ and $(*6b)$, they are very similar, $r_{in}$ and $-r_{out}$ are the two roots of the same equation in $R$.

$$\left( \frac{1}{R} + \sum_{k=1}^{n-1} \frac{1}{r_k} \right)^2 = d \left( \frac{1}{R^2} + \sum_{k=1}^{n-1} \frac{1}{r_k^2} \right)$$

If the two roots of this equation has different sign, the positive root will be the inner Soddy radius $r_{in}$, the negative root will be $-r_{out}$, the negative of the outer Soddy radius.


In high-school, my geometry teacher mentioned the Soddy-Gosset theorem and I had been curious about its proof ever since. About 5 years ago, I decided to attempt to prove it and came up with a nice proof and a corollary that, according to my investigation, was unknown, at least in the form that I have.

The first section below is taken from a paper I wrote at that time. The later sections are a clarification of the later parts of that paper.


Tangent Spheres

Suppose that the sphere with radius $r_1$ and center $c_1$, $(c_1,r_1)$, and the sphere with radius $r_2$ and center $c_2$, $(c_2,r_2)$, are externally tangent. Then we have $$ \begin{align} |c_1-c_2|^2&=(r_1+r_2)^2\tag{1a}\\ |c_1|^2-2c_1\cdot c_2+|c_2|^2&=r_1^2+2r_1r_2+r_2^2\tag{1b} \end{align} $$ For a given sphere, $(c,r)$, the radius of the inverted sphere, $\bar{r}$, satisfies $$ \begin{align} \frac{|c|^2}{r^2}-1&=\frac1{r\bar{r}}\tag{2a}\\ |c|^2-r^2&=\frac r{\bar{r}}\tag{2b} \end{align} $$ We don't actually use any of the properties of spherical inversion, so we can simply take $(2)$ as a definition of $\bar{r}$.

Thus, if we subtract $r_1^2+r_2^2$ from both sides of $\mathrm{(1b)}$ and use equation $\mathrm{(2b)}$, we get $$ \begin{align} \frac{r_1}{\bar{r}_1}+\frac{r_2}{\bar{r}_2}-2c_1\cdot c_2&=2r_1r_2\tag{3a}\\ \frac{c_1}{r_1}\cdot\frac{c_2}{r_2}-\frac{1}{2}(\frac{1}{\bar{r}_1r_2}+\frac{1}{r_1\bar{r}_2})&=-1\tag{3b} \end{align} $$ Furthermore, from equation $\mathrm{(2a)}$, we get $$ \frac{|c|^2}{r^2}-\frac{1}{r\overline{r}}=1\tag{4} $$ Therefore, if, for a sphere in $\mathbb{R}^n$, $(c,r)$, we define the row vector in $\mathbb{R}^{n+2}$ $$ v=\left\langle\frac cr,\frac1r,\frac1{\bar{r}}\right\rangle\tag{5} $$ and make a square matrix, $V$, from $n+2$ of these row vectors for $n+2$ mutually tangent spheres, then we get from equation $\mathrm{(3b)}$, for the off-diagonal elements, and equation $(4)$, for the diagonal elements, that $$ V\left[\begin{array}{c c c} I_n & 0_{n\times 1} & 0_{n\times 1} \\ 0_{1\times n} & 0 & -\frac{1}{2} \\ 0_{1\times n} & -\frac{1}{2} & 0 \end{array}\right]V^T = 2I_{n+2}-C_{n+2}\tag{6} $$ where $C_{n+2}$ is an $(n{+}2)\times(n{+}2)$ matrix of all $1$s. Because $C_{n+2}^2=(n+2)C_{n+2}$, it is easily verified that $$ (2I_{n+2}-C_{n+2})^{-1}=\frac{1}{2}I_{n+2}-\frac{1}{2n}C_{n+2}\tag{7} $$ Invert equation $(6)$ and move $V$ to the right side to get $$ \left[\begin{array}{c c c} I_n & 0_{n\times 1} & 0_{n\times 1} \\ 0_{1\times n} & 0 & -2 \\ 0_{1\times n} & -2 & 0 \end{array}\right] = V^T(\frac{1}{2}I_{n+2}-\frac{1}{2n}C_{n+2})V\tag{8} $$ Looking at the $(n+1,n+1)$ element of the matrix equation $(8)$, we get the equation $$ n\sum_{k=1}^{n+2}\frac{1}{r_k^2}=\left(\sum_{k=1}^{n+2}\frac{1}{r_k}\right)^2\tag{9} $$ If one of the spheres encompasses the rest, so that all the other spheres are internally tangent to it, then equation $\mathrm{(1a)}$ becomes $|c_1-c_2|^2=(r_1-r_2)^2$. If we consider the radius of the encompassing sphere to be negative, then this becomes $|c_1-c_2|^2=(r_1+r_2)^2$ as before. Therefore, we will consider the radius of an encompassing sphere to be negative.

Thus, we have shown

The Soddy-Gosset Theorem: Given $n+2$ mutually tangent spheres in $\mathbb{R}^n$, their radii, $\{r_k:1\le k\le n+2\}$ satisfy equation $(9)$ if the radius of an encompassing sphere is considered negative.

Furthermore, if we look at the $(j,n+1)$ elements of the matrix equation $(8)$, for $1\le j\le n$, we get the vector equation $$ n\sum_{k=1}^{n+2}\frac{c_k}{r_k^2}=\sum_{k=1}^{n+2}\frac{c_k}{r_k}\;\sum_{k=1}^{n+2}\frac{1}{r_k}\tag{10} $$ If we divide equation $(10)$ by equation $(9)$, we get $$ \frac{\displaystyle\sum_{k=1}^{n+2}\frac{c_k}{r_k^2}}{\displaystyle\sum_{k=1}^{n+2}\frac{1}{r_k^2}} =\frac{\displaystyle\sum_{k=1}^{n+2}\frac{c_k}{r_k}}{\displaystyle\sum_{k=1}^{n+2}\frac{1}{r_k}}\tag{11} $$ Thus, we have derived the following

Corollary: Given $n+2$ mutually tangent spheres in $\mathbb{R}^n$ with radii $\{r_k:1\le k\le n+2\}$, the mean of the centers weighted by $\frac1{r_k^2}$ is equal to the mean of the centers weighted by $\frac1{r_k}$.


Conclusions from Soddy-Gosset

Given the radii of $n+1$ mutually tangent spheres in $\mathbb{R}^n$, the Soddy-Gosset Theorem allows us to compute the radius of two possible spheres that are tangent to the other $n+1$. That is, we can solve $(9)$ for $\frac1{r_{n+2}}$: $$ n\left(\color{#C00000}{\frac1{r_{n+2}^2}}+\sum_{k=1}^{n+1}\frac1{r_k^2}\right) =\color{#C00000}{\frac1{r_{n+2}^2}}+2\color{#C00000}{\frac1{r_{n+2}}}\sum_{k=1}^{n+1}\frac1{r_k}+\left(\sum_{k=1}^{n+1}\frac1{r_k}\right)^2\tag{12a} $$ $$ (n-1)\color{#C00000}{\frac1{r_{n+2}^2}}-2\color{#C00000}{\frac1{r_{n+2}}}\sum_{k=1}^{n+1}\frac1{r_k}+n\sum_{k=1}^{n+1}\frac1{r_k^2}-\left(\sum_{k=1}^{n+1}\frac1{r_k}\right)^2=0\tag{12b} $$ applying the quadratic formula to $(12)$, we get $$ \frac{n-1}{r_{n+2}}=\sum\limits_{k=1}^{n+1}\frac1{r_k}\pm\sqrt{n}\left[\left(\sum\limits_{k=1}^{n+1}\frac1{r_k}\right)^2-(n-1)\sum\limits_{k=1}^{n+1}\frac1{r_k^2}\right]^{1/2}\tag{13} $$ Note that the discriminant in $(13)$ vanishes precisely when the first $n+1$ spheres satisfy the Soddy-Gosset Theorem on $\mathbb{R}^{n-1}$. This implies that $\{c_k:1\le k\le n+1\}$ lie in an $n-1$ dimensional hyperplane, and the two possibilities for $(c_{n+2},r_{n+2})$ are reflections of each other across this hyperplane.


Conclusions from the Corollary

Given the radii of $n+2$ mutually tangent spheres and the centers of the first $n+1$ of them, the Corollary allows us to determine $c_{n+2}$, except in the case where the coefficients of $c_{n+2}$ cancel in $(11)$. That is, when $$ \frac{\dfrac1{r_{n+2}^2}}{\displaystyle\sum_{k=1}^{n+2}\frac{1}{r_k^2}} =\frac{\dfrac1{r_{n+2}}}{\displaystyle\sum_{k=1}^{n+2}\frac{1}{r_k}}\tag{14} $$ which is equivalent to $$ \begin{align} \frac1{r_{n+2}} &=\frac{\displaystyle\sum_{k=1}^{n+2}\frac1{r_k^2}}{\displaystyle\sum_{k=1}^{n+2}\frac1{r_k}} =\frac{\displaystyle\sum_{k=1}^{n+1}\frac1{r_k^2}}{\displaystyle\sum_{k=1}^{n+1}\frac1{r_k}}\tag{15a,15b}\\ &=\frac1n\sum_{k=1}^{n+2}\frac1{r_k} =\frac1{n-1}\sum_{k=1}^{n+1}\frac1{r_k}\tag{15c,15d} \end{align} $$ Explanation:
$\mathrm{(15a)}$: restatement of $(14)$
$\mathrm{(15b)}$: multiply $\mathrm{(15a)}$ by $\sum\limits_{k=1}^{n+2}\frac1{r_k}$, subtract $\frac1{r_{n+2}^2}$, then divide by $\sum\limits_{k=1}^{n+1}\frac1{r_k}$
$\mathrm{(15c)}$: apply $(9)$ to $\mathrm{(15a)}$
$\mathrm{(15d)}$: multiply $\mathrm{(15c)}$ by $n$, subtract $\frac1{r_{n+2}}$, then divide by $n-1$

Note that $\mathrm{(15b)}$ and $\mathrm{(15d)}$ say that the first $n+1$ spheres satisfy the Soddy-Gosset Theorem in $\mathbb{R}^{n-1}$. Again, this implies that $\{c_k:1\le k\le n+1\}$ lie in an $n-1$ dimensional hyperplane.

Furthermore, $\mathrm{(15d)}$ matches $(13)$ when the discriminant vanishes.


Locating $\boldsymbol{c_{n+2}}$ in the Exceptional Case

By rotation and translation, we can arrange that the $n-1$ dimensional hyperplane containing $\{c_k:1\le k\le n+1\}$ is $x_n=0$. That is, if $1\le k\le n+1$, $$ c_{k,n}=0\tag{16} $$ and $c_{n+2,n}$ is the perpendicular distance of $c_{n+2}$ above or below that hyperplane. For clarity of presentation, let $\hat{p}$ be the projection of $p$ onto the hyperplane $x_n=0$; that is, $\hat{p}_k=p_k$ for $1\le k\le n-1$ and $\hat{p}_n=0$.

Looking at the $(j,n)$ elements of matrix equation $(8)$ for $1\le j\le n-1$, we get the vector equation $$ \begin{align} \frac{\hat{c}_{n+2}}{r_{n+2}} &=\frac1n\sum_{k=1}^{n+2}\frac{\hat{c}_k}{r_k}\tag{17a}\\ &=\frac1{n-1}\sum_{k=1}^{n+1}\frac{c_k}{r_k}\tag{17b} \end{align} $$ which, in light of $\mathrm{(15d)}$, gives the vector equation $$ \hat{c}_{n+2}=\frac{\displaystyle\sum_{k=1}^{n+1}\frac{c_k}{r_k}}{\displaystyle\sum_{k=1}^{n+1}\frac1{r_k}}\tag{18} $$ Looking at the $(n,n)$ element of matrix equation $(8)$ we get $$ \begin{align} 1 &=\frac12\frac{c_{n+2,n}^2}{r_{n+2}^2}-\frac1{2n}\frac{c_{n+2,n}^2}{r_{n+2}^2}\\ &=\frac{n-1}{2n}\frac{c_{n+2,n}^2}{r_{n+2}^2}\tag{19} \end{align} $$ which gives us $$ c_{n+2,n}=\pm r_{n+2}\sqrt{\frac{2n}{n-1}}\tag{20} $$ We can now remove the rotation, and use $(18)$ and $(20)$ to get that $$ c_{n+2}=\frac{\displaystyle\sum_{k=1}^{n+1}\frac{c_k}{r_k}}{\displaystyle\sum_{k=1}^{n+1}\frac1{r_k}}\pm ur_{n+2}\sqrt{\frac{2n}{n-1}}\tag{21} $$ or equivalently, using Soddy-Gosset and the Corollary in $\mathbb{R}^{n-1}$ and $\mathrm{(15b)}$, $$ \frac{c_{n+2}}{r_{n+2}}=\frac{\displaystyle\sum_{k=1}^{n+1}\frac{c_k}{r_k}\frac1{r_k}}{\displaystyle\sum_{k=1}^{n+1}\frac1{r_k}}\pm u\sqrt{\frac{2n}{n-1}}\tag{22} $$ where $u$ is a unit vector perpendicular to the plane containing $\{c_k:1\le k\le n+1\}$.