We introduce and study an extension of the heat equation relevant to relativistic energy formula involving square root of differential operators. We furnish exact solutions of corresponding Cauchy (initial) problem using the operator formalism invoking one-sided Lévy stable distributions. We note a natural appearance of Bessel polynomials which allow one the obtention of closed form solutions for a number of initial conditions. The resulting relativistic diffusion is slower than the non-relativistic one, although it still can be termed a normal one. Its detailed statistical characterization is presented in terms of exact evaluation of arbitrary moments and is compared with the non-relativistic case. In this note we develop an alternative point of view which takes advantage of the concepts of anomalous diffusion [7] governed by the evolution equations involving square root pseudo-differential operators. We postulate and study in this note the following one-dimensional differential equation for the distribution function ψ(ξ, τ ) which we believe would be appropriate in the relativistic (R) situation: which is formally of infinite order in ∂ 2 x . For the purpose of this work we shall refer to Eq. (3) as a relativistic heat equation. In Quantum Mechanics the difficulty to treat the square root in Eqs.(1) and (3) (so-called pseudodifferential operators [11,12]) is usually resolved by introducing the wave functions with several components, like in Pauli and in Dirac equations [11,12]. However the interest for equations involving pseudo-differential operators arose even before the inception of the Dirac equation. It can be traced back to H. Weyl [13]. The fact that Eqs. (1) and (3) correctly reproduce the NR limit is reminiscent of the R extension of the analogy between the Schrödinger and the heat equations, carried out in 1984 in [4], by identifying the underlying Poisson proarXiv:1610.05605v1 [math-ph]