Heat Kernel Asymptotics on Manifold with Boundary This is crosspost from math.stackexchange https://math.stackexchange.com/questions/311213/heat-kernel-asymptotics-on-manifold-with-boundary where the question did not yield any answer
On a closed Riemannian manifold $M$, the heat kernel $k_t(x, y)$ of the Laplace-Beltrami operator (or more general of any generalized symmetric Laplace-type operator acting on sections of a vector bundle) admits an asymptotic expansion of the form
$$  k_t(x, y) \sim \exp\left( -\frac{1}{4t}d(x, y)^2\right) \sum_{j=0}^\infty t^j \Phi_j(x, y) $$
where $d(x, y)$ denotes the Riemannian distance and $\Phi_j$ are appropriate smooth functions, not depending on $t$. This is meant in the sense that for each $N \in \mathbb{N}$, there exists a constant $C>0$ such that for all $x, y \in M$, we have
$$ \left| k_t(x, y) -  \chi(x, y)\exp\left( -\frac{1}{4t}d(x, y)^2\right) \sum_{j=0}^Nt^j \Phi_j(x, y) \right| < C t^{N+1}$$
where $\chi(x, y)$ is an appropriate cutoff function that is $\equiv 1$ near the diagonal.
In the case that $M$ is still compact but has a boundary, in many books there can be found an asymptotic expansion of the trace, but I could not find an asymptotic expansion of the kernel itself, uniform on $M \times M$. Is there such an expansion?
Remark 1: When $M$ has a totally geodesic boundary, such an expansion is easy to get with help of the "Riemannian double". But of course, this case is quite unlikely (it is not even fulfilled for domains in Euclidean space).
Remark 2: I did not specify any boundary conditions, but one can assume that we are in the simplest case, i.e. the Laplace-Beltrami operator acting on functions with either Dirichlet or Neumann boundary conditions, whichever you like.

Edit: Clearly, the same asymptotics as above cannot hold on a manifold with boundary. One could however expect results like
$$\Bigl|k_t(x, y) - (2\pi t)^{-n/2}e^{-\frac{d(x, y)^2}{4t}} \pm (2\pi t)^{-n/2}e^{-\frac{\sigma(x, y)^2}{4t}} \Bigr| \leq C t^{-n/2+\varepsilon},$$
where $\sigma(x, y)$ is the distance "over the boundary" (i.e. the length of the shortest path from $x$ to $y$ that touches the boundary somewhere) and the sign $\pm$ depends on the boundary conditions (+ for Dirichlet, - for Neumann). (Addendum: The whole statement is: For every compact set $K$ in the $M\times M$ minus cutpoints (which may include boundary points), there exists a constant $C$ such that the above statement is true for. $x,y \in K$).
Results like this come from physical arguments (like reflecting paths etc.) and are used with handwavy arguments in physics literature, but I do not know if something like this is ever rigourously proved (probably this is even wrong). 
Edit: This result does hold in the case that the boundary is totally geodesic, but will most definitely be wrong otherwise.
 A: Uniformity on compact sets interior to the domain
The boundary condition has no real  impact on the small-time asymptotics of the heat kernel for interior points. More precisely, consider a Riemannian manifold $\mathbb{M}$ and $\Omega$ a relatively compact smooth open domain in $\mathbb{M}$. Denote by $p_\Omega$ the heat kernel associated with a self-adjoint extension of the Laplace-Beltrami operator on $\Omega$ and $p_\mathbb{M}$ the usual heat kernel on $\mathbb{M}$. If $K$ is a compact subset of $\Omega$, then we have for $ t \le T$, $x,y \in K$
$
| p_\mathbb{M}(t,x,y) - p_\Omega(t,x,y) | \le C_1 e^{-C_2/t} \quad  (\star)  
$
where $C_1$ depends on $T$ and $C_2$ on $K$. As a consequence the small time asymptotics of  $p_\mathbb{M}$ leads to a small-time asymptotics of $ p_\Omega$ and we obtain the following theorem:

Let $K$ be any compact subset of $\Omega$.  For every $N \ge 0$, there exist $\delta> 0$, smooth functions $\Phi_0,  \Phi_1, \cdots, \Phi_N$ and a constant $C$ such that for every $x,y \in K$ with $d(x,y) \le \delta$, $t \in [0,1]$,
  $
\left| p_\Omega(t,x,y) -\frac{e^{-\frac{d^2(x,y)}{4t}} }{t^{n/2}} \left( \sum_{k=0}^N  \Phi_k (x,y) t^k \right) \right| \le C \frac{e^{-\frac{d^2(x,y)}{4t}} }{t^{n/2-(N+1)}}
$

This is how we can prove the inequality $ (\star)  $.  Fix $x \in \Omega$, and consider a smooth domain $\tilde{\Omega}$ that contains $x,y$. Denote by $p^D_{\tilde{\Omega}}$ the Dirichlet heat kernel on $\tilde{\Omega}$. We have
$
p_\Omega(t,x,y)=p^D_{\tilde{\Omega}}(t,x,y)+\mathbb{E}_x \left( 1_{ t \ge T_{\tilde{\Omega}}} p_\Omega(t-T_{\tilde{\Omega}},X_{T_{\tilde{\Omega}}},y) \right)
$
where $X$ is the diffusion process associated with $p_\Omega$ and $T_{\tilde{\Omega}}$ the hitting time of the boundary of $\tilde{\Omega}$.
We easily deduce from that
$
0\le p_\Omega (t,x,y) - p^D_{\tilde{\Omega}}(t,x,y) \le C_1 \mathbb{P}_x ( t \ge T_{\tilde{\Omega}} ) \le C_1 e^{-C_2/t}$
Similarly we obtain
$
0 \le p_\mathbb{M}(t,x,y) -p^D_{\tilde{\Omega}}(t,x,y)   \le C'_1 e^{-C'_2/t}
$
and we get the inequality  $ (\star)  $ by combining the two above inequalities.
Uniformity on compact sets intersecting the boundary
In general, we can not hope for any estimate  of the type
$
\left|  p_\Omega (t,x,y) -\frac{e^{-\frac{d^2(x,y)}{4t}} }{t^{n/2}} \Phi (x,y)\right| \le \frac{C}{t^{n/2-\varepsilon}}
$
which would be uniform on a compact $K$ intersecting the boundary of $\Omega$. Indeed, take for $p_\Omega$ the Dirichlet heat kernel. If the above estimate is satisfied we would have
$
\left| t^{n/2} p_\Omega (t,x,x) - \Phi (x,x)\right| \le C t^{\varepsilon}
$
If $x \in \partial \Omega \cap K$, then $ p_\Omega (t,x,x) =0$, so we deduce $\Phi (x,x)=0$. If $x \in \Omega \cap K$, $\lim_{t \to 0} t^{n/2} p_\Omega (t,x,x) =\frac{1}{(4\pi)^{n/2}}$, so $\Phi(x,x)=\frac{1}{(4\pi)^{n/2}}$. 
So $\Phi(x,x)$ is not continuous, which is a contradiction because $t^{n/2} p_\Omega (t,x,x) $ is continuous and supposed to converge uniformly to $\Phi(x,x)$.
However as OP observed, in some cases like when the boundary is totally geodesic and the boundary condition is Neumann or Dirichlet, there exists a modified asymptotics 
$
\left|  p_\Omega (t,x,y) -\left( \frac{e^{-\frac{d^2(x,y)}{4t}} }{t^{n/2}} \pm \frac{e^{-\frac{\sigma^2(x,y)}{4t}} }{t^{n/2}}  \right) \Phi (x,y)\right| \le \frac{C}{t^{n/2-\varepsilon}}
$
which holds uniformly on a compact intersecting the boundary. 
We can not   hope for such "reflection principle" in  great generality for the Neumann kernel  because under some convexity assumption  of the boundary we have for $x,y \in \partial{\Omega}$
$
p_\Omega (t,x,y)\simeq \frac{\Phi_1 (x,y)}{t^{n/2+1/6}}  e^{-\frac{d^2(x,y)}{4t}-\frac{\Phi_2( d(x,y))}{t^{1/3}}} $
where $\Phi_1$ and $\Phi_2$ are directly related to the normal curvature of the boundary.
The result was proved by Ikeda and Hsu, using probabilistic methods
SHORT-TIME ASYMPTOTICS OF THE HEAT KERNEL ON A CONCAVE BOUNDARY
