Let \(\Omega \subseteq \mathbb {R}^d\) be open, A a complex uniformly strictly accretive \(d\times d\) matrix-valued function on \(\Omega \) with \(L^\infty \) coefficients, and V a locally integrable function whose negative part is subcritical. We consider the operator \({{\mathscr {L}}}= -\textrm{div}(A\nabla ) + V\) with mixed boundary conditions on \(\Omega \) . We extend the bilinear embedding of Carbonaro and Dragičević (Publ Mat 69(1):83–108, 2025), established for nonnegative potentials under the p-ellipticity of the matrices, by introducing a novel condition on the coefficients that reduces to p-ellipticity when V is nonnegative. As a consequence, the solution to the parabolic problem \(u^\prime (t) + {{\mathscr {L}}}u(t) = f(t)\) with \(u(0)=0\) has maximal regularity on \(L^p(\Omega )\) . Moreover, under this new condition, we study mapping properties of the semigroup generated by \(-{{\mathscr {L}}}\) , thereby extending classical results for the Schrödinger operator \(-\Delta + V\) on \(\mathbb {R}^d\) .