Higher regularity in congested traffic dynamics

In this paper, we consider minimizers of integral functionals of the type F(u):=∫Ω[1p(|Du|-1)+p+f·u]dx\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\mathcal {F}}(u) := \int _\Omega \big [\tfrac{1}{p} \big (|Du|-1)^p_+ + f\cdot u\big ]\mathrm {d}x\nonumber \end{aligned}$$\end{document}for p>1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p>1$$\end{document} in the vectorial case of mappings u:Rn⊃Ω→RN\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u:{\mathbb {R}}^n\supset \Omega \rightarrow {\mathbb {R}}^N$$\end{document} with N≥1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N\ge 1$$\end{document}. Assuming that f belongs to Ln+σ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^{n+\sigma }$$\end{document} for some σ>0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma >0$$\end{document}, we prove that H(Du)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {H}}(Du)$$\end{document} is continuous in Ω\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega $$\end{document} for any continuous function H:RNn→RNn\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {H}}:{\mathbb {R}}^{Nn}\rightarrow {\mathbb {R}}^{Nn}$$\end{document} vanishing on {ξ∈RNn:|ξ|≤1}\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{\xi \in {\mathbb {R}}^{Nn} : |\xi |\le 1\}$$\end{document}. This extends previous results of Santambrogio and Vespri (Nonlinear Anal 73:3832–3841, 2010) when n=2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n=2$$\end{document}, and Colombo and Figalli (J Math Pures Appl (9) 101(1):94–117, 2014) for n≥2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n\ge 2$$\end{document}, to the vectorial case N≥1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N\ge 1$$\end{document}.


Introduction and main result
In this paper, we study C 1 -regularity of minimizers of integral functionals of the Calculus of Variations with widely degenerate convex integrands of the form where ⊂ R n , n ≥ 2, is a bounded domain and u : → R N , N ≥ 1, a possibly vector valued function. We concentrate ourself on the study of the prototype integrand for some p > 1. The datum f is required to belong to L n+σ for some σ > 0. The functional F with the specific integrand F from (1.2) is the prototype for a class of more general functionals where F is a convex function vanishing inside some convex set, and satisfying specific growth and ellipticity assumptions. For sake of clarity, the results in this paper are stated and proved for the functionals F(u) as in (1.1)-(1.2). However, we expect our techniques to apply to a general class of integrands with a widely degenerate structure as well. The functional F(u) and its associated Euler-Lagrange system div |Du| − 1 naturally arise in problems of optimal transport with congestion effects. In fact, minimizing (1.1) with N = 1 and the integrand from (1.2) is equivalent to the dual minimization problem min H(σ ) dx : σ ∈ L q ( , R n ), div σ = f , σ · ν ∂ = 0 , (1.4) where the integrand with H (t) = t + 1 q t q and 1 p + 1 q = 1 is the convex conjugate of F, or equivalently F = H * , and σ represents the traffic flow. The function g(t) = H (t) models the congestion effect. Note that σ → g(σ ) is increasing and g(0) = 1 > 0, so that moving in an empty street has nonzero cost. As shown in [4] the unique minimizer σ (x) of (1.4) is given by D ξ F(Du(x)). We refer to [2-4, 6, 7, 36] and the references therein for detailed motivations and for the physical meaning of the regularity of minimizers. It would be interesting to investigate if there are applications of the vectorial problem. However, our main motivation to consider the very degenerate system (1.3) was from a mathematical point of view.
In connection with congested traffic dynamic problems the regularity of minimizers, as well as the regularity of weak solutions of the associated autonomous Euler-Lagrange system has been an active field of research in recent years. For instance, in [2,4,9,24] Lipschitz regularity of minimizers has been established under suitable assumptions on the datum f .
At this point it is worthwhile to observe that, in general, no more than Lipschitz regularity can be expected for solutions of equations or systems as in (1.3). Indeed when f = 0, every 1-Lipschitz continuous function solves (1.3). On the other hand, in the scalar case N = 1 assuming f ∈ L n+σ for some σ > n, it was shown by Santambrogio and Vespri [36] for n = 2 and Figalli and Colombo [10,11] for n ≥ 2, that the composition of an arbitrary continuous function vanishing on the set {|∇u| ≤ 1} with ∇u is continuous.
Our aim in this paper is to investigate C 1 -regularity of minimizers in the vectorial case N ≥ 1. In general, regularity in the vectorial case is much more delicate and minimizers may be irregular although the integrand is smooth, cf. [15,38]. In this respect, regularity can be expected only for integrands with special structure. The first result in this direction has been obtained by Uhlenbeck [40] for the p-Laplace system when p ≥ 2. She proved that weak solutions are of class C 1,α . The scalar case had previously been established by Ural'ceva [41], while the case p ∈ (1, 2) was obtained by Tolksdorf [39]. As already mentioned, we cannot hope for a C 1,α -regularity result for the elliptic system (1.3), since Lipschitz continuity is optimal. However, we are able to establish in the vectorial setting that the composition K(Du) is continuous for any continuous function K : R N n → R vanishing on {ξ ∈ R N n : |ξ | ≤ 1}. This phenomenon is somewhat reminiscent of comparable results for the Stefan problem, in which the continuity of the energy cannot been shown, but the temperature shows a logarithmic type continuity, cf. [16,31].

Statement of the main result
Before formulating the main regularity result, we need to introduce a few notations. The natural energy space to deal with (local) minimizers of the integral functional F is the Sobolev space W 1, p ( , R N ). Then, (local) minimizers in W By carefully tracing the dependence of constants on the parameter δ in the proof of Theorem 3.6, one could determine an explicit modulus of continuity of G(Du), where G is defined in (2.2). However, it is not clear if G(Du) is Hölder continuous in general. For a different very degenerate elliptic equation a counterexample to Hölder continuity is provided in [11]. Theorem 1.1 can be regarded as the vectorial analog of the regularity results of Santambrogio and Vespri [36,Theorem 11] and of Figalli and Colombo [11,Theorem 1.1] as far as the model type integral functional is considered. The vectorial case cannot be treated with the methods from [11,36], since these are tailored to the scalar case. Nevertheless, some steps in our proof are similar, for example, the approximation procedure by some sequence of uniformly elliptic problems. The main difference in our proof is that we are establishing a variant of DiBenedetto's and Manfredi's proofs of C 1,α -regularity of minimizers to p-energy type functionals [17,32]. It is also inspired by the arguments from DiBenedetto and Friedman's pioneering results on C 1,α -regularity for parabolic p-Laplacian systems [18,20]. Roughly speaking, our strategy is the adaptation of De Giorgi's approach to the level of gradients in combination with Campanato type comparison arguments. The past has shown that De Giorgi's approach is extremely flexible. Therefore, we expect that our approach can be transferred to larger classes of widely degenerate functionals in the vectorial case. However, and in order to keep the individual steps as simple as possible, we limit ourselves to treating the model case.

Strategy of the proof
Concerning the overall strategy of proof a few words are in order. First, we observe that weak solutions of (1.5) are Lipschitz continuous. This has been proved in [2,4,9]. Moreover, functionals as in (1.1) fit into the broader context of asymptotically convex functionals, i.e. functionals having a p-Laplacian type structure only at infinity. This class of functionals has been widely studied, since the local Lipschitz regularity result by Chipot and Evans [8]. In particular we mention generalizations allowing super-and sub-quadratic growth [28,30,35], lower order terms [34]. Extensions to various other settings can be found in the non-complete list [12-14, 23-27, 37].
The proof of Theorem 1.1 is divided into several steps and starts by an approximation procedure. Indeed, by replacing h(t) by h ε (t) := h(t) + ε for ε > 0 and considering instead of (1.5) the Dirichlet-problem on a ball compactly contained in associated to the regularized coefficients h ε and with Dirchlet boundary datum u we obtain a sequence of approximating more regular mappings u ε . In particular, u ε has second weak derivatives in L 2 loc . In Sect. 3.1 we summarize the most important properties, i.e. uniform energy bounds, uniform quantitative interior L ∞ -gradient bounds, uniform quantitative higher differentiability W 2,2 -estimates, and finally strong L pconvergence of G δ (Du ε ) → G δ (Du) in the limit δ ↓ 0. The nonlinear mapping G δ : R N n → R N n with δ ∈ (0, 1] is defined by Observe that G δ vanishes on the larger set {|ξ | ≤ 1 + δ}. The reason for considering G δ is that on the complement of {|ξ | > 1 + δ} the system (1.5) behaves non-degenerate in the sense that the vector field A admits a uniform ellipticity bound from below, of course, with constants depending on δ. This point of view has already been exploited in [11,36]. As a first main result we prove that G δ (Du ε ) is Hölder continuous uniformly with respect to ε. However, the constants in the quantitative estimate, i.e. the Hölder exponent and the Hölder norm, may blow up when δ ↓ 0. We distinguish between two different regimes: the degenerate and non-degenerate regime. The degenerate regime is characterized by the fact that the measure of those points in a ball in which |G δ (Du ε )| is far from its supremum is large, while the non-degenerate regime is characterized by the opposite. In the non-degenerate regime we compare u ε with a solution of a linearized system. This allows us to derive a quantitative L 2 -excess-improvement for G 2δ (Du ε ) on some smaller ball (see Proposition 3.4). This step utilizes a suitable comparison estimate and the higher integrability of u ε . On the smaller ball we are again in the non-degenerate regime, so that the argument can be iterated yielding a Campanato-type estimate for the L 2 -excess of G 2δ (Du ε ). In the degenerate regime we establish that U ε := (|Du ε | − 1 − δ) 2 + is a subsolution to a linear uniformly elliptic equation with measurable coefficients; of course the ellipticity constants depend on δ and blow up as δ ↓ 0. At this stage a De Giorgi type argument allows a reduction of the modulus of G δ (Du ε ) on some smaller ball (see Proposition 3.5). However, on this smaller scale it is not clear whether or not we are in the degenerate or non-degenerate regime. Therefore one needs to distinguish between these two regimes again. In the non-degenerate regime we can conclude as above, while in the degenerate regime the reduction of the modulus of G δ (Du ε ) applies again. This argument can be iterated as long as we stay in the degenerate regime. However, if at a certain scale the switching from degenerate to non-degenerate occurs, the above Campanato type decay applies. If no switching occurs, we have at any scale of the iteration process a reduction of the modulus of G δ (Du ε ). This, however, shows that the supremum of |G δ (Du ε )|and hence also the one of |G 2δ (Du ε )|-on shrinking concentric balls converges to 0. Altogether this leads to a quantitative Hölder estimate for G 2δ (Du ε ) which remains stable under the already established convergence u ε → u as ε ↓ 0. The final step consits in passing to the limit δ ↓ 0 and conclude that G(Du) := (|Du|−1) + |Du| Du is continuous. This can be achieved by an application of Ascoli-Arzela's theorem. It is here where we loose control on the quantitative Hölder exponent. At this point the continuity of K(Du) for any continuous function K vanishing on {|ξ | ≤ 1} is an immediate consequence.

Notation
For the open ball of radius > 0 and center If the center is clear from the context we omit the reference to the center and write B respectively (v) for short. For the standard scalar product on Euclidean spaces R k as well as the space R kn of k × n matrices, we use the notation ξ · η. Finally, we use the notion ∇u for the gradient of a scalar function u, while we use Du for a vector field u.
Throughout this paper we abbreviate and Moreover, for δ ∈ (0, 1] we define and note that G 0 ≡ G. Generic constants are denoted by c. They may vary from line to line. Relevant dependencies on parameters and special constants will be suitably emphasized using parentheses or subscripts.

Algebraic inequalities
In this section, we summarize the relevant algebraic inequalities that will be needed later on. The first lemma follows from an elementary computation.
The next lemma can be deduced as in [29,Lemma 8.3].

Lemma 2.4
There exists a constant c = c( p) such that for any a > 1 and b ≥ 0 we have Proof We apply Lemma 2.2 with α = p − 1 > 0 to obtain This proves the claim.

Lemma 2.5
For a > 1 we have Moreover, for a, b > 1 there holds Proof By direct computation we have for a > 1 that , (2.4) from which the first claim immediately follows. We now turn our attention to the second claim. We may assume that 1 < a < b; otherwise we interchange the role of a and b. In view of (2.4) we find For the first term we have For the second term we estimate This completes the proof of the lemma.

Lemma 2.6
For any t ∈ R + we have Proof The first assertion can be achieved, since for t ∈ R + we have For the second claim, we first observe that both sides are zero for t ≤ 1. Therefore, it remains to consider t > 1. Recalling (2.4) we compute Moreover, we have so that Dividing the previous inequality by h(t) + h (t)t we infer that holds, proving the second claimed inequality.

Bilinear forms
For ε ∈ [0, 1] we define and Observe that all forms are symmetric in the arguments η and ζ . Due to the special structure of h ε and h ε the compositions A ε (Dv), B ε (Dv) and C ε (Dv) are well definined for v ∈ W 1, p as integrable functions. Therefore integral calculations involving these quantities make sense. The next Lemma provides the relevant ellipticity and boundedness properties of the bilinear forms A ε (ξ ), B ε (ξ ) and C ε (ξ ). The following abbreviations and λ(t) = 0 = (t) for 0 ≤ t ≤ 1 prove to be useful in the formulation of the Lemma.

Lemma 2.7
Let ε ∈ [0, 1] and ξ ∈ R N n \ {0}. The bilinear form A ε (ξ ) defined above satisfies The analogous estimates hold for the bilinear form B ε and any η, ζ ∈ R N n , as well as for C ε and any η, ζ ∈ R n .
Proof If |ξ | ≤ 1 the inequality holds trivially. Therefore, it remains to consider the case |ξ | > 1. We first establish the lower bound. If h ε (|ξ |) ≥ 0 we omit the second term in the definition of A ε and obtain while for h ε (|ξ |) < 0 we use Cauchy-Schwarz inequality and (2.5) to conclude Now, we turn our attention to the upper bound. If h ε (|ξ |) ≥ 0 we have This proves the claim for the bilinear form A ε . The corresponding estimates for B ε and C ε follow in the same way.
It should also be mentioned that the coercive symmetric bilinear forms fulfill Cauchy-Schwarz inequality. In particular, we have for any η, ζ ∈ R n .
In the next Lemma we put together the monotonicity and growth properties of the vector field A ε .
Proof The first inequality results from the following chain of inequalities From the third last line to the second last we used Lemma 2.4. For the proof of the second inequality we abbreviate Keeping this in mind we compute In turn we used Lemma 2.7 and the elementary inequality (|ξ s |−1) + ≤ |ξ s | ≤ |ξ |+|ξ |. Now, we distinguish whether or not |ξ | ≤ |ξ |. If |ξ | ≤ |ξ |, then For s ∈ 0, |ξ |−1 4|ξ | this implies a bound from below in the form holds true. In the case that |ξ | > |ξ |, we estimate |ξ s | from below by Therefore, for s ∈ 3|ξ |+1 4|ξ | , 1 we obtain This Inserting this above, we obtain the second claim of the Lemma.
In the following Lemma we quantify the remainder term in the linearization of A ε . In the application, it can be assumed that the linearization only takes place in points ξ with |ξ | > 1 in a quantifiable way. The precise statement is Proof We distinguish two cases. We start with the case |ξ −ξ | ≤ 1 Similarly to the computations in the proof of Lemma 2.8 we have This allows us to re-write We decompose and estimate the integrand appearing on the right-hand side and obtain in this way which immediately implies For the second term we use Lemma 2.1, Lemma 2.5, the assumptions on ξ,ξ and (2.9) to obtain Inserting the preceding estimates above, we find that At this stage it remains to consider the case |ξ −ξ | > 1 8 μ. Note that Therefore, by Lemma 2.7 and the assumption μ < 8|ξ − ξ | we obtain Joining both cases yields the claim.
Summing with respect to α ∈ {1, . . . , n}, we obtain In the preceding inequality we replace the term g(|Dv|) 2 by h(|Dv|)(|Dv| − 1) p + , which is possible by an application of the first inequality in Lemma 2.6. Moreover, we would also like to replace g (|Dv|) 2 |Dv| 2 by h (|Dv|)|Dv|(|Dv| − 1) p + . To this aim we have to distinguish two cases. If h (|Dv|) ≥ 0, we use the second inequality in Lemma 2.6 to replace g (|Dv|) 2 Thereby, we may omit the positive term g(|Dv|) 2 on the left-hand side. Inserting this above and using also Kato's inequality in the form |∇|Dv|| ≤ |D 2 v|, we obtain Otherwise, if h (|Dv|) < 0 we use Kato's inequality twice and again the second inequality in Lemma 2.6 to obtain This proves the asserted inequality also in the second case.

Proof of Theorem 1.1
In this section we will prove Theorem 1.1 under the hypothesis that Propositions 3.4 and 3.5 below are true. The remainder of the paper is then devoted to the proof of those two propositions.
Here and in the following we denote by u ∈ W 1, p ( , R N ) a weak solution of (1.5). We first observe that u ∈ W 1,∞ loc ( , R N ); cf [2,4] for the scalar case and [9] for the vectorial case. Therefore, we may always assume that Du is locally bounded in .

Regularization
The first step in the proof consists in the construction of more regular approximating solutions. To this aim we consider a fixed ball . We let ε ∈ (0, 1] and p := max{ p, 2} and recall the definition of the regularized coefficients A ε from (2.6).
The weak formulation of (3.1) is This can be retrieved similarly as in [33].
Our first observation is a uniform energy bound for Du ε .

Lemma 3.2 There exists a constant c = c(n, p) such that for any
Proof The desired estimate can be deduced with a standard argument by testing the weak formulation (3.2) with the test-function ϕ := u ε − u. Indeed, we have By Young's inequality we obtain for the first integral on the right-hand side The second integral on the right-hand side is estimated with Hölder's and Sobolev's inequality, so that where c = c(n, p). We insert these inequalities above and reabsorb the terms containing Du ε from the right-hand side into the left. In this way, we get with a constant c = c(n, p). The desired uniform energy bound can easily be deduced from the preceding inequality.
The next lemma ensures strong convergence of the approximating solutions, in the sense that G δ (Du ε ) strongly converges to G δ (Du) in L p .
Using Lemma 2.9 and Young's inequality we find Re-absorbing the first integral from the right into the left-hand side, we obtain The integral on the right-hand side is finite, since Du ∈ L ∞ loc ( , R N n ). Therefore, the preceding inequality implies strong convergence G δ (Du ε ) → G δ (Du) in L p (B R , R N n ) as ε ↓ 0.

Hölder-continuity of G ı (Du )
We recall the definition of G δ from (2.3). In this subsection we will prove that G δ (Du ε ) is locally Hölder continuous in B R for any δ ∈ (0, 1]. This will be achieved in Theorem 3.6. Thereby, it is essential that all constants are independent of ε. The proof of Theorem 3.6 relies on the distinction between two different regimesthe degenerate and non-degenerate regime. In the non-degenerate regime we will prove an excess-decay estimate for G δ (Du ε ) (see Proposition 3.4), while in the degenerate regime we establish a reduction of the modulus of G δ (Du ε ) (see Proposition 3.5). The precise setup is as follows. We consider a ball B R and denote by u ε the unique weak solution of the Cauchy-Dirichlet problem (3.1). For 0 < r o < R we let for some μ > 0 such that Next, for ν ∈ (0, 1) we define the super-level set of |Du ε | by The definition of the super-level set allows us to distinguish between the degenerate regime which is characterized by the measure condition and the non-degenerate regime which is characterized by the reversed inequality. Roughly speaking, in the degenerate regime the set of points with |Du ε | small has large measure, while in the non-degenerate regime the set of points with |Du ε | small is small in measure. We start with the latter one. In the following we abbreviate β := σ n+σ .
are satisfied, then the limit exists, and the excess decay estimate holds true. Moreover, we have The statement for the degenerate regime is as follows.
In turn, substituting 2δ by δ, this proves the claim of the proposition. We proceed in two steps.
Step 1: We prove that the limit holds true for a constant c = c(n, p, σ, f n+σ , M, δ). To this aim, we define for i ∈ N 0 radii i := 2 −i and μ i := κ i μ and observe that Now, suppose that assumption (3.10) holds on B (x o ). Then, Proposition 3.5 yields that either Note that the first alternative cannot happen, since it would imply Hence, we conclude that (3.4) holds on B 1 (x o ) with μ = μ 1 . If the measure condition (3.10) is satisfied with = 1 and μ = μ 1 , then a second application of Proposition 3.5 yields that either μ 2 As before, the first alternative cannot happen, since it would imply Assume now that (3.10) is satisfied for i = 1, . . . , i o − 1 up to some i o ∈ N, i.e. that (3.10) holds true on the balls B i (x o ) with μ = μ i . Then, we iteratively conclude that Now assume that (3.10) fails to hold for some i o ∈ N 0 . If μ i o > δ, the hypothesis of Proposition 3.4 are satisfied on B io (x o ) and we conclude that the limit Moreover, we have (3.14) Therefore, we obtain from the preceding inequality and (3.12) that holds true for any 0 < r ≤ i o . For a radius r ∈ ( i o , ] there exists i ∈ {0, . . . , i o } such that i+1 < r ≤ i . Using (3.13), (3.14) and (3.12) we obtain Combining the preceding two inequalities, we have shown (3.11) provided μ i o > δ.
Step 2: Here, we prove that the Lebesgue representative x → x of G 2δ (Du ε ) is Hölder continuous in B r o . The proof is standard once the excess decay (3.11) is established. For convenience of the reader we give the details. We consider x 1 , x 2 ∈ B r o . If r := |x 1 − x 2 | ≤ * we definex := 1 2 (x 1 + x 2 ) and obtain from Step 1 that This can be re-written in the form Otherwise, if r = |x 1 − x 2 | > * , we trivially have Together, this establishes that the Lebesgue representative x → x of G 2δ (Du ε ) is Hölder continuous in B r o with Hölder exponent α. Note that α admits the same dependencies as κ, i.e. α = α(n, N , p, σ, f n+σ , M, δ). This finishes the proof of the theorem.

Continuity of G(Du)
In this subsection it is important that all estimates are independent of ε. More precisely, constants might depend on δ, but are independent of ε.
Combining the estimates from above, we end up with Now, we consider the case |Du(x)| > 1 + √ ε. Here, Lemma 2.3 and (3.16) imply Hence, K(Du) is continuous on B r . Since B r B R were arbitrary, we have shown that K(Du) is continuous in . This completes the proof of the theorem.
As mentioned above, we have now finished the proof of the main theorem Theorem 1.1 under the condition that Propositions 3.4 and 3.5 are true. The rest of the paper is now devoted to the proof of those two propositions.

The main integral inequality for second derivatives
Throughout this subsection we assume as a general requirement that u ε : B R → R N is a weak solution to the regularized system (3.1). Instead of u ε , we write u for the sake of simplicity. In contrast, we will continue to use the subscript ε in the notation for the coefficients h ε and its associated bilinear forms, such as C ε . We recall that the bilinear forms have been defined in Sect. 2.3.
For some index β = 1, . . . , n we differentiate the regularized system (3.1) with respect to x β and obtain In (4.1) we choose the testing function ϕ = ζ φ(|Du|)D β u, where ζ ∈ C 1 0 (B R ) is non-negative and φ ∈ W 1,∞ loc (R ≥0 , R ≥0 ) is non-decreasing. Note that The resulting equations are then summed with respect to β from 1 to n. This leads to Now, we compute the right-hand side.
with the obvious meaning of R i . For the first term, we have by Young's inequality for any τ ∈ (0, 1). Similarly, we get Inserting this above and re-absorbing the terms containing second derivatives from the right-hand side into the left, we obtain for any non-negative function ζ ∈ C 1 0 (B R ). In the preceding inequality the parameter τ is at our disposal. We choose τ = 1 2 min{1, p − 1}. For the first factor, i.e. the term [. . . ], in the integrand of the first integral on the left-hand side we have In fact, if h ε (|Du|) ≥ 0 the inequality is obvious. If otherwise h ε (|Du|) < 0, which can only happen if p < 2 and |Du| > 1, the result follows by an application of Kato's inequality and (2.4). Indeed For the term in brackets of the second integral on the left-hand side a similar computation applies. The result of the calculation is Using the last and second last inequality above, we obtain an inequality which can be interpreted in two ways. On the one hand it can be seen as an energy inequality for the second derivatives of u. On the other hand-by discarding on the left-hand side the two non-negative terms containing second derivatives-the inequality implies that |Du| is a subsolution of an elliptic equation with measurable coefficients.

Subsolution to an elliptic equation
We start by showing that (|Du ε | − 1 − δ) 2 + is a sub-solution of a certain elliptic equation. More precisely holds true for any non-negative test function ζ ∈ C 1 0 (B (x o )) and with a universal constant c = c ( p, M, δ). The coefficients A αβ are given by where C ε is the bilinear form defined in (2.8).
Proof We apply Lemma 4.1 with φ(t) = (t − 1 − δ) + . Due to Lemma 2.7 the first two integrals on the left-hand side are non-negative and therefore can be discarded. In this way, we obtain with c = 4 min{1, p−1} that holds true for any non-negative ζ ∈ C 1 0 (B (x o )). To proceed further we compute The above inequalities together with the general assumptions (3.4), (3.5) allow to estimate Moreover, we have φ(|Du ε |)|Du ε | ≤ Mμ. In this way, we obtain for the right-hand side in (4.4) the estimate where c = c( p, M, δ). Now, we consider the left-hand side in (4.4). Observe that ∇U ε = 2φ(|Du ε |)∇|Du ε |, so that by the linearity of C ε (Du ε ) with respect to the first variable we have For the left-hand side this has the consequence that holds true. Here, we have taken into account the definition of the coefficients A αβ . Altogether we have shown the claim (4.3).

The coefficients A αβ in (4.3) can be explicitly written as
They are only degenerate elliptic due to the factor h ε (|Du ε |)|Du ε | which, for ε = 0, vanishes on the set {|Du| ≤ 1}. On the other hand U ε has its support in the set B R ∩{|Du ε | ≥ 1+δ}. This allows us to modify the coefficients on B R ∩{|Du ε | ≤ 1+δ}. This idea will lead us to an energy estimate for U ε in the next lemma.  function defined in (4.2). Then, for any k > 0 and any τ ∈ (0, 1) we have where c = c(n, p, M, δ) and β = σ n+σ .
Proof We let A αβ be the coefficients defined in Lemma 4.2 and extend them from The new coefficients A αβ are thus defined by From this definition and Lemma 4.2 we observe that U ε is a weak sub-solution also with the modified coefficients. More precisely, we have that for any non-negative ζ ∈ C 1 0 (B R ). We now investigate the upper bound and ellipticity of the coefficients A αβ . We will show that there exist 0 < λ ≤ < ∞ both depending at most on p, M and δ such that for any x ∈ B r o and ζ ∈ R n . We start with the former one. On the set where |Du ε | ≤ 1 + δ the upper bound holds with = 1, while on the set where |Du ε | > 1 + δ we have from Lemma 2.7 that A αβ ζ α ζ β = |Du ε |C ε (Du ε )(ζ, ζ ) which proves the claim with = ( p, M, δ) = M + pM p δ . Similarly, on the set where |Du ε | ≤ 1 + δ the ellipticity holds with λ = 1, while on the set where |Du ε | > 1 + δ we have from Lemma 2.7 that Now, the claimed energy estimate follows in a standard way by choosing in (4.5) a test-function of the form ζ = η 2 (U ε − k) + with a cut-off function η ∈ C 1 0 (B (x o )) with η ≡ 1 on B τ (x o ) and |∇η| ≤ 2 τ , cf. [19, Chapter 10.1].

Energy estimates
Here, we assume that the hypothesis of Proposition 3.4 are in force. Our starting point is again Lemma 4.1. This time we keep the two non-negative terms containing the quadratic forms A ε and C ε on the left-hand side. For any non-decreasing function φ ∈ W 1,∞ loc (R ≥0 , R ≥0 ) and any non-negative function ζ = η 2 ∈ C 1 0 (B (x o )) we have with the obvious meaning of I-III. For III we have Moreover, we estimate the integral I by Cauchy-Schwarz inequality and obtain From the second to last inequality we used the upper bound from Lemma 2.7 to estimate the second integral. Inserting the results above and re-absorbing the first term from the right into the left, we find that holds true with a constant c = c( p). We now choose φ(t) : For t ∈ [0, 1 + 2μ] we compute In turn we used δ ≤ μ ≤ M from (3.5) and (3.6), which implies on the one hand t ≤ 1 + 2μ ≤ 2 + 1 δ )μ, and on the other hand μ 2 ≤ max{δ 2− p , M 2− p }μ p . Next, we compute Due to assumptions (3.4) and (3.6) we know that |Du ε | ≤ 1 + 2μ on B (x o ). This allows us to use the preceding estimates in (4.6) to bound the right-hand side from above. Moreover, by Lemma 2.11 the left-hand side in (4.6) can be estimated from below. Proceeding in this way we obain for any η ∈ C 1 0 (B (x o )). The constant c depends only on p, M, and δ. Different concrete choices ofφ in (4.7) result in two important energy inequalities. The first one is for some universal constant c = c(n, p, M, δ).

This leads us to
with a constant c = c(n, p, M, δ), which is the claimed energy estimate.
The second energy estimate is for a constant c = c(n, p, M, δ).

Proof This time we choosẽ
in inequality (4.7), and obtain since ν ≤ 1 4 . Again, we choose η ∈ C 1 0 (B (x o )) to be a non-negative cut-off function with η ≡ 1 in B τ (x o ), 0 ≤ η ≤ 1, and |∇η| ≤ 2 (1−τ ) . This, together with the fact that |Du ε | ≤ 1 + 2μ on B (x o ), allows us to estimate the right-hand side in the above inequality. Indeed, we have Therefore it remains to estimate the left-hand side from below. The integral has to be taken only on the set of points We shrink this set to those points satisfying the stronger condition |Du ε ( On this set we have Inserting this above we conclude that holds true. This proves the claim.

The non-degenerate regime
The aim of this section is to prove Proposition 3.4. Throughout this section we presume the following general assumptions. For given ε ∈ (0, 1] we denote by u ε ∈ W 1, p (B R , R N ) the unique weak solution of the regularized system (3.1). Moreover, we assume that for some δ ∈ (0, 1] and μ > δ and a ball B 2 (x o ) ⊂ B r 1 with ≤ 1 assumptions (3.4)-(3.6) are in force. We denote by i.e. the L 2 -mean square deviation of Du ε from its mean value (Du ε ) x o , .

Higher integrability
An ingredient in the proof of Proposition 3.4 is the following higher integrability result.
We test the weak form (3.2) of the elliptic system by the testing function We use the monotonicity of A ε from Lemma 2.8 in order to estimate the first term from below. Due to our assumption on ξ and (3.4) we have |ξ | + |Du ε | ≤ 5|ξ | and therefore obtain To bound the second integral from above we use the structural upper bound from Lemma 2.8. Moreover, we observe that (|Du ε | − 1) + ≤ 4(|ξ | − 1) due to our assumption on ξ . This allows us to estimate for a constant c = c( p). Re-absorbing terms on the left-hand side, we find that where c = c( p). Due to the assumption on ξ , the assumption (3.5), and the particular choice of η, we conclude with an application of Hölder's and Sobolev-Poincaré's inequality a reverse Hölder inequality of the form

Comparison with a linear system
In this section we will consider the weak solution v ∈ u ε + W 1,2 0 (B /2 (x o ), R N ) of the linear elliptic system for any ϕ ∈ W 1,2 0 (B /2 (x o ), R N ) as comparison function to our solution u ε of the regularized elliptic system (3.1). Recall that B ε has been defined in (2.7).
Proof Throughout the proof we omit the reference to the center x o and write B instead of B (x o ). Moreover, we abbreviate ξ := (Du ε ) . Using the weak form (3.2) of the elliptic system we obtain for any ϕ ∈ W 1,2 0 (B /2 , R N ). Using also the fact that v is a weak solution of the linear elliptic system (5.2), we find that Here we used from the second to last line Lemma 2.10. This is possible since (3.4) and (5.3) are in force. Since u ε −v ∈ W 1,2 0 (B /2 , R N ), we may choose the testing function ϕ = u ε − v. Together with the bound from below from Lemma 2.7 and Hölder's and Poincaré's inequality this leads us to for a constant γ = γ (p, δ) > 0. We divide both sides by square the result and finally take means. This implies with a constant c = c(n, p, M, δ). Here we have also used δ < μ ≤ M. At this stage arrived, we want to reduce the integrability exponent on the right-hand side from 4 to 2(1 + ϑ), where ϑ = ϑ(n, p, M, δ) ∈ (0, min{ 1 2 , n+σ 2 − 1}] is the integrability exponent from the higher integrability Lemma 5.1. This is possible since |Du ε | and |ξ | are bounded by (2 + 1 δ )μ on account of (3.4) and (3.6). Then, the application of the higher integrability lemma yields where c = c(n, p, M, δ). Inserting this inequality above and noting that ≤ 1 finishes the proof of the lemma.
The following a priori estimate for solutions to linear elliptic systems can be inferred from [5] once the ellipticity conditions for the quadratic form B ε (Du ε ) x o , are established; see also [21,Theorem 2.3].  N , p, δ) such that for any τ ∈ (0, 1 2 Proof As mentioned before, the a priori estimate is standard. The constant c o depends on the dimensions n, N and the ellipticity constant and the upper bound of the quadratic form B ε (Du ε ) x o , . Due to assumption (5.3) and Lemma 2.7 these quantities only depend on p and δ.

Exploiting the measure theoretic information
The aim of this subsection is to convert the measure theoretic information (3.7) into a lower bound for the mean value of Du ε and smallness of the excess.
To prove (5.4) 1 , we first observe that the measure theoretic assumption (3.7) implies Hence, due to the definition of the set E ν , we obtain On the other hand, due to (5.4) 2 , we have Due to the choice of ν and the fact that n ≥ 2 we have ν ≤ ( 1 3 θ) . Together with the assumptions δ ≤ μ and θ ≤ 1 64 we obtain Inserting this above yields the claim (5.4) 1 and finishes the proof of the lemma.

Proof of Proposition 3.4
Our aim in this subsection is to prove Proposition 3.4. We start with an excess-decay estimate for the excess (x o , ) of Du ε .
Note that the constant c * depends on n, N , p, f L n+σ (B R ) , M and δ. For the proof of (II) 1 we use (5.6) 2 and τ n+2 ϑ ≤ τ n+2 to obtain Together with (5.6) 1 this implies (II) 1 . Now, we consider i > 1 and prove (I) i and (II) i assuming that (I) i−1 and (II) i−1 hold. From (I) i−1 and (II) i−1 we observe that the assumptions of Lemma 5.6 as formulated

The degenerate regime
Our aim in this section is to prove Proposition 3.5, which treats the degenerate regime. The proof relies on a De Giorgi type reduction argument reducing the supremum of U ε = (|Du ε | − 1 − δ) 2 + under the measure theoretic assumption (3.10). The starting point is the energy estimate for U ε from Lemma 4.3.
As in Section 5, we first formulate the general assumptions. For ε ∈ (0, 1] we denote by u ε ∈ W 1, p (B R , R N ) the unique weak solution to the Dirichlet problem (3.1) associated to the regularized system. We assume that (3.4) is in force for some μ, δ > 0 on some ball B 2 (x o ) ⊂ B r 1 B R . Let U ε := (|Du ε | − 1 − δ) 2 + denote the function defined in (4.2). Note that (3.4) implies Moreover, we set β := σ n+σ ∈ (0, 1). We start by a De Giorgi type lemma for U ε , which can for instance be deduced as in [19, Chap. 10, Proposition 4.1] by the use of the energy estimate from Lemma 4.3. For the readers convenience we provide the proof in the appendix Sect. 7. The proof of the next Lemma can be deduced as in [19, Chap. 10, Proposition 5.1] utilizing the energy estimate from Lemma 4.3; see also Sect. 7. Lemma 6.2 Assume that the general assumptions of Sect. 6 are in force and assume that (3.10) is satisfied for some ν ∈ (0, 1). Then, for any i * ∈ N we either have for a constant c * = c * (n, p, f n+σ , M, δ). Now, we have all the prerequisites at hand to provide the Proof of Proposition 3.5 Letν ∈ (0, 1) and c * be the constants from Lemmas 6.1 and 6.2. Note that both depend on n, p, f n+σ , M and δ. Choose i * ∈ N such that i * ≥ c * νν 2 .
Here we used that U ε − k i ≥ k i+1 − k i on A i+1 . Joining the preceding two inequalities and recalling the definition of Y i shows If μ 2 < β /θ , the lemma is proved. Otherwise, the term in brackets on the right-hand side is bounded by 2, so that where c = c(n, p, f n+σ , M, δ). Therefore, Lemma 7.1 ensures that Y i → 0 in the limit i → ∞, provided that Y 0 ≤ 2 −n 2 c − n 2 =:ν and hence U ε ≤ (1 − 1 2 θ)μ 2 on B /2 (x o ).
Proof of Lemma 6.2 For convenience in notation we omit the reference to the center x o . For i ∈ {0, . . . , i * } we consider the super-level sets where k i := (1 − 2 −i ν)μ 2 .
Then, due to assumption (3.10) we have |B \ A i | ≥ ν|B | for any i ∈ {0, . . . , i * }. Applying Lemma 7.2 to U ε on B with k = k i and = k i+1 , we obtain with a constant c = c(n) that In view of the energy estimate from Lemma 4.3 and the definition of k i we have with a constant c = c(n, p, f n+σ , M, δ). If μ 2 < 2 i β /ν, the lemma is proved. Otherwise, the term in brackets on the right-hand side is bounded by 2. Inserting this above yields