Measure of overlap between two arbitrary ellipses on a sphere.

Various packing problems and simulations of hard and soft interacting particles, such as microscopic models of nematic liquid crystals, reduce to calculations of intersections and pair interactions between ellipsoids. When constrained to a spherical surface, curvature and compactness lead to non-trivial behaviour that finds uses in physics, computer science and geometry. A well-known idealized isotropic example is the Tammes problem of finding optimal non-intersecting packings of equal hard disks. The anisotropic case of elliptic particles remains, on the other hand, comparatively unexplored. We develop an algorithm to detect collisions between ellipses constrained to the two-dimensional surface of a sphere based on a solution of an eigenvalue problem. We investigate and discuss topologically distinct ways two ellipses may touch or intersect on a sphere, and define a contact function that can be used for construction of short- and long-range pair potentials.

Introduction

It comes as no surprise that packing of ellipses and ellipsoids is a very thoroughly researched topic that appears in many different fields of research, both in experimental realizations and in numerical models used to study them. Ellipsoids appear in Gay–Berne (anisotropic Lennard–Jones) models [1] of liquid crystals as a coarse-grained replacement for the full molecular structure [2-4], in colloidal dispersions with an anisotropic dispersed phase [5-7], and in granular and jammed matter [8-11], where random and optimal packings are of particular interest [12,13]. All these examples are, however, Euclidean—yet many experimental systems call for a confinement of particles to a curved surface, often that of a sphere. Recent examples include packings of rods [14] and ellipsoids [15], spherocylinder simulations of nematics [16] and proteins adsorbed on vesicles [17,18]. This calls for an adaptation of ellipse–ellipse intersection algorithms for use on a spherical surface. Such an algorithm would also allow answering the question of optimal packing: while the well-researched Tammes problem [19,20] considers optimal packings of circles on a sphere, a generalization from circles to ellipses of arbitrary aspect ratios can provide us with the packing fraction for hard ellipses, which so far remains an open question. Furthermore, an algorithm which can be applied to ellipses of different sizes and aspect ratios opens up the possibility to consider polydisperse systems.

The bread-and-butter of computing ellipse–ellipse interactions lies in detecting collisions and overlaps in simulations of hard particles [21], and, for long-range interactions, measuring the closest distance between them [22]. One of the widely used and cited algorithms developed by Perram et al. [23,24] has been used, optimized and adapted in numerous ways and for various applications—in two dimensions (for ellipses) [25,26], three dimensions (for ellipsoids) [27-29] and was even generalized to hyperellipsoids [30]. However, all these algorithms are limited to Euclidean space and cannot be applied to the spherical case without modification.

In this work, we present a new algorithm that tackles the previously unsolved question of computing the distance and detecting overlap of ellipses confined to the two-dimensional surface of a sphere. Spherical confinement poses interesting challenges to the algorithm. Stretching is not a linear operation on a sphere, and two ellipses can interact in topologically different ways—if they interact at all. These situations differ strongly from the Euclidean case. We explain the intricacies of the spherical ellipse–ellipse interaction with examples, discuss the performance of the numerical algorithm and conclude by showing a few packing solutions.

Numerical algorithm

Problem formulation

First, we define what constitutes an ellipse on the surface of a sphere. We adopt the conventional definition of an ellipse as the set of points with a constant sum of distances to the foci. To generalize it to a sphere, we require a constant sum of geodesic distances to the foci : . Without loss of generality, we can set . To convert the trigonometric relations into an algebraic form, we work with the cosine of the focal property , which introduces an extraneous solution—another ellipse at the opposite side of the sphere due to . With further manipulations, we obtain a quadratic form representing an elliptic cylinder, which is oriented in the plane. However, on the unit sphere, the set of solutions is invariant to addition of any multiple of , which gives a whole family of quadratic forms that specify the same ellipse pair. The most natural representation among these is an elliptic cylinder in the plane with zero term, which simply projects a planar ellipse onto the sphere along the axis through the centre of the spherical ellipse

We will thus define a spherical ellipse as an intersection of the unit sphere and an elliptical cylinder represented as a degenerate positive semidefinite quadratic form where is a rotation matrix that will not be explicitly needed. We assume from now on that is a given quantity which can be computed from any representation of the ellipses, such as by rotating the focal representation (equation 2.2), or from centre vectors and major semiaxis directions. Looking for an intersection of two arbitrary spherical ellipses is therefore equivalent to looking for an intersection of two quadratic forms and the unit sphere Quadratic forms are invariant to inversion and produce a pair of antipodal ellipses when intersected with the unit sphere. This poses an additional challenge for the collision detection algorithm, as we must specify which ellipse is the correct one and which collisions to ignore. The correct ellipses can be specified by vectors and corresponding to the centres of the ellipses—signed eigenvectors corresponding to the zero eigenvalue of the quadratic forms and . The dot product between the ellipse centre and any point on the ellipse is positive for the correct ellipse and negative for the antipode.

Unlike in the Euclidean case, scaling the semiaxes of the quadratic form has an important effect on the topology of its intersection with the unit sphere. When the semiaxes are small compared with the radius of the sphere, the ellipses are similar to Euclidean ellipses. If the semiaxes are scaled to be comparable with the sphere radius, the apexes become sharper and converge to a ‘lemon wedge’ shape in the limit where the large semiaxis of the quadratic form matches the sphere radius. In this configuration, the antipodes touch at two ‘poles’, forming two intersecting great circles. Beyond this size, the intersection with the sphere splits again into a new pair of ellipses, but now their centres are directed along the shorter of the quadratic form semiaxes. At this crossover, the former antipodal pair recombines, and no longer corresponds to elliptical particles centred at . These cases with inverted ellipses will play a role in our theoretical analysis, but have no physical significance.

The goal of our algorithm is to detect when two ellipses are tangent or overlapping by defining a contact function and to obtain the contact point . If forces at the contact point are required, the direction of the force should be along the normal to the ellipse, which is given by the gradient of the quadratic form (magnitudes can be normalized—here we halve the expression to simplify notation) From the force and the intersection point, we can also compute torques acting on the ellipse, which is useful for molecular dynamics simulations.

Solving for ellipse contacts

Following the same steps as Perram and Wertheim [11,23,24], we define a linear interpolation of the quadratic forms, with the parameter , so that on the whole sphere. Solution for contacts of spherical ellipses, i.e. level sets and , is based on finding the minimal values of this interpolation at each . Constraining the solutions to the unit sphere, the problem can be restated in terms of finding the stationary points of the Lagrangian function Equation reduces to solving the eigenvalue problem , with solutions and eigenvectors that satisfy (for ). Let be the eigenvector corresponding to the smallest eigenvalue . Plugging back into expression (2.8), we get the minimum value of on the sphere,

Consider that the value of the quadratic form is always greater than 1 in the part of the sphere that is outside both ellipses, as it is an interpolation of two values greater than 1. As is varied from 0 to 1, will trace a continuous path on the sphere from to where . If the ellipses do not overlap, this means that will have to cross the region outside both ellipses for some , where consequently . On the other hand, for overlapping ellipses, will always be smaller than 1 in the intersection region, meaning that for all . It follows from equation (2.10) that ellipses and intersect on the unit sphere if and only if the smallest eigenvalue of never exceeds 1 on the interval . We denote this extremum and the corresponding eigenvector , Positive definiteness ensures there are always three non-negative real eigenvalues, corresponding to the casus irreducibilis of the cubic equation, which is solvable in closed form through trigonometry. To find the maximum , any one-dimensional maximization algorithm can be used, such as the golden section search. We can rely on this function being anticonvex with a single maximum, which ensures reliable and fast convergence.

To gain additional understanding of the relation between eigenvalue extrema and ellipse contacts, one can consider that the geometric representation of the unconstrained three-dimensional level set is a generic ellipsoid (or possibly a degenerate elliptical cylinder when one eigenvalue of is zero). The eigenvalues of correspond to inverse squares of its semiaxes. This ellipsoid is thus completely contained within the unit sphere if all its eigenvalues are greater than 1 and intersects the unit sphere if this is not the case. For non-overlapping ellipses, where , the level set , constrained to the surface of the sphere, will therefore be empty for some . This supports the fact that the space of allowed level set locations on the sphere is discontinuous if the ellipses do not overlap and the level set therefore cannot evolve continuously from to as is increased from 0 to 1. Conversely, if the ellipses intersect, the level set on the sphere is always non-empty as , with the intersection points of and a part of the level set for each . Examples of disjoint, touching and overlapping ellipses, and the level sets of , are shown in figure 1a–c.

Figure 1

Cases of (a) disjoint, (b) touching and (c) overlapping ellipses on a sphere, with ellipses stretched to achieve tangency shown by dashed coloured lines. Level sets of (equation (2.8)) at different are shown by dashed grey lines. Black dashed lines represent the lines described by the eigenvector corresponding to the smallest eigenvalue of when runs from 0 to 1, and marks the intersection point found when the eigenvalue is maximized. (d) Smallest eigenvalue with respect to for the three cases in (a–c), showing that the smallest eigenvalue exceeds 1 when the ellipses are disjoint. (Online version in colour.)

The value of has a clear geometric meaning: if we observe the intersection with a sphere instead of the unit sphere, the ellipses and touch at a single point of tangency, given by the appropriately scaled eigenvector . Scaling the system back to the unit sphere by a factor of , we see that is the factor by which the orthogonal projected area of both ellipses must be grown to make them tangent (scaling the semiaxes by ). Values of signify non-overlapping ellipses that become tangent when grown, and values of overlapping ellipses that become tangent when shrunk. This property makes an appropriate choice for a contact function, with the same meaning it has in the Euclidean case (see the work of Perram and Wertheim [11,23,24]). However, without additional tests, the value of does not distinguish between the two antipodal ellipses represented by the same quadratic form and thus signals an overlap even when the ellipses in question are on the opposite sides of the sphere. For a usable algorithm, collisions with the antipodes of the ellipses represented by the quadratic forms must be ignored. This is handled in the following section.

Solution branches and secondary contacts

Points of tangency of ellipses on the sphere can be defined in terms of the full intersection set of two elliptical cylinders in three dimensions, The ellipses, obtained as intersections of and with the sphere of radius , are intersecting at points on at radius and are tangent in critical points on with locally extremal distance from the origin. The maximized smallest eigenvalue , which we derived in the previous section, simply corresponds to the critical point of farthest from the origin; but this is just one of the critical points.

Degenerate cases aside, the set consists of an antipodal pair of two disjoined loops. Each loop can have at most four critical points—two with locally maximal and two with locally minimal distance to the origin, corresponding to four values (figure 2a). Depending on the relative orientation and size of the ellipses, there may be only two critical points, (figure 2b). At the transition between these two regimes, the critical points and merge into an inflection point before disappearing. In other borderline cases with zero measure, the intersection set may be a ‘basket’ with two fourfold junctions, or may have whole arcs at constant distance from the origin. These can all be understood as limiting cases with degenerate maxima and minima.

Figure 2

(a) A generic intersection set of two obliquely intersecting elliptical cylinders. The intersection consists of two antipodal loops, with two points of maximal distance and two points of minimal distance from the origin. These represent the four tangency cases; only and are relevant for our analysis. (b) In exceptional cases, the two loops might be joined in a four-way junction. (Online version in colour.)

The maximal critical points correspond to the tangency with appearance of two new intersections when ellipses are stretched past the tangency condition. The minimal critical points correspond to the disappearance of intersections when stretching ellipses past the tangency condition. Only the maxima—the critical points —are relevant for detecting ellipse contacts. The remaining two critical points involve inverted ellipses, as they describe points on with locally minimal distance to the origin and are thus closer than at least one of the quadratic form semiaxes.

The antipodal doubling of ellipses means that the tangency at may correspond to the contact with the antipode of the second ellipse, so it might not be the one we are looking for. If there are only two critical points, there is no other possible contact. If there are four critical points, growing the ellipses further makes them touch again at the next locally maximal critical point (). This contact might be between the correct pair of ellipses, or it could be between the same pair of ellipses as the critical point, in which case it is not a candidate for a true contact either.

As already discussed, the maximum of the lowest eigenvalue, , solves for the first contact. The rest of the contacts can also be tied to extrema of the eigenvalues of over . The values and correspond to the minimum and the maximum of the middle eigenvalue, and to the minimum of the largest eigenvalue (figure 3). Unlike the lowest eigenvalue of , which is guaranteed to have a local maximum between and , the remaining eigenvalues can have extrema outside the interval , or none at all. In these cases, constrained minimization returns one of the edge points of the interval.

Figure 3

Eigenvalue spectrum of a generic case, with all four extrema occurring inside the interval . If the first contact is between the antipodes (lower left inset), the true contact and thus the correct value of the contact function is found by the minimum of the second eigenvalue. Observe that the green ellipse (a = 1) still intersects the red (b = 1) antipode ellipse, which we are ignoring. If the first contact is between the correct ellipses, then the lowest eigenvalue is the correct solution—we need information about the correct antipode to test for that. The upper two extrema correspond to inverted contacts (the shaded area lies above the lowest non-zero eigenvalue of and ). (Online version in colour.)

If there are only two critical points on each loop of the intersection manifold , then the middle eigenvalue has no local extrema, neither inside the interval nor anywhere else on the real line, and are undefined. If there are four critical points, the local extrema may lie outside the interval . This corresponds to a second contact between the same pair of ellipses as , meaning that either both critical points signify contact between the true ellipses or both signify contact with the antipode, in which case there is no contact (figure 4). This is convenient, as simply checking for the existence of a minimum of the middle eigenvalue inside the interval includes all cases in which the critical point can constitute a real contact. Finally, if the resulting or exceed any of the eigenvalues of or (which coincide with the non-zero eigenvalues of and ), it signifies a contact where at least one ellipse is inverted. We can test this by finding the minimum non-zero eigenvalue of and , corresponding to the largest semiaxis of the largest ellipse. Critical points that exceed this value, , do not correspond to valid contacts, nor can their values be unambiguously used as an analitical continuation of the contact function, because mixing of antipodes into the inverted ellipse makes the choice between the branches impossible.

Figure 4

Eigenvalue spectrum of a case where the first two contacts are between the same ellipses, signifying that either the first contact is correct or neither of them is (as in the latter case, both contacts are between antipodes). Such situations are characterized by the middle eigenvalue not having a local minimum in the interval . Top two contacts are not depicted. The shaded area corresponds to inverted ellipses. (Online version in colour.)

Contact function

To define a well-behaved contact function to use as a test for ellipse–ellipse intersections, the correct eigenvalue must be selected. This is done with the help of the eigenvectors, which correspond to critical points. We denote with and the eigenvectors corresponding to and , respectively, and and are the true centres of the ellipses on the unit sphere, their signs picking the correct ellipse of the antipodal pair. If the true ellipse collides with the antipode of the second one, the projections of the intersection vector onto the vectors of ellipse centres are of opposite signs, and vice versa. If there is no contact, or the contact is with an inverted ellipse, assigning the value makes the function continuous under variations of the relative position of the ellipses. The full algorithm for computing the contact function is described in algorithm 1.

The contact function , which is according to the above criterion equal to , , or , can be used either to directly detect when ellipses overlap () or to construct a pair potential. Instead of a hard core repulsion, a soft repulsion potential for overlapping cases can be defined based on the value of , such as , , or , the last being a soft potential of finite strength at complete overlap. On the other hand, long-range values of could act as a distance metric, e.g. in a Lennard–Jones-like potential, as they do in Euclidean space [22]. Setting the function to in cases for which the ellipses cannot intersect no matter the stretch factor, ensures a constant potential and zero force on the particles for that entire region, and makes the function well-behaved for use in methods that require a potential (e.g. Monte Carlo methods). Even though there is no correspondence between such an artificially fashioned potential and any physical phenomena we know of, such an academic exercise could provide a reasonable approximation to medium-range behaviour that could match empirical observations in certain physical systems.

Examples

Intersection of unequal circles

The simplest example that can be used for interpretation of the contact function is a pair of unequal circles. Define the following pair of quadratic forms: with and the angular separation of the circle centres. In this case, the extremal eigenvalues (without applying the restriction to ) have a relatively simple closed form, and the contact function can be expressed as The function’s behaviour with respect to is depicted in figure 5 for a few combinations of circle sizes and . At , we have , corresponding to the second contact, as the first contact is with the antipode. We observe that the crossover between the branches is continuously differentiable. However, with the exception of equal circles, we see that the function reaches a maximum at and then goes back to zero at . This part of the contact function corresponds to the second collision also being with the antipode. The collision is internal (non-facing normals), and the interpolation parameter at minimal middle eigenvalue is . In our algorithm, we assign these collisions .

Figure 5

Contact function for circles of different relative radii separated by angle . Transition to the antipodal contact at is continuous; the insets show that the secondary contact is the correct one, while the contact with the antipode (dashed circle outline without infill) is ignored. If the circles are of the same size, the contact function is monotonously increasing, but if they are of different sizes, the decreasing part (shaded below the curve) corresponds to the case when the ellipses cannot be made to touch by stretching, and the corresponding eigenvalue detects the second contact between the ‘wrong’ pair of ellipses. In this region, the value of is set to (horizontal line), which corresponds to the inverse square radius of the larger circle. Parameters and correspond to inverse square radii (see equation (3.1)). (Online version in colour.)

Computational cost of the algorithm

The algorithm itself is fast, as the eigenvalue calculation can be expressed in a closed form, although it uses trigonometric functions that are slower than simple multiplications. One-dimensional minimization and maximization routines are available in any number of numerical libraries. We implemented two such routines, the golden section search (GSS) and the Brent method (GSS with quadratic interpolation), and compared both the numbers of eigenvalue evaluations to achieve the desired accuracy (tolerance of in ) as well as calculation times. In figure 6, we show the results for a pair of ellipses on a unit sphere with major and minor semiaxes and , respectively (aspect ratio ). At a given angular separation , the contact function and the computational cost required to determine it with the Brent method depend on orientations of both ellipses as shown for in figure 6a and 6b for the number of first and second eigenvalue evaluations, respectively. The number of evaluations for mostly lies between 10 and 20, with the exceptions of diagonals with fewer evaluations and two loops with that correspond to the cases near the and crossover. For relatively high aspect ratios , as is the case in the demonstrated example, the second derivative close to the crossover becomes large, which is unfavourable for the Brent minimization. Inside these loops, the second eigenvalue becomes relevant for the contact function, as indicated in figure 6b ( is only evaluated in regions where the eigenvector test fails, see algorithm 1). Additionally, closer to the diagonals, the local minimum of inside the interval disappears and the algorithm returns . Note again that despite the algorithm branch changes, the contact function is continuous in the whole configuration space.

Figure 6

Number of eigenvalue evaluations and contact function calculation times for a pair of ellipses with and . Number of evaluations needed to determine (a) and (b) for the Brent method at angular separation in the whole orientational domain. For , in a large part of the domain where is the correct eigenvalue for determining the contact function (outside of the white contours in (b)). Around this eigenvalue crossover, the number of evaluations is increased. The white contour line in (a) surrounds the region where even with the ET enabled. The red contours in (b) show the border where the contact function transitions to the constant value of . (c) Number of contact evaluations for GSS and Brent line minimizators. The increase in evaluations is a consequence of growing regions where is not the correct eigenvalue. ET (dashed lines) significantly decreases the number of evaluations. The grey area corresponds to distances where the overlap of ellipses depends on their orientations; on the left side of this region, overlap is guaranteed, while there can be no overlap on the right side of the region. (d) Comparison of contact function calculation times for GSS and Brent methods. ET results are relevant only for , as they can only be used to determine overlap/no overlap. (Online version in colour.)

We evaluate the necessary computational cost to determine the contact function both for the GSS and Brent methods. The results with respect to the angular separation are shown in figure 6c, where the value at each distance represents the average number of eigenvalue evaluations over the whole orientational domain ( points, figure 6a,b). The number of evaluations with the GSS method remains (almost) constant for all distances, as a fixed number of interval divisions is necessary to achieve the desired precision. This number is also markedly higher compared with the Brent method, which shows that quadratic interpolation is highly effective for this problem (this could be expected from eigenvalue curves in figures 3 and 4). Note that the number of evaluations is symmetric around , as elliptical cylinder configurations are invariant to coordinate transformation and only the antipode interpretations for the correct/wrong ellipse are exchanged. The number of evaluations does not show this symmetry. At small angular separations, the first eigenvalue will always be the correct one and only for higher does the second eigenvalue evaluation become necessary in parts of the orientational space (figure 6b). These regions become larger as is increased (at some point, they consume the whole orientational domain), which in turn increases the average evaluation numbers.

In some situations, e.g. for simulations of hard particles, the calculation of the exact contact function is not needed. The optimization algorithm can be terminated immediately after a value of is encountered, as that means no overlap. The average number of evaluations with this early termination (ET) condition is shown in figure 6c with dashed lines and leads to a sharp decrease of the necessary calculations in a large part of the plot. As shown in (a) for , more than one evaluation is necessary only inside the white contour which grows/shrinks for smaller/larger distances. Additionally, the grey region in the plot highlights the distances where the overlap appears only for certain ellipse orientations ()—on the left side of this region, ellipses overlap for all orientations and on the right, overlap is not possible as they are too distant and the eigenvalue calculation can be skipped entirely.

Finally, figure 6d shows the average calculation time to evaluate the contact function. The results are on the order of and closely follow the combined number of and evaluations from (c), with the increase in calculation times corresponding to additional evaluations needed to determine the second eigenvalue at larger distances. If ET is enabled, the efficiency of the calculation is significantly improved.

Dense packings of spherical ellipses

To demonstrate the use case of the proposed algorithm in multiparticle simulations, we calculated dense packings of spherical ellipses with for both monodispersed and bidispersed systems (figure 7). We employed an energy minimization-based approach similar to the scheme used by Mailman et al. [31] where the system is randomly initialized at a packing fraction far from the jamming point, with subsequent iterative increases of particle sizes and relaxations to remove all overlaps. As angular separation between the centres of neighbouring (touching) ellipses is smaller than for our system parameters ( and ), it is sufficient to calculate only the minimal first eigenvalue to determine the contacts—possible cases with antipodal contacts can be excluded based on ellipse separation alone.

Figure 7

Examples of monodispersed and bidispersed dense packings of spherical ellipses with aspect ratio . In the bidispersed case, half of all ellipses are smaller by a factor of 1.4. (Online version in colour.)

Discussion

Depending on the requirements, the algorithm can be optimized further. For example, with SIMD instructions, evaluations at multiple could be performed with minimal overhead, allowing for faster determination of the correct eigenvalue branch and narrower initial bracket for the optimization algorithm. For the purposes of collision-driven molecular dynamics, the expensive complexity of evaluating pair interactions for a large number of particles can be alleviated by keeping track of nearest neighbours (e.g. by adapting pre-existing methods that make use of the contact functions [32]). Tracking and changing particle positions and orientations while keeping their shapes constant requires keeping track of the rotation matrices in a numerically stable form, which can be done either by tracking the ellipse centre vectors and the vector of its principal component (e.g. through Euler angles) or by using unit quaternions.

Our algorithm is largely based on the algorithm of [23] but has some important differences due to the differences between spherical and Euclidean geometry. On the one hand, spherical geometry of the problem makes it simpler, because in the Euclidean space, translations and rotations have to be considered separately, while on the sphere, the only parameter for the position and orientation of the ellipse is a single rotation matrix. Similarities can be partially restored by handling the Euclidean case in homogeneous (projective) coordinates, but then the confinement surface is a plane, not a sphere, resulting in a different algorithm. Due to this difference, our algorithm requires solving an eigenvalue problem and not a linear system of equations. In general, the eigenvalue problem is numerically more expensive, but for matrices, a closed-form solution is available.

From the aspect of finding the correct solutions, the spherical version of the algorithm is more involved, as the configuration space of possible intersections is topologically non-trivial and splits into different parts based on the behaviour of eigenvalue bands with respect to the parameter . The antipodal doubling means we need additional information to treat different branches of the solution differently. However, as shown in our work, this can be done with a few trivial tests, with the only caveat that the long-range contact function is spliced and undefined (clipped to ) in parts of the configuration space.

Conclusion

The simplicity and speed of the presented algorithm makes it a viable workhorse for future simulations on a sphere, be it interactions of hard particles or general long-range interactions where distances are needed, although the concept of the contact function as a distance metric must be considered with care. Collision detection and generalized distance can be used for Monte Carlo simulations, while molecular dynamics can make use of the intersection vector and the normal vector to the surface as well. Elongation of particles is known to affect optimal packing fraction of random packings in Euclidean space [9,33], and with the presented algorithm, related questions can be answered for packings on a sphere.

Simulations can also be augmented with other potentials that do not use the contact function—for example, multipolar interactions, which may account for elliptical magnetic particles or electrostatically charged macromolecules. The algorithm is viable for particles of different aspect ratios and sizes, so it can be used for simulations of polydisperse particle systems. Another important use case is in representation of arbitrarily shaped objects as isosurfaces of Gaussian sums (called blobs or metaballs in three-dimensional graphics). A product of Gaussians, resulting directly in addition of quadratic forms when constraned to a sphere, also resembles posterior Bayesian update when handling probability models for directional or geographical data, which may be relevant in data processing and machine learning.

Finally, more fundamental questions can also be tackled. Recall that both the Tammes problem and its long-range potential cousin, the Thomson problem, have been well studied not only by physicists but also from the perspective of fundamental and applied mathematics and computer science. Generalization to an anisotropic case is a richer example, which without doubt hides many undiscovered facts about spherical packings.