Solution (source code)

= Solution

For a <plane wave>, $\partial_t$ acts as $-i\omega$ and $\partial_x$ as $ik$. Substitution gives
$$
\omega^2(\omega-\omega_n)-c^2k^2(\omega-\omega_d)=0,
$$
so, for $\omega\ne\omega_d$, the <dispersion relation> is
$$
\boxed{k^2=\frac{\omega^2(\omega-\omega_n)}{c^2(\omega-\omega_d)}.}
$$
Assume $\omega_n>\omega_d>0$ and real <angular frequency>. The ratio is negative exactly on $\omega_d<\omega<\omega_n$, giving imaginary <wavenumbers> and an <evanescent wave>. Thus \b[the propagating frequency gap is $(\omega_d,\omega_n)$]. At $\omega_d$ there is no finite-$k$ solution, while $\omega_n$ has $k=0$; neither endpoint is an ordinary finite-wavelength propagating mode.

Choose the branch $k=(\omega/c)\sqrt{(\omega-\omega_n)/(\omega-\omega_d)}$, so the <phase velocity> is positive in the allowed bands. Direct differentiation gives the <phase velocity> and <group velocity>
$$
v_p=c\sqrt{\frac{\omega-\omega_d}{\omega-\omega_n}},\qquad
v_g=\frac{v_p}{1+\dfrac{\omega(\omega_n-\omega_d)}{2(\omega-\omega_n)(\omega-\omega_d)}}.
$$
The opposite spatial branch reverses both signs. On the chosen branch, the requested <limits> are
$$
\boxed{\begin{array}{c|cc}
\text{frequency limit}&v_p&v_g\\\hline
\omega\to0&c\sqrt{\omega_d/\omega_n}&c\sqrt{\omega_d/\omega_n}\\
|\omega|\to\infty&c&c\\
\omega\to\omega_n^+&+\infty&0^+\\
\omega\to\omega_d^-&0^+&0^+
\end{array}}
$$
The band-edge limits are taken from their propagating sides. More precisely, with $\epsilon=\omega-\omega_n\downarrow0$ and $\delta=\omega_d-\omega\downarrow0$,
$$
v_p\sim c\sqrt{\frac{\omega_n-\omega_d}{\epsilon}},\quad
v_g\sim\frac{2c}{\omega_n}\sqrt{(\omega_n-\omega_d)\epsilon};\qquad
v_p\sim c\sqrt{\frac\delta{\omega_n-\omega_d}},\quad
v_g\sim\frac{2c\delta^{3/2}}{\omega_d\sqrt{\omega_n-\omega_d}}.
$$

For the final case, set $K=k_0c>0$. Substitution of $T(t)e^{ik_0x}$ at $\omega_d=0$ gives $T^{(3)}+(3iK/2)T''+K^2T'=0$. Its <characteristic polynomial> factors as $r(r+2iK)(r-iK/2)$, hence
$$
T=A+Be^{-2iKt}+Ce^{iKt/2}.
$$
The initial data give $A+B+C=1$, $-2B+C/2=0$ and $-4B-C/4=-5$. Thus $B=1$, $C=4$, $A=-4$, and
$$
\boxed{\phi(x,t)=e^{ik_0x}\left[-4+e^{-2ik_0ct}+4e^{ik_0ct/2}\right].}
$$
The time-independent contribution is essential: <static modes lost by cancelling frequency> would be missed if one divided the $\omega_d=0$ <dispersion relation> by $\omega$ and retained only its two nonzero-frequency roots.