On the one hand, from the product representation (we disregard $f \equiv 0$, where the existence of a square root is trivial)
$$f(z) = e^{a(z)} z^{2\nu_0} \prod \left(1-\tfrac{z}{b_k}\right)^{2\nu_k}e^{p_k(z)}\tag{1}$$
of $f$, we directly obtain what should be a product representation of a square root of $f$:
$$g(z) = e^{a(z)/2} z^{\nu_0} \prod\left(1-\tfrac{z}{b_k}\right)^{\nu_k} e^{p_k(z)/2},\tag{2}$$
and then it remains to see that the product in $(2)$ converges locally uniformly on all of $\mathbb{C}$. That is maybe a little tedious if done rigorously.
We can avoid that tedium if we use the Weierstraß product theorem in a different way:
There is an entire function $g_1$ that has zeros of order $\nu_k$ in the $b_k$ (with $b_0 = 0$ and possibly $\nu_0 = 0$), and no other zeros. Then the function
$$q(z) = \frac{f(z)}{g_1(z)^2}$$
is entire and has no zeros. Hence $q$ has a logarithm, $q(z) = e^{g_2(z)}$. Now it is clear that $g(z) = g_1(z) e^{g_2(z)/2}$ is an entire square root of $f$.
An alternative way to establish the existence of a square root of $f$ considers the logarithmic derivative
$$h(z) = \frac{f'(z)}{f(z)}.$$
$h$ is an entire meromorphic function, with simple poles in the zeros of $f$, where the residue in $b_k$ is $2\nu_k$, and no other poles. Therefore the residue of $\tilde{h} = \frac{1}{2}\cdot h$ in all poles is an integer, and hence
$$g(z) = \exp \left(c_0 + \int_{z_0}^z \tilde{h}(\zeta)\,d\zeta\right),$$
where
- $z_0$ is an arbitrary point such that $f(z_0) \neq 0$,
- $c_0$ is chosen so that $e^{2c_0} = f(z_0)$,
- and the integral is over an arbitrary piecewise differentiable path from $z_0$ to $z$ that does not pass through any zero of $f$,
is a well-defined holomorphic function on $\mathbb{C}\setminus f^{-1}(\{0\})$ with $g(z)^2 \equiv f(z)$ there. From that, it follows that the zeros of $f$ are removable singularities of $g$, and the identity $g(z)^2 \equiv f(z)$ holds on all of $\mathbb{C}$ after removing the removable singularities.
Note: the same arguments work on any simply connected domain, but on a domain that is not simply connected, not all holomorphic functions with all zeros of even order have a square root. Also, the arguments, mutatis mutandis, also work for any $m$-th root of $f$ if the order of all zeros of $f$ is a multiple of $m$.
Best Answer
To add on to what Makuasi stated an entire function has a pole of order $n$ at $\infty$ iff the function is a polynomial of order $n$. So if $f$ does not have a pole at $\infty$ $f$ is either bounded or has an essential singularity at $\infty$. $f$ cannot be bounded by Liouville's Thm. So $f$ must have an essential singularity at $\infty$. By Casorati–Weierstrass theorem there is a sequence, $(z_{n})$, such that $|z_{n}| \rightarrow \infty$, and $f(z_{n})\rightarrow 0$