[guided]The goal is to extend the real-time formula $S(t) = e^{-tA}$ from $t > 0$ to complex time $z$ with positive real part. The spectral calculus reduces this operator-norm question to scalar estimates for the functions $\lambda \mapsto e^{-z\lambda}$ on the spectrum of $A$, which is contained in $[0,\infty)$ because $A$ is nonnegative.
Let
\begin{align*}
\mathbb{C}_+ := \{z \in \mathbb{C} : \operatorname{Re} z > 0\}.
\end{align*}
For each $z \in \mathbb{C}_+$, define $m_z: [0,\infty) \to \mathbb{C}$ by $m_z(\lambda)=e^{-z\lambda}$.
This is the correct scalar multiplier because $S(t)$ is $e^{-tA}$ for real $t > 0$. The condition $\operatorname{Re} z > 0$ is exactly what keeps this multiplier bounded on the nonnegative spectrum:
\begin{align*}
|m_z(\lambda)| = |e^{-z\lambda}| = e^{-(\operatorname{Re} z)\lambda} \leq 1
\end{align*}
for every $\lambda \geq 0$. The spectral calculus therefore defines a bounded operator $m_z(A) \in \mathcal{L}(H)$, and its norm satisfies
\begin{align*}
\|m_z(A)\|_{\mathcal{L}(H)} \leq \sup_{\lambda \in [0,\infty)} |m_z(\lambda)| \leq 1.
\end{align*}
We set $S: \mathbb{C}_+ \to \mathcal{L}(H)$ by $S(z)=m_z(A)$.
It remains to prove holomorphicity in the operator norm, not merely in scalar matrix coefficients. Fix $z_0 \in \mathbb{C}_+$, and choose $r > 0$ with $r < \operatorname{Re} z_0$. If $h \in \mathbb{C}$ and $0 < |h| < r/2$, then $z_0+h \in \mathbb{C}_+$, and every point $z_0 + sh$ with $s \in [0,1]$ satisfies
\begin{align*}
\operatorname{Re}(z_0+sh) \geq \operatorname{Re} z_0-r/2.
\end{align*}
Define $b_{z_0}: [0,\infty) \to \mathbb{C}$ by $b_{z_0}(\lambda)=\lambda e^{-z_0\lambda}$. Since $\operatorname{Re} z_0 > 0$, the scalar function $b_{z_0}$ is bounded on $[0,\infty)$, with
\begin{align*}
|b_{z_0}(\lambda)| \leq \lambda e^{-(\operatorname{Re} z_0)\lambda} \leq (e\operatorname{Re} z_0)^{-1}.
\end{align*}
Thus the spectral calculus defines $A S(z_0)=b_{z_0}(A)$ as a bounded operator on $H$.
To compute the derivative, define the difference-quotient error multiplier $d_h: [0,\infty) \to \mathbb{C}$ by
\begin{align*}
d_h(\lambda)=\frac{e^{-(z_0+h)\lambda}-e^{-z_0\lambda}}{h}+\lambda e^{-z_0\lambda}.
\end{align*}
The spectral multiplier estimate gives
\begin{align*}
\left\|\frac{S(z_0+h)-S(z_0)}{h}+A S(z_0)\right\|_{\mathcal{L}(H)} \leq \sup_{\lambda \in [0,\infty)} |d_h(\lambda)|.
\end{align*}
Therefore the operator-norm derivative will follow once this scalar supremum tends to $0$.
We use Taylor's formula with integral remainder for the one-variable holomorphic function $\zeta \mapsto e^{-\zeta\lambda}$ along the segment $\zeta = z_0+sh$, where $s \in [0,1]$. For each $\lambda \geq 0$,
\begin{align*}
e^{-(z_0+h)\lambda} = e^{-z_0\lambda}-h\lambda e^{-z_0\lambda}+h^2\lambda^2\int_0^1 (1-s)e^{-(z_0+sh)\lambda}\,d\mathcal{L}^1(s).
\end{align*}
Substituting this identity into the definition of $d_h$ gives
\begin{align*}
|d_h(\lambda)| \leq |h|\lambda^2 e^{-(\operatorname{Re} z_0-r/2)\lambda}\int_0^1 (1-s)\,d\mathcal{L}^1(s).
\end{align*}
The one-dimensional integral equals $1/2$. If
\begin{align*}
\beta := \operatorname{Re} z_0-r/2,
\end{align*}
then $\beta > 0$, and the function $\lambda \mapsto \lambda^2 e^{-\beta\lambda}$ on $[0,\infty)$ has maximum $4e^{-2}\beta^{-2}$. Consequently,
\begin{align*}
\sup_{\lambda \in [0,\infty)} |d_h(\lambda)| \leq 2e^{-2}\beta^{-2}|h|.
\end{align*}
The right-hand side tends to $0$ as $h \to 0$, so
\begin{align*}
\lim_{h\to 0}\left\|\frac{S(z_0+h)-S(z_0)}{h}+A S(z_0)\right\|_{\mathcal{L}(H)} = 0.
\end{align*}
Thus $S$ is complex differentiable in operator norm at $z_0$, with derivative $S'(z_0)=-A S(z_0)$. Since $z_0$ was arbitrary, $S: \mathbb{C}_+ \to \mathcal{L}(H)$ is holomorphic and extends the real heat semigroup.[/guided]