Why do bees create hexagonal cells ? (Mathematical reasons)

There are two principles at play here: a mathematical principle that favors hexagonal networks, and a physical principle that favors a network with straight walls.

The mathematical principle that prefers hexagonal planar networks is Euler's theorem applied to the two-torus $\mathbb{T}^2$ (to avoid boundary effects), $$V-E+F=0,$$ with $V$ the number of vertices, $E$ the number of edges, and $F$ the number of cells. Because every vertex is incident with three edges$^\ast$ and every edge is incident with two vertices, we have $2E = 3V$, hence $E/F=3$. Since every edge is adjacent to two cells, the average number of sides per cell is 6 --- hence a uniform network must be hexagonal.

_{$^\ast$ A vertex with a higher coordination number than 3 is mechanically unstable, it will split as indicated in this diagram to lower the surface energy.}

_{blue: total edge length for the left diagram (diagonals of a unit square), gold: total edge length for the right diagram, as a function of the length $x$ of the splitting.}

Euler's theorem still allows for curved rather than straight walls of the cells. The physical principle that prefers straight walls is the minimization of surface area.

_{source: Honeybee combs: how the circular cells transform into rounded hexagons}

An experiment that appears to be directly relevant for honeybee combs is the transformation of a closed-packed bundle of circular plastic straws into a hexagonal pattern on heating by conduction until the melting point of the plastic. Likewise, the honeybee combs start out as such a closed-packed bundle of circular cells (panel a). The wax walls of the cells are heated to the melting point by the bees and then become straight to reduce the surface energy (panel b).

There is this theorem of Thomas Hales from 1999, which proves the Honeycomb Conjecture:

Theorem. Let $\Gamma$ be a locally finite graph in $\mathbb{R}^2$, consisting of smooth curves, and such that $\mathbb{R}^2\setminus \Gamma$ has infinitely many bounded connected components, all of unit area. Let $C$ be the union of these bounded components. Then $$ \limsup_{r \to \infty} \frac{ \mathrm{perim}\, (C \cap B(0, r))}{\mathrm{area}\, (C \cap B(0, r))} \geq \sqrt[4]{12} $$ Equality is attained for the regular hexagonal tiling.

So basically it is optimal way to partition the plane into cells of equal area using least amount of perimeter. This doesn't account for the fact that honeycomb lattice is 3d and not exactly cylindrical with hexagonal cross section.

The paper introduction has a bit of discussion https://arxiv.org/abs/math/9906042

Isn't it just the 2d sphere packing? If one assumes that the larvae needs a disc of fixed radius to grow up to an adult form and that the bees want to have as many cells as possible then the hexagonal lattice is the optimal one.