Given an elliptic PDE with Neumann boundary condition \begin{align} \left\{ \begin{aligned} -\sum_{i,j=1}^N\partial_i(a_{ij}\partial_j u)+cu&=f &&\mbox{in}\,\,\,\Omega, \\ \sum_{i,j=1}^Na_{ij}\partial_j un_i&=0 &&\mbox{on}\,\,\,\partial\Omega , \end{aligned} \right. \end{align} in a bounded smooth domain $\Omega$, where $a_{ij}\in W^{1,\infty}(\Omega)$ and $c\in L^\infty(\Omega)$, and the following strong ellipticity condition holds for some constant $K$: \begin{align} &K^{-1}\sum_{i=1}^N|\xi_i|^2\le \sum_{i,j=1}^Na_{ij}(x)\xi_i\xi_j \le K\sum_{i=1}^N|\xi_i|^2 , \quad \forall\, \xi_i\in{\mathbb R},\,\, \forall\, x\in\Omega,\\ &c(x)\ge c_0\,\,\,\mbox{for some constant $c_0>0$}, \end{align}

Is there any reference which says that \begin{align} \|u\|_{W^{2,p}(\Omega)}\le C_p\|f\|_{L^p(\Omega)},\quad\forall\, 1<p<\infty \,\, ? \end{align} All the references I found are only about Dirichlet boundary condition.