Why do reverse Bessel polynomials give rise to flat group delay in Bessel filters
It's not the polynomials that make the filter, it's the filter that needed those polynomials in order for it to be a constant group delay.
The theory goes something like this: there is a need for a maximally flat group delay filter. That translates into a very simple Laplace expression:
$$H(s)=\mathrm{e}^{-s\tau}=\frac{1}{\mathrm{e}^{s\tau}}=\frac{1}{\sinh{s\tau}+\cosh{s\tau}}\tag{1}$$
Or the definition for an ideal delay. The starting point is the generic formula for an all-pole Laplace transfer function:
$$H(s)=\prod_{k=0}^N{\frac{1}{a_ks^k}}$$
which is expanded in frequency domain as:
$$H(j\omega)=\frac{1}{\sum_{k=0}^N{(-1)^ka_{2k}\omega^{2k}}+j\omega\sum_{k=0}^N{(-1)^ka_{2k+1}\omega^{2k}}}=\frac{1}{\Re{H(j\omega)}+j\Im{H(j\omega)}}$$
From which the phase and the group delay are:
$$
\angle{H(j\omega)}=\arctan\frac{\Im{H(j\omega)}}{\Re{H(j\omega)}} \\
\delta(\omega)=-\frac{\mathrm{d}\angle{H(j\omega)}}{\mathrm{d}\omega} $$
To solve this, consider the fraction from the two equations above as \$\frac{f(x)}{g(x)}\$, and the derivative becomes:
$$
\frac{\mathrm{d}}{\mathrm{d}x}\frac{f(x)}{g(x)}=\frac{1}{f(x)^2+g(x)^2}\left(\frac{\mathrm{d}g(x)}{\mathrm{d}x}f(x)-\frac{\mathrm{d}f(x)}{\mathrm{d}x}g(x)\right)
$$
To avoid continuing with generic expansions and sums, or with Taylor series, let's consider the case for a 3rd order filter, which has 4 terms in the denominator, thus two real and two imaginary. For a group delay to be flat (and unity), the derivative of the phase at \$\omega=0\$ must be zero, thus the phase must be linear, or \$\omega\$ itself. Since the delay will be the derivative of the phase, the first term of the \$\omega\$ term for the phase will be \$1\$ for the derivative. So, from the (non-monic) transfer function we end up with:
$$
H_3(s)=\frac{1}{a_3s^3+a_2s^2+a_1s+1} \\
H_3(s\tau)=\frac{1}{a_3(s\tau)^3+a_2(s\tau)^2+s\tau+1}
$$
The group delay, considering the derivations from above, becomes:
$$
\delta(\omega)=-\frac{\mathrm{d}}{\mathrm{d}\omega}\arctan{\left(\frac{a_1\omega-a_3\omega^3}{1-a_2\omega^2}\right)}=\frac{a_2a_3\omega^4+(a_2-3a_3)\omega^2+1}{a_3^2\omega^6+(a_2^2-2a_3)\omega^4+(1-2a_2)\omega^2+1} 
$$
Flat group delay implies the coefficients of the same powers of \$\omega\$ from both numerator and denominator need to be made equal. Then the following system of equations can be formed:
$$
\left\{
\begin{aligned}
a_2-3a_3&amp;=1-2a_2 \\
a_2a_3&amp;=a_2^2-2a_3
\end{aligned}
\right. \\
\Rightarrow \\
a_3=\frac{1}{15}\,,\quad a_2=\frac25 \\
\Rightarrow \\
H_3(s\tau)=\frac{1}{\frac{1}{15}(s\tau)^3+\frac25(s\tau)^2+s\tau+1}=\frac{15}{(s\tau)^3+6(s\tau)^2+15s\tau+15}
$$
The denominator is recognized now as a Bessel polynomial. The same derivation can be applied for any other order, and the transfer function will have a Bessel polynomial in the denominator, and a DC scaling term in the numerator.
This is why I said in the beginning that the Bessel polynomials are not chosen because they have some special property that gives a flat group delay, it's because the derivation of a filter that needs to have a flat group delay leads to a transfer function having the Bessel polynomials.