[proofplan]
Fix $t>0$ and write the heat semigroup through the heat-kernel convolution formula. The proof has two parts: first, Tonelli's theorem and the normalization of the heat kernel show that the convolution is absolutely integrable in $x$, so $u(t,\cdot)\in L^1(\mathbb{R}^n)$. Second, [Fubini's theorem](/theorems/2961) permits interchange of the $x$- and $y$-integrals; the inner heat-kernel integral equals $1$, leaving exactly the integral of $f$.
[/proofplan]
[step:Declare the heat-kernel convolution and its normalization]
Fix $t>0$. Define the kernel
\begin{align*}
\Gamma_t: \mathbb{R}^n \to (0,\infty), \qquad z \mapsto (4\pi t)^{-n/2}\exp\left(-\frac{|z|^2}{4t}\right).
\end{align*}
The heat-kernel normalization is
\begin{align*}
\int_{\mathbb{R}^n}\Gamma_t(z)\,d\mathcal{L}^n(z)=1.
\end{align*}
Define
\begin{align*}
u_t: \mathbb{R}^n \to \mathbb{R}, \qquad x \mapsto \int_{\mathbb{R}^n}\Gamma_t(x-y)f(y)\,d\mathcal{L}^n(y),
\end{align*}
where $u_t$ denotes $u(t,\cdot)$.
[guided]
We fix a time $t>0$ because the assertion is pointwise in time. The heat semigroup is represented by convolution with the heat kernel, so the relevant kernel is the map
\begin{align*}
\Gamma_t: \mathbb{R}^n \to (0,\infty), \qquad z \mapsto (4\pi t)^{-n/2}\exp\left(-\frac{|z|^2}{4t}\right).
\end{align*}
The constant $(4\pi t)^{-n/2}$ is chosen precisely so that the heat kernel has total mass one:
\begin{align*}
\int_{\mathbb{R}^n}\Gamma_t(z)\,d\mathcal{L}^n(z)=1.
\end{align*}
For this fixed $t$, define the function
\begin{align*}
u_t: \mathbb{R}^n \to \mathbb{R}, \qquad x \mapsto \int_{\mathbb{R}^n}\Gamma_t(x-y)f(y)\,d\mathcal{L}^n(y).
\end{align*}
This is the convolution $u_t=\Gamma_t*f$, which is the meaning of $u(t,\cdot)=e^{t\Delta}f$ in the statement. The rest of the proof shows that this integral defines an $L^1$ function and that its total integral agrees with the total integral of $f$.
[/guided]
[/step]
[step:Use Tonelli to prove that the convolution is integrable]
Consider the non-negative measurable function
\begin{align*}
H_t: \mathbb{R}^n \times \mathbb{R}^n \to [0,\infty), \qquad (x,y) \mapsto \Gamma_t(x-y)|f(y)|.
\end{align*}
By Tonelli's theorem and the change of variables $z=x-y$, under which $d\mathcal{L}^n(z)=d\mathcal{L}^n(x)$ and $\mathbb{R}^n$ is mapped onto $\mathbb{R}^n$, we have
\begin{align*}
\int_{\mathbb{R}^n}\int_{\mathbb{R}^n}\Gamma_t(x-y)|f(y)|\,d\mathcal{L}^n(x)\,d\mathcal{L}^n(y)=\int_{\mathbb{R}^n}|f(y)|\,d\mathcal{L}^n(y).
\end{align*}
The right-hand side is finite because $f\in L^1(\mathbb{R}^n)$. Hence the function
\begin{align*}
G_t: \mathbb{R}^n \times \mathbb{R}^n \to \mathbb{R}, \qquad (x,y) \mapsto \Gamma_t(x-y)f(y)
\end{align*}
belongs to $L^1(\mathbb{R}^n\times\mathbb{R}^n)$. Therefore, for $\mathcal{L}^n$-a.e. $x\in\mathbb{R}^n$, the integral defining $u_t(x)$ is absolutely convergent. Moreover,
\begin{align*}
\int_{\mathbb{R}^n}|u_t(x)|\,d\mathcal{L}^n(x)\leq \int_{\mathbb{R}^n}\int_{\mathbb{R}^n}\Gamma_t(x-y)|f(y)|\,d\mathcal{L}^n(y)\,d\mathcal{L}^n(x).
\end{align*}
By the same Tonelli computation, the right-hand side equals $\|f\|_{L^1(\mathbb{R}^n)}$. Thus $u_t\in L^1(\mathbb{R}^n)$.
[guided]
We must first prove that the convolution is an $L^1$ function. The natural object to estimate is the absolute-value integrand, so define the non-negative measurable function
\begin{align*}
H_t: \mathbb{R}^n \times \mathbb{R}^n \to [0,\infty), \qquad (x,y) \mapsto \Gamma_t(x-y)|f(y)|.
\end{align*}
The function is measurable because $(x,y) \mapsto x-y$ is continuous, $\Gamma_t$ is continuous, and $f$ is Lebesgue measurable as an $L^1(\mathbb{R}^n)$ function. Since $H_t \geq 0$, Tonelli's theorem applies without any prior integrability assumption.
For each fixed $y\in\mathbb{R}^n$, use the translation change of variables $z=x-y$. Under this substitution, [Lebesgue measure](/page/Lebesgue%20Measure) is preserved, so $d\mathcal{L}^n(z)=d\mathcal{L}^n(x)$, and the domain $\mathbb{R}^n$ in the variable $x$ is mapped onto the domain $\mathbb{R}^n$ in the variable $z$. Therefore the inner integral is
\begin{align*}
\int_{\mathbb{R}^n}\Gamma_t(x-y)|f(y)|\,d\mathcal{L}^n(x)=|f(y)|\int_{\mathbb{R}^n}\Gamma_t(z)\,d\mathcal{L}^n(z).
\end{align*}
By the heat-kernel normalization, the last integral equals $1$. Tonelli's theorem therefore gives
\begin{align*}
\int_{\mathbb{R}^n}\int_{\mathbb{R}^n}\Gamma_t(x-y)|f(y)|\,d\mathcal{L}^n(x)\,d\mathcal{L}^n(y)=\int_{\mathbb{R}^n}|f(y)|\,d\mathcal{L}^n(y).
\end{align*}
The right-hand side is finite because $f\in L^1(\mathbb{R}^n)$. Hence the signed integrand
\begin{align*}
G_t: \mathbb{R}^n \times \mathbb{R}^n \to \mathbb{R}, \qquad (x,y) \mapsto \Gamma_t(x-y)f(y)
\end{align*}
satisfies $|G_t|=H_t$ and belongs to $L^1(\mathbb{R}^n\times\mathbb{R}^n)$.
This product-integrability has two consequences. First, for $\mathcal{L}^n$-a.e. $x\in\mathbb{R}^n$, the section
\begin{align*}
y \mapsto \Gamma_t(x-y)f(y)
\end{align*}
is absolutely integrable over $\mathbb{R}^n$, so the integral defining $u_t(x)$ is absolutely convergent for such $x$. Second, the pointwise triangle inequality gives
\begin{align*}
|u_t(x)|\leq \int_{\mathbb{R}^n}\Gamma_t(x-y)|f(y)|\,d\mathcal{L}^n(y)
\end{align*}
for every $x$ for which the defining integral is absolutely convergent. Integrating this inequality in $x$ and applying Tonelli to the non-negative right-hand side gives
\begin{align*}
\int_{\mathbb{R}^n}|u_t(x)|\,d\mathcal{L}^n(x)\leq \int_{\mathbb{R}^n}\int_{\mathbb{R}^n}\Gamma_t(x-y)|f(y)|\,d\mathcal{L}^n(y)\,d\mathcal{L}^n(x).
\end{align*}
The same Tonelli and translation computation evaluates the right-hand side as
\begin{align*}
\int_{\mathbb{R}^n}|f(y)|\,d\mathcal{L}^n(y)=\|f\|_{L^1(\mathbb{R}^n)}.
\end{align*}
Thus $\int_{\mathbb{R}^n}|u_t(x)|\,d\mathcal{L}^n(x)<\infty$, which proves $u_t\in L^1(\mathbb{R}^n)$.
[/guided]
[/step]
[step:Interchange the order of integration and collapse the heat kernel]
Since $G_t\in L^1(\mathbb{R}^n\times\mathbb{R}^n)$, Fubini's theorem applies to $G_t$. Therefore
\begin{align*}
\int_{\mathbb{R}^n}u_t(x)\,d\mathcal{L}^n(x)=\int_{\mathbb{R}^n}\int_{\mathbb{R}^n}\Gamma_t(x-y)f(y)\,d\mathcal{L}^n(y)\,d\mathcal{L}^n(x).
\end{align*}
Interchanging the order of integration gives
\begin{align*}
\int_{\mathbb{R}^n}u_t(x)\,d\mathcal{L}^n(x)=\int_{\mathbb{R}^n}f(y)\left(\int_{\mathbb{R}^n}\Gamma_t(x-y)\,d\mathcal{L}^n(x)\right)d\mathcal{L}^n(y).
\end{align*}
For each fixed $y\in\mathbb{R}^n$, use the translation change of variables $z=x-y$. This preserves Lebesgue measure, so $d\mathcal{L}^n(z)=d\mathcal{L}^n(x)$, and the domain $\mathbb{R}^n$ remains $\mathbb{R}^n$. Hence
\begin{align*}
\int_{\mathbb{R}^n}\Gamma_t(x-y)\,d\mathcal{L}^n(x)=\int_{\mathbb{R}^n}\Gamma_t(z)\,d\mathcal{L}^n(z)=1.
\end{align*}
Substitution into the previous identity yields
\begin{align*}
\int_{\mathbb{R}^n}u_t(x)\,d\mathcal{L}^n(x)=\int_{\mathbb{R}^n}f(y)\,d\mathcal{L}^n(y).
\end{align*}
[guided]
Now we compute the total mass of $u_t$. The only delicate point is that we want to interchange the order of integration, so we must first verify absolute integrability. That was done in the previous step: the function
\begin{align*}
G_t: \mathbb{R}^n \times \mathbb{R}^n \to \mathbb{R}, \qquad (x,y) \mapsto \Gamma_t(x-y)f(y)
\end{align*}
belongs to $L^1(\mathbb{R}^n\times\mathbb{R}^n)$. Therefore Fubini's theorem applies to $G_t$.
Using the definition of $u_t$, we have
\begin{align*}
\int_{\mathbb{R}^n}u_t(x)\,d\mathcal{L}^n(x)=\int_{\mathbb{R}^n}\int_{\mathbb{R}^n}\Gamma_t(x-y)f(y)\,d\mathcal{L}^n(y)\,d\mathcal{L}^n(x).
\end{align*}
Fubini's theorem allows us to integrate first in $x$ and then in $y$:
\begin{align*}
\int_{\mathbb{R}^n}u_t(x)\,d\mathcal{L}^n(x)=\int_{\mathbb{R}^n}f(y)\left(\int_{\mathbb{R}^n}\Gamma_t(x-y)\,d\mathcal{L}^n(x)\right)d\mathcal{L}^n(y).
\end{align*}
The inner integral is the total mass of a translated heat kernel. For a fixed $y\in\mathbb{R}^n$, define the change of variables $z=x-y$. Translation preserves Lebesgue measure, so $d\mathcal{L}^n(z)=d\mathcal{L}^n(x)$, and as $x$ ranges over $\mathbb{R}^n$, so does $z$. Hence
\begin{align*}
\int_{\mathbb{R}^n}\Gamma_t(x-y)\,d\mathcal{L}^n(x)=\int_{\mathbb{R}^n}\Gamma_t(z)\,d\mathcal{L}^n(z).
\end{align*}
By the normalization of the heat kernel, this last integral is $1$. Therefore
\begin{align*}
\int_{\mathbb{R}^n}u_t(x)\,d\mathcal{L}^n(x)=\int_{\mathbb{R}^n}f(y)\,d\mathcal{L}^n(y).
\end{align*}
This is exactly the desired conservation of total mass at the fixed time $t>0$.
[/guided]
[/step]
[step:Return from the fixed time notation to the theorem statement]
Since $t>0$ was arbitrary and $u_t=u(t,\cdot)=e^{t\Delta}f$, the preceding steps prove that $u(t,\cdot)\in L^1(\mathbb{R}^n)$ and
\begin{align*}
\int_{\mathbb{R}^n}u(t,x)\,d\mathcal{L}^n(x)=\int_{\mathbb{R}^n}f(x)\,d\mathcal{L}^n(x)
\end{align*}
for every $t>0$.
[/step]
