Consider a symmetric random walk with density $p_t(x)$ with time $t$ in discrete steps. Then the probability for the next step in $t$ is \[ p_{t + \Delta t}(x) = \frac{1}{2} p_t(x + \Delta x) + \frac{1}{2} p_t(x - \Delta x) \] \[ \partial_t u_t(x) = \Delta u_t(x) \] If $\frac{(\Delta x)^2}{\Delta t} = 1$, the random walk approximately solves the heat equation. With initial value $p_0$ , we have that \[ \frac{p_{t + \Delta t}(x) - p_{t + \Delta t}(x)}{\Delta t} = \frac{\frac{1}{2} p_t(x + \Delta x) - \frac{1}{2} p_t(x) + \frac{1}{2} p_t(x - \Delta x) - \frac{1}{2} p_t(x)}{(\Delta x)^2} \] is approximately \[ \partial_t p_t(x) = \frac{1}{2} \Delta p_t(x). \] Now, we seek to solve the heat equation by taking a limit of a properly scaled random walk and this will be called Brownian motion.
Interpolate a 1-dimensional random walk with time steps $\frac{1}{n}$ and space steps $\frac{1}{\sqrt{n}}$. This is \[ S_{nt} = \sum_{k = 1}^{\lfloor nt \rfloor} \frac{1}{\sqrt{n}} Z_k + (nt - \lfloor nt \rfloor) \frac{1}{\sqrt{n}} Z_{\lfloor nt \rfloor + 1} \] where $Z_k$ are independent standard normal random variables. Now, $S_{nt}$ is normal with $Var (S_{nt}) = t$ and $E S_{nt} = 0$ and the finite dimensional distributions $(S_{nt_1}, \ldots, S_{nt_k})$ are normal with covariance $E S_{nt_i} S_{nt_j} = t_i \wedge t_j$. Because this is normal the tails are controlled with \[ E|S_{nt} - S_{ns}|^4 = 3 (E|S_{nt} - S_{ns}|^2)^2 = 3 |t - s|^2. \] Thus, a Theorem by Kolmogorov shows this has a continuous limit with the same finite dimensional distributions. This is Brownian motion. Everything can be extended from the 1-dimensional case, so we only need to prove the 1-dimensional case.
Definition (Brownian Motion). A d-dimensional Brownian Motion is a process $W_t$ on a probability space $(\Omega, \F, P)$ over [0, T] where $t \mapsto W_t$ is continuous and $W_{t_1}, \ldots, W_{t_k}$ is Normally distributed with $EW_{t_i} W_{t_j} = t_i \wedge t_j \delta_{i j}$.
Thereom (Existence of Brownian Motion). $W_t$ is a limit of $S_{nt}$.
Proof. By Chebyshev, \[ P(|S_{nt} - S_{ns}| \ge \epsilon) \le E|S_{nt} - S_{ns}|^{4} \epsilon^{-4} = 3|t - s|^2 \epsilon^{-4}. \] Since $[0, T]$ is compact, we can partition it with $t_1 < \cdots < t_{2^n} = T$. Then \[ P(|S_{nt} - S_{ns}| \ge \epsilon) = P(\max_{0 \le k \le 2^n} \max_{\frac{T k}{2^n} \le t, s \le \frac{T(k + 1)}{2^n}} |S_{nt} - S_{ns}| \ge \epsilon) \le 2^{n + 1} 3 T 2^{-2n} \epsilon^{-4} \] Now, choose $\epsilon = 2^{-\frac{n}{8}}$ and \[ P(\max_{0 \le s, t \le T} |S_{nt} - S_{ns}| \ge \epsilon) \le 6 T 2^{-\frac{n}{2}}. \] By Borel-Cantelli, we have that \[ \max_{0 \le s, t \le T} |S_{nt} - S_{ns}| \] is Cauchy a.e.. Since $t \mapsto S_{nt}$ is continuous by construction and $C[0, T]$ is complete, there exists a continuous a continuous limit $S^*$ on $[0, T] $a.e. By continuity, the finite dimensional distributions must coincide. Define the measure $PS^*$ as the measure for the Brownian motion process $(W_t)_{t \in [0, T]}$. $\blacksquare$
Finally, let's solve the heat equation. Now, $W_t$ is Markov (Easy to show). Let $f$ be $C^2$. Because it is Markov, we have \[ E_x f(W_{t + h}) = E_x E_{W_t} f(W_h) = E_x f(W_{t}) + \frac{h}{2} \Delta E_x f(W_t) + C. \] Therefore, with initial value $f(x)$, we can solve the heat equation \[ \partial_t E_x f(W_t) = \frac{1}{2} \Delta E_x f(W_t) \] with Brownian motion. Very rad!
References.
1. Stroock, Daniel W.; Varadhan, S. R. S. Diffusion processes. Proceedings of the Sixth Berkeley Symposium on Mathematical Statistics and Probability, Volume 3: Probability Theory, 361--368, University of California Press, Berkeley, Calif., 1972. https://projecteuclid.org/euclid.bsmsp/12005143462. Applied Stochastic Analysis. https://cims.nyu.edu/~holmes/teaching/asa2017.html
3. Koralov, Leonid, Sinai, Yakov G. Theory of Probability and Random Processes. 2007. Springer.