[proofplan]
We prove both directions of the equivalence. For the forward direction ($H^k z \to 0$ for all $z$ implies $\rho(H) < 1$), we argue by contrapositive: if some eigenvalue $\lambda$ satisfies $|\lambda| \geq 1$, we produce a real vector $z$ for which $H^k z \not\to 0$ by extracting real or imaginary parts from a (possibly complex) eigenvector. For the reverse direction ($\rho(H) < 1$ implies $H^k z \to 0$), we use the Jordan normal form to express $H^k$ in terms of powers $\lambda_j^k$ and polynomial factors in $k$, both of which decay to zero when $|\lambda_j| < 1$.
[/proofplan]
[step:Show that $H^k z \to 0$ for all $z$ forces every eigenvalue to satisfy $|\lambda| < 1$]
We prove the contrapositive: suppose there exists an eigenvalue $\lambda$ of $H$ with $|\lambda| \geq 1$. We produce a real vector $z \in \mathbb{R}^n$ for which $H^k z \not\to 0$.
Let $w \in \mathbb{C}^n \setminus \{0\}$ satisfy $Hw = \lambda w$. Then $H^k w = \lambda^k w$ for all $k \geq 0$, so
\begin{align*}
\|H^k w\|_\infty = |\lambda|^k \|w\|_\infty \geq \|w\|_\infty > 0 \quad \text{for all } k \geq 0,
\end{align*}
since $|\lambda| \geq 1$.
**Case 1:** $\lambda \in \mathbb{R}$. Then $w$ can be chosen real (the eigenequation $(H - \lambda I)w = 0$ has a real solution since $H$ and $\lambda$ are real). Setting $z = w$ gives $\|H^k z\|_\infty \geq \|w\|_\infty > 0$, so $H^k z \not\to 0$.
**Case 2:** $\lambda \notin \mathbb{R}$. Write $w = u + iv$ with $u, v \in \mathbb{R}^n$. Since $H$ has real entries, $H^k w = H^k u + iH^k v$. By the triangle inequality:
\begin{align*}
\|H^k u\|_\infty + \|H^k v\|_\infty \geq \|H^k u + iH^k v\|_\infty = \|H^k w\|_\infty \geq \|w\|_\infty > 0.
\end{align*}
(Here we use $\|a + ib\|_\infty \leq \|a\|_\infty + \|b\|_\infty$ for real vectors $a, b$, applied in the reverse direction.) If both $H^k u \to 0$ and $H^k v \to 0$, then $H^k w = H^k u + iH^k v \to 0$, contradicting $\|H^k w\|_\infty \geq \|w\|_\infty > 0$. Therefore at least one of $z = u$ or $z = v$ satisfies $H^k z \not\to 0$.
In both cases, $H^k z \not\to 0$ for some $z \in \mathbb{R}^n$, so the hypothesis "$H^k z \to 0$ for all $z$" fails. Taking the contrapositive: if $H^k z \to 0$ for all $z \in \mathbb{R}^n$, then $|\lambda| < 1$ for every eigenvalue $\lambda$ of $H$, i.e., $\rho(H) < 1$.
[guided]
The subtlety in this direction is that eigenvalues of a real matrix can be complex, but our hypothesis involves only real vectors $z \in \mathbb{R}^n$. We need to bridge the gap from "eigenvector power does not decay" (which may involve a complex eigenvector) to "some real vector does not decay."
If $\lambda$ is real, the eigenvalue equation $(H - \lambda I)w = 0$ is a real linear system, so it has a nonzero real solution. In that case $z = w$ works directly: $H^k z = \lambda^k z$, which does not tend to zero because $|\lambda|^k \geq 1$.
If $\lambda$ is complex, the eigenvector $w$ is necessarily complex. We decompose $w = u + iv$ into real and imaginary parts. Since $H$ is a real matrix, $H^k$ is also real, so $H^k$ maps real vectors to real vectors and $H^k w = H^k u + i H^k v$.
Now we use the contrapositive of the triangle inequality. For any complex vector $c = a + ib$ with $a, b \in \mathbb{R}^n$:
\begin{align*}
\|c\|_\infty = \max_j |a_j + ib_j| \leq \max_j (|a_j| + |b_j|) \leq \|a\|_\infty + \|b\|_\infty.
\end{align*}
Applying this with $a = H^k u$ and $b = H^k v$: if both $H^k u \to 0$ and $H^k v \to 0$, then $\|H^k u\|_\infty + \|H^k v\|_\infty \to 0$, which forces $\|H^k w\|_\infty \to 0$. But we established $\|H^k w\|_\infty \geq \|w\|_\infty > 0$ for all $k$, a contradiction. So at least one of the real vectors $u$ or $v$ witnesses the failure of $H^k z \to 0$.
[/guided]
[/step]
[step:Show that $\rho(H) < 1$ implies $H^k z \to 0$ for all $z$ via the Jordan normal form]
Assume $\rho(H) < 1$. The Jordan normal form provides an invertible matrix $P \in \mathbb{C}^{n \times n}$ such that
\begin{align*}
H = PJP^{-1}, \qquad J = \operatorname{diag}(J_1, \ldots, J_s),
\end{align*}
where each $J_l$ is a Jordan block of size $m_l$ with eigenvalue $\lambda_l$:
\begin{align*}
J_l = \lambda_l I_{m_l} + N_l, \qquad (N_l)_{ij} = \begin{cases} 1 & \text{if } j = i+1, \\ 0 & \text{otherwise}. \end{cases}
\end{align*}
The nilpotent part satisfies $N_l^{m_l} = 0$. Taking the $k$-th power of a single Jordan block via the binomial theorem (which is valid because $\lambda_l I$ and $N_l$ commute):
\begin{align*}
J_l^k = (\lambda_l I + N_l)^k = \sum_{j=0}^{\min(k, m_l - 1)} \binom{k}{j} \lambda_l^{k-j} N_l^j,
\end{align*}
where the sum terminates at $j = m_l - 1$ because $N_l^{m_l} = 0$. Each term satisfies
\begin{align*}
\left|\binom{k}{j} \lambda_l^{k-j}\right| = \binom{k}{j} |\lambda_l|^{k-j}.
\end{align*}
Since $|\lambda_l| \leq \rho(H) < 1$ and $j \leq m_l - 1$ is bounded, the polynomial growth $\binom{k}{j} \leq k^j / j!$ is dominated by the geometric decay $|\lambda_l|^{k-j}$:
\begin{align*}
\binom{k}{j} |\lambda_l|^{k-j} \leq \frac{k^j}{j!} \cdot |\lambda_l|^{k-j} \to 0 \quad \text{as } k \to \infty,
\end{align*}
because $k^j |\lambda_l|^k \to 0$ for any fixed $j \geq 0$ when $|\lambda_l| < 1$ (exponential decay beats polynomial growth). Therefore $J_l^k \to 0$ for each block, hence $J^k \to 0$, and
\begin{align*}
H^k = PJ^kP^{-1} \to P \cdot 0 \cdot P^{-1} = 0.
\end{align*}
Since $H^k \to 0$ as a matrix (i.e., every entry tends to zero), $H^k z \to 0$ for every $z \in \mathbb{R}^n$.
[guided]
The Jordan normal form is the standard tool for analysing the asymptotic behaviour of matrix powers. Every matrix over $\mathbb{C}$ has a Jordan form $H = PJP^{-1}$, and $H^k = PJ^kP^{-1}$, so the behaviour of $H^k$ reduces to the behaviour of powers of individual Jordan blocks.
A Jordan block $J_l = \lambda_l I + N_l$ decomposes into a scalar part (the eigenvalue $\lambda_l$) and a nilpotent part ($N_l$, the superdiagonal of ones). Since $\lambda_l I$ and $N_l$ commute, the binomial theorem applies:
\begin{align*}
J_l^k = \sum_{j=0}^{m_l - 1} \binom{k}{j} \lambda_l^{k-j} N_l^j.
\end{align*}
The sum terminates at $j = m_l - 1$ because $N_l^{m_l} = 0$ (each application of $N_l$ shifts the superdiagonal up by one row; after $m_l$ applications, everything has been shifted past the top of the matrix).
Now we ask: does each term $\binom{k}{j} \lambda_l^{k-j}$ tend to zero? The binomial coefficient $\binom{k}{j}$ grows polynomially in $k$ (as $k^j / j! + \text{lower order terms}$), while $|\lambda_l|^{k-j} = |\lambda_l|^{-j} \cdot |\lambda_l|^k$ decays exponentially in $k$ since $|\lambda_l| < 1$. Exponential decay dominates polynomial growth: for any fixed $j$ and any $r \in (0,1)$, $k^j r^k \to 0$ as $k \to \infty$ (for instance, by the ratio test: $\frac{(k+1)^j r^{k+1}}{k^j r^k} = r \cdot (1 + 1/k)^j \to r < 1$). Therefore every entry of $J_l^k$ tends to zero, and so $J^k \to 0$.
The conjugation $H^k = PJ^k P^{-1}$ preserves convergence to zero because $P$ and $P^{-1}$ are fixed matrices:
\begin{align*}
\|H^k\| = \|PJ^kP^{-1}\| \leq \|P\| \cdot \|J^k\| \cdot \|P^{-1}\| \to 0.
\end{align*}
Applying $H^k$ to any fixed $z$: $\|H^k z\| \leq \|H^k\| \cdot \|z\| \to 0$.
What if $|\lambda_l| = 1$? Then $|\lambda_l|^{k-j} = 1$ and the polynomial factor $\binom{k}{j}$ is no longer beaten down by exponential decay. For a $1 \times 1$ block, $J_l^k = \lambda_l^k$ has $|J_l^k| = 1$ and does not tend to zero (though it stays bounded). For larger blocks, the polynomial growth makes $J_l^k$ unbounded. This is why the strict inequality $\rho(H) < 1$ is essential.
[/guided]
[/step]