1. Introduction
The aim of this paper is to present some regularity results of vectorial minimizers of variational integrals of the form
\begin{align}
\mathcal{F}(u,\Omega):=\int_\Omega \left[ f(x,Du)- g \cdot u \right] dx .
\end{align} Here, 
$\Omega \subset \mathbb{R}^n$ is a bounded open subset, 
$n\geq 2$, 
$u : \Omega \to \mathbb{R}^N$, 
$N \ge 1$, 
$f:\Omega\times\mathbb{R}^{N\times n}\to [0,+\infty)$, 
$\xi \mapsto f(x,\xi) \in \mathcal{C}^1(\mathbb{R}^{N\times n})$ and satisfies the so-called Uhlenbeck structure, i.e. f is represented in the form 
$f (x, \xi) = \tilde{f}(x, |\xi|)$ for a given function 
$\tilde{f} :\Omega\times[0,+\infty)\to [0,+\infty)$. Besides, we assume that there exist positive constants 
$\nu, L, M, l$ and exponents 
$2 \leq p \lt q \lt \infty$ such that
\begin{align}
\nu|\xi|^p \leq f(x,\xi) \leq L(1+|\xi|^q)
\end{align}
\begin{align}
\langle D_\xi f(x, \xi)-D_\xi f(x, \eta),\xi-\eta\rangle \ge M (1 +|\xi|^2+|\eta|^2)^\frac{p-2}{2}|\xi - \eta|^2
\end{align}
\begin{align}
|D_\xi f(x, \xi)-D_\xi f(x, \eta)| \le l (1 +|\xi|^2+|\eta|^2)^\frac{q-2}{2}|\xi - \eta|
\end{align} for a.e. 
$x \in \Omega$ and all 
$\xi,\eta \in \mathbb{R}^{N \times n}$. Concerning the dependence on the x-variable, we assume that there exists a non-negative function 
$k: \Omega \to [0, +\infty)$ such that
\begin{align}
|D_\xi f(x,\xi)-D_\xi f(y, \xi)| \leq |x-y|^{\gamma} (k(x)+k(y)) (1 +|\xi|^2)^\frac{q-1}{2}
\end{align} for a.e. 
$x,y \in \Omega$ and every 
$\xi \in \mathbb{R}^{N \times n}$, where 
$\gamma \in (0,1)$.
It is worth noting that assumption (Equation F4) with 
$k \in L^r_{loc}(\Omega)$ states that the partial map 
$x \mapsto D_\xi f(x,\xi)$ belongs to the Besov space 
$B_{r,\infty}^{\gamma}(\Omega, \mathbb{R}^{N \times n})$ (for precise definition and properties of Besov spaces see 
$\S$ 2.2).
On the other hand, we say that assumption (Equation F5) is satisfied if there exists a sequence of measurable non-negative functions 
$g_k \in L^{r}_{loc}(\Omega)$ such that for every 
$\Omega' \Subset \Omega$
\begin{align*}\displaystyle\sum_{k=1}^{\infty} \Vert g_k \Vert^{s}_{L^{r}(\Omega')} \lt \infty,\end{align*} for some exponent 
$s \ge 1$, and at the same time
\begin{align}
|D_{\xi}f(x,\xi)-D_{\xi}f(y, \xi)| \leq |x-y|^{\gamma} (g_k(x)+g_k(y)) (1 +|\xi|^2)^{\frac{q-1}{2}}
\end{align} for a.e. 
$x,y \in \Omega$ such that 
$2^{-k} \text{diam}(\Omega) \leq |x-y| \lt 2^{-k+1}\text{diam}(\Omega)$ and for every 
$\xi \in \mathbb{R}^{N \times n}$. As before, we remark that assumption (Equation F5) implies that the partial map 
$x \mapsto D_\xi f(x,\xi)$ belongs to the Besov space 
$B^\gamma_{r,s}(\Omega, \mathbb{R}^{N \times n})$.
Let us recall the definition of local minimizer.
Definition 1.1. A function 
$u\in W_{\rm loc}^{1,1}(\Omega, \mathbb{R}^N)$ is a local minimizer of (Equation 1.1) if, for every open subset 
$\tilde{\Omega} \Subset \Omega$, we have 
$\mathcal{F}(u, \tilde{\Omega}) \lt \infty$ and 
$\mathcal{F}(u;\tilde{\Omega})\le \mathcal{F}(\varphi;\tilde{\Omega})$ holds for all 
$\varphi\in u+W_0^{1,1}(\tilde{\Omega},\mathbb{R}^N)$.
We will say that a function F has (p, q)-growth conditions if assumption (Equation F1) is in force. The study of regularity properties of local minima to functionals with general growth started with the pioneering papers by Marcellini [Reference Marcellini29], see also [Reference Marcellini30, Reference Marcellini31].
When referring to (p, q)-growth conditions (Equation F1), we call the quantity 
$q/p \gt 1$ the gap ratio of the integrand F, or simply, the gap. A main point is that in order to get regularity of minimizers to non-autonomous functionals with non-standard growth conditions, even the local boundedness, a restriction between p and q need to be imposed, usually expressed in the form
\begin{align*}
q \le c(n,r)p, \qquad \text{with } c(n,r ) \to 1 \ \text{as} \ r \to n,
\end{align*} where here r is the exponent appearing in (Equation F4)-(Equation F5) and usually denotes the degree of regularity of the map 
$x \mapsto D_\xi f(x, \xi)$. We refer to [Reference Giaquinta17, Reference Marcellini30] for counterexamples.
In order to explain the difficulties of our problem, we now recall some class of functionals included in our setting and some known regularity properties of the related minimizers.
It is well known that no extra differentiability properties for minimizers to
\begin{align*}
\mathcal{G}(u,\Omega):=\int_\Omega f(x,Du) \ dx
\end{align*} can be expected, unless some assumption is given on the coefficients of the operator 
$D_\xi f(x,\xi)$. A 
$W^{1,r}$ Sobolev regularity, with 
$r \ge n$, on the partial map 
$x \mapsto D_\xi f(x,\xi)$ is a sufficient condition for the higher differentiability of minima (see [Reference Eleuteri and Passarelli di Napoli10, Reference Eleuteri and Passarelli di Napoli11, Reference Gavioli14–Reference Gentile16, Reference Giova18, Reference Giova and Passarelli di Napoli19, Reference Grimaldi21, Reference Passarelli di Napoli33–Reference Russo35]). In the case of standard p-growth, the higher fractional differentiability of minimizers has been proved assuming that the coefficients of 
$ D_\xi f(x,\xi)$ belong to the Besov space 
$B^\alpha_{\frac{n}{\alpha},s}$ (see [Reference Baisón, Clop, Giova, Orbitg and Passarelli di Napoli1, Reference Baisón, Clop and Orobitg2, Reference Clop, Faraco and Ruiz5, Reference Cruz, Mateu and Orobitg7, Reference Kristensen and Mingione28]). Regarding the non-standard growth, in [Reference Grimaldi and Ipocoana22, Reference Grimaldi and Ipocoana23] this analysis has been carried out in the setting of obstacle problems, assuming that the map 
$x \mapsto D_\xi f(x,\xi)$ satisfies a 
$B^\alpha_{r,s}$-regularity for 
$r \gt \frac{n}{\alpha}$ (see [Reference Grimaldi and Ipocoana24] for the unconstrained case). It is worth noticing that solutions u to obstacle problems solve elliptic equations of the form
\begin{align}
\text{div} D_\xi f(x,Du) = \text{div} D_\xi f(x,D \psi)=: \bar{g},
\end{align} where ψ is the obstacle function and 
$u \ge \psi$. In [Reference Grimaldi and Ipocoana22, Reference Grimaldi and Ipocoana23], the higher differentiability has been obtained assuming that 
$D \psi \in W^{1,2q-p}(\Omega)$, that implies 
$\bar{g} \in L^{2q-p}(\Omega)$. Here, in Theorem 1.2 below, we will prove some higher differentiability properties of minimizers under weaker assumptions on the right-hand side.
Actually, the research in this field is so intense that it is almost impossible to give an exhaustive and comprehensive list of references; in addition to those mentioned in this paper and the references therein, we confine ourselves to refer the interested reader to the survey [Reference Mingione and Radulescu32].
In particular in [Reference Clop, Giova and Passarelli di Napoli6], the authors proved the extra fractional differentiability of weak solutions of the following nonlinear elliptic equations in divergence form
\begin{align*}\text{div} \mathcal{A}(x,Du)= \text{div } G \quad \text{in } \Omega,\end{align*} where the operator 
$\mathcal{A}$ satisfies (Equation F2)–(Equation F4) with p = q. It has been proved that a higher differentiability on the partial map 
$x \mapsto \mathcal{A}(x, \cdot) $ and on the right-hand side G transfer to the function 
$V_p(Du)$ (see 
$\S$ 2 for the definition).
The main novelty of our results is that we consider (p, q)-growth, which means we are able to extend the results of [Reference Clop, Giova and Passarelli di Napoli6] to the context of non-standard growth. Moreover, we show that the gradient of minimizers to (Equation 1.1) posseses some fractional differentiability properties, without assuming any differentiability properties on the function g, but only a suitable order of integrability in the setting of Lebesgue spaces.
Our first result is the following
Theorem 1.2 Let f satisfy (Equation F1)–(Equation F4) and let 
$g \in L^{m}_{loc}(\Omega, \mathbb{R}^N)$, with 
$p' \leq m \leq 2$, for exponents 
$2\leq p \lt \frac{n}{\lambda} \lt r$, 
$\textstyle p \lt q $ such that
\begin{align}
\dfrac{q}{p} \lt 1+ \dfrac{\lambda}{n}- \frac{1}{r} ,
\end{align} where 
$\lambda:= \min \{\gamma, \frac{m}{2} \}$. Let 
$u\in W^{1,p}_{loc}(\Omega, \mathbb{R}^N)$ be a local minimizer of (Equation 1.1). Then, 
$V_p(Du) \in B^\lambda_{2,\infty,\textrm{\text{loc}}}(\Omega, \mathbb{R}^{N \times n})$ and the following estimate holds
\begin{align}
\int_{B_{R/4}} |\tau_h V_p (Du)|^2 dx \leq c |h|^{2 \lambda}& \left( 1+\int_{B_{R}} (1+f(x,Du)) dx + \Vert u \Vert_{W^{1,p}(B_R)}\right. \nonumber\\ & \qquad \left. + \Vert k \Vert_{L^r (B_{R})} + \Vert g \Vert_{L^{m}(B_R)} \vphantom{\int_{B_{R}}} \right)^\sigma,
\end{align} for every concentric balls 
$B_{R/4} \subset B_{R} \Subset \Omega$, where 
$c = c(n,p,q,m,\gamma,r, M,l,R)$ and 
$\sigma= \sigma (n,p,q,m,\gamma,r)$ are positive constants.
One of the main motivations to the study of functionals with (p, q)-growth comes from the applications, for instance to the theory of elasticity for strongly anisotropic materials (see Zhikov [Reference Zhikov37, Reference Zhikov38]). A very well-known model is given by the so-called double phase functional defined by
\begin{align}
\int_\Omega \left[ |Du|^p+a(x)|Du|^q \right] dx,
\end{align}where a(x) is non-negative and Hölder continuous with exponent α. The regularity properties of local minimizers to such functional have been widely investigated in [Reference Baroni, Colombo and Mingione3, Reference Esposito, Leonetti and Mingione12, Reference Fonseca, Maly and Mingione13] with the aim of identifying sufficient and necessary conditions on the relation between p, q and α to establish the Hölder continuity of the gradient of the local minimizers.
Very recently, in [Reference De Filippis and Mingione8, Reference De Filippis and Mingione9] it has been proved the local gradient Hölder continuity of minimizers of non-uniformly elliptic integrals, that are not necessarily equipped with a Euler-Lagrange equation due to the mere local Hölder continuity of coefficients, under the gap
\begin{align}
\dfrac{q}{p} \lt 1 + \dfrac{\alpha}{n}.
\end{align} We observe that in the model case (Equation 1.5), the inequality (Equation 1.3) gives back (Equation 1.6). Indeed, by Lemma 2.14 below, if 
$a(x) \in B^\lambda_{r,\infty}(\Omega)$ then 
$a(x) \in \mathcal{C}^{0,\alpha}(\Omega)$ with exponent
\begin{align*}
\alpha=\lambda-\dfrac{n}{r}.
\end{align*} Moreover, in [Reference De Filippis and Mingione8] the authors proved the Hölder continuity of Du assuming that 
$g \in L^{n/\alpha,1/2}(\Omega)$. However, in Theorem 1.2 we impose a weaker summability on the function g, that is 
$g \in L^{m}_{loc}(\Omega, \mathbb{R}^N)$ for 
$p' \leq m \leq 2$, and so our results are not covered by those in [Reference De Filippis and Mingione8].
Now assuming (Equation F5) instead of (Equation F4) and arguing in a similary way, we obtain the following theorem.
Theorem 1.3 Let f satisfy (Equation F1)–(Equation F3), (Equation F5) and let 
$g \in L^{m}_{loc}(\Omega, \mathbb{R}^N)$, with 
$p' \leq m \leq 2$, for exponents 
$2\leq p \lt \frac{n}{\lambda} \lt r$, 
$\textstyle p \lt q $ such that (Equation 1.3) holds. Let 
$u\in W^{1,p}_{loc}(\Omega, \mathbb{R}^N)$ be a local minimizer of (Equation 1.1). Then, 
$V_p(Du) \in B^\lambda_{2,s,\textrm{\text{loc}}}(\Omega, \mathbb{R}^{N \times n})$ and the following estimate holds
\begin{align*}
\displaystyle\int_{B_{1}} \biggl( \displaystyle\int_{B_{R/4}} \dfrac{|\tau_hV_p (Du)|^{2}}{|h|^{2 \lambda }} dx \biggr)^{\frac{s}{2}} \dfrac{dh}{|h|^{n}} \leq c\Biggl( & 1+ \int_{B_{R}} (1 + f(x,Du)) dx+ \Vert u \Vert_{W^{1,p}(B_R)} \notag\\
& \qquad+ \sum_{k=1}^{\infty }\Vert g_k \Vert_{L^r (B_{R})} + \Vert g \Vert_{L^{m}(B_R)} \Biggr)^\sigma, \notag
\end{align*} for every concentric balls 
$B_R \subset B_{2R} \Subset \Omega$, where 
$c = c(n,p,q,m,\gamma,r, M,l,R, s)$ and 
$\sigma= \sigma (n,p,q,m,\gamma,r, s)$ and 
$\lambda:= \min \{\gamma, \frac{m}{2} \}$ are positive constants.
We briefly describe the structure of the paper. After summarizing some known results and fixing few notation, we recall an approximation lemma in 
$\S$ 3. Section 4 is devoted to the proof of some a priori estimates for solutions to a family of approximating problems. Next, in 
$\S$ 5, we pass to the limit in the approximating problems.
2. Preliminaries
In this section, we introduce some notations and collect several results that we shall use to establish our main result.
We will denote by c or C a general constant that may vary on different occasions, even within the same line of estimates. Relevant dependencies on parameters and special constants will be suitably emphasized using parentheses or subscripts.
In what follows, 
$B(x,r)=B_{r}(x)= \{y \in \mathbb{R}^{n} : |y-x | \lt r \}$ will denote the ball centred at x of radius r. We shall omit the dependence on the centre and on the radius when no confusion arises.
We define an auxiliary function by
\begin{align*}
V_{p}(\xi):=( 1 +|\xi|^{2})^\frac{p-2}{4} \xi
\end{align*} for all 
$\xi\in \mathbb{R}^{m}$, 
$m \in \mathbb{N}$. For the function Vp, we recall the following estimate (see e.g. [Reference Giusti20, Lemma 8.3]).
Lemma 2.1. Let 
$1 \lt p \lt +\infty$. There exists a constant 
$c=c(n,p) \gt 0$ such that
\begin{align*}
c^{-1}(1+|\xi|^{2}+|\eta|^{2})^{\frac{p-2}{2}} \leq \dfrac{|V_{p}(\xi)-V_{p}(\eta)|^{2}}{|\xi-\eta|^{2}} \leq c(1+|\xi|^{2}+|\eta|^{2})^{\frac{p-2}{2}}
\end{align*} for any 
$\xi, \eta \in \mathbb{R}^{n}$, 
$\xi \neq \eta$.
Now we state a well-known iteration lemma (see [Reference Giusti20, Lemma 6.1] for the proof).
Lemma 2.2. Let 
$\Phi : [\frac{R}{2},R] \rightarrow \mathbb{R}$ be a bounded nonnegative function, where R > 0. Assume that for all 
$\frac{R}{2} \leq r \lt s \leq R$ it holds
\begin{align*}\Phi (r) \leq \theta \Phi(s) +A + \dfrac{B}{(s-r)^2}+ \dfrac{C}{(s-r)^{\gamma}}\end{align*} where 
$\theta \in (0,1)$, A, B, 
$C \geq 0$ and γ > 0 are constants. Then there exists a constant 
$c=c(\theta, \gamma)$ such that
\begin{align*}\Phi \biggl(\dfrac{R}{2} \biggr) \leq c \biggl( A+ \dfrac{B}{R^2}+ \dfrac{C}{R^{\gamma}} \biggr).\end{align*}2.1. Difference quotient
We recall some properties of the finite difference quotient operator that will be needed in the sequel.
Definition 2.3. Let F be a function defined in an open set 
$\Omega \subset \mathbb{R}^n$ and let 
$h \in \mathbb{R}^n$. We call the difference quotient of F with respect to h the function
\begin{align*}
\tau_{h}F(x) :=F(x+h)-F(x) .\end{align*} The function 
$\tau_{h}F$ is defined in the set
\begin{align*}\tau_{h}\Omega := \{ x \in \Omega : x+h \in \Omega \},\end{align*}and hence in the set
\begin{align*}\Omega_{|h|}: = \{x \in \Omega : \mathrm{dist}(x,\partial \Omega) \gt |h| \}.\end{align*}We start with the description of some elementary properties that can be found, for example, in [Reference Giusti20].
Proposition 2.4. Let 
$F \in W^{1,p}(\Omega)$, with 
$p \geq1$, and let 
$G:\Omega \rightarrow \mathbb{R}$ be a measurable function. Then
(i) 
$\tau_{h}F \in W^{1,p}(\Omega_{|h|})$ and
\begin{align*}D_{i}(\tau_{h}F)=\tau_{h}(D_{i}F).\end{align*} (ii) If at least one of the functions F or G has support contained in 
$\Omega_{|h|}$, then
\begin{align*}\displaystyle\int_{\Omega}F \tau_h G dx = \displaystyle\int_{\Omega} G \tau_{-h}F dx.\end{align*}(iii) We have
\begin{align*}\tau_{h} (FG)(x)= F(x+h)\tau_{h} G(x)+G(x) \tau_{h} F(x).\end{align*}The next result about the finite difference operator is a kind of integral version of Lagrange Theorem.
Lemma 2.5. If 
$0 \lt \rho \lt R,$ 
$|h| \lt \frac{R-\rho}{2},$ 
$1 \lt p \lt +\infty$ and 
$F\in W^{1,p}(B_{R})$, then
\begin{align*}
\displaystyle\int_{B_{\rho}} |\tau_{h}F(x)|^{p} dx \leq c(n,p)|h|^{p} \displaystyle\int_{B_{R}} |DF(x)|^{p} dx.
\end{align*}Moreover,
\begin{align*}
\displaystyle\int_{B_{\rho}} |F(x+h)|^{p} d x \leq \displaystyle\int_{B_{R}} |F(x)|^{p}d x.
\end{align*}We recall a fractional version of Sobolev embedding property (see [Reference Esposito, Leonetti and Mingione12, Lemma 3]).
Lemma 2.6. Let 
$F \in L^2(B_R)$. Suppose that there exist 
$\rho \in (0,R)$, 
$0 \lt \alpha \lt 1$ and M > 0 such that
\begin{align*} \displaystyle\int_{B_{\rho}}|\tau_{h}F(x)|^2 dx \leq M^2 |h|^{2 \alpha}, \end{align*} for every h such that 
$|h| \lt \frac{R-\rho}{2}$. Then 
$F \in L^{\frac{2n}{n-2 \beta}}(B_{\rho})$ for every 
$\beta \in (0, \alpha)$ and
\begin{align*}\Vert F \Vert _{L^{\frac{2n}{n-2 \beta}}(B_{\rho})} \leq \frac{c}{(R-\rho)^{2\beta+2\alpha+2}} (M+ \Vert F \Vert _{L^2(B_R)}), \end{align*} with 
$c=c(n, \alpha, \beta)$.
2.2. Besov spaces
We give the definition of Besov spaces as done in [Reference Haroske26, Section 2.5.12].
Definition 2.7. Let 
$1 \le p \lt +\infty$ and 
$0 \lt \alpha \lt 1 $. Given 
$1 \leq s \lt +\infty$, we say that a function 
$v : \mathbb{R}^n \rightarrow \mathbb{R}^m$, 
$m \in \mathbb{N}$, belongs to the Besov space 
$B^{\alpha}_{p,s}(\mathbb{R}^{n})$ if, and only if, 
$v \in L^{p}(\mathbb{R}^{n})$ and
\begin{align*}
[v]_{B^{\alpha}_{p,s}(\mathbb{R}^{n})}: = \biggl( \displaystyle\int_{\mathbb{R}^{n}} \biggl( \displaystyle\int_{\mathbb{R}^{n}} \dfrac{|\tau_hv(x)|^{p}}{|h|^{\alpha p}} dx \biggr)^{\frac{s}{p}} \dfrac{dh}{|h|^{n}} \biggr)^{\frac{1}{s}} \lt + \infty .
\end{align*} Equivalently, we could simply say that 
$v \in L^{p}(\mathbb{R}^{n})$ and 
$\frac{\tau_{h}{v}}{|h|^{\alpha}} \in L^{s}\bigl( \frac{dh}{|h|^{n}}; L^{p}(\mathbb{R}^{n}) \bigr)$.
When 
$s = + \infty$, the Besov space 
$B^{\alpha}_{p,\infty}(\mathbb{R}^{n})$ consists of functions 
$v \in L^{p}(\mathbb{R}^{n})$ such that
\begin{align*}
[v]_{B^{\alpha}_{p, \infty}(\mathbb{R}^{n})} := \displaystyle\sup_{h \in \mathbb{R}^{n}} \biggl( \displaystyle\int_{\mathbb{R}^{n}} \dfrac{|\tau_hv(x)|^{p}}{|h|^{\alpha p}} dx \biggr)^{\frac{1}{p}} \lt +\infty .
\end{align*} Accordingly, for 
$1 \le s \le + \infty$, the Besov space 
$B^{\alpha}_{p,s}(\mathbb{R}^{n})$ is normed with
\begin{align*}
\Vert v \Vert_{B^{\alpha}_{p,s}(\mathbb{R}^{n})}: = \Vert v \Vert_{L^{p}(\mathbb{R}^{n})} + [v]_{B^{\alpha}_{p,s}(\mathbb{R}^{n})} .
\end{align*} Observe that, if 
$1 \le s \lt + \infty$, integrating for 
$h \in B(0, \delta)$ for a fixed δ > 0 then an equivalent norm is obtained, because
\begin{align*}
\biggl( \displaystyle\int_{\{|h| \geq \delta\}} \biggl( \displaystyle\int_{\mathbb{R}^{n}} \dfrac{|\tau_h v(x)|^{p}}{|h|^{\alpha p}} dx \biggr)^{\frac{s}{p}} \dfrac{dh}{|h|^{n}} \biggr)^{\frac{1}{s}} \leq c(n, \alpha,p,s, \delta) \Vert v \Vert_{L^{p}(\mathbb{R}^{n})} .
\end{align*} In the case 
$s = + \infty$, one can simply take supremum over 
$|h| \leq \delta$ and obtain an equivalent norm. By construction, one has 
$B^{\alpha}_{p, s}(\mathbb{R}^{n}) \subset L^{p}(\mathbb{R}^{n})$. One also has the following version of Sobolev embeddings (a proof can be found at [Reference Haroske26, Proposition 7.12]).
Lemma 2.8. Suppose that 
$0 \lt \alpha \lt 1$.
(a) If 
$1 \lt p \lt \frac{n}{\alpha}$ and 
$1 \leq s \leq p^{*}_{\alpha} = \frac{np}{n- \alpha p}$, then there is a continuous embedding 
$B^{\alpha}_{p, s}(\mathbb{R}^{n}) \subset L^{p^{*}_{\alpha}}(\mathbb{R}^{n})$.
(b) If 
$p = \frac{n}{\alpha}$ and 
$1 \leq s \leq + \infty$, then there is a continuous embedding 
$B^{\alpha}_{p, s}(\mathbb{R}^{n}) \subset BMO(\mathbb{R}^{n})$,
where BMO denotes the space of functions with bounded mean oscillations [Reference Giusti20, Chapter 2].
We recall the following inclusions between Besov spaces ([Reference Haroske26, Proposition 7.10 and Formula (7.35)]).
Lemma 2.9. Suppose that 
$0 \lt \beta \lt \alpha \lt 1 $.
(a) If 
$1 \lt p \lt +\infty$ and 
$1 \leq s \leq t \leq + \infty$, then 
$B^{\alpha}_{p, s}(\mathbb{R}^{n}) \subset B^{\alpha}_{p, t}(\mathbb{R}^{n})$.
(b) If 
$1 \lt p \lt +\infty$ and 
$1 \leq s , t \leq + \infty$, then 
$B^{\alpha}_{p, s}(\mathbb{R}^{n}) \subset B^{\beta}_{p, t}(\mathbb{R}^{n})$.
(c) If 
$1 \le s \le + \infty$, then 
$B^{\alpha}_{\frac{n}{\alpha}, s}(\mathbb{R}^{n}) \subset B^{\beta}_{\frac{n}{\beta}, s}(\mathbb{R}^{n})$.
Combining Lemmas 2.8 and 2.9, we get the following Sobolev type embedding theorem for Besov spaces 
$B^\alpha_{p,\infty}(\mathbb{R}^n)$.
Lemma 2.10. Suppose that 
$0 \lt \alpha \lt 1$ and 
$1 \lt p \lt \frac{n}{\alpha}$. There is a continuous embedding 
$B^\alpha_{p,\infty}(\mathbb{R}^n) \subset L^{p^*_\beta}(\mathbb{R}^n)$, for every 
$0 \lt \beta \lt \alpha$. Moreover, for every function 
$F \in B^\alpha_{p,\infty}(\mathbb{R}^n)$ the following local estimate
\begin{align*}
\Vert F \Vert _{L^{\frac{np}{n-\beta p}}(B_\varrho)} \leq \dfrac{c}{(R-\varrho)^{\delta}} (\Vert F \Vert_{L^{p}(B_R)} + [F]_{B^{\alpha}_{p,\infty}(B_R)}) \notag
\end{align*} holds for every ball 
$B_\varrho \subset B_R$, with 
$c=c(n,p,\alpha,\beta)$ and 
$\delta=\delta(n,p)$.
Given a domain 
$\Omega \subset \mathbb{R}^{n}$, we say that a function 
$v: \mathbb{R}^n \rightarrow \mathbb{R}^m$ belongs to the local Besov space 
$ B^{\alpha}_{p, s,\text{loc}}$ if 
$\varphi v \in B^{\alpha}_{p, s}(\mathbb{R}^{n})$ whenever 
$\varphi \in \mathcal{C}^{\infty}_{0}(\Omega)$. It is worth noticing that one can prove suitable version of Lemmas 2.8 and 2.9, by using local Besov spaces.
We have the following lemma (we refer to [Reference Baisón, Clop, Giova, Orbitg and Passarelli di Napoli1, Lemma 7] for the proof).
Lemma 2.11. Let 
$1 \le p \lt + \infty$, 
$1 \le s \le + \infty$ and 
$0 \lt \alpha \lt 1$. A function 
$v \in L^{p}_{\text{loc}}(\Omega)$ belongs to the local Besov space 
$B^{\alpha}_{p,s,\text{loc}}$, if, and only if,
\begin{align*}
\biggl\Vert \dfrac{\tau_{h}v}{|h|^{\alpha}} \biggr\Vert_{L^{s}\bigl(\frac{dh}{|h|^{n}};L^{p}(B)\bigr)} \lt + \infty,
\end{align*} for any ball 
$B\subset2B\subset\Omega$ with radius rB. Here the measure 
$\frac{dh}{|h|^n}$ is restricted to the ball 
$B(0,r_B)$ on the h-space.
It is known that Besov spaces of fractional order 
$\alpha \in (0,1)$ can be characterized in pointwise terms. We give the following definition according to [Reference Koskela, Yang and Zhou27].
Definition 2.12. Given a measurable function 
$v:\mathbb{R}^{n} \rightarrow \mathbb{R}^m$, a fractional α-Hajlasz gradient for v is a sequence 
$\{g_{k}\}_{k}$ of measurable, non-negative functions 
$g_{k}:\mathbb{R}^{n} \rightarrow \mathbb{R}$, together with a null set 
$A\subset\mathbb{R}^{n}$, such that the inequality
\begin{align*}
|v(x)-v(y)|\leq (g_{k}(x)+g_{k}(y))|x-y|^{\alpha}
\end{align*} holds whenever 
$k \in \mathbb{Z}$ and 
$x,y \in \mathbb{R}^{n}\setminus A$ are such that 
$2^{-k} \leq|x-y| \lt 2^{-k+1}$. We say that 
$\{g_{k}\}_{k} \in l^{s}(\mathbb{Z};L^{p}(\mathbb{R}^{n}))$ if
\begin{align*}
\Vert \{g_{k}\}_{k} \Vert_{l^{s}(L^{p})}=\biggl(\displaystyle\sum_{k \in \mathbb{Z}}\Vert g_{k} \Vert^{s}_{L^{p}(\mathbb{R}^{n})} \biggr)^{\frac{1}{s}} \lt +\infty.
\end{align*}The following result was proved in [Reference Koskela, Yang and Zhou27].
Theorem 2.13 Let 
$0 \lt \alpha \lt 1,$ 
$1 \leq p \lt +\infty$ and 
$1\leq s \leq +\infty $. Let 
$v \in L^{p}(\mathbb{R}^{n})$. One has 
$v \in B^{\alpha}_{p,s}(\mathbb{R}^{n})$ if, and only if, there exists a fractional α-Hajlasz gradient 
$\{g_{k}\}_{k} \in l^{s}(\mathbb{Z};L^{p}(\mathbb{R}^{n}))$ for v. Moreover,
\begin{align*}
\Vert v \Vert_{B^{\alpha}_{p,s}(\mathbb{R}^{n})}\simeq \inf \Vert \{g_{k}\}_{k} \Vert_{l^{s}(L^{p})},
\end{align*}where the infimum runs over all possible fractional α-Hajlasz gradients for v.
We conclude this section with the following embedding result (see [Reference Triebel36, Section 2.7.1, Remark 2.]).
Lemma 2.14. Let 
$0 \lt \alpha \lt 1$ and 
$ 1 \le s \le + \infty$. If 
$p \gt \frac{n}{\alpha}$, then we have the continuous embedding 
$B^\alpha_{p,s}(\mathbb{R}^n) \subset \mathcal{C}^{0,\alpha-\frac{n}{p}}(\mathbb{R}^n)$.
3. An approximation lemma
The main tool to prove Theorem 1.2 is the following approximation lemma, which allows us to approximate from below the function f with a sequence of functions 
$(f_j)$ satisfying standard p-growth conditions and the same regularity properties of f w.r.t. the x-variable (see [Reference Grimaldi and Ipocoana23]).
Lemma 3.1. Let 
$f : \Omega \times \mathbb{R}^{N \times n} \rightarrow [0,+\infty), f = f(x, \xi)$, be a Carathéodory function satisfying the Uhlenbeck structure and assumptions (Equation F1), (Equation F2), (Equation F3), (Equation F4). Then, there exists a sequence 
$(f_j)$ of Carathéodory functions 
$f_j: \Omega \times \mathbb{R}^{N \times n} \rightarrow [0,+\infty)$, convex with respect to the last variable, monotonically convergent to f, such that for a.e. 
$x,y \in \Omega$ and every 
$\xi,\eta \in \mathbb{R}^{N \times n}$
(i)
$f_j(x,\xi) = \tilde{f}_j (x,|\xi|)$;(ii)
$f_j(x,\xi)\leq f_{j+1}(x, \xi) \leq f(x,\xi)$,
$\forall j \in \mathbb{N}$;(iii)
$ \langle D_\xi f_j(x, \xi)-D_\xi f_j(x, \eta),\xi-\eta\rangle \ge M_1 (1 +|\xi|^2+|\eta|^2)^\frac{p-2}{2}|\xi - \eta|^2 $, with
$M_1=M_1(M,p) \gt 0$;(iv)
$ |D_\xi f_j(x, \xi)-D_\xi f_j(x, \eta)| \le l_1 (1 +|\xi|^2+|\eta|^2)^\frac{q-2}{2}|\xi - \eta|$, with
$l_1=l_1(l,p) \gt 0$;(v) there exist
$L_1 \gt 0$, independent of j, and
$\bar{L}_1$, depending on j, such that
\begin{align*}
& 1/L_1(|\xi|^p-1) \leq f_j(x,\xi)\leq L_1(1 + |\xi|^2)^\frac{q}{2},\\
&f_j (x,\xi)\leq \bar{L}_1(j)(1 + |\xi|^2)^\frac{p}{2};
\end{align*}(vi) there exists a constant
$C(j) \gt 0$ such that
\begin{align*}
&|D_\xi f_j(x,\xi)-D_\xi f_j(y,\xi)| \leq |x-y|^{\gamma} {(k(x)+k(y))}(1 +|\xi|^2)^\frac{q-1}{2},\\
&|D_\xi f_j(x,\xi)-D_\xi f_j(y,\xi)| \leq C(j)|x-y|^{\gamma} {(k(x)+k(y))}(1 +|\xi|^2)^\frac{p-1}{2}.
\end{align*}
Moreover, supposing (Equation F5) instead of (Equation F4), statement (vi) would change as follows:
(vi) there exists a constant
$C(j) \gt 0$ such that
\begin{align*}
&|D_\xi f_j(x,\xi)-D_\xi f_j(y,\xi)| \leq |x-y|^{\gamma} (g_k(x)+g_k(y))(1 +|\xi|^2)^\frac{q-1}{2},\\
&|D_\xi f_j(x,\xi)-D_\xi f_j(y,\xi)| \leq C(j)|x-y|^{\gamma} (g_k(x)+g_k(y))(1 +|\xi|^2)^\frac{p-1}{2}
\end{align*}for a.e.
$x,y \in \Omega$ such that
$2^{-k} \text{diam}(\Omega) \leq |x-y| \lt 2^{-k+1}\text{diam}(\Omega)$ and for every
$\xi \in \mathbb{R}^{N \times n}$.
4. A priori estimate
We consider the following approximating functionals
\begin{align}
\int_{\Omega}\left( f_j (x, Dw) - g \cdot w \right) dx ,
\end{align} where 
$(f_j)$ is the sequence of functions obtained applying Lemma 3.1 to the integrand f of the functional at (Equation 1.1). We denote by uj a local minimizer of functional (Equation 4.1). Since fj satisfies standard p-growth, then uj solves the Euler-Lagrange system
\begin{align}
\int_{\Omega} \langle D_{\xi}f_j (x, Du_{j}), D \varphi \rangle \, dx = \int_{\Omega} g \cdot \varphi \, dx,
\end{align} for every 
$ \varphi \in W^{1,p}_{0}(\Omega, \mathbb{R}^N).$
The main aim of this section is to establish the following
Theorem 4.1 Let 
$g \in L^{m}_{loc}(\Omega, \mathbb{R}^N)$, with 
$p' \leq m \leq 2$. Let 
$u_j \in W^{1,p}(\Omega, \mathbb{R}^{N })$ be a local minimizer of (Equation 4.1) under assumptions (Equation F1)–(Equation F4), for some exponents 
$2\leq p \lt \frac{n}{{\lambda}} \lt r$, where 
$\lambda:= \min \{\gamma, \frac{m}{2} \}$, 
$p \lt q$ such that (Equation 1.3) is in force. If we a priori assume that
\begin{align*}V_p(Du_j) \in B^{{\lambda}}_{2,\infty,\text{loc}}(\Omega,\mathbb{R}^{N \times n}),\end{align*}then the following estimates
\begin{align*}
\left( \int_{B_{R/4}}|Du_j|^\frac{np}{n-2 \beta} dx \right)^\frac{n-2 \beta}{n} \leq c_1 \left( 1+\int_{B_{R}} (1 + |Du_j|^p) dx + \Vert k \Vert_{L^r (B_{R})} + \Vert g \Vert_{L^{m}(B_R)} \right)^{\sigma_1}
\end{align*}for every 
$0 \lt \beta \lt \lambda$ and
\begin{align*}
\int_{B_{R/4}} |\tau_h V_p (Du_j)|^2 dx \leq c_2 {|h|^{2 \lambda} }\left( 1+\int_{B_{R}} (1 + |Du_j|^p) dx + \Vert k \Vert_{L^r (B_{R})} + \Vert g \Vert_{L^{m}(B_R)} \right)^{\sigma_2}
\end{align*} hold for every concentric balls 
$B_{R/4} \subset B_{R} \Subset \Omega$, where 
$ c_1 = c_1(n,p,q,m,\gamma,r, M,l,R, \beta)$ and 
$\sigma_1= \sigma_1 (n,p,q,m,\gamma,r, \beta)$ are positive constants independent of j. The same happens to c 2 and σ 2, but they do not depend on β.
Proof. In the following, we write u instead of uj for simplicity of notation. For further needs, we notice that from the hypothesis
\begin{align*}V_p(Du) \in B^\lambda_{2,\infty,{loc}}(\Omega,\mathbb{R}^{N \times n})\end{align*}and Lemma 2.10, we have
\begin{align}
Du \in L^{\frac{np}{n-2 \beta}}_{loc}(\Omega ,\mathbb{R}^{N \times n}),
\end{align} for all 
$0 \lt \beta \lt \lambda$. One can easily check that (Equation 1.3) yields
\begin{align}
\dfrac{r(2q-p)}{r-2} \lt \dfrac{np}{n-2 \beta} ,
\end{align} for every 
$\beta \in (\frac{\lambda n r}{nr+2(\lambda r -n)}, \lambda)$.
Fix a ball 
$B_R \Subset \Omega$ and consider radii 
$\frac{R}{4} \lt \rho \lt t' \lt R$. We divide the interval 
$(\rho, t')$ into three equal parts by means of s and t, so that 
$\frac{R}{4} \lt \rho \lt s \lt t \lt t' \lt \frac{R}{2}$, with
\begin{align*}s-\rho = t-s=t'-t=\dfrac{t'-\rho}{3}.\end{align*} Let us take 
$t^{\prime\prime}=2t'-t$ so that 
$t' \lt t^{\prime\prime} \lt R$ and 
$t^{\prime\prime}-t'=t'-t$. Let us consider a cut-off function 
$\eta \in C_0^{\infty}(B_t)$, with η = 1 on Bs, 
$0 \leq \eta \leq 1$, 
$|D \eta | \leq \frac{c}{t-s}$ and 
$|h|\leq \frac{t' - t}{2}$. We test (Equation 4.2) with the function
\begin{align*}\varphi = \tau_{-h} \left( \eta^2 \tau_h u \right) \end{align*}thus obtaining
\begin{align*}
\int_{\Omega} \langle D_{\xi}f_j (x, Du), \tau_{-h} D \left( \eta^2 \tau_h u \right) \rangle \, dx = \int_{\Omega} g \cdot \tau_{-h} \left( \eta^2 \tau_h u \right) \, dx,
\end{align*}hence by (ii) of Proposition 2.4
\begin{align}
\int_{\Omega} \langle \tau_{h} D_{\xi}f_j (x, Du), D \left( \eta^2 \tau_h u \right) \rangle \, dx = \int_{\Omega} g \cdot \tau_{-h} \left( \eta^2 \tau_h u \right) \, dx =: B.
\end{align}Now we rewrite the integral in the left-hand side of the previous equality as follows
\begin{align}
\int_{\Omega}& \langle D_{\xi}f_j (x+h, Du(x+h))- D_{\xi}f_j (x, Du(x)) , 2 \eta D \eta \tau_h u \rangle \, dx \notag \\
&+ \int_{\Omega} \langle D_{\xi}f_j (x+h, Du(x+h))- D_{\xi}f_j (x, Du(x)) , \eta^2 \tau_h Du \rangle \, dx \notag \\
&= \int_{\Omega} \langle D_{\xi}f_j (x+h, Du(x+h))- D_{\xi}f_j (x+h, Du(x)) , 2 \eta D \eta \tau_h u \rangle \, dx \notag \\
&+ \int_{\Omega} \langle D_{\xi}f_j (x+h, Du(x))- D_{\xi}f_j (x, Du(x)) , 2 \eta D \eta \tau_h u \rangle \, dx \notag \\
&+ \int_{\Omega} \langle D_{\xi}f_j (x+h, Du(x+h))- D_{\xi}f_j (x+h, Du(x)) , \eta^2 \tau_h Du \rangle \, dx \notag \\
&+ \int_{\Omega} \langle D_{\xi}f_j (x+h, Du(x))- D_{\xi}f_j (x, Du(x)) , \eta^2 \tau_h Du \rangle \, dx \notag \\
& =: A_1 + A_2 + A_3 + A_4.
\end{align}So we have from Equation 4.5
\begin{align}
A_3 \leq |A_1| + |A_2| + |A_4| + |B| .
\end{align} We start by considering the term 
$|A_1|$. By using property (iv) of Lemma 3.1, Young’s and Hölder’s equalities, we get
\begin{align}
|A_1| &= \left| \int_{\Omega} \sum_{i, \alpha} \left( f_{\xi_{i}^{\alpha}}(x+h, Du(x+h))-f_{\xi_{i}^{\alpha}}(x+h, Du(x))\right) 2 \eta \, \eta_{x_i} \tau_h u^{\alpha}dx \right| \notag \\
& \leq 2 l_1\int_{\Omega} \eta \, |D\eta| \left( 1 + |Du(x+h)|^2 + |Du(x)^2| \right)^{\frac{q-2}{2}} \, |\tau_h
Du| \, |\tau_h u| dx \notag \\
& \leq \varepsilon \int_{\Omega} \eta^2 \, \left( 1 + |Du(x+h)|^2 + |Du(x)^2| \right)^{\frac{p-2}{2}} \, |\tau_h
Du|^2 dx \notag \\
& \qquad + c_\varepsilon \int_{\Omega} |D \eta|^2 \left( 1 + |Du(x+h)|^2 + |Du(x)^2| \right)^{\frac{2q-p-2}{2}} \, |\tau_h u|^2 dx\notag \\
& \leq \varepsilon \int_{\Omega} \eta^2 \, \left( 1 + |Du(x+h)|^2 + |Du(x)^2| \right)^{\frac{p-2}{2}} \, |\tau_h
Du|^2 dx \notag \\
& \qquad + \frac{c_\varepsilon}{(t-s)^2} \left( \int_{B_{t'}}(1+ |Du|^2)^{\frac{2q-p}{2}} dx \right)^{\frac{2q-p-2}{2q-p}} \, \left( \int_{B_{t}} |\tau_h u|^{2q-p} dx \right)^{\frac{2}{2q-p}},
\end{align} where in the last inequality we used that 
$|D \eta| \leq \frac{c}{t-s}$. Next applying Lemma 2.5, we have
\begin{align*}
|A_1| & \leq \varepsilon \int_{\Omega} \eta^2 \, \left( 1 + |Du(x+h)|^2 + |Du(x)^2| \right)^{\frac{p-2}{2}} \, |\tau_h
Du|^2 dx \notag \\
& \qquad + \frac{c_\varepsilon}{(t-s)^2} \, |h|^2 \, \left( \int_{B_{t'}}(1+ |Du|)^{2q-p} dx \right)^{\frac{2q-p-2}{2q-p}} \, \left( \int_{B_{t'}} |Du|^{2q-p}dx \right)^{\frac{2}{2q-p}} \notag \\
& \leq \varepsilon \int_{\Omega} \eta^2 \, \left( 1 + |Du(x+h)|^2 + |Du(x)^2| \right)^{\frac{p-2}{2}} \, |\tau_h
Du|^2 dx \notag \\
& \qquad + \frac{c_\varepsilon}{(t-s)^2} \, |h|^2 \, \int_{B_{t'}}(1+ |Du|)^{2q-p} dx ,\notag
\end{align*}where the last integral is finite by virtue of (Equation 4.3) and (Equation 4.4).
The ellipticity assumption (iii) of Lemma 3.1 implies
\begin{align}
A_3 \geq M_1 \int_{\Omega} \eta^2 \, \left( 1 + |Du(x+h)|^2 + |Du(x)^2| \right)^{\frac{p-2}{2}} \, |\tau_h
Du|^2 dx.
\end{align}From condition (vi) of Lemma 3.1, Lemma 2.5, the properties of η and Hölder’s inequality, we derive
\begin{align}
|A_2| &= \left| \int_{\Omega} \sum_{i, \alpha} \left( f_{\xi_{i}^{\alpha}}(x+h, Du(x))-f_{\xi_{i}^{\alpha}}(x, Du(x)) \right)2 \eta \, \eta_{x_i} \tau_h u^{\alpha} dx \right| \notag \\
& \leq 2 \int_{\Omega} \eta |D \eta||\tau_h u| |h|^{\gamma}(k(x+h)+ k(x))(1 + |Du|^2)^{\frac{q-1}{2}} dx \notag \\
& \leq \frac{c}{t-s}|h|^{\gamma} \int_{B_t} |\tau_h u| (k(x+h)+ k(x))(1 + |Du|^2)^{\frac{q-1}{2}} dx \notag \\
& \leq \frac{c}{t-s} \, |h|^{\gamma} \, \left( \int_{B_{R}} k^{r} dx \right)^{\frac{1}{r}} \, \left( \int_{B_{t}}|\tau_h u|^{\frac{r}{r-1}} (1 +|Du|)^{\frac{r(q-1)}{r-1}} dx \right)^{\frac{r-1}{r}} \notag \\
& \leq \frac{c}{t-s} \, |h|^{\gamma} \, \left( \int_{B_{R}} k^{r} dx \right)^{\frac{1}{r}} \, \left( \int_{B_{t}}|\tau_h u|^{\frac{rq}{r-1}}dx\right)^{\frac{r-1}{rq}} \left( \int_{B_t} (1 +|Du|)^{\frac{rq}{r-1}} dx \right)^{\frac{(r-1)(q-1)}{rq}} \notag \\
&\leq \frac{c}{t-s} \, |h|^{1+ \gamma} \, \left( \int_{B_{R}} k^{r} dx \right)^{\frac{1}{r}} \, \left( \int_{B_{t'}} (1 +|Du|)^{\frac{rq}{r-1}} dx \right)^{\frac{r-1}{r}},
\end{align} where we also used that 
$Du \in L^{\frac{rq}{r-1}}_{loc}(\Omega)$ by virtue of (Equation 4.3) and (Equation 4.4).
Similarly as for the estimate of A 2, we obtain
\begin{align}
|A_4| &= \left| \int_{\Omega} \sum_{i, \alpha} \left( f_{\xi_{i}^{\alpha}}(x+h, Du(x))-f_{\xi_{i}^{\alpha}}(x, Du(x)) \right) \eta^2 \, \tau_h u_{x_i}^{\alpha} dx \right| \notag \\
& \leq |h|^{\gamma} \int_{\Omega} \eta^2 |\tau_h Du| (k(x+h)+ k(x))( 1 + |Du|^2)^{\frac{q-1}{2}} dx \notag \\
& \leq \varepsilon \int_{\Omega} \eta^2 \, \left( 1 + |Du(x+h)|^2 + |Du(x)^2| \right)^{\frac{p-2}{2}} \, |\tau_h Du|^2 dx \notag \\
& \qquad + c_{\varepsilon} |h|^{2 \gamma} \int_{\Omega} \eta ^2 (k(x+h)+ k(x))^2 ( 1 + |Du|^2)^{\frac{2q-p}{2}} dx \notag \\
& \leq \varepsilon \int_{\Omega} \eta^2 \, \left( 1 + |Du(x+h)|^2 + |Du(x)^2| \right)^{\frac{p-2}{2}} \, |\tau_h Du|^2 dx \notag \\
& \qquad + c_{\varepsilon} |h|^{2 \gamma} \left( \int_{B_{R}} k^{r} dx \right)^{\frac{2}{r}} \, \left( \int_{B_{t}} (1 +|Du|)^{\frac{r(2q-p)}{r-2}} dx \right)^{\frac{r-2}{r}} ,
\end{align}where the last integral is finite by virtue of (Equation 4.3) and (Equation 4.4).
Now, we take care of B. Recalling that 
$p' \le m \le 2$, the Hölder’s inequality and Lemma 2.5 give that
\begin{align*}
\left| B\right|
& \le \left( \int_{B_{t'}} |g|^{m}dx \right)^\frac{1}{m} \left( \int_{B_{t'}} |\tau_{-h}(\eta^2 \tau_h u)|^{m'}\right)^\frac{1}{m'} \notag\\
&\leq c|h|\, \left( \int_{B_R} |g|^{m}dx \right)^\frac{1}{m}\,\left( \int_{B_{t^{\prime\prime}}} |D(\eta^2 \tau_h u)|^{m'} dx \right)^{\frac{1}{m'}} \notag \\
&\leq c|h|\, \left( \int_{B_R} |g|^{m}dx \right)^\frac{1}{m}\,\left( \int_{B_{t^{\prime\prime}}} \eta^{m'} |D \eta|^{m'} |\tau_h u|^{m'} dx \right)^{\frac{1}{m'}} \notag \\
& \ \ \ \ + c|h|\, \left( \int_{B_R} |g|^{m}dx \right)^\frac{1}{m}\,\left( \int_{B_{t^{\prime\prime}}} \eta^{2m'} |\tau_h Du|^{m'} dx \right)^{\frac{1}{m'}}. \notag
\end{align*} Using again Lemma 2.5, Young’s inequality and that 
$|D \eta | \le \frac{c}{t-s}$, 
$0 \le \eta \le 1$ and 
$2 \leq m' \leq p$, we have
\begin{align}
|B| &\leq \dfrac{c}{t-s}|h|^2\, \left( \int_{B_R} |g|^{m}dx \right)^\frac{1}{m}\,\left( \int_{B_{R}} |D u|^{m'} dx \right)^{\frac{1}{m'}} \notag \\
& \ \ \ \ + c_\varepsilon \ |h|^{m}\, \int_{B_R} |g|^{m}dx \notag\\
& \ \ \ \ + \varepsilon\int_{\Omega} \eta^{2} \left( 1+ |Du(x+h)|^2+|Du(x)|^2\right)^\frac{p-2}{2}|\tau_h Du|^2 dx .
\end{align}Inserting (Equation 4.8), (Equation 4.10), (Equation 4.9), (Equation 4.11) and (Equation 4.12) in (Equation 4.7), we infer
\begin{align*}
M_1 & \int_{\Omega} \eta^2 \, \left( 1 + |Du(x+h)|^2 + |Du(x)^2| \right)^{\frac{p-2}{2}} \, |\tau_h
Du|^2 dx \notag \\
& \qquad \leq 3 \varepsilon \int_{\Omega} \eta^2 \, \left( 1 + |Du(x+h)|^2 + |Du(x)^2| \right)^{\frac{p-2}{2}} \, |\tau_h
Du|^2 dx \notag \\
& \qquad +\frac{c_\varepsilon}{(t-s)^2} |h|^2 \int_{B_{t'}}(1 + |Du|)^{2q-p} dx \notag \\
& \qquad + \frac{c}{t-s} \, |h|^{1+ \gamma} \, \left( \int_{B_{R}} k^{r} dx \right)^{\frac{1}{r}} \, \left( \int_{B_{t'}} (1 +|Du|)^{\frac{rq}{r-1}} dx \right)^{\frac{r-1}{r}} \notag \\
& \qquad +c_\varepsilon \, |h|^{2 \gamma} \, \left( \int_{B_{R}} k^{r} dx \right)^{\frac{2}{r}} \, \left( \int_{B_{t'}} (1 +|Du|)^{\frac{r(2q-p)}{r-2}} dx \right)^{\frac{r-2}{r}} \notag \\
& \qquad + \dfrac{c}{t-s}|h|^2\, \left( \int_{B_R} |g|^{m}dx \right)^\frac{1}{m}\,\left( \int_{B_{R}} |D u|^p dx \right)^{\frac{1}{p}} \notag \\
& \qquad + c_\varepsilon \ |h|^{m}\, \int_{B_R} |g|^{m}dx. \notag
\end{align*} Choosing 
$\varepsilon = \frac{M_1}{6}$ and reabsorbing the first integral in the right-hand side by the left-hand side, we get
\begin{align*}
\int_{\Omega} \eta^2 \,& \left( 1 + |Du(x+h)|^2 + |Du(x)^2| \right)^{\frac{p-2}{2}} \, |\tau_h Du|^2 dx \notag \\
&\leq \frac{c}{(t-s)^2} |h|^2 \int_{B_{t'}}(1 + |Du|)^{2q-p} dx \notag \\
& \, + \frac{c}{t-s} \, |h|^{1+ \gamma} \, \left( \int_{B_{R}} k^{r} dx \right)^{\frac{1}{r}} \, \left( \int_{B_{t'}} (1 +|Du|)^{\frac{rq}{r-1}} dx \right)^{\frac{r-1}{r}} \notag \\
& \, +c \, |h|^{2 \gamma} \, \left( \int_{B_{R}} k^{r} dx \right)^{\frac{2}{r}} \, \left( \int_{B_{t'}} (1 +|Du|)^{\frac{r(2q-p)}{r-2}} dx \right)^{\frac{r-2}{r}} \notag \\
& \ + \dfrac{c}{t-s}|h|^2\, \left( \int_{B_R} |g|^{m}dx \right)^\frac{1}{m}\,\left( \int_{B_{R}} |D u|^p dx \right)^{\frac{1}{p}} \notag \\
& \ + c \ |h|^{m}\, \int_{B_R} |g|^{m}dx .\notag
\end{align*}Since q > p and r > 2, it holds
\begin{align*}2q-p \leq \frac{r(2q-p)}{r-2} \quad \text{and} \quad \frac{rq}{r-1} \leq \frac{r(2q-p)}{r-2}.\end{align*}Therefore, we have
\begin{align}
\int_{\Omega} \eta^2 \, & \left( 1 + |Du(x+h)|^2 + |Du(x)^2| \right)^{\frac{p-2}{2}} \, |\tau_h Du|^2 dx \notag \\
&\leq \frac{c}{(t-s)^2} |h|^2 \left( \int_{B_{t'}}(1 + |Du|)^\frac{r(2q-p)}{r-2} dx \right)^\frac{r-2}{r}\notag \\
& \, + \frac{c}{t-s} \, |h|^{1+ \gamma} \, \left( \int_{B_{R}} k^{r} dx \right)^{\frac{1}{r}} \, \left( \int_{B_{t'}} (1 +|Du|)^{\frac{r(2q-p)}{r-2}} dx \right)^{\frac{q(r-2)}{(2q-p)r}} \notag \\
& \, +c \, |h|^{2 \gamma} \, \left( \int_{B_{R}} k^{r} dx \right)^{\frac{2}{r}} \, \left( \int_{B_{t'}} (1 +|Du|)^{\frac{r(2q-p)}{r-2}} dx \right)^{\frac{r-2}{r}} \notag \\
& \ + \dfrac{c}{t-s}|h|^2\, \left( \int_{B_R} |g|^{m}dx \right)^\frac{1}{m}\,\left( \int_{B_{R}} |D u|^p dx \right)^{\frac{1}{p}} \notag \\
& \ + c \ |h|^{m}\, \int_{B_R} |g|^{m}dx .
\end{align}We now consider the interpolation inequality
\begin{align}
\left( \int_{B_{\rho}} |h|^\frac{r(2q-p)}{r-2}dx \right)^{\frac{r-2}{r(2q-p)}} \leq \left( \int_{B_\rho} |h|^p dx \right)^\frac{\pi}{p} \, \, \left( \int_{B_\rho} |h|^\frac{np}{n-\beta p} dx \right)^{\frac{(n-2 \beta)(1-\pi)}{np}},
\end{align} where 
$0 \lt \pi \lt 1$ is defined by
\begin{align*}
\dfrac{r-2}{r(2q-p)}= \dfrac{\pi}{p}+ \dfrac{(1-\pi)(n-2 \beta)}{np} \notag
\end{align*}which implies
\begin{align}
\pi= \dfrac{nr(p-q)-np+\beta r(2q-p)}{\beta r(2q-p)}, \ \ \ 1- \pi= \dfrac{n[r(q-p)+p]}{\beta r(2q-p)}.
\end{align}Applying Lemma 2.1 in the left-hand side of (Equation 4.13) and the interpolation inequality (Equation 4.14) in the right-hand side, we get
\begin{align*}
& \int_{B_{s}}|\tau_h V_p (Du)|^2 dx \le
\int_{B_{t}}\eta ^2 |\tau_h V_p (Du)|^2 dx \notag\\ & \leq \quad
\frac{c}{(t-s)^2} |h|^2 \left( \int_{B_{t'}}(1 + |Du|)^p dx \right)^{\frac{(2q-p) \pi}{p}} \notag\\ & \quad\cdot \left( \int_{B_{t'}}(1 + |Du|)^{\frac{np}{n-2 \beta}} dx \right)^{\frac{n-2 \beta}{np}(2q-p)(1- \pi)} \notag \\
& \quad \, + \frac{c}{t-s} \, |h|^{1+ \gamma} \, \left( \int_{B_{R}} k^{r} dx \right)^{\frac{1}{r}} \, \left( \int_{B_{t'}} (1 +|Du|)^p dx \right)^{\frac{\pi q}{p}} \cdot \notag \\
& \quad \qquad \qquad \cdot \left( \int_{B_{t'}}(1 + |Du|)^{\frac{np}{n-2 \beta}} dx \right)^{\frac{n-2 \beta}{np}q(1- \pi)} \notag \\
& \quad \, + c \, |h|^{2 \gamma} \, \left( \int_{B_{R}} k^{r} dx \right)^{\frac{2}{r}} \, \left( \int_{B_{t'}} (1 +|Du|)^p dx \right)^{\frac{\pi(2 q-p)}{p}} \cdot \notag \\
& \quad \qquad \qquad \cdot \left( \int_{B_{t'}}(1 + |Du|)^{\frac{np}{n-2 \beta}} dx \right)^{\frac{n-2 \beta}{np}(2 q- p)(1-\pi)} \notag \\
& \quad \ + \dfrac{c}{t-s}|h|^2\, \left( \int_{B_R} |g|^{m}dx \right)^\frac{1}{m}\,\left( \int_{B_{R}} |D u|^p dx \right)^{\frac{1}{p}} \ + \ c \ |h|^{m}\, \int_{B_R} |g|^{m}dx. \notag
\end{align*} where in the first inequality we used that η = 1 on Bs. Then, when 
$|h| \le 1$, since 
$\lambda:= \min \{\gamma, \frac{m}{2} \}$, we obtain
\begin{align}
\int_{B_{s}}|\tau_h V_p (Du)|^2 dx \le &
|h|^{2 \lambda } \Biggl\{
\frac{c}{(t-s)^2} \left( \int_{B_{t'}}(1 + |Du|)^p dx \right)^{\frac{(2q-p) \pi}{p}} \notag \\
& \, \cdot \left( \int_{B_{t'}}(1 + |Du|)^{\frac{np}{n-2 \beta}} dx \right)^{\frac{n-2 \beta}{np}(2q-p)(1- \pi)} \notag \\
& \, + \frac{c}{t-s} \, \left( \int_{B_{R}} k^{r} dx \right)^{\frac{1}{r}} \, \left( \int_{B_{t'}} (1 +|Du|)^p dx \right)^{\frac{\pi q}{p}} \cdot \notag \\
& \qquad \qquad \cdot \left( \int_{B_{t'}}(1 + |Du|)^{\frac{np}{n-2 \beta}} dx \right)^{\frac{n-2 \beta}{np}q(1- \pi)} \notag \\
& \, + c \, \left( \int_{B_{R}} k^{r} dx \right)^{\frac{2}{r}} \, \left( \int_{B_{t'}} (1 +|Du|)^p dx \right)^{\frac{\pi(2 q-p)}{p}} \cdot \notag \\
& \qquad \qquad \cdot \left( \int_{B_{t'}}(1 + |Du|)^{\frac{np}{n-2 \beta}} dx \right)^{\frac{n-2 \beta}{np}(2 q- p)(1-\pi)} \notag \\
& \ + \dfrac{c}{t-s}\, \left( \int_{B_R} |g|^{m}dx \right)^\frac{1}{m}\,\left( \int_{B_{R}} |D u|^p dx \right)^{\frac{1}{p}} \ + \ c \, \int_{B_R} |g|^{m}dx \Biggr\} .
\end{align} On the other hand, when 
$1 \lt |h| \le \frac{t'-t}{2}$, we have
\begin{align*}\int_{B_s} |\tau_hV_p(Du)|^2 dx & \le c \int_{B_{t'}} |V_p(Du)|^2 dx \le c \int_{B_{t'}}(1+|Du|)^p dx \nonumber\\ & \le c |h|^{2 \lambda} \int_{B_{t'}} (1+|Du|)^pdx.\end{align*} Combining Lemma 2.6 with the previous estimate, we infer, since 
$s-\rho=t'-t$,
\begin{align}
& \left( \int_{B_{\rho}}|Du|^\frac{np}{n-2 \beta} dx \right)^\frac{n-2 \beta}{n} \notag\\
& \quad \leq
\frac{c}{(t-s)^2(s-\rho)^{2(2\lambda+2\beta+2)}} \left( \int_{B_{t'}}(1 + |Du|)^p dx \right)^{\frac{(2q-p) \pi}{p}} \nonumber\\ & \quad \cdot \left( \int_{B_{t'}}(1 + |Du|)^{\frac{np}{n-2 \beta}} dx \right)^{\frac{n-2 \beta}{np}(2q-p)(1- \pi)} \notag \\
&\quad \, + \frac{c}{(t-s)(s-\rho)^{2(2\lambda+2\beta+2)}} \, \left( \int_{B_{R}} k^{r} dx \right)^{\frac{1}{r}} \, \left( \int_{B_{t'}} (1 +|Du|)^p dx \right)^{\frac{\pi q}{p}} \cdot \notag \\
& \qquad \qquad \cdot \left( \int_{B_{t'}}(1 + |Du|)^{\frac{np}{n-2 \beta}} dx \right)^{\frac{n-2 \beta}{np}q(1- \pi)} \notag \\
& \quad \, + \dfrac{c}{(s-\rho)^{2(2\lambda+2\beta+2)}} \, \left( \int_{B_{R}} k^{r} dx \right)^{\frac{2}{r}} \, \left( \int_{B_{t'}} (1 +|Du|)^p dx \right)^{\frac{\pi(2 q-p)}{p}} \cdot \notag \\
& \qquad \qquad \cdot \left( \int_{B_{t'}}(1 + |Du|)^{\frac{np}{n-2 \beta}} dx \right)^{\frac{n-2 \beta}{np}(2 q- p)(1-\pi)} \notag \\
& \quad \ + \dfrac{c}{(t-s)(s-\rho)^{2(2\lambda+2\beta+2)}}\, \left( \int_{B_R} |g|^{m}dx \right)^\frac{1}{m}\,\left( \int_{B_{R}} |D u|^p dx \right)^{\frac{1}{p}} \notag \\
&\quad \ + \dfrac{c}{(s-\rho)^{2(2\lambda+2\beta+2)}} \ \, \int_{B_R} |g|^{m}dx \,+\, \dfrac{c}{(s-\rho)^{2(2\lambda+2\beta+2)}} \int_{B_R}(1+|Du|)^p dx.
\end{align}We note that, since p < q, it holds
\begin{align*} \frac{p}{q(1- \pi)} \gt \frac{p}{(2q-p)(1- \pi)}.\end{align*}Moreover, recalling the definition of π at (Equation 4.15), we have
\begin{align*} \frac{p}{(2q-p)(1- \pi)} \gt 1
\Longleftrightarrow
\dfrac{p \beta r}{n [r(q-p)+p]} \gt 1
\Longleftrightarrow \frac{q}{p} \lt 1+ \frac{\beta}{n}-\frac{1}{r},\end{align*} which is true for every 
$\beta \in \left(n \left(\frac{q}{p}-1+\frac{1}{r}\right), \lambda\right)$, since condition (Equation 1.3) is in force. Hence, since both exponents 
$\frac{p}{q(1-\pi)}$ and 
$\frac{p}{(2q-p)(1-\pi)}$ are bigger than 1, we can apply Young’s inequality in the right-hand side of (Equation 4.17) and, assuming 
$\tilde{q}$ as the conjugate of 
$ \frac{p}{q(1-\pi)}$, 
$\tilde{p}$ as the conjugate of 
$\frac{p}{(2q-p)(1- \pi)}$, we obtain for a fixed 
$\varepsilon \in (0,1)$
\begin{align}
& \left( \int_{B_{\rho}}|Du|^\frac{np}{n-2 \beta} dx \right)^\frac{n-2 \beta}{n} \notag\\ & \qquad \leq \, \frac{\varepsilon}{3} \left( \int_{B_{t'}}(1+ |Du|)^{\frac{np}{n-2 \beta}} dx \right)^{\frac{n-2 \beta}{n}} \notag \\
& \qquad + \frac{c_\varepsilon}{(t-s)^{2\tilde{p}}(s-\rho)^{2(2\lambda+2\beta+2)\tilde{p}}} \left( \int_{B_{2R}}(1+ |Du|)^p dx\right)^{\frac{(2q-p)\pi}{p}\tilde{p}} \notag \\
& \qquad + \frac{\varepsilon}{3} \left( \int_{B_{t'}}(1+ |Du|)^\frac{np}{n-2 \beta} dx\right)^\frac{n-2 \beta}{n} \notag \\
& \qquad + \frac{c_\varepsilon}{(t-s)^{\tilde{q}}(s-\rho)^{2(2\lambda+2\beta+2)\tilde{q}}}\left( \int_{B_{R}} k^r dx \right)^{\frac{\tilde{q}}{r}} \left( \int_{B_{t'}}(1+ |Du|)^p dx \right)^{\frac{q \pi \tilde{q}}{p}} \notag \\
& \qquad + \frac{\varepsilon}{3} \left( \int_{B_{t'}}(1+ |Du|)^\frac{np}{n-2 \beta} dx \right)^\frac{n-2 \beta}{n} \notag \\
& \qquad + \dfrac{c_\varepsilon}{(s-\rho)^{2(2\lambda+2\beta+2)\tilde{p}}}\left( \int_{B_{R}} k^r dx \right)^{\frac{2\tilde{p}}{r}} \left( \int_{B_{t'}}(1+ |Du|)^p dx \right)^{\frac{(2q-p) \pi \tilde{p}}{p}} \notag \\
& \qquad + \dfrac{c}{(t-s)(s-\rho)^{2(2\lambda+2\beta+2)}}\, \left( \int_{B_R} |g|^{m}dx \right)^\frac{1}{m}\,\left( \int_{B_{R}} |D u|^p dx \right)^{\frac{1}{p}} \notag \\
& \qquad + \dfrac{c}{(s-\rho)^{2(2\lambda+2\beta+2)}} \, \int_{B_R} |g|^{m}dx \, + \, \dfrac{c}{(s-\rho)^{2(2\lambda+2\beta+2)}} \int_{B_R}(1+|Du|)^p dx,
\end{align} where the constants c and cɛ are independent of the radii 
$\rho, s, t, t', t^{\prime\prime}$.
Since
\begin{align*}t-s=s-\rho= \dfrac{t'-\rho}{3},\end{align*}estimate (Equation 4.18) becomes
\begin{align*}
& \left( \int_{B_{\rho}}|Du|^\frac{np}{n-2 \beta} dx \right)^\frac{n-2 \beta}{n} \nonumber\\ \quad & \leq \varepsilon \left( \int_{B_{t'}}(1+ |Du|)^\frac{np}{n-2 \beta} dx \right)^\frac{n-2 \beta}{n} \notag + \frac{c_\varepsilon}{(t'-\rho)^{2(2\lambda+2\beta+4)\tilde{p}}} \left( \int_{B_{R}}(1+ |Du|)^p dx\right)^{\frac{(2q-p)\pi}{p} \tilde{p}} \notag \\
& \qquad + \frac{c_\varepsilon}{(t'-\rho)^{(4\lambda+4\beta+5)\tilde{q}}}\left( \int_{B_{R}} k^r dx \right)^{\frac{\tilde{q}}{r}} \left( \int_{B_{R}}(1+ |Du|)^p dx \right)^{\frac{q \pi \tilde{q}}{p}} \notag \\
& \qquad + \dfrac{c_\varepsilon}{(t'-\rho)^{2(2\lambda+2\beta+2)\tilde{p}}} \left( \int_{B_{R}} k^r dx \right)^{\frac{2 \tilde{p}}{r}} \left( \int_{B_{R}}(1+ |Du|)^p dx \right)^{\frac{(2q-p) \pi \tilde{p}}{p}} \notag \\
& \qquad + \dfrac{c}{(t'-\rho)^{4\lambda+4\beta+5}}\, \left( \int_{B_R} |g|^{m}dx \right)^\frac{1}{m}\,\left( \int_{B_{R}} |D u|^p dx \right)^{\frac{1}{p}} \notag \\
& \qquad + \dfrac{c}{(t'-\rho)^{2(2\lambda+2\beta+2)}} \ \int_{B_R} |g|^{m}dx \ +\ \dfrac{c}{(t'-\rho)^{2(2\lambda+2\beta+2)}} \int_{B_R}(1+|Du|)^pdx . \notag
\end{align*}Applying Lemma 2.2
\begin{align*}
\left( \int_{B_{R/4}}|Du|^\frac{np}{n-2 \beta} dx \right)^\frac{n-2 \beta}{n} & \leq c_1 \left( 1+\int_{B_{R }} (1 + |Du|^p) dx + \Vert k \Vert_{L^r (B_{R})} + \Vert g \Vert_{L^{m}(B_R)}\right)^{\sigma_1}
\end{align*} which is finite since 
$u \in W^{1,p}_{loc}(\Omega, \mathbb{R}^N)$, 
$k \in L^r(\Omega)$ and 
$g \in L^{m}_{loc}(\Omega, \mathbb{R}^N)$ and so
\begin{align}
\int_{B_{R/4}} |\tau_h V_p (Du)|^2 dx \leq c_2|h|^{2 \lambda} \left( 1+\int_{B_{R}} (1 + |Du|^p) dx + \Vert k \Vert_{L^r (B_{R})} + \Vert g \Vert_{L^{m}(B_R)} \right)^{\sigma_2},
\end{align} for positive constants 
$c_1= c_1(n,p,q,m,\gamma,r, M,l,R, \beta)$ and 
$\sigma_1= \sigma_1 (n,p,q,m,\gamma,r, \beta)$ independent of j. The same happens to c 2 and σ 2, but they do not depend on β.
5. Proof of Theorem 1.2
In the case of p-growth conditions, i.e. when p = q, we have the following higher differentiability result.
Theorem 5.1 Let 
$g \in L^{m}_{loc}(\Omega, \mathbb{R}^N)$, with 
$p' \leq m \leq 2$. Moreover, assume that f satisfies (Equation F1)–(Equation F4) for exponents 
$2\leq p =q \lt \frac{n}{{\lambda}}$, where 
$\lambda:= \min \{\gamma, \frac{m}{2} \}$, and for a function 
$k \in L^\infty_{loc}(\Omega)$. If 
$u\in W^{1,p}_{loc}(\Omega, \mathbb{R}^N)$ is a local minimizer of (Equation 1.1), then 
$V_p(Du) \in B^\lambda_{2,\infty,\textrm{\text{loc}}}(\Omega, \mathbb{R}^{N \times n})$.
Proof. The proof goes similarly as the one of Theorem 4.1. Indeed, one can easily check that, under the stronger assumptions p = q and 
$k \in L^\infty_{loc}(\Omega)$, the exponent 
$\frac{r(2q-p)}{r-2}$ at (Equation 4.4) reduces to p. Therefore, we are able to obtain that 
$V_p(Du) \in B^\lambda_{2,\infty,\textrm{\text{loc}}}(\Omega, \mathbb{R}^{N \times n})$ with the estimate (Equation 4.19).
For a given ball 
$B_R \Subset \Omega$, we consider the variational problems
\begin{align}
\inf \left\{ \int_{B_R} [f_j(x,Dv)-g \cdot v] dx \ : \ v \in W^{1,p}_0(B_R, \mathbb{R}^N)+u \right\}, \quad j \in \mathbb{N},
\end{align} where 
$(f_j)_j$ is the sequence of functions defined in Lemma 3.1. The lower semicontinuity and the p-growth conditions of the integrand fj give the existence of a unique solution 
$u_j \in W^{1,p}_0(B_R, \mathbb{R}^N)+ u $ of the problem (Equation 5.1).
Now, fix a non-negative smooth kernel 
$\phi \in \mathcal{C}^\infty_0(B_1(0))$ such that 
$\int_{B_1(0)} \phi =1$ and consider the corresponding family of mollifiers 
$(\phi_\varepsilon)_{\varepsilon \gt 0}$. We set
\begin{align*}k_\varepsilon=k * \phi_\varepsilon\end{align*}and
\begin{align}
f_{j}^\varepsilon(x,\xi)= \int_{B_1(0)} \phi(y) f_j(x+\varepsilon y, \xi) dy
\end{align} on BR, for every 
$\varepsilon \lt \text{dist}(B_R,\Omega)$. One can easily check that 
$f_j^\varepsilon$ satisfies assumptions (Equation F1)–(Equation F3) and the set of conditions:
\begin{align}
1/L_1(|\xi|^p-1) \leq f_{j}^\varepsilon(x,\xi) \leq \bar{L}_1(j)(1+|\xi|^p),
\end{align}
\begin{align}
\langle D_\xi f_{j}^\varepsilon(x, \xi)-D_\xi f_{j}^\varepsilon(x, \eta),\xi-\eta\rangle \ge M_1 (1 +|\xi|^2+|\eta|^2)^\frac{p-2}{2}|\xi - \eta|^2,
\end{align}
\begin{align}
|D_\xi f_{j}^\varepsilon(x, \xi)-D_\xi f_{j}^\varepsilon(x, \eta)| \le l(j) (1 +|\xi|^2+|\eta|^2)^\frac{p-2}{2}|\xi - \eta|.
\end{align} By virtue of assumption (Equation F4), we have that for a.e. 
$x,y \in B_R$ and every 
$\xi \in \mathbb{R}^{N \times n}$,
\begin{align}
|D_\xi f_{j}^\varepsilon(x,\xi)-D_\xi f_{j}^\varepsilon(y, \xi)| \leq |x-y|^{\gamma} (k_\varepsilon(x)+k_\varepsilon(y)) (1 +|\xi|^2)^\frac{q-1}{2} ,
\end{align}
\begin{align}
|D_\xi f_{j}^\varepsilon(x,\xi)-D_\xi f_{j}^\varepsilon(y, \xi)| \leq \tilde{c}(j) |x-y|^{\gamma} (k_\varepsilon(x)+k_\varepsilon(y)) (1 +|\xi|^2)^\frac{p-1}{2}.
\end{align} Let us define the boundary value problem in 
$B_R \Subset \Omega$
\begin{align}
\inf \left\{ \int_{B_R} [f_j^\varepsilon(x,Dv)-g \cdot v] dx \ : \ v \in W^{1,p}_0(B_R, \mathbb{R}^N)+u_j \right\}.
\end{align} By the strict convexity of 
$f_j^\varepsilon$, there exists a unique minimum 
$u_j^\varepsilon \in W^{1,p}_0(B_R, \mathbb{R}^N)+ u_j $ of the problem (Equation 5.3). The p-growth from above in (Equation F1*) implies that 
$u_j^\varepsilon$ satisfies the Euler-Lagrange system
\begin{align}
\int_{B_R} \langle D_\xi f_j^\varepsilon(x,Du_j^\varepsilon),D \varphi \rangle dx = \int_{B_R} g \cdot \varphi dx, \quad \text{for every } \varphi \in W^{1,p}_0(B_R, \mathbb{R}^N).
\end{align} Exploiting the ellipticity assumption (Equation F2*), choosing 
$\varphi= u_j^\varepsilon -u_j$ as test function in (Equation 5.4) and (Equation 5.1), we have
\begin{align*}
M_1 & \int_{B_R} (1+|Du_j|^2+|Du_{j}^\varepsilon|^2)^{\frac{p-2}{2}}
|Du_{j}^\varepsilon-Du_j|^2 dx \notag\\
& \le \int_{B_R} \langle D_\xi f_j^\varepsilon(x,Du_j^\varepsilon) - D_\xi f_j^\varepsilon(x,Du_j),D (u_j^\varepsilon - u_j) \rangle dx
\notag\\
& = \underbrace{\int_{B_R} \langle D_\xi f_j^\varepsilon(x,Du_j^\varepsilon) ,D (u_j^\varepsilon - u_j) \rangle dx}_{= \ \int_{B_R} g \cdot (u_j^\varepsilon -u_j)dx} \notag\\
& \ \ \ \ - \int_{B_R} \langle D_\xi f_j^\varepsilon(x,Du_j),D (u_j^\varepsilon - u_j) \rangle dx \notag\\
& \ \ \ \ + \int_{B_R} \langle D_\xi f_j(x,Du_j),D (u_j^\varepsilon - u_j) \rangle dx \notag\\
& \ \ \ \ -\underbrace{\int_{B_R} \langle D_\xi f_j(x,Du_j),D (u_j^\varepsilon - u_j) \rangle dx}_{= \ \int_{B_R} g \cdot (u_j^\varepsilon -u_j)dx} \notag
\end{align*}and so
\begin{align*}
& \int_{B_R} (1+|Du_j|^2+|Du_{j}^\varepsilon|^2)^{\frac{p-2}{2}}
|Du_{j}^\varepsilon-Du_j|^2 dx \\
& \le c \int_{B_R} \langle D_\xi f_j(x,Du_j) - D_\xi f_j^\varepsilon(x,Du_j),D (u_j^\varepsilon - u_j) \rangle dx,
\end{align*} for a constant 
$c:=c(M,p)$. Applying Hölder’s inequality in the right hand side of the previous estimate, we infer
\begin{align*}
& \int_{B_R}
|Du_{j}^\varepsilon-Du_j|^p dx \notag\\
& \le c \int_{B_R} (1+|Du_j|^2+|Du_{j}^\varepsilon|^2)^{\frac{p-2}{2}}
|Du_{j}^\varepsilon-Du_j|^2 dx \notag\\
& \le c\left( \int_{B_R} | D_\xi f_j(x,Du_j) - D_\xi f_j^\varepsilon(x,Du_j)|^\frac{p}{p-1} dx \right)^\frac{p-1}{p} \left( \int_{B_R} |D (u_j^\varepsilon - u_j) |^p dx \right)^\frac{1}{p}, \notag
\end{align*} where the first inequality follows from the fact that 
$p \ge 2$. Next, we divide both sides by
\begin{align*} \left( \int_{B_R} |D (u_j^\varepsilon - u_j) |^p dx \right)^\frac{1}{p}\end{align*}thus getting
\begin{align}
\int_{B_R}
|Du_{j}^\varepsilon-Du_j|^p dx \le c \int_{B_R} | D_\xi f_j(x,Du_j) - D_\xi f_j^\varepsilon(x,Du_j)|^\frac{p}{p-1} dx .
\end{align}Assumptions (Equation F1*) and (Equation F2*) yields
\begin{align*}|D_\xi f_j^\varepsilon (x,Du_j)| \le c(j) (1+|Du_j|^{p-1}) \quad \text{a.e.\ in } B_R\end{align*} with c(j) independent of ɛ. Moreover, 
$|Du_j|^{p-1} \in L^{\frac{p}{p-1}}(B_R)$ and, by the very definition of 
$f_j^\varepsilon$, we have
\begin{align*}D_\xi f_j^\varepsilon (x,Du_j) \to D_\xi f_j (x,Du_j) \quad \text{a.e.\ in } B_R \text{\ \ \ as } \varepsilon \to 0^+\end{align*}and so the Dominated convergence theorem gives that
\begin{align*}\lim_{\varepsilon \to 0^+} \int_{B_R} | D_\xi f_j(x,Du_j) - D_\xi f_j^\varepsilon(x,Du_j)|^\frac{p}{p-1} dx =0. \end{align*} Therefore, passing to the limit as 
$\varepsilon \to 0^+$ in (Equation 5.5),
\begin{align}
u_j^\varepsilon \to u_j \text{\ strongly in \ } W^{1,p}(B_R, \mathbb{R}^N)
\end{align}and, up to subsequences,
\begin{align}
D u^\varepsilon_j \to Du_j \text{\ a.e.\ in \ } B_R .
\end{align} By the regularity Theorem 5.1, we have that 
$V_p(Du_j^\varepsilon) \in B^\lambda_{2,\infty,{{loc}}}(B_R, \mathbb{R}^{N \times n})$. Then, by virtue of Theorem 4.1, 
$u_j^\varepsilon$ satisfies the following estimate
\begin{align}
\int_{B_{\rho}} |\tau_h V_p (Du^\varepsilon_j)|^2 dx \leq C {|h|^{2 \lambda} }\left( {1+}\int_{B_{R}} (1 + |Du^\varepsilon_j|^p) dx + \Vert k \Vert_{L^r (B_{R})} + \Vert g \Vert_{L^{m}(B_R)} \right)^\sigma ,
\end{align} for all balls Bρ such that 
$B_{4\rho} \Subset B_R$, with constants 
$C = C(n,p,q,m,\gamma,r, M,l,R)$ and 
$\sigma= \sigma (n,p,q,m,\gamma,r)$, both independent of j and ɛ. Hence, by virtue of (Equation 5.7), we can apply Fatou’s Lemma, together with the continuity of the function 
$\xi \mapsto V_p(\xi)$ and of the operator τh, in estimate (Equation 5.8), and by using (Equation 5.6) and the strong convergence of kɛ to k in 
$L^r_{loc}(B_R)$, we obtain
\begin{align}
\int_{B_{\rho}} |\tau_h V_p (Du_j)|^2 dx \le &\liminf_{\varepsilon}\int_{B_{\rho}} |\tau_h V_p (Du^\varepsilon_j)|^2 dx\notag\\ \leq & \ C {|h|^{2 \lambda} }\left( {1+}\int_{B_{R}} (1 + |Du_j|^p) dx + \Vert k \Vert_{L^r (B_{R})} + \Vert g \Vert_{L^{m}(B_R)} \right)^\sigma.
\end{align}Now, by the growth assumption at (v) of Lemma 3.1, the minimality of uj and using u as test function, we get
\begin{align}
\int_{B_R} (|Du_j|^p-1)dx \le & \
L_1 \int_{B_R} f_j(x,Du_j)dx \notag \\
\le & \ L_1 \int_{B_R} [f_j(x,Du_j)-g \cdot u_j+ g \cdot u_j]dx \notag\\
\le & \ L_1 \int_{B_R} [f_j(x,Du)-g \cdot u]dx +L_1 \int_{B_R} g \cdot u_jdx \notag\\
\le & \ L_1 \int_{B_R} [f(x,Du)-g \cdot u]dx + L_1\int_{B_R} g \cdot u_jdx,
\end{align} where in the last inequality we used the monotonicity of 
$(f_j)$. From Poincaré inequality, we derive that
\begin{align*}
\Vert u_j -u\Vert_{L^p(B_R)} \leq c \Vert Du_j - Du \Vert_{L^p(B_R)}, \notag
\end{align*}that implies
\begin{align}
\Vert u_j \Vert_{L^p(B_R)} \le c(n,p,R) \left( \Vert Du_j \Vert_{L^p(B_R)}+\Vert u \Vert_{W^{1,p}(B_R)}\right)
\end{align}and so the second integral in the right hand side of (Equation 5.10) can be estimated as follows
\begin{align}
\int_{B_R} g \cdot u_jdx \le & \ \Vert g \Vert_{L^{p'}(B_R)} \Vert u_j \Vert_{L^p(B_R)} \notag\\
\le & \ c(n,p,R) \Vert g \Vert_{L^{p'}(B_R)} \left( \Vert Du_j \Vert_{L^p(B_R)}+\Vert u \Vert_{W^{1,p}(B_R)}\right) \notag\\
\le & \ \dfrac{1}{2L_1} \Vert Du_j \Vert_{L^p(B_R)}^p+ c(n,p,L_1,R) \left( \Vert g \Vert_{L^{p'}(B_R)}^{p'}+ \Vert u \Vert_{W^{1,p}(B_R)}^p \right),
\end{align} where we also used Hölder’s and Young’s inequalities. Putting the previous estimate in (Equation 5.10) and reabsorbing the term 
$\int_{B_R} |Du_j|^p dx$ from the right hand side by the left hand side, we get
\begin{align}
\int_{B_R} |Du_j|^pdx \le c \int_{B_R} [f(x,Du)-g \cdot u]dx + c \left( 1+\Vert g \Vert_{L^{p'}(B_R)}^{p'}+ \Vert u \Vert_{W^{1,p}(B_R)}^p \right),
\end{align} for a positive constant c independent of j. Therefore, up to subsequences, 
$(u_j)$ converges weakly in 
$W^{1,p}(B_R, \mathbb{R}^N)$ as 
$j \to \infty$ to a function 
$v \in W^{1,p}_0(B_R, \mathbb{R}^N)+u$.
The convergence of 
$(f_j)$ to f, the weak lower semicontinuity of 
$v \mapsto \int_{B_R} [f_j(x,Dv)-g \cdot v] dx$ and the minimality of uj give that
\begin{align*}
\int_{B_R} [f(x,Dv)-g \cdot v]dx \le & \liminf_{j_0} \int_{B_R} [f_{j_0}(x,Dv)-g \cdot v]dx \notag\\
\le & \liminf_{j_0} \liminf_j\int_{B_R} [f_{j_0}(x,Du_j)-g \cdot u_j]dx \notag\\
\le & \liminf_j\int_{B_R} [f_{j}(x,Du_j)-g \cdot u_j]dx \notag\\
\le & \liminf_j\int_{B_R} [f_{j}(x,Du)-g \cdot u]dx \notag\\
\le & \int_{B_R} [f(x,Du)-g \cdot u]dx,
\end{align*}which implies v = u a.e. in BR for the strict convexity of f.
Now, since uj satisfies (Equation 5.9) and estimate (Equation 5.13) holds, we have
\begin{align}
\int_{B_{\rho}} |\tau_h V_p (Du_j)|^2 dx \leq C |h|^{2 \lambda} & \left( {1+}\int_{B_{R}} (1 +f (x, D u))dx + \Vert k \Vert_{L^r (B_{R})} + \Vert g \Vert_{L^{m}(B_R)}\right. \nonumber\\ & \qquad \left. +\Vert u \Vert_{W^{1,p}(B_R)}^p\right)^\sigma,
\end{align} Using [Reference Hamburger25, Lemma 2.1] and 
$p \ge 2$, we get
\begin{align*}|H(\xi)-H(\eta)|^2 \le c(|\xi|^2+|\eta|^2)^\frac{p-2}{2}|\xi-\eta|^2 \le \tilde{c}|V_p(\xi)-V_p(\eta)|^2\end{align*}which implies, together with (Equation 5.14),
\begin{align*}
\int_{B_{\rho}} |\tau_h H(Du_j)|^2 dx \leq C |h|^{2 \lambda} & \left( 1+\int_{B_{R}} (1 +f (x, D u))dx + \Vert k \Vert_{L^r (B_{R})} + \Vert g \Vert_{L^{m}(B_R)} \right. \nonumber\\ & \qquad \left.+\Vert u \Vert_{W^{1,p}(B_R)}^p\right)^\sigma, \notag
\end{align*} where we denoted by 
$H(\xi):=|\xi|^\frac{p-2}{2}\xi$. The sequence 
$(H (Du_j))_j$ is bounded in 
$B_{2, \infty}^{\lambda}(B_{\rho},\mathbb{R}^{N \times n})$ for all balls Bρ such that 
$B_{4\rho} \Subset B_R$ and, therefore, is bounded in 
$B_{2, \infty,loc}^{\lambda}(B_{R}, \mathbb{R}^{N \times n})$ thanks to a standard covering argument. Besov spaces can be seen as interpolation spaces between Lp and 
$W^{1,p}$ (see [Reference Bergh and Löfström4, Chapter 6]), so applying classical interpolation theory ([Reference Bergh and Löfström4, Theorem 3.8.1]), there exists a function w such that
\begin{align*}H (Du_j) \rightharpoonup w \qquad \text{weakly in }
B_{2, \infty,loc}^{\lambda}(B_{R},\mathbb{R}^{N \times n})\end{align*}and
\begin{align*}H (Du_j) \to w \qquad \text{strongly in } L^2_{loc}(B_R,\mathbb{R}^{N \times n}) \text{and a.e.\ in } B_R.\end{align*}Recalling that
\begin{align*}|H(Du_j)|^2 = |Du_j|^p,\end{align*} we deduce that Duj converges strongly in 
$L^p_{loc}(B_R,\mathbb{R}^{N \times n})$. Hence, since
\begin{align*}u_j\rightharpoonup u \qquad \text{weakly in }
W_0^{1,p}(B_R,\mathbb{R}^{N })+u,\end{align*} we get, by uniqueness of weak limit, that 
$w= H(Du)$ and
\begin{align*}
Du_j \to Du \qquad \text{strongly in } L^p_{loc}(B_R,\mathbb{R}^{N \times n})
\end{align*}and, up to subsequences,
\begin{align}
D u_j \to Du \text{\ a.e.\ in \ } B_R .
\end{align} Then, by virtue of (Equation 5.15), we can apply Fatou’s Lemma, together with the continuity of the function 
$\xi \mapsto V_p(\xi)$ and of the operator τh, in estimate (Equation 5.14), thus getting
\begin{align*}
\int_{B_{\rho}} |\tau_h V_p (Du)|^2 dx &\leq \liminf_j \int_{B_{R/4}} |\tau_h V_p (Du_j)|^2 dx \notag \\
&\leq C |h|^{2 \lambda} \Big( {1+}\int_{B_{R}} (1 +f (x, D u)) dx + \Vert k \Vert_{L^r (B_{R})} \nonumber\\ & \qquad + \Vert g \Vert_{L^{m}(B_R)} +\Vert u \Vert_{W^{1,p}(B_R)}^p\Big)^\sigma, \notag
\end{align*} for every ball Bρ such that 
$B_{4 \rho} \Subset B_R \Subset \Omega$. Finally, we can conclude by a standard covering argument that 
$V_p(Du) \in B^\lambda_{2,\infty,loc}(\Omega,\mathbb{R}^{N \times n})$.
Acknowledgements
The authors sincerely thank Prof. A. Passarelli Di Napoli for proposing the research topic, engaging in insightful discussions, providing meticulous feedback, and generously dedicating her time. The authors wish to acknowledge the anonymous referee for her/his valuable suggestions. The authors are members of the Gruppo Nazionale per l’Analisi Matematica, la Probabilità e le loro Applicazioni (GNAMPA) of the Istituto Nazionale di Alta Matematica (INdAM). A.G. Grimaldi has been partially supported through the INdAM−GNAMPA 2024 Project “Interazione ottimale tra la regolarità dei coefficienti e l’anisotropia del problema in funzionali integrali a crescite non standard” (CUP: E53C23001670001). The authors have been supported through the INdAM - GNAMPA 2025 Project “Regolarità di soluzioni di equazioni paraboliche a crescita nonstandard degeneri” (CUP: E5324001950001). In addition S. Russo has also been supported through the project: Sustainable Mobility Center (Centro Nazionale per la Mobilità Sostenibile - CNMS) - SPOKE 10. In addition, A.G. Grimaldi has also been supported through the project: Geometric-Analytic Methods for PDEs and Applications (GAMPA) - funded by European Union - Next Generation EU within the PRIN 2022 program (D.D. 104 - 02/02/2022 Ministero dell’Università e della Ricerca). This manuscript reflects only the authors’ views and opinions and the Ministry cannot be considered responsible for them.
