[proofplan]
The proof proceeds in three stages. First, we establish that the convergence exponent of the zero sequence is at most $\rho$, which guarantees that the canonical product $P(z) = \prod_{k=1}^\infty E_p(z/z_k)$ converges and defines an entire function of order at most $\rho$ with exactly the same zeros as $f$ (with multiplicity). Second, we form the quotient $h(z) = f(z) / (z^m P(z))$, which is entire and zero-free, hence of the form $h(z) = e^{g(z)}$ for some entire function $g$. Third, we show that $g$ must be a polynomial of degree at most $p = \lfloor \rho \rfloor$ by proving that $h$ has order at most $\rho$, which constrains the growth of $g$.
[/proofplan]
[step:Establish that the convergence exponent of the zeros is at most $\rho$]
Let $n(r)$ denote the number of zeros of $f$ (counted with multiplicity) in the closed disc $\overline{B}(0, r)$. Since $f$ has finite order $\rho$, for every $\varepsilon > 0$ there exists $R_\varepsilon > 0$ such that
\begin{align*}
|f(z)| \le e^{|z|^{\rho + \varepsilon}} \quad \text{for all } |z| \ge R_\varepsilon.
\end{align*}
We claim that $n(r) \le C_\varepsilon\, r^{\rho + \varepsilon}$ for $r$ sufficiently large and some constant $C_\varepsilon > 0$. To see this, apply Jensen's formula: for $r > 0$ such that no zeros lie on $|z| = r$,
\begin{align*}
\sum_{|z_k| \le r} \log \frac{r}{|z_k|} = \frac{1}{2\pi} \int_0^{2\pi} \log |f(re^{i\theta})| \, d\mathcal{L}^1(\theta) - \log |c_m|,
\end{align*}
where $c_m$ is the leading coefficient of the Taylor expansion of $f$ at the origin (after factoring out $z^m$). Since $\log |f(re^{i\theta})| \le r^{\rho + \varepsilon}$ for $r \ge R_\varepsilon$, the right-hand side is bounded by $r^{\rho + \varepsilon} + C'$ for a constant $C'$. Restricting the sum on the left to zeros in $\overline{B}(0, r/2)$, each term satisfies $\log(r/|z_k|) \ge \log 2$, so
\begin{align*}
n(r/2) \log 2 \le r^{\rho + \varepsilon} + C'.
\end{align*}
Hence $n(r) = O(r^{\rho + \varepsilon})$ as $r \to \infty$, for every $\varepsilon > 0$.
It follows that the convergence exponent $\lambda := \inf\!\left\{\alpha > 0 : \sum_{k=1}^\infty |z_k|^{-\alpha} < \infty\right\}$ satisfies $\lambda \le \rho$. Indeed, for $\alpha > \rho$, choose $\varepsilon > 0$ with $\rho + \varepsilon < \alpha$. Then $n(r) \le C r^{\rho + \varepsilon}$, and a Stieltjes integration gives
\begin{align*}
\sum_{k=1}^\infty \frac{1}{|z_k|^\alpha} = \alpha \int_0^\infty \frac{n(r)}{r^{\alpha + 1}} \, d\mathcal{L}^1(r) \le C \alpha \int_{r_0}^\infty r^{\rho + \varepsilon - \alpha - 1} \, d\mathcal{L}^1(r) < \infty,
\end{align*}
since $\rho + \varepsilon - \alpha < 0$. Therefore $\lambda \le \alpha$, and since $\alpha > \rho$ was arbitrary, $\lambda \le \rho$.
In particular, since $p = \lfloor \rho \rfloor \ge \lambda - 1$, we have $p + 1 > \lambda$, so $\sum_{k=1}^\infty |z_k|^{-(p+1)} < \infty$.
[guided]
The order $\rho$ of $f$ controls how fast $|f(z)|$ can grow: by definition, $\rho = \inf\{\alpha \ge 0 : |f(z)| \le e^{|z|^\alpha} \text{ for } |z| \text{ large enough}\}$. For every $\varepsilon > 0$, we therefore have $|f(z)| \le e^{|z|^{\rho+\varepsilon}}$ for $|z|$ sufficiently large.
Why does this constrain the density of zeros? The intuition is that a function with many zeros in a disc must be small somewhere on the disc (by the intermediate value theorem for modulus), and if the function also grows slowly, it cannot have too many zeros. The precise tool is Jensen's formula, which relates the zeros of $f$ inside $B(0, r)$ to the average of $\log |f|$ on $\partial B(0, r)$:
\begin{align*}
\sum_{|z_k| \le r} \log \frac{r}{|z_k|} = \frac{1}{2\pi} \int_0^{2\pi} \log |f(re^{i\theta})| \, d\mathcal{L}^1(\theta) - \log |c_m|,
\end{align*}
where $c_m$ is the leading nonzero Taylor coefficient. The right-hand side is at most $r^{\rho+\varepsilon} + C'$ for $r$ large. To extract a count of zeros from the left-hand side, restrict to zeros in $\overline{B}(0, r/2)$. For such zeros, $|z_k| \le r/2$, so $\log(r/|z_k|) \ge \log 2$. Therefore
\begin{align*}
n(r/2) \log 2 \le \sum_{|z_k| \le r/2} \log \frac{r}{|z_k|} \le \sum_{|z_k| \le r} \log \frac{r}{|z_k|} \le r^{\rho + \varepsilon} + C',
\end{align*}
giving $n(r) = O(r^{\rho + \varepsilon})$ after replacing $r/2$ by $r$.
The convergence exponent $\lambda$ of the zero sequence is defined as the infimum of all $\alpha > 0$ for which $\sum_k |z_k|^{-\alpha} < \infty$. The bound $n(r) = O(r^{\rho + \varepsilon})$ is converted to a bound on the Dirichlet series via Stieltjes integration:
\begin{align*}
\sum_{k=1}^\infty \frac{1}{|z_k|^\alpha} = \alpha \int_0^\infty \frac{n(r)}{r^{\alpha+1}} \, d\mathcal{L}^1(r).
\end{align*}
For $\alpha > \rho + \varepsilon$, the integrand decays as $r^{\rho + \varepsilon - \alpha - 1}$ with exponent less than $-1$, so the integral converges. Since $\varepsilon > 0$ was arbitrary, $\lambda \le \rho$.
The consequence $p + 1 > \lambda$ (where $p = \lfloor \rho \rfloor$) ensures that the Weierstrass elementary factors $E_p(z/z_k)$ will produce a convergent infinite product.
[/guided]
[/step]
[step:Construct the canonical product $P(z)$ and verify it has the correct zeros and finite order]
Define the canonical product
\begin{align*}
P: \mathbb{C} &\to \mathbb{C} \\
z &\mapsto \prod_{k=1}^\infty E_p\!\left(\frac{z}{z_k}\right).
\end{align*}
Since $\sum_{k=1}^\infty |z_k|^{-(p+1)} < \infty$, the standard estimate on Weierstrass elementary factors guarantees that for each compact subset $K \subset \mathbb{C}$,
\begin{align*}
\sum_{k=1}^\infty \left|1 - E_p\!\left(\frac{z}{z_k}\right)\right| < \infty \quad \text{uniformly for } z \in K.
\end{align*}
Therefore the infinite product converges uniformly on compact subsets of $\mathbb{C}$, and $P$ is entire. The zeros of $P$ are exactly $z_1, z_2, \ldots$ with the correct multiplicities (each factor $E_p(z/z_k)$ vanishes to order $1$ at $z = z_k$ and nowhere else). Moreover, the order of $P$ is at most $\lambda \le \rho$.
[guided]
The Weierstrass elementary factor $E_p(u) = (1-u)\exp(u + u^2/2 + \cdots + u^p/p)$ is designed so that $E_p(u) = 1 + O(|u|^{p+1})$ as $u \to 0$. More precisely, $\log E_p(u) = -\sum_{j=p+1}^\infty u^j/j$ for $|u| < 1$, and consequently $|1 - E_p(u)| \le C|u|^{p+1}$ for $|u| \le 1/2$. This means the convergence of the product $\prod_k E_p(z/z_k)$ is controlled by $\sum_k |z/z_k|^{p+1}$. For $z$ in a compact set $K$, we have $|z| \le R$ for some $R$, and for $k$ large enough $|z_k| > 2R$, so $|z/z_k| < 1/2$ and
\begin{align*}
\sum_{k \text{ large}} |1 - E_p(z/z_k)| \le C R^{p+1} \sum_{k \text{ large}} |z_k|^{-(p+1)} < \infty.
\end{align*}
The convergence $\sum_k |z_k|^{-(p+1)} < \infty$, established in the previous step, is precisely what is needed. By the theory of infinite products, uniform convergence of $\sum_k |1 - E_p(z/z_k)|$ on compact sets implies that the product defines an entire function whose zeros are exactly the $z_k$ with the correct multiplicities.
The order of $P$ is bounded by the convergence exponent $\lambda$ of the zero sequence (this is a standard result in the theory of entire functions, following from the estimate $\log |E_p(u)| \le C_p |u|^{p+1}$ for $|u| \le 1/2$ and $\log |E_p(u)| \le C_p |u|^p \log |u|$ for $|u| > 1/2$). Since $\lambda \le \rho$, the product $P$ has order at most $\rho$.
[/guided]
[/step]
[step:Form the quotient $f(z)/(z^m P(z))$ and express it as $e^{g(z)}$]
Define
\begin{align*}
h: \mathbb{C} &\to \mathbb{C} \\
z &\mapsto \frac{f(z)}{z^m P(z)},
\end{align*}
where $m \ge 0$ is the order of the zero of $f$ at the origin. The numerator $f$ vanishes at $z = 0$ to order $m$ and at each $z_k$ to order equal to the multiplicity of $z_k$ as a zero of $f$. The denominator $z^m P(z)$ has the same zeros with the same multiplicities. Therefore $h$ extends to an entire function (the removable singularities at the zeros are removed by construction), and $h$ is zero-free on $\mathbb{C}$.
Since $h$ is entire and zero-free on the simply connected domain $\mathbb{C}$, there exists an entire function $g: \mathbb{C} \to \mathbb{C}$ such that
\begin{align*}
h(z) = e^{g(z)} \quad \text{for all } z \in \mathbb{C}.
\end{align*}
This gives the factorisation $f(z) = z^m e^{g(z)} P(z) = z^m e^{g(z)} \prod_{k=1}^\infty E_p(z/z_k)$.
[guided]
The point of introducing $h$ is to separate the "zero structure" of $f$ (captured by $z^m P(z)$) from its "growth beyond what the zeros dictate" (captured by $e^{g(z)}$).
Why is $h$ entire? At each zero $z_k$ of $f$, the factor $E_p(z/z_k)$ in $P(z)$ contributes a simple zero. If $z_k$ appears with multiplicity $\nu$ in the list (i.e., $z_k = z_{k+1} = \cdots = z_{k+\nu-1}$), then $P(z)$ vanishes at $z_k$ to order $\nu$, matching the multiplicity of $z_k$ as a zero of $f$. At $z = 0$, the factor $z^m$ cancels the $m$-fold zero. So $f(z)/(z^m P(z))$ has removable singularities at each zero, and extends to an entire function by Riemann's removable singularity theorem.
Why is $h$ zero-free? If $h(z_0) = 0$ for some $z_0 \in \mathbb{C}$, then $f(z_0) = 0$ but $z_0$ is not a zero of $z^m P(z)$, contradicting the fact that $z^m P(z)$ accounts for all zeros of $f$.
Why can we write $h = e^g$ for an entire $g$? A zero-free holomorphic function on a simply connected domain always admits a holomorphic logarithm. Since $\mathbb{C}$ is simply connected and $h$ is zero-free, we define $g(z) = \int_0^z h'(w)/h(w) \, dw + \log h(0)$ (choosing any branch of $\log h(0)$). Then $g$ is entire and $h = e^g$.
[/guided]
[/step]
[step:Show that $g$ is a polynomial of degree at most $p = \lfloor \rho \rfloor$]
It remains to show that $g$ is a polynomial of degree at most $p$. Since $f$ has order $\rho$ and $z^m P(z)$ has order at most $\rho$ (as established above), the quotient $h = f/(z^m P)$ also has order at most $\rho$. That is, for every $\varepsilon > 0$ and $|z|$ sufficiently large:
\begin{align*}
|h(z)| \le e^{|z|^{\rho + \varepsilon}}.
\end{align*}
Since $h = e^g$, this gives $\operatorname{Re}(g(z)) = \log |h(z)| \le |z|^{\rho + \varepsilon}$ for $|z|$ large. By the Borel-Carath\'eodory theorem, a bound on $\operatorname{Re}(g)$ on circles implies a bound on $|g|$: for $|z| = r$ and $R = 2r$,
\begin{align*}
|g(z)| \le \frac{2R}{R - r} \sup_{|w| = R} \operatorname{Re}(g(w)) + \frac{R + r}{R - r} |g(0)| \le 4 \sup_{|w| = 2r} \operatorname{Re}(g(w)) + 3|g(0)|.
\end{align*}
Since $\operatorname{Re}(g(w)) \le |w|^{\rho + \varepsilon} \le (2r)^{\rho + \varepsilon}$ for $|w| = 2r$ sufficiently large, we obtain $|g(z)| \le C_\varepsilon r^{\rho + \varepsilon}$ for $|z| = r$. This means $g$ is an entire function of order at most $\rho$.
By the [Generalised Liouville Theorem](/theorems/3342), an entire function satisfying $|g(z)| \le C_\varepsilon(1 + |z|^{\rho + \varepsilon})$ for every $\varepsilon > 0$ is a polynomial. More precisely: since $g$ is entire and of order at most $\rho$, the Taylor coefficients $b_k = g^{(k)}(0)/k!$ satisfy $|b_k| \le M_R / R^k$ by the [Cauchy Estimates](/theorems/2571). Taking $M_R \le e^{R^{\rho + \varepsilon}}$ and optimising $R$ shows that $b_k = 0$ for $k > p = \lfloor \rho \rfloor$. Therefore $g$ is a polynomial of degree at most $p$.
[guided]
This is the most delicate step. We need to show that the "leftover" entire function $g$ (the logarithm of the zero-free part) is actually a polynomial, and moreover of the right degree.
The chain of reasoning is:
1. $f$ has order $\rho$ and $z^m P(z)$ has order at most $\rho$, so $h = f/(z^m P)$ has order at most $\rho$.
2. Since $h = e^g$, we have $\operatorname{Re}(g(z)) = \log|h(z)|$, so $\operatorname{Re}(g)$ grows at most as fast as $|z|^{\rho + \varepsilon}$.
3. The Borel-Carath\'eodory theorem converts a bound on $\operatorname{Re}(g)$ into a bound on $|g|$, showing $g$ has order at most $\rho$.
For point 1: more precisely, the order of a quotient of entire functions satisfies $\rho(f_1/f_2) \le \max(\rho(f_1), \rho(f_2))$ when $f_1/f_2$ is entire. We apply this with $f_1 = f$ (order $\rho$) and $f_2 = z^m P$ (order at most $\rho$).
For point 3: the Borel-Carath\'eodory theorem states that if $g$ is holomorphic on $\overline{B}(0, R)$, then for $|z| = r < R$,
\begin{align*}
|g(z)| \le \frac{2r}{R - r} \sup_{|w| = R} \operatorname{Re}(g(w)) + \frac{R + r}{R - r}|g(0)|.
\end{align*}
Taking $R = 2r$, the coefficients become $2r/r = 2 \cdot 2 = 4$ and $(2r + r)/r = 3$, giving
\begin{align*}
|g(z)| \le 4 \sup_{|w| = 2r} \operatorname{Re}(g(w)) + 3|g(0)| \le 4(2r)^{\rho + \varepsilon} + 3|g(0)| \le C_\varepsilon r^{\rho + \varepsilon}
\end{align*}
for $r$ sufficiently large. So $g$ is entire of order at most $\rho$.
Now, why must an entire function of order at most $\rho$ that is also a polynomial in disguise (since $e^g$ has finite order) actually be a polynomial of degree at most $\lfloor \rho \rfloor$? If $g(z) = \sum_{k=0}^\infty b_k z^k$ were not a polynomial, then $e^g$ would have infinite order (since $e^{z^N}$ has order $N$ for integer $N$, and the presence of infinitely many nonzero terms in $g$ would produce at least that growth). More precisely, if $g$ has degree $d$ as a polynomial, then $e^g$ has order exactly $d$. Since $e^g = h$ has order at most $\rho$, we need $d \le \rho$, hence $d \le \lfloor \rho \rfloor = p$.
[/guided]
[/step]
[step:Assemble the Hadamard factorisation]
Combining the results of the previous steps, we have
\begin{align*}
f(z) = z^m e^{g(z)} \prod_{k=1}^\infty E_p\!\left(\frac{z}{z_k}\right),
\end{align*}
where $m \ge 0$ is the multiplicity of the zero at the origin, $g$ is a polynomial of degree at most $p = \lfloor \rho \rfloor$, and the infinite product converges uniformly on compact subsets of $\mathbb{C}$. This is the Hadamard factorisation.
[/step]