Consider the Cauchy problem for the heat equation:
\begin{align*}
\left\{ \begin{array}{r l}
\frac{\partial u}{\partial t} - \Delta u = 0 & \text{in} \, \, \mathbb{R}^{d} \times (0,\infty) \\
u(x,0) = u_{0} & \text{on} \, \, \mathbb{R}^{d} \end{array} \right.
\end{align*}
where $u_{0} : \mathbb{R}^{d} \to \mathbb{R}$ is a "nice" function. (A convenient choice is $u_{0} \in \mathscr{S}(\mathbb{R}^{d})$.)
The beauty of the Fourier transform in this context is $\Delta$ becomes an algebraic "multiplier." Indeed, if I multiply all sides by $e^{-i 2 \pi \langle \xi, x \rangle}$ and integrate with respect to $x$, then I obtain
\begin{align*}
\int_{\mathbb{R}^{d}} \frac{\partial u}{\partial t}(x,t) e^{-i 2 \pi \langle \xi, x \rangle} \, dx &= \frac{\partial}{\partial t} \left(\int_{\mathbb{R}^{d}} u(x,t) e^{- i 2 \pi \langle \xi, x \rangle} \, dx \right) = \frac{\partial \hat{u}}{\partial t}(\xi,t) \\
\int_{\mathbb{R}^{d}} \frac{\partial^{2} u}{\partial x_{j}^{2}}(x,t) e^{-i 2 \pi \langle \xi, x \rangle} \, dx &= (- i2 \pi \xi_{j})^{2} \int_{\mathbb{R}^{d}} u(x,t) e^{-i 2 \pi \langle \xi, x \rangle} \, dx = - 4 \pi^{2} \xi_{j}^{2} \hat{u}(\xi,t),
\end{align*}
where in the second line I integrated by parts twice and assumed that $\lim_{|x| \to \infty} u(x,t) = 0$ uniformly (for each fixed $t$). (Note: Here $\langle \xi, x \rangle = \sum_{j = 1}^{d} \xi_{j} x_{j}$, that is, it's a fancy way of writing the dot product.) Thus, summing over $j$, we obtain
\begin{align*}
\left\{ \begin{array}{r l}
\frac{\partial u}{\partial t}(\xi,t) + 4 \pi^{2} |\xi|^{2} \hat{u}(\xi,t) = 0 & \text{in} \, \, \mathbb{R}^{d} \times (0,\infty) \\
\hat{u}(\xi,0) = \hat{u}_{0}(\xi) & \text{on} \, \, \mathbb{R}^{d} \end{array} \right.
\end{align*}
For each fixed $\xi \in \mathbb{R}^{d}$, this is an ODE. In particular,
$$\hat{u}(\xi,t) = \hat{u}_{0}(\xi) e^{- 4 \pi^{2} |\xi|^{2} t}.$$
Therefore, by Fourier inversion,
\begin{align*}
u(x,t) = \int_{\mathbb{R}^{d}} \hat{u}_{0}(\xi) e^{-4 \pi^{2} |\xi|^{2} t + i 2 \pi \langle \xi, x \rangle} \, d\xi
\end{align*}
This isn't entirely what we want: we would like to express the answer in terms of $u_{0}$, not $\hat{u}_{0}$.
Recall the following fact about the Fourier transform: if $f,g$ are "nice" functions (e.g. $f,g \in \mathscr{S}(\mathbb{R}^{d})$), then
$$(f * g)(x) = \int_{\mathbb{R}^{d}} \hat{f}(\xi) \hat{g}(\xi) e^{i 2 \pi \langle \xi, x \rangle} \, d\xi.$$
In other words, "convolution in space corresponds to multiplication in frequency." Thus, let $G_{t}$ denote the function satisfying
$$\hat{G}_{t}(\xi) = e^{-4 \pi^{2} |\xi|^{2} t}.$$ By the last remark, we are led to guess that
$$\int_{\mathbb{R}^{d}} u_{0}(y) G_{t}(x - y) \, dy = \int_{\mathbb{R}^{d}} \hat{u}_{0}(\xi) \hat{G}_{t}(\xi) e^{i 2 \pi \langle \xi, x \rangle} \, d\xi.$$
This leads naturally to the question: what is $G_{t}$? By Fourier inversion,
$$G_{t}(x) = \int_{\mathbb{R}^{d}} e^{- 4 \pi^{2} |\xi|^{2} t + i 2 \pi \langle \xi ,x \rangle} \, d\xi.$$
The RHS can be evaluated using complex analysis (or found in a Fourier transform table). The results is:
$$G_{t}(x) = (4 \pi t)^{-\frac{d}{2}} e^{- \frac{|x|^{2}}{4t}}.$$
Since $G_{t}$ is a very nice function (i.e. $G_{t} \in \mathscr{S}(\mathbb{R}^{d})$), our guess is justified and we obtain
$$u(x,t) = (4 \pi t)^{-\frac{d}{2}} \int_{\mathbb{R}^{d}} u_{0}(y) e^{- \frac{|y - x|^{2}}{4t}} \, dy.$$
In the case when $u_{0} = \delta_{0}$, we find $u(x,t) = G_{t}(x)$. Therefore, we showed $\tilde{G}(x,t) = G_{t}(x)$ satisfies
\begin{align*}
\left\{\begin{array}{r l}
\frac{\partial \tilde{G}}{\partial t} - \Delta \tilde{G} = 0 & \text{in} \, \, \mathbb{R}^{d} \times (0,\infty) \\
\tilde{G}(dx,0) = \delta_{0}(dx) & \text{on} \, \, \mathbb{R}^{d}.
\end{array} \right.
\end{align*}
The function $\tilde{G}$ is called the heat kernel.
A classical question is now: how can we solve inhomogeneous problem, i.e.
\begin{align*}
\left\{ \begin{array}{r l}
\frac{\partial v}{\partial t} - \Delta v = f(x,t) & \text{in} \, \, \mathbb{R}^{d} \times (0,\infty) \\
v(x,0) = u_{0} & \text{on} \, \, \mathbb{R}^{d} \end{array} \right.
\end{align*}
By linearity, it suffices to consider the case when $u_{0} = 0$. (Otherwise, add the solution with $0$ initial data to the solution of the homogeneous equation with initial data $u_{0}$.) If we take the Fourier transform in space again, we find
\begin{align*}
\left\{ \begin{array}{r l}
\frac{\partial \hat{v}}{\partial t}(\xi,t) + 4 \pi^{2} |\xi|^{2} \hat{v}(\xi,t) = \hat{f}(\xi,t) & \text{in} \, \, \mathbb{R}^{d} \times (0,\infty) \\
\hat{v}(\xi,0) = 0 & \text{on} \, \, \mathbb{R}^{d}.
\end{array} \right.
\end{align*}
Recall that this linear ODE can be solved using Duhamel's formula:
\begin{align*}
\hat{v}(\xi,t) = \int_{0}^{t} \hat{f}(\xi,s) e^{- 4 \pi^{2} |\xi|^{2} (t - s)} \, ds.
\end{align*}
If we apply the inverse Fourier transform and Fubini's Theorem, we find
\begin{align*}
v(x,t) &= \int_{0}^{t} \left( \int_{\mathbb{R}^{d}} \hat{f}(\xi,s) e^{- 4 \pi^{2} |\xi|^{2} (t - s) + i 2 \pi \langle \xi, x \rangle} \, d \xi \right) \, ds \\
&= \int_{0}^{t} \left( \int_{\mathbb{R}^{d}} \hat{f}(\xi,s) \hat{G}_{t - s}(\xi) e^{i 2 \pi \langle \xi, x \rangle} \, d \xi \right) \, ds \\
&= \int_{0}^{t} \left(\int_{\mathbb{R}^{d}} f(y,s) G_{t - s}(x - y) \, dy \right) \, ds.
\end{align*}
In the case when $f(dx,dt) = \delta_{0}(dx) \delta_{0}(dt)$, the solution is
$$v(x,t) = G_{t}(x),$$
which shows that the heat kernel is Green's function for the heat equation.
