Notice that a weak solution of a inhomogenous $p$-Laplace equation is by definition $W^{1,p}_{loc}$. So, by Sobolev’s embedding theorem we have $u \in C^{0,1-n/p}$. It is know that is possible to prove that solutions of $$ \Delta_p u = f,\quad f \in L^q,\; n/p <q<n $$ are $C^{0,\frac{p-n/q}{p-1}}$. So, we have a gain of regularity only when $$ \frac{p-n/q}{p-1} > 1-n/p. \label{1}\tag{1} $$ That is, only if $q> \frac{p}{np-n+p}$. In particular, we always have gains of regularity when $q=p/(p-1)$ the conjugate exponet of $p$. In fact \eqref{1} become equivalent to $$ 1>(1-p)n. $$ Since we are assuming $p>1$ the inequality above is satisfied.
I'd like to know how is the respective problem with Orlicz Spaces. More specifically, the $p$-Laplace case is the case when the parameters $g_0$ and $g_1$ coincide with $p-1$. In the $p$-Laplacian case, we saw that we have gain of regularity when $q$ is large. What can we do in the Orlicz setting? What would be the limits of $q$ in this case (since we have two parameters involved)? That is,
I think that if someone can answer this then they can also answer this question.