  Let $\Omega\subseteq \mathbb{R}^{d}$ be open and $A$ a complex uniformly
strictly accretive $d\times d$ matrix-valued function on $\Omega$ with
$L^{\infty}$ coefficients. Consider the divergence-form operator ${\mathscr
L}^{A}=-{\rm div}(A\nabla)$ with mixed boundary conditions on $\Omega$. We
extend the bilinear inequality that we proved in [16] in the special case when
$\Omega=\mathbb{R}^{d}$. As a consequence, we obtain that the solution to the
parabolic problem $u^{\prime}(t)+{\mathscr L}^{A}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 we introduced in [16]. This range of
exponents is optimal for the class of operators we consider. 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. The methods of [16]
do not apply directly to the present case and a new argument is needed.
