Let \(D\subset\subset \mathbb R^m\) be a domain. For \(l\in\{0,\dots,m\}\), \( i,j\in\{1,\dots,m\}\) let \(A^l_{i j}\in \mathbb L^2(\mathbb R^m)\) be such that \(A^l_{ij}=-A^l_{ji}\). If \(\mathcal R_l\) denotes the \(l\)th Riesz transform, \(l=1,\dots, m\), and \(\mathcal R_0\) is the identity in \(\mathbb R^m\), set \[ \Omega_{ij} =\sum_{l=0}^mA^l_{i j}\mathcal R_l. \]
The main result of this paper is the following: Let \(1<p<\infty\). If \(w \in L^2(\mathbb R^m)\) is a solution to the system \[ \int w_i|\nabla|^\mu\varphi=-\int\Omega_{ik}[w_k]\varphi,\quad\text{for all } \varphi\in C^\infty_0(D), \] where \(\mu\leq \min\{1,m/2\}\) or \(\mu= m/2\), and \(w\) also satisfies \[ \sup_{B_\rho\subset D}\rho^{\frac(2\mu-m){2}}||w||_{2,B_\rho}<\infty \] then \(w\in L^p_{\mathrm{loc}}(D)\), provided that all the quantities \[ \sup_{B_\rho(x),x\in D}\rho^{\frac{(2\mu-m)}{2}}||A^l_{ij}| |_{2,B_\rho} \] are sufficiently small.