\section{Introduction}

While uniformly elliptic partial differential equations are in general well understood, degenerate elliptic PDEs arising in the case when the lower bound for the ellipticity constant is not bounded away from zero present some challenges. They are often studied from the perspective that they \emph{close} to a uniformly elliptic PDE. In this paper we propose a different approach for 
operators that fail to be uniformly elliptic but which are $\gamma$-homogeneous. More precisely we focus on operators which arise naturally in the study of anisotropic geometric variational problems. Namely we study:


\begin{equation} 
\label{e:homo1} \int_{\Omega} \mathfrak{D}(u,\varphi) = 0 \quad \varphi \in W^{1,\gamma^{\prime}}_{0}(\Omega),
\end{equation}
where $\mathfrak{D} : W^{1,\gamma}(\Omega) \times W^{1,\gamma^{\prime}}_{0}(\Omega) \to W^{1,1}(\Omega)$ is defined by
\begin{equation}  \label{e:integrand}
\mathfrak{D}(u, \varphi) = \left \langle \rho(x,Du)^{\gamma-1}(D \rho(x,\cdot))(Du), D \varphi \right \rangle,
\end{equation}
for some function so norm $\rho(x, \cdot)$ is $1$-homogeneous and strictly convex and $1 < \gamma < \infty$. 


A Moser iteration approach allows us to exploit the $\gamma$-homogeneity of the PDE above to prove a Harnack inequality and a Bernstein type theorem for global solutions. These two results arise as corollaries of the theorems described below.


Throughout this paper, for any real number $\alpha > 1$ we let $\alpha^{\prime}$ denotes the Holder conjugate of $\alpha$. Also, whenever $\rho(x,\cdot)$ is a $1$-homogeneous function $\rho_{*}(x,\cdot)$ is its convex dual defined by $\rho_{*}(x, \zeta) = \sup_{\rho(x,\xi) < 1} \zeta \cdot \xi$. Throughout this paper we consider the PDE defined for some $\vec{F} \in L^{q}_{\operatorname{loc}}(\Omega), f \in L^{\frac{q}{\gamma^{\prime}}}_{\operatorname{loc}}(\Omega)$, by
\begin{equation} \label{e:inhomo}
\int_{\Omega} \left \langle \rho(x,(Du))^{\gamma-1} (D \rho(x, \cdot))(Du) , D \varphi  \right \rangle = \int_{\Omega} \left \langle \vec{F}, D\varphi \right \rangle + f \varphi.
\end{equation}

\begin{theorem} Suppose that for some $1 < \gamma < n$, $u \in W^{1,\gamma}(\Omega)$ is a non-negative subsolution of \eqref{e:inhomo}, that $\rho: \Omega \times \mathbb{R}^n \setminus \{0\} \to (0,\infty)$ satisfies \eqref{e:meas}, \eqref{e:p1h}, \eqref{e:cx}, and \eqref{e:eb}. Assume $\vec{F} \in L^{q}_{\operatorname{loc}}(\Omega)$ and $f \in L^{\frac{q}{\gamma^{\prime}}}_{\operatorname{loc}}(\Omega)$ for some $q > \frac{n}{\gamma-1}$. If $0 < r < R < 1$ and $B_{R} \subset \Omega$, then for all $0 < p < \infty$, there exists $C = C(n,\gamma,\nu,p)$ so that
\begin{equation} \label{e:supbound}
\sup_{B_{r}} u \le C \left[ (R-r)^{-\frac{n}{p}} \|u\|_{L^{p}(B_{R})} + R^{\delta} \|\rho_{*}(x,\vec{F})\|^{\frac{1}{\gamma-1}}_{L^{q}(B_{R})} + R^{\gamma^{\prime} \delta} \|f\|_{L^{\frac{q}{\gamma^{\prime}}}}^{\frac{1}{\gamma-1}} \right]
\end{equation}
where $\delta = 1 - \frac{n}{q(\gamma-1)} > 0$.
\end{theorem}

In the case when $\gamma = 2 = \gamma^{\prime}$, which includes the elliptic case, \eqref{e:supbound} recovers the local-boundedness results of De Giorgi, Nash, and Moser for elliptic divergence form equations. One can verify this includes the elliptic case with symmetric coefficients by considering $\rho(x,\xi) = |S(x) \xi|$ where $S(x) = A(x)^{1/2}$, which ensures $\langle A(x) \xi, \xi \rangle = |S(x) \xi|^{2}$. The next example demonstrates that the $\gamma = 2$ is not restricted to the uniformly elliptic setting. 

\begin{example} \label{x:p2}
If $\gamma = 2$ and $\rho = \| \cdot \|_{\ell^{p}}$ for some $p > 2$, then \eqref{e:homo1} with integrand \eqref{e:integrand} becomes the PDE $ \divr \left(  D( \| \cdot \|_{\ell^{p}}^{2})(Du) \right) = 0$. A quick computation verifies that 
$$
(\partial_{i} \partial_{j} \| \cdot \|_{\ell^{p}}^{2} )(x) = 2(p-1) \left[ \delta_{ij} \left( \frac{x_{i}}{\|x\|_{p}} \right)^{p-2} - \frac{x_{i} |x_{i}|^{p-2} x_{j} |x_{j}|^{p-2}}{\|x\|_{p}^{2(p-1)}} \right].
$$
Notably, $(D_{2} \| \cdot \|_{\ell^{p}}^{2})(e) \equiv 0$ whenever $e$ is a standard basis vector.
\end{example}

The fact that Example \ref{x:p2} is handled with the methods in this paper demonstrates that the growth rate is what truly matters when studying this form of PDE. The convexity conditions typically imposed when studying elliptic equations are merely forcing a certain growth rate with respect to the Euclidean norm, whereas enforcing a certain growth rate with respect to \emph{any} $C^{1}$-norm satisfying \eqref{e:cx} is just as good.

\begin{example} \label{x:anisotropic}
We revisit the previous example in more generality. Consider $\rho( \cdot) = \| A (\cdot) \|_{\ell^{p}}$ for $p > 2$, where $A$ is an orthogonal matrix. This time, we consider any $\gamma > 2$. The PDE under consideration takes the form
$$
\divr( D(\| A (\cdot ) \|_{\ell^{p}}^{\gamma})(Du)) = 0.
$$
We note that when $A = I$ and $\gamma = p$ this is precisely the so-called pseudo $p$-Laplacian or anisotropic $p$-Laplacian.  By the results here-in we prove a new Harnack inequality for solutions to the pseudo $p$-Laplacian.

By the previous example, whenever $A(Du) \in \{ \pm e_{1}, \pm e_{2}, \dots, \pm e_{n}\}$ the Hessian, $D^{2} \| A( \cdot )\|_{\ell^{p}}^{\gamma}(Du)$, is zero. Performing a Taylor series expansion of $\| A( \cdot )\|_{\ell^{p}}^{\gamma}$ around an arbitrary $\xi_{1} \in \mathbb{R}^n \setminus \{0\}$, you find for $\xi_{2} \in \xi_{1}^{\perp}$ small enough,
$$
\| A( \xi_{1} + \xi_{2}) \|^{\gamma}_{p} - \|A(\xi_{1})\|^{\gamma}_{p} = O \left( \sum_{e_{i} :  A(\xi_{1}) \cdot e_{i} = 0} |\xi_{2} \cdot e_{i}|^{p} \right) + O \left( \sum_{e_{i} : A(\xi_{1}) \cdot e_{i} \neq 0 } |\xi_{2} \cdot e_{i}|^{2} \right).
$$
whenever $|\xi|$ is small enough. This seems to demonstrate that the gains and losses in energy of functional
$$
\int_{\Omega} \| A( Du) \|_{\ell^{p}}^{p}
$$
by small perturbations behaves similarly to a ``mixed growth'' problem. Meanwhile, the corresponding PDE for minimizers is covered as a special case of the work here.

Finally, without straying too far from this same family of PDEs, we note that the techniques below even cover the case when $\rho(x, \xi) = \| A(x)[\xi] \|_{\ell^{p(x)}}$ so long as $p(x) \in (1, \infty)$ and $A(x)$ is reasonably bounded.
\end{example}



The second main result of this paper is for non-negative supersolutions.


\begin{theorem}
Suppose $1 < \gamma < n$ and $u \in W^{1,\gamma}(\Omega)$ is a non-negative supersolution to \eqref{e:inhomo}, for some $\rho$ satisfying \eqref{e:meas} - \eqref{e:eb}. Assume $\vec{F} \in L^{q}_{\operatorname{loc}}(\Omega)$, $f \in L^{\frac{q}{\gamma^{\prime}}}_{\operatorname{loc}}(\Omega)$ for some $q > \frac{n}{\gamma-1}$. If $0 < r < R < 1$ and $B_{R} \subset \Omega$, then for all $0 < p < \frac{n(\gamma-1)}{n-\gamma}$ and all $0 < \theta < \tau < 1$, there exists $C = C(n, \gamma, \nu, \Lambda, q, p, \theta, \tau) > 0$ so that

$$
\inf_{B_{\theta R}} u + R^{\delta} \| \rho_{*}(\vec{F})\|_{L^{q}(B_{R})} + R^{\gamma^{\prime} \delta} \| f\|_{L^{\frac{q}{\gamma^{\prime}}}(B_{R})} \ge C R^{-\frac{n}{p}} \|u\|_{L^{p}(B_{\tau R})},
$$
where $\delta = 1 - \frac{n}{q(\gamma-1)}$.
\end{theorem}

In the case $\gamma = 2 = \gamma^{\prime}$ this recovers Moser's weak Harnack inequality. Moreover, we want to emphasize that the above theorem does not require any \emph{a priori} boundedness of $u$.

As in the classical case, and addressed in Section \ref{s:theorems}, these two main theorems can be combined to achieve a Harnack inequality. Hence as corollaries they prove $C^{\alpha}$-regularity and a Liouville-type theorem for bounded solutions on $\mathbb{R}^{n}$.


There has been substantial work already done regarding minimizers in the setting of multi-phase growth, i.e., minimizers of functionals that originally take the form
$$
\int |Du|^{p} + a(x) |Du|^{q},
$$
for appropriate choices of $p,q$ depending on $n$ and the amount of Holder regularity you may want for your solution. See for instance \cite{acerbi1991partial,marcellini1991regularity, colombo2015bounded,byun2017global,ok2018harnack,baroni2018regularity}. Despite many similarities, all of these works in some way handle a new difficulty not present in the others. Likewise, the class of PDEs studied here have similarities, but to the best of the author's knowledge handle difficulties not covered in the literature to date. For instance, difficulties addressed here that ensure these results are not covered by the aforementiond works are: anisotropy, no type of uniform strict quasiconvexity, $\vec{F}$ being unbounded, the lack of a priori regularity of the solution $u$, and achieving a Harnack inequality with bounds independent of the solution's \emph{a priori} regularity.

Other bodies of literature on similar subjects appear in the likes of variable-exponent growth equations, see for instance, \cite{harjulehto2008harnack, harjulehto2009harnack}, whose results include wildly varying growth rates at different points, while taking advantage of having the same growth rate in all directions at each point. It would be interesting, but outside the scope of this paper, to see if the techniques here can also be applied in a variable-exponent setting. That is: 

{\bf Question: }
Do analogous statements to Theorem \ref{t:sup} and Theorem \ref{t:inf} hold when $\gamma$ is replaced with a function depending on $x\in \Omega$? 

The paper \cite{cupini2009regularity} handles many types of anisotropies than can arise in an integrand, and handles this difficulties with carefully balancing the different exponents while also not multiplying partial derivatives from different directions. Example \ref{x:anisotropic} produces a PDE that falls outside the structure discussed there.

In the preparation of these notes the author learned about the very recent paper \cite{benyaiche2020weak}. After discovering that paper we learned there is also a wealth of literature about minimizers with Orlicz-type growth conditions. See, for instance, \cite{harjulehto2017holder,benyaiche2019harnack,arriagada2018harnack}, and the citations therein. We believe that combining the work of \cite{benyaiche2019harnack,benyaiche2020weak} is how one comes closest to finding the results here-in in the pre-existing literature. 

In the works of Benyaiche, et. al., they study the setting when what we call $\vec{F}, f$ are zero, and using Orlicz spaces they address additional difficulties including including different upper and lower (almost)-homogeneities of the integrand, while proving similar bounds to what we recovered here when $f, \vec{F}$ are both constantly zero. In the isotropic setting, \cite{arriagada2018harnack} shows a Harnack inequality similar to the \ref{t:iwh}. Meanwhile \cite{toivanen2012harnack} proves a general Harnack inequality, while assuming that the terms like $\vec{F}, f$ are in $L^{\infty}$.

Finally, also during preparation of these notes \cite{fazly2020partial} proved a Liouville type theorem when $\gamma = 2$ and $\rho$ has strictly positive Hessian in the sense that 

$$
\lambda^{2} |\zeta|^{2} \le (\partial_{\xi_{i}} \partial_{\xi_{j}} \rho)(\xi)\zeta_{i} \zeta_{j} \le \Lambda |\zeta|^{2} \qquad \forall \zeta \in \xi^{\perp}.
$$
Their techniques made use of an interesting monotonicity formula. 



\section{Preliminaries}

Throughout we will suppose $\rho: \Omega \times \mathbb{R}^n \setminus \{0\} \to (0, \infty)$ is so that 
\begin{equation} \label{e:meas}
\begin{cases}
\rho(x, \cdot ) \in C^{1}(\mathbb{R}^n \setminus \{0\}) & \forall x \in \Omega \\ 
\rho(\cdot, \xi) \in L^{\infty}(\Omega)  & \forall \xi \in \mathbb{R}^n \setminus \{0\}.
\end{cases}
\end{equation}
and that $\rho$ positively $1$-homogeneous function in its second variable, i.e.,
\begin{equation} \label{e:p1h}
\rho( x, \lambda \xi) = \lambda \rho(x, \xi) \qquad \forall \lambda > 0, \quad  \forall x \in \Omega, \quad \forall \xi \in \mathbb{R}^n \setminus \{0\}.
\end{equation}
We further assume that $\rho(x, \cdot)$ is strictly convex in the sense that
\begin{equation} \label{e:cx}
\rho(x, \xi_{1} + \xi_{2}) < \rho(x, \xi_{1}) + \rho(x, \xi_{2}) \quad \forall x \in \Omega, \quad \forall \xi_{1} \neq \lambda \xi_{2} \in \mathbb{R}^n \setminus \{0\} .
\end{equation}
Finally, we assume there exists $0 < \nu \le \lambda < \infty$ independent of $x$ so that
\begin{equation} \label{e:eb}
 \nu \le \rho(x, \xi) \le \Lambda \qquad \forall |\xi| = 1, \qquad \forall x \in \Omega.
\end{equation}

We will let $\rho_{x}$ denote $\rho(x, \cdot)$ to simplify notation. Namely, $(D \rho_{x})(Du) = (D \rho(x, \cdot))(Du)$. We say that $u \in W^{1,\gamma}(\Omega)$ is a subsolution (supersolution) to \eqref{e:inhomo} if

\begin{equation*} 
\int_{\Omega} \left \langle \rho_{x}(Du)^{\gamma-1} (D \rho_{x})(Du), D \varphi \right \rangle \le (\ge) \int_{ \Omega} \left \langle \vec{F}, D \varphi \right \rangle + f \varphi 
\end{equation*}
for all non-negative $\varphi \in W^{1,\gamma^{\prime}}_{0}(\Omega)$. We say that $u$ is a solution if it is both a subsolution and supersolution.



Given a real number, say $\alpha$, in $(1, \infty)$ or $[1,n)$, we will respectively always let $\alpha^{\prime}$ and $\alpha^{*}$ denote the Holder and Sobolev exponents. That is,
$$
\frac{1}{\alpha} + \frac{1}{\alpha^{\prime}} = 1 \qquad \alpha^{*} = \frac{n \alpha}{n-\alpha}.
$$



\begin{theorem}[Gagliado-Nirenberg-Sobolev]
If $u \in W^{1,\gamma}_{0}(\Omega)$ then there exists $C = C(n,\gamma) > 0$ so that
$$
\| u\|_{L^{\gamma^{*}(\Omega)}} \le C \|Du\|_{L^{\gamma}(\Omega)} 
$$
\end{theorem}



\begin{theorem}[Poincare in a ball]
For each $1 \le \gamma < n$ there exists a $C = C(n,\gamma)$ so that
$$
\left( r^{-n} \int_{B(x,r)} |f-(f)_{x,r}|^{\gamma^{*}} dy \right)^{\frac{1}{\gamma^{*}}} \le C_{2} r^{1- \frac{n}{\gamma}} \left( \int_{B(x,r)} |Df|^{\gamma} dy \right)^{\frac{1}{\gamma}}
$$
\end{theorem}

\begin{remark} \label{r:sobolev}
Since $\rho$ will always satisfies \eqref{e:meas}, \eqref{e:p1h}, and \eqref{e:cx}, if it also solves the first inequality in \eqref{e:eb} then the Sobolev embedding theorem can be re-written as 
$$
\|u\|_{L^{\gamma^{*}}(\Omega)} \le C \nu^{-1} \| \rho_{x}(Du)\|_{L^{\gamma}(\Omega)}.
$$
Similarly Poincare in a ball can be re-written as
$$
\left( r^{-n} \int |f-(f)_{x,r}|^{\gamma^{*}} dy \right)^{\frac{1}{\gamma^{*}}} \le C_{2} \nu^{-1} r^{1- \frac{n}{\gamma}} \left( \int_{B(x,r)} \rho_{x}(Df)^{\gamma} dy \right)^{\frac{1}{\gamma}}.
$$
\end{remark}

Given $\rho : \Omega \times \mathbb{R}^n \setminus \{0\} \to (0,\infty)$ as in \eqref{e:meas}, \eqref{e:p1h}, \eqref{e:cx}, define $\rho_{*} : \Omega \times \mathbb{R}^n \setminus \{0\} \to (0, \infty)$ by
$$
\rho_{*}(x, \xi^{*}) = \sup_{\rho(x,\xi) < 1} \xi \cdot \xi^{*}.
$$
That is, $\rho_{*}(x, \cdot)$ is the convex dual of $\rho(x, \cdot)$ for all $x \in \Omega$. We record for the reader's convenience a few facts that are frequently used:
\begin{proposition}
Let $a,b,c, \epsilon > 0$, $\alpha \in (1, \infty)$, $\rho$ as in \eqref{e:meas}, \eqref{e:p1h}, and \eqref{e:cx}. Suppose $\xi_{1}, \xi_{2} \in \mathbb{R}^n \setminus \{0\}$. 
\begin{itemize}
\item Young's inequality says
$$
abc \le a \epsilon^{\alpha} \frac{b^{\alpha}}{\alpha} + a \epsilon^{-\alpha^{\prime}} \frac{c^{\alpha^{\prime}}}{\alpha^{\prime}}.
$$
\item Fenchel's inequality guarantees
$$
\xi_{1} \cdot \xi_{2} \le \rho(\xi_{1}) \rho_{*}(\xi_{2})
$$
\item It holds,
\begin{equation} \label{e:eq1}
\rho_{*}(x,D \rho(\xi_{1})) \equiv 1.
\end{equation}
\item If $\rho$ satisfies \eqref{e:p1h} and \eqref{e:cx} then so does $\rho_{*}$.

\item If $\rho$ satisfies \eqref{e:eb} then for all $|\xi| = 1$,
$$
\Lambda^{-1} \le \rho_{*}(\xi) \le \nu^{-1}.
$$
In particular, $\vec{F} \in L^{q}(\Omega)$ if and only if $\rho_{*}(\cdot,\vec{F}) \in L^{q}(\Omega)$
\end{itemize}
\end{proposition}




A warm-up exercise is the following Cacciopolli type inequality.

\begin{lemma}
If $u$ is a subsolution to \eqref{e:inhomo} with $\vec{F},f \equiv 0$ and $1 < \gamma < \infty$, then
$$
\| \eta \rho_{x}(D u) \|_{L^{\gamma}(\Omega)} \le C(n,\gamma) \| u \rho_{x}(D \eta)\|_{L^{\gamma}(\Omega)}.
$$
\end{lemma}


\begin{proof}
Consider $\varphi = \eta^{\gamma} u$. Then $D \varphi = \gamma \eta^{\gamma-1} u D \eta+ \eta^{\gamma} Du$. So, using the $1$-homogeneity of $\rho$ and Fenchel's inequality, 
\begin{align*}
& \left \langle \rho_{x}(Du)^{\gamma-1} (D \rho_{x})(Du) , D \varphi \right \rangle \\ & \ge- \gamma (\eta \rho_{x}(Du))^{\gamma-1} \rho_{*}(x,(D \rho_{x})(Du)) u \rho_{x}(D \eta) + \eta^{\gamma} \rho_{x}(Du)^{\gamma}.
\end{align*}
Using $u$ is a subsolution with $\vec{F},f \equiv 0$, and \eqref{e:eq1} yields
\begin{align*}
\int_{\Omega} \eta^{\gamma} \rho_{x}(Du)^{\gamma} & \le \gamma \int_{\Omega} ( \eta \rho_{x}(Du))^{\gamma-1} u \rho_{x}(D \eta) \\
& \le \left( \int ( \eta \rho(Du))^{\gamma} \right)^{1 - \frac{1}{\gamma}} \left( \int_{\Omega} u^{\gamma} \rho(D \eta)^{\gamma} \right)^{\frac{1}{\gamma}}.
\end{align*}
Dividing completes the proof.
\end{proof}

Finally, we recall a technical lemma for later use. 
\begin{lemma} \label{l:ibs}
Let $\omega, \sigma$ be non-decreasing functions in $(0,R]$. Suppose there exists $0 < \tau, \tilde \delta < 1$ so that for all $r \le R$,
$$
\omega( \tau r) \le \tilde \delta \omega(r) + \sigma(r).
$$
Then for any $\mu \in (0,1)$ and $r \le R$
$$
\omega(r) \le C \left\{ \left( \frac{r}{R} \right)^{\alpha} \omega(R) + \sigma(r^{\mu} R^{1-\mu}) \right\}
$$
where $C = C(\tilde \delta, \tau)$ and $\alpha = \alpha( \tilde \delta, \tau, \mu)$.
\end{lemma}


\section{Main results} \label{s:theorems}

In this section we focus on functions $u$ that solve
\begin{equation} \label{e:dt1}
\int_{\Omega} \left \langle \rho_{x}(Du)^{\gamma-1} (D \rho_{x})(Du) , D \varphi \right \rangle dx = \int_{\Omega} \left \langle \vec{F}, D \varphi \right \rangle + f \varphi  \qquad \forall \varphi \in W^{1, \gamma^{\prime}}_{0}(\Omega).
\end{equation}
As a corollary of these results, we can answer further questions about functions $u$ that solve 
\begin{equation} \label{e:homo}
\int_{\Omega} \left \langle \rho_{x}(Du)^{\gamma-1} (D \rho_{x})(Du) , D \varphi \right \rangle dx = 0.
\end{equation}

We begin with a Caccioppoli inequality when $\vec{F},f \not\equiv 0$.

\begin{theorem} \label{t:crhs}
Suppose $u \in W^{1, \gamma}(\Omega)$ is a subsolution of \eqref{e:dt1} with $1 < \gamma < \infty$ and $\rho : \Omega \times \mathbb{R}^n \setminus \{0\}$ satisfies \eqref{e:meas}, \eqref{e:p1h}, \eqref{e:cx}, and \eqref{e:eb}. Assume $\vec{F},f \in L^{\tilde \gamma}(\Omega)$ where $\tilde \gamma = \max \{\gamma, \gamma^{\prime} \}$. If $B_{2R} \subset \Omega$ and $0 < R \le 10$ then,
\begin{equation} \label{e:crhs}
\| \rho_{x}(Du)\|_{L^{\gamma}(B_{R})} \le c_{\gamma,\Lambda} \left[ R^{-1} \|u\|_{L^{\gamma}}(B_{2R}) + \| \rho_{*}(x,\vec{F}) \|_{L^{\gamma^{\prime}}(B_{2R})}^{\frac{1}{\gamma-1}} + R^{\frac{1}{\gamma-1}} \| f \|_{L^{\gamma^{\prime}}(B_{2R})}^{\frac{1}{\gamma-1}} \right].
\end{equation}
\end{theorem}

\begin{remark}
Note, it is respectively necessary for $\vec{F}, f \in L^{\gamma}$ ($\vec{F}, f \in L^{\gamma^{\prime}}$) for the equation (conclusion) to make sense\footnote{Technically, by Sobolev embedding we only need $f \in L^{(\frac{\gamma^{\prime} n}{n-\gamma^{\prime}} )^{\prime}} \cap L^{\gamma^{\prime}}$.}. 
\end{remark}

\begin{proof}
Suppose without loss of generality, $R = 1$. Consider the test function $\varphi = \eta^{\gamma} u$ for some $\eta \in C^{1}_{0}(B_{2})$ to be chosen later. Note, 
$$
D \varphi = \eta^{\gamma} Du + \gamma u \eta^{\gamma-1} D \eta.
$$ 
Taking advantage of the $1$-homogeneity of $\rho_{x}$, Fenchel's inequality, Young's inequality, and \eqref{e:eq1} we choose $\epsilon > 0$ so that $(\gamma-1) \epsilon^{\gamma^{\prime}} = \frac{1}{2} \rho( \eta Du)^{\gamma}$ and compute
\begin{align} 
\nonumber \langle \rho_{x}(Du)^{\gamma-1} (D\rho_{x})(Du), D \varphi \rangle & \ge \rho_{x}( \eta Du)^{\gamma} -  \gamma \rho_{x}(\eta Du)^{\gamma-1} \rho_{*}(x,(D \rho)(Du)) \rho_{x}(u D \eta) \\
\nonumber & \ge \rho_{x}( \eta Du)^{\gamma} - \gamma \left[ \frac{ \epsilon^{\gamma^{\prime}} \rho_{x}( \eta D u)^{\gamma}}{\gamma^{\prime}} + \frac{ \rho_{x}(u D \eta)^{\gamma}}{\epsilon^{\gamma} \gamma} \right] \\
\label{e:c1} & \ge \frac{\rho_{x}( \eta Du)^{\gamma}}{2} - c_{\gamma} \rho_{x}(u D \eta)^{\gamma}.
\end{align}
On the other hand, for $\epsilon> 0$ chosen so that $\gamma^{-1} \epsilon^{\gamma} = 1/4$, Fenchel's and Young's inequalities ensure
\begin{align}
\nonumber \langle \vec{F}, D \varphi \rangle & \le \gamma \left(\rho_{*}(x,\vec{F})\eta^{\gamma-1}\right) \rho_{x}(u D \eta) + \eta^{\gamma} \left( \rho_{*}(\vec{F})  \rho(Du) \right) \\
\nonumber & \le \gamma \left[ \frac{\rho_{*}(x,\vec{F})^{\gamma^{\prime}} \eta^{\gamma}}{ \gamma^{\prime}} +  \frac{\rho_{x}(u D \eta)^{\gamma}}{\gamma} \right] + \eta^{\gamma} \left[ \frac{\rho_{*}( x,\vec{F})^{\gamma^{\prime}}}{\epsilon^{\gamma^{\prime}} \gamma^{\prime}} + \epsilon^{\gamma} \frac{ \rho_{x}(Du)^{\gamma}}{\gamma} \right] \\
& \label{e:c2}= \frac{\rho_{x}( \eta Du)^{\gamma}}{4} + c_{\gamma} \rho_{*}(x,\vec{F})^{\gamma^{\prime}} \eta^{\gamma} + \rho_{x}(u D \eta)^{\gamma}.
\end{align}
Combining \eqref{e:c1}, \eqref{e:c2}, and \eqref{e:dt1} yields,
\begin{align} 
\nonumber \int \rho_{x}(\eta Du)^{\gamma} & \le c_{\gamma} \left[ \int \rho_{*}(x,\vec{F})^{\gamma^{\prime}} \eta^{\gamma} + \int \rho_{x}(u D \eta)^{\gamma} + \int f \eta^{\gamma} u \right] \\
\label{e:c4} & \le c_{\gamma} \left[ \int \rho_{*}(x,\vec{F})^{\gamma^{\prime}} \eta^{\gamma} + \int u^{\gamma}(\eta^{\gamma} + \rho_{x}( D \eta)^{\gamma}) +  \int  \eta^{\gamma } f^{\gamma^{\prime}}  \right].
\end{align}
Choosing $0 \le \eta \le 1$, $\eta \equiv 1$ on $B_{1}$, $\eta \equiv 0$ on $B_{2}^{c}$ and $|D \eta| \le 2$ we find
\begin{equation*}
\| \rho_{x}( Du)\|_{L^{\gamma}(B_{R})} \le c_{\gamma, \Lambda,} \left[ \|u\|_{L^{\gamma}(B_{2R})} + \| \rho_{*}(x,\vec{F}) \|_{L^{\gamma^{\prime}}(B_{2R})}^{\gamma^{\prime}-1} + \| f \|_{L^{\gamma^{\prime}}(B_{2R})}^{\gamma^{\prime}-1} \right].
\end{equation*}
Equation \eqref{e:crhs} is recovered by scaling.
\end{proof}



We note that in Theorem \ref{t:crhs}, the fact that $\rho$ can depend on $x$ never needs to be dealt with separately. This is unsurprising because conditions \eqref{e:p1h}, \eqref{e:cx}, and Fenchel's inequality are used at a pointwise level while \eqref{e:eb} are used at a global level, see \ref{r:sobolev}. Hence, to further simplify the notation, we only explicitly write-out the dependence of $\rho$ on $x$ in the statements of theorems and suppress this dependence on $x$ throughout the proofs.


\begin{theorem} \label{t:sup}
Suppose $\rho : \Omega \times \mathbb{R}^n \setminus \{0\}$ satisfies \eqref{e:meas}, \eqref{e:p1h}, \eqref{e:cx},  and \eqref{e:eb}. If $u$ is a subsolution to \eqref{e:dt1}, $\vec{F} \in L^{q}(\Omega)$ and $f \in L^{\frac{q}{\gamma^{\prime}}}(\Omega)$ for some $q > \frac{n}{\gamma-1}$, $1 < \gamma < n$, $0 < r < R < 1$, and $x \in \Omega$ with $\overline{B_{R}} \subset \Omega$ then there exists some $C = C(n,\nu, \Lambda,\gamma,q,p)$ and $\delta = 1 - \frac{n}{q(\gamma-1)} > 0$ so that 
\begin{align*}
\sup_{B_{r}} u^{+} & \le C  \bigg[ \frac{\| u^{+}\|_{L^{p}(B_{R})}}{(R-r)^{\frac{n}{p}} } + R^{\delta} \|\rho_{*}(x,\vec{F})\|_{L^{q}(B_{R})}^{\frac{1}{\gamma-1}} +  R^{\gamma^{\prime} \delta}\|f\|_{L^{\frac{q}{\gamma^{\prime}}}(B_{R})}^{\frac{1}{\gamma-1}} \bigg]
\end{align*}
\end{theorem}

\begin{proof}
We consider the test function $\varphi = \eta^{\gamma} v^{\beta} \bar{u}$ for $\beta \ge 0$ where $\bar{u} = u^{+} + k$ and $v = \min\{u,m\}$ for some $0 < k < m < \infty$, $k$ to be chosen later. Notice 
$$
D \varphi = \gamma \eta^{\gamma-1} v^{\beta} \bar{u} D \eta + \beta \eta^{\gamma} v^{\beta} \bar{U} Dv + \eta^{\gamma} v^{\beta} D \bar{u}.$$

We wish to expand out \eqref{e:dt1} with this choice of $\varphi$. To this end, first observe $1$-homogeneity, i.e., \eqref{e:p1h} ensures
\begin{align} \label{e:dt2}
\nonumber \big \langle \rho(Du)^{\gamma-1} (D \rho)(Du), & \beta v^{\beta-1} \bar{u} \eta^{\gamma} D v + \eta^{\gamma} v^{\beta} D \bar{u} \big \rangle \\
& = \beta \rho(Dv)^{\gamma} v^{\beta} \eta^{\gamma} + \rho(D\bar{u})^{\gamma} \eta^{\gamma} v^{\beta}.
\end{align}
Next we apply Fenchel's and Young's inequalities in combination with \eqref{e:eq1} for some $\epsilon = \epsilon(\gamma) > 0$ to be chosen immediately after \eqref{e:dt3},
\begin{align} 
\nonumber \big \langle \rho(Du)^{\gamma-1} (D \rho)(Du),&  \gamma \eta^{\gamma-1} v^{\beta} \bar{u} D \eta \big \rangle  \\
\nonumber & \ge - \gamma v^{\beta} \left( \rho(D \bar{u}) \eta) \right)^{\gamma-1} \rho_{*}(D\rho(Du)) \left( \rho(D \eta) \bar{u} \right)      \\
\label{e:dt3}  & \ge - \gamma v^{\beta} \left[  \epsilon^{\gamma^{\prime}} \frac{\rho(D \bar{u})^{\gamma} \eta^{\gamma}}{\gamma^{\prime}} + \frac{ \rho(D \eta)^{\gamma} \bar{u}^{\gamma}}{ \epsilon ^{\gamma} \gamma} \right].
\end{align}
Choose $\epsilon$ so that $ (\gamma-1) \epsilon^{\gamma^{\prime}} = 1/2$. Since $\frac{\gamma}{\gamma^{\prime}} = \gamma-1$, this choice of $\epsilon$ ensures \eqref{e:dt3} becomes
\begin{align} \label{e:dt4}
 \left \langle \rho(Du)^{\gamma-1} (D \rho)(Du), \gamma \eta^{\gamma-1} v^{\beta} \bar{u} D \eta \right \rangle & \ge  -  \frac{v^{\beta} \rho(D \bar{u})^{\gamma} \eta^{\gamma}}{2} - c_{\gamma} v^{\beta} \rho(D\eta)^{\gamma} \bar{u}^{\gamma}, 
\end{align}
where $c_{\gamma}$ may change depending on the line, but depends only on $\gamma$.

Now we look at the righthand side. We split this into two pieces and treat the first piece in much the same fashion as above. 
\begin{align}
\nonumber \big\langle \vec{F}, & \beta v^{\beta-1} \bar{u} \eta^{\gamma} Dv + \eta^{\gamma} v^{\beta} D \bar{u} \big\rangle  \le  \rho_{*}(\vec{F}) \left[ \beta v^{\beta} \eta^{\gamma} \rho(Dv) \right] + \rho_{*}(\vec{F}) \left[ \eta^{\gamma} v^{\beta} \rho(D \bar{u}) \right] \\
\nonumber & = (\eta^{\gamma} v^{\beta}) \left[ \left(\beta\right) \left( \rho_{*}(\vec{F}) \right) \left(  \rho(Dv) \right) + \rho_{*}(\vec{F}) \rho(D \bar{u}) \right] \\
\nonumber & \le   \eta^{\gamma} v^{\beta} \left[ \left( \beta \frac{\rho_{*}(\vec{F})^{\gamma^{\prime}}}{\epsilon_{1}^{\gamma^{\prime}} \gamma^{\prime}} + \beta \epsilon_{1}^{\gamma} \frac{ \rho(D v)^{\gamma}}{\gamma} \right) + \left( \frac{ \rho_{*}(\vec{F})^{\gamma^{\prime}}}{ \epsilon_{2}^{\gamma^{\prime}} \gamma^{\prime}} + \epsilon_{2}^{\gamma} \frac{ \rho(D \bar{u})^{\gamma}}{\gamma} \right) \right] \\
& \label{e:dt5} =  \eta^{\gamma} v^{\beta} \rho_{*}(\vec{F})^{\gamma^{\prime}}\left( \frac{\beta}{\epsilon_{1}^{\gamma^{\prime}} \gamma^{\prime}} + \frac{1}{\epsilon_{2}^{\gamma^{\prime}} \gamma^{\prime}} \right) + \frac{\beta \epsilon_{1}^{\gamma} }{\gamma} \eta^{\gamma} v^{\beta} \rho(Dv)^{\gamma} + \frac{\epsilon_{2}^{\gamma}}{2} \eta^{\gamma} v^{\beta} \rho(D \bar{u})^{\gamma}.
\end{align}
Since we want to absorb the last two terms of \eqref{e:dt5} into \eqref{e:dt2} we choose $\epsilon_1, \epsilon_{2}$ so that $\gamma^{-1} \beta \epsilon_{1}^{\gamma} = \frac{\beta}{2}$ and $\gamma^{-1} \epsilon_{2}^{\gamma} = \frac{1}{4}$. The need for choosing $1/4$ for the $\epsilon_{2}$ coefficient is due to the fact that we'll also be absorbing the $\rho(D \bar{u})$-term from \eqref{e:dt4} into \eqref{e:dt2}. We note both $\epsilon_{1}$ and $\epsilon_{2}$ depend solely on $\gamma$. All-in-all this allows us to re-write \eqref{e:dt5} as
\begin{align} \label{e:dt6}
 \left \langle  \vec{F}, \beta v^{\beta-1} \bar{u} \eta^{\gamma} Dv + \eta^{\gamma} v^{\beta} D \bar{u} \right \rangle & \le  c_{\gamma}(1 + \beta) \eta^{\gamma} v^{\beta} \rho_{*}(\vec{F})^{\gamma^{\prime}} + \frac{\beta}{2} \eta^{\gamma} v^{\beta} \rho(Dv)^{\gamma} \\ 
 \nonumber & + \frac{1}{4} \eta^{\gamma} v^{\beta} \rho(D \bar{u})^{\gamma}.
\end{align}

We now deal with the final term via Fenchel and Young's inequalities
\begin{align}
\nonumber \left \langle \vec{F}, \gamma \eta^{\gamma-1} v^{\beta} \bar{u} D \eta \right \rangle & \le  \gamma \rho_{*}(\vec{F}) \eta^{\gamma-1} v^{\beta} \bar{u} \rho(D \eta) \\
\nonumber & \le  \gamma v^{\beta} \left[ \frac{ \rho_{*}(\vec{F})^{\gamma^{\prime}} \eta^{\gamma} }{\gamma^{\prime}} + \frac{ \bar{u}^{\gamma} \rho(D \eta)^{\gamma}}{\gamma} \right] \\
\label{e:dt7} & =  (\gamma-1) \rho_{*}(\vec{F})^{\gamma^{\prime}} \eta^{\gamma} v^{\beta} + \bar{u}^{\gamma} \rho(D \eta)^{\gamma} v^{\beta}.
\end{align}

Finally, plugging \eqref{e:dt2}, \eqref{e:dt4}, \eqref{e:dt6}, and \eqref{e:dt7} into \eqref{e:dt1} yields  
\begin{align} 
\nonumber \frac{\beta}{2}& \int_{\Omega} \rho(Dv)^{\gamma} v^{\beta} \eta^{\gamma} + \frac{1}{4} \int_{\Omega} \rho(D \bar{u})^{\gamma} \eta^{\gamma} v^{\beta} \\
\nonumber &  \le  c_{\gamma} \left[ \int_{\Omega} v^{\beta} \rho(D \eta)^{\gamma} \bar{u}^{\gamma} + \int_{\Omega} \eta^{\gamma} v^{\beta} \rho_{*}(\vec{F})^{\gamma^{\prime}} + \int_{\Omega} f \eta^{\gamma} v^{\beta} \bar{u} \right] \\
 \label{e:dt8} & \le c_{\gamma} \left[ \int_{\Omega} v^{\beta} \rho(D \eta)^{\gamma} \bar{u}^{\gamma} + (1+\beta)\int_{\Omega} \eta^{\gamma} v^{\beta} \frac{\bar{u}^{\gamma}}{k^{\gamma}} \rho_{*}(\vec{F})^{\gamma^{\prime}} + \int_{\Omega}  \frac{f}{k^{\gamma-1}} \eta^{\gamma} v^{\beta} \bar{u}^{\gamma} \right].
\end{align}

The final inequality used $\bar{u} \ge k$. Set $w= v^{\beta/\gamma} \bar{u}$, we note
$$
Dw = \frac{\beta}{\gamma} v^{\frac{\beta}{\gamma}-1} \bar{u} Dv + v^{\frac{\beta}{\gamma}} D \bar{u} = \frac{\beta}{\gamma} v^{ \frac{\beta}{\gamma}} Dv + v^{\frac{\beta}{\gamma}} D \bar{u}.
$$
Since $\rho(\xi_{1} + \xi_{2}) \le \rho(\xi_{1}) + \rho(\xi_{2})$ for all $\xi_{1}, \xi_{2}$ it follows that
\begin{align*}
\eta^{\gamma} \rho(D w)^{\gamma} &\le \eta^{\gamma} \left( \frac{\beta}{\gamma} v^{\frac{\beta}{\gamma}} \rho(Dv) + v^{\frac{\beta}{\gamma}} \rho(D \bar{u}) \right)^{\gamma} \\
& \le 2^{\gamma-1} \left( \left( \frac{\beta}{\gamma} \right)^{\gamma} \eta^{\gamma} v^{\beta} \rho(Dv)^{\gamma} + \eta^{\gamma} v^{\beta} \rho(D \bar{u})^{\gamma} \right).
\end{align*}
In particular, this guarantees that for some $c_{\gamma}$
\begin{equation} \label{e:dt9}
\int_{\Omega} \eta^{\gamma} \rho(Dw)^{\gamma} \le c_{\gamma} (1 + \beta^{\gamma-1}) \left[ \frac{\beta}{2} \int_{\Omega} \rho(Dv)^{\gamma} v^{\beta} \eta^{\gamma} + \frac{1}{4} \int_{\Omega} \rho(D \bar{u})^{\gamma} \eta^{\gamma} v^{\beta} \right].
\end{equation}
Combining \eqref{e:dt8} and \eqref{e:dt9} yields
\begin{equation} \label{e:dt10} 
\int \eta^{\gamma} \rho(Dw)^{\gamma} \le c_{\gamma}(1 + \beta^{\gamma}) \left[ \int_{\Omega} w^{\gamma} \rho(D \eta)^{\gamma}  +\int_{\Omega} (\eta w)^{\gamma} \left( \frac{\rho_{*}(\vec{F})^{\gamma^{\prime}}}{k^{\gamma}} + \frac{f}{k^{\gamma-1}} \right) \right].
\end{equation}
Due to the observation that $\rho(D(\eta w))^{\gamma} \le 2^{\gamma-1} (w^{\gamma} \rho(D \eta)^{\gamma} + \eta^{\gamma} \rho(D w)^{\gamma})$, \eqref{e:dt10} implies
\begin{equation} \label{e:dt10.1}
\int_{\Omega} \rho(D( \eta w))^{\gamma} \le c_{\gamma}(1 + \beta^{\gamma}) \left[ \int_{\Omega} w^{\gamma} \rho(D \eta)^{\gamma} + \int_{\Omega} (\eta w)^{\gamma} \left(
\frac{\rho_{*}(\vec{F})^{\gamma^{\prime}}}{k^{\gamma}} + \frac{f}{k^{\gamma-1}} \right) \right].
\end{equation}

Our next goal is to deal with the term $\int (\eta w)^{\gamma} \frac{\rho_{*} (\vec{F})^{\gamma^{\prime}}}{k^{\gamma}}$. We roughly explain in words how we do this. We  first use Holder's inequality to make the $L^{q}$ norm of $\vec{F}$ appear. To this end, we will introduce the parameter $\alpha_{1} = \alpha(q,\gamma) = \frac{q}{q-\gamma^{\prime}}$ which is equivalent to $\alpha_{1}^{\prime} = \frac{q}{\gamma^{\prime}}$. This is precisely where the requirement $f \in L^{\frac{q}{\gamma^{\prime}}}$ comes from.

Next, by choosing $k$ appropriately, we will make the term with $\rho(\vec{F}) + f$ be absorbed into a $1$. At this point, the remaining term with $\eta w$ will have a strange power. So, we use interpolation, and the fact that $ \gamma < \gamma \alpha_{1} < \frac{n \gamma}{n- \gamma}$, to re-write our strange power as a linear combination of the $L^{\gamma}$ and $L^{\gamma^{*}}$ norms of $\eta w$.  The necessary upper-bound on $\alpha_{1}$ is satisfied so long as $q > \frac{n}{\gamma-1}$. 

By making the coefficient of the $L^{\gamma}$ norm of $\eta w$ larger, we can choose the $L^{\gamma^{*}}$ norm of $\eta w$ to be arbitrarily small. This is necessary, as we will apply the Gagliardo-Nirenberg-Sobolev inequality to turn this latter norm into an estimate on the $L^{\gamma}$ norm of $\rho(D (\eta w))$ which we can finally absorb into the left hand side. When we apply Young's inequality in order to make this weighted-linear combination of norms appear, we will choose $\alpha_{2} = \alpha(n,q,\gamma) = \theta_{1}^{-1}$, where $\theta_{1}$ is the interpolation power. This conveniently makes all exponents outside of integrals disappear, allowing the desired simplifications to all occur.

We begin the process outlined above by applying Holder's inequality,
\begin{align} \nonumber
\int_{\Omega}(\eta w)^{\gamma} \left(\frac{\rho_{*}(\vec{F})^{\gamma^{\prime}}}{k^{\gamma}} + \frac{f}{k^{\gamma-1}} \right)& \le \| (\eta \omega)^{\gamma}\|_{L^{\alpha_{1}}(\Omega)}  \left\| \frac{\rho_{*}(\vec{F})^{\gamma^{\prime}}}{k^{\gamma}} +\frac{ f }{k^{\gamma-1}} \right\|_{L^{\alpha_{1}^{\prime}}(\Omega)} \\
\label{e:dt11} & \le C(\gamma,q) \| \eta w \|_{L^{\gamma \alpha_{1}}(\Omega)}^{\gamma}
\end{align}
where $\alpha_{1}^{\prime} = \frac{q}{\gamma^{\prime}}$ is as above, and $k$ is chosen so that $k = k_{1} + k_{2}$ where $k_{1}^{\gamma} = \|\rho_{*}(\vec{F})\|_{L^{q}(\Omega)}^{\gamma^{\prime}}$ and $k_{2}^{\gamma-1} = \|f\|_{L^{\frac{q}{\gamma^{\prime}}}(\Omega)}$. If $\vec{F},f \equiv 0$ choose $k > 0$ arbitrary, and you can later send $k$ to zero. Next, we define $\theta_{1} = \theta_{1}(q, n,\gamma) \in (0,1)$ so that $\frac{1}{\gamma \alpha_{1}} = \frac{\theta_{1}}{\gamma} + \frac{(1- \theta_{1})(n-\gamma)}{n \gamma}$. Note that if $\alpha_{2} = \alpha(q,n,\gamma) = \theta_{1}^{-1}$ then $\alpha_{2}^{\prime} = (1-\theta_{1})^{-1}$.  Riesz-Thorin interpolation applied to \eqref{e:dt11}, Young's inequality, and the Gagliardo-Nirenberg-Sobolev inequality consecutively ensure
\begin{align} \nonumber
\int_{\Omega}(\eta w)^{\gamma} \frac{\rho_{*}(\vec{F})^{\gamma^{\prime}}}{k^{\gamma}}& \le \left( \| \eta w\|_{L^{\gamma}(\Omega)}^{\theta_{1}} \| \eta w\|_{L^{\gamma^{*}}(\Omega)}^{1-\theta}\right)^{\gamma} \\
& \nonumber \le  \frac{1}{\epsilon^{\alpha_{2}} \alpha_{2}} \|\eta w\|_{L^{\gamma}(\Omega)}^{\gamma \theta_{1} \alpha_{2}} + \frac{\epsilon^{\alpha_{2}^{\prime}}}{\alpha_{2}^{\prime} } \|\eta w\|_{L^{\gamma^{*}}(\Omega)}^{\gamma(1 - \theta_{1}) \alpha_{2}^{\prime}} \\
& \nonumber =  \frac{1}{\epsilon^{\alpha_{2}} \alpha_{2}} \|\eta w\|_{L^{\gamma}(\Omega)}^{\gamma} + \frac{\epsilon^{\alpha_{2}^{\prime}}}{\alpha_{2}^{\prime} } \|\eta w\|_{L^{\gamma^{*}}(\Omega)}^{\gamma} \\
\label{e:dt11.1} & \le \epsilon^{- \alpha_{2}} \alpha_{2}^{-1} \| \eta w \|_{L^{\gamma}(\Omega)}^{\gamma} + C(q,n,\gamma,\rho) \epsilon^{\alpha_{2}^{\prime}}  \| \rho(D(\eta w)) \|_{L^{\gamma}(\Omega)}^{\gamma}.
\end{align}

See Remark \ref{r:sobolev} for our non-standard application of Gagliardo-Nirenberg-Sobolev inequality. Next, we want to plug \eqref{e:dt11.1} into \eqref{e:dt10.1} and subtract over the $\rho(D(\eta w))$ term, so we choose $\epsilon$ so that the coefficient of $\| \rho(D(\eta w))\|_{L^{\gamma}(\Omega)}^{\gamma}$ is $1/2$. That is, choose $\epsilon = \epsilon(q,n,\gamma,\rho) >0$ by 
$$
c_{\gamma}(1+\beta^{\gamma}) C(q,n,\gamma,\nu) \epsilon^{\alpha^{\prime}_{2}} = \frac{1}{2}. 
$$
Then using our choice of $\epsilon$ and plugging \eqref{e:dt11.1} into \eqref{e:dt10.1} yields
\begin{align} \nonumber
\int_{\Omega} \rho(D(\eta w))^{\gamma} & \le c_{\gamma}(1 + \beta^{\gamma}) \left[ \int_{\Omega} w^{\gamma} \rho(D \eta)^{\gamma} + C_{q,n,\gamma,\nu} (1 + \beta^{\gamma})^{\alpha_{2}-1} \int_{\Omega} w^{\gamma} \eta^{\gamma} \right] \\
 \label{e:dt12} & \le  C_{q,n,\gamma,\nu} (1 + \beta^{\gamma})^{\alpha_{2}} \left[ \int_{\Omega} w^{\gamma} (\rho(D \eta)^{\gamma} + \eta^{\gamma}) \right],
\end{align}
where as always, $\alpha_{2} = \alpha(n,q,\gamma) > 0$. Finally, the Gagliardo-Nirenberg-Sobolev inequality applied to \eqref{e:dt12} implies
\begin{equation} \label{e:dt13}
\| \eta w \|_{\gamma \chi} \le C_{q,n,\gamma,\nu}(1 + \beta^{\gamma})^{\frac{\alpha_{2}}{\gamma}} \left[ \int_{\Omega} w^{\gamma}(\rho(D \eta)^{\gamma} + \eta^{\gamma}) \right]^{\frac{1}{\gamma}},
\end{equation}
where $\chi = \frac{n}{n-\gamma}$. Choose the cut-off function $\eta$ so that with $0 < r < R <1$ and some $B_{R}(x) \subset \subset \Omega$, $\eta \in C^{1}_{0}(B_{R})$ and
\begin{equation} \label{e:dt13.1}
\eta \equiv 1 \text{ in } B_{r} \quad \text{and} \quad |D \eta| \le \frac{2}{R-r}.
\end{equation}


Then \eqref{e:dt13} guarantees
\begin{equation*} 
\left(\int_{B_{r}(x)} w^{\gamma \chi}\right)^{\frac{1}{\chi}} \le C_{n,\gamma,q, \nu, \Lambda} \frac{(1+ \beta^{\gamma-1})^{\alpha_{2}}}{(R-r)^{\gamma}} \int_{B_{R}} w^{\gamma}
\end{equation*}
where our constant gained a dependence on $\Lambda$, which arises by applying \eqref{e:dt13.1} in the form $\rho(D \eta) \le \Lambda |D \eta| \le  2\Lambda(R-r)^{-1}$. Recall the definition of $w$ and use $v \le \bar{u}$ to discover
\begin{equation} \label{e:dt14} 
\left(\int_{B_{r}(x)} v^{(\beta + \gamma)\chi}  \right)^{\frac{1}{\chi}} \le C_{n,\gamma,q, \nu, \Lambda } \frac{(1+ \beta^{\gamma-1})^{\alpha_{2}}}{(R-r)^{\gamma}} \int_{B_{R}} \bar{u}^{\beta + \gamma}.
\end{equation}
 Taking $m \uparrow \infty$ in \eqref{e:dt14} yields
\begin{equation} \label{e:dt15}
\| \bar{u} \|_{L^{(\beta + \gamma)\chi}B_{r}(x)}  \le \left( C_{n,\gamma,q, \nu , \Lambda } \frac{(1+ \beta^{\gamma-1})^{\alpha_{2}}}{(R-r)^{\gamma}}\right)^{\frac{1}{\beta + \gamma}} \| \bar{u}\|_{L^{\beta + \gamma}(B_{R}(x))},
\end{equation}
so long as the right hand side is finite. Note the constant on the righthand side is independent of $\beta$, suggesting we iterate. We begin with $\beta = \beta_{0} = 0$ and for $k =1,2,3 \dots$ define $\beta_{k} = \gamma( \chi^{k-1}-1)$. For $k = 0, 1,2, \dots$ we consider $r_{k} = r + \frac{R-r}{2^{k+1}}$. Note, $(\beta_{k} + \gamma)\chi = \gamma \chi^{k} = (\beta_{k+1}+\gamma)$ and $r_{k-1} - r_{k} = 2^{-(k+1)} (R-r)$.

Writing $C = C_{n,\gamma,q, \nu, \Lambda}$, for $k \ge 0$ this choice of $\beta_{k}, r_{k}$ \eqref{e:dt15} reads
\begin{align*} 
\| \bar{u}\|_{\displaystyle L^{\beta_{k} + \gamma}(B_{r_{k}}(x))} & \le \left( C \frac{ (1 + \beta_{k}^{\gamma-1})^{\alpha_{2}}}{(R-r)^{\gamma} } \right)^{\frac{1}{\beta_{k}+\gamma}} 2^{\frac{k+1}{\gamma \chi^{k-1}}} \| \bar{u} \|_{\displaystyle  L^{\beta_{k-1} + \gamma}(B_{r_{k-1}}(x))}.
\end{align*}
Recalling $\beta_{k} = \gamma(\chi^{k-1}-1)$ we observe that for $k \ge 2$,
\begin{align*} 
\left( \frac{ (1 + \beta_{k}^{\gamma-1})^{\alpha_{2}}}{(R-r)^{\gamma} } \right)^{\frac{1}{\beta_{k}+\gamma}} & \le (R-r)^{- \chi^{k}} (C (1+\gamma \chi))^{\frac{k (\gamma-1) \alpha_{2}}{\gamma \chi^{k-1} -1}} \\
& \le (R-r)^{-\chi^{k}} C^{ \frac{k}{\chi^{k}}}.
\end{align*}
Consequently, after iterating out the first few cases by hand, for all $k$,
$$
\| \bar{u}\|_{\displaystyle L^{\beta_{k} + \gamma}(B_{r_{k}}(x))} \le (R-r)^{-\sum_{i=0}^{\infty} \chi^{i}}C^{ \sum_{i=0}^{\infty} \frac{i}{\chi^{i-1}}} \| \bar{u}\|_{L^{\beta_{0}+\gamma }(B_{r_{1}}(x))},
$$
where still $C= C(n,\gamma, q, \nu, \Lambda)$. Since $u \in W^{1,\gamma}(\Omega)$ and $\beta_{0} = 0$, the right hand side is seen to be finite, and is independent of $k$. Taking $k \to \infty$ on the left hand side, and recalling $\bar{u} = u^{+} + k$ yields 
$$
 \| u^{+} \|_{L^{\infty}(B_{r}(x))} \le C (R-r)^{-\frac{n}{\gamma}}\left[ \|u^{+}\|_{L^{\gamma}(B_{R})} + k \right].
$$
Recalling $k= \| \rho_{*}(\vec{F})\|_{q}^{\frac{\gamma^{\prime}}{\gamma}} + \|f\|_{\frac{q}{\gamma^{\prime}}}^{\frac{1}{\gamma-1}}$ and noting $\frac{1}{1- \chi^{-1}} = \frac{n}{\gamma}$ yields the result for the case $p \ge \gamma$. A classical scaling argument covers the case $0 < p < \gamma$. See, for instance, \cite[p. 75]{han2011elliptic}. 
\end{proof}

\begin{theorem} \label{t:inf}
Let $\rho : \Omega \times \mathbb{R}^n \setminus \{0\}$ satisfy \eqref{e:meas}, \eqref{e:p1h}, \eqref{e:cx}, and \eqref{e:eb}. Suppose $u \in W^{1,\gamma}(\Omega)$ is a supersolution to \eqref{e:dt1} in $\Omega$ and $\vec{F} \in L^{q}_{\operatorname{loc}}(\Omega)$, $f \in L^{\frac{q}{\gamma^{\prime}}}_{\operatorname{loc}}(\Omega)$ for some $q > \frac{n}{\gamma-1}$. Suppose that $B_{2 R} \subset \Omega$ and $0 < p < \frac{(\gamma-1) n}{n- \gamma}$. Then for any $0 < \theta < \tau < 1$, there exists $C = C(n,\gamma,q,p,\nu,\Lambda,\theta,\tau)$ and $\delta = 1 - \frac{n}{q(\gamma-1)} > 0$ so that
$$
\inf_{B_{\theta R}} u^{+} + R^{\delta} \|\rho_{*}(\vec{F})\|_{L^{q}(B_{R})}^{\frac{1}{\gamma-1}} + R^{\gamma^{\prime}\delta} \| f\|_{L^{\frac{q}{\gamma^{\prime}}}(B_{R})}^{\frac{1}{\gamma-1}} \ge C R^{-\frac{n}{p}}\| u^{+}\|_{L^{p}(B_{\tau R})}
$$
and the dependence on $\rho$ only appears as a dependence on $\sup_{|\xi|=1} \rho(\xi)$ and $\inf_{|\xi|=1} \rho(\xi)$.
\end{theorem}


An immediate corollary of Theorem \ref{t:sup} and Theorem \ref{t:inf} is the weak-Harnack inequality.

\begin{theorem} \label{t:iwh}
Suppose $\rho : \Omega \times \mathbb{R}^n \setminus \{0\} \to (0, \infty)$ satisfies \eqref{e:meas}, \eqref{e:p1h}, \eqref{e:cx}, and \eqref{e:eb}. If $1 < \gamma < n$ and $u \in W^{1,\gamma}(\Omega)$ is a nonnegative solution to \eqref{e:dt1} in $\Omega$ for some $\vec{F} \in L^{q}_{\operatorname{loc}}(\Omega)$ and $f \in L^{\frac{q}{\gamma^{\prime}}}_{\operatorname{loc}}(\Omega)$, $q > \frac{n}{\gamma-1}$ and $\overline{B_{3R}} \subset \Omega$, then there exists some $C = C(n, \gamma, q, \nu, \Lambda)$ and $\delta = 1 - \frac{n}{q(\gamma-1)} > 0$ so that
\begin{equation} \label{e:iwh}
\sup_{B_{R}} u \le C_{n, \gamma,\nu, \Lambda,q} \left[ \inf_{B_{2R}} u+ R^{\delta} \|\rho_{*}(\vec{F})\|_{L^{q}(B_{3R})}^{\frac{1}{\gamma-1}} + R^{\gamma^{\prime}\delta} \| f\|_{L^{\frac{q}{\gamma^{\prime}}}(B_{3R})}^{\frac{1}{\gamma-1}} \right] .
\end{equation}
\end{theorem}

The corollary follows readily from Theorem \ref{t:sup} and Theorem \ref{t:inf} by choosing, for instance, $p =  \frac{(\gamma-1)n}{2(n-\gamma)}$ in both theorems. In the case that $\vec{F}, f \equiv 0$ this recovers an inequality of the same form as the classical Harnack inequality. 


\begin{proof}

Let $k > 0$ to be chosen later, $\bar{u} = u^{+} + k$ , $v = \min\{ \bar{u}, m\}$. Consider the test function $\varphi = \eta^{\gamma} v^{-\beta} \bar{u}$ for $\beta \ge \beta_{p} > 1$, where $\beta_{p} = \gamma - \frac{p}{\chi}$, $\chi = \frac{n}{n - \gamma}$.

We now proceed in a similar fashion to the proof of Theorem \ref{t:sup}. Observe,
\begin{align}
\nonumber \big \langle \rho(Du)^{\gamma-1} (D \rho)(Du),& -\beta v^{-\beta-1} \bar{u} \eta^{\gamma} D v + \eta^{\gamma} v^{-\beta} D \bar{u} \big \rangle  \\
\nonumber & = -\beta \rho(Dv)^{\gamma} v^{-\beta} \eta^{\gamma} + \rho(D\bar{u})^{\gamma} \eta^{\gamma} v^{-\beta} \\
\label{e:lb0} & = (1- \beta) \rho(Dv)^{\gamma} v^{-\beta} \eta^{\gamma}.
\end{align}
We note that $(\beta-1) \ge 1 - \beta_{p}$; this observation will simplify our computations later as, when we iterate $\beta$, we may now keep constants independent of $\beta$. They will instead depend on $\beta_{p}$, which depends on $n, \gamma, p$. 

Analogous to \eqref{e:dt3} we compute

\begin{equation} \label{e:lb1}
\left \langle \rho(Du)^{\gamma-1}(D \rho)(Du), \gamma \eta^{\gamma-1} v^{-\beta} \bar{u} D \eta \right \rangle \le \gamma v^{-\beta} \left[ \epsilon^{\gamma^{\prime}} \frac{ \rho(D \bar{u})^{\gamma} \eta^{\gamma}}{\gamma^{\prime}} + \frac{\rho(D \eta)^{\gamma} \bar{u}^{\gamma}}{\epsilon^{\gamma} \gamma}\right].
\end{equation}
Choosing $\epsilon$ so that $(\gamma-1)\epsilon^{\gamma^{\prime}} = -(1- \beta)/2 >(1 - \gamma)/2 > 0$ yields
\begin{align} 
\nonumber \bigg \langle \rho(Du)^{\gamma-1}(D \rho)(Du), & \gamma \eta^{\gamma-1} v^{-\beta} \bar{u} D \eta \bigg \rangle \\ 
\label{e:lb2} & \le  \frac{- (1-\beta) }{2} v^{-\beta}\eta^{\gamma}  \rho(D v)^{\gamma} + c_{n,\gamma,p} v^{-\beta} \rho(D \eta)^{\gamma} \bar{u}^{\gamma}.
\end{align}
Next we treat the $\vec{F}$ term. Proceeding as in \eqref{e:dt5} we discover
\begin{align} 
\nonumber \big\langle \vec{F}, &  -\beta v^{-\beta-1} \bar{u} \eta^{\gamma} Dv + \eta^{\gamma} v^{\beta} D \bar{u} \big\rangle  \ge - (\beta-1)  \eta^{\gamma}v^{-\beta} \left[  \left(\rho_{*}(\vec{F}) \right) \left(  \rho(Dv) \right) \right] \\
\nonumber & \ge (1- \beta) \left[ \frac{ \eta^{\gamma} v^{-\beta} \rho_{*}(\vec{F})^{\gamma^{\prime}}}{\epsilon^{\gamma^{\prime}} \gamma^{\prime}} + \epsilon^{\gamma} \frac{ \eta^{\gamma} v^{-\beta} \rho(Dv)^{\gamma}}{\gamma} \right] \\
\label{e:lb3} & \ge  \frac{(1-\beta)}{4} \eta^{\gamma} v^{-\beta} \rho(Dv)^{\gamma} -  c_{n,\gamma,p}  \eta^{\gamma} v^{-\beta} \rho_{*}(\vec{F})^{\gamma^{\prime}}
\end{align}
where we chose $\epsilon > 0$ so that $\gamma^{-1} \epsilon^{\gamma}= 1/4$. To bound the final piece, we compute

\begin{align}
\nonumber \left \langle \vec{F}, \gamma \eta^{\gamma-1} v^{-\beta} \bar{u} D \eta \right \rangle & \ge - \gamma \rho_{*}(\vec{F}) \eta^{\gamma-1} v^{-\beta} \bar{u} \rho(D \eta) \\
\nonumber & \ge -  \gamma v^{-\beta} \left[ \frac{ \rho_{*}(\vec{F})^{\gamma^{\prime}} \eta^{\gamma} }{\gamma^{\prime}} + \frac{ \bar{u}^{\gamma} \rho(D \eta)^{\gamma}}{\gamma} \right] \\
\label{e:lb4} & = - (\gamma-1) \rho_{*}(\vec{F})^{\gamma^{\prime}} \eta^{\gamma} v^{-\beta} - \bar{u}^{\gamma} \rho(D \eta)^{\gamma} v^{-\beta}.
\end{align}

Using that $u$ is a supersolution and plugging \eqref{e:lb0}, \eqref{e:lb2}, \eqref{e:lb3}, and \eqref{e:lb4} into \eqref{e:dt1} achieves
\begin{align*}
& \frac{-(\beta-1)}{4}  \int_{\Omega} \rho(Dv)^{\gamma} v^{-\beta} \eta^{\gamma}\\
 & \ge - c_{n,\gamma,p} \int_{\Omega} v^{-\beta} \rho(D \eta)^{\gamma} \bar{u}^{\gamma} - c_{n,\gamma,p} \int_{\Omega} \eta^{\gamma} v^{-\beta} \rho_{*}(\vec{F})^{\gamma^{\prime}} - \int_{\Omega} \eta^{\gamma} v^{-\beta+1} f
\end{align*}
or after recalling $\bar{u} \ge k$, $\beta-1 \ge \beta_{p} -1 >0$ and multiplying by $-1$,
\begin{align} \label{e:lb5}
\nonumber \int_{\Omega} \rho(Dv)^{\gamma} v^{-\beta} \eta^{\gamma} &\le c_{n,\gamma,p} \bigg[ \int_{\Omega} v^{-\beta} \bar{u}^{\gamma} \rho(D \eta)^{\gamma} + \\
& + \int_{\Omega} \eta^{\gamma} v^{-\beta} \bar{u}^{\gamma} \frac{\rho_{*}(\vec{F})^{\gamma^{\prime}}}{k^{\gamma}} + \int_{\Omega} \eta^{\gamma} v^{-\beta} \bar{u}^{\gamma} \frac{f}{k^{\gamma-1}} \bigg].
\end{align}

Note the striking similarity to \eqref{e:dt8}. The main difference being the benefits that \eqref{e:lb5} has no constants depending on $\beta$, which indicates we should follow the process done in the proof of Theorem \ref{t:sup} and choose 
$$
k = \| \rho_{*}(\vec{F})\|_{L^{q}}^{\frac{1}{\gamma-1}} + \| f \|_{L^{\frac{q}{\gamma^{\prime}}}}^{\frac{1}{\gamma-1}}.
$$ 
Recalling $0 < \theta < \tau < 1$, setting $w = v^{- \beta/\gamma} \bar{u}$ and proceeding as in the proof of Theorem \ref{t:sup} leads to
\begin{equation} \label{e:lb6}
\left( \int_{B_{R \theta}} \bar{u}^{(\gamma-\beta)\chi} \right)^{\frac{1}{\chi}} \le C_{n, \gamma,\nu, \Lambda, p,q, R \theta, R\tau} \int_{B_{R \tau}} \bar{u}^{\gamma-\beta}.
\end{equation}
The dependence on $R \theta, R\tau$ comes from the magnitude of the derivative of the cut-off function. Moreover, in \eqref{e:lb6}, the dependence on $p$ that is not present in \eqref{e:dt13} is due to having assumed $\beta \ge \beta_{p}$. To improve readability, we write $C = C_{n,\gamma,\nu,\Lambda,p,q,R \theta,R\tau}$ and suppose $R=1$ through the end of \eqref{e:lb8.4} . 
Now, whenever $\gamma- \beta < 0$ \eqref{e:lb6} ensures
\begin{equation} \label{e:lb7}
C  \| \bar{u}\|_{L^{\gamma-\beta}(B_{\tau})} \le \| \bar{u}\|_{L^{(\gamma-\beta)\chi}(B_{\theta})}.
\end{equation} 
Since $\inf v^{+} = \lim_{p \to - \infty} \| v\|_{L^{p}}$, iterating as in Theorem \ref{t:sup}, for any $\beta > \gamma$ \eqref{e:lb7} implies,
\begin{equation} \label{e:lb8}
C\| \bar{u}\|_{L^{\gamma-\beta}(B_{\tau})} \le \inf_{B_{\theta}} \bar{u}.
\end{equation}

On the other hand, whenever $0 < \beta < \gamma-1$, \eqref{e:lb6} guarantees
\begin{equation} \label{e:lb8.1}
\| \bar{u}\|_{L^{(\gamma-\beta)\chi}(B_{\theta})} \le C  \|\bar{u}\|_{L^{\gamma-\beta}(B_{\tau})}.
\end{equation}
Recalling $- \beta < -1$, we can iterate \eqref{e:lb8.1} finitely many times depending on $p,p_{0}$ so that when $0 < p_{0} < p < (\gamma-1)\chi$, choosing $\gamma-\beta = p_{0}$ and finitely many iterations of \eqref{e:lb8.1} ensures
\begin{equation} \label{e:lb8.2}
\| \bar{u}\|_{L^{p}(B_{\theta})} \le C \cdot C(p_{0}) \| \bar{u}\|_{L^{p_{0}}(B_{\tau})}.
\end{equation}


We claim the proof is complete once we show that there exists $p_{0} > 0$ so that
\begin{equation} \label{e:lb8.3}
\int_{B_{\tau}} \bar{u}^{-p_{0}} \int_{B_{\tau}} \bar{u}^{p_{0}} \le C.
\end{equation}

Indeed, choosing $\beta = \gamma + p_{0}$ in \eqref{e:lb8} combined with \eqref{e:lb8.2} and \eqref{e:lb8.3} yields
\begin{align} \label{e:lb8.4}
\inf_{B_{1}} \bar{u} & \ge C \cdot C(p_{0}) \| \bar{u}\|_{L^{-p_{0}}} \\
\nonumber & =C \cdot C(p_{0}) \left( \| \bar{u}\|_{L^{-p_{0}}} \| \bar{u}\|_{L^{p_{0}}}^{-1} \right) \| \bar{u}\|_{L^{p_{0}}} \\
\nonumber & \ge C \cdot C(p_{0}) \| \bar{u}\|_{L^{p}}.
\end{align}
It happens that $p_{0} = p_{0}(n,\gamma,\nu,\Lambda)$ so that the constant $C(p_{0})$ can be absorbed into our universal constant. Recalling that $\bar{u} = u^{+} + k$ and using scaling, this says
$$
\inf_{B_{R \theta}} u^{+} + R^{\delta} \|\rho_{*}(\vec{F})\|_{L^{q}(B_{R})}^{\frac{1}{\gamma-1}} + R^{\gamma^{\prime} \delta} \| f\|_{L^{\frac{q}{\gamma^{\prime}}}(B_{R})}^{\frac{1}{\gamma-1}} \ge C_{n, \gamma,  \nu, \Lambda, q, p, \theta, \tau} \| u^{+}\|_{L^{p}(B_{R \tau})},
$$
as desired.


Hence, it only remains to show \eqref{e:lb8.3}. We consider the test function $\varphi = \frac{\eta^{\gamma}}{\bar{u}^{\gamma-1}}$. Note,
$$
D \varphi = \gamma \left( \frac{\eta}{\bar{u}} \right)^{\gamma-1} D \eta - (\gamma-1) \left( \frac{\eta}{\bar{u}} \right)^{\gamma} D \bar{u}.
$$
We estimate as usual:

\begin{equation} \label{e:lb9}
\left \langle \rho(D \bar{u})^{\gamma-1} (D \rho) (D \bar{u}), - (\gamma-1) \left( \frac{\eta}{\bar{u}} \right)^{\gamma} D \bar{u} \right \rangle= - (\gamma-1) \rho\left( \frac{ \eta D \bar{u}}{\bar{u}}\right)^{\gamma} .
\end{equation}
and choosing $\epsilon > 0$ so that $(\gamma-1) \epsilon^{\gamma^{\prime}} = (\gamma-1)/2$.
\begin{align} \label{e:lb10}
\nonumber \left \langle \rho(D \bar{u})^{\gamma-1} (D \rho)(D \bar{u}), \gamma \left( \frac{\eta}{\bar{u}} \right)^{\gamma-1} D \eta \right \rangle & \le   \gamma \left[\frac{ \epsilon^{\gamma^{\prime}}}{\gamma^{\prime}} \rho \left( \frac{ \eta D \bar{u}}{\bar{u}} \right)^{\gamma} + \frac{ \rho(D \eta)^{\gamma}}{\epsilon^{\gamma} \gamma} \right] \\
& \le \frac{\gamma-1}{2} \rho \left( \frac{ \eta D \bar{u}}{\bar{u}} \right)^{\gamma} + c_{\gamma} \rho(D \eta)^{\gamma}.
\end{align}


On the other hand, since $\bar{u} \ge k$
\begin{align}
\nonumber \left \langle \vec{F}, \gamma \left( \frac{\eta}{\bar{u}} \right)^{\gamma-1} D \eta \right \rangle & \ge - \gamma \rho_{*}(\vec{F})^{\gamma^{\prime}} \left( \frac{ \eta}{\bar{u}} \right)^{\gamma-1} \rho(D \eta) \\
\label{e:lb11}& \ge - \gamma  \left[ \frac{\rho_{*}(\vec{F})^{\gamma^{\prime}} \eta^{\gamma}}{k^{-\gamma} \gamma^{\prime}} + \frac{ \rho(D \eta)^{\gamma}}{\gamma} \right]
\end{align}
and choosing $\epsilon > 0$ so that $(\gamma^{\prime})^{-1} \epsilon^{\gamma} = 1/4$ guarantees
\begin{align}
\nonumber \langle \vec{F}, -(\gamma-1) \left( \frac{\eta}{\bar{u}} \right)^{\gamma} D \bar{u} \rangle & \ge - (\gamma-1) \left( \frac{\eta}{\bar{u}} \right)^{\gamma} \left[ \frac{ \epsilon^{\gamma} \rho(D \bar{u})^{\gamma}}{\gamma} + \frac{\rho_{*}(\vec{F})^{\gamma^{\prime}}}{\epsilon^{\gamma^{\prime}} \gamma^{\prime}} \right]  \\
\label{e:lb12}& \ge - \frac{(\gamma-1)}{4} \rho \left( \frac{ \eta D \bar{u}}{\bar{u}} \right)^{\gamma} - c_{\gamma} \eta^{\gamma} k^{-\gamma} \rho_{*}(\vec{F})^{\gamma^{\prime}}.
\end{align}
Finally, note
\begin{equation} \label{e:lb13}
\int f \eta^{\gamma} \bar{u}^{-(\gamma-1)} \ge - k^{-(\gamma-1)} \| f \|_{\frac{n}{\gamma}} \| \eta \|_{L^{\frac{n \gamma}{n-\gamma}}}^{\gamma} \ge - c_{\gamma,n,\nu} k^{-(\gamma-1)} \| f\|_{\frac{n}{\gamma}} \|D \eta \|_{L^{\gamma}}^{\gamma}
\end{equation}

Combining \eqref{e:lb9} - \eqref{e:lb13}, with the fact that $u$ is a supersolution yields
\begin{align}
\nonumber - \int_{\Omega} &  \eta^{\gamma} \rho \left(  D (\log \bar{u}) \right)^{\gamma} \ge - c_{n,\gamma, \nu} \left[ \int \rho( D \eta)^{\gamma} \left( 1 + \frac{\|f\|_{\frac{n}{\gamma}}}{k^{\gamma-1}} \right) + \int \frac{ \rho_{*}(\vec{F})^{\gamma^{\prime}}}{k^{\gamma}} \eta^{\gamma} \right] \\
& \ge  - c_{n,\gamma, \nu} \left[ \int \rho(D \eta)^{\gamma} \left( 1 + k^{-(\gamma-1)} \|f\|_{\frac{n}{\gamma}} + k^{-\gamma} \|\rho_{*}(\vec{F})\|_{L^{\frac{n}{\gamma-1}}}^{\gamma^{\prime}} \right) \right] ,
\end{align}
where the 2nd inequality follows from the first by Holder and Sobolev embedding similar to \eqref{e:lb13}. 
Choosing $k = \| f \|_{\frac{n}{\gamma}}^{\frac{1}{\gamma-1}} + \|\rho_{*}(\vec{F})\|_{L^{\gamma^{\prime}}}^{\frac{1}{\gamma-1}}$ this can be written
\begin{equation}  \label{e:lb14}
\int_{\Omega} \eta^{\gamma} \rho(  D ( \log \bar{u}))^{\gamma} \le c_{n, \gamma, \nu} \int \rho(D \eta)^{\gamma}.
\end{equation}
Now, for all $B_{2r}^{\prime} \subset B_{2R}$, we can choosing $\eta \equiv 1$ on $B_{r}^{\prime}$ and $\eta \equiv 0$ on $\Omega \setminus B_{2r}^{\prime}$ with $|D \eta| \le \frac{2}{r}$, \eqref{e:lb14} says

\begin{equation*}
\int_{B_{r}^{\prime}} |D ( \log \bar{u})|^{\gamma} \le c_{n,\gamma,\nu} \int_{B_{r}^{\prime}} \rho(D(\log \bar{u}))^{\gamma} \le c_{n,\gamma,\nu, \Lambda} r^{n - \gamma}
\end{equation*}
or taking the $\gamma$-root and multiplying by $r^{1 - \frac{n}{\gamma}}$ that is
\begin{equation} \label{e:lb14.1}
r^{1-\frac{n}{\gamma}} \left( \int_{B_{r}^{\prime}} |D(\log \bar{u})|^{\gamma} \right)^{\frac{1}{\gamma}} \le C_{n,\gamma,\nu,\Lambda}. 
\end{equation}

Noticing $(r^{-n})^{\frac{1}{\gamma^{*}}} = r^{1- \frac{n}{\gamma}}$, we consecutively apply Jensen's inequality, the Poincare inquality in a ball, and \eqref{e:lb14.1} to achieve
\begin{align} \label{e:lb15}
\frac{1}{r^{n}} \int_{B_{r}^{\prime}} |\log \bar{u} - (\log \bar{u})_{B_{r}^{\prime}}|  &\le C \left( \frac{1}{r^{n}} \int_{B_{r}^{\prime}} |\log \bar{u}  - (\log \bar{u})_{B_{r}^{\prime}}|^{\gamma^{*}} \right)^{\frac{1}{\gamma^{*}}} \\
\nonumber & \le C_{n, \gamma, \nu,\Lambda}.
\end{align}
Since this holds uniformly for all $B_{2r}^{\prime} \subset B_{2R}$, and hence \eqref{e:lb15} holds for all $B_{r}^{\prime} \subset B_{R}$, this ensures $\log \bar{u} \in BMO(B_{R})$ and consequently, by the John-Nirenberg lemma there exists $0< p_{0}$ depending only on $n$ and the constant in \eqref{e:lb15} so that
\begin{equation} \label{e:lb16}
\sup_{B\subset B_{3}} \frac{1}{|B|} \int_{B} e^{p_{0} |\log \bar{u} - (\log \bar{u})_{B}|} < \infty.
\end{equation}
Finally, since $e^{|a-b|} \ge 1$ we note that $e^{p_{1}|a-b|} \le e^{p_{0}|a-b|}$ if $p_{1} \le p_{0}$. In particular, without loss of generality, suppose $p_{0} < \frac{(\gamma-1)n}{n-\gamma}$. But then, \eqref{e:lb16} implies \eqref{e:lb8.3}.
\end{proof}


\begin{theorem} [Improvement of Oscillation] \label{t:ioo}
Suppose $u, \rho$, and $\gamma$ are as in Theorem \ref{t:iwh} for some $\vec{F} \in L^{q}_{\operatorname{loc}}(\Omega), f \in L^{\frac{q}{\gamma^{\prime}}}_{\operatorname{loc}}(\Omega)$ some $q > \frac{n}{\gamma-1}$. If $\overline{B_{3R}} \subset \Omega$, then for all $0 < \theta < 1$,
$$
\operatorname{osc}_{B_{\theta R}} u \le C_{n,\rho,\gamma,q}  \theta^{\alpha} \left[ \operatorname{osc}_{B_{R}} u  +  \| \rho_{*}(\vec{F})\|^{\frac{1}{\gamma-1}}_{L^{q}(B_{2R})} + \|f\|_{L^{\frac{q}{\gamma^{\prime}}}(B_{2R})}^{\frac{1}{\gamma-1}} \right].
$$
\end{theorem}

\begin{proof}
For $0 < s < 2R$ let 
$$
M_{s} = \sup_{B_{s}} u \qquad \text{and} \qquad m_{s} = \inf_{B_{s}} u.
$$
Then $v \in \{M_{2R} - u, u - m_{2r}\}$ is a positive solution of \eqref{e:dt1} on $B_{2R}$.

Consider $p = 1$ and $\theta = 1/2 < \tau < 1$ in Theorem \ref{t:inf}. That is, for a positive solution $v$,
\begin{equation} \label{e:ioo1}
\inf_{B_{R/2}} v + G(R) \ge C R^{-n} \int v
\end{equation}
where 
$$
G(R) = R^{\delta} \| \rho_{*}(\vec{F})\|^{\frac{1}{\gamma-1}}_{L^{q}(B_{R})}   + R^{\gamma^{\prime} \delta} \| f\|_{L^{\frac{q}{\gamma^{\prime}}}(B_{2R})}^{\frac{1}{\gamma-1}}.
$$
Note, $\inf_{B_{R/2}} M_{2R} - u = M_{2R} - M_{R/2}$ and $\inf_{B_{R/2}} u - m_{2R} = m_{r/2} - m_{2r}$. Therefore, applying \eqref{e:ioo1} to $M_{2R} - u$ and $u - m_{2R}$ yields
\begin{align*}
M_{2R} - M_{R/2} + G(R) &\ge  C R^{-n} \int M_{2R} - u \\
m_{R/2} - m_{2R} + G(R) & \ge C R^{-n} \int u - m_{2R}.
\end{align*}
Adding yields
$$
(M_{2R} - m_{2R}) - (M_{R/2} - m_{R/2}) + 2 G(R)\ge C (M_{2R} - m_{2R} ),
$$
or equivalently,
$$
\operatorname{osc}_{B_{R/2}} u \le (1 - C) \operatorname{osc}_{B_{2R}} u + 2\left[ R^{\delta} \| \rho_{*}(\vec{F})\|^{\frac{1}{\gamma-1}}_{L^{q}(B_{R})}   + R^{\gamma^{\prime}\delta} \| f\|_{L^{\frac{q}{\gamma^{\prime}}}(B_{2R})}^{\frac{1}{\gamma-1}} \right].
$$
Lemma \ref{l:ibs} verifies the result by choosing the parameters from Lemma \ref{l:ibs} by $\tau = 1/4$,$\tilde \delta = (1-C)$, and $\mu (1 - \frac{n}{q(\gamma-1)}) > \alpha$. The latter can be done by making $\alpha$ smaller if necessary. Notably $\alpha = \alpha(1/4,\tilde \delta)$ so it has the expected dependencies.
\end{proof}



Holder regularity is a classic result of the improvement of oscillation in Theorem \ref{t:ioo}.
\begin{corollary} 
Suppose $\rho, \gamma,$ and $u$ are as in Theorem \ref{t:iwh}, for some for $\vec{F} \in L^{q}_{\operatorname{loc}}(\Omega), f \in L^{\frac{q}{\gamma^{\prime}}}_{\operatorname{loc}}(\Omega)$ some $q > \frac{n}{\gamma-1}$. Then $u \in C_{\operatorname{loc}}^{\alpha}(\Omega)$ for some $\alpha = \alpha(n, \gamma, \nu, \Lambda, q)$.
\end{corollary}

Last, we conclude with a Liouville or Bernstein type theorem in the case that $f, \vec{F} \equiv 0$.

\begin{theorem} [Bernstein Theorem]
Let $\rho, \gamma,$ and $u$ be as in Theorem \ref{t:iwh}. Suppose additionally that $f, \vec{F} \equiv 0$. If $\Omega = \mathbb{R}^n$ then $u$ is bounded if and only if it is constant.
\end{theorem}

\begin{proof}
Since $f, \vec{F} \equiv 0$, iterating Theorem \ref{t:ioo} says there exists $0 < \gamma < 1$ so that
$$
\operatorname{osc}_{B_{R}} u \le \gamma^{k} \operatorname{osc}_{B_{2^{k} R}} \le  \gamma^{k} \left( 2 \|u \|_{L^{\infty}} \right).
$$
Taking $k$ and $r$ to infinity consecutively completes the proof.
\end{proof}

\bibliographystyle{alpha}
\bibdata{references}
