There are many places in the literature where the positivity of some semigroups is treated. However I did not know anyone which states and proves the strong positivity even for the basic semigroups like the Neumann laplacian semigroup.
Here is a simplified mathematical problem:
$$ \begin{cases}\dfrac{\partial u}{\partial t}(t,x)-\Delta u(t,x)+u(t,x)=0, & (t,x)\in (0,T)\times\Omega \\ \dfrac{\partial u}{\partial\nu}(t,x)=0, & (t,x)\in (0,T)\times\partial\Omega \\ u(0,x)=f(x), & x\in\Omega \end{cases}, $$
where $\Omega\subset\mathbb{R}^N$ is an open, bounded set with smooth boundary, and $f\in L^{\infty}(\Omega)^+=\{g\in L^{\infty}(\Omega)\ |\ g(x)\geq 0,\ \text{a.e. on}\ \Omega\}$.
If we denote by $S(t)$ the semigroup generated by $-\Delta+I$ with Neumann b.c. on $L^2(\Omega)$ then $S(t)$ is a positive semigroup, i.e. $u(t,\cdot)=S(t)f\in L^{\infty}(\Omega)^+$ for any $t\in [0,T]$. See for example W. Arendt - Heat Kernels (Theorem 3.3.1) https://www.uni-ulm.de/fileadmin/website_uni_ulm/mawi.inst.020/arendt/downloads/internetseminar.pdf.
In some articles I found, without proof or references that if $f\in L^{\infty}(\Omega)^+$ and $f\neq 0$ then $u(t,\cdot)\in \text{int}(L^{\infty}(\Omega)^+)$, i.e. there is some constant $c(t,u_0)>0$ such that $u(t,x)>c(t,u_0)$ a.e. on $\Omega$. How can we prove that?
It looks like a parabolic Harnack-type inequality is needed here...