  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, $b$ and $c$ two $d$-dimensional vector-valued functions on
$\Omega$ with $L^{\infty}$ coefficients and $V$ a locally integrable nonegative
function on $\Omega$. Consider the operator ${\mathscr L}^{A,b,c,V}=-{\rm
div}\,(A\nabla) + \left\langle \nabla , b \right\rangle - {\rm div}\,(c \,
\cdot) + V $ with mixed boundary conditions on $\Omega$. We extend the bilinear
inequality that Carbonaro and Dragi\v{c}evi\'c proved in the special cases when
$b=c = 0$. As a consequence, we obtain that the solution to the parabolic
problem $u^{\prime}(t)+{\mathscr L}^{A,b,c,V}u(t)=f(t)$, $u(0)=0$, has maximal
regularity in $L^{p}(\Omega)$, for all $p>1$ such that $A$ satisfies the
$p$-ellipticity condition that Carbonaro and Dragi\v{c}evi\'c introduced in
arXiv:1611.00653 and $b,c,V$ satisfy another condition that we introduce in
this paper. Roughly speaking, $V$ has to be ``big'' with respect to $b$ and
$c$. We do not impose any conditions on $\Omega$, in particular, we do not
assume any regularity of $\partial\Omega$, nor the existence of a Sobolev
embedding.
