  The paper focuses on the $L^{p}$-Positivity Preservation property ($L^{p}$-PP
for short) on a Riemannian manifold $(M,g)$. It states that any $L^p$ function
$u$ with $1<p<+\infty$, which solves $(-\Delta + 1)u\ge 0$ on $M$ in the sense
of distributions must be non-negative. Our main result is that the $L^{p}$-PP
holds if (the possibly incomplete) $M$ has a finite number of ends with respect
to some compact domain, each of which is $q$-parabolic for some, possibly
different, values $2p/(p-1) < q \leq +\infty$. When $p=2$, since
$\infty$-parabolicity coincides with geodesic completeness, our result settles
in the affirmative a conjecture by M. Braverman, O. Milatovic and M. Shubin in
2002. On the other hand, we also show that the $L^{p}$-PP is stable by removing
from a complete manifold a possibly singular set with Hausdorff co-dimension
strictly larger than $2p/(p-1)$ or with a uniform Minkowski-type upper estimate
of order $2p/(p-1)$. The threshold value $2p/(p-1)$ is sharp as we show that
when the Hausdorff co-dimension of the removed set is strictly smaller, then
the $L^{p}$-PP fails. This gives a rather complete picture. The tools developed
to carry out our investigations include smooth monotonic approximation and
consequent regularity results for subharmonic distributions, a manifold version
of the Brezis-Kato inequality, Liouville-type theorems in low regularity,
removable singularities results for $L^{p}$-subharmonic distributions and a
Frostman-type lemma.
  Since the seminal works by T. Kato, the $L^{p}$-PP has been linked to the
spectral theory of Schr\"odinger operators with singular potentials $\Delta -
V$. Here we present some applications of the main results of this paper to the
case where $V\in L^p_{loc}$, addressing the essential self-adjointness of the
operator when $p=2$ and whether or not $C^\infty_c(M)$ is an operator core for
$\Delta-V$ in $L^p$.
