Mean value formula for fractional heat equation For the solution $u(z) = u(t,x)$ of the heat equation $u_t -\Delta u = 0$ we have
$$u(z_0) = \int_{\Omega_r(z_0)}u(z) K_r(z_0-z) dz,$$
where $$\Omega_r(z_0) = \left\{z \in \mathbb{R}^{N+1}: \Gamma(z_0-z) > 1/r\right\} \quad K_r(z) = K_r(x,t) = \frac{|x|^2}{4rt^2}$$
and $\Gamma(z) = \Gamma(x,t)$ is the fundamental solution at the origin of the heat equation.
Does something similar hold for the fractional heat equation $u_t + (-\Delta)^su=0$?
 A: The answer depends on what you mean by "something similar". :-)

The expression for the classical heat equation has the following probabilistic interpretation. Let $X_t$ be the $d$-dimensional Brownian motion, started at a given point $x$, and let $\mathbb P^x$ denote the corresponding probability. As a consequence of Itô's lemma, if $u$ is a solution of the heat equation, then the process
$$ M_t = u(t, X_t) $$
is a martingale. Thus, for any "reasonable" Markov time $\tau$ we have
$$ u(0, x) = M_0 = \mathbb E^x M_\tau = \mathbb E^x u(\tau, X_\tau) . $$
Now if $\tau$ is the first exit time from a set $\Omega$, then $(\tau, X_\tau)$ is concentrated on the boundary of $\Omega$, and we get a measure $\pi_\Omega$ such that
$$ u(0, x) = \int_{\partial \Omega} u(t, y) \pi_\Omega(dt, dy) . $$
Considering a family $\Omega_r$ of sets contained in $\Omega$ (for example, the ones defined in the question for $z_0 = (0, x)$) and averaging appropriately over $r$, we get
$$ u(0, x) = \int_\Omega u(t, y) \pi(dt, dy) $$
for a measure $\pi$ that is the mixture of measures $\pi_{\Omega_r}$. If we do this properly, $\pi$ will have a "nice" density.
(Of course, actually carrying out the above procedure to get an explicit formula for $\pi(dt, dy)$ is a different story.)

For the fractional heat equation, one can carry out the same procedure, replacing the Brownian motion by an appropriate stable Lévy process. However, in this case $(\tau, X_\tau)$ is no longer distributed on the boundary of $\Omega$: its support is typically the complement of $\Omega$ in the smallest strip $\mathbb R^d \times [0, T]$ that contains $\Omega$.
As a consequence, the measure $\pi_\Omega$ is no longer concentrated on the boundary, and so the support of $\pi$ is not compact. So we do have a similar formula, but with an integral over $[0, T] \times \mathbb R^d$ rather than over $\Omega$ — which seems pretty useless in applications, as we can express $u$ in terms of the fractional heat kernel straight away, with no effort.

On the other hand, one should be able to get a similar formula with an integral over a compact set, but with a signed kernel $K$. That said, I do not know whether this has been studied, let alone an explicit expression for the kernel $K$.
