If we have nonnegative $V \in L^1_{\textrm{loc}}(\mathbb{R}^{n})$, then the operator $H = -\Delta + V$ can be defined on $L^{2}(\mathbb{R}^{n})$ via quadratic form methods. This is done by, for example, E.B. Davies in Heat Kernels and Spectral Theory. He also proves that $H$ is the infinitesimal generator of an ultracontractive symmetric Markov semigroup $e^{-Ht}$ which is given by integration against a kernel $p(x,y,t)$:
$$e^{-Ht}[f] = \int_{\mathbb{R}^{n}} p(x,y,t)f(y)\,dy \qquad f \in L^2(\mathbb{R}^{n})$$
Now I have seen several authors use, directly or indirectly, the fact that for $y$ fixed, $p(\cdot,y,t)$ is a weak solution of the equation $(\partial_{t} + H)u = 0$. But how exactly do we know this is true? I can't find anything about this in Davies.