Squared magnitude of the Z-transform
I understand your confusion. What we are looking for is a function of $z$ which equals the squared magnitude $|H(e^{j\omega})|^2$ of a Fourier transform $H(e^{j\omega})=\mathcal{F}\{h[n]\}$ when evaluated on the unit circle $|z|=1$. Note that the inverse discrete-time Fourier transform of $|H(e^{j\omega})|^2$ is given by the convolution $h[n]\star h^*[-n]$, where $^*$ denotes complex conjugation. Since the $\mathcal{Z}$-transform of $h^*[-n]$ is given by $H^*(1/z^*)$ we see that the $\mathcal{Z}$-transform of $h[n]\star h^*[-n]$ is given by
$$\mathcal{Z}\big\{h[n]\star h^*[-n]\big\}=H(z)H^*\left(\frac{1}{z^*}\right)\tag{1}$$
If $h[n]$ is real-valued we obviously have $h^*[-n]=h[-n]$, and $(1)$ is equivalent to
$$\mathcal{Z}\big\{h[n]\star h[-n]\big\}=H(z)H\left(\frac{1}{z}\right)\tag{2}$$
It's important to understand that the equations
$$\big|H(z)\big|^2=H(z)H^*\left(\frac{1}{z^*}\right)\tag{3}$$
and, for real-valued $h[n]$,
$$\big|H(z)\big|^2=H(z)H\left(\frac{1}{z}\right)\tag{4}$$
are only valid on the unit circle $|z|=1$.
In general we of course have
$$\big|H(z)\big|^2=H(z)H^*(z)\tag{5}$$
but $(5)$ is no valid $\mathcal{Z}$-transform of any time domain sequence, because it is not analytic. On the other hand, $(3)$ and $(4)$ are valid $\mathcal{Z}$-transforms, namely of the sequence $h[n]\star h^*[-n]$. On the unit circle they equal the squared magnitude of the Fourier transform of $h[n]$, i.e.,
$$H(z)H^*\left(\frac{1}{z^*}\right)=\big|H(e^{j\omega})\big|^2,\qquad |z|=1\tag{6}$$
Let's look at a simple example. With $h[n]=\delta[n]+\delta[n-1]$, we have $H(z)=1+z^{-1}$. Since $h[n]$ is real-valued we can use Eq. $(4)$:
$$H(z)H\left(\frac{1}{z}\right)=(1+z^{-1})(1+z)=z+2+z^{-1}\tag{7}$$
The actual squared magnitude of $H(z)$ is
$$\big|H(z)\big|^2=|1+z^{-1}|^2=1+2\textrm{Re}\{z^{-1}\}+|z|^{-2}\tag{8}$$
On the unit circle $z=e^{j\omega}$ both $(7)$ and $(8)$ equal
$$\big|H(e^{j\omega})\big|^2=2+2\cos(\omega)\tag{9}$$
However, the inverse $\mathcal{Z}$-transform of $(7)$ equals the inverse Fourier transform of $|H(e^{j\omega})|^2$, whereas $(8)$ is no $\mathcal{Z}$-transform.