Weighted Orlicz gradient estimates for the class of singular p-Laplace system

Let $n \in \{2, 3, 4, \ldots\}$, $N \in \{1, 2, 3, \ldots\}$ and $p \in \big(1, 2-\frac{1}{n}\big]$. Let $\beta \in (1,\infty)$ be such that $$\frac{np}{n-p}<\beta'<\frac{n}{n(2-p)-1}$$ and $f \in L^{\beta}(\mathbb R^n;\mathbb R^N)$. Consider the $p$-Laplace system $$-\Delta_p u=-\operatorname{div}\big(|Du|^{p-2}Du\big)=f \quad in \quad \mathbb R^n.$$ We obtain a weighted gradient estimate for distributional solutions of this system.