Group delay of multiple feedback filter
What you're talking about is the time constant, usually denoted \$\tau=RC\$. It serves in calculating the frequency of the filter. The group delay (not the best article on Wikipedia) is calculated based on the whole transfer function, not just \$\omega_n\$, or \$\omega_c\$. It's also not something that has a fixed value. As you can see in the 2nd graph you posted, it's a function of frequency. For a generic allpass, :
$$\begin{align}
H(s)&amp;=K\dfrac{s^2-\dfrac{\omega_p}{Q}s+\omega_p^2}{s^2+\dfrac{\omega_p}{Q}s+\omega_p^2} \tag{1} \\
\Phi(\omega)&amp;=\arctan{\Biggl(\dfrac{\Im{\bigl(H(j\omega)\bigr)}}{\Re{\bigl(H(j\omega)\bigr)}}\Biggr)} \tag{2} \\
\tau_{gd}(\omega)&amp;=-\dfrac{\mathrm{d}}{\mathrm{d}\omega}\Phi(\omega) \\
&amp;=\dfrac{2\omega_p\omega^6+2\omega_p^3\Biggl(1-\dfrac{1}{Q^2}\Biggr)\omega^4+2\omega_p^5\Biggl(1-\dfrac{1}{Q^2}\Biggr)\omega^2+2\omega_p^7}{\omega^8-\dfrac{2\omega_p^2}{Q^2}\omega^6-\omega_p^4\Biggl(2-\dfrac{4}{Q^2}-\dfrac{1}{Q^4}\Biggr)\omega^4-\dfrac{2\omega_p^6}{Q^2}\omega^2+\omega_p^8} \tag{3}
\end{align}$$
Please let me know if I made any mistakes, I tried to prettify the monster. As you can see, an allpass has a double order for the group delay. This is to be expected, since what you have there is 1-2*BP (not because of the 2*BP, but because of the 1-...). A 2nd order transfer function has a 4th order group delay, normally. This can be easily tested:

As you can see, the two traces overlap. The expression for B1 is too long, so here it is, tamed:
(2*omega_p^7-((2*Q^2-2)*(s/sqrt(-1))^2*omega_p^5)/Q^2-((2*Q^2-2)*(s/sqrt(-1))^4*omega_p^3)/Q^2+2*(s/sqrt(-1))^6*omega_p)/(omega_p^8-(2*(s/sqrt(-1))^2*omega_p^6)/Q^2-((2*Q^4-4*Q^2-1)*(s/sqrt(-1))^4*omega_p^4)/Q^4-(2*(s/sqrt(-1))^6*omega_p^2)/Q^2+(s/sqrt(-1))^8)

It's divided by \$j\$ because the expression, itself, is a function of \$\omega\$, not \$j\omega\$, and the gain is 1, for simplicity.
For the case where you might want to think about the group delay at a certain frequency, you need to use the formula above (or whichever you have derived). If you have a flat(-tish) group delay, similar to what you're showing there, your interest might lie in the lower section of the frequency, so you might consider the value at DC -- in the same way a Bessel filter might be considered, since its group delay is flat across the passband. For this case, it's enough to divide the single terms in the numerator, since they are not attached to any \$\omega\$ term and, thus, determine the value at DC. For the case above, \$\omega_p=1\$:
$$\tau_{gd}(0)=\dfrac{2\omega_p^7}{\omega_p^8}=2\;\mathrm{s}$$
which is what can be read in the test, above.

I realize the formula above can look intimidating, especially since the 8th order is actually a 4th order, in disguise. It turned out this way because I expanded the fraction and because of the rather complicated terms I used. So let me try again, with two generic 2nd order transfer functions for an allpass and a bandpass. You'll see much clearer where the doubling comes from, through comparison.
First, instead of writing the transfer functions with \$\omega_p\$ and \$Q\$, let's replace the terms with \$a=\omega_p/Q\$ and \$b=\omega_p^2\$:
$$\begin{align}
A(s)&amp;=\dfrac{s^2-as+b}{s^2+as+b} \tag{4} \\
B(s)&amp;=\dfrac{as}{s^2+as+b} \tag{5}
\end{align}$$
This looks much cleaner. For the group delay, the derivative is applied to the phase which, in turn, is a function of the real and imaginary parts of the transfer function:
$$\begin{align}
A(j\omega)&amp;=\dfrac{(b-\omega^2)^2-a^2\omega^2\;-\;j2a\omega(b-\omega^2)}{(b-\omega^2)^2+a^2\omega^2} \tag{6} \\
B(j\omega)&amp;=\dfrac{a^2\omega^2\;+\;ja\omega(b-\omega^2)}{(b-\omega^2)^2+a^2\omega^2} \tag{7} \\
\end{align}$$
Already there is a hint of what is to come: the imaginary components of \$A(j\omega)\$ have a \$2\$ term in front, compared to \$B(j\omega)\$. The phase would be:
$$\begin{align}
\Phi_A(\omega)&amp;=\arctan{\Biggl(-\dfrac{2a\omega(b-\omega^2)}{(b-\omega^2)^2+a^2\omega^2}\Biggr)} \tag{8} \\
\Phi_A(\omega)&amp;=\arctan{\Biggl(\dfrac{a\omega(b-\omega^2)}{a^2\omega^2}\Biggr)}=\arctan{(\dfrac{b-\omega^2}{a\omega})} \tag{9}
\end{align}$$
The group delay means the negative of the derivative of the phase, and since the argument is the arctangent of a fraction, it follows the generic:
$$\arctan'{\dfrac{f}{g}}=\dfrac{f'g-g'f}{f^2+g^2}$$
I won't practise the math, you can use symbolab, for example, but the results are:
$$\begin{align}
\tau_A(\omega)&amp;=2\dfrac{a(\omega^2+b)}{\omega^4+(a^2-2b)\omega^2+b^2} \tag{10} \\
\tau_B(\omega)&amp;=\dfrac{a(\omega^2+b)}{\omega^4+(a^2-2b)\omega^2+b^2} \tag{11}
\end{align}$$
Clearly, the only difference between the two is a factor of 2.