Suppose $M$ is a Riemannian manifold with some compact quotient under isometries.
Associated with the Riemannian metric one has the Laplace-Beltrami operator $\Delta$ and the heat kernel $p(t,x,y)$ which is positive on $(0,+\infty) \times M \times M$ and satisfies:
$\partial_t p(t,x,\cdot) = \Delta p(t,x,\cdot)$
$\lim_{t \to 0}f(y)p(t,x,y)\mathrm{d}y = f(x)$
among other properties (e.g. $p(t,x,y) = p(t,y,x)$).
Suppose one replaces the Laplacian by a second order differential operator which is given in local coordinates by: $\sum_{i,j}a_{ij}\partial_{i}\partial_{j} + \sum_{k}b_k \partial_k$ where the coeficients $a_{ij},b_k$ vary smoothly and the matrix $(a_{ij})$ is symmetric and positive definite.
My question is: What are the possible pitfalls if one tries to generalize properties from the Laplace-Beltrami case to this case?
I know (e.g. from Candel "The harmonic measures of Lucy Garnett") that one still has a "heat kernel" $p(t,x,y)$ but that it may not be symmetric (i.e. $p(t,x,y) \neq p(t,y,x)$). It is claimed in said reference that estimates such as that of Cheng-Li-Yau also hold for these types of kernels. Furthermore "Harnack-type" estimates on what would be the analog of harmonic functions also hold as seen in Serrin's "On the Harnack inequality for linear elliptic equations.".
I'd be interested in a reference where differences and similarities with the Laplace-Beltrami case are discussed and worked out.
