There appears to be a connection between the heat kernel and Maxwell-Boltzmann distribution, but I have not seen this in the literature before. I'd appreciate any kind comments or corrections/suggestions on this line of thinking (EDIT: especially with respect to, Eq. \eqref{1}, \eqref{2}, \eqref{3}).
The heat kernel is the fundamental solution to the heat equation $\text{d} K/\text{d} t = D\nabla^2 K$ on a particular domain with a given set of boundary conditions, and in the case of a single spatial dimension for $t>0$, it gives the probability of a particle being displaced by a position $\Delta x \equiv x-x_0 = x(t)-x(t_0)$, with $\Delta t=t-t_0$, such that \begin{align} \tilde K(x,t)\ \text{d} x=\sqrt{\frac{1}{4\pi D \Delta t}}\exp\left[-\frac{(x-x_0)^2}{4D \Delta t}\right] \text{d} x, \end{align} where $D$ is the diffusion coefficient.
By its Markov property, an arbitrary heat kernel can be defined as the product of kernels, defined over successive intervals of displacement and time, such that \begin{align} \tilde K(x_c,t_c) = \tilde K(x_b, t_b) \tilde K(x_a, t_a),\label{eq} \end{align} where $x_c>x_b>x_a$ and $t_c>t_b>t_a$. Therefore, without loss of generality, one can define a kernel for a displacement in an infinitesimal time interval, from which arbitrary kernels can be defined. As such, we consider the limiting case as $\Delta t \to 0$: \begin{align} \lim_{\Delta t \to 0} \tilde K(x,t)\ \text{d} x &=\lim_{\Delta t \to 0} \sqrt{\frac{1}{4\pi D \Delta t}}\exp\left[-\frac{(x-x_0)^2}{4D \Delta t}\right] \text{d} x. \end{align} Letting $x_0=0$ and relabeling $x$ as $\Delta x$, we have \begin{align} \lim_{\Delta t \to 0} \tilde K(\Delta x,t)\ \text{d}(\Delta x) &=\lim_{\Delta t \to 0} \sqrt{\frac{1}{4\pi D \Delta t}}\exp\left[-\frac{1}{4D}\frac{(\Delta x)^2}{\Delta t}\right] \text{d} (\Delta x). \end{align} In order to make the connection with the Maxwell-Boltzmann distribution, we introduce a Lebesgue measure on the infinitesimal space of position and time, and we interpret $K(x,t)$ as a probability density function in this space. Now, by replacing the limit with the formal substitution $\Delta t \to \text{d} t$ and $\Delta x \to \text{d} x$, we can write \begin{align} \label{1}\tag{1} \tilde K(\text{d} x,t)\ \text{d}(\text{d} x) &=\sqrt{\frac{1}{4\pi D \text{d}t}}\exp\left[-\frac{1}{4D}\frac{(\text{d} x)^2}{\text{d} t}\right] \text{d} (\text{d} x). \end{align} We now turn to the Maxwell-Boltzmann distribution. For a canonical ensemble of particles of mass $m$ and temperature $T$, the Maxwell-Boltzmann distribution gives the probability of finding a particle with an instantaneous velocity $v$. In terms of a single degree of freedom in an element of the velocity phase space, the distribution takes the form \begin{align} f(v) \ \text{d} v &=\sqrt{\frac{m}{2\pi k_BT}}\exp\left[-\frac{mv^2}{2k_BT}\right] \ \text{d} v, \end{align} where the mean square speed is $\langle v^2 \rangle= k_BT/m$ and $k_B$ is the Boltzmann constant. Making the substitution $v= \text{d} x /\text{d} t$ and introducing factors of $\text{d} t$, we have \begin{align} f(v) \ \text{d} v &=\sqrt{\frac{m}{2\pi k_BT(\text{d} t)^2}}\exp\left[-\frac{m}{2k_BT}\left(\frac{\text{d} x}{\text{d} t}\right)^2\right] \ \text{d} t \ \text{d} v. \end{align} Now, by virtue of the fact that $\text{d} v / \text{d} t = \text{d}^2 x / \text{d} t^2$, such that $\text{d} t \ \text{d} v = \text{d}^2 x = \text{d}(\text{d} x)$, we find \begin{align} \label{2}\tag{2} f(v) \ \text{d} v &=\sqrt{\frac{m}{2\pi k_BT(\text{d} t)^2}}\exp\left[-\frac{m}{2k_BT\text{d} t}\frac{(\text{d} x)^2}{\text{d} t}\right] \ \text{d}(\text{d} x), \end{align} such that we can equate \eqref{1} and \eqref{2}. Given this equality, the diffusion coefficient can be related to the mean square speed as \begin{align} \label{3}\tag{3} 2D&= \lim_{\Delta t \to 0} \ \frac{k_BT \Delta t} {m} = \lim_{\Delta t \to 0} \ \langle v^2 \rangle \Delta t. \end{align} Moreover, as the mean square displacement, $\langle d^2 \rangle$, in one dimension satisfies the identity \begin{align} \langle d^2 \rangle = \langle (x(t) - x_0)^2\rangle = 2D \Delta t, \end{align} it therefore follows that \begin{align} \langle v^2 \rangle = \lim_{\Delta t \to 0} \frac{\langle d^2 \rangle}{\Delta t^2}. \end{align} As this is true for an infinitesimal time-step, by the product property of kernels, it follows for arbitrarily large time intervals, as well. Integrating over all possible displacements for each kernel factor, one can write the resulting kernel as \begin{align} \tilde K(x,t) &= \lim_{\Delta t \to 0} \int_{-\infty}^{\infty} \dots \int_{-\infty}^{\infty} \left(\frac{1}{4\pi D \Delta t}\right)^{n/2}\prod_{j=1}^n \exp\left[ -\frac{(x_j-x_{j-1})^2}{4D \Delta t} \right] \ \text{d} x_{n-1} \dots \text{d} x, \end{align} which is precisely the Euclidean path integral definition of the kernel. In each of the $n$ infinitesimal intervals, the identity $\langle v^2 \rangle= \lim_{\Delta t \to 0} \langle d^2 \rangle/\Delta t^2$ holds. Hence, the identity holds for a finite time interval in general: \begin{align} \langle v^2 \rangle = \frac{\langle d^2 \rangle}{\Delta t^2}. \end{align}