Derivation of Bessel filters
The optimality criterion for Bessel (low pass) filters is that they have a maximally flat group delay at $\omega=0$, as correctly pointed out in a comment by Peter K.. The derivation is slightly tortuous, but you can find it for the continuous-time case in these lecture notes, which are taken from the book Introduction to Circuit Synthesis and Design by Temes and La Patra. I can recommend this book to anybody who wants to gain a deep understanding of electrical network theory.

The interesting fact is that unlike for other standard filters (Butterworth, Chebyshev, Cauer), the discrete-time version of a Bessel filter cannot be derived from its continuous-time counterpart via the bilinear transform. The reason is the frequency warping inherent in the bilinear transform, which doesn't affect the optimality conditions of the other standard filters (maximal flatness of the magnitude, or equiripple behavior of the magnitude, or a combination thereof), but it does not preserve the maximally flat group delay. A discrete-time filter derived from an analog Bessel filter via the bilinear transfer will not be optimal in any way, and it will not have a maximally flat group delay. This means that a discrete-time Bessel filter must be derived directly in the discrete-time domain. This has been done in this famous paper by Thiran. Again, the derivation is slightly painful but I'm afraid there is no easier route.

What you can take away:


Bessel filters have a maximally flat group delay.
Discrete-time Bessel filters cannot be obtained from continuous-time Bessel filters via the bilinear (or any other) transform.