Calculating the density of optimally packed spheres

[IridescentRain]

Hello.

I wish to calculate the density of optimally packed spheres (i.e. the fraction of space which is occupied by the spheres). This number is approximately 0.74, but I know this because I have looked on the net, not because I have calculated it. I know the essentials of how to calculate it, but I am stuck in my calculations.

All of the following analysis and calculations were done by me; I've checked my results (the density of optimally packed circles, a regular tetrahedron's height, etc) with the known expressions found on the internet and in books.

Circles

I'll first derive the density of optimally packed circles on a plane, which is more illustrative and easier to understand, and then I'll proceed to the problem of spheres in space.

Consider an infinite array of optimally-packed identical circles on a plane:

We may trace an (equilateral) triangle with its vertices at the centres of three mutually adjacent circles. If we do this with all the circles, then the (infinite) set of all such triangles covers the entire plane, so the fraction of a single triangle occupied by circle fragments is equal to the fraction of the plane occupied by all the circles. Let us call this the "volume fraction" of the circles and denote it by $f$.

Let $a$ be the radius of the circles and $D=2\,a$ be the length of the triangle's sides.

Because the triangle is equilateral, its angles all equal $\pi/3$ radians. Therefore, the fraction of each circle that lies within the triangle is $1/6$. The volume fraction is then
$$f=3\frac{\pi\,a^2/6}{A_T}=\frac{\pi\,a^2}{2\,A_T},$$where $A_T$ is the triangle's area and the 3 multiplying the first fraction comes from the fact that there are 3 circles which lie at least partially within the triangle.

By Heron's formula,
$$A_T=\sqrt{s\,(s-D)^3},$$where $s=3\,D/2$. Therefore,
$$A_T=\sqrt{\frac{3}{2}D\,\left(\frac{3}{2}-1\right)^3\,D^3}=\frac{\sqrt{3}}{4}D^2,$$whereby
$$f=\frac{\pi\,a^2}{(\sqrt{3}/2)\,D^2}=\frac{2\,\pi}{\sqrt{3}}\frac{a^2}{D^2}=\frac{\pi}{2\,\sqrt{3}}\approx0.9069.$$
Spheres

We generalise the problem to spheres as follows: consider optimally-packed identical spheres, and take a regular tetrahedron with its vertices at the centres of four mutually adjacent spheres. The (infinite) set of all such tetrahedra again fills all of space, so, if $V_S$ equals the volume of each sphere which lies within the tetrahedron and $V_T$ is the tetrahedron's volume, we have
$$f=4\frac{V_S}{V_T}.$$

Since the tetrahedron is regular, its faces are equilateral triangles. Let $a$ be the radius of the spheres, and let $D=2\,a$ once again be the length of the tetrahedron's edges.

The implicit equation of a sphere of radius $a$ is, of course,
$$x^2+y^2+z^2-a^2=0.$$
Consider a regular tetrahedron with edge length $a'$ placed such that one of its vertices lies at the origin and the centroid of the face opposite this vertex lies on the $z$ axis. The position of said face's centroid is then of the form $(0,0,z_0)$. We proceed to calculate $z_0$:

The top face of the tetrahedron (that is, the face which lies at a height $z_0$) is, of course, an equilateral triangle with side length $a'$. Let $B$ be the distance from its centroid to any of its vertices and $b$ be the distance from its centroid to the middle point of any of its sides:

If the coordinates of the triangle's vertices are $(0,y_1),(-x_2,y_2),(x_2,y_2)$ (omitting the $z$ coordinate for clarity), then, by symmetry,
$$y_1^2=x_2^2+y_2^2=B^2.$$Also,
$$x_2^2+(y_1-y_2)^2=a'^2.$$
Now, $x_2=a'/2$, whereby the above equations become
$$y_1^2=\frac{a'^2}{4}+y_2^2,$$$$\frac{a'^2}{4}+(y_1-y_2)^2=a'^2.$$The latter equation implies that
$$y_1=\sqrt{\frac{3}{4}a'^2}+y_2=\frac{\sqrt{3}}{2}a'+y_2.$$Substituting this in the former equation, we obtain
$$\left(\frac{\sqrt{3}}{2}a'+y_2\right)^2=\frac{3}{4}a'^2+\sqrt{3}\,a'\,y_2+y_2^2=\frac{a'^2}{4}+y_2^2$$$$\Rightarrow\sqrt{3}\,a'\,y_2=-\frac{1}{2}a'^2$$$$\Rightarrow y_2=-\frac{1}{2\,\sqrt{3}}a'$$$$\Rightarrow B=\sqrt{\left(\frac{a'}{2}\right)^2+\left(\frac{1}{2\,\sqrt{3}}a'\right)^2}=\sqrt{\frac{1}{4}+\frac{1}{12}}\,a'=\frac{a'}{\sqrt{3}}.$$
To calculate $b$, we divide the triangle into four identical equilateral triangles with side length $a'/2$:

The central triangle's centroid lies at the same location as that of the original triangle by symmetry. $b$ is then the distance from the centroid of an equilateral triangle with side length $a"=a'/2$ to any of that triangle's vertices. Using the expression we obtained for $B$, we have
$$b=\frac{a"}{\sqrt{3}}=\frac{a'}{\sqrt{12}}.$$
The height of the original triangle is
$$h=B+b=\left(1+\frac{1}{2}\right)\frac{a'}{\sqrt{3}}=\frac{\sqrt{3}}{2}a'.$$
Now consider the entire tetrahedron (with edge length $a$). If we place a new coordinate system $x'y'z'$ so that the $z'$ axis lies along the centre of one of the three faces which are not horizontal in the $xyz$ coordinate system, then the "height" $z'_0$ along $z'$ of that face is
[tex]z'_0=\frac{\sqrt{3}}{2}a.

The (non-equilateral) triangular region enclosed by the $z$ and $z'$ axes and the tetrahedron's top face is shown in the following figure:

Notice that
$$\sin(\varphi)=\frac{a/2\,\sqrt{3}}{z'_0}=\frac{a}{2\,\sqrt{3}}\frac{2}{\sqrt{3}\,a}=\frac{1}{3},$$$$\tan(\varphi)=\frac{a/2\,\sqrt{3}}{z_0}.$$Therefore,
$$z_0 = \frac{a}{2\,\sqrt{3}\,\tan(\varphi)} = \frac{a}{2\,\sqrt{3}\,\tan \left[ \arcsin(1/3) \right]} = \frac{a\,\sqrt{1-(1/3)^2}}{2\,\sqrt{3}/3} = \frac{a\,\sqrt{8/9}}{2/\sqrt{3}} = \sqrt{\frac{2}{3}}\,a.$$
Having finally calculated the tetrahedron's height $z_0$, we may consider the tetrahedron to be composed of a collection of equilateral-triangular prisms (with differential height $\text{d}z$ with side length $a'(z)$ stacked on top of each other. $a'$ is a linearly increasing function of $z$ with the following boundary conditions:
$$a'(0)=0,$$$$a'(z_0)=a'\left(\sqrt{\frac{2}{3}}a\right)=a.$$Therefore,
$$a'=\sqrt{\frac{3}{2}}z.$$
This expression for $a'$ in terms of $z$ --along with the expressions at the beginning of this section which relate $x_2$, $y_1$ and $y_2$ to $a'$-- may be used to construct an implicit equation for the tetrahedron (which I will not do at the moment). The volume $V_S$ of the portion of each sphere which lies within the ($D$-sided) tetrahedron is, then,
$$V_S=\int\int\int\,\text{d}x\,\text{d}y\,\text{d}z,$$where the integration limits are given by the intersection of the regions described by the implicit equations for the sphere and the (infinite, with side length $a'(z)$) tetrahedron.

My question is how to write these limits. The limits of the integral over $z$ are trivial: $z$ runs from $0$ to $a$. The limits for $x$ and $y$, however, are absolutely not trivial.

Of course, I may write the integration limits as $\min(A,B)$, where $A$ and $B$ are the corresponding integration limits for the sphere and the tetrahedron, respectively, but this is useless because it does not allow me to calculate $V_S$.

Notes

I know there is an easier way to calculate $f$ which uses a hexagonal prism and 10 or so spheres (as opposed to a tetrahedron and 4 spheres). There's a good reason for this, but I shall not go into it in detail here unless asked.

The procedure I have used is probably not the simplest one possible, but it also has a reason: since it ultimately requires nothing more than the implicit equations for a sphere (the packed particle in question) and a tetrahedron (the space-filling figure which allows me to take a finite number (4) of particles rather than having to deal with an infinite amount of them), whichever expression for $f$ I obtain at the end is trivially generalisable to particles of any shape (provided I can describe said shape via an implicit equation).

Thank you! :)

 

[mfb]

You cannot fill all space with tetrahedrons. The ideal sphere packing has other shapes inside, too.

 

[IridescentRain]

In fact, I can; the "cannonball arrangement" does just that. Take any four mutually adjacent spheres (which, of course, form a tetrahedron). Now take any three of those spheres and a new sphere which is mutually adjacent to all three; they will form a second tetrahedron which has a face in common with the first tetrahedron (and no overlap, of course). You can continue doing this until you fill all of three-dimensional space.

There are other arrangements which are also optimal and, as you say, have other shapes inside.

Regardless of whether I can fill all of space with tetrahedra, though, my problem concerns how to express the integration limits of a volume integral which takes place in the intersection of two regions which I can describe using implicit equations:
$$f_1(x,y,z)=0,$$$$f_2(x,y,z)=0.$$

 

[mfb]

Consider the foremost ball in the second layer from the bottom, and the two spheres behind that: They have a sphere on top, but no sphere below them (touching all 3 of them). The point in the middle, below those 3 spheres is the center of a triangle of 6 spheres (of the bottom layer), which is not the position of a sphere.

Any 4 spheres touching each other form a tetrahedron, of course, but you cannot fill space by adding just those tetrahedron shapes everywhere.

Regardless of whether I can fill all of space with tetrahedra, though, my problem concerns how to express the integration limits of a volume integral which takes place in the intersection of two regions which I can describe using implicit equations:

I'm not sure if that integral is relevant (or sufficient) for the original question of the optimal sphere packing.

 

[IridescentRain]

Ah, you're right. I misinterpreted the picture. Thanks.

I'll check what the unit cell for the "cannonball arrangement" is and use that instead of the tetrahedron. I suppose I'll lose generalisability, though...