In $\mathbb R^d$, $d \geq 3$, consider the divergence and the non-divergence form operators \begin{equation} \tag{$i$} -Δ- \nabla \cdot (a-I) \cdot \nabla + b \cdot \nabla, \end{equation} \begin{equation} \tag{$ii$} - Δ- (a-I) \cdot \nabla^2 + b \cdot \nabla, \end{equation} where the second order perturbations are given by the matrix $$a-I=c|x|^{-2}x \otimes x, \quad c>-1.$$ The vector field $b:\mathbb R^d \rightarrow \mathbb R^d$ is form-bounded with the form-bound $δ>0$ (this includes a sub-critical class $[L^d + L^\infty]^d$, as well as vector fields having critical-order singularities). We characterize quantitative dependence on $c$ and $δ$ of the $L^q \rightarrow W^{1,qd/(d-2)}$ regularity of the resolvents of the operator realizations of ($i$), ($ii$) in $L^q$, $q \geq 2 \vee ( d-2)$ as (minus) generators of positivity preserving $L^\infty$ contraction $C_0$ semigroups.