Let $u$ be the weak solution on a smooth bounded domain $\Omega \subset \mathbb{R}^n$ (for $n \leq 3$) of $$u_t - \Delta u = f$$ $$u(0) = u_0$$ $$\partial_\nu u = 0 \quad\text{on $\partial\Omega$}$$ for $u_0 \in L^\infty(\Omega)$ (non-negative) and $f \in L^\infty(0,T;L^2(\Omega))$. We know that $$u(t,x) \leq C(\lVert u_0 \rVert_\infty + \lVert f \rVert_{L^\infty(0,T;L^2)})$$ is satisfied almost every $(x,t)$ by De Giorgi's method (eg. Theorem 4.2.2 in "Elliptic and Parabolic Equations" book by Wu, Yin and Wang).

If instead we had the equation $$u_t - \Delta u +\alpha u= f$$ for a constant $\alpha > 0$ (and all else is the same as above), can I expect some kind of decay on the $L^\infty$ estimate involving $\alpha$, eg. $$u(t,x) \leq \frac{C}{\alpha}(\lVert u_0 \rVert_\infty + \lVert f \rVert_{L^\infty(0,T;L^2)})?$$ It is essential for me that $f$ is only in $L^\infty(0,T;L^2)$.

Of course I tried using an integrating factor, but this gives a bad contribution on the term involving the $f$ norm. I wonder if it can be true.