Solution 1:

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.