Is the group velocity times the phase velocity always the velocity of light squared?

For a classical relativistic scalar field, we have $$E(p) = sqrt{(pc)^2 + left( mc^2 right)^2},$$ and using $E = hbar omega$ and $p = hbar k$ we find that (in units with $hbar = c = 1$) $$omega(k) = sqrt{k^2 + m^2},$$ so a wave with wave number $k$ has phase velocity $$frac{omega}{k} = sqrt{1 + left(frac{m}{k} right)^2}$$ and group velocity $$frac{domega}{dk} = frac{1}{sqrt{1 + left(frac{m}{k}right)^2}}.$$

So indeed, such a wave has that the product of the group and phase velocities equals $c^2$ for any wave number $k$.

But this only holds for relativistic scalar fields; in general waves can have pretty much any dispersion relation $omega(k)$ and so any relation between the group and phase velocities. (For example, waves like sounds waves or water waves obviously don't have any direct connection to the speed of light.)

$$c^2~stackrel{?}{=}~v_pv_g~=~frac{E}{p}frac{d E}{d p}~=~frac{d (E^2)}{d (p^2)}qquadLeftrightarrow qquad E^2-c^2p^2 ~=~{rm const} $$ is true for relativistic scalar matter, but is violated e.g. for non-relativistic matter waves $$ E~=~frac{p^2}{2m}qquadRightarrow qquad v_g~=~2v_p .$$

No.

$v_g v_p=c^2$ is true for waves in waveguides. But that's just one area where $v_g$ and $v_p$ are different, and that happens whenever the velocity depends on the wavelength (or frequency). In other areas - matter waves, and even the basic dispersion of colours by a prism - it does not apply.

This is quite a common trap for students. If they first meet group velocity when learning about waveguides they assume that $v_g v_p=c^2$ is fundamental. It isn't: it happens to be true in this case but not in general.