Let $\Omega\subset \mathbb{R}^d$ be a Lipschitz bounded domain. In $L_2(\Omega;{\mathbb C}^n)$, we consider a strongly elliptic second order operator $A$ with constant coefficients. Namely, it is assumed that $A$ is given by the differential expression $b(D)^* b(D)$ with Neumann boundary condition. Here $b(D)= \sum_{l=1}^d b_l D_l$ is a first order differential operator, where $b_l$ are constant $(m\times n)$-matrices. It is assumed that $m \ge n$ and the symbol $b(\xi) = \sum_{l=1}^d b_l \xi_l$ is subject to the condition $$ {\rm rank}\, b(\xi) =n,\ \ 0\ne \xi \in \mathbb{C}^d, $$ which ensures the coercivity. Let $u \in H^1(\Omega;\mathbb{C}^n)$ be a (weak) solution of the equation $$ (A u)(x) = \sum_{j=1}^d D_j F_j(x) $$ with Neumann (natural) boundary condition and $\int_\Omega u(x)\,dx=0$. This is understood as the integral identity $$ \intop_\Omega \langle b(D) u(x), b(D) \eta(x) \rangle\, dx = \sum_{j=1}^d \intop_\Omega \langle F_j(x), D_j \eta(x) \rangle\, dx, \quad \forall \eta \in H^1(\Omega;\mathbb{C}^n). $$
QUESTION: Assume that $F_j \in L_p(\Omega;\mathbb{C}^n)$ with some $2\le p < \infty$. Is this true that $u \in W^1_p(\Omega;\mathbb{C}^n)$ and we have an estimate $$ \| u \|_{W_p^1(\Omega)} \le C \sum_{j=1}^d\|F_j\|_{L_p(\Omega)} ? $$ Is it true for any Lipschitz bounded domain or for more smooth domains? What is the reference?
