In the free case one can compute the resolvents of the Laplacian $-\Delta$ in many cases explicitly, in the sense that they are given by an integral operator. Often, one uses the Hille-Yosida theorem or Fourier transform to do so.

Also, the Schwartz-kernel theorem tells us that in a distributional sense there is always an "integral kernel" associated to the resolvent of a Schrodinger operator.

However, examples in the free case suggest that the resolvent of $H:=-\frac{d^2}{dx^2} +V$ ( I am now in dimension $1$ as I would just like to get some conceptual insight) where $V$ is a $C^{\infty}$ and bounded potential with suitable boundary conditions such that this operator is self-adjoint, gives rise to a resolvent

$$(H-z)^{-1}\psi(x) = \int K_z(x,y)\psi(y) dy$$ where $K_z$ has, for fixed $z$ in the resolvent set, nice properties (almost always smooth) in $x$ and $y$.

See for example this calculation: of the one-dimensional resolvent or this one of the three-dimensional on math.stackexchange

I would like to ask why this is the case? Which theorem guarantees us that such nice $K_z$ kernels exist?