Symmetry break in the eight bubble compaction

1.   Introduction

Often used for children's enjoyment, soap bubbles are the simplest physical example of a lot of mathematical problems. They are solutions of the minimal surface problem ^[1]: the connection between soap films and minimal surfaces was established by Gauss in 1830. Indeed, working on capillarity problems, he found that at equilibrium any liquid surface minimizes the potential energy caused by molecular forces. In the specific case of soap films, such a energy is proportional to the area. Hence, soap films are nothing but the physical representation of stable minimal surfaces. Then, soap bubbles are used to study two dimensional hydrodynamics and turbulence ^[2]: their use is motivated by the fact that atmospheric turbulence at large scales displays two dimensional features due to the small thickness of the Earth's atmosphere ^[3]. Indeed, to mimic the emergence of hurricanes or cyclones in the atmosphere ^[4], which are the manifestation of turbulence, a half bubble is heated at the equator giving rise to thermal convection which leads, under some physical conditions, the appearance of the searched long lived isolated vortices ^[5]. Finally, they solve an optimization problem ^[6]: when two or more bubbles cluster together, their configuration obeys some physical optimization rules. Indeed, assembling several bubbles traps pockets of gas in a liquid and results in foam: the surfactants added to the liquid stabilize the bubbles by reducing the surface tension and by arranging themselves at the liquid/gas interfaces ^[7].

Regarding the behavior of a single soap bubble, everything is known. What is not completely understood is the geometrical and mechanical properties of a cluster of many bubbles, known as a foam. For instance, its optimal rearrangement in space is still matter of debate.

In $2$D there are more results: Hales proved the honeycomb conjecture, which states that the partition of the plane into regular hexagons of equal area has least perimeter, i.e., it minimizes the perimeter fixing the area ^[8]. In this context, it has been proved that the optimal (i.e., minimal) configuration exists for $N$ clusters ^[9,10] and Cox et al., obtained their numerical visualizations up to $N = 200$ ^[11]. Due to the non linearity of the problem, in a lot of physical situations, the equilibrium solution is only stable with respect to small displacements, i.e., it is not a global minimum of the system. This aspect leads to mechanical instabilities which break the symmetry of the system ^[12,13,14]. These instabilities should be classified into two groups ^[15]: the first one is primarily due to the presence of an asymmetry in the system and it can be studied though a purely static model. For the second one, the system jumps away from this initial configuration to settle in another which is not close to it and its description may require a dynamic model. Just to mention some famous instabilities, we recall the $T1$ process in $2$-dimension ^[16], and the $2D$ flower cluster ^[13].

If we move to $3$D, the problem is minimizing the area functional and few exact result exist. The only rigorous one is the proof of the Double Bubble Conjecture ^[17], which states that the standard double bubble provides the least-area enclose and separates two regions of prescribed volume in $\mathbb{R} ^3$. As regarding the numerical results, Kelvin ^[18] proposed an optimal candidate structure with identical cells, which has been numerically refuted in ^[19]: by numerical calculations, an agglomerate of two different types of bubbles has less perimeter, for fixed area.

The aim of this work is to study the symmetry break of a eight bubble compaction in $3$D, i.e., we want to investigate the mechanical cues driving the three-dimensional packing of $8$ bubbles displayed in a recalling foam structure. To circumvent the difficulty of a variational approach, i.e., the minimization of an area functional satisfying some geometrical constraints, we follow a different strategy. First of all, in Section 2, we fix the geometrical arrangement of the eight bubbles as the solution of the Tammes' problem ^[20]: we exploit a geometrical principle of maximal mutual distance between neighbor points on a spherical surface obtaining seven symmetrical peripheral spheres tangent to the central one ^[21]. Then, compaction is produced by packing the outer bubbles along the radial direction of the aggregate. The obtained agglomerate recalls the foam structure ^[22]: the central sphere is completely covered, while the peripheral ones have a free-curved surface. While our construction does not ensure that the obtained final configuration is the minimal one, we will prove that the our tessellation on the central sphere is the optimal one among all the possible ^[23] according to physical assumptions: the liquid/liquid interface is favored versus the liquid/gas one and it maximizes the volume ^[7]. Fixing the geometrical arrangement, in Section 3, we look for the balance of forces that produces such a configuration ^[6]. Enforcing the mechanical equilibrium and the conservation of volume, we derive the tensional balance laws on this geometrically optimal packing. Finally, in Section 4, we discuss the results and we add few concluding remarks.

2.   Geometrical principle

2.1.   Spatial arrangement

In this section, we introduce a geometrical principle that we exploit to describe the spatial arrangement of the $8$-bubbles configuration, to determine the position of the seven bubbles surrounding the central one. The coordinates of the peripheral bubbles centroids are given as solution of the classical Tammes's Problem ^[20]: determine the arrangement of $n$ points on the surface of a sphere maximizing the minimum distance between nearest points (maxmin principle). This is equivalent to determine (up to rigid rotations) the $n$ unit vectors $\{{\bf{r}}_i\}$ such that

is maximum, where the limit $m \rightarrow + \infty$ selects the distance among closest points only.

In our case $n = 7$; we want to find the position of seven points on a sphere with minimum distance from their nearest neighbors.

Here, we exploit the graph theory ^[24] to find this maximizing configuration. A set of $n$ points on a sphere forms a graph $G$ of $n$ points connected by arcs of great circles of length $a$ ^[25]. The maximal spatial arrangement of seven points can be obtained by the projective argument presented in ^[25].

Consider a frame of reference centered in $O = (0, 0, 0)$ and the coordinates on the spherical surface $(r, \theta, \phi)$, where $\theta \in [0, 2\pi)$ (longitude) and $\phi \in [0, \pi]$ (latitude) on $S^2$, namely

where $r$ is the radius of the sphere $N$ and $S$ the North and South Pole, respectively. Three points $\left\{A, B, C\right\}$ are placed at the same latitude on the surface of the central sphere, such that they are connected by arcs of length $a$ and they form an equilateral triangle centered in $S$, see the triangle $\widehat{ABC}$ with in Figure 1a. Three more identical triangles are then created, adjacent to the former ones, with vertices $D, E$ and $F$: they share the same latitude too, see the triangle $\widehat{ABF}$ in Figure 1b. The final step is then to connect $D, E$ and $F$ with $N$ and vary the radius $r$ (for fixed $a$) until also the latter arcs have length $a$^*, (see Figures 1c and 2). The associated extremal graph defines four triangles and three quadrangles on the spherical surface, as illustrated in Figure 3.

^* In a fully equivalent way, one can fix the radius $r$ and vary the chord length $a$.

Fundamental relations of spherical trigonometry tell us that the internal angle of an equilateral spherical triangle is $\alpha = \frac{4\pi}{9}$, while the arc angle $\beta$ with respect to the centre of the sphere $\beta$ is given by ^[26]

as illustrated in Figure 4.

This allows to find explicitly the linear relation between arc length and radius, i.e., $a = \beta r$. The length $\ell$ of the chord between closest points is

2.2.   Coordinates

By construction, the spherical distance on $S^2$ of the points $D, E$ and $F$ from $N$ is equal to $a$, so their latitude is the angle $\phi = \beta$. For the triplet of points $\{A, B, C\}$, the calculations are a little bit more elaborated. Let $r_{\phi}$ be the radius of the circumference defined by the intersection of the sphere and the plane where $A, B$ and $C$ are. The following relation holds

where ${{\phi }_{abc}}$ is the latitude of the points $A$, $B$ and $C$.

Since the chord length is the same for the spherical arc and for the in–plane circle, we find

By combining Eqs. (2.3) and (2.4) we get

Summarizing, the coordinates of the seven Tammes points depicted in Fig. 2 are

This configuration, given by the maxmin principle (2.1), is here adopted as the ideal reference bubble arrangement: seven spheres are tangent to the former one in the Tammes' points, as illustrated in Figure 5a.

2.3.   "Bubble compaction": tiling the central sphere

The tessellation of the spherical surface illustrated in the previous sections, is composed by four equilateral triangles and three quadrilaterals, see Figure 3a. However it can be made of triangles only by connecting points $\{D, E, F\}$, see Figure 3b. The triangle $\{D, E, N\}$ is not equilateral, since the distance between $D$ and $E$ is equal to $1.34 \, a$. On the basis of such a triangular tessellation we can produce a dual tessellation connecting the circumcenters of the triangles: the locus where the axis of the edges cross each other (Figure 6).

The Tammes' points are the centroids of the polygons that define the dual tessellation (see Figure 6 and for more mathematical details Appendix A). The bubble packing is obtained ideally moving each peripheral bubble, initially tangent to the central one in the Tammes' points, towards the origin $O$ along the radial direction, as illustrated in Figure 5b, while enforcing the volume conservation. In other words, to pack the bubbles aggregate we generate a collection of flat surfaces of contact among bubbles starting from the maxmin distribution of the tangent points: each peripheral bubble adheres to the central one moving centripetally, see Figure 5b. At the same time we shuffle the peripheral and the central spheres to preserve the initial volume $\mathscr{V}$. The contact surfaces between central and peripheral bubbles obtained by such a dive and shuffle procedure are nothing but the polygons obtained connecting the points of the dual tessellation defined above.

The final configuration is a tiling of the spherical surface of the central bubble with seven polygons, as illustrated in Figure 7:

● an equilateral triangle centered in $N$ (area $\simeq 0.66 \, a^2$),

● three quadrilaterals centered in $A, B, C$ (area $\simeq 0.7 \, a^2$),

● three pentagons centered in $D, E, F$ (area $\simeq 0.74\, a^2$).

2.4.   Geometrical optimality

While the produced polygonal surface covers the central sphere, it naturally arises the question if such a tiling is optimal according to some suitable criterion. The problem to cover a spherical surface with polygons is old, rigorous results dating back to Cauchy ^[27]. In the present context, all the bubbles are identical and it is therefore tempting the idea to cover the central bubble with identical polygons. Unfortunately, this is not allowed by Euler's Polyhedron Formula^†.

^†Euler's Polyhedron Formula has been proved by Cauchy ^[27] and it gives a relation among the number of faces, edges and vertices of a polyhedron, such as

The number of vertices and edges are related with the number of faces faces $F$ as follows

where $n$ is the number of edges of the polygon at hand, while $m$ is the number of faces which insist on the same vertex. We are interested in the case $F = 7$. By elementary calculations one can easily see that there is no $n$ for which a suitable integer $m$ exists.

There is no regular heptahedron. In order to prove if our tiling is optimal, we can construct all the convex polyhedra with $7$ faces which can be inscribed into a sphere of fixed radius $r$. By the software Plantri ^[23], we find $34$ convex polytopes with seven faces. They can be classified in terms of number of edges and number of vertices, see Table 1.

In general, an optimal packing rearrangement of soap bubbles is the one which saves energy, i.e., we expect the chosen one is the equilibrium configuration which should minimize a suitable assigned energy functional. However, due to the huge number of variables involved, we are not able to write down such a functional, hence we use some geometrical and physical considerations to deduce the optimality and maybe the minimality of the obtained configuration.

From a physical point of view, it is known that creating an interface liquid/liquid energetically costs less than one liquid/air ^[7]. Since the number of edges represents the intersection between nearest bubbles placed on the Tammes' points which creates a liquid/liquid interface, maximizing this number implies that we are requiring that the final polytope on the central sphere has more liquid/liquid intersections then liquid/air ones, i.e., we are looking for a polyhedron $\widetilde{P}$ such that

where $E_{\mathcal{P}_{i}}$ is the number of edges the $i$-polyhedron. In this way, we can reduce the number of polyhedra: we pass from $34$ convex polytopes with $7$ faces to just $5$ in which our tiling is included. In Appendix B, we show these $5$ polyhedra obtained and drawn by the software Plantri.

We need another condition to select just one configuration. Using the fact that the shape of a single soap bubble is spherical, which easily derives from the isoperimetric inequality, i.e., the sphere is the only solid that, fixing the area, it maximizes the volume, we search among the five configurations obtained by Plantri, the one that has the maximum volume, indeed we look for a polyhedron $\overline{\mathcal{P}}$ such that

where $\mathcal{L}^3$ is the volume measure. By numerically computing the five volumes, we find that the tiling obtained as the dual of the Tammes' one is the one with maximum volume, i.e., the optimal one according to our criteria.

Remark 1. Since each polyhedron is not regular, we do not have an explicit formula to compute the volume. However, each polyhedron can be divided into $7$ pyramids, where the basis is a face. Using this geometrical argument, the total volume can be computed as the sum of the volumes of the pyramids.

2.5.   Surfaces, edges and vertices

The dual tessellation defines ten nodes and it belongs to $\mathcal{C}_{3z}(1, 3, 3, 3)$, see Figure 7b^‡. We denote the nodes of the dual tessellation on the basis of the vertices of the Tammes' triangles they belong to, such as

^‡ $\mathcal{C}_{nz}$ is the group of a cyclic symmetry after a rotation $2\pi / n$ with respect to the axis $z$ ^[28,29]: the configuration is invariant for rotations of an angle $2 \pi/3$ around the $z$ axis. The notation $(1, 3, 3, 3)$ denotes how many nodes of the tessellation share the same longitude.

Because of the symmetry of the problem, there are only four independent nodes, as depicted in Figure 7b.

The dual tessellation is formed by 15 edges, only 4 of them being independent. Each edge is identified by the polygons it belongs to on the tiled surface and the boldface denotes the unit vector parallel to the edge. Therefore, ${\bf{c}}_1 = \overline{V_1V_2}$ denotes the edge between two quadrilaterals, ${\bf{c}}_2 = \overline{V_2V_3}$ separates a quadrilateral and a pentagon, ${\bf{c}}_3 = \overline{V_3V_4}$ separates a pentagon and a pentagon and ${\bf{c}}_4 = \overline{V_2V_3}$ is between a pentagon and the triangle.

We generate a three dimensional structure projecting radially the dual tessellation, by an height to be fixed later on the basis of volume conservation arguments. Each vertex of the tessellation on the central bubble has therefore a corresponding outer one that we denote by $V_i^h$, $i = 1, 2, 3, 4$. The connection between inner, outer and side surfaces is defined by two classes of edges:

● ${\bf{c}}_5 = \overline{V_1V_1^h}$, ${\bf{c}}_6 = \overline{V_2V_2^h}$, ${\bf{c}}_7 = \overline{V_3V_3^h}$, ${\bf{c}}_8 = \overline{V_4V_4^h}$ point radially,

● ${\bf{c}}_9 = \overline{V_1^hV_2^h}$, ${\bf{c}}_{10} = \overline{V_2^hV_3^h}$, ${\bf{c}}_{11} = \overline{V_3^hV_4^h}$, ${\bf{c}}_{12} = \overline{V_2^hV_3^h}$, are parallel to the ones on the tessellation of the central bubble.

At this stage the geometrical characterization of the 8-bubbles configuration derived on the basis of a maximum-minimum distance of the centroids of the peripheral bubbles is completed. The inner bubble has no free surface: it is surrounded by contact interfaces with other bubbles only. The external ones have the shape of a pyramidal frustum covered by a laterally cut spherical cap: the lower basis is the polygon generated by the adhesion with the central bubble, lateral sides are flat too, their edges being radially oriented, the upper basis of the frustum is a radial projection of the lower one. The upper geometrical structure is a spherical vault on a polygonal frustum, intriguingly known since the Middle Age in Sicilian architecture ^[30]. The radius of the spherical cap and the height of the frustum are to specified on the basis of balance and conservation arguments discussed below.

3.   Mechanical balance

In this section, we compute the surface tensions that make the geometrical packing mechanically equilibrated. We remark that the central sphere has only bubble–bubble contact interfaces, while the peripheral ones also possess a traction–free surface. Each bubble–bubble interface and each free surface is characterized by a tension $\tau_i$, defined as the energy density per unit area of the liquid/liquid or liquid/air interfaces ^[31]. Thus, for the strong symmetry of the problem, we have ten unknown independent tensions $\tau_i$, one variable for each portion of surface. Precisely, they are: three on the central bubble, four on lateral bubble-bubble interfaces and three at the free surface denoted by

where the subscript identifies the surface of the polygon it applies to and the superscript $s$ specifies the tensions at the free surfaces.

3.1.   Tensional balance

First, we enforce the mechanical equilibrium imposing that the surface tensions are balanced on each independent edge ${\bf{c}}_i$, where $i = 1, ..., 12$ (see Figure 8a). Three (flat or curved) surfaces are attached to each edge, their local orientation being denoted by the normal unit vectors ${\bf{n}}^i_j$, $j = 1, 2, 3$. The balance of tensions on each edge is defined by the sum of the tensions, oriented orthogonally to the edge and in-plane with the corresponding interface (see Figure 8a). Therefore it must hold

Equation (3.2) defines 36 scalar equations, 12 of them being trivially null because all the summed vectors are in the plane orthogonal to the edge under consideration. With the help of a symbolic software^§, we eventually find that, given the unit vectors ${\bf{n}}^i_j$ only 10 of them are independent. The equations are detailed in Appendix D in Tables 2–4.

^§We used Mathematica (Wolfram Inc., Version 12).

The orientation of the normal unit vectors is not defined according any specific rule because it is expected to affect only the sign of tension, that we know to be positive.

The linear system Eq (3.2) is however not closed because, while the direction orthogonal to the flat surfaces is uniquely defined, the edge contribution of the tension defined on the free surface depends on the curvature of the surface itself. Curvature, tension and pressure gap on the free surface of the peripheral bubbles obey the Young-Laplace equation ^[32] in the following way

where $R_i$ is the radius of curvature of the free surface of the $i$-th bubble and $\Delta p$ is the difference between the outer and the inner of pressure and there is an extra factor $2$ since the surface of the bubble is composed by two leaflets. Since we have only three independent types of polygons on the tessellation, Eq (3.3) gives three independent equations.

3.2.   Volume conservation

Finally, we have to impose the conservation of bubble volume under compaction. While the central (packed) bubble is bounded by flat interfaces, the peripheral ones have the shape of a pyramidal frustum covered with a spherical vault (see Figure 8b). The basis of the pyramidal frustum are

● the interface with the central bubble;

● its radially directed homothetic projection, by a factor $\frac{r+h}{r}$, where $h$ is the radial height of the intersection surface among three adjacent cells (see the yellow segment in Figure 9).

The value of $h$ has to be fixed on the basis of volume conservation arguments: the sum of the pyramidal and apsal volumes must be equal to the common volume of all bubbles.

Details about the calculation of these volumes are given in Appendix C. As the area of each polygonal basis and the curvature radius of the apse are different, the radial height $h$ of the cells is not actually the same; however differences are below $1\%$.

3.3.   Results

We can solve the system of $16$ equations given by Eqs (3.2) and (3.3), constrained to volume conservation, with respect the $16$ unknowns: 10 tensions, 3 curvature radii and 3 heights. Using some experimental data coming from the foam literature, ^[22], we can fix the pressure difference $\Delta p = 50 \, {\rm Pa}$ and the radius of the round bubble as $r = 1 \, {\rm mm}$, before the compaction process. Numerical solution of the nonlinear system of equations predicts the following surface tensions

The computed radii of curvature are

The obtained surface tensions are consistent with experimental results ^[22]: to create a soap bubble, the surface tension has to be less than the one of water, which is $\tau_{\rm water} \simeq 73 \, m{\rm N}/{\rm m}$, otherwise the bubble cannot exist. The obtained field of forces Eq (3.4) is the one at the equilibrium. We immediately notice that there is an anisotropy in the distribution of the field of forces: surface tensions of bubble-bubble interfaces with normal oriented in the circumferential direction of bubbles aggregate (second line of Eq (3.4)) are larger than the ones with normal unit vector pointing radially out of the aggregate (first and third line of Eq (3.4)).

This result supports our conjecture: the anisotropy in the mechanical cues may be the cause of the symmetry break, i.e., there might be a preferential direction of the next topological instability ^[13]. Indeed, from experiments and numerical results, it is known that a similar aggregate, due to some physical involved parameters, can develop an asymmetry or a topological transition ^[15]. The study of the stability of this configuration is out of this paper. This result wants just to show that the distribution of forces in the equilibrium configuration itself is not symmetric, hence we can state that any small perturbations can change the rearrangement of forces inside the system and can develop a topological transition which breaks the starting symmetrical structure.

4.   Final remarks

In this work, we studied the symmetry break of a particular configuration of $8$ spheres, showing the anisotropy in the distribution of the field of forces of the equilibrium position, that might possible originate a topological transition and break the symmetric structure of the starting agglomerate ^[13,15]. This result can be applied to a physical situation, i.e., the study of foam, since the selected rearrangement remembers the one of soap bubbles in a single module of the foam structure ^[22].

We considered $7$ identical spheres symmetrically surrounding a central one: their initial position is dictated by the solution of the Tammes' problem ^[20]. Neglecting any dynamical process, the final configuration is obtained by a compaction process which results into a full tiling of the central sphere. By introducing physical criteria of optimality dictated by the energy minimality, Eq (2.7), and by the volume maximality, Eq (2.8), we deduced that our polyhedra is the optimal one among all the $34$ convex polytopes inscribed into a sphere with radius $r$ ^[23], since due to Euler Polyhedra Formula no regular heptahedrons exist ^[27].

Fixing this geometrical arrangement, we looked for the force balance that realizes such a configuration: we computed balance of forces on every edge, we forced the conservation of volume (by experimental evidences soap bubbles can be assumed to be incompressible ^[33]) and we imposed the Laplace law on the possibly curved free surface. We obtained a force field, Eq (3.4), which fulfills an acceptable physical range ^[22], but it shows an anisotropy in its orientation, see second line in Eq (3.4). This result suggests that a difference in tension, generated by a purely mechanical principle, might be crucial for next topological transitions and for the development of anisotropies ^[15]. In this respect, we conjecture that the expulsion of a single bubble, dictated by any small perturbations, in the flower cluster in $2$D ^[13] can be replicated in higher dimensions, breaking the symmetrical structure. Indeed, our configuration represents a $3$D asymmetrical flower cluster. Hence, up to construct a suitable dynamical problem to study the evolution in time of such a configuration, a possible solution might be the ejection of one bubble from the agglomerate.

Future efforts will be to devoted first to try to understand more rigorously if the obtained dual tessellation realizes the minimum of a suitable energy functional, then to reproduce this system both in a laboratory and numerically to study the dynamical evolution of this agglomerate.

The main novelty of this work is the introduction of a new non-destructive technique to study the geometrical rearrangement and the change of shape of a configuration by using mechanical and geometrical considerations. Indeed, the used mathematical method here is applied to a specific context, i.e., the symmetry break of a $8$-bubble compaction, but it can be extended to better understand many different physical phenomena. For instance, it can be used to design new meta-materials, where it is fundamental to know a priori the balance of forces, or to study the mitosis of cells. Indeed, just by knowing their geometrical rearrangement at a fixed stage, we can determine if the distribution of forces has an anisotropy which can favour the duplication process along a particular direction ^[34]. This means that our method can give an insight for a more detail comprehension on the mechanics of morphogenesis of a variety of tissues.

Acknowledgments

The author thanks Davide Ambrosi, Pasquale Ciarletta, Luca Lussardi, Alfredo Marzocchi and Maurizio Paolini for helpful suggestions and fruitful discussions. This work has been partially supported by INdAM$-$GNFM.

Conflict of interest

The author declares no conflict of interest.

Appendix A.   Dual tessellation of the central bubble

A visual three-dimensional representation of the dual tessellation produced on the central sphere by the bubble packing is depicted in Figures 10 and 11. The polygons representing the flat bubble–bubble interfaces are here plotted inside the original spherical bubble of radius $r$ before reshuffling the polyhedron to recover the original bubble volume.

As we can see from Figures 10 and 11, the projection of Tammes' points along the radial direction represents the centroid of each polygon on the dual tessellation. Viceversa, the vertices $V_i$ with $i = 1, \dots, 10$ are the centroids of the triangles of the modified Tammes' tessellation, see Figure 3b.

Appendix B.   Polyedra with seven faces

In this appendix, we want to show the five polyhedra among all the $34$ convex polytopes with seven faces inscribed into a sphere of a fixed radius $r$, which satisfies the first optimal criterium, i.e., Eq (2.7). By the software Plantri, we can draw them and their graphical representation is presented in Figure 12.

Appendix C.   Calculation of volumes

In order to compute the volume of the peripheral bubbles, we have to mathematical describe the shape of a pyramidal frustum covered by a vault. The vault is obtained radially cutting slices of a spherical cap drawn on the vertices of the external basis of the frustum, by prolongation of the lateral surfaces of the frustum itself (see Figures 8b and 14a). From a mathematical point of view, this calculation is a bit technical since the final solid is not a known or a common one.

First of all, we have to better clarify which are the involved unknowns. As it concerns the packing parameter, which is used to define the position of the vertices on the free surface and it is obtained though the interaction among three surfaces, we a priori have four different values of packing parameter, i.e.,

where the letters into brackets denote the interaction between three polygons, i.e., $P$ for pentagons, $Q$ for quadrilaterals and $T$ for the triangle.

We denote with the symbol $\, \tilde{}\,$ the centroids of the dual tessellation polygon: $\tilde A, \tilde B, \tilde C, \dots$, and with $C_i$, with $i = A, B, C, D, E, F, N$, the center of the sphere passing through the vertices of the homotetically projected polygons. Since we have three type of polygons, we make the calculations only on a representative of each class, such as on the quadrilateral $A$, on the pentagon $D$ and on the triangle $N$. We define $h_P = \tilde D - C_D$, $h_Q = \tilde A- C_A$ and $h_T = \tilde N - C_N$ (see Figure 13) the distance between the centre of the sphere associated with the vault curvature and the centroid of the polygonal inner basis of the frustum (the projection of the Tammes' points on the central bubble interfaces). These distances have to be fixed in order to enforce conservation of volumes of the peripheral bubbles. Depending on the type of adjacent polygons (pentagon-pentagon, pentagon-triangle, ....), a different packing parameter $h_i$ with $i = 1, 2, 3, 4$ is expected (see Eq (C.1)). For illustrative purposes, here below we only show how the radial height $h_2$, obtained by the intersection among the sphere constructed on the triangle and on the two adjacent pentagons. We are going to show that it can be rewritten as a function of $h_P$ (Figure 13).

The sphere of radius $R_P$, centered in $C_{D}$ is defined by the equation

where the coordinates of $C_{D}$ are

Analogously, we construct the sphere on the adjacent pentagon $E$, so that

The same compation parameter $d$ scales the coordinates $\left(x_{C_{E}}, y_{C_{E}}, z_{C_{E}}\right)$ because both $D$ and $E$ are pentagons.

The intersection between the two spheres is a circumference, and the edge ${\bf{c}}_3 = \overline{V_3V_4}$ lies on it. In the same way, we can also consider another circumference obtained by the intersection of, for instance, the sphere centered in $C_{D}$ and the adjacent constructed on the triangle, i.e., centered in $C_{N}$. The length of the side of the frustum $h_2$ is the distance between the node $V_4$ of the dual tessellation and the outer intersection point of the two circumferences defined above. The parametric representation of the radial line passing through $V_4 = (x_{V_{4}}, y_{V_{4}}, z_{V_{4}})$ is given by

where $t$ is a positive real parameter. The intersection between the two spherical surfaces Eqs (C.2) and (C.3) with the line Eq (C.4) is the point ${\bf x}_0$ on the free surface with the following coordinates (where $z_{V_{4}}\neq 0$ by construction)

Hence, we define

This procedure can be repeated on all the lateral interfaces in order to calculate all the $h_i$ with $i = 1, 2, 3, 4$.

The second and final step is to write down the volume of the solid (pyramidal frustum plus spherical vault) as a function of $h_P$, $h_T$ and $h_Q$ calculated as in Eq (C.5). The volume of the pyramidal frustum is

where $A_{\rm bottom}^{i}$ is the area of the $i$-th polygon on the dual tessellation, $A_{\rm top}^{i} = \frac{r+h}{r} A_{\rm bottom}^{i}$ is the area of the upper basis, and $H_{\rm py}^{i}$ is the height of the pyramidal frustum.

The volume of the spherical vault is nothing but the volume of the laterally cut spherical cap. For the sake of simplicity, we consider here the spherical vault based on the triangle $N$. We use local coordinates with origin in $C_N$, the centroid of the triangle is in $\tilde{N} = (0, 0, h_T)$.

In terms of the local coordinates $(x', y', z')$, the volume of the spherical cap is

and its measure is obtained by standard volume integration. We have to subtract to the measure of $\Omega'$ the volume of the slices obtained by prolongation of the flat interface defined by the vertices $V_3V_7$ (and so on). It is worth to remark that the plane attached to $V_3V_7$ is not perpendicular to basis of the spherical cap. So, first of all, the expression of the area of a circular segment of radius $\rho$ and chord $b$ is

Here $\rho$ and $b$ are functions of the quote $z' \in [h+h_T, h_{\rm max}]$ and $h_{\rm max}$ has to be determined, see Figure 14a. The value of $h_{\rm max}$ is fixed at the $z'$-level such that the homothetic projection of the inner triangle is circumscribed into the circumference. Namely, a plane at given $z'$ crosses the plane defined by $V_3V_7$ along a line, that we represent by its equation $ax'+by'+c = 0$. To obtain the upper integration bounds in $z'$ we impose

where $\bar x' = 0$, $\bar y' = 0$ and the left-hand-side of Eq (C.9) is exactly the radius $\rho$ of the circumference at fixed $z'$.

By solving (C.9), we get

where $a = 2.48-2.96h_{T}$, $b = 0$, $c = 1$ and $d = 0.42-h_{T}$ and $h_{\rm max}$ is the positive value, since it belongs to $\Omega'$.

Finally, the volume of the bubble constructed on the triangle is given by

where $|\cdot|$ denotes the measure of the volume of $\Omega'$. All the calculations are computed numerically by using the Newton's method with the software Mathematica 12 (Wolfram Research, Champaign, IL, USA). For the other polygons, the computation is similar up to a rotation which has to be done before computing the translation of the centre of the frame of reference, see Figure 14b.

Appendix D.   Equations of balance of tensions

The detailed expressions of Eq (3.2) and their geometrical representation are listed in Tables 2–4. We introduce different superscripts to distinguish the different directions. The unit vectors denoting the direction of the force on the side edge is denoted by the symbol $\mathit{\boldsymbol{t}}$, $\mathit{\boldsymbol{t}}^{r}$ for the ones oriented in the radial direction, $\mathit{\boldsymbol{t}}^s$ on the free surface.

In the following tables, the unit vector normal to the free surface applied in the homothetical vertex of the dual Tammes' tessellation are not reported due to absence of space. We collect them below, i.e.,

where $h_P$, $h_Q$ and $h_T$ are the distance between the center of the sphere associated with the vault curvature and the centroid of the polygonal inner basis of the frustum, for more details see Appendix C.