We give two proofs. The first is a direct classical verification using the Leibniz integral rule. The second uses the space-time distributional convolution $E * \tilde{f}$, where $E$ is the causal extension of $\Phi$ satisfying $(\partial_t - \Delta)E = \delta$ in $\mathcal{D}'(\mathbb{R}^{n+1})$.

---

## Proof 1: Classical verification

**Proof plan.** For each fixed $s \in [0,t)$, define $v^s(x,t) := \int_{\mathbb{R}^n}\Phi(x-y,t-s)\,f(y,s)\,d\mathcal{L}^n(y)$, which solves the homogeneous heat equation for $t > s$ with initial data $f(\cdot, s)$ at time $s$ by the [solution of the Cauchy problem](/theorems/54). The candidate is $u(x,t) = \int_0^t v^s(x,t)\,d\mathcal{L}^1(s)$. We differentiate in $t$ using the Leibniz rule (Claim 1): the boundary term at $s = t$ contributes $v^t(x,t) = f(x,t)$ by the delta-limit property (Claim 2), while the integral term contributes $\int_0^t \partial_t v^s\,d\mathcal{L}^1(s) = \int_0^t \Delta_x v^s\,d\mathcal{L}^1(s) = \Delta u$ (Claim 3). Combining gives $\partial_t u = f + \Delta u$.

**Step 1 (The auxiliary functions).**

For each $s \geq 0$, define the auxiliary function

\begin{align*} v^s: \mathbb{R}^n \times (s, \infty) &\to \mathbb{R} \\ (x,t) &\mapsto \int_{\mathbb{R}^n}\Phi(x - y, t - s)\,f(y,s)\,d\mathcal{L}^n(y). \end{align*}

Since $f(\cdot, s) \in C^2(\mathbb{R}^n) \cap L^\infty(\mathbb{R}^n)$ (bounded because $f$ has compact support), the [solution of the Cauchy problem](/theorems/54) applied with initial data $g := f(\cdot, s)$ and time origin shifted to $s$ gives: $v^s \in C^\infty(\mathbb{R}^n \times (s,\infty))$, $v^s$ satisfies $\partial_t v^s - \Delta_x v^s = 0$ for $t > s$, and

\begin{align*} \lim_{\substack{(x,t) \to (x_0, s) \\ t > s}} v^s(x,t) = f(x_0, s) \quad \text{for every } x_0 \in \mathbb{R}^n. \end{align*}

The candidate solution is then

\begin{align*} u(x,t) = \int_0^t v^s(x,t)\,d\mathcal{L}^1(s). \end{align*}

**Step 2 (Differentiating in $t$: the Leibniz rule).**

[claim:Leibniz Rule For The Duhamel Integral] For all $(x,t) \in \mathbb{R}^n \times (0,\infty)$, \begin{align*} \partial_t u(x,t) = v^t(x,t) + \int_0^t \partial_t v^s(x,t)\,d\mathcal{L}^1(s). \end{align*} [/claim]

[proof] The function $u(x,t) = \int_0^t v^s(x,t)\,d\mathcal{L}^1(s)$ has both the integrand and the upper limit depending on $t$. The Leibniz integral rule states that if $F(t) = \int_0^{a(t)} h(s,t)\,d\mathcal{L}^1(s)$ with $a(t) = t$, and if $h$ and $\partial_t h$ are continuous on the region $\{(s,t) : 0 \leq s \leq t\}$, then

\begin{align*} F'(t) = h(t,t)\cdot a'(t) + \int_0^{a(t)} \partial_t h(s,t)\,d\mathcal{L}^1(s) = h(t,t) + \int_0^t \partial_t h(s,t)\,d\mathcal{L}^1(s). \end{align*}

We verify the hypotheses. Set $h(s,t) := v^s(x,t)$ for fixed $x$. For $0 \leq s < t$, the function $(s,t) \mapsto v^s(x,t)$ is continuous: continuity in $t$ follows from $v^s \in C^\infty$ for $t > s$, and continuity in $s$ follows from the smoothness of $f$ in its time argument (the integral $\int_{\mathbb{R}^n}\Phi(x-y,t-s)\,f(y,s)\,d\mathcal{L}^n(y)$ depends continuously on $s$ by dominated convergence, since $f$ has compact support in the spatial variable and varies continuously in $s$). Similarly, $\partial_t v^s(x,t) = \int_{\mathbb{R}^n}\partial_t\Phi(x-y,t-s)\,f(y,s)\,d\mathcal{L}^n(y)$ is continuous in $(s,t)$ for $0 \leq s < t$ (the differentiation under the integral sign is justified exactly as in the proof of the [Cauchy problem](/theorems/54), Claim 1, using the Gaussian decay of $\partial_t\Phi$ and boundedness of $f$).

The remaining point is the behaviour as $s \to t$. For $s$ near $t$, the kernel $\Phi(x-y,t-s)$ concentrates near $y = x$ as $t - s \to 0$, and $v^s(x,t) \to f(x,s)$ by the delta-limit property. Since $f$ is continuous in $s$, the function $s \mapsto v^s(x,t)$ extends continuously to $s = t$ with value $f(x,t)$. This gives $h(t,t) = f(x,t)$, and the Leibniz rule applies. [/proof]

**Step 3 (The boundary term).**

[claim:Boundary Term Evaluation] For all $(x,t) \in \mathbb{R}^n \times (0,\infty)$, the boundary term satisfies $v^t(x,t) = f(x,t)$. [/claim]

[proof] This is a direct application of the delta-limit property (property 4 of the [fundamental solution](/theorems/53)). Formally, the value $v^t(x,t)$ is the limit

\begin{align*} v^t(x,t) = \lim_{s \to t^-} v^s(x,t) = \lim_{s \to t^-} \int_{\mathbb{R}^n}\Phi(x - y, t - s)\,f(y,s)\,d\mathcal{L}^n(y). \end{align*}

As $s \to t^-$, the time parameter $\tau := t - s \to 0^+$, and $f(y,s) \to f(y,t)$ uniformly in $y$ (since $f$ is continuous in $s$ with compact spatial support). Therefore

\begin{align*} v^t(x,t) &= \lim_{\tau \to 0^+}\int_{\mathbb{R}^n}\Phi(x - y, \tau)\,f(y, t - \tau)\,d\mathcal{L}^n(y). \end{align*}

Write the integrand as $\Phi(x-y,\tau)\,f(y,t) + \Phi(x-y,\tau)\,(f(y,t-\tau) - f(y,t))$. The first term converges to $f(x,t)$ by the delta-limit property applied to the continuous bounded function $y \mapsto f(y,t)$. The second term is bounded in absolute value by $\Phi(x-y,\tau)\cdot\sup_{y}|f(y,t-\tau) - f(y,t)|$. Since $f$ has compact spatial support and is continuous in $t$, the supremum tends to $0$ as $\tau \to 0$, while $\int\Phi\,d\mathcal{L}^n = 1$. Therefore the second term vanishes and $v^t(x,t) = f(x,t)$. [/proof]

**Step 4 (The integral term equals $\Delta u$).**

[claim:Interchange Of Laplacian And Time Integral] For all $(x,t) \in \mathbb{R}^n \times (0,\infty)$, \begin{align*} \int_0^t \partial_t v^s(x,t)\,d\mathcal{L}^1(s) = \Delta u(x,t). \end{align*} [/claim]

[proof] Since $v^s$ solves the homogeneous heat equation for $t > s$ (Step 1), we have $\partial_t v^s(x,t) = \Delta_x v^s(x,t)$ for each fixed $s \in [0,t)$. Therefore

\begin{align*} \int_0^t \partial_t v^s(x,t)\,d\mathcal{L}^1(s) = \int_0^t \Delta_x v^s(x,t)\,d\mathcal{L}^1(s). \end{align*}

It remains to show that $\Delta_x$ passes through the $s$-integral. For each $s \in [0,t)$ and each coordinate index $i$, the second spatial derivative is

\begin{align*} \partial_{x_i x_i} v^s(x,t) = \int_{\mathbb{R}^n}\partial_{x_i x_i}\Phi(x-y,t-s)\,f(y,s)\,d\mathcal{L}^n(y). \end{align*}

This is justified by the dominating function argument from the [Cauchy problem](/theorems/54), Claim 1: $|\partial_{x_i x_i}\Phi(x-y,\tau)| \leq C\,\tau^{-n/2-1}\exp(-|x-y|^2/(8\tau))$ for $\tau = t - s > 0$, and $f$ is bounded with compact spatial support.

To interchange $\int_0^t$ and $\partial_{x_i x_i}$, we need the joint integrability of $\partial_{x_i x_i} v^s(x,t)$ in $s$ near $s = t$, where $\tau = t - s \to 0$. Using the bound above:

\begin{align*} |\partial_{x_i x_i} v^s(x,t)| &\leq C\,\|f\|_{L^\infty}\,(t-s)^{-n/2-1}\int_{\mathbb{R}^n}\exp\!\left(-\frac{|x-y|^2}{8(t-s)}\right)d\mathcal{L}^n(y) = C'\,(t-s)^{-1}, \end{align*}

where the Gaussian integral contributes $(t-s)^{n/2}$, cancelling the $(t-s)^{-n/2}$ factor and leaving $(t-s)^{-1}$. This is integrable on $[0, t-\delta]$ for any $\delta > 0$ but only barely fails at $s = t$. However, the actual integrand $\partial_{x_i x_i}v^s(x,t)$ is better behaved because $f$ has compact spatial support: for $s$ close to $t$, the function $f(\cdot, s)$ is concentrated in a compact set, and $\Phi(\cdot, t-s)$ concentrates near the origin, so the convolution $v^s(x,t)$ and its derivatives remain bounded as $s \to t$ (they converge to $\partial_{x_i x_i}f(x,t)$, which exists since $f \in C^{2,1}$). More precisely, writing $\partial_{x_i x_i}v^s(x,t) = \int \Phi(x-y,t-s)\,\partial_{y_i y_i}f(y,s)\,d\mathcal{L}^n(y)$ (integrating by parts twice, with boundary terms vanishing since $f$ has compact support) shows that $|\partial_{x_i x_i}v^s(x,t)| \leq \|\partial_{x_i x_i}f\|_{L^\infty}$ uniformly in $s$. This uniform bound makes $s \mapsto \partial_{x_i x_i}v^s(x,t)$ bounded on $[0,t]$, and Fubini's theorem (or differentiation under the integral sign via dominated convergence) gives

\begin{align*} \int_0^t \Delta_x v^s(x,t)\,d\mathcal{L}^1(s) = \Delta_x \int_0^t v^s(x,t)\,d\mathcal{L}^1(s) = \Delta u(x,t). \end{align*} [/proof]

**Step 5 (Conclusion).**

Combining Claims 1–3:

\begin{align*} \partial_t u(x,t) &= v^t(x,t) + \int_0^t \partial_t v^s(x,t)\,d\mathcal{L}^1(s) = f(x,t) + \Delta u(x,t). \end{align*}

Rearranging gives $\partial_t u - \Delta u = f$ in $\mathbb{R}^n \times (0,\infty)$.

For the initial condition, at $t = 0$ the domain of integration $[0,0]$ has measure zero:

\begin{align*} u(x,0) = \int_0^0 v^s(x,0)\,d\mathcal{L}^1(s) = 0. \end{align*}

---

## Proof 2: Distributional space-time convolution

**Proof plan.** Extend $\Phi$ to all of $\mathbb{R}^{n+1}$ by setting it to zero for $t \leq 0$, obtaining a distribution $E \in \mathcal{D}'(\mathbb{R}^{n+1})$ satisfying $(\partial_t - \Delta)E = \delta$ (Claim 1). Extend $f$ similarly to obtain $\tilde{f}$ with compact support in $\mathbb{R}^{n+1}$. The Duhamel formula is exactly the space-time convolution $u = E * \tilde{f}$ (Claim 2). Since $\tilde{f}$ has compact support, the heat operator commutes with convolution: $(\partial_t - \Delta)(E * \tilde{f}) = ((\partial_t - \Delta)E) * \tilde{f} = \delta * \tilde{f} = \tilde{f}$ (Claim 3).

**Step 1 (The causal fundamental solution).**

[claim:Causal Fundamental Solution Satisfies Delta Equation] Define $E: \mathbb{R}^{n+1} \to \mathbb{R}$ by \begin{align*}

E(x,t) := \begin{cases} \Phi(x,t) & \text{if } t > 0, \\ 0 & \text{if } t \leq 0. \end{cases} \end{align*} Then $E \in L^1_{\mathrm{loc}}(\mathbb{R}^{n+1})$ and, as an element of $\mathcal{D}'(\mathbb{R}^{n+1})$, it satisfies $(\partial_t - \Delta)E = \delta$ where $\delta$ is the Dirac distribution at the origin in $\mathbb{R}^{n+1}$. [/claim] [proof] First, $E \in L^1_{\mathrm{loc}}(\mathbb{R}^{n+1})$: for any compact set $K \subset \mathbb{R}^{n+1}$, the integral $\int_K |E(x,t)|\,d\mathcal{L}^{n+1}(x,t)$ is finite because $\Phi(x,t) \leq C\,t^{-n/2}$ and $t^{-n/2}$ is locally integrable in $t$ near $t = 0$ when integrated against the Gaussian spatial factor (in fact, $\int_{\mathbb{R}^n}\Phi(x,t)\,d\mathcal{L}^n(x) = 1$ for all $t > 0$, so the spatial integral is uniformly bounded). To verify $(\partial_t - \Delta)E = \delta$, let $\psi \in C_c^\infty(\mathbb{R}^{n+1})$ be an arbitrary test function. We must show \begin{align*} \langle (\partial_t - \Delta)E, \psi \rangle = \psi(0,0). \end{align*} By definition of distributional derivatives, $\langle (\partial_t - \Delta)E, \psi \rangle = -\langle E, \partial_t \psi \rangle - \langle E, \Delta \psi \rangle$, but since $E$ is a regular distribution and $\partial_t$ acts on $E$ distributionally, we write \begin{align*} \langle (\partial_t - \Delta)E, \psi \rangle &= -\int_0^\infty\int_{\mathbb{R}^n}\Phi(x,t)\,\partial_t\psi(x,t)\,d\mathcal{L}^n(x)\,d\mathcal{L}^1(t) - \int_0^\infty\int_{\mathbb{R}^n}\Phi(x,t)\,\Delta\psi(x,t)\,d\mathcal{L}^n(x)\,d\mathcal{L}^1(t). \end{align*} For the second integral, since $\Phi(\cdot,t)$ is smooth and decays as a Gaussian for each $t > 0$, integration by parts in $x$ (with all boundary terms vanishing) gives $\int_{\mathbb{R}^n}\Phi\,\Delta\psi\,d\mathcal{L}^n = \int_{\mathbb{R}^n}(\Delta\Phi)\,\psi\,d\mathcal{L}^n$. For the first integral, integrate by parts in $t$ on $[\varepsilon, T]$ (where $T$ is large enough that $\psi$ vanishes for $t > T$): \begin{align*} -\int_\varepsilon^T\int_{\mathbb{R}^n}\Phi\,\partial_t\psi\,d\mathcal{L}^n\,d\mathcal{L}^1 &= \int_\varepsilon^T\int_{\mathbb{R}^n}(\partial_t\Phi)\,\psi\,d\mathcal{L}^n\,d\mathcal{L}^1 - \left[\int_{\mathbb{R}^n}\Phi(x,t)\,\psi(x,t)\,d\mathcal{L}^n(x)\right]_{t=\varepsilon}^{t=T}. \end{align*} The boundary term at $t = T$ vanishes because $\psi(\cdot, T) = 0$. The boundary term at $t = \varepsilon$ is $-\int_{\mathbb{R}^n}\Phi(x,\varepsilon)\,\psi(x,\varepsilon)\,d\mathcal{L}^n(x)$. Combining and using $\partial_t\Phi - \Delta\Phi = 0$ for $t > 0$ (property 2 of the [fundamental solution](/theorems/53)): \begin{align*} \langle (\partial_t - \Delta)E, \psi \rangle &= \lim_{\varepsilon \to 0^+}\left(\int_\varepsilon^T\int_{\mathbb{R}^n}\underbrace{(\partial_t\Phi - \Delta\Phi)}_{= 0}\,\psi\,d\mathcal{L}^n\,d\mathcal{L}^1 + \int_{\mathbb{R}^n}\Phi(x,\varepsilon)\,\psi(x,\varepsilon)\,d\mathcal{L}^n(x)\right). \end{align*} The first integral vanishes. For the second, since $\psi$ is continuous and bounded, and $\psi(x,\varepsilon) \to \psi(x,0)$ uniformly in $x$ as $\varepsilon \to 0$, the delta-limit property (property 4 of the [fundamental solution](/theorems/53)) gives \begin{align*} \lim_{\varepsilon \to 0^+}\int_{\mathbb{R}^n}\Phi(x,\varepsilon)\,\psi(x,\varepsilon)\,d\mathcal{L}^n(x) = \psi(0,0). \end{align*} Therefore $\langle (\partial_t - \Delta)E, \psi \rangle = \psi(0,0) = \langle \delta, \psi \rangle$. [/proof] **Step 2 (The Duhamel formula as space-time convolution).** [claim:Duhamel Formula Equals Space Time Convolution] Define $\tilde{f}: \mathbb{R}^{n+1} \to \mathbb{R}$ by $\tilde{f}(x,t) := f(x,t)$ for $t \geq 0$ and $\tilde{f}(x,t) := 0$ for $t < 0$. Then $\tilde{f}$ has compact support in $\mathbb{R}^{n+1}$ and \begin{align*} u(x,t) = (E * \tilde{f})(x,t) \end{align*} for all $(x,t) \in \mathbb{R}^n \times [0,\infty)$. [/claim] [proof] Since $f$ has compact support in $\mathbb{R}^n \times [0,\infty)$, $\tilde{f}$ has compact support in $\mathbb{R}^{n+1}$. The space-time convolution is \begin{align*} (E * \tilde{f})(x,t) &= \int_{-\infty}^{\infty}\int_{\mathbb{R}^n}E(x - y, t - s)\,\tilde{f}(y,s)\,d\mathcal{L}^n(y)\,d\mathcal{L}^1(s). \end{align*} The integrand is nonzero only when both $E(x-y,t-s) \neq 0$ and $\tilde{f}(y,s) \neq 0$. The first condition requires $t - s > 0$, i.e. $s < t$. The second requires $s \geq 0$. Therefore the $s$-integral reduces to the interval $[0,t]$, and $E(x-y,t-s) = \Phi(x-y,t-s)$, $\tilde{f}(y,s) = f(y,s)$: \begin{align*} (E * \tilde{f})(x,t) = \int_0^t\int_{\mathbb{R}^n}\Phi(x - y, t - s)\,f(y,s)\,d\mathcal{L}^n(y)\,d\mathcal{L}^1(s) = u(x,t). \end{align*} [/proof] **Step 3 (Applying the heat operator to the convolution).** [claim:Heat Operator Passes Through Convolution] $(\partial_t - \Delta)u = f$ in $\mathbb{R}^n \times (0,\infty)$ and $u(\cdot, 0) = 0$. [/claim] [proof] Since $\tilde{f}$ has compact support and $E \in L^1_{\mathrm{loc}}(\mathbb{R}^{n+1})$, the convolution $E * \tilde{f}$ is well-defined and yields a $C^{2,1}$ function (this follows from the smoothing property of convolution: the regularity of $E$ for $t > 0$ combined with the compact support and $C^{2,1}$ regularity of $\tilde{f}$ ensures that $u$ inherits $C^{2,1}$ regularity). For distributional derivatives of convolutions, the general identity states: if $T \in \mathcal{D}'(\mathbb{R}^{n+1})$ and $\varphi$ has compact support, then $D^\alpha(T * \varphi) = (D^\alpha T) * \varphi$ for any partial derivative $D^\alpha$. Applying this with $T = E$ and $\varphi = \tilde{f}$, and using the linearity of $\partial_t - \Delta$: \begin{align*} (\partial_t - \Delta)(E * \tilde{f}) = ((\partial_t - \Delta)E) * \tilde{f} = \delta * \tilde{f}, \end{align*} where the last equality uses Claim 1. Since convolution with $\delta$ is the identity ($\delta * \tilde{f} = \tilde{f}$): \begin{align*} (\partial_t - \Delta)u = \tilde{f}. \end{align*} For $t > 0$, $\tilde{f}(x,t) = f(x,t)$, giving $\partial_t u - \Delta u = f$ in $\mathbb{R}^n \times (0,\infty)$. For the initial condition: when $t \leq 0$, the convolution $(E * \tilde{f})(x,t)$ involves integrating over $s \in (-\infty, t] \cap [0,\infty)$. For $t \leq 0$, this set is at most $\{0\}$, which has measure zero, so $u(x,t) = 0$ for $t \leq 0$. In particular, $u(\cdot, 0) = 0$. [/proof]