Boundary regularity results for minimisers of convex functionals with $(p,q)$-growth

We prove improved differentiability results for relaxed minimisers of vectorial convex functionals with $(p, q)$-growth, satisfying a H\"older-growth condition in $x$. We consider both Dirichlet and Neumann boundary data. In addition, we obtain a characterisation of regular boundary points for such minimisers. In particular, in case of homogeneous boundary conditions, this allows us to deduce partial boundary regularity of relaxed minimisers on smooth domains for radial integrands. We also obtain some partial boundary regularity results for non-homogeneous Neumann boundary conditions.


Introduction and results
We study minimisation problems of the form min u∈W 1,p (Ω,R m ) and min u∈W 1,p g D (Ω,R m ) Here f , g N and g D are sufficiently regular data. In the Neumann case (N) we will assume the compatibility conditionˆΩ where the integral is to be understood as being applied componentwise. In order to comment on our results, we state our assumptions on F precisely. Let 0 < α ≤ 1, n ≥ 2, m ≥ 1, and let Ω ⊂ R n be a bounded Lipschitz domain. We assume that F = F (x, z) : Ω × R m×n → R is measurable in x, continuously differentiable in z and moreover satisfies natural growth conditions and a Hölder-continuity assumption in x of the form for some µ, ν, Λ > 0, α ∈ (0, 1], all z, w ∈ R m×n and almost every x, y ∈ Ω, where 1 < p ≤ q. At times we additionally assume In particular, F (x, z) is convex. Further properties of such integrands will be discussed in Section 2.5. It is well-known that a major obstruction to regularity in the setting of (p, q)-growth is the possible occurrence of the Lavrentiev phenomenon. This describes the possibility that inf u∈W 1,p g (Ω) F (u) < inf u∈W 1,q g (Ω) F (u).
This phenomenon was first described in [53]. In the context of (p, q)-growth functionals the theory was expanded in [78], [76], [77]. Adapting the viewpoint and terminology of [15], we will not deal with pointwise minimisers of (N) and (D), but rather with relaxed minimisers in the following sense: Definition 1.1. We say u ∈ W 1,p (Ω) is a W 1,q -relaxed minimiser (usually referred to as a relaxed minimiser) of F N if u minimises the relaxed functional amongst all v ∈ X = W 1,p (Ω) where Y = W 1,q (Ω). Similarly, we say u ∈ W 1,p g D (Ω) is a W 1,q -relaxed minimiser of F D , if u minimises the relaxed functional amongst all v ∈ X = W 1,p g D (Ω) where Y = W 1,q (Ω) ∩ W 1,p g D (Ω). Note that if u ∈ W 1,q (Ω), then the relaxed functional agrees with the pointwise definition, justifying the terminology.
There is an extensive literature on the Lavrentiev phenomenon, an overview of which can be found in [14], [36]. The phenomenon also arises in nonlinear elasticity [37]. If the Lavrentiev phenomenon can be excluded, then minimisers of the relaxed and pointwise functional agree. In general, under the assumptions of this paper, it is not known whether the Lavrentiev phenomenon can be excluded. We refer however to [50] for conditions under which Lavrentiev can be excluded, see also [32,31] for local versions of these conditions. We further remark that while our improved differentiability results for minimisers of (N) or (D) are stated for relaxed minimisers, a closer look reveals that they imply an improved differentiability result for pointwise minimisers with an argument or assumption excluding Lavrentiev. Finally, we point out that there is a different approach to the relaxed functional, which involves studying measure representations of it. We refer to [34,1] for results and further references in this direction.
The study of regularity theory for minimisers in the case p < q started with the seminal papers [54,55]. We do not aim to give a complete overview of the theory here, and we will instead focus on results directly relevant to this paper. We refer to [59] for a good overview and further references. In general, the study has been almost completely concentrated on (D) and in particular, unless explicitly mentioned, all references pertain to this setting. A particular focus of research have been the special cases of the doublephase functional F (x, z) = |z| p + a(x)|z| q and functionals with p(x)-growth. For an introduction and further references with regards to these special cases we refer to the introduction of [6] and [28,62], respectively. Already in the scalar m = 1 autonomous case F (x, z) ≡ F (z) counterexamples show that in order to prove regularity of minimisers p and q may not be too far apart [38,54,44]. We stress that the counterexamples only apply to pointwise minimisers. We list the to our knowledge best available W 1,q locregularity results for general autonomous convex functionals with (p, q)-growth (when n ≥ 2): Under natural growth conditions it suffices to assume q < np n−1 [18] in order to obtain W 1,q loc -regularity of minimisers. To obtain the same conclusion under controlled growth conditions the gap may be widened to q < p 1 + 2 n−1 [65], and under controlled duality growth conditions it suffices to take q < np n−2 (if n = 2 it suffices to take q < ∞) [23]. We note that in all three cases higher integrability goes hand in hand with a higher differentiability result. We refer to [8] for results and references in the case of parabolic systems with (p, q)-growth.
The global theory is less developed and the only general results available extend the results of [18] up to the boundary in [49]. Additionally in [13], Lipschitz regularity up to the boundary is obtained for minimisers of scalar autonomous functionals satisfying nonstandard growth conditions and the structure condition F ≡ F 0 (|z|). The growth conditions considered include (p, q)-growth.
We now turn to the case of non-autonomous functionals F (x, z), convex and with (p, q)-growth in z, while satisfying a uniform α-Hölder condition in x. For n ≥ 2 counterexamples to W 1,q regularity with 1 < p < n < n + α < q are due to [32], see also [35]. Recent work suggests that the condition p < n < q may be removed [5]. If q < (n+α)p n , it was proven in [32] for many standard examples that minimisers enjoy W 1,q loc -regularity. Using [31], the result may be extended to functionals satisfying an additional condition on the x-dependence. W 1,q loc regularity is in general not known if q = (n+α)p n . An exception are functionals modeled on the double-phase functional [6], see also [24].
Concerning the theory of global improved differentiability and W 1,q -regularity re-sults for non-autonomous functionals in [49], the results of [31] are extended up to the boundary. For functionals satisfying a structural assumption inspired by the doublephase functional Caldéron-Zygmund estimates valid up to the boundary are obtained in [16]. Hölder-regularity up to the boundary for double-phase functionals is studied in [71].
If additional structure assumptions such as F ≡ F 0 (x, |z|) are imposed or if it is assumed that minimisers are bounded it is possible to improve on the results listed so far. Without going into further detail we refer to [11], [17] for results and further references in these directions. Local boundedness of minimisers for convex non-autonomous (p, q)-growth functionals under natural growth conditions and the additional assumption F (x, 2z) 1 + F (x, z) is studied in [43].
A secondary problem is that of partial C 1,α regularity for solutions, which is the best we can hope for in light of classical examples [26,57]. In the setting of (p, q)-growth, the first partial regularity results were obtained in [3] for functionals with a specific structure, and later in [61] for F ≡ F (z) subject to a controlled growth condition for exponents satisfying q ≥ p ≥ 2 and q < min p + 1, np n − 1 . (1. 2) The latter authors used the smoothing operator introduced in [34], and these ideas were developed further in [66] to obtain a Caccioppoli-type inequality in the quasiconvex setting. Building upon these techniques, the case of convex integrands F (x, z) was obtained in [27] assuming the exponents satisfy (1.2). While the above results only considered pointwise minimisers, the case of W 1,q -relaxed minimisers was treated in [68] through a delicate analysis of the relaxed functional. Recently, this was improved in [41] to the range 1 ≤ p < q < min p + 1, np n−1 . Integrands satisfying a p-Laplace type degeneration at the origin were treated in [67], and recently in this setting partial gradient continuity in the presence of a forcing term f ∈ L(n, 1) was obtained in [22,25].
The aforementioned results only consider interior regularity however, which raises the question of the existence of regular boundary points. For convex problems with regular growth this was settled in [51], where the integrand can additionally depend on u. This relied on a characterisation of regular boundary points obtained in [52] for quasiconvex integrands, however without an improved differentiability result it is unclear whether this characterisation is satisfied at any x 0 ∈ ∂Ω. To our knowledge this has not been considered in the case of (p, q)-growth, nor in the setting of Neumann boundary.
Concerning maximal improved differentiability subject to a forcing term, regarding the p-Laplace system, there is a wealth of information available, c.f. [60,12,19] and the references therein. For the purposes of this paper, we remark that a classical result in [69,70] establishes interior B s,p ∞ regularity with s = 1 + min p − 1, 1 with f ∈ L p ′ , and that this is sharp for p ≥ 2. Further, up to the boundary, [70] establishes that it is possible to take s = min (p − 1) 2 , 1 (p−1) 2 . Our results indicate that this likely is not sharp. If f ∈ W 1,p ′ , the obtained regularity improves to V p (Du) ∈ W 1,2 ; this was noted to hold globally in [30]. Finally, we remark that maximal improved Besovregularity for the p-Laplace equation was studied in [21].

Overview of results
Our first main result extends the results of [48] to the Neumann case: , and that u is a relaxed minimiser of F N in the class W 1,p (Ω). Then u ∈ W 1,q (Ω), and for any β < α the estimate holds for some γ > 0, where the implicit constant depends on Ω, λ, Λ, n, β, and p. In fact, we have If we additionally assume (H4), then we can obtain a further improvement in differentiability. Roughly speaking, once solutions lie in W 1,q they achieve the same improvement in differentiability as solutions to integrands with (p, p)-growth, up to possibly the endpoint scale. The argument uses second order difference quotients and a delicate iteration argument, and we will state it in the interior case where the result is already new. Theorem 1.3. Suppose 1 < p ≤ q < ∞, α ∈ (0, 1], and suppose F satisfies (H1)-(H4). Let q < n+α n p and f ∈ B α−1,p ′ ∞ (Ω). Suppose u is a relaxed minimiser for F . Then for all for β ∈ [α, 2], then we can take δ < min p ′ β 2 , α . Additionally, if a-priori u ∈ W 1,q (Ω), it suffices to ensure q < np n−α and (H3) can be omitted.
Restricting to radial integrands F ≡ F 0 (x, |z|) and homogeneous boundary conditions, adapting techniques developed in [30,29] we can extend this to hold globally.
In the Neumann case we can also infer some higher differentiability when g N ≡ 0; see Theorem 4.8.
Remark 1.5. The main result in [29] is stated to hold for non-radial integrands, however we believe that the proof is not correct. When restricting to radial integrands F ≡ F 0 (x, |z|), the proof can be corrected along the lines of the argument we give in Section 4.2 and the following statement can be shown to be true: (Ω) and V p,µ (Du) ∈ W 1,2 (Ω). The issue with the proof given in [29] lies in the proof of Lemma 3 of the paper. In the notation used there, it is not the case that Since v is the even extension of u, the normal derivative ∂ n v is odd, whereas ∂ n u was defined to be an even extension of ∂ n u. This also causes issues with the claim that J 03 = −J 04 . It is precisely these issues that force us to restrict to radial integrands in the above theorem.
Our results can also be adapted to piecewise C 1,1 domains, and to the case of mixed Dirichlet-Neumann boundary conditions, which we investigate in Section 4.4. We also remark that in the autonomous setting, it is possible to improve the range of q by employing arguments from [49] that are based on [65], extending the results of [18] from the interior case. Corollary 1.6. If F is independent of x, Theorems 1.2, 1.3 and 1.4 hold under the assumption 2 ≤ p ≤ q < max np n−1 , p + 1 . We remark that for the p-Laplace operator with forcing term f ∈ W −1+s,p ′ , 0 ≤ s ≤ 1, in [69], it is shown that solutions may in general fail to lie in W 1+ s max{2,p}−1 +ε,p loc (Ω) for any ε > 0. For p ≥ 2, this example has been adapted to Besov spaces in [75] showing in particular that, if f ∈ B −1+s,p ′ ∞ (Ω), solutions in general do not lie in B (Ω). An example in [10] shows that if p ≥ 2 and f ∈ W s,p ′ , then solutions may in general fail to satisfy V p,µ (Du) ∈ W p 2 s+1 p−1 +ε,2 loc (Ω). Thus our results are (up to excluding the edge case) optimal.
As a consequence of our improved differentiability results, we can establish the existence of regular boundary points. For this we will need to characterise the singular set of relaxed minimisers via a suitable ε-regularity result, which extends the corresponding interior result from [27].
, and the exponents satisfy We will assume further that Ω ⊂ R n is a bounded C 1,α domain and f ∈ L n 1−α (Ω, R m ). For the boundary function we will assume g D ∈ C 1,α (Ω, R m ) in the case of Dirichlet boundary, and g N ∈ C 0,α (Ω, R m ) in the case of Neumann boundary, in which case we also assume the compatibility condition (1.1). Suppose u is a relaxed minimiser of (N) or (D). Then for each M > 0 and 0 < β < α, there is ε > 0 and As a direct consequence of Theorem 1.7, we obtain the following characterisation of regular points. Here the singular set is understood to be the set of points x 0 where u fails to be C 1 in every neighbourhood of x 0 . Corollary 1.8 (Characterisation of regular points). Under the assumptions of Theorem 1.7, the singular set of u in Ω is given by The following is now immediate.
Finally, in the Neumann case, by Theorem 4.8 we obtain the existence of regular boundary points for inhomogeneous boundary conditions.
, then H n−1 -almost every boundary point is regular for u.
It would be interesting to obtain a statement akin to Theorem 1.11 in the case of Dirichlet boundary conditions. Remark 1.12. If a-priori u ∈ W 1,q (Ω), it suffices to ensure q < min np n−α , p + 1 in the statements of Corollary 1.9, Theorem 1.10 and Theorem 1.11. In particular in the case of homogeneous integrands F ≡ F 0 (|z|), 2 ≤ p ≤ q < min np n−1 , p + 1 suffices.
The outline of the paper is follows. In Section 2 we collect our notation and a number of preliminary results, concerning in particular Besov spaces and the relation of the relaxed functional to appropriate regularised versions of F . In Section 3 we consider a boundary version of a smoothing operator introduced in [34], where an additivity property for the relaxed function is proved allowing us to localise and flatten the boundary in our analysis.
In Section 4 we prove our differentiability results, starting with the basic result in the Neumann case in Section 4.1. We then consider the improved results, proving the interior statement in Section 4.2 where the key iteration argument is described. The corresponding global version is considered in Section 4.3 following the argument of [30,29], and in the Neumann case we must show an even extension of u preserves the desired fractional differentiability; this follows from Lemmas 4.6 and 4.7. We also explain how to apply these results to mixed boundary conditions in Section 4.4.
Finally in Section 5, we prove the characterisation of regular boundary points, by establishing a Caccioppoli-type inequality which we combine with a blow-up argument.

Notation
In this section we introduce our notation. Ω will always denote a open, bounded domain in R n . Given ω ⊂ R n , ω will denote its closure and χ ω the indicator function of ω. For subsets of R n we will use both L n and |·| to denote the Lebesgue measure, H s for the Hausdorff measure, and dim H for the Hausdorff dimension. We also denote the half-space in R n by R n + = {x ∈ R n : x n > 0}. For x ∈ R n and r > 0, B r (x) is the open Euclidean ball of radius r, centered at x, while S n−1 is the unit sphere in R n . We also write Ω r (x) = Ω ∩ B r (x). Further we denote by K r (x) the cube of side-length r, centered at x. We denote the cone of height ρ, aperture θ and axis in direction n by C ρ (θ, n). That is Here and elsewhere |·| denotes the Euclidean norm of a vector in R n and likewise the Euclidean norm of a matrix A ∈ R n×n . Id is the identity matrix in R n×n . Given an open set Ω, for λ > 0 set Here d(x, ∂Ω) = inf y∈∂Ω |x − y| is the distance of x from the boundary of Ω. Abusing notation, we write Ω h = Ω |h| for h ∈ R n .
If p ∈ [1, ∞], denote by p ′ = p p−1 the Hölder conjugate. The symbols a ∼ b and a b mean that there is some constant C > 0, depending only on n, m, p, q, Ω, µ, ν, Λ and α, and independent of a and b such that C −1 a ≤ b ≤ Ca and a ≤ Cb, respectively.
We will often find it useful to write for a function v defined on R n and a vector h ∈ R n , v h (x) = v(x+h). We also pick a family { ρ ε } of radially symmetric, non-negative mollifiers of unitary mass. We denote convolution of a locally integrable function u with For bounded and open sets Ω ⊂ R n , we will also use the notation − Ω f dx = 1 |Ω|´Ω f dx for averaged integrals. Function spaces will be denoted X(Ω) = X(Ω, R m ) where the target space may be dropped if it is clear from context, and for f ∈ X, g ∈ X ′ we will write´Ω f · g dx to be understood as the duality pairing.
At times, we need to consider versions of F N , F D and F N , F D localised to ω ⊂ R n . These are obtained by restricting the domain of integration and definition of the functions involved to ω and denoted F N (·, ω), F D (·, ω), F N (·, ω) and F D (·, ω), respectively. For the case of mixed boundary, and when localising in the Neumann case, we will need the relaxed functional defined for ω ⊂ R n , γ ⊂ ∂ω by

WB coverings, Lipschitz domains and extensions of F
We define Whitney-Besicovitch coverings, which combine properties of Whitney and Besicovitch coverings and were introduced in [46]. A nice presentation of the theory is given in [47]. To be precise: A family of cubes { K i } i∈I is called a Whitney-Besicovitch-covering (WB-covering) of Ω if there is a triple (δ, M, ε) of positive numbers such that The existence of a Whitney-covering for Ω is classical. The refinement to a WBcovering can be found in [47]: It will be of crucial importance to us that there exists a partition of unity associated to a WB-covering. Theorem 2.3 (cf. Theorem 3.19, [47]). Suppose the cubes { K i } i∈I form a WB-covering of Ω with constants (δ, M, ε). Then there is a family { ψ i } i∈I of infinitely differentiable functions that form a partition of unity on Ω with the following properties: It is well-known that Lipschitz domains satisfy a uniform interior and exterior cone condition. To be precise: Let Ω be a Lipschitz domain. Then (see e.g. [42, Section 1.2.2]): there are ρ 0 , θ 0 > 0 and a map n : R n → S n−1 such that for every x ∈ R n Finally, we extend F outside Ω. Let Ω ⋐ B(0, R). We extend F (x, z) to a functional on B(0, R) × R n×m still denoted F by setting It is straightforward to check that if F satisfies any of the properties (H1)-(H3) on Ω, F satisfies them on B(0, R) (after potentially increasing Λ). In Section 4.3 we will also need our extension to preserve (H4), however there we will flatten the boundary and take an even reflection in x instead.

Function spaces
We recall some basic properties of Sobolev and Besov spaces following the exposition in [64], alternatively the theory can be found in [72].
Denote by [·, ·] s,q the real interpolation functor. Let s ∈ (0, 1) and p, q ∈ [1, ∞]. We define (Ω) and that for 1 ≤ q < ∞, B s,p q (Ω) embeds continuously in B s,p ∞ (Ω). For α < 0, we define W α,p (Ω) as the dual of W −α,p ′ (Ω), and for s ∈ (0, 1) we can also define Note we also have W −s,p (Ω) = B −s,p p (Ω) by the duality theorem for real interpolation. We will use a characterisation of Besov spaces B s,p ∞ in terms of difference quotients. Let D be a set generating R n , star-shaped with respect to 0. For s ∈ (0, 1), p ∈ [1, ∞], consider Moreover there are positive constants C 1 , C 2 > 0 depending only on s, p, D, Ω such that If Ω = B r (x 0 ), then C 1 , C 2 are unchanged by replacing D with QD, where Q is an orthonormal matrix. In particular, when D = C ρ (θ, n) is a cone, they are independent of the choice of n.
Recall also that B s,p q (Ω) may be localised for s ∈ (0, 2) and 1 ≤ p, q ≤ ∞: We recall a well-known embedding theorem, see e.g. [73]: If q > q 1 , the embedding remains true provided the inequality (2.4) is strict.
We also recall the trace theorem in the following form, see e.g. [4]. Finally we recall the following well-known result (see [33] for the ingredients of the proof) which will justify extending u by extensions of g. We will later collect results concerning extensions preserving Besov-regularity in Section 4.3.2.

V -functions and inequalities
For p ∈ [1, ∞) and µ ≥ 0, we will introduce the V -functions If µ is clear from context, we will write V p,µ = V p . The fundamental property these functions satisfy is that noting that |V p,0 (z)| 2 = |z| p when µ = 0. This modular quantity naturally arises instead of the p-norm, which is the primary complication in the subquadratic (p ≤ 2) case. We will record several useful and well-known properties and estimates involving these functions.
We will later need to compare V -functions with different exponents, for which it will be useful to record that for p ≤ q. This follows by using (2.5) and distinguishing between the case of small and large |z| separately. We will also use the following Poincaré-Sobolev inequality adapted to V -functions. To apply these near the boundary it will be necessary to ensure the associated constants do not change, which is recorded in the following result.
In both cases the constant C depends on n, p, r/R only.
where we have used (2.9) in the last line, and the fact that |D| 1 n ∼ R ∼ r where the implicit constant depends on n and r/R. It remains to show the estimate Proof. Setting θ = s − r, consider the cones and choose finitely many We can also ensure that j χ U j ≤ M with M independent of s, t. Now we wish to show the corresponding Sobolev inequality on each U j , from which the result will follow. If B s (x 0 ) ⊂ Ω, this can be done for instance by parametrising B 1 by spherical coordinates, and noting this parametrisation restricts to a diffeomorphism on each U j by our choice of angle θ, and applying Lemma 2.8 there (noting this diffeomorphism maps each U j to a convex set). In the boundary case, our restriction on (x 0 ) n ensures each U j contains a coneŨ i associated to R n ) and becomes convex under this parametrisation. Thus we can also apply Lemma 2.8 in this case.
We will frequently pass from differentiability of the V -functional to differentiability of the function itself. The following is elementary, but will be useful in the sequel.
Proof. The case p ≥ 2 follows by noting that | · | p ≤ |V p,µ (·)| 2 using (2.5). If p < 2, for h ∈ R n we have by Hölder and (2.8) that from which the conclusion follows.

Bounds on the integrand
We collect some basic properties satisfied by the integrand F = F (x, z), that we will use extensively in what follows. We first note that for z ∈ R n×m and almost every x ∈ Ω, (2.13) By convexity of F, we have (H2) implies that (2.14) We also show that (H1) implies a more familiar ellipticity condition.
Lemma 2.11. Let 1 < p < ∞, µ ≥ 0, and suppose f : R m×n → R is C 1 and satisfies (H1). Then we have This follows from an explicit calculation. Now using the fundamental theorem of calculus we infer that for z, w ∈ R N n with tz and by approximation this holds for all z, w ∈ R m×n . For general f in C 1 , we will fix a standard mollifier ρ δ on R m×n and apply the above to F δ = F (x, ·) * ρ δ . Then we have and for |z| ≥ 2δ, we have where the implicit constant is independent of δ. Therefore by the C 2 case proved above, Sending δ → 0 we deduce that (2.15) holds for f and all z, w ∈ R m×n .

Properties of the relaxed and the regularised functional
In this section, we collect some results regarding the relaxed functional F N and its relation to a regularised version of F . For local versions of F N these results are well-known, c.f. Section 6 in [54] and [32]. Moreover, for global versions incorporating Dirichlet boundary conditions, similar results were obtained in [48]. The results in this section are minor modifications of those obtained in [48]. Define for u ∈ W 1,p (Ω) and ε ∈ (0, 1], is measurable for any z ∈ R n and F (x, ·) is continuously differentiable for almost every x ∈ Ω, it is easy to see that the same properties hold for F ε . Moreover, if F satisfies (H1)-(H3), then so does F ε with bounds that are independent of ε.
Minimisers of F ε (·) and F N (·) are related as follows: Lemma 2.12 (c.f. Lemma 6.4. in [54]). Suppose F (x, z) satisfies (H1) and (H2), and Proof. Note that by our assumptions´Ω f · u dx and´∂ Ω g N · u dH n−1 as defined for u ∈ W 1,q (Ω). Hence existence and uniqueness of u ε follows from the direct method and strict convexity, respectively. We further note that To prove the reverse implication note that for any v ∈ W 1,q (Ω), By definition of F N the inequality above extends to all v ∈ W 1,p (Ω). In particular, it holds with the choice v = u. Thus F ε (u ε ) → F N (u).
Using (H1) we may extract a (non-relabelled) subsequence of u ε , so that u ε ⇀ v weakly in W 1,p (Ω) for some v ∈ W 1,p (Ω). Note using our calculations above that v is a relaxed minimiser of F (·) in the class W 1,p g (Ω). Using (H1) it is easy to see that for Using the definition of F N and weak lower semicontinuity of norms, we see that this estimate extends to w 1 , w 2 ∈ W 1,p g (Ω). In particular, F N is convex and so u = v. Moreover the choice w 1 = u, w 2 = u ε in the estimate shows that u ε → u in W 1,p (Ω).
A useful criterion for establishing the equality F N (u) = F (u) is the following: We close this section by showing that relaxed minimisers are very weak solutions of the Euler-Lagrange system. We first note an elementary bound. Lemma 2.14. Suppose F satisfies (H1) and (H2). Then Proof. Introduce the Fenchel-conjugate of F , By the equality case of the Fenchel-young inequality ).

An elementary calculation shows that (H2) implies
This concludes the proof, recalling that we assumed also (H1).
We now adapt an argument from [23].
Proof. In the Neumann case, by the direct method we obtain u ε ∈ W 1,q (Ω) minimising F ε (·) in the class W 1,q (Ω) for each ε > 0. Without loss of generality we may assume that´Ω u ε dx = 0. Denote σ ε = ∂ z F ε (·, Du ε ) and µ ε = |Du ε | q−2 Du ε . Then u ε satisfies the Euler-Lagrange system Choose ψ = u ε and use Lemma 2.14, Hölder and Young's inequality to find that for any δ > 0, Hence choosing δ sufficiently small, we may rearrange to conclude Since ∂ z F (x, ·) is Carathéodory and Du ε → Du in measure by Lemma 2.12, we also have σ ε → F z (·, Du) in measure. Finally note that σ ε + εµ ε ⇀ F z (·, Du) in L 1 (Ω), so the latter is row-wise divergence free and since it is also in L q ′ (Ω), the result is established. The Dirichlet case is analogous, the only difference is that we test the equation against ψ = u ε − g.

Smoothing and relaxation
We will establish some technical results concerning the smoothing operator introduced in [34].
for any affine map a : R n → R n , where the implicit constants depend on n, p and Ω only. If in addition Ω = R n + and u ∈ W 1,p 0 (Ω), then w can be chosen to also lie in W 1,p 0 (Ω).
Proof. We will closely follow the proof of [34, Lemma 2.4], with a slight modification to the extension operator defined in [34,Lemma 2.2]. In our setting we can take η(x) = |x|, so for s ′ , t ′ to be determined we define By approximation this property extends in the sense of traces to all u ∈ W 1,p (Ω). We can also see that T preserves affine maps, that is T a = a. Now by the uniform cone condition, there is R 0 > 0 such that |Ω r (x)| ∼ |B r (x)| for all r ≤ R 0 and x ∈ Ω. Using this we can now argue exactly as in the interior case to obtain the corresponding estimates for p-norms. For the claimed estimates (3.1), (3.2) involving V -functions, we can argue as in [66, Lemma 6.5] when p ≤ 2, and if p ≥ 2 we simply note that V p,µ (z) ∼ |z| 2 + |z| p (see also [61,Lemma 2.3]).
We now assume that u ∈ W 1,p 0 (Ω). Using the first part, we obtain (3.2) are satisfied. It remains to modify w ′ so it vanishes on ∂Ω.
Let {Q i } be a WB-covering for A ′ using Theorem 2.2, and {ψ i } i∈I the associated partition of unity guaranteed by Theorem 2.3. We set We analyse the boundary behaviour of w on ∂A ′ , for which we assume Ω = R n + . By replacing B s \ B r by either B s \ B s+r 2 or B s+r 2 \ B r , we can assume Lemma 2.9 holds.
By a density argument (which is allowed since q ≤ np n−1 ) we may assume that v ∈ C(Ω).
In particular, extending w by u and using Lemma 2.6, we obtain w ∈ W 1,p 0 (Ω). Now using the additivity property (2.7), the finite overlap property (2.2) and (3.1), in A ′ , which follows by noting that a = i∈I a⋆ρ δ i ψ i . Similarly to above we can estimatê To estimate the latter term, we apply Lemma 2.8 followed by Jensen's inequality to n−1 , noting that ϕ(t) ∼ e 1,p (t). This gives,

Now combining the above estimates and summing we havê
where we have used Hölder's inequality for sequences and the finite overlap property in the last two lines. To complete the estimate, we apply Lemma 2.9 which gives and so (3.1) follows by the corresponding estimates for w ′ . For (3.2) a similar argument and so the result follows by the estimates for w ′ .
Remark 3.2. Inspecting the proof of Lemma 2.4 in [34], it holds that lim sup We now adapt the additivity property for relaxed functions established in [68, Lemma 7.7] to hold near the boundary. For the Neumann case we will use the mixed version F M as defined in (2.1).
Suppose Ω ⊂ R n is a bounded domain and that F satisfies (H2) and (2.13). Suppose u ∈ W 1,p (Ω) and letũ be a W 1,p -extension of u to R n such that ũ W 1,p (R n ) u W 1,p (Ω) . If the boundary condition holds, then we have with u k ⇀ u in W 1,p (Ω s (x 0 )). Extend each u k to a W 1,q function on R n denotedũ k in such a way that Due to [68,Lemma 7.7], there is a sequence Lemma 7.7] concerns autonomous integrands of the form F = F (z), an inspection of the proof reveals that we only need to assume the growth bounds (H2), (2.13) are satisfied. Hence abusing notation and using w k to denote the restriction of w k to Ω s (x 0 ), we deduce w k ⇀ũ weakly in W 1,p (Ω s (x 0 )) and further By a similar argument and using [68,Lemma 7 Composing the sequences u k and v k , we deduce The other direction is immediate from the definition.
We now turn to the Dirichlet case, where we need to modify the sequences {w k } and {v k } slightly in order to preserve the boundary condition w k = g D on ∂Ω ∩ B s (x 0 ) and v k = g D on ∂Ω \ B s (x 0 ) respectively. We only illustrate how to do so in the case of {w k }. For {v k } the argument is similar. Thus we assume that u ∈ W 1,p g D (Ω). Let Let {Q i } i∈I be a WB-partition of B s (x 0 ) and {ψ i } i∈I be a partition of unity associated to it. Set Arguing as in the proof of Lemma 3.1, we see that w k,ε ∈ W 1,q u (Ω s (x 0 )) and w k,ε → w k in W 1,q (Ω s (x 0 )) as ε → 0. Moreover, due to continuity of F in W 1,q , it holds that F (w k,ε , Ω s (x 0 )) → F (w k ) as ε → 0. By a diagonal subsequence argument, we thus obtain a sequence {w k } ⊂ W 1,q u (Ω s (x 0 )) such that w k ⇀ u weakly in W 1,p (Ω s (x 0 )) and This concludes the proof.

Reduction to a flat boundary
For many of our results we will assume that our boundary is locally flat, and consider the case where Ω = U + = U ∩ R n + for some bounded open subset U ⊂ R n containing B 1 (0). In this section we will outline how to reduce to this setting, so let Ω ⊂ R n be a bounded C 1,α domain, and let F = F (x, z) satisfy (H1)-(H4).
For x 0 ∈ ∂Ω let R > 0 be sufficiently small so there exists a bi-C 1,α mapping Ψ : here D τ denotes the vector of tangential derivatives along ∂Ω. The Dirichlet case is analogous, except that we drop the final term on the boundary.
Note that the associated integrand F (x, z) = F Φ(x), Dv · DΦ(x) −1 |det DΦ| inherits the properties (H1)-(H4) (up to possibly changing the constants appearing within them), and the Besov regularity of f is preserved. Additionally if F = F 0 (x, |z|) is radial, the same will hold for F .
We now claim that if u is a relaxed minimiser of F N or F D , by shrinking R if necessary we haveũ is a relaxed minimiser of F on U + . To achieve this we choose R > 0 so that (3.3) holds, using Remark 3.2. In the Neumann case, using (3.4) we have u minimises the mixed problem F M (u, Ω s (x 0 ), Ω ∩ ∂B s (x 0 )), whereas in the Dirichlet case we use (3.5) to see that u minimises F D (u, Ω s (x 0 )). Now composing the minimisers u ε of the regularised functional (2.16) with a flattening map Φ = Ψ −1 , we haveũ with F (v) as in (3.6). Now sending ε → 0 and noting v → v • Φ is weak-to-weak continuous in W 1,p , we deduce thatũ is a relaxed minimiser of F N , F D respectively on U + .
∂Ω) be such that (1.1) is satisfied. Suppose F satisfies (H1)-(H3) with 1 < p ≤ q < (n+α)p n , and that u is a relaxed minimiser of F N in the class W 1,p (Ω). Then u ∈ W 1,q (Ω), and for any β < α the estimate holds for some γ > 0, where the implicit constant depends on Ω, λ, Λ, n, β, and p. In fact, we have The proof follows the argument in [48] for the Dirichlet case closely, and we focus on the key differences.
Proof. We will write g = g N to simplify notation. The key is to prove the following a-priori estimate: Let v ε be the minimiser of F ε (·) in the class W 1,q (Ω). Note that v ε exists by the direct method and using the strict convexity deriving from (H1). Then for any 0 ≤ β < α, we claim there is γ > 0 such that the estimate holds, with the implicit constant independent of ε and γ. Let ρ 0 > 0 and n : R n → S n−1 be so that the uniform cone property (2.3) holds. Possibly reducing ρ 0 , assume without loss of generality that Ω + B 3ρ 0 (x) ⊂ B(0, R) for all x ∈ Ω. Here B(0, R) ⋑ Ω is the ball defined in Section 2.2. Given x 0 ∈ Ω, let 0 ≤ φ = φ x 0 ,ρ 0 ≤ 1 be a smooth cut-off supported in B 2ρ 0 (x 0 ) with φ(x) = 1 in B ρ 0 (x 0 ) and |D k φ(x)| ≤ C k ρ 0 −k for some C k > 0 and k ∈ N. Given functions v defined on R n and h ∈ R n introduce As in the Dirichlet case, for h ∈ C ρ 0 (θ 0 , −n(x 0 )), we note using Lemma 2.11 the lower bound The proof will conclude as in the Dirichlet case once we show that for every x 0 ∈ R n there is a constant C = C(n, ρ 0 , Λ, Ω) such that for all v ∈ W 1,q (Ω) We note that if p < 2, we combine (4.1), (4.2) and use the improved differentiability at the level of V -functions as is done in [32].
Indeed, given v ∈ W 1,q (Ω), letṽ be a W 1,q -extension of v to R n . Then considering h ∈ C ρ 0 (θ 0 , n(x 0 )), evidently T hṽ ∈ W 1,q (Ω) and so we can write The terms A 1 and A 2 are estimated exactly as in the Dirichlet case. For A 3 we make a small refinement and estimate noting this holds when α = 0, 1 and interpolating f → |A 3 | in [·, ·] 1−α,∞ for the general case. Concerning A 4 , we claim that for α ∈ (0, 1), we have We will establish this by similarly interpolating between α = 0, 1. When α = 0 we have where we have used Lemma 2.5. On the other hand, if α = 1, letg ∈ B 1,q ′ ∞ (Ω) be an extension of g to Ω. Denote for σ ∈ ∂Ω ∩ B 2ρ 0 , h σ (t) = σ − tn and put t max (σ) = max{t > 0 : σ − tn ∈ B 2ρ 0 }. Then, using the fundamental theorem of calculus, the coarea formula and the interior cone condition we have Note that we may apply the previous two estimates to v − (v) Ω instead, in order to replace the W 1,q norm by Dv L q (Ω) using the Poincaré inequality.
Improving the argument by optimising the cut-off function as demonstrated in [48], we obtain the following improved version in the autonomous case.
is autonomous and satisfies (H1), (H2). Then if u is a relaxed minimiser of F N in the class W 1,p (Ω), we have u ∈ W 1,q (Ω) and for any β < 1 the estimate (Ω) γ holds for some γ > 0. In fact, we have Remark 4.2 (Excluding Lavrentiev). We point out that under the assumptions of Theorem 1.2 (or the corresponding analog in the Dirichlet case, see [48]), we can rule out the Lavrentiev phenomenon provided we assume the following structural condition: we assume there is ε 0 > 0 such that for all 0 < ε < ε 0 and x ∈ Ω, there isŷ ∈ Ω ε (x) such that against all y ∈ Ω ε (x) and z ∈ R m×n . This condition was introduced in [31] in the context of functionals with (p, q)-growth and is similar to Assumption 2.3 in [77]. It has also been used in the context of lower semi-continuity in [2]. If this holds, we have that F N (u) = F N (u) for all u ∈ W 1,p (Ω) such that F (x, Du) ∈ L 1 (Ω) by a straightforward adaption of the argument for the Dirichlet case in [50].

Improved interior regularity
We now present an improved differentiability result in the interior case, which goes beyond the α/p differentiability obtained in [32] (and correspondingly Theorem 1.2 and [48] in the global case).
Remark 4.3. In light of the results in [69,75] our results are essentially sharp, however in the first part the limiting factor is the regularity of the forcing term f. Assuming regularity beyond f ∈ L p ′ gives further differentiability when p > 2, up until β = p − 1 where p ′ (β−2) 2 = 1. Moreover, we remark that in general the regularity in this case is also sharp, c.f. [10].
Suppose for 0 < r < s, B r ⋐ B s ⋐ Ω. Set r 1 = r + s−r 3 , r 2 = r + 2(s−r) 3 . Let φ be a radial cut-off with φ = 1 in B r , supported on B r 1 , and such that |∇ k φ| (s − r) −1 . Further set . We will split the proof into several steps by first establishing a basic estimate using difference quotients, describing an iteration procedure separately in the cases when p ≤ 2 and p > 2.
Step 1: Basic estimate: Consider h ∈ R n with |h| ≤ (s − r)/3. Testing the Euler-Lagrange equation for u with the test function −φ 2 ∆ 2 h u we find Concerning A, we note using discrete integration by parts Using (H1) via (2.14) and (H4) we have whereas for A 2 , A 3 by (2.14) and (H4) we have Similarly, we estimate B 1 as To estimate B 2 , as in Theorem 1.2 we consider the endpoint estimates in α. The case α = 0 follows by the duality estimate whereas in the case α = 1 we use the difference quotient estimate .
We will repeatedly use (4.5) to iteratively infer an improvement in differentiability, which will involve carefully estimating C 2 and C 3 at each step. Assume that We know by (a local version of) Theorem 1.2 that we can take δ = α 2 . By Lemma 2.10, we can estimate For the remaining two terms, suppose further that there is τ 1 , τ 2 ≥ 1 and ε > 0, σ > 0 such that B δ,2 Then by Lemma 2.10, since V p,µ (Du) ∈ B δ,2 ∞ (B s ) this implies that and so by Hölder we can estimate We can similarly bound C 3 , except we use that so we arrive at In particular, we can estimate (4.5) as Hence this shows that with the associated estimate This idea will be to iterate this procedure for a carefully chosen set of parameters δ, σ, τ 1 , τ 2 , ε. To do this, given 0 < r < s < 1 we set ρ k = r + 2 −k (s − r) for each k ≥ 0. We will apply the claim with ρ k+1 , ρ k in place of r, s for each k. We set δ 0 = α 2 and show by induction that V p,µ (Du) ∈ B δ k ,p ∞ (B ρ tk ), for some t = t(n, p, q, α) ≥ 1 where which tends to αp 2(p−1) as k → ∞. By Theorem 2.4, we can take Note that since q < np n−α and δ k ≥ α 2 , in the equality case, by shrinking ε > 0 we can assume that p 2 (τ 1 − ε) > q, so that we can take pτ 2 2 > 1. We will consider two cases depending on the possible values of τ 2 ; note that if τ 2 ≤ 2 then we may take σ = δ, and distinguishing between these two cases leads to a two-step iteration process. We consider the two possible cases at each step.
Case (a): Suppose that for some δ k we have τ 2 > 2 for all ε > 0 sufficiently small, which is equivalent to the condition Given ε > 0 to be determined, defining τ 1 , τ 2 to satisfy (4.6), (4.8) with equality, we will choose σ to satisfy for some κ 0 > 0, noting that q p < n+α n < n n−α . Thus we can choose ε > 0 so that Then by (4.7) we deduce that for some small c 0 > 0. We can now replace δ k = δ k,0 by δ k,1 and apply this iteratively, Since the second term is divergent, there is some j 0 = j 0 (α, n, p, q) ∈ N for which either δ k,j 0 = δ k+1 or δ k,j 0 ≥ n α(q−1) In the former case we are done, otherwise we can turn to case (b).
Step 3: Iteration, case p < 2. We will employ a similar argument, and we will only highlight the differences. As before assume there is δ > 0, k ≥ 0 such that V p,µ (Du) ∈ B 1+δ,2 ∞ (B s ) with the estimate (4.5). We will estimate (4.5) using this, where similarly as in the p ≥ 2 we use Lemma 2.10 to bound For C 2 we first estimate where we have used (2.8). We seek to estimate C 3 similarly, but we will use the fundamental theorem of calculus to write where we have used the fact that p − 2 ≤ 0 in the last line. Now suppose there is τ 1 , τ 2 ≥ 1 and ε, σ > 0 such that we have the embedding ∞ (B s ) as before, but we now require For the embeddings to hold, by Theorem 2.4 we require Now by Hölder we can estimate and there is γ > 0 such that To iterate this, we set δ 0 = α 2 and show by induction that which tends to α as k → ∞.
Choose τ 1 so that the first inequality in (4.10) holds with equality, then since δ ≥ α 2 and q < np n−α we have τ 1 ≥ 2n n−α . Therefore for ε > 0 sufficiently small we get so this ensures that we can choose 1 < τ 2 < 2q p so that (4.9) holds. As before, if the equality case gives τ 2 ≤ 2, then we can take τ 2 = 2 and σ = δ. Otherwise we have for all ε > 0 small. Determining τ 2 by the equality case of (4.9), by (4.10) 2 we choose σ = δ k + n 1 Since q < np n−α we have κ 0 > 0, and so we can choose ε > 0 sufficiently small, so that Then analogously as in the case p ≥ 2, we can use this iteratively to show that we have We can now iterate this exactly as in the case p ≥ 2.
Step 4: Case of regular f . Now suppose f ∈ B β−1,p ′ ∞ (Ω) with β ≥ α, and suppose that V p,µ (Du) ∈ B δ,2 ∞ (Ω). We will assume p > 2, as we otherwise see no improvement. We argue analogously as in Step 1, however we will use a different estimate for B 2 ; if β ≤ 1 we can argue as in (4.4) to estimate using Lemma 2.10 for the second inequality. Now assuming β > 1, we will write We can then estimate Now letting τ = min{1, δ + β − 1}, by interpolation along σ ∈ (0, 1) we will show that If σ = 0, by definition of the Besov seminorm we have If σ = 1 note that for τ ∈ (0, 1) by interpolation. Now interpolating in σ we deduce (4.11) for σ ∈ (0, 1) and τ ∈ [0, 1]. Note that we do not need the estimate in the endpoint case σ = 1, as it suffices to run this argument with σ = p−2 2 < 1 as noted in Remark 4.3. Also we have assumed f ∈ B β−1,p ′ 1 (Ω), however for any choice of δ < β p−1 we can findβ < β such that δ <β p−1 , and we can run the above argument knowing that f ∈ Bβ −1,p ′ 1 (Ω). Hence estimating the other terms in the same way as before, we arrive at the estimate as in the end of step 1, assuming that δ + β − 1 ≤ 1 (decreasing β if necessary). We also estimate the latter two terms slightly differently. For C 2 we use (2.8) to estimate For C 3 , we will write so we can estimate Now given V p,µ (Du) ∈ B δ,2 ∞ (B s ), suppose τ 1 , τ 2 ≥ 1 and σ, ε > 0 are such that the embedding B δ,2 ∞ (B s ) ֒→ L τ 1 −ε (B s ) ∩ B σ,τ 2 ∞ (B s ) and (4.9) holds, exactly as in Step 3. Then we can estimate Now we employ the same iteration argument as in the p ≥ 2 case from Theorem 1.3 to conclude. In fact, this gives rise to the iteration Now we claim we can take the iteration This follows by arguing as in Step 3, where an intermediate iteration allows us to replace σ with δ. Thus, taking k → ∞, we obtain the desired regularity improvement.

Improved boundary regularity for radial integrands
Our next goal is to obtain a version of Theorem 1.3 applicable up to the boundary, by adapting the approach in [30,29]. Our results will be restricted to radial integrands F ≡ F 0 (|z|), and to homogeneous boundary conditions in the Dirichlet case. More precisely, we will prove the following.
Following the procedure in Section 3.1 we can reduce to the case of a flat boundary. While this procedure only requires ∂Ω to be of class C 1,α , the C 1,1 regularity is assumed to ensure the fractional differentiability is preserved under this transformation. In particular, we will write Ω t = Ω ∩ B t (0), Γ t = B t ∩ ∂Ω, and suppose Ω = B + s , Γ t = B t ∩{x n = 0}. In the case of the Laplacian it is known that solutions in C 1 domains are in general no better than B 3 2 ,2 2 (Ω) [20], so additional regularity of ∂Ω is necessary, however we have not attempted to optimise this.

The boundary difference quotient argument
The key ingredient is that we require a version of (4.5) valid up to the boundary, and to achieve this we will flatten the boundary.
Here r 1 = r + s−r 3 and r 2 = r + 2(s−r) 3 . In the Dirichlet case, the term involving g N may be dropped. Moreover, if a-priori u ∈ W 1,q (Ω), it suffices to ensure q < np n−α .
A crucial part in the proof will be played by the following elementary lemma: Proof of Lemma 4.4. We denote a(x, t) = ∂ z F (x, t) and a 0 (x, t) = ∂ t F 0 (x, t)/t. As in Theorem 1.3, let φ be a radial cut-off with φ = 1 in B r , supported in B r 1 .
Step 1: Tangential directions. We will consider difference quotients in the direction h = he i with i ∈ {1, . . . , n − 1} and |h| ≤ s−r 3 . Suppose h = he i with i ∈ {1, . . . , n − 1} and |h| ≤ s−r 3 . Noting that φ 2 ∆ 2 h u is an admissible test function in both the Dirichlet and mixed case, there is no problem in applying the proof of Theorem 1.3 in order to obtain (4.5).
The only difference is that in the Neumann case we obtain an extra term which, analogously as in (4.3) from the proof of Theorem 1.2, we estimate as Step 2: Normal direction. Let now h = he n with |h| ≤ s−r 3 . Similarly as in the first step, we will test the system against φ∆ 2 hũ . In the Dirichlet case, we observe this is admissible owing to the fact thatũ is odd and that g D = 0. In both cases this gives which we can write as understanding that g N ≡ 0 in the Dirichlet case. For this we observe the identity Unlike in the tangential case the last two terms do not necessarily vanish, however by exploiting the symmetries of the system we will show that they cancel; the regions of integration of the respective terms are shown in Figure 1. For 1 ≤ i ≤ n, let a i denote the i th column of a, so that a i (x, Du) = a 0 (x, |Du|)∂ i u, as then we can write (4.12) We apply Lemma 4.5 with for each 1 ≤ i ≤ n and x ′ ∈ B n−1 s . Across the x n -axis we have ∂ i u is odd or even depending on whether i = n and the boundary condition, and hence a 0 (x, |Du|) is even. Also φ is even along the x n -axis since it is a radial cutoff, so it follows that σ and τ have the same parity. Hence, using this lemma in (4.12) and summing over n, we deduce that Figure 1: Now we can proceed as in the tangential case, where we write so by (2.14), (2.15), (H4) we can estimate and similarly which gives us our estimates for For the remaining terms we can as before estimate Therefore collecting our estimates we have shown the desired estimate.

Odd and even extensions
Lemma 4.4 indicates that we want to construct W 1,q -extensions of solutions that preserve fractional differentiability. Under zero Dirichlet boundary conditions, we will take an odd extension given bỹ Observe that since u vanishes on Γ s , by Lemma 2.6 it is weakly differentiable in B s . Moreover by [74,Theorem 4.5.2] for all δ ∈ (0, 1) we have the B 1+δ,p ∞ -regularity is preserved with the corresponding estimate [Du] B δ,p ∞ (Ωs) . In the case of zero Neumann boundary, we can instead take an even extension given bỹ (4.14) Note this preserves fractional differentiability up to order δ < 1 p , as illustrated by the following lemma.
. This can be deduced from [73, Section 2.8.7], which asserts that the characteristic function χ {xn>0} is a multiplier for B δ,p ∞ (R n ) provided δ < 1 p . Here we instead present an elementary proof, based on a reduction to the one-dimensional case.
Proof. It suffices to verify fractional differentiability in the direction e n . For h > 0, we haveˆB Evidently . For the second term note, that by Fubini we havê , since, by considering a countable dense set of ε ∈ (0, h), we can show that Then for any such x ′ , by Sobolev embedding we have f (x ′ , ·) ∈ B δ,p ∞ ((0, h)) ֒→ L p 1−δp ((0, h)) provided δp < 1. In this case, by Hölder we can splitˆh , so integrating over x ′ ∈ B n−1 1−r we deduce that , from which the result follows.
In the Neumann case, we can use this to perform our iteration to deduce differentiability given by some β > 1 2 for V p,µ (Du) as we will see in Theorem 4.8 below. However beyond this point our extension will no longer preserve the regularity of u, so our iteration scheme will be forced to terminate. To go beyond this range we will need to assume g N = 0 and apply the following lemma, which will be used in the proof of Theorem 1.4.
for all ϕ ∈ W 1,q (Ω), where a 0 : Ω×[0, ∞) → R is continuous satisfying the growth bounds Then, if we additionally have V q,µ (Dv) ∈ B δ,2 ∞ (B + s ) for some δ > 1 2 and all s < 1, we have Proof. By our assumption on V q,µ (Dv), we have Dv admits a trace on Γ s for each s < 1 by [74,Section 4.4.2], so by taking the precise representative this will be understood as Dv. Let denotes the unit ball in R n−1 . In particular, there is s < 1 such that ϕ is supported in B n−1 s . Then for ε sufficiently small we have φ(x ′ ) (ε − x n )1 {xn≤ε} is a valid test function, and using this gives The second integral, we can bound using (2.14) by Since the integral vanishes as ε → 0, it follows that Since Dv ∈ L q (Γ s ) as noted above, we can pass to the limit to get 0 =ˆΓ 1 a 0 (x, |Dv|)Dv · e n ⊗ φ(x ′ ) dx ′ In particular, we deduce that a 0 (x, |Dv|)Dv = 0 H n−1 -almost everywhere on Γ 1 . By (4.15), we deduce that |Dv| q−2 |∂ n v| = 0 H n−1almost everywhere on Γ 1 , and in particular this implies that ∂ n v = 0 on Γ 1 in the sense of traces.

Proof of boundary differentiability
We first prove the improved differentiability in the Neumann case, in the case of non-homogeneous boundary conditions. Theorem 4.8. Suppose 1 < p ≤ q < ∞, α ∈ (0, 1], q < n+α n p, Ω is a bounded C 1,1 -domain and F ≡ F 0 (x, |z|) satisfies (H1)-(H4). Then there is δ 0 > 1 2 , such that the following holds for all δ < min , and u is a relaxed minimiser of F N , then V p,µ (Du) ∈ B β,2 ∞ (Ω) and we have the corresponding estimate The same conclusion holds for some δ > 1 2 (with a modified estimate) if p ≥ 2, α > 1 2 , By a standard straightening of the boundary described in Section 3.1, it suffices to consider the case where Ω = B + 1 and u is a relaxed minimiser for F M with γ = ∂B + 1 \Γ 1 . Note in particular, that it is possible to localise the relaxed functional due to Lemma 3.3. We need to prove improved differentiability in Ω r = B + r for some r < 1. This is achieved by combining Lemma 4.4 with the even extension provided by (4.14). Note that by Lemma 4.6, this preserves B 1+s,p ∞ -differentiability provided s < 1 p . Employing this extension and making the obvious changes to the domains of integration, we can employ the iteration technique of Theorem 1.3 without change using the bound in Lemma 4.4 as a starting point in order to deduce the claimed differentiability V p,µ (Du) ∈ B s,2 ∞ (Ω) for all s < min α max(2,p)−1) , 1 2 . If α max(2,p)−1 > 1 2 , to conclude, we run the iteration step once more with s sufficiently close to 1 2 , to find β 0 > 1 2 such that V p,µ (Du) ∈ B β 0 ,2 ∞ (Ω). The case when f ∈ B δ,p ′ ∞ (Ω) follows by arguing as in Step 4 of Theorem 1.3 with δ = 0, noting that p ′ 2 ≥ 1 2 .
Using this and the results from the previous sections, we can now complete the proof of Theorem 1.4.
Proof of Theorem 1.4. As in Theorem 4.8, we can reduce to the case when Ω = B + 1 . In the case of Dirichlet boundary, taking an odd extension (4.13) we can apply Lemma 4.4 and iterate as in the proof of Theorem 1.3 to conclude. In the Neumann case, the result follows by Theorem 4.8 when α ≤ 1 max{2,p ′ } . Otherwise there is δ > 1 2 such that V p,µ (Du) ∈ B δ,p ∞ (Ω r ) for all r < 1, and we use the fact that g N = 0 to iterate further. Claim: If α > 1 max{2,p ′ } , then we have ∂ n u = 0 on Γ. We will establish this for minimisers u ε of the relaxation from Section 2.6, which in this localised form minimises over v ∈ W 1,q (B + 1 ) such that v = u on γ = ∂B + 1 \ Γ 1 . Then each F ε satisfies the hypotheses of Theorem 4.8 with p = q, so given our assumption on α we know that V p,µ (Du ε ) ∈ B δ,2 ∞ (B + s ) for some δ > 1 2 and s < 1. Hence by Lemma 4.7 it follows that ∂u ε = 0 on Γ s for each s < 1. Now passing to the limit, noting u ε converges to u in W 1,p (B + 1 ), the claim follows. Hence, applying [74, Section 4.5.2], we have an even extension of u preserves the B s,2 ∞ -norm of V p,µ (Du) for all s < 1 since ∂ n u is odd and vanishes on Γ 1 . Now continuing to iterate as in Theorem 1.3, the result follows. Remark 4.9. Due to Corollary 4.1 and the corresponding result in the Dirichlet case in [49], if F is autonomous, Theorem 1.4 applies if 2 ≤ p ≤ q < min p + 1, np n−1 .

Mixed boundary problems
We will conclude this section with a brief discussion of mixed boundary problems of the form min For convex functionals with p-growth, techniques using first-order difference quotients have been developed in [64], and we will employ similar arguments. The relaxed functional to use in this case is where v ∈ X = {u ∈ W 1,p (Ω) : Tr u = g D on Γ D } and Y = W 1,q (Ω) ∩ X. Using this definition it is straightforward to adapt the results of Section 2.6 to this setting. We wish to employ similar arguments as in the pure Dirichlet and Neumann cases, however we will need to ensure our difference quotients preserve the respective boundary conditions. This will involve considering domains with piecewise regular boundary, such that the respective boundary conditions are prescribed on the regular portions. More precisely we will consider C k,α -domains with corners, in that every x 0 ∈ Ω admits a C k,α -diffeomorphism to a neighbourhood of the closure of the model space We write B ++ = R n ++ ∩ B 1 (0). Moreover we will impose further restrictions on the boundary components, and assume one of the following two cases are satisfied at each x 0 ∈ Γ D ∩ Γ N . We say that we are in the flat case, if after localising, we may use a C 1,α -flattening map such that Γ D , Γ N are mapped to {x n−1 > 0, x n = 0} and {x n−1 < 0, x n = 0}, respectively. We say we are in the corner case, if after localising and flattening, we may work in B k,++ := B 1 (0) ∩ {x i > 0, ∀k ≤ i ≤ n} and moreover, Γ D is mapped to ℓ≤i≤n ∂B ++ ∩ {x i = 0} for some k, ℓ with 1 ≤ k ≤ n and k ≤ ℓ ≤ n + 1 (understanding that Γ D = ∅ when ℓ = n + 1). The two cases are illustrated in Figure 2 below. Note that in the corner case we may take k = ℓ or ℓ = n + 1, which corresponds to prescribing purely Dirichlet and Neumann boundary conditions in piecewise regular domains respectively.
Under these assumptions, we obtain the following analogue to Theorem 1.2.
Let Ω be a C 1,1 -domain with corners, and let Γ D , Γ N ⊂ ∂Ω be such that locally we are either in the flat or corner case. Let F satisfy (H1)-(H3), with In particular, u ∈ W 1,q (Ω).
Since our estimates are local in nature, it suffices to consider the case x 0 ∈ Γ D ∩ Γ N where the two boundary conditions meet. The only aspect in which the proof differs from [48] and Theorem 1.2 is in the construction of the difference quotients; as such we will only outline the necessary changes to treat the mixed case. Informally we need to ensure the difference quotients do not move the unconstrained regions (the Neumann part) into the constrained regions (the Dirichlet part). This is achieved by flattening the boundary to reduce to the model case, however it is unclear whether this is really necessary; in particular we must assume that the boundary is C 1,1 -regular to preserve both the Hölder regularity of the coefficients and Besov-regularity of functions under the flattening map, but it will be interesting to understand whether this can be relaxed.
Proof of Theorem 4.10. We will derive estimates to the minimia u ε ∈ W 1,q (Ω) associated to F ε (·) for each ε > 0, and pass to the limit at the end. Upon flattening the boundary, we will consider the flat and corner cases separately. In both cases, we will take a B 1+α,p ∞ extensiong of g to B 1 (0), and consider v ε = u ε −g. We will need to extend this to the full ball B 1 (0), and the respective procedures are illustrated in Figure 2.  In the flat case, we will then take a zero extension of v ε to {x n−1 > 0, x n < 0}, before taking a W 1,q extensionṽ ε to the full ball B 1 (0). By constructionṽ agrees with v on B + , and vanishes on Γ D (and also on {x n−1 = 0, x n < 0}). From here we observe that we can take difference quotients in directions h ∈ R n such that h n−1 < 0 and h n > 0 (see Figure 2(a)). Since v h vanishes on Γ D we have T h v is a valid competitor, and we can argue in the same way as before.
In the corner case, we will take a first-order reflection along each x i for k ≤ i < ℓ (as in Theorem 4.5.2 in [74], which is a W 1,q -extension). This gives a W 1,q function on B ℓ,++ vanishing on ∂B k,++ ∩ B 1 (0),, so we can extend it to the full ball by means of a zero extension to obtainṽ. Now we can take difference quotients in directions h ∈ R n such that h i > 0 for each ℓ ≤ i ≤ n (see Figure 2(b) for the case k = n − 1, ℓ = n), and argue as before to conclude.
Assuming F ≡ F 0 (x, |z|) is a radial integrand and homogeneous boundary conditions, we may also improve on Theorem 1.4 as in Section 4.2. In addition to assuming (H4), we also found it was necessary to work on piecewise regular domains. Here the idea is that the reflection procedure in Section 4.3 can be carried out separately in each x i direction, taking odd and even reflections for the case of Dirichlet and Neumann boundaries respectively. We point out that this result also shows our regularity results hold for piecewise C 1,1 domains, even in the purely Dirichlet and Neumann cases. (Ω) for some β ∈ [α, 2] then V p,µ (Du) ∈ B δ,2 ∞ (Ω) for all δ < min p ′ β 2 , α . Moreover, if a-priori u ∈ W 1,q (Ω), it suffices to ensure q < np n−α .
Proof. By flattening we can reduce one again to the corner ball B k,++ , where Γ D is mapped to ℓ≤i≤n {x i = 0} and Γ N is mapped to k≤i<ℓ {x i = 0}. Here Γ D , Γ D are understood to be empty if ℓ = n+1, k = ℓ respectively. We know that u ∈ B 1+ α p ,p ∞ (B k,++ ) due to Theorem 4.10. Now we extend u to the full ball by taking an even reflection along x i for k ≤ i < ℓ, and an odd reflection along x i for ℓ ≤ i ≤ n. From Section 4.3.2 we see the Besov regularity is preserved by these extensions, adapting Lemmas 4.6 and 4.7 for the case of an even extension. Therefore we may iteratively apply the difference quotient argument from Lemma 4.4 in each direction, and argue as in the proof of Theorem 1.4 to conclude.
Remark 4.12. In both Theorem 4.10 and Theorem 4.11, if F is independent of x, due to Corollary 4.1, it suffices to assume 2 ≤ p ≤ q < min np n−1 , p + 1 .
Remark 4.13. In order to obtain Theorem 4.11, we genuinely must assume that Γ D and Γ N meet at a corner. This was already observed in [63] for linear equations, where the function u(x, y) = ℑ(x + iy) 1 2 is seen to be harmonic on the upper-half plane {ℑz > 0}, and satisfies a mixed Dirichlet-Neumann condition on the boundary (u vanishes on {x = 0, y > 0}, and ∂ y u vanishes on {x = 0, y < 0}). Note that u does not lie in W 2,r loc (R 2 + ) for any r ≥ 4 3 , whereas the above results, if they apply, would imply u lies in W 2,r loc (R 2 + ) for all 1 ≤ r < 2. Similarly, we have already seen that Γ D must be sufficiently regular to deduce higher differentiability, so it must map to one of the corner faces under the flattening map.

Characterisation of regular boundary points
In this section, we will establish the following ε-regularity result.
Theorem 1.7. Assume F : Ω × R m×n → R is C 2 in z, satisfies (H1)-(H4) with µ = 1, and the exponents satisfy We will assume further that Ω ⊂ R n is a bounded C 1,α domain and f ∈ L n 1−α (Ω, R m ). For the boundary function we will assume g D ∈ C 1,α (Ω, R m ) in the case of Dirichlet boundary, and g N ∈ C 0,α (Ω, R m ) in the case of Neumann boundary, in which case we also assume the compatibility condition (1.1). Suppose u is a relaxed minimiser of (N) or (D). Then for each M > 0 and 0 for some x 0 ∈ Ω and 0 < R < R 0 , then u is Throughout this section we will assume the above assumptions hold. Additionally for z 0 ∈ R m×n , we introduce the shifted integrand This satisfies the estimates for all x 1 , x 2 ∈ Ω and z ∈ R m×n , which follows by distinguishing between when |z| ≤ 1 and |z| > 1. Here the constants depend on M > 0 where |z 0 | ≤ M.
Following [27], we will employ a blow-up argument to establish a suitable excess decay estimate. In the case of Dirichlet boundary, one can also argue directly by means of an A-harmonic approximation similarly as in [52] and we intend to return to direct arguments, also in the case of Neumann boundary, in future work.
Remark 5.1. While our ε-regularity results are local in nature, for our arguments to hold it is necessary to impose only one type of boundary condition on ∂Ω; in particular this does not apply to the problems considered in Section 4.4. This is because the proof will make use of the translation invariance u → u + c of the Neumann problem, which is not available in the mixed setting. It is unclear whether this is merely a technical difficulty, but it highlights that we crucially use a global property in the Neumann case.
We will locally flatten the boundary, following the procedure in Section 3.1. In what follows, we will denote B + = B 1 (0)∩R n + and Γ = B 1 (0) ∩ {x n = 0}. Finally, for x 0 ∈ B + and R > 0 we will denote B +

Caccioppoli-type inequalities
We now prove the following boundary Caccioppoli estimate in the Neumann case.
, and the implicit constant depends on M, F, Ω and other structural constants indicated in Section 2.1.
Proof. We will putF = F Da for the shifted integrand (5.1), as well asũ = u − a. Note thatũ is a relaxed minimiser of the problem where we have written g = g N . Indeed, if we let F 0,N denote the relaxed functional associated to v →´ΩF (x, Dv)), by definition of the relaxation we see that Let R/2 < r < s < R/2, applying Lemma 3.1 to u gives w and r < r ′ < s ′ < r. Let φ be a smooth cut-off supported on B s ′ , with φ = 1 on B r ′ and |Dφ| ≤ c s−r . We will set w = w − a and ψ = (1 − φ)w. Then by Remark 3.2 and Lemma 3.3 applied to F 0,N , we have and a similar statement holds for ψ. Further, by lower-semicontinuity of the energý Ω |V p (Dv)| 2 dx and (5.4) we havê Hence using the above and minimality ofũ we havê We find using (5.2) and (2.10), Hence using (3.1), (3.2) we can bound For I 2 , we note thatũ − ψ vanishes on Ω ∩ ∂B s ′ , allowing us to use a Poincaré-Sobolev inequality to estimate This is justified by integrating along one of the tangential directions. To estimate I 3 , we use the translation invariance of the problem (due to the compatibility condition (1.1)) to replaceũ byũ + b, which in turn replacesw byw + b. Then we can choose b ∈ R so that 0 =ˆB (5.6) Using this we can estimate For I 4 , we use (5.5), the divergence theorem, (5.6), and the assumption that ∂Ω is C 1,α to bound Now using Young's inequality (2.9) applied to |V p (z)| 2 we have Collecting estimates we have shown that Filling the hole, applying a standard iteration argument (see for instance [40, Lemma 6.1]) and choosing δ > 0 sufficiently small to absorb the final gradient term, we deduce the result.
The Dirichlet case is essentially the same, and we will only outline the necessary modifications.
Lemma 5.3. In the setting of Theorem 1.7, suppose that u ∈ W 1,p g D (Ω) is a relaxed minimiser to (D) and let x 0 ∈ Ω, 0 < R < R 0 . Then for each M > 0, if a : R n → R m is an affine map such that |Da| ≤ M , , and the implicit constant depends on M, F, Dg D L ∞ (Ω) , Ω and other structural constants.
Proof. We will choose R 0 > 0 sufficiently small so B R 0 (x 0 ) can be flattened by a C 1,αdiffeomorphism, which we will use to define the smooth operator from Lemma 3.1.
Note that the shifted estimates (5.2)-(5.4) continue to hold, however the constants depend on M ≥ |z 0 | + Dg L ∞ (Ω) . For the Hölder bound we can use (5.3) and also local uniform estimates for ∂ 2 z F to obtain Now for R 2 ≤ r < s < R, applying Lemma 3.1 to u − g gives v and r ≤ r ′ < s ′ ≤ s for which the claimed estimates hold. We then setũ = u − a −g andw = v − a − (g −g), which satisfies similar estimates as we have only shifted v by an affine map. As before let φ be a cutoff supported in B s ′ such that χ ≡ 1 on B r ′ and |Dφ| ≤ C s ′ −r ′ . Setting ψ = (1 − φ)w, using (5.4), Lemma 3.3, Remark 3.2 and minimality of u we havê We can estimate I 1 and I 2 as before. For I 3 we can use the Hölder estimate (5.7) and the fact thatũ − ψ vanishes on ∂Ω s ′ (x 0 ) to show that from which the rest follows as in the Neumann case.

Excess decay estimates
We will obtain, by means of a blow-up argument, estimates for the excess energy where β < α is fixed. Throughout this section we will assume we are in the case of a flat boundary where Ω R (x 0 ) = B + R (x 0 ).
Lemma 5.4. In the setting of Theorem 1.7, assume u is a relaxed minimiser of (N) or Proof in the Neumann case. We will employ a blow-up argument, which involves several steps. As before write g = g N .
Step 1 (Blow up): Suppose otherwise, so there exists x j ∈ B + 1/2 and R j > 0 such that with each λ j ≤ 1, and, for C M > 0 to be chosen appropriately, we have To simplify notation, put

and introduce in addition
so in particular we have Note that B + ⊂ B + j ⊂ B for any j and (v j ) B + j = 0, so by Lemma 2.8 we deduce that {v j } is bounded in W 1,min{p,2} (B + ). Also, by passing to a subsequence we can assume that x j → x 0 , and moreover that x 0 ∈ Γ; indeed since R 2β for j sufficiently large, and we can apply the interior argument from [27] to derive a contradiction. Also by passing to a subsequence we have A j → A 0 for some A 0 .
Step 2 (Extremality of v j ): We introduce the rescaled integrands where F A j is given by (5.1). It is convenient to introduce f j (y) = R j λ −1 j f (x j + R j y) and g j (y) = λ −1 j g(x j + R j y). Then by Lemma 2.15, v j satisfieŝ We will additionally requiré Γ j ϕ j dH n−1 = 0, so we can estimate where we have used the α-Hölder continuity of F. For the second term we apply the Sobolev inequality to estimate More precisely we used Lemma 2.8, noting this gives an extra term which we estimate by , integrating along one of the tangential directions. Finally, for I 3 , we use the cancellation condition and Gagliardo's trace theorem (noting Γ j is flat) to estimate Combining the estimates we conclude that Step 3 (Limit is harmonic): Passing to a subsequence, we obtain a weak limit v j ⇀ v in W 1,min{2,p} (B + , R N ). We will show this satisfies a constant-coefficient equation. Let . We claim that we may extend ϕ to ϕ j defined on B + j so that ϕ j = 0 on B + j \ Γ j and moreover ϕ j as well as´ Γ j ϕ j = 0, where the implicit constant is independent of j. In addition, we may ensure that ϕ j → ϕ in W 1,∞ (B + ). We postpone the construction of ϕ j to Lemma 5.5 after the proof. We now wish to send j → ∞ in (5.8), for which we split By Markov's inequality we have So far we have assumed that ϕ ∈ W 1,∞ (B + ), but we claim we can take ϕ ∈ W 1,p (B + ) by density. We will prove this in Lemma 5.6 after the proof. Now, by standard energy methods, we deduce from (5.9) that v is C 1 (B + τ ) for each τ ∈ (0, 1) with the associated estimate − One needs to be slightly careful in the choice of test function however; we typically take ϕ h = ∆ k −he i (η 2 ∆ k he i v) for k ≥ 1, 1 ≤ i ≤ n − 1 and a suitable cutoff η. However, u is only defined up to a constant. In particular, changing u → u + λ j R j b, changes v j → v j + b and hence v → v + b. Note that these changes do not impact any of our estimates. Thus, with the choice b h = −(∆ k −he i (η 2 ∆ k he i v) Γ ), ϕ h is a valid test function. The rest of the argument is standard; see for instance [40,Chapter 10].
Step 4 (Conclusion): We now consider the affine maps · y approximating u and v j respectively, where B + j,2τ = R −1 j (B + 2τ R j (x j ) − x j ). Then by the Caccioppoli-type inequality (Lemma 5.2) we have Note that R 2β j λ 2 j ≤ 1, and sending j → ∞ we have Also by the Poincaré inequality (Lemma 2.8) we can estimate for 1 ≤ r ≤ n n−1 . Taking r = 1 and noting that λ 2 p (q−p) j → 0 as j → 0, it follows that lim sup To pass to the limit in the first term, observe by (5.11) that the sequence is uniformly bounded in L n n−1 (B 2τ ) and hence is uniformly integrable. Now by strong L p convergence we have v j → v a.e. in B + 2τ andã j,2τ →ã(y) = (v) B + (τ ) + (Dv j ) B + τ · y uniformly in B + 2τ , so by Vitali's convergence theorem we deduce that lim sup Hence there is C M > 0 for which for j sufficiently large, and so choosing C M > C M gives the desired contradiction.
We now prove the following Lemma, which guarantees the existence of maps {φ j } with the properties claimed in Step 3 of the previous proof. Moreover, ϕ j → ϕ in W 1,∞ (B + ).
Proof. Consider the cone C j , which is uniquely determined by having sections Γ and Γ j . Denote by C ′ j = C j ∩ {x n ≤ 1}. We construct a family of invertible, uniformly bi-C 2 -maps f j : and further f j → id in C 2 (B + 1 ) f −1 j → id in C 2 (f j (B + 1 )). (5.12) Once, {f j } is constructed, we extend ϕ by 0 to R n−1 + and set Using the properties of {f j } (see Fig.3), we find that with implicit constants independent of j. It remains to construct {f j }. If B + j = B + , we may set f j = id. Hence we may assume h j = x j,n > 0. Note that then B + j = B 1 ∩ {x n > −h j }. It is straightforward to check that C j is centered at −a where a = Note that as h → 0, a → ∞. Thus, ensuring that h j < c 0 , we may assume that a−h > 1. Then x n + a > 1 and a + y n > 1, so that f j and f −1 j are smooth. Moreover, using again that a → ∞ as h → 0, it is a straightforward calculation to check that the desired convergence (5.12) holds.
It remains to show the density statement we used in the proof of Lemma 5.4.
is dense in with respect to the W 1,p norm.
We will conclude by sketching the necessary modifications to the proof of Lemma 5.4 in the Dirichlet case.
Proof of Lemma 5.4 in the Dirichlet case. We will establish a decay estimate for the modified excess We suppose otherwise, then we can find x j ∈ B + 1/2 and R j > 0 such that with each λ j ≤ 1, and for C M > 0 to be chosen we have Letg(x) = g(x) − g(x 0 ) − Dg(x 0 ) · (x − x 0 ), we will consider the rescaled sequence v j (y) = u(x j + R j y) − a j (y) −g(x j + R j y) λ j R j .
Arguing as in the Neumann case, we see that v j is bounded in W 1,min{2,p} (B + ), and we can assume that x j → x 0 ∈ Γ and that R j → 0. We now use the fact that each v j is extremal with respect to the shifted integrand F j (y, z) = λ −2 j F A j +Dg(x) (x j + R j y, λ j z), to deduce the estimate 0 ≤ˆ B + j ∂ z F j (y, Dv j ) · Dϕ j dx + CR α j λ jˆ B + j |Dϕ j | dx for any ϕ j ∈ W 1,∞ 0 ( B + j ). Here we are using Lemma 2.15. This will involve estimating 1 λ 2 jˆ B + j ∂ z F (x 0 + R j y, λ j Da + λ j Dg(x)) · Dϕ j dx ≤ CR α j λ jˆ B + j |Dϕ j | dx, using (5.7) and the fact that ϕ j vanishes on ∂ B + j . The remainder of the proof is analogous to the Neumann case. We note that for ϕ ∈ W 1,∞ 0 (B + ), we can simply extend it by zero to B + j to obtain ϕ j , from which we obtain the limit map v is harmonic. Moreover since v is affine on Γ, we deduce the same the decay estimate (5.10) as in the Neumann case using results in [40,Chapter 10]. Finally applying Lemma 5.3 we can argue analogously as in Step 4 to obtain a contradiction. Hence we infer a decay estimate of the form The corresponding estimate for E(x 0 , R) follows by (2.7) and noting that − B R |V p (Dg)| dx ≤ C|V p (R β )| 2 ≤ CR 2β for R ≤ 1.

Iteration of excess and conclusion
We can now conclude by a standard iteration argument. Note that we can treat the Dirichlet and Neumann problems simultaneously, since we established the same excess decay estimate.
By shrinking ε > 0 further, if necessary, we claim that Noting that e −1 p (t) ∼ √ t for t > 0 sufficiently small, we can estimate which is less than 2 n+2 M if ε > 0 is sufficiently small. This shows that E(x, r) ≤ C r R β for all x ∈ B + R/2 (x 0 ) and 0 < r < R, verifying the Campanato-Meyers characterisation of Hölder continuity (see for instance [40,Theorem 2.9]).
To deduce Corollary 1.9, straightening the boundary and employing a covering argument, it suffices to consider the case where Ω = B + 1 . We note that if v ∈ W θ,p (Ω, R N ), we have, c.f. [58], dim H x ∈ Ω : lim sup