Analytical way to generate an IIR that is the square root response of another IIR
The only time you'll be able to do this exactly is when you know that the poles and zeros of $H(z)$ are in complex conjugate pairs inside and outside of the unit circle.

Note that for systems with real-valued impulse responses, this means you need four of each: two inside and two outside the unit circle.

The keyword to search on is spectral factorization.

This means that if you can write $H(z)$ as:
$$
H(z) = G(z) G(z^{-1})
$$
where $G(z)$ is a minimum phase system, then $G$ just takes all the poles and zeros of $H$ that are inside the unit circle (and half those that are on the unit circle).

There exist approximate algorithms for finding a solution, but I am not well-versed in those.