Squared magnitude of System Function H(z)
It's important to understand that the equation
$$|H(z)|^2=H(z)H\left(\frac{1}{z}\right)\tag{1}$$
is only valid on the unit circle, i.e., for $z=e^{j\omega}$. For complex-valued systems the general form of $(1)$ is
$$|H(z)|^2=H(z)H^*\left(\frac{1}{z^*}\right)\tag{2}$$
For real-valued systems we have $H(z)=H^*(z^*)$ and, consequently, $(1)$ is valid for real-valued systems. Note that both $(1)$ and $(2)$ are only valid on the unit circle.
The reason why we use $(1)$ or $(2)$ as a generalization of $|H(e^{j\omega})|^2$ to the whole complex plane is that we want an analytic function that is actually a $\mathcal{Z}$-transform of some sequence. And the function $H(z)H^*(z)$ does not satisfy that requirement.
Note that the IDFT of $H^*(e^{j\omega})$ is $h^*[-n]$. So if we find the $\mathcal{Z}$-transform of $h^*[-n]$ we have what we need:
$$\sum_nh^*[-n]z^{-n}=\left(\sum_nh[n](z^*)^{n}\right)^*=H^*\left(\frac{1}{z^*}\right)\tag{3}$$