[proofplan]
We use the Gaussian convolution formula for the heat semigroup. First we verify that, for each fixed $t>0$ and $x \in \mathbb{R}^n$, the convolution integral defining $(e^{t\Delta}f)(x)$ is finite by pairing the heat kernel with $f$ using Hölder's inequality, with the endpoint cases included. Then we use the pointwise positivity of the heat kernel and the a.e. nonnegativity of $f$ to show that the integrand is nonnegative up to a null set, so its integral is nonnegative. Since this produces a nonnegative representative of $e^{t\Delta}f$, the $L^p$ equivalence class is nonnegative almost everywhere.
[/proofplan]
[step:Fix a nonnegative representative outside a null set]
Fix $t>0$. Since $f \ge 0$ $\mathcal{L}^n$-a.e., there exists a Borel set $N \subset \mathbb{R}^n$ such that $\mathcal{L}^n(N)=0$ and $f(y) \ge 0$ for every $y \in \mathbb{R}^n \setminus N$, after modifying $f$ on a null set if necessary. This modification does not change the element of $L^p(\mathbb{R}^n)$ represented by $f$.
For each $x \in \mathbb{R}^n$, define the translated heat kernel $K_x: \mathbb{R}^n \to (0,\infty)$ by
\begin{align*}
K_x(y) := \Gamma_t(x-y).
\end{align*}
The formula for $\Gamma_t$ gives $K_x(y)>0$ for every $x,y \in \mathbb{R}^n$.
[/step]
[step:Verify that the heat convolution is finite for every spatial point]
Let $x \in \mathbb{R}^n$. We prove that
\begin{align*}
u_t(x) := \int_{\mathbb{R}^n} K_x(y) f(y)\, d\mathcal{L}^n(y)
\end{align*}
is finite as an absolutely convergent integral.
If $1 \le p < \infty$, let $q \in [1,\infty]$ be the Hölder conjugate exponent, so $1/p+1/q=1$, with $q=\infty$ when $p=1$. The [translation invariance of Lebesgue measure](/theorems/4911) gives $K_x \in L^q(\mathbb{R}^n)$ because $\Gamma_t \in L^q(\mathbb{R}^n)$. Applying Hölder's inequality to the functions $K_x$ and $f$ gives
\begin{align*}
\int_{\mathbb{R}^n} |K_x(y) f(y)|\, d\mathcal{L}^n(y) \le \|K_x\|_{L^q(\mathbb{R}^n)} \|f\|_{L^p(\mathbb{R}^n)} < \infty.
\end{align*}
If $p=\infty$, then $f \in L^\infty(\mathbb{R}^n)$ and $K_x \in L^1(\mathbb{R}^n)$. Hence
\begin{align*}
\int_{\mathbb{R}^n} |K_x(y) f(y)|\, d\mathcal{L}^n(y) \le \|f\|_{L^\infty(\mathbb{R}^n)} \int_{\mathbb{R}^n} K_x(y)\, d\mathcal{L}^n(y) < \infty.
\end{align*}
Thus $u_t: \mathbb{R}^n \to \mathbb{R}$ is defined pointwise by the heat convolution formula.
[guided]
Fix a point $x \in \mathbb{R}^n$. Before using positivity, we must first know that the convolution integral is a genuine finite real number. Define $u_t(x)$ by
\begin{align*}
u_t(x) := \int_{\mathbb{R}^n} K_x(y) f(y)\, d\mathcal{L}^n(y),
\end{align*}
where $K_x: \mathbb{R}^n \to (0,\infty)$ is the map $K_x(y)=\Gamma_t(x-y)$.
Assume first that $1 \le p < \infty$. Let $q \in [1,\infty]$ be determined by
\begin{align*}
\frac{1}{p}+\frac{1}{q}=1.
\end{align*}
The heat kernel $\Gamma_t$ belongs to $L^q(\mathbb{R}^n)$ for every $1 \le q \le \infty$: for finite $q$ this follows from Gaussian integrability, and for $q=\infty$ it follows from the bound $\Gamma_t(z) \le (4\pi t)^{-n/2}$. Since $K_x$ is a translate of $\Gamma_t$, translation invariance of $\mathcal{L}^n$ gives
\begin{align*}
\|K_x\|_{L^q(\mathbb{R}^n)} = \|\Gamma_t\|_{L^q(\mathbb{R}^n)} < \infty.
\end{align*}
Now apply Hölder's inequality to the product of $K_x \in L^q(\mathbb{R}^n)$ and $f \in L^p(\mathbb{R}^n)$. This gives
\begin{align*}
\int_{\mathbb{R}^n} |K_x(y) f(y)|\, d\mathcal{L}^n(y) \le \|K_x\|_{L^q(\mathbb{R}^n)} \|f\|_{L^p(\mathbb{R}^n)} < \infty.
\end{align*}
For the endpoint $p=\infty$, the correct pairing is $L^1$ with $L^\infty$. Since $K_x$ is a translate of the heat kernel, $K_x \in L^1(\mathbb{R}^n)$, and therefore
\begin{align*}
\int_{\mathbb{R}^n} |K_x(y) f(y)|\, d\mathcal{L}^n(y) \le \|f\|_{L^\infty(\mathbb{R}^n)} \int_{\mathbb{R}^n} K_x(y)\, d\mathcal{L}^n(y) < \infty.
\end{align*}
Thus the heat convolution is absolutely integrable for every spatial point $x$.
[/guided]
[/step]
[step:Integrate a nonnegative representative of the heat convolution]
For every $x \in \mathbb{R}^n$ and every $y \in \mathbb{R}^n \setminus N$, we have $K_x(y)>0$ and $f(y)\ge 0$, hence
\begin{align*}
K_x(y) f(y) \ge 0.
\end{align*}
Since $\mathcal{L}^n(N)=0$ and $K_x f \in L^1(\mathbb{R}^n)$ by the previous step, the integral over $\mathbb{R}^n$ is the same as the integral over $\mathbb{R}^n \setminus N$. Therefore
\begin{align*}
u_t(x) = \int_{\mathbb{R}^n} K_x(y) f(y)\, d\mathcal{L}^n(y) \ge 0.
\end{align*}
Thus the convolution representative $u_t$ of $e^{t\Delta}f$ is nonnegative for every $x \in \mathbb{R}^n$.
Since equality of representatives in $L^p(\mathbb{R}^n)$ is understood up to $\mathcal{L}^n$-null sets, this pointwise nonnegative representative implies
\begin{align*}
e^{t\Delta}f \ge 0
\end{align*}
$\mathcal{L}^n$-a.e. on $\mathbb{R}^n$. Because $t>0$ was arbitrary, the conclusion holds for every $t>0$.
[/step]
